PhysicallyMapped: implement point_evaluation (#144)

pbrubeck · web-flow · commit ce61877bcb94 · 2025-04-10T11:30:13.000+01:00
* PhysicallyMapped: implement point_evaluation
diff --git a/finat/cube.py b/finat/cube.py
@@ -77,11 +77,12 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
 
         return self.product.basis_evaluation(order, ps, self._unflatten[entity])
 
-    def point_evaluation(self, order, point, entity=None):
+    def point_evaluation(self, order, point, entity=None, coordinate_mapping=None):
         if entity is None:
             entity = (self.cell.get_spatial_dimension(), 0)
 
-        return self.product.point_evaluation(order, point, self._unflatten[entity])
+        return self.product.point_evaluation(order, point, self._unflatten[entity],
+                                             coordinate_mapping)
 
     @property
     def dual_basis(self):
diff --git a/finat/direct_serendipity.py b/finat/direct_serendipity.py
@@ -122,7 +122,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
 
         return result
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, point, entity=None, coordinate_mapping=None):
         raise NotImplementedError("Not done yet, sorry!")
 
     def mapping(self):
diff --git a/finat/discontinuous.py b/finat/discontinuous.py
@@ -68,8 +68,8 @@ def fiat_equivalent(self):
     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         return self.element.basis_evaluation(order, ps, entity, coordinate_mapping=coordinate_mapping)
 
-    def point_evaluation(self, order, refcoords, entity=None):
-        return self.element.point_evaluation(order, refcoords, entity)
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
+        return self.element.point_evaluation(order, refcoords, entity, coordinate_mapping)
 
     @property
     def dual_basis(self):
diff --git a/finat/enriched.py b/finat/enriched.py
@@ -125,7 +125,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
                    for element in self.elements]
         return self._compose_evaluations(results)
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         '''Return code for evaluating the element at an arbitrary points on
         the reference element.
 
@@ -136,7 +136,7 @@ def point_evaluation(self, order, refcoords, entity=None):
                           free indices are arbitrary.
         :param entity: the cell entity on which to tabulate.
         '''
