Interpolation

Firedrake offers various ways to interpolate expressions onto fields (Functions). Interpolation is often used to set up initial conditions and/or boundary conditions. The basic syntax for interpolation is:

# create new function f on function space V
f = interpolate(expression, V)

# alternatively:
f = Function(V).interpolate(expression)

# setting the values of an existing function
f.interpolate(expression)

Warning

Interpolation currently only works if all nodes of the target finite element are point evaluation nodes.

The recommended way to specify the source expression is UFL. UFL produces clear error messages in case of syntax or type errors, yet UFL expressions have good run-time performance, since they are translated to C interpolation kernels using TSFC technology. Moreover, UFL offers a rich language for describing expressions, including:

  • The coordinates: in physical space as SpatialCoordinate, and in reference space as ufl.geometry.CellCoordinate.
  • Firedrake Functions, derivatives of Functions, and Constants.
  • Literal numbers, basic arithmetic operations, and also mathematical functions such as sin, cos, sqrt, abs, etc.
  • Conditional expressions using UFL conditional.
  • Compound expressions involving any of the above.

Here is an example demonstrating some of these features:

# g is a vector-valued Function, e.g. on an H(div) function space
f = interpolate(sqrt(3.2 * div(g)), V)

Interpolation from external data

Unfortunately, UFL interpolation is not applicable if some of the source data is not yet available as a Firedrake Function or UFL expression. Here we describe a recipe for moving external to Firedrake fields.

Let us assume that there is some function mydata(X) which takes as input an \(n \times d\) array, where \(n\) is the number of points at which the data values are needed, and \(d\) is the geometric dimension of the mesh. mydata(X) shall return a \(n\) long vector of the scalar values evaluated at the points provided. (Assuming that the target FunctionSpace is scalar valued, although this recipe can be extended to vector or tensor valued fields.) Presumably mydata works by interpolating the external data source, but the precise details are not relevant now. In this case, interpolation into a target function space V proceeds as follows:

# First, grab the mesh.
m = V.ufl_domain()

# Now make the VectorFunctionSpace corresponding to V.
W = VectorFunctionSpace(m, V.ufl_element())

# Next, interpolate the coordinates onto the nodes of W.
X = interpolate(m.coordinates, W)

# Make an output function.
f = Function(V)

# Use the external data function to interpolate the values of f.
f.dat.data[:] = mydata(X.dat.data_ro)

This will also work in parallel, as the interpolation will occur on each process, and Firedrake will take care of the halo updates before the next operation using f.

C string expressions

Warning

This is a deprecated feature, but it remains supported for compatibility with FEniCS.

The Expression class wraps a C string expression, e.g. Expression("sin(x[0]*pi)"), which is then copy-pasted into the interpolation kernel. Consequently, C syntax rules apply inside the string, with the following “environment”:

  • The physical spatial coordinate is available as “vector” x (array in C), i.e. the coordinates x, y, and z are accessed as x[0], x[1], and x[2].
  • The mathematical constant \({\pi}\) is available as pi.
  • Mathematical functions from the C header math.h.

It is possible to augment this environment. For example, Expression('sin(x[0]*t)', t=t) takes the value from the Python variable t, and uses that value for t inside the C string.

Since C string expressions are deprecated, below are a few examples on how to replace them with UFL expressions:

# Expression:
f = interpolate(Expression("sin(x[0]*pi)"), V)

# UFL equivalent:
x = SpatialCoordinate(V.mesh())
f = interpolate(sin(x[0] * math.pi), V)

# Expression with a Constant parameter:
f = interpolate(Expression('sin(x[0]*t)', t=t), V)

# UFL equivalent:
x = SpatialCoordinate(V.mesh())
f = interpolate(sin(x[0] * t), V)

Python expression classes

Warning

This is a deprecated feature, but it remains supported for compatibility with FEniCS.

One can subclass Expression and define a Python method eval on the subclass. An example usage:

class MyExpression(Expression):
    def eval(self, value, x):
        value[:] = numpy.dot(x, x)

    def value_shape(self):
        return ()

f.interpolate(MyExpression())

Here the arguments value and x of eval are numpy arrays. x contains the physical coordinates, and the result of the expression must be written into value. One must not reassign the local variable value, but overwrite its content.

Since Python Expression classes expressions are deprecated, below are a few examples on how to replace them with UFL expressions:

# Python expression:
class MyExpression(Expression):
    def eval(self, value, x):
        value[:] = numpy.dot(x, x)

    def value_shape(self):
        return ()

f.interpolate(MyExpression())

# UFL equivalent:
x = SpatialCoordinate(f.function_space().mesh())
f.interpolate(dot(x, x))