UFL

From The FEniCS project

UFL (Unified Form Language) is a unified language for definition of variational forms intended for finite element discretization. More precisely, it defines a fixed interface for choosing finite element spaces and defining expressions for weak forms in a notation close to mathematical notation. Examples of form compilers that are planned to use the UFL language are FFC and SyFi.

At this stage of the project we are defining what we want to include in the language. All wanted features should be added as tests in the repository. Currently the prototype of the UFL language consists of a python package.

Finite element definitions

TODO: Define all variations of finite element choices that we want

Element families:

• family = "Lagrange" | "CG"
• family = "Discontinuos Lagrange" | "DG"
• family = "Crouzeix-Raviart" | "CR"
• family = "Brezzi-Douglas-Marini" | "BDM"
• family = "Brezzi-Douglas-Fortin-Marini" | "BDFM"
• family = "Raviart-Thomas" | "RT"
• family = "Nedelec" (FIXME: how to handle 1st/second kind here?)
• family = "Arnold-Winther" | "AW"
• family = "Bubble"
• familu = "Quadrature" | "Q" (Kristian's quadrature elements in FFC)
• family = "..." # TODO: someone extend this list

Supported polygons (same as UFC):

• polygon = "interval"
• polygon = "triangle"
• polygon = "quadrilateral"
• polygon = "tetrahedron"
• polygon = "hexahedron"

Cell definitions (including high order geometries): # TODO: make this general and clean

• cell = polygon
• cell = Cell(polygon, order, ...)

Alternatives for basis function mapping:

• mapping = "affine"
• mapping = "isoparametric"
• mapping = "Piola" (FIXME: how to handle covariant/contravariant/symmetric Piola here?)

Element definitions:

• element = FiniteElement(family, cell, order, mapping="affine")
• element = VectorElement(family, cell, order, mapping="affine", size=None)
• element = TensorElement(family, cell, order, mapping="affine", shape=None, symmetric=False)

Element definitions with mapping:

• element = FiniteElement(family, cell, order, mapping="affine")
• element = VectorElement(family, cell, order, mapping="affine", size=None)
• element = TensorElement(family, cell, order, mapping="affine", shape=None, symmetric=False)

Mixed elements:

• mixed_element = MixedElement(element1, element2[, element3, ...])
• mixed_element = element1 + element2 # equivalent with MixedElement(element1, element2)

Terminal operands

Operands that depend only on element spaces and decorative names.

Basis functions

• u = BasisFunction(element)
• v = TrialFunction(element)
• u = TestFunction(element)
• u, p = BasisFunctions(mixed_element): basisfunctions over the subelements of a mixed element
• u, p = TrialFunctions(mixed_element): trialfunctions over the subelements of a mixed element
• v, q = TestFunctions(mixed_element): testfunctions over the subelements of a mixed element

Coefficient functions

These suggestions differs a bit from the current FFC language:

• f = Function(element, name): name could be used to generate easier to read code, or latex, and we can add a function to get this name from the compiled UFC form?
• c = Constant(polygon, name): the polygon could technically be deduced from the UFL Form by the form compiler, but having it here allows a more rigorous error checking during expression definition.

Predefined objects

• I = Identity(): TODO: Take the rank as argument? This will usually be the rank 2 identity matrix, but may mean other things.
• n = FacetNormal()
• h = MeshSize(): TODO: Define what this means, maybe we need several different mesh size variables?
• default indices i-s: i = Index("i") etc.
• default time t
• default coordinates x[0], x[1], x[2]
• default integral operators dx, ds, dS, dx#, ds#, dS#: dx0 = IntegralOperator("cell", 0) etc.

Language operators to include

TODO: Define all language operators we want.

Scalar operators

• a+b: scalar addition
• a-b: scalar subtraction
• a*b: scalar multiplication
• a/b: scalar division
• a**b: scalar a to the power of scalar b

Vector operators

• v+w: vector addition
• v-w: vector subtraction
• v*a, a*v: multiplication of all components in vector v with scalar a
• v/a: division of all components in vector v with scalar a
• v[i]: component access where i is an integer or Index

Matrix operators

• A+B: matrix addition
• A-B: matrix subtraction
• A*a, a*A: multiplication of all components in matrix A with scalar a
• A/a: division of all components in matrix A with scalar a
• A[i,j]: component access where i and j are integers or Index objects
• A[i,:]: row access, rank(A[i,:]) == 1
• A[:,j]: column access, rank(A[:,j]) == 1

Index notation

TODO: This needs more work!

• ii = Index("ii"): define a new index
• Einstein summation convention: TODO: define this convention here
• rank(f): the number of free indices in f
• f[i]: component access where i is an index tuple, where rank(f) == len(i), and rank(f) == 0
• f[..., i]: component access where i is an index tuple, the rank r of f[i] is rank(f) - len(i), the first r indices of f are free
• f[i, ...]: component access where i is an index tuple, the rank r of f[i] is rank(f) - len(i), the last r indices of f are free
• delta(i, j): Kronecker delta / identity matrix

Compound operators

These are operators that can be defined in terms of the basic operators.

• dot(a, b): contraction of the last index of a and the first index of b, rank(dot(a,b)) == rank(a)+rank(b)-2
• outer(a, b): outer product of a and b, rank(outer(a,b)) == rank(a) + rank(b)
• inner(a, b): contraction of all indices of a and b, they must have the exact same dimensions, rank(inner(a,b)) == 0
• transpose(A): transpose of matrix A, transpose(A)[i,j] == A[j,i]
• trace(A): trace of matrix A
• det(A): determinant of matrix A
• dev(A): deviatoric part of matrix A
• inverse(A): inverse of matrix A (limited to 1x1, 2x2 and 3x3 tensors)
• cofactor(A): cofactor of matrix A, cofactor(A) == det(A)*inverse(A) (limited to 1x1, 2x2 and 3x3 tensors)

Differential operators

• v.dx(i): partial derivative in coordinate direction i
• Dx(v, i): alternative notation for partial derivative # NB!: changed this from D to Dx, is this ok?
• v.dt(): time differentiation
• Dt(v): alternative notation for time differentiation
• oD(v, s): ordinary derivative of f w.r.t. s # TODO: choose a better name than oD, and define what s can be, parametric symbol of some kind
• pD(v, s): partial derivative of f w.r.t. s # TODO: choose a better name than pD, and define what s can be, parametric symbol of some kind

• grad(f): gradient of expression f
• div(f): divergence of expression f
• curl(f): curl of expression f
• should we have curl for scalars and rot for vectors in 2D?
• wedge(f, g): wedge product of expressions f and g

Basic functions

These nonlinear scalar functions should be self explaining. More functions like these can be added ad lib, f.ex. matching the functions found in the standard C library.

• abs(f)
• sqrt(f)
• pow(f, a)
• exp(f)
• ln(f)
• log(f,base)
• sin(f)
• cos(f)
• tan(f)

Integration

• (...) * dx: integrate (...) over all cells
• (...) * ds: integrate (...) over all exterior facets (the boundary)
• (...) * dS: integrate (...) over all interior facets (all facets not on the boundary)

• (...) * dxN where N is in (0-9): integrate (...) over all cells in cell subdomain N (subset of domain)
• (...) * dsN where N is in (0-9): integrate (...) over all exterior facets in exterior facet subdomain N (subset of the boundary)
• (...) * dSN where N is in (0-9): integrate (...) over all interior facets in interior facet subdomain N (subset of facets not on the boundary)

Subdomains 0-9 of cells, exterior facets, and interior facets have no correlation, ie subdomain k of the exterior facets does not have to be the boundary of subdomain k of the cells.

Transformations of entire forms

Automated computation of the derivative of a form with respect to a coefficient: If f is a rank 0 form that is nonlinear in coefficient w, F is the rank 1 derivative of f with respect to w and J is the rank 2 Jacobi of F with respect to w.

Example:

 w = Function(element, "w")
 f = (w**2)/2 * dx
 F = Derivative(f, w)
 J = Derivative(F, w)

is equivalent to:

 w = Function(element, "w")
 f = (w**2)/2 * dx

 v = BasisFunction(element)
 F = w*v * dx

 u = BasisFunction(element)
 F = u*v * dx

Automated Action computation:

 L = Action(a)

If a is a rank 2 form used to assemble the matrix A, L is a rank 1 form that can be used to assemble the vector b = Ax. Useful to define both the form of a matrix and the form of its action, and for the action of a Jacobi matrix computed using Derivative.

Extraction of left- and right-hand sides:

If

 v = TestFunction(element)
 u = TrialFunction(element)
 u0 = Function(element)

then

 F = v*(u - u0)*dx + k*dot(grad(v), grad(0.5*(u0 + u)))*dx
 a = lhs(F)
 L = rhs(F)

should be equivalent to

 a = v*u*dx + 0.5*k*dot(grad(v), grad(u))*dx
 L = v*u0*dx - 0.5*k*dot(grad(v), grad(u0))*dx

Forms

• Determine whether or not a bilinear form is self-adjoint (symmetric)

Low priority ideas

A couple of fairly complicated low priority ideas here. These likely won't be part of UFL in the first version, perhaps never.

• A way to plug in low level user functions at the source code level?
• Conditionals: depending on some condition, assign an expression to a variable that can be used in further expressions.

Conditionals are used f.ex. in a lot of turbulence models, but can often be handled as Functions that are computed elsewhere. Alternatively, functions like min(a, b) and max(a, b) should cover some uses if conditionals. Both of these make features like differentiation difficult or impossible, and are probably only usable with quadrature.

Main authors

Martin Alnæs, Kent-Andre Mardal, Anders Logg

License

The UFL specification is released into the public domain. FIXME: But the implementation is GPL?

Dependencies and requirements

The UFL python code does not specify any additional dependencies.