## 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 `Function`

s.

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
# [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.

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 `Interpolator`

s and interpolation see the
interpolation section.

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

or

```
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
`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.

### Points outside the domain¶

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