Firedrake’s first API for evaluating functions at arbitrary points,
at(), is designed for simple interrogation of a function with
a few points.
# evaluate f at a 1-dimensional point
f.at(0.3)
# evaluate f at two 1-dimensional points, or at one 2-dimensional point
# (depending on f's geometric dimension)
f.at(0.2, 0.4)
# evaluate f at one 2-dimensional point
f.at([0.2, 0.4])
# evaluate f at two 2-dimensional point
f.at([0.2, 0.4], [1.2, 0.5])
# evaluate f at two 2-dimensional point (same as above)
f.at([[0.2, 0.4], [1.2, 0.5]])
While in these examples we have only shown lists, other iterables
such as tuples and numpy arrays are also accepted. The following
are equivalent:
f.at(0.2, 0.4)
f.at((0.2, 0.4))
f.at([0.2, 0.4])
f.at(numpy.array([0.2, 0.4]))
For a single point, the result is a numpy array, or a tuple of
numpy arrays in case of mixed functions. When evaluating
multiple points, the result is a list of values for each point.
To summarise:
Single point, non-mixed: numpy array
Single point, mixed: tuple of numpy arrays
Multiple points, non-mixed: list of numpy arrays
Multiple points, mixed: list of tuples of numpy arrays
When any point is outside the domain of the function,
PointNotInDomainError exception is raised. If
dont_raise=True is passed to at(), the result is
None for those points which fall outside the domain.
mesh = UnitIntervalMesh(8)
f = mesh.coordinates
f.at(1.2) # raises exception
f.at(1.2, dont_raise=True) # returns None
f.at(0.5, 1.2) # raises exception
f.at(0.5, 1.2, dont_raise=True) # returns [0.5, None]
If you move the mesh, by changing the mesh coordinates, then the bounding box tree that Firedrake
maintains to ensure fast point evaluation must be rebuilt. To do
this, after moving the mesh, call
clear_spatial_index() on the mesh you have just
moved.
There is limited support for at() when running Firedrake
in parallel. There is no special API, but there are some restrictions:
Point evaluation is a collective operation.
Each process must ask for the same list of points.
Each process will get the same values.
If RuntimeError: Point evaluation gave different results across processes.
is raised, try lowering the mesh tolerance.
Firedrake’s principal API for evaluating functions at arbitrary points,
interpolation onto a VertexOnlyMesh(), is designed for evaluating a
function at many points, or repeatedly, and for creating expressions which
contain point evaluations. It is parallel-safe. Whilst at()
produces a list of values, cross-mesh interpolation onto
VertexOnlyMesh() gives Firedrake Functions.
This is discussed in detail in [NHSCH23] but, briefly,
the idea is that the VertexOnlyMesh() is a mesh that represents a
point cloud domain. Each cell of the mesh is a vertex at a chosen location in
space. As usual for a mesh, we represent values by creating functions in
function spaces on it. The only function space that makes sense for a mesh
whose cells are vertices is the space of piecewise constant functions, also
known as the Polynomial degree 0 Discontinuous Galerkin (P0DG) space.
Our vertex-only meshes are immersed in some ‘parent’ mesh. We perform point
evaluation of a function \(f\) defined in a function space
\(V\) on the parent mesh by interpolating into the P0DG space on the
VertexOnlyMesh(). For example:
parent_mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(parent_mesh, "CG", 2)
# Create a function f on the parent mesh to point evaluate
x, y = SpatialCoordinate(parent_mesh)
f = Function(V).interpolate(x**2 + y**2)
# 3 points (i.e. vertices) at which to point evaluate f
points = [[0.1, 0.1], [0.2, 0.2], [0.3, 0.3]]
vom = VertexOnlyMesh(parent_mesh, points)
# P0DG is the only function space you can make on a vertex-only mesh
P0DG = FunctionSpace(vom, "DG", 0)
# Interpolation performs point evaluation
f_at_points = interpolate(f, P0DG)
print(f_at_points.dat.data_ro)
will print [0.02, 0.08, 0.18] when running in serial, the values of
\(x^2 + y^2\) at the points \((0.1, 0.1)\), \((0.2, 0.2)\) and
\((0.3, 0.3)\). For details on viewing the outputs in parallel, see the
section on the input ordering property.
Note that f_at_points is a Function which takes
on all the values of f evaluated at points. The cell ordering of a
VertexOnlyMesh() follows the ordering of the list of points it is given
at construction. In general VertexOnlyMesh() accepts any numpy array of
shape (num_points, point_dim) (or equivalent list) as the set of points to
create disconnected vertices at.
The operator for evaluation at the points specified can be
created by making an Interpolator acting on a
TestFunction()
u = TestFunction(V)
Interpolator(u, P0DG)
For more on Interpolators and interpolation see the
interpolation section.
When interpolating from vector or tensor valued function spaces, the P0DG
function space on the vertex-only mesh must be a
VectorFunctionSpace() or TensorFunctionSpace()
respectively. For example:
V = VectorFunctionSpace(parent_mesh, "CG", 2)
or
V = FunctionSpace(parent_mesh, "N1curl", 2)
each require
vom = VertexOnlyMesh(parent_mesh, points)
P0DG_vec = VectorFunctionSpace(vom, "DG", 0)
for successful interpolation.
In parallel the points given to VertexOnlyMesh() are assumed to be
the same on each MPI process and are taken from rank 0. To let different ranks
provide different points to the vertex-only mesh set the keyword argument
redundant = False
# Default behaviour
vom = VertexOnlyMesh(parent_mesh, points, redundant = True)
# Different points on each MPI rank to add to the vertex-only mesh
vom = VertexOnlyMesh(parent_mesh, points, redundant = False)
In this case, points will redistribute to the mesh partition where they are
located. This means that if rank A has points \(\{X\}\) that are not
found in the mesh cells owned by rank A but are found in the mesh cells owned
by rank B then they will be moved to rank B.
If the same coordinates are supplied more than once, they are always assumed to
be a new vertex: this is true for both redundant = True and
redunant = False. So if we have the same set of points on all MPI processes
and switch from redundant = True to redundant = False we will get point
duplication.
Be default points outside the domain by more than the specified
tolerance will generate
a VertexOnlyMeshMissingPointsError. This can be switched to a
warning or switched off entirely:
# point (1.1, 1.0) is outside the mesh
points = [[0.1, 0.1], [0.2, 0.2], [1.1, 1.0]]
# This will raise a VertexOnlyMeshMissingPointsError
vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour='error')
# This will generate a warning and the point will be lost
vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour='warn')
# This will cause the point to be silently lost
vom = VertexOnlyMesh(parent_mesh, points, missing_points_behaviour=None)
Integrating over a vertex-only mesh is equivalent to summing over
it. So if we have a vertex-only mesh \(\Omega_v\) with \(N\) vertices
at points \(\{x_i\}_{i=0}^{N-1}\) and we have interpolated a function
\(f\) onto it giving a new function \(f_v\) then
\[\int_{\Omega_v} f_v \, dx = \sum_{i=0}^{N-1} f(x_i).\]
These equivalent expressions for point evaluation
\[\sum_{i=0}^{N-1} f(x_i) = \sum_{i=0}^{N-1} \int_\Omega f(x) \delta(x - x_i) \, dx\]
where \(N\) is the number of points, \(x_i\) is the \(i\)th point,
\(\Omega\) is a ‘parent’ mesh, \(f\) is a function on that mesh,
\(\delta\) is a dirac delta distribition can therefore be written in
Firedrake using VertexOnlyMesh() and interpolate() as
omega = UnitSquareMesh(100, 100, quadrilateral=True)
V = FunctionSpace(omega, "CG", 2)
f = Function(V)
points = [[0.1, 0.1], [0.2, 0.2], [1.0, 1.0]]
# Create a vertex-only mesh at the points
vom = VertexOnlyMesh(omega, points)
# Create a P0DG function space on the vertex-only mesh
P0DG = FunctionSpace(vom, "DG", 0)
# Interpolating f into the P0DG space on the vertex-only mesh evaluates f at
# the points
expr = assemble(interpolate(f, P0DG)*dx)
Aside from input_ordering, we can use
interpolate() to interact with external data to, for example,
compare a PDE solution with the point data. The \(l_2\) error norm
(euclidean norm) of a function \(f\) (which may be a PDE solution)
evaluated against a set of point data \(\{y_i\}_{i=0}^{N-1}\) at points
\(\{x_i\}_{i=0}^{N-1}\) is defined as
\[\sqrt{ \sum_{i=0}^{N-1} (f(x_i) - y_i)^2 }.\]
We can express this in Firedrake as
error = sqrt(assemble((interpolate(f, P0DG) - y_pts)**2*dx))
# or equivalently
error = errornorm(interpolate(f, P0DG), y_pts)
We can then use the input_ordering vertex-only mesh
to safely check the values of error at the points
\(\{x_i\}_{i=0}^{N-1}\).
If points are outside the mesh domain but ought to still be found a
tolerance parameter can be set. The tolerance is relative to the size of
the mesh cells and is a property of the mesh itself
parent_mesh = UnitSquareMesh(100, 100, quadrilateral=True)
# point (1.1, 1.0) is outside the mesh
points = [[0.1, 0.1], [0.2, 0.2], [1.1, 1.0]]
# This prints 0.5 - points can be up to around half a mesh cell width away from
# the edge of the mesh and still be considered inside the domain.
print(parent_mesh.tolerance)
# This changes the tolerance and will cause the spatial index of the mesh
# to be rebuilt when first performing point evaluation which can take some
# time
parent_mesh.tolerance = 20.0
# This will now include the point (1.1, 1.0) in the mesh since each mesh
# cell is 1.0/100.0 wide.
vom = VertexOnlyMesh(parent_mesh, points)
# Similarly .at will not generate an error
V = FunctionSpace(parent_mesh, 'CG', 2)
Function(V).at((1.1, 1.0))
Alternatively we can modify the tolerance by providing it as an argument to the
vertex-only mesh. This will modify the tolerance property of the parent mesh.
Warning
To avoid confusion, it is recommended that the tolerance be set for a mesh
before any point evaluations are performed, rather than making use of these
keyword arguments.
# The point (1.1, 1.0) will still be included in the vertex-only mesh
vom = VertexOnlyMesh(parent_mesh, points, tolerance=30.0)
# The tolerance property has been changed - this will print 30.0
print(parent_mesh.tolerance)
# This doesn't generate an error
Function(V).at((1.1, 1.0), tolerance=20.0)
# The tolerance property has been changed again - this will print 20.0
print(parent_mesh.tolerance)
try:
# This generates an error
Function(V).at((1.1, 1.0), tolerance=1.0)
except PointNotInDomainError:
# But the tolerance property has still been changed - this will print 1.0
print(parent_mesh.tolerance)
Note that since our tolerance is relative, the number of cells in a mesh
dictates the point loss behaviour close to cell edges. So the mesh
UnitSquareMesh(5, 5, quadrilateral = True) will include the point
\((1.1, 1.0)\) by default.
Changing mesh tolerance only affects point location after it has been changed.
To apply the new tolerance to a vertex-only mesh, a new vertex-only mesh must
be created. Any existing immersed vertex-only mesh will have been created
using the previous tolerance and will be unaffected by the change.
UFL reserves the function call operator for evaluation:
will evaluate \(f\) at \((0.2, 0.4)\). UFL does not accept
multiple points at once, and cannot configure what to do with a point
which is not in the domain. The advantage of this syntax is that it
works on any Expr, for example:
will evaluate \(f \cdot \sin(f)\) at \((0.2, 0.4)\).
Note
The expression itself is not translated into C code. While the
evaluation of a function uses the same infrastructure as the
Firedrake APIs, which use generated C code, the expression tree is
evaluated by UFL in Python.
