Fix expressions on facets and interpolation on facet submeshes

src/scifem/bcs.py

@@ -3,10 +3,56 @@
 import ufl
 import numpy.typing as npt
 import numpy as np
+from packaging.version import Version
+from ffcx.ir.elementtables import (
+    permute_quadrature_interval,
+    permute_quadrature_triangle,
+    permute_quadrature_quadrilateral,
+)
 
 __all__ = ["interpolate_function_onto_facet_dofs"]
 
 
+def build_quadrature_permutations(facet_type, points):
+    if Version(dolfinx.__version__) < Version("0.11.0.dev0"):
+        # In older versions of dolfinx, the permutation is handled internally
+        # in the C++ code, so we can just return the original points
+        if facet_type == basix.CellType.interval:
+            num_permutations = 2
+        elif facet_type == basix.CellType.triangle:
+            num_permutations = 6
+        elif facet_type == basix.CellType.quadrilateral:
+            num_permutations = 8
+        else:
+            raise ValueError(f"Unsupported {facet_type=}")
+        return [points for _ in range(num_permutations)]
+    else:
+        if facet_type == basix.CellType.interval:
+            return [permute_quadrature_interval(points, ref) for ref in range(2)]
+        elif facet_type == basix.CellType.triangle:
+            perms = []
+            # FFCx order: rot is outer loop, ref is inner loop
+            for rot in range(3):
+                for ref in range(2):
+                    # Counteract the mapping with the inverse permutation
+                    rot_inv = (3 - rot) % 3 if ref == 0 else rot
+                    perms.append(permute_quadrature_triangle(points, ref, rot_inv))
+            return perms
+
+        elif facet_type == basix.CellType.quadrilateral:
+            perms = []
+            # FFCx order: rot is outer loop, ref is inner loop
+            for rot in range(4):
+                for ref in range(2):
+                    # Counteract the mapping with the inverse permutation
+                    rot_inv = (4 - rot) % 4 if ref == 0 else rot
+                    perms.append(permute_quadrature_quadrilateral(points, ref, rot_inv))
+            return perms
+
+        else:
+            raise ValueError(f"Unsupported {facet_type=}")
+
+
 def interpolate_function_onto_facet_dofs(
     Q: dolfinx.fem.FunctionSpace,
     expr: ufl.core.expr.Expr,
@@ -49,17 +95,17 @@ def interpolate_function_onto_facet_dofs(
     )
     ref_top = ref_cmap.reference_topology
     ref_geom = ref_cmap.reference_geometry
+    facet_type = facet_types.pop()
     facet_cmap = basix.ufl.element(
         "Lagrange",
-        facet_types.pop(),
+        facet_type,
         1,
         shape=(domain.topology.dim,),
         dtype=np.float64,
     )
     facet_cel = dolfinx.fem.CoordinateElement(
         dolfinx.cpp.fem.CoordinateElement_float64(facet_cmap.basix_element._e)
     )
-
     reference_facet_points = None
     for i, points in enumerate(interpolation_points[fdim]):
         geom = ref_geom[ref_top[fdim][i]]
@@ -71,9 +117,10 @@ def interpolate_function_onto_facet_dofs(
         else:
             assert np.allclose(reference_facet_points, ref_points)
     assert reference_facet_points is not None
-    # Create expression for BC
-    normal_expr = dolfinx.fem.Expression(expr, reference_facet_points)
-
+    facet_points = build_quadrature_permutations(facet_type, reference_facet_points)
+    expressions = []
+    for i, perm in enumerate(facet_points):
+        expressions.append(dolfinx.fem.Expression(expr, perm))
     points_per_entity = [sum(ip.shape[0] for ip in ips) for ips in interpolation_points]
     offsets = np.zeros(domain.topology.dim + 2, dtype=np.int32)
     offsets[1:] = np.cumsum(points_per_entity[: domain.topology.dim + 1])
@@ -85,14 +132,22 @@ def interpolate_function_onto_facet_dofs(
     all_connected_cells = dolfinx.mesh.compute_incident_entities(
         domain.topology, facets, domain.topology.dim - 1, domain.topology.dim
     )
-    values = np.zeros(len(all_connected_cells) * offsets[-1] * domain.geometry.dim)
+    expr_value_size = expressions[0].value_size
+
+    # Update array allocations to use the exact expression value size
+    values_per_entity = np.zeros((offsets[-1], expr_value_size), dtype=dolfinx.default_scalar_type)
+    values = np.zeros(len(all_connected_cells) * offsets[-1] * expr_value_size)
+
     domain.topology.create_connectivity(domain.topology.dim, fdim)
     c_to_f = domain.topology.connectivity(domain.topology.dim, fdim)
     num_facets_on_process = (
         domain.topology.index_map(fdim).size_local + domain.topology.index_map(fdim).num_ghosts
     )
     is_marked = np.zeros(num_facets_on_process, dtype=np.int8)
     is_marked[facets] = 1
+    domain.topology.create_entity_permutations()
+    num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(domain.topology.cell_type, fdim)
+    facet_permutations = domain.topology.get_facet_permutations().reshape(-1, num_facets_per_cell)
     for i, cell in enumerate(all_connected_cells):
         values_per_entity[:] = 0.0
         local_facets = c_to_f.links(cell)
@@ -102,24 +157,24 @@ def interpolate_function_onto_facet_dofs(
             insert_pos = offsets[fdim] + reference_facet_points.shape[0] * j
             # Backwards compatibility
             entity = np.array([[cell, j]], dtype=np.int32)
+            perm = facet_permutations[cell, j]
             try:
-                normal_on_facet = normal_expr.eval(domain, entity)
+                normal_on_facet = expressions[perm].eval(domain, entity)
             except (AttributeError, AssertionError):
-                normal_on_facet = normal_expr.eval(domain, entity.flatten())
+                normal_on_facet = expressions[perm].eval(domain, entity.flatten())
             # NOTE: evaluate within loop to avoid large memory requirements
             values_per_entity[insert_pos : insert_pos + reference_facet_points.shape[0]] = (
-                normal_on_facet.reshape(-1, domain.geometry.dim)
+                normal_on_facet.reshape(-1, expr_value_size)
             )
-        values[
-            i * offsets[-1] * domain.geometry.dim : (i + 1) * offsets[-1] * domain.geometry.dim
-        ] = values_per_entity.reshape(-1)
+        values[i * offsets[-1] * expr_value_size : (i + 1) * offsets[-1] * expr_value_size] = (
+            values_per_entity.reshape(-1)
+        )
 
     qh = dolfinx.fem.Function(Q)
     if hasattr(qh._cpp_object, "interpolate_f"):
         interpolate_func = qh._cpp_object.interpolate_f
     else:
         interpolate_func = qh._cpp_object.interpolate
-    interpolate_func(values.reshape(-1, domain.geometry.dim).T.copy(), all_connected_cells)
+    interpolate_func(values.reshape(-1, expr_value_size).T.copy(), all_connected_cells)
     qh.x.scatter_forward()
-
     return qh

src/scifem/interpolation.py

@@ -337,7 +337,10 @@ def interpolate_to_surface_submesh(
     """
     Interpolate a function `u_volume` into the function `u_surface`.
     Note:
-        Does not work for DG as no dofs are associated with the facets
+        Does not work for DG as no dofs are associated with the facets in versions of DOLFINx
+        prior to https://github.com/FEniCS/dolfinx/pull/4140, which is included in version
+        0.11.0 and later.
+
     Args:
         u_volume: Function to interpolate data from
         u_surface: Function to interpolate data to
@@ -362,18 +365,21 @@ def interpolate_to_surface_submesh(
 
     submesh.topology.create_entity_permutations()
     mesh.topology.create_entity_permutations()
-    cell_info = mesh.topology.get_cell_permutation_info()
     ft = V_surf.element.basix_element.cell_type
 
-    for i in range(integration_entities.shape[0]):
-        perm = np.arange(data.shape[1], dtype=np.int32)
-        V_vol.element.basix_element.permute_subentity_closure_inv(
-            perm,
-            cell_info[integration_entities[i, 0]],
-            ft,
-            int(integration_entities[i, 1]),
-        )
-        data[i] = data[i][perm]
+    if Version(dolfinx.__version__) < Version("0.11.0.dev0"):
+        # Before the introduction of https://github.com/FEniCS/dolfinx/pull/4140
+        # one needed to permute the data according to the facet permutations.
+        cell_info = mesh.topology.get_cell_permutation_info()
+        for i in range(integration_entities.shape[0]):
+            perm = np.arange(data.shape[1], dtype=np.int32)
+            V_vol.element.basix_element.permute_subentity_closure_inv(
+                perm,
+                cell_info[integration_entities[i, 0]],
+                ft,
+                int(integration_entities[i, 1]),
+            )
+            data[i] = data[i][perm]
 
     if len(data.shape) == 3:
         # Data is now (num_cells, value_size,num_points)

tests/test_bcs.py

@@ -26,9 +26,9 @@ def right_facets(x):
 def test_normal_enforcement(cell_type: dolfinx.mesh.CellType):
     tdim = dolfinx.mesh.cell_dim(cell_type)
     if tdim == 2:
-        mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
+        mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, cell_type=cell_type)
     elif tdim == 3:
-        mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3)
+        mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3, cell_type=cell_type)
 
     V = dolfinx.fem.functionspace(mesh, ("BDM", 1))
 
@@ -59,9 +59,9 @@ def test_normal_enforcement(cell_type: dolfinx.mesh.CellType):
 def test_tangent_enforcement(cell_type: dolfinx.mesh.CellType):
     tdim = dolfinx.mesh.cell_dim(cell_type)
     if tdim == 2:
-        mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
+        mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, cell_type=cell_type)
     elif tdim == 3:
-        mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3)
+        mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 8, 10, 3, cell_type=cell_type)
     else:
         raise ValueError(f"Unsupported {cell_type=}")
 

tests/test_interpolation.py

@@ -231,13 +231,16 @@ def test_discrete_curl(degree, use_petsc, cell_type):
 
 
 @pytest.mark.parametrize("degree", [1, 2, 3])
-@pytest.mark.parametrize("family", ["Lagrange"])
+@pytest.mark.parametrize("family", ["Lagrange", "DG"])
 @pytest.mark.skipif(
     Version(dolfinx.__version__) < Version("0.10.0"), reason="Requires DOLFINx version >0.10.0"
 )
 def test_interpolate_to_interface_submesh(family, degree):
     # Create a unit square
     comm = MPI.COMM_WORLD
+
+    if Version(dolfinx.__version__) < Version("0.11.0.dev0") and family == "DG":
+        pytest.skip("Interpolation to surface submesh does not work for DG in DOLFINx < 0.11.0")
     domain = dolfinx.mesh.create_unit_square(
         comm, 48, 48, ghost_mode=dolfinx.mesh.GhostMode.shared_facet
     )