Merge branch 'main' into pbrubeck/fix/integral_type_map

Author: pbrubeck

firedrake/__init__.py

@@ -67,7 +67,7 @@ def init_petsc():
 from firedrake.cofunction import Cofunction, RieszMap  # noqa: F401
 from firedrake.constant import Constant  # noqa: F401
 from firedrake.deflation import DeflatedSNES, Deflation  # noqa: F401
-from firedrake.exceptions import ConvergenceError  # noqa: F401
+from firedrake.exceptions import ConvergenceError, MismatchingDomainError  # noqa: F401
 from firedrake.function import (  # noqa: F401
     Function, PointNotInDomainError,
     CoordinatelessFunction, PointEvaluator

firedrake/exceptions.py

@@ -1,3 +1,4 @@
+from tsfc.exceptions import MismatchingDomainError  # noqa: F401
 
 
 class ConvergenceError(Exception):

firedrake/interpolation.py

@@ -34,7 +34,7 @@
 from tsfc.driver import compile_expression_dual_evaluation
 from tsfc.ufl_utils import extract_firedrake_constants, hash_expr
 
-from firedrake.utils import IntType, ScalarType, known_pyop2_safe, tuplify
+from firedrake.utils import IntType, ScalarType, cached_property, known_pyop2_safe, tuplify
 from firedrake.tsfc_interface import extract_numbered_coefficients, _cachedir
 from firedrake.ufl_expr import Argument, Coargument, action
 from firedrake.mesh import MissingPointsBehaviour, VertexOnlyMeshMissingPointsError, VertexOnlyMeshTopology, MeshGeometry, MeshTopology, VertexOnlyMesh
@@ -155,6 +155,57 @@ def options(self) -> InterpolateOptions:
         """
         return self._options
 
+    @cached_property
+    def _interpolator(self):
+        """Access the numerical interpolator.
+
+        Returns
+        -------
+        Interpolator
+            An appropriate :class:`Interpolator` subclass for this
+            interpolation expression.
+        """
+        arguments = self.arguments()
+        has_mixed_arguments = any(len(arg.function_space()) > 1 for arg in arguments)
+        if len(arguments) == 2 and has_mixed_arguments:
+            return MixedInterpolator(self)
+
+        operand, = self.ufl_operands
+        target_mesh = self.target_space.mesh()
+
+        try:
+            source_mesh = extract_unique_domain(operand) or target_mesh
+        except ValueError:
+            raise NotImplementedError(
+                "Interpolating an expression with no arguments defined on multiple meshes is not implemented yet."
+            )
+
+        try:
+            target_mesh = target_mesh.unique()
+            source_mesh = source_mesh.unique()
+        except RuntimeError:
+            return MixedInterpolator(self)
+
+        submesh_interp_implemented = (
+            all(isinstance(m.topology, MeshTopology) for m in [target_mesh, source_mesh])
+            and target_mesh.submesh_ancesters[-1] is source_mesh.submesh_ancesters[-1]
+            and target_mesh.topological_dimension == source_mesh.topological_dimension
+        )
+        if target_mesh is source_mesh or submesh_interp_implemented:
+            return SameMeshInterpolator(self)
+
+        if isinstance(target_mesh.topology, VertexOnlyMeshTopology):
+            if isinstance(source_mesh.topology, VertexOnlyMeshTopology):
+                return VomOntoVomInterpolator(self)
+            if target_mesh.geometric_dimension != source_mesh.geometric_dimension:
+                raise ValueError("Cannot interpolate onto a VertexOnlyMesh of a different geometric dimension.")
+            return SameMeshInterpolator(self)
+
+        if has_mixed_arguments or len(self.target_space) > 1:
+            return MixedInterpolator(self)
+
+        return CrossMeshInterpolator(self)
+
 
 @PETSc.Log.EventDecorator()
 def interpolate(expr: Expr, V: WithGeometry | BaseForm, **kwargs) -> Interpolate:
@@ -353,46 +404,7 @@ def get_interpolator(expr: Interpolate) -> Interpolator:
         An appropriate :class:`Interpolator` subclass for the given
         interpolation expression.
     """
-    arguments = expr.arguments()
-    has_mixed_arguments = any(len(arg.function_space()) > 1 for arg in arguments)
-    if len(arguments) == 2 and has_mixed_arguments:
-        return MixedInterpolator(expr)
-
-    operand, = expr.ufl_operands
-    target_mesh = expr.target_space.mesh()
-
-    try:
-        source_mesh = extract_unique_domain(operand) or target_mesh
-    except ValueError:
-        raise NotImplementedError(
-            "Interpolating an expression with no arguments defined on multiple meshes is not implemented yet."
-        )
-
-    try:
-        target_mesh = target_mesh.unique()
-        source_mesh = source_mesh.unique()
-    except RuntimeError:
-        return MixedInterpolator(expr)
-
-    submesh_interp_implemented = (
-        all(isinstance(m.topology, MeshTopology) for m in [target_mesh, source_mesh])
-        and target_mesh.submesh_ancesters[-1] is source_mesh.submesh_ancesters[-1]
-        and target_mesh.topological_dimension == source_mesh.topological_dimension
-    )
-    if target_mesh is source_mesh or submesh_interp_implemented:
-        return SameMeshInterpolator(expr)
-
-    if isinstance(target_mesh.topology, VertexOnlyMeshTopology):
-        if isinstance(source_mesh.topology, VertexOnlyMeshTopology):
-            return VomOntoVomInterpolator(expr)
-        if target_mesh.geometric_dimension != source_mesh.geometric_dimension:
-            raise ValueError("Cannot interpolate onto a VertexOnlyMesh of a different geometric dimension.")
-        return SameMeshInterpolator(expr)
-
-    if has_mixed_arguments or len(expr.target_space) > 1:
-        return MixedInterpolator(expr)
-
-    return CrossMeshInterpolator(expr)
+    return expr._interpolator
 
 
 class DofNotDefinedError(Exception):

firedrake/mesh.py

@@ -3123,12 +3123,11 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F
             pyop2_index.extend(cell_node_map.values[ngidx])
 
         # Find the correct coordinate permutation for each cell
-        # NB: Coordinates must be cast to real when running Firedrake in complex mode
         permutation = find_permutation(
             physical_space_points,
-            new_coordinates.dat.data[pyop2_index].reshape(
+            new_coordinates.dat.data[pyop2_index].real.reshape(
                 physical_space_points.shape
-            ).astype(np.float64, copy=False),
+            ),
             tol=permutation_tol
         )
 

tests/firedrake/regression/test_interpolate.py

@@ -592,9 +592,12 @@ def test_interpolator_reuse(family, degree, mode):
         u = Function(V.dual())
         expr = interpolate(TestFunction(V), u)
 
-    I = get_interpolator(expr)
+    Iorig = get_interpolator(expr)
 
     for k in range(3):
+        I = get_interpolator(expr)
+        assert I is Iorig
+
         u.assign(rg.uniform(u.function_space()))
         expected = u.dat.data.copy()
 

tests/firedrake/regression/test_multiple_domains.py

@@ -41,10 +41,10 @@ def test_mismatching_meshes_indexed_function(mesh1, mesh3):
     with pytest.raises(NotImplementedError):
         project(d1, target)
 
-    with pytest.raises(ValueError):
+    with pytest.raises(MismatchingDomainError):
         assemble(inner(d1, TestFunction(V2))*dx(domain=mesh3))
 
-    with pytest.raises(ValueError):
+    with pytest.raises(MismatchingDomainError):
         assemble(inner(d1, TestFunction(V2))*dx(domain=mesh1))
 
 
@@ -177,29 +177,29 @@ def test_multi_domain_assemble():
 
     for i, j in [(0, 1), (1, 0)]:
         a1 = inner(u[i], v[j])*dx(domain=mesh1)
-        with pytest.raises(ValueError):
+        with pytest.raises(MismatchingDomainError):
             assemble(a1)
         a2 = inner(u[i], v[j])*dx(domain=mesh2)
-        with pytest.raises(ValueError):
+        with pytest.raises(MismatchingDomainError):
             assemble(a2)
         l1 = inner(f[i], v[j])*dx(domain=mesh1)
-        with pytest.raises(ValueError):
+        with pytest.raises(MismatchingDomainError):
             assemble(l1)
         l2 = inner(f[i], v[j])*dx(domain=mesh2)
-        with pytest.raises(ValueError):
+        with pytest.raises(MismatchingDomainError):
             assemble(l2)
 
     for i, j in [(0, 0), (1, 1)]:
         a = inner(u[i], v[j])*dx(domain=mesh1)
         if i == 1:
-            with pytest.raises(ValueError):
+            with pytest.raises(MismatchingDomainError):
                 assemble(a)
             continue
         A = assemble(a)
         assert A.M.values.shape == (V.dim(), V.dim())
 
     a = inner(u[0], v[0])*dx(domain=mesh1) + inner(u[0], v[1])*dx(domain=mesh2)
-    with pytest.raises(ValueError):
+    with pytest.raises(MismatchingDomainError):
         assemble(a)
 
     a = inner(u[0], v[0])*dx(domain=mesh1) + inner(u[1], v[1])*dx(domain=mesh2)

tests/firedrake/regression/test_projection_zany.py

@@ -80,7 +80,8 @@ def run_convergence_test(mh, el, degree, convrate):
                           ('Hermite', 3, 3.8),
                           ('Bell', 5, 4.7),
                           ('Argyris', 5, 5.8),
-                          ('Argyris', 6, 6.7)])
+                          ('Argyris', 6, 6.7),
+                          ('Nonconforming Robust Wu-Xu', 7, 3.8)])
 def test_projection_zany_convergence_2d(hierarchy_2d, el, deg, convrate):
     run_convergence_test(hierarchy_2d[2:], el, deg, convrate)
 
