Interpolation from a DG0 space on a VertexOnlyMesh

[nbouziani]

Describe the bug

I would like to interpolate a function from a DG0 space on a VertexOnlyMesh, P0DG, to a CG1 space on another mesh, U. This is presently possible, in the context of point forcing, by using the adjoint interpolation of the Interpolator mapping from U to P0DG. In other words, we can have the following interpolation: P0DG* -> CG1* (see. https://gist.github.com/ReubenHill/b600e56e8530b79117fa4b12caf1d0ff). However, I would like to directly interpolate between the primal spaces, i.e. P0DG -> CG1, which doesn't work at the moment, see the below example:

I am not sure if there are valid reasons for which only the adjoint interpolation should be possible or if it is just a limitation of the code.

PS: As a workaround, I am currently using the adjoint interpolation and taking the l2 Riesz representer.

Steps to Reproduce

import numpy as np
from firedrake import *

# Define U
mesh = UnitSquareMesh(20, 20)
U = FunctionSpace(mesh, "CG", 1)

# Define P0DG
coords = np.asarray([[0.5, 0.25], [0.5, 0.75], [0.25, 0.5], [0.75, 0.5]])
vom = VertexOnlyMesh(mesh, coords)
P0DG = FunctionSpace(vom, "DG", 0)

# Interpolator from P0DG to U
I = Interpolator(TestFunction(P0DG), U)

# Interpolate
f_p0dg = Function(P0DG).assign(1.)
f_u = Function(U)
I.interpolate(f_p0dg, output=f_u)

Expected behavior

I should be able to populate f_u with the result of interpolating f_p0dg into U.

Error message

NotImplementedError                       Traceback (most recent call last)
Cell In[4], line 14
     11 P0DG = FunctionSpace(vom, "DG", 0)
     13 # Interpolator from P0DG to U
---> 14 I = Interpolator(TestFunction(P0DG), U)
     16 # Interpolate
     17 f_p0dg = Function(P0DG).assign(1.)

File ~/firedrake/src/pyadjoint/pyadjoint/tape.py:109, in no_annotations.<locals>.wrapper(*args, **kwargs)
    106 @wraps(function)
    107 def wrapper(*args, **kwargs):
    108     with stop_annotating():
--> 109         return function(*args, **kwargs)

File ~/firedrake/src/firedrake/firedrake/interpolation.py:412, in CrossMeshInterpolator.__init__(self, expr, V, subset, freeze_expr, access, bcs, allow_missing_dofs)
    410 dest_node_coords = f_dest_node_coords.dat.data_ro
    411 try:
--> 412     self.vom_dest_node_coords_in_src_mesh = firedrake.VertexOnlyMesh(
    413         src_mesh,
    414         dest_node_coords,
    415         redundant=False,
    416         missing_points_behaviour=missing_points_behaviour,
    417     )
    418 except VertexOnlyMeshMissingPointsError:
    419     raise DofNotDefinedError(src_mesh, dest_mesh)

File petsc4py/PETSc/Log.pyx:188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File ~/firedrake/src/firedrake/firedrake/mesh.py:2840, in VertexOnlyMesh(mesh, vertexcoords, missing_points_behaviour, tolerance, redundant, name)
   2838 if pdim != gdim:
   2839     raise ValueError(f"Mesh geometric dimension {gdim} must match point list dimension {pdim}")
-> 2840 swarm, original_swarm, n_missing_points = _pic_swarm_in_mesh(
   2841     mesh, vertexcoords, tolerance=tolerance, redundant=redundant, exclude_halos=False
   2842 )
   2844 missing_points_behaviour = MissingPointsBehaviour(missing_points_behaviour)
   2846 if missing_points_behaviour != MissingPointsBehaviour.IGNORE:

File ~/firedrake/src/firedrake/firedrake/mesh.py:3114, in _pic_swarm_in_mesh(parent_mesh, coords, fields, tolerance, redundant, exclude_halos)
   3102 tdim = parent_mesh.topological_dimension()
   3103 gdim = parent_mesh.geometric_dimension()
   3105 (
   3106     coords_local,
   3107     global_idxs_local,
   3108     reference_coords_local,
   3109     parent_cell_nums_local,
   3110     owned_ranks_local,
   3111     input_ranks_local,
   3112     input_coords_idxs_local,
   3113     missing_global_idxs,
-> 3114 ) = _parent_mesh_embedding(
   3115     parent_mesh,
   3116     coords,
   3117     tolerance,
   3118     redundant,
   3119     exclude_halos,
   3120     remove_missing_points=False,
   3121 )
   3122 visible_idxs = parent_cell_nums_local != -1
   3123 if parent_mesh.extruded:
   3124     # need to store the base parent cell number and the height to be able
   3125     # to map point coordinates back to the parent mesh

File ~/firedrake/src/firedrake/firedrake/mesh.py:3607, in _parent_mesh_embedding(parent_mesh, coords, tolerance, redundant, exclude_halos, remove_missing_points)
   3534 """Find the parent mesh cells containing the given coordinates.
   3535 
   3536 Parameters
   (...)
   3603 
   3604 """
   3606 if isinstance(parent_mesh.topology, VertexOnlyMeshTopology):
-> 3607     raise NotImplementedError(
   3608         "VertexOnlyMeshes don't have a working locate_cells_ref_coords_and_dists method"
   3609     )
   3611 import firedrake.functionspace as functionspace
   3612 import firedrake.constant as constant

NotImplementedError: VertexOnlyMeshes don't have a working locate_cells_ref_coords_and_dists method

Environment:

---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |4066827   |False     |
|PyOP2               |master                        |d017d594  |False     |
|fiat                |master                        |d571bcc   |False     |
|firedrake           |master                        |2ae2ce8c6 |False     |
|h5py                |firedrake                     |6cc4c912  |False     |
|libspatialindex     |master                        |4768bf3   |True      |
|libsupermesh        |master                        |dbe226b   |False     |
|loopy               |main                          |8158afdb  |False     |
|petsc               |firedrake                     |34704c14dc|False     |
|pyadjoint           |master                        |1711274   |False     |
|pytest-mpi          |main                          |a478bc8   |False     |
|tsfc                |master                        |5515b2a   |False     |
|ufl                 |master                        |3f337c83  |False     |
---------------------------------------------------------------------------

Additional Info
Add any other context about the problem here.

[dham]

What you are asking for makes no mathematical sense. Interpolating a Function on a point space into a FunctionSpace defined on a mesh means evaluating the dual functions of the FunctionSpace on the Function. This requires the Function to be defined at all the points that the dual functions evaluate. This won't be true (except by vanishing chance) for a point space.

Please ask a question about the system you are actually attempting to simulate.

[nbouziani]

I have an algorithm that predicts the values of a field $u$ at a set of points $(x_i)_{i}$. At the moment, I can use those predictions in Firedrake to define a DG0 function on a point cloud. However, I would like, in addition, to be able to evaluate a residual form using this Function. Because this residual may involve the gradient of the reconstructed $u$, a DG0 function will not be enough.

[dham]

Please write up the maths somewhere (e.g. in a discussion) and we can talk about the right way to do it.