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
# 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)
Interpolation is supported for most, but not all, of the elements that Firedrake provides. In particular, higher-continuity elements such as Argyris and Hermite do not presently support interpolation.
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:
Literal numbers, basic arithmetic operations, and also mathematical functions such as
Conditional expressions using UFL
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)
This also works as expected when interpolating into a a space defined on the facets of the mesh:
# where trace is a trace space on the current mesh: f = interpolate(expression, trace)
Firedrake is also able to generate reusable
objects which provide caching of the interpolation operation. The
following line creates an interpolator which will interpolate the
current value of expression into the space V:
interpolator = Interpolator(expression, V)
If expression does not contain a
the interpolation can be performed with:
f = interpolator.interpolate()
Alternatively, one can use the interpolator to set the value of an existing
f = Function(V) interpolator.interpolate(output=f)
w = TestFunction(W) interpolator = Interpolator(w, V)
Here, interpolator acts as the interpolation matrix from the
FunctionSpace() W into the
FunctionSpace() V. Such that if f is a
Function in W then g = interpolator.interpolate(f) is its
interpolation into a function g in V. As before, the output parameter can
be used to write into an existing
Function. Passing the
transpose=True option to
cause the transpose interpolation to occur. This is equivalent to the
multigrid restriction operation which interpolates assembled 1-forms
in the dual space to V to assembled 1-forms in the dual space to
The interpolation API supports interpolation between meshes where the target function space has finite elements (as given in the list of supported elements)
Lagrange/CG (also known a Continuous Galerkin or P elements),
Q (i.e. Lagrange/CG on lines, quadrilaterals and hexahedra),
Discontinuous Lagrange/DG (also known as Discontinuous Galerkin or DP elements) and
DQ (i.e. Discontinuous Lagrange/DG on lines, quadrilaterals and hexahedra).
Vector, tensor and mixed function spaces can also be interpolated into from other meshes as long as they are constructed from these spaces.
The list of supported elements above is only for target function spaces. Function spaces on the source mesh can be built from most of the supported elements.
There are few constraints on the meshes involved: the target mesh can have a different cell shape, topological dimension, or resolution to the source mesh. There are many use cases for this: For example, two solutions to the same problem calculated on meshes with different resolutions or cell shapes can be interpolated onto one another, or onto a third, finer mesh, and be directly compared.
The target mesh for a cross-mesh interpolation need not cover the full domain
of the source mesh. Volume, surface and line integrals can therefore be
calculated by interpolating onto the mesh or
immersed manifold which defines the volume,
surface or line of interest in the domain. The integral itself is calculated
assemble() on an approriate form over the target mesh
# Start with a simple field exactly represented in the function space over # the unit square domain. m = UnitSquareMesh(2, 2) V = FunctionSpace(m, "CG", 2) x, y = SpatialCoordinate(m) f = Function(V).interpolate(x * y) # We create a 1D mesh immersed 2D from (0, 0) to (1, 1) which we call "line". # We want to calculate the line integral of f along it. To do this we # create a function space on the line mesh... V_line = FunctionSpace(line, "CG", 2) # ... and interpolate our function f onto it. f_line = interpolate(f, V_line) # The integral of f along the line is then a simple form expression which # we assemble: assemble(f_line * dx) # this outputs sqrt(2) / 3
For more on forms, see this section of the manual.
Interpolation from high-order meshes is currently not supported.
If the target mesh extends outside the source mesh domain, then cross-mesh
interpolation will raise a
# These meshes only share some of their domain src_mesh = UnitSquareMesh(2, 2) dest_mesh = UnitSquareMesh(3, 3, quadrilateral=True) dest_mesh.coordinates.dat.data_wo[:] *= 2 # We consider a simple function on our source mesh... x_src, y_src = SpatialCoordinate(src_mesh) V_src = FunctionSpace(src_mesh, "CG", 2) f_src = Function(V_src).interpolate(x_src**2 + y_src**2) # ... and want to interpolate into a function space on our target mesh ... V_dest = FunctionSpace(dest_mesh, "Q", 2) # ... but get a DofNotDefinedError if we try f_dest = interpolate(f_src, V_dest) # raises DofNotDefinedError
This can be overriden with the optional
f_dest = interpolate(f_src, V_dest, allow_missing_dofs=True) f_dest.at(0.5, 0.5) # returns 0.5**2 + 0.5**2 = 0.5
In this case, the missing degrees of freedom (DoFs, the global basis function coefficients which could not be set) are, by default, set to zero:
f_dest.at(1.5, 1.5) # returns 0.0
If we specify an output
Function then the missing DoFs are
We can optionally specify a value to use for our missing DoFs. Here
we set them to be
nan (‘not a number’) for easy identification:
f_dest = interpolate( f_src, V_dest, allow_missing_dofs=True, default_missing_val=np.nan ) f_dest.at(1.5, 1.5) # returns np.nan
If we specify an output
Function, this overwrites the missing
argument is set at construction:
interpolator = Interpolator(f_src, V_dest, allow_missing_dofs=True)
default_missing_val keyword argument is then set whenever we call
f_dest = interpolator.interpolate(default_missing_val=np.nan)
If we supply an output
Function and don’t set
default_missing_val then any missing DoFs are left as they were prior to
x_dest, y_dest = SpatialCoordinate(dest_mesh) f_dest = Function(V_dest).interpolate(x_dest + y_dest) f_dest.at(0.5, 0.5) # returns x_dest + y_dest = 1.0 interpolator.interpolate(output=f_dest) f_dest.at(0.5, 0.5) # now returns x_src^2 + y_src^2 = 0.5 # Similarly, using the interpolate method on a function will not overwrite # the pre-existing values if default_missing_val is not set f_dest.interpolate(f_src, allow_missing_dofs=True)
Unfortunately, UFL interpolation is not applicable if some of the
source data is not yet available as a Firedrake
or UFL expression. Here we describe a recipe for moving external to
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
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
proceeds as follows:
# First, grab the mesh. m = V.mesh() # 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
For interaction with external point data, see the corresponding manual section.
C string expressions were a FEniCS compatibility feature which has now been removed. Users should use UFL expressions instead. This section only remains to assist in the transition of existing code.
Here are a couple of old-style C string expressions, and their modern replacements.
# Expression: f = interpolate(Expression("sin(x*pi)"), V) # UFL equivalent: x = SpatialCoordinate(V.mesh()) f = interpolate(sin(x * math.pi), V) # Expression with a Constant parameter: f = interpolate(Expression('sin(x*t)', t=t), V) # UFL equivalent: x = SpatialCoordinate(V.mesh()) f = interpolate(sin(x * t), V)
Python expression classes were a FEniCS compatibility feature which has now been removed. Users should use UFL expressions instead. This section only remains to assist in the transition of existing code.
Expression classes expressions are
deprecated, below are a few examples on how to replace them with UFL
# 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))
randomfunctiongen module wraps the external numpy package numpy.random,
which gives Firedrake users an easy access to many stochastically sound random number generators,
SFC64, which are parallel-safe.
All distribution methods defined in numpy.random,
are made available, and one can pass a
FunctionSpace to most of these methods
to generate a randomised
mesh = UnitSquareMesh(2,2) V = FunctionSpace(mesh, "CG", 1) # PCG64 random number generator pcg = PCG64(seed=123456789) rg = RandomGenerator(pcg) # beta distribution f_beta = rg.beta(V, 1.0, 2.0) print(f_beta.dat.data) # produces: # [0.56462514 0.11585311 0.01247943 0.398984 0.19097059 0.5446709 0.1078666 0.2178807 0.64848515]