@@ -97,7 +98,9 @@ def test_projection_zany_convergence_3d(hierarchy_3d, el, deg, convrate):
                           ('HCT', 3),
                           ('HCT', 4),
                           ('Argyris', 5),
-                          ('Argyris', 6)])
+                          ('Argyris', 6),
+                          ('Nonconforming Wu-Xu', 4),
+                          ('Alfeld C2', 5)])
 def test_mass_conditioning(element, degree, hierarchy_2d):
     mass_cond = []
     for msh in hierarchy_2d[1:4]:

tsfc/__init__.py

@@ -1,5 +1,6 @@
 from tsfc.driver import compile_form, compile_expression_dual_evaluation  # noqa: F401
 from tsfc.parameters import default_parameters  # noqa: F401
+from tsfc.exceptions import MismatchingDomainError  # noqa: F401
 
 
 def register_citations():

tsfc/driver.py

@@ -23,6 +23,8 @@
 from tsfc.parameters import default_parameters, is_complex
 from tsfc.ufl_utils import apply_mapping, extract_firedrake_constants
 import tsfc.kernel_interface.firedrake_loopy as firedrake_interface_loopy
+from tsfc.exceptions import MismatchingDomainError
+
 
 # To handle big forms. The various transformations might need a deeper stack
 sys.setrecursionlimit(3000)
