TSFC: support same-mesh ufl.Interpolate in Form kernels by pbrubeck · Pull Request #5258 · firedrakeproject/firedrake

pbrubeck · 2026-07-16T08:51:16Z

Description

Support same-mesh Interpolate nodes in normal TSFC Form kernels.
Fuse standalone Interpolate and Form(Interpolate) codegen.
Corollary: native mat_type="aij" support for same-mesh Interpolation across MixedFunctionSpaces
Corollary: interpolation into VertexOnlyMesh also goes through the unified TSFC codegen.
Kept cross-mesh interpolation on BaseFormAssembler.

Depends on
FEniCS/ufl#497
firedrakeproject/fiat#261

pbrubeck · 2026-07-16T09:42:38Z

    value_dict = {}
    for alpha, table in per_derivative.items():
        table_qi = gem.Indexed(table, beta + zeta)
+        if not hasattr(vec_beta, "index_ordering"):


hasattr is an antipattern, check isinstance instead

pbrubeck · 2026-07-16T09:50:47Z

 preserve_geometry_types = (CellVolume, FacetArea)


+class InterpolateMapper(MultiFunction):


DAGTraverser

pbrubeck · 2026-07-16T11:21:11Z

This PR should reimplement SameMeshInterpolator as a dispatcher to {Zero|One|Two}FormAssembler. We want to unify the TSFC codepaths for Form assembly and same-mesh interpolation.

pbrubeck · 2026-07-16T20:48:13Z

            (k, subdomain_id) for k in kernels for subdomain_id in k.kinfo.subdomain_id
        )

+    def _compile_interpolate(self):


This could be fused into the rest of the form assembly

pbrubeck · 2026-07-16T20:53:59Z

        if len(self._form.arguments()) == 2 and not self._diagonal:
            if not self._bcs:
                return None
-


put this back

pbrubeck · 2026-07-17T09:29:15Z

+
        summands = []
+        argument_multiindices = ctx.argument_multiindices
+        if hasattr(argument_multiindices, "values"):


Is this to check if we have a dict or None?

Ah, right we have mismatching conventions: Form uses a tuple and Inteprolate uses a dict

pbrubeck · 2026-07-17T14:49:43Z

-        integrand_exprs = builder.compile_integrand(integral.integrand(), params, ctx)
-        integral_exprs = builder.construct_integrals(integrand_exprs, params)
-        builder.stash_integrals(integral_exprs, params, ctx)
+        params["mode"] = "vanilla"


Is this hiding a bug?

pbrubeck · 2026-07-17T15:26:58Z

+        from firedrake.assemble import (
+            OneFormAssembler, TwoFormAssembler, ZeroFormAssembler,
+        )


TODO: let get_assembler do the dispatch

pbrubeck · 2026-07-17T15:30:51Z

+                for mesh, subelement in zip(
+                    domain.iterable_like(element), element.sub_elements
+                )
+            ]


Some of this could live in compute_form_data

pbrubeck · 2026-07-22T17:09:01Z

+        if (
+            isinstance(tensor, Function | Cofunction)
+            and any(set(tensor.dat).intersection(set(input_.dat))
+                    for input_ in inputs)
+        ):
+            output = tensor
+            tensor = assembler.allocate()
+            copyout = (partial(tensor.dat.copy, output.dat),)


I think this case is already handled by OneFormAssembler

pbrubeck · 2026-07-22T17:10:48Z

+            for op in ufl.algorithms.extract_base_form_operators(
+                integral.integrand()
+            ):
+                if (
+                    not isinstance(op, ufl.Interpolate)
+                    or not set(extract_domains(op)) <= valid_domains
+                ):


change to all()

pbrubeck · 2026-08-10T13:57:14Z

@@ -0,0 +1,172 @@
+Plate Bending Using Mixed Interpolation of Tensorial Components (MITC)


Make this demo more general by including the modified Morely method from https://www.jstor.org/stable/43694071

