Clement interpolant #4244

jwallwork23 · 2025-04-22T07:48:00Z

I've been meaning to upstream some functionality from Animate that might be useful for the wider Firedrake community, so am starting with this. Let me know if it's of interest!

This PR adds functionality for the Clément interpolant, which transfers fields from P1 space to P0 space via local volume averaging. We've found this particularly useful for gradient/Hessian recovery in P1 space (can upstream those later).

I've included tests on simplices in 1D, 2D, and 3D, and for scalar-, vector-, and matrix-valued functions.

A limitation of the approach is that the local averaging breaks down on the boundary, but it does a good job in the domain interior. I also have an implementation of local facet averaging on the boundary that I can provide, if it'd be of use.

pbrubeck · 2025-04-22T07:50:44Z

firedrake/interpolation.py

+    For arguments, see :class:`.Interpolator`.
+    """
+
+    def __new__(cls, expr, V, **kwargs):


Should this check go to __init__ ?

The reason I put it here is that the class inherits from Interpolator, which defines __new__ to account for the different implementations in SameMeshInterpolator and CrossMeshInterpolator. I could overwrite __new__ to just pass and then put this logic into __init__ if that'd be preferred - just let me know.

__new__ should only be overloaded (and very rarely!) to dispatch to different class constructors. It isn't needed here so I wouldn't define a __new__ at all.

Well I figured it'd make sense for ClementInterpolator to inherit from Interpolator, so something needs to be done to account for its implementation there, right? Or did I miss something?

Could you not just leave __new__ unimplemented, and move the check to __init__ ?

If I do that I get

# Projecting into P0 space and then applying Clement interpolation should recover # the original function x_P0 = assemble(project(expr, P0)) interpolator = ClementInterpolator(TestFunction(P0), P1) > x_P1 = interpolator.interpolate(x_P0) tests/firedrake/regression/test_interpolate.py:604: _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func ??? petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func ??? firedrake/interpolation.py:418: in interpolate interp = self._interpolate_future(*function, adjoint=adjoint, _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ self = <firedrake.interpolation.SameMeshInterpolator object at 0x72d277a4b7d0> transpose = None, adjoint = False, default_missing_val = None function = (Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x72d277bb3440>, TensorElement(FiniteEl...ape=(3, 3), symmetry={}), name=None), Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 45)), 69),) def _interpolate_future(self, *function, transpose=None, adjoint=False, default_missing_val=None): """Define the :class:`Interpolate` object corresponding to the interpolation operation of interest. Parameters ---------- *function: firedrake.function.Function or firedrake.cofunction.Cofunction If the expression being interpolated contains an argument, then the function value to interpolate. transpose : bool Deprecated, use adjoint instead. adjoint: bool Set to true to apply the adjoint of the interpolation operator. default_missing_val: bool For interpolation across meshes: the optional value to assign to DoFs in the target mesh that are outside the source mesh. If this is not set then the values are either (a) unchanged if some ``output`` is specified to the :meth:`interpolate` method or (b) set to zero. This does not affect adjoint interpolation. Ignored if interpolating within the same mesh or onto a :func:`.VertexOnlyMesh`. Returns ------- firedrake.interpolation.Interpolate or ufl.action.Action or ufl.adjoint.Adjoint The symbolic object representing the interpolation operation. Notes ----- This method is the default future behaviour of interpolation. In a future release, the ``Interpolator.interpolate`` method will be replaced by this method. """ > V = self.V E AttributeError: 'SameMeshInterpolator' object has no attribute 'V' firedrake/interpolation.py:336: AttributeError

Oh god. I bet that type(ClementInterpolator(TestFunction(P0), P1)) returns SameMeshInterpolator because of the __new__ magic in Interpolator.

I think that this level of the code really needs cleaning up, but for now I think you can just have a minimal __new__:

def __new__(cls, *args, **kwargs): return object.__new__(ClementInterpolator)

with a comment explaining why it is necessary.

Thanks - done in db99911.

pbrubeck · 2025-04-22T08:46:13Z

firedrake/interpolation.py

+        if rank != len(self.V.value_shape):
+            raise ValueError(f"Rank-{rank} input inconsistent with target space.")
+        mesh = self.V.mesh()
+        dim = mesh.topological_dimension()


This dim should be the number of copies of the element V.block_size. I think a VectorFunctionSpace and TensorFunctionSpace could be indexed with a double index (flat_dof_id, component_id). So rank could possibly be unrestricted and there wouldn't have to be special cases if you treat the rank-0 as the general case.

Addressed in 7ec5197.

pbrubeck · 2025-04-22T08:47:10Z

tests/firedrake/regression/test_interpolate.py

+
+@pytest.mark.parametrize("rank", range(3))
+@pytest.mark.parametrize("dim", range(1, 4))
+def test_clement_interpolator_simplex(dim, rank):


We should test a vector with 3 components in a 2D mesh.

Good suggestion - refactoring for this helped to address some limitations in the implementation. Addressed in 7ec5197.

pbrubeck · 2025-04-22T08:53:05Z

firedrake/interpolation.py

+        domain = "{[i]: 0 <= i < patch.dofs}"
+        instructions = "patch[i] = patch[i] + vol[0]"
+        keys = {"vol": (volume, op2.READ), "patch": (patch_volume, op2.RW)}
+        parloops.par_loop((domain, instructions), dX, keys)


The patch volume here could become a cached_property so it is not recomputed every time. Maybe we could attach this property to the mesh, so it can be reused more widely across different Interpolators.

Good idea. I've addressed this in 34ceb67, using the cell_sizes cached property as a template.

jwallwork23 · 2025-04-28T09:16:25Z

I added Clement1975 to docs/source/_static/references.bib but notice the docs build gives

/__w/firedrake/firedrake/venv/lib/python3.12/site-packages/firedrake/interpolation.py:docstring of firedrake.interpolation.ClementInterpolator:4: WARNING: could not find bibtex key "Clement1975" [bibtex.key_not_found]

Is there somewhere else that I should be putting this reference?

firedrake/interpolation.py

firedrake/mesh.py

Co-authored-by: Pablo Brubeck <[email protected]>

connorjward · 2025-04-28T10:48:26Z

firedrake/interpolation.py

+        if subset:
+            raise NotImplementedError("subset not implemented")
+        if freeze_expr:
+            raise NotImplementedError("freeze_expr not implemented")
+        if access != op2.WRITE:
+            raise NotImplementedError("access other than op2.WRITE not implemented")
+        if bcs:
+            raise NotImplementedError("bcs not implemented")


This sort of thing makes me wonder whether this is the best place to implement this functionality. Might it be better to simply have a free function called interpolate_clement that does the parloop?

As in instead of having the Interpolator at all?

Yes. I think we need @dham to weigh in here. Is this the right abstraction?

This doesn't look great. Interpolator objects should be dying out of the public API of Firedrake. They are really just the implementation backend to interpolate. Having users directly instantiate Interpolator objects will, for example, break differentiability.

I think the actually correct way to do this is to inherit from AbstractExternalOperator so that you get all the right symbolic properties for free. The syntax would then be assemble(interpolate_clement(source, target)). If you haven't yet implemented the operator (i.e. the matrix) or the adjoint then these will raise NotImplementedError and at least the right interface will be there for the future.

It might be nice to have assemble(interpolate(source, target, variant="clement")) but that might be tricky/unwieldy.

I think we'd have to think harder about whether that is actually the Right Thing. For a finite element space with a known dual basis, there is a single canonical interpolation operation. Other interpolation operations are meaningfully different: they are just other projections that don't directly arise from the definition of the space.

The variant option at least requires heavier surgery because interpolate is defined in UFL.

pbrubeck · 2025-04-28T10:56:09Z

tests/firedrake/regression/test_interpolate.py

+
+@pytest.mark.parametrize("tdim,shape", [(1, tuple()), (2, tuple()), (3, tuple()),
+                                        (1, (1,)), (2, (2,)), (2, (3,)), (3, (3,)),
+                                        (1, (1, 1)), (2, (2, 2)), (3, (3, 3))],


Suggested change

(1, (1, 1)), (2, (2, 2)), (3, (3, 3))],

(1, (1, 2)), (2, (2, 3)), (3, (2, 3))],

We would have to create expr with the right shape. Maybe something like:

P1 = fs_constructor(mesh, "CG", 1) x = SpatialCoordinate(mesh) expr = as_tensor(np.reshape([dot(Constant(tuple(range(i+1, i+1+tdim))), x) for i in range(P1.block_size)], shape))

Addressed in cd20a5e and simplified the notation while I was at it.

jwallwork23 · 2025-05-02T08:01:06Z

The most recent test fails don't seem to be related to the changes in this PR.

dham · 2025-05-03T10:25:02Z

firedrake/mesh.py

+        from firedrake.functionspace import FunctionSpace
+        DG0 = FunctionSpace(self, "Discontinuous Lagrange", 0)
+        volume = Function(DG0)
+        return volume.interpolate(ufl.CellVolume(self))


Is this actually mathematically correct for higher order meshes? Shouldn't it be a projection?

As long as CellVolume is correctly integrated then this should be fine. Recall that CellVolume is intrinsically a DG0 function on any mesh.

jwallwork23 added 3 commits April 21, 2025 14:04

Implement Clement interpolant

560bad2

Test Clement interpolant on simplex meshes

88d2748

Use parloop instead of interpolate call; clearer var names

b4dfbd9

jwallwork23 added the enhancement label Apr 22, 2025

jwallwork23 self-assigned this Apr 22, 2025

pbrubeck reviewed Apr 22, 2025

View reviewed changes

jwallwork23 added 3 commits April 22, 2025 18:12

Add reference for Clement1975

651b01f

Introduce cached properties for cell_volume and patch_volume

34ceb67

Extend Clement to account for vector and tensor sizes other than tdim

7ec5197

jwallwork23 requested a review from pbrubeck April 28, 2025 09:05

Lint

6e67568

pbrubeck requested changes Apr 28, 2025

View reviewed changes

firedrake/interpolation.py Outdated Show resolved Hide resolved

firedrake/interpolation.py Outdated Show resolved Hide resolved

firedrake/interpolation.py Outdated Show resolved Hide resolved

firedrake/interpolation.py Outdated Show resolved Hide resolved

pbrubeck reviewed Apr 28, 2025

View reviewed changes

firedrake/interpolation.py Outdated Show resolved Hide resolved

pbrubeck reviewed Apr 28, 2025

View reviewed changes

firedrake/mesh.py Outdated Show resolved Hide resolved

pbrubeck reviewed Apr 28, 2025

View reviewed changes

firedrake/mesh.py Outdated Show resolved Hide resolved

jwallwork23 and others added 2 commits April 28, 2025 11:05

Apply suggestions from @pbrubeck code review

63b5cb7

Co-authored-by: Pablo Brubeck <[email protected]>

Lint again

9323bcf

jwallwork23 requested review from pbrubeck and connorjward April 28, 2025 10:22

connorjward reviewed Apr 28, 2025

View reviewed changes

pbrubeck reviewed Apr 28, 2025

View reviewed changes

jwallwork23 added 2 commits May 2, 2025 08:03

Simplify __new__

db99911

Extend tensor cases; simplify notation

cd20a5e

jwallwork23 requested review from pbrubeck and connorjward May 2, 2025 07:18

dham reviewed May 3, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Clement interpolant #4244

Clement interpolant #4244

jwallwork23 commented Apr 22, 2025 •

edited

Loading

pbrubeck Apr 22, 2025

jwallwork23 Apr 28, 2025

connorjward Apr 28, 2025

jwallwork23 Apr 28, 2025

pbrubeck Apr 28, 2025

jwallwork23 Apr 28, 2025

connorjward Apr 28, 2025

connorjward Apr 28, 2025

jwallwork23 May 2, 2025

pbrubeck Apr 22, 2025

jwallwork23 Apr 28, 2025

pbrubeck Apr 22, 2025 •

edited

Loading

jwallwork23 Apr 28, 2025

pbrubeck Apr 22, 2025

jwallwork23 Apr 25, 2025

jwallwork23 commented Apr 28, 2025

connorjward Apr 28, 2025

jwallwork23 May 2, 2025

connorjward May 2, 2025

dham May 2, 2025

connorjward May 2, 2025

dham May 2, 2025 •

edited

Loading

dham May 2, 2025

pbrubeck Apr 28, 2025

pbrubeck Apr 28, 2025 •

edited

Loading

jwallwork23 May 2, 2025

jwallwork23 commented May 2, 2025

dham May 3, 2025

pbrubeck May 4, 2025

	(1, (1, 1)), (2, (2, 2)), (3, (3, 3))],
	(1, (1, 2)), (2, (2, 3)), (3, (2, 3))],

Clement interpolant #4244

Are you sure you want to change the base?

Clement interpolant #4244

Conversation

jwallwork23 commented Apr 22, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pbrubeck Apr 22, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jwallwork23 commented Apr 28, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dham May 2, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pbrubeck Apr 28, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jwallwork23 commented May 2, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jwallwork23 commented Apr 22, 2025 •

edited

Loading

pbrubeck Apr 22, 2025 •

edited

Loading

dham May 2, 2025 •

edited

Loading

pbrubeck Apr 28, 2025 •

edited

Loading