@@ -147,7 +149,7 @@ def compile_integral(integral_data, form_data, prefix, parameters, *, diagonal=F
         integral_type=integral_data.integral_type,
         subdomain_id=integral_data.subdomain_id,
         domain_number=domain_number,
-        domain_integral_type_map=integral_data.domain_integral_type_map,
+        domain_integral_type_map={mesh: integral_data.domain_integral_type_map.get(mesh, None) for mesh in all_meshes},
         arguments=arguments,
         coefficients=coefficients,
         coefficient_split=coefficient_split,
@@ -188,19 +190,19 @@ def validate_domains(form):
         domain = itg.ufl_domain()
         for other_domain in itg.extra_domain_integral_type_map():
             if domain.submesh_youngest_common_ancester(other_domain) is None:
-                raise ValueError("Assembly of forms over unrelated meshes is not supported. "
-                                 "Try using Submeshes or cross-mesh interpolation.")
+                raise MismatchingDomainError("Assembly of forms over unrelated meshes is not supported. "
+                                             "Try using Submeshes or cross-mesh interpolation.")
 
         # Check that all Arguments and Coefficients are defined on the valid domains
         valid_domains = set(itg.extra_domain_integral_type_map())
         valid_domains.add(domain)
 
         itg_domains = set(extract_domains(itg))
         if len(itg_domains - valid_domains) > 0:
-            raise ValueError("Argument or Coefficient domain not found in integral. "
-                             "Possibly, the form contains coefficients on different meshes "
-                             "and requires measure intersection, for example: "
-                             'Measure("dx", argument_mesh, intersect_measures=[Measure("dx", coefficient_mesh)]).')
+            raise MismatchingDomainError("Argument or Coefficient domain not found in integral. "
+                                         "Possibly, the form contains coefficients on different meshes "
+                                         "and requires measure intersection, for example: "
+                                         'Measure("dx", argument_mesh, intersect_measures=[Measure("dx", coefficient_mesh)]).')
 
 
 def preprocess_parameters(parameters):

tsfc/exceptions.py

@@ -0,0 +1,4 @@
+
+
+class MismatchingDomainError(Exception):
+    """Error raised for unsupported multidomain problems"""