From b01a5239e838f8c9f7a8b1add9c729a1f7837ace Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 16 Jul 2026 21:36:39 +0100 Subject: [PATCH] Interpolation into MixedElement --- finat/mixed.py | 15 +++++++++++++++ test/finat/test_dual_basis.py | 24 ++++++++++++++++++++++++ 2 files changed, 39 insertions(+) diff --git a/finat/mixed.py b/finat/mixed.py index 1fef33055..d47effea1 100644 --- a/finat/mixed.py +++ b/finat/mixed.py @@ -95,6 +95,21 @@ def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=Non core_eval = self.element.point_evaluation(order, refcoords, entity) return self._transform_evaluation(core_eval) + def dual_evaluation(self, fn, coordinate_mapping=None): + """Evaluate the dual basis against this mixed component.""" + def component(points): + value = fn(points) + values = numpy.empty(self.element.value_shape, dtype=object) + for j, zeta in enumerate(numpy.ndindex(self.element.value_shape)): + values[zeta] = gem.Indexed(value, (self.offset + j,)) + if self.element.value_shape: + return gem.ListTensor(values) + return values.item() + + return self.element.dual_evaluation( + component, coordinate_mapping=coordinate_mapping + ) + @property def mapping(self): return self.element.mapping diff --git a/test/finat/test_dual_basis.py b/test/finat/test_dual_basis.py index 1737730a4..fbb1ee460 100644 --- a/test/finat/test_dual_basis.py +++ b/test/finat/test_dual_basis.py @@ -298,3 +298,27 @@ def summands(e): assert set(Q.free_indices) <= set(x.indices) assert len(element.dual_basis[1].points) \ == sum(len(e.dual_basis[1].points) for e in elements) + + +def test_mixed_subelement_dual_evaluation(): + # A mixed element is the direct sum of its components, so each component + # dual evaluates against its own slice of the flattened value vector. + cell = ufc_simplex(2) + scalar = finat.Lagrange(cell, 1) + vector = finat.TensorFiniteElement(scalar, (2,)) + mixed = finat.MixedElement([scalar, vector]) + vector_component = mixed.elements[1] + + def fn(points): + coordinates = points.expression + return gem.ListTensor([ + gem.Literal(0.0), + gem.Indexed(coordinates, (0,)), + gem.Indexed(coordinates, (1,)), + ]) + + expression, point_indices, basis_indices = vector_component.dual_evaluation(fn) + result, = evaluate([gem.ComponentTensor(gem.IndexSum(expression, point_indices), + basis_indices)]) + + assert numpy.allclose(result.arr, cell.vertices)