Interpolation into MixedSubElement

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

test/finat/test_dual_basis.py

@@ -306,3 +306,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)