Warning

You are reading a version of the website built against the unstable main branch. This content is liable to change without notice and may be inappropriate for your use case. You can find the documentation for the current stable release here.

Point evaluation ¶

Firedrake can evaluate Functions at arbitrary physical points. This feature can be useful for the evaluation of the result of a simulation, or for creating expressions which contain point evaluations. Three APIs are offered to this feature: two Firedrake-specific ones, and one from UFL.

Firedrake convenience function ¶

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 points
f.at([0.2, 0.4], [1.2, 0.5])

# evaluate f at two 2-dimensional points (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

Points outside the domain ¶

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]

Evaluation on a moving mesh ¶

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.

Evaluation with a distributed mesh ¶

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.

Primary API: Interpolation onto a vertex-only mesh ¶

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 [NHSCH24] 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
# [test_vertex_only_mesh_manual_example 2]
f_at_points = assemble(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.

Vector and tensor valued function spaces ¶

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)

V = FunctionSpace(parent_mesh, "N1curl", 2)

each require

vom = VertexOnlyMesh(parent_mesh, points)
P0DG_vec = VectorFunctionSpace(vom, "DG", 0)

for successful interpolation.

Parallel behaviour ¶

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 redundant = 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.

Points outside the domain ¶

By 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="ignore")

`PointEvaluator` convenience object ¶

The PointEvaluator class performs point evaluation using vertex-only meshes, as described above. This is a convenience object for users who want to evaluate functions at points without writing the vertex-only mesh boilerplate code each time.

First, create a PointEvaluator object by passing the parent mesh and the points to evaluate at:

point_evaluator = PointEvaluator(mesh, points)

Internally, this creates a vertex-only mesh at the given points, immersed in the given mesh. To evaluate a Function defined on the parent mesh at the given points, we use evaluate():

f_at_points = point_evaluator.evaluate(f)

Under the hood, this creates the appropriate P0DG function space on the vertex-only mesh and performs the interpolation. The points are then reordered to match the input ordering, as described in the section on the input ordering property. The result is a Numpy array containing the values of f at the given points, in the order the points were provided to the PointEvaluator constructor.

If redundant=True (the default) was used when creating the PointEvaluator, then only the points given to the constructor on rank 0 will be evaluated. The result is then broadcast to all ranks. Use this option if the same points are given on all ranks. If redundant=False was used when creating the PointEvaluator, then each rank will evaluate the points it was given. Use this option if different points are given on different ranks, for example when using external point data.

The parameters missing_points_behaviour and tolerance (discussed here and here respectively) can be set when creating the PointEvaluator and will be passed to the VertexOnlyMesh() it creates internally.

If the coordinates or the tolerance of the parent mesh are changed after creating the PointEvaluator, then the vertex-only mesh will be reconstructed on the new mesh the next time evaluate() is called.

Expressions with point evaluations ¶

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 distribution 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
f_at_points = assemble(interpolate(f, P0DG))
expr = assemble(f_at_points * dx)

Interacting with external point data ¶

Using the input ordering property ¶

Any set of points with associated data in our domain can be expressed as a P0DG function on a VertexOnlyMesh(). The recommended way to import data from an external source is via the input_ordering property: this produces another vertex-only mesh which has points in the order and MPI rank that they were specified when first creating the original vertex-only mesh. For example:

# We have a set of points with corresponding data from elsewhere which vary
# from rank to rank
vom = VertexOnlyMesh(parent_mesh, point_locations_from_elsewhere, redundant=False)
P0DG = FunctionSpace(vom, "DG", 0)

# Create a P0DG function on the input ordering vertex-only mesh
P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0)
point_data_input_ordering = Function(P0DG_input_ordering)

# We can safely set the values of this function, knowing that the data will
# be in the same order and on the same MPI rank as point_locations_from_elsewhere
point_data_input_ordering.dat.data_wo[:] = point_data_values_from_elsewhere

# Interpolate puts this data onto the original vertex-only mesh
point_data = assemble(interpolate(point_data_input_ordering, P0DG))

This is entirely parallel safe.

Similarly, we can use input_ordering to get data out of a vertex-only mesh in a parallel-safe way. If we return to our example from the section where we introduced vertex only meshes, we had

f_at_points = assemble(interpolate(f, P0DG))

print(f_at_points.dat.data_ro)

In parallel, this will print the values of f at the given points list after the points have been distributed over the parent mesh. If we want the values of f at the points list before the points have been distributed we can use input_ordering as follows:

# Make a P0DG function on the input ordering vertex-only mesh again
P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0)
f_at_input_points = Function(P0DG_input_ordering)

# We interpolate the other way this time
f_at_input_points.interpolate(f_at_points)

print(f_at_input_points.dat.data_ro)  # will print the values at the input points

Note

When a vertex-only mesh is created with redundant = True (which is the default when creating a VertexOnlyMesh()) the input_ordering method will return a vertex-only mesh with all points on rank 0.

If we ran the example in parallel, the above code would print [0.02, 0.08, 0.18] on rank 0 and [] on all other ranks. If we set redundant = False when creating the vertex-only mesh, the above code would print [0.02, 0.08, 0.18] on all ranks and we would have point duplication.

If any of the specified points were not found in the mesh, the value on the input ordering vertex-only mesh will not be affected by the interpolation from the original vertex-only mesh. In the above example, the values would be zero at those points. To make it more obvious that those points were not found, it’s a good idea to set the values to nan before the interpolation:

import numpy as np
f_at_input_points.dat.data_wo[:] = np.nan
f_at_input_points.interpolate(f_at_points)

print(f_at_input_points.dat.data_ro)  # any points not found will be NaN

More ways to interact with external data ¶

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}\).

Mesh tolerance ¶

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))

Keyword arguments ¶

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 API ¶

UFL reserves the function call operator for evaluation:

f([0.2, 0.4])

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:

(f*sin(f))([0.2, 0.4])

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.