In this way we illustrate two compositions interpolate(grad) vs grad(interpolate)

pbrubeck · 2026-08-24T16:02:25Z

+  Pi_beta = interpolate(beta, R)
+  Pi_theta = interpolate(theta, R)
+
+  a_shear = G_shear * inner(grad(w) - Pi_beta, grad(v) - Pi_theta) * dx


This is the main feature exposed by this PR

Route a same-mesh Interpolate through the ordinary Form assembly path, so one code generator serves both: the kernel builder compiles an Interpolate the way it compiles an integrand, and the assembler allocates and applies boundary conditions for it the way it does for a Form.

function_space_match folded two questions into one predicate: whether a boundary condition belongs to a form's space, and which block of a mixed space it belongs to. Its callers then had to pre-index the space they passed, so the block comparison happened in two places at once and the predicate had to reconstruct what the caller had already discarded. Give function_space a `parent` argument, which walks up through indexed and component subspaces to the top-level space, and let each caller ask the two questions separately: the top-level spaces must be equal, and for a mixed space the block index must match the row or column. This surfaces the TypeError that ExplicitMatrixAssembler already raises for a boundary condition on an unrelated space, so the two tests that assert on it now name that type. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

prepare_arguments looks an argument's domain up in the integral type map. An argument that the form compiler has not split, which is what compile_form(split=False) hands it for a mixed space, carries a MeshSequence rather than one of the meshes the map is keyed on, so the lookup raises KeyError. PatchPC builds its patch operator that way, so every mixed-space Vanka or star patch failed there. Fall back to the meshes the sequence holds. They are all iterated by the same integral, so a single integral type has to come back, and the tuple unpacking says so. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

A point cloud target has no element to evaluate a dual basis against, so the interpolation compiles against a quadrature element on the source cell whose points are only known at run time. That element reached apply_mapping but not the kernel builder, which rebuilt the target element from the dual argument and tabulated the source on the vertex cell's empty point set. Carry it on the interpolation data instead. The domain numbering was inconsistent either side of the kernel: the point cloud counts as a domain of the interpolation, so it takes a slot in the cell sizes it cannot fill, and extract_domains has to reach the argument slots to see it at all. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

An Interpolate takes its test function from the dual space, and a dual space compares unequal to its primal. Once same-mesh interpolation assembles through ExplicitMatrixAssembler, that made two comparisons in the boundary condition path answer the wrong question. Asking which argument a condition belongs to, a condition named on the target space matched neither, and assembly raised a TypeError instead of applying it. Asking whether the matrix is square, an interpolation from a space to itself looked rectangular, so its boundary rows were zeroed without ever taking the diagonal entry that keeps them solvable. Take the primal space of each argument once, where the arguments are read, and let both comparisons stand as they are. The matrix a mixed interpolation assembles is now the one Form assembly builds, so the tests and the manual describe that matrix: every block of a nest is allocated, the off-diagonal ones assembling to zero rather than staying NULL; a blocked space defaults its nest blocks to baij; an unmixed space collapses its one-block nest instead of refusing it; and a mixed space rejects baij with the ValueError a Form raises. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

A Real-space argument is replaced by the literal 1 before a form reaches TSFC, because the Real element has a single global degree of freedom whose basis function is constant. That reasoning applies to a primal argument the form integrates against. An Interpolate's Coargument instead names the space the result lands in, and TSFC dual-evaluates it, so replacing it with 1 hands Interpolate.__init__ an IntValue where a space belongs. Mangle before the empty check, as the components are split before they are mangled and an empty component carries no Real argument either way.

An interpolation whose output is the mesh coordinate field reads those same coordinates through the kernel's coordinate argument, so it has to compute into a tensor of its own and copy out. The coordinates never appear among the form's coefficients, so the check for that overlap missed them and the interpolation overwrote, cell by cell, the coordinates it was still reading. Building the octahedral sphere is the case that notices: it remaps its own coordinates twice, and the second remap collapsed cells it had already written. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>