-        results = [element.point_evaluation(order, refcoords, entity)
+        results = [element.point_evaluation(order, refcoords, entity, coordinate_mapping)
                    for element in self.elements]
         return self._compose_evaluations(results)
 
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
@@ -156,7 +156,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
                 result[alpha] = expr
         return result
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         '''Return code for evaluating the element at an arbitrary points on
         the reference element.
 
@@ -166,6 +166,9 @@ def point_evaluation(self, order, refcoords, entity=None):
                           a vector with the correct dimension, its
                           free indices are arbitrary.
         :param entity: the cell entity on which to tabulate.
+        :param coordinate_mapping: a
+           :class:`~.physically_mapped.PhysicalGeometry` object that
+           provides physical geometry callbacks (may be None).
         '''
         if entity is None:
             entity = (self.cell.get_dimension(), 0)
diff --git a/finat/finiteelementbase.py b/finat/finiteelementbase.py
@@ -165,7 +165,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         '''
 
     @abstractmethod
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         '''Return code for evaluating the element at an arbitrary points on
         the reference element.
 
@@ -175,6 +175,9 @@ def point_evaluation(self, order, refcoords, entity=None):
                           a vector with the correct dimension, its
                           free indices are arbitrary.
         :param entity: the cell entity on which to tabulate.
+        :param coordinate_mapping: a
+           :class:`~.physically_mapped.PhysicalGeometry` object that
+           provides physical geometry callbacks (may be None).
         '''
 
     @property
diff --git a/finat/hdivcurl.py b/finat/hdivcurl.py
@@ -74,7 +74,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         core_eval = self.wrappee.basis_evaluation(order, ps, entity)
         return self._transform_evaluation(core_eval)
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         core_eval = self.wrappee.point_evaluation(order, refcoords, entity)
         return self._transform_evaluation(core_eval)
 
diff --git a/finat/mixed.py b/finat/mixed.py
@@ -91,7 +91,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         core_eval = self.element.basis_evaluation(order, ps, entity, coordinate_mapping=coordinate_mapping)
         return self._transform_evaluation(core_eval)
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         core_eval = self.element.point_evaluation(order, refcoords, entity)
         return self._transform_evaluation(core_eval)
 
diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py
@@ -264,7 +264,7 @@ def basis_transformation(self, coordinate_mapping):
         :arg coordinate_mapping: Object providing physical geometry."""
         pass
 
-    def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
+    def map_tabulation(self, tabulation, coordinate_mapping):
         assert coordinate_mapping is not None
 
         M = self.basis_transformation(coordinate_mapping)
@@ -281,13 +281,15 @@ def matvec(table):
             val = gem.ListTensor(expressions)
             return gem.optimise.aggressive_unroll(val)
 
-        result = super().basis_evaluation(order, ps, entity=entity)
+        return {alpha: matvec(tabulation[alpha]) for alpha in tabulation}
 
-        return {alpha: matvec(table)
-                for alpha, table in result.items()}
+    def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
+        result = super().basis_evaluation(order, ps, entity=entity)
+        return self.map_tabulation(result, coordinate_mapping)
 
-    def point_evaluation(self, order, refcoords, entity=None):
-        raise NotImplementedError("TODO: not yet thought about it")
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
+        result = super().point_evaluation(order, refcoords, entity=entity)
+        return self.map_tabulation(result, coordinate_mapping)
 
 
 class DirectlyDefinedElement(NeedsCoordinateMappingElement):
diff --git a/finat/quadrature_element.py b/finat/quadrature_element.py
@@ -149,7 +149,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         sd = self.cell.get_spatial_dimension()
         return {(0,) * sd: gem.ComponentTensor(delta, basis_indices)}
 
-    def point_evaluation(self, order, refcoords, entity=None):
+    def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         raise NotImplementedError("QuadratureElement cannot do point evaluation!")
 
     @property
diff --git a/finat/runtime_tabulated.py b/finat/runtime_tabulated.py
@@ -94,7 +94,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
                 result[alpha] = gem.partial_indexed(gem.Variable(name, shape), ps.indices)
         return result
 
-    def point_evaluation(self, order, point, entity=None):
+    def point_evaluation(self, order, point, entity=None, coordinate_mapping=None):
         raise NotImplementedError("Point evaluation not supported for runtime tabulated elements")
 
     @property
diff --git a/finat/tensor_product.py b/finat/tensor_product.py
@@ -143,7 +143,7 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
 
         return self._merge_evaluations(factor_results)
 
-    def point_evaluation(self, order, point, entity=None):
+    def point_evaluation(self, order, point, entity=None, coordinate_mapping=None):
         entities = self._factor_entity(entity)
         entity_dim, _ = zip(*entities)
 
diff --git a/finat/tensorfiniteelement.py b/finat/tensorfiniteelement.py
@@ -119,9 +119,9 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         scalar_evaluation = self._base_element.basis_evaluation
         return self._tensorise(scalar_evaluation(order, ps, entity, coordinate_mapping=coordinate_mapping))
 
-    def point_evaluation(self, order, point, entity=None):
+    def point_evaluation(self, order, point, entity=None, coordinate_mapping=None):
         scalar_evaluation = self._base_element.point_evaluation
-        return self._tensorise(scalar_evaluation(order, point, entity))
+        return self._tensorise(scalar_evaluation(order, point, entity, coordinate_mapping))
 
     def _tensorise(self, scalar_evaluation):
         # Old basis function and value indices
diff --git a/test/finat/test_point_evaluation.py b/test/finat/test_point_evaluation.py
@@ -0,0 +1,74 @@
+import pytest
+
+import FIAT
+import finat
+import gem
+import numpy
+
+
+@pytest.fixture(params=[1, 2, 3])
+def cell(request):
+    dim = request.param
+    return FIAT.ufc_simplex(dim)
+
+
+@pytest.mark.parametrize('degree', [1, 2])
+def test_cellwise_constant(cell, degree):
+    dim = cell.get_spatial_dimension()
+    element = finat.Lagrange(cell, degree)
+    index = gem.Index()
+    point = gem.partial_indexed(gem.Variable('X', (17, dim)), (index,))
+
+    order = 2
+    for alpha, table in element.point_evaluation(order, point).items():
+        if sum(alpha) < degree:
+            assert table.free_indices == (index,)
+        else:
+            assert table.free_indices == ()
+
+
+@pytest.mark.parametrize("element,degree", [
+    (finat.HsiehCloughTocher, 3),
+    (finat.Argyris, 5),
+    (finat.MardalTaiWinther, 3),
+])
+def test_point_evaluation_zany(ref_to_phys, element, degree):
+    dim = 2
+    ref_cell = ref_to_phys[dim].ref_cell
+    phys_cell = ref_to_phys[dim].phys_cell
+    A, b = FIAT.reference_element.make_affine_mapping(ref_cell.vertices, phys_cell.vertices)
+
+    ref_pt = numpy.array([0.2, 0.3])
+    phys_pt = numpy.dot(A, ref_pt) + b
+
+    finat_kwargs = {}
+    if element in {finat.HsiehCloughTocher, finat.Argyris}:
+        finat_kwargs["avg"] = True
+
+    order = 0
+    point = gem.Literal(ref_pt)
+    ref_element = element(ref_cell, degree, **finat_kwargs)
+    result = ref_element.point_evaluation(order, point, coordinate_mapping=ref_to_phys[dim])
+
+    phys_element = element(phys_cell, degree, **finat_kwargs).fiat_equivalent
+    expected = phys_element.tabulate(order, phys_pt)
+
+    num_dof = ref_element.space_dimension()
+    for alpha in result:
+
+        ref_val, = gem.interpreter.evaluate([result[alpha]])
+
+        if phys_element.mapping()[0] == "covariant piola":
+            val = numpy.tensordot(ref_val.arr, A, (-1, 0))
+        elif phys_element.mapping()[0] == "contravariant piola":
+            detA = numpy.linalg.det(A)
+            val = numpy.tensordot(ref_val.arr, A/detA, (-1, 1))
+        else:
+            val = ref_val.arr
+
+        assert numpy.allclose(val, expected[alpha][:num_dof])
+
+
+if __name__ == '__main__':
+    import os
+    pytest.main(os.path.abspath(__file__))
diff --git a/test/finat/test_point_evaluation_ciarlet.py b/test/finat/test_point_evaluation_ciarlet.py
