Codim 2 coupling

ffcx/codegeneration/C/form.py

@@ -87,7 +87,7 @@ def generator(ir: FormIR, options):
     integral_ids = []
     integral_offsets = [0]
     # Note: the order of this list is defined by the enum ufcx_integral_type in ufcx.h
-    for itg_type in ("cell", "exterior_facet", "interior_facet"):
+    for itg_type in ("cell", "exterior_facet", "interior_facet", "ridge"):
         unsorted_integrals = []
         unsorted_ids = []
         for name, id in zip(ir.integral_names[itg_type], ir.subdomain_ids[itg_type]):

ffcx/codegeneration/C/integrals.py

@@ -90,5 +90,4 @@ def generator(ir: IntegralIR, options):
         tabulate_tensor_complex64=code["tabulate_tensor_complex64"],
         tabulate_tensor_complex128=code["tabulate_tensor_complex128"],
     )
-
     return declaration, implementation

ffcx/codegeneration/access.py

@@ -43,6 +43,7 @@ def __init__(self, entity_type: str, integral_type: str, symbols, options):
             ufl.geometry.FacetEdgeVectors: self.facet_edge_vectors,
             ufl.geometry.CellEdgeVectors: self.cell_edge_vectors,
             ufl.geometry.CellFacetJacobian: self.cell_facet_jacobian,
+            ufl.geometry.CellRidgeJacobian: self.cell_ridge_jacobian,
             ufl.geometry.ReferenceCellVolume: self.reference_cell_volume,
             ufl.geometry.ReferenceFacetVolume: self.reference_facet_volume,
             ufl.geometry.ReferenceCellEdgeVectors: self.reference_cell_edge_vectors,
@@ -70,7 +71,6 @@ def get(
                 if isinstance(e, k):
                     handler = self.call_lookup[k]
                     break
-
         if handler:
             return handler(mt, tabledata, quadrature_rule)  # type: ignore
         else:
@@ -246,6 +246,18 @@ def cell_facet_jacobian(self, mt, tabledata, num_points):
         else:
             raise RuntimeError(f"Unhandled cell types {cellname}.")
 
+    def cell_ridge_jacobian(self, mt, tabledata, num_points):
+        """Access a cell ridge jacobian."""
+        cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
+        if cellname in ("tetrahedron", "prism", "hexahedron"):
+            table = L.Symbol(f"{cellname}_cell_ridge_jacobian", dtype=L.DataType.REAL)
+            ridge = self.symbols.entity("ridge", mt.restriction)
+            return table[ridge][mt.component[0]][mt.component[1]]
+        elif cellname in ["triangle", "quadrilateral"]:
+            raise RuntimeError("The ridge jacobian doesn't make sense for 2D cells.")
+        else:
+            raise RuntimeError(f"Unhandled cell types {cellname}.")
+
     def reference_cell_edge_vectors(self, mt, tabledata, num_points):
         """Access a reference cell edge vector."""
         cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname()
@@ -376,10 +388,10 @@ def facet_edge_vectors(self, mt, tabledata, num_points):
 
         # Get edge vertices
         facet = self.symbols.entity("facet", mt.restriction)
-        facet_edge = mt.component[0]
+        edge_ridge = mt.component[0]
         facet_edge_vertices = L.Symbol(f"{cellname}_facet_edge_vertices", dtype=L.DataType.INT)
-        vertex0 = facet_edge_vertices[facet][facet_edge][0]
-        vertex1 = facet_edge_vertices[facet][facet_edge][1]
+        vertex0 = facet_edge_vertices[facet][edge_ridge][0]
+        vertex1 = facet_edge_vertices[facet][edge_ridge][1]
 
         # Get dofs and component
         component = mt.component[1]

ffcx/codegeneration/definitions.py

@@ -68,6 +68,7 @@ def __init__(self, entity_type: str, integral_type: str, access, options):
             ufl.geometry.FacetEdgeVectors: self.pass_through,
             ufl.geometry.CellEdgeVectors: self.pass_through,
             ufl.geometry.CellFacetJacobian: self.pass_through,
+            ufl.geometry.CellRidgeJacobian: self.pass_through,
             ufl.geometry.ReferenceCellVolume: self.pass_through,
             ufl.geometry.ReferenceFacetVolume: self.pass_through,
             ufl.geometry.ReferenceCellEdgeVectors: self.pass_through,

ffcx/codegeneration/geometry.py

@@ -17,6 +17,8 @@ def write_table(tablename, cellname):
         return facet_edge_vertices(tablename, cellname)
     if tablename == "cell_facet_jacobian":
         return cell_facet_jacobian(tablename, cellname)
+    if tablename == "cell_ridge_jacobian":
+        return cell_ridge_jacobian(tablename, cellname)
     if tablename == "reference_cell_volume":
         return reference_cell_volume(tablename, cellname)
     if tablename == "reference_facet_volume":
@@ -64,6 +66,14 @@ def cell_facet_jacobian(tablename, cellname):
     return L.ArrayDecl(symbol, values=out, const=True)
 
 
+def cell_ridge_jacobian(tablename, cellname):
+    """Write a reference facet jacobian."""
+    celltype = getattr(basix.CellType, cellname)
+    out = basix.cell.edge_jacobians(celltype)
+    symbol = L.Symbol(f"{cellname}_{tablename}", dtype=L.DataType.REAL)
+    return L.ArrayDecl(symbol, values=out, const=True)
+
+
 def reference_cell_volume(tablename, cellname):
     """Write a reference cell volume."""
     celltype = getattr(basix.CellType, cellname)

ffcx/codegeneration/integral_generator.py

@@ -198,6 +198,7 @@ def generate_geometry_tables(self):
         ufl_geometry = {
             ufl.geometry.FacetEdgeVectors: "facet_edge_vertices",
             ufl.geometry.CellFacetJacobian: "cell_facet_jacobian",
+            ufl.geometry.CellRidgeJacobian: "cell_ridge_jacobian",
             ufl.geometry.ReferenceCellVolume: "reference_cell_volume",
             ufl.geometry.ReferenceFacetVolume: "reference_facet_volume",
             ufl.geometry.ReferenceCellEdgeVectors: "reference_cell_edge_vectors",
@@ -221,7 +222,6 @@ def generate_geometry_tables(self):
         for i, cell_list in cells.items():
             for c in cell_list:
                 parts.append(geometry.write_table(ufl_geometry[i], c))
-
         return parts
 
     def generate_element_tables(self):

ffcx/codegeneration/symbols.py

@@ -108,6 +108,8 @@ def entity(self, entity_type, restriction):
                 return self.entity_local_index[0]
         elif entity_type == "vertex":
             return self.entity_local_index[0]
+        elif entity_type == "ridge":
+            return self.entity_local_index[0]
         else:
             logging.exception(f"Unknown entity_type {entity_type}")
 
@@ -136,7 +138,9 @@ def x_component(self, mt):
     def J_component(self, mt):
         """Jacobian component."""
         # FIXME: Add domain number!
-        return L.Symbol(format_mt_name("J", mt), dtype=L.DataType.REAL)
+        return L.Symbol(
+            format_mt_name(f"J{mt.expr.ufl_domain().ufl_id()}", mt), dtype=L.DataType.REAL
+        )
 
     def domain_dof_access(self, dof, component, gdim, num_scalar_dofs, restriction):
         """Domain DOF access."""

ffcx/compiler.py

@@ -122,5 +122,4 @@ def compile_ufl_objects(
     cpu_time = time()
     code_h, code_c = format_code(code)
     _print_timing(4, time() - cpu_time)
-
     return code_h, code_c

ffcx/element_interface.py

@@ -52,3 +52,18 @@ def map_facet_points(
         ],
         dtype=np.float64,
     )
+
+
+def map_edge_points(
+    points: npt.NDArray[np.float64], edge: int, cellname: str
+) -> npt.NDArray[np.float64]:
+    """Map points from a reference edge to a physical edge."""
+    geom = np.asarray(basix.geometry(_CellType[cellname]))
+    edge_vertices = [geom[i] for i in basix.topology(_CellType[cellname])[-3][edge]]
+    return np.asarray(
+        [
+            edge_vertices[0] + sum((i - edge_vertices[0]) * j for i, j in zip(edge_vertices[1:], p))
+            for p in points
+        ],
+        dtype=np.float64,
+    )

ffcx/ir/elementtables.py

@@ -141,11 +141,10 @@ def get_ffcx_table_values(
     # Extract arrays for the right scalar component
     component_tables = []
     component_element, offset, stride = element.get_component_element(flat_component)
-
     for entity in range(num_entities):
         if codim == 0:
             entity_points = map_integral_points(points, integral_type, cell, entity)
-        elif codim == 1:
+        elif codim == 1 or codim == 2:
             entity_points = points
         else:
             raise RuntimeError("Codimension > 1 isn't supported.")
@@ -200,7 +199,7 @@ def generate_psi_table_name(
     if any(derivative_counts):
         name += "_D" + "".join(str(d) for d in derivative_counts)
     name += {None: "", "cell": "_AC", "facet": "_AF"}[averaged]
-    name += {"cell": "", "facet": "_F", "vertex": "_V"}[entity_type]
+    name += {"cell": "", "facet": "_F", "vertex": "_V", "ridge": "_E"}[entity_type]
     name += f"_Q{quadrature_rule.id()}"
     return name
 
@@ -247,8 +246,12 @@ def get_modified_terminal_element(mt) -> typing.Optional[ModifiedTerminalElement
     assert (mt.averaged is None) or not (ld or gd)
     # Change derivatives format for table lookup
     tdim = domain.topological_dimension()
-    local_derivatives: tuple[int, ...] = tuple(ld.count(i) for i in range(tdim))
-
+    local_derivatives: tuple[int, ...]
+    # Special handling if the element is defined on a vertex
+    if tdim == 0:
+        local_derivatives = (0,)
+    else:
+        local_derivatives = tuple(ld.count(i) for i in range(tdim))
     return ModifiedTerminalElement(element, mt.averaged, local_derivatives, fc)
 
 
@@ -357,89 +360,124 @@ def build_optimized_tables(
         tdim = cell.topological_dimension()
         codim = tdim - element.cell.topological_dimension()
         assert codim >= 0
-        if codim > 1:
-            raise RuntimeError("Codimension > 1 isn't supported.")
+        if codim > 2:
+            raise RuntimeError("Codimension > 2 isn't supported.")
 
         # Only permute quadrature rules for interior facets integrals and for
         # the codim zero element in mixed-dimensional integrals. The latter is
         # needed because a cell may see its sub-entities as being oriented
-        # differently to their global orientation
+        # differently to their global orientation#
         if integral_type == "interior_facet" or (is_mixed_dim and codim == 0):
-            if tdim == 1 or codim == 1:
-                # Do not add permutations if codim-1 as facets have already gotten a global
-                # orientation in DOLFINx
-                t = get_ffcx_table_values(
-                    quadrature_rule.points,
-                    cell,
-                    integral_type,
-                    element,
-                    avg,
-                    entity_type,
-                    local_derivatives,
-                    flat_component,
-                    codim,
-                )
-            elif tdim == 2:
-                new_table = []
-                for ref in range(2):
-                    new_table.append(
-                        get_ffcx_table_values(
-                            permute_quadrature_interval(quadrature_rule.points, ref),
-                            cell,
-                            integral_type,
-                            element,
-                            avg,
-                            entity_type,
-                            local_derivatives,
-                            flat_component,
-                            codim,
-                        )
+            if entity_type == "facet":
+                if tdim == 1 or codim == 1 or codim == 2:
+                    # Do not add permutations if codim-1 as facets have already gotten a global
+                    # orientation in DOLFINx
+                    t = get_ffcx_table_values(
+                        quadrature_rule.points,
+                        cell,
+                        integral_type,
+                        element,
+                        avg,
+                        entity_type,
+                        local_derivatives,
+                        flat_component,
+                        codim,
                     )
-
-                t = new_table[0]
-                t["array"] = np.vstack([td["array"] for td in new_table])
-            elif tdim == 3:
-                cell_type = cell.cellname()
-                if cell_type == "tetrahedron":
+                elif tdim == 2:
                     new_table = []
-                    for rot in range(3):
-                        for ref in range(2):
-                            new_table.append(
-                                get_ffcx_table_values(
-                                    permute_quadrature_triangle(quadrature_rule.points, ref, rot),
-                                    cell,
-                                    integral_type,
-                                    element,
-                                    avg,
-                                    entity_type,
-                                    local_derivatives,
-                                    flat_component,
-                                    codim,
-                                )
+                    for ref in range(2):
+                        new_table.append(
+                            get_ffcx_table_values(
+                                permute_quadrature_interval(quadrature_rule.points, ref),
+                                cell,
+                                integral_type,
+                                element,
+                                avg,
+                                entity_type,
+                                local_derivatives,
+                                flat_component,
+                                codim,
                             )
+                        )
+
                     t = new_table[0]
                     t["array"] = np.vstack([td["array"] for td in new_table])
-                elif cell_type == "hexahedron":
-                    new_table = []
-                    for rot in range(4):
-                        for ref in range(2):
-                            new_table.append(
-                                get_ffcx_table_values(
-                                    permute_quadrature_quadrilateral(
-                                        quadrature_rule.points, ref, rot
-                                    ),
-                                    cell,
-                                    integral_type,
-                                    element,
-                                    avg,
-                                    entity_type,
-                                    local_derivatives,
-                                    flat_component,
-                                    codim,
+                elif tdim == 3:
+                    cell_type = cell.cellname()
+                    if cell_type == "tetrahedron":
+                        new_table = []
+                        for rot in range(3):
+                            for ref in range(2):
+                                new_table.append(
+                                    get_ffcx_table_values(
+                                        permute_quadrature_triangle(
+                                            quadrature_rule.points, ref, rot
+                                        ),
+                                        cell,
+                                        integral_type,
+                                        element,
+                                        avg,
+                                        entity_type,
+                                        local_derivatives,
+                                        flat_component,
+                                        codim,
+                                    )
+                                )
+                        t = new_table[0]
+                        t["array"] = np.vstack([td["array"] for td in new_table])
+                    elif cell_type == "hexahedron":
+                        new_table = []
+                        for rot in range(4):
+                            for ref in range(2):
+                                new_table.append(
+                                    get_ffcx_table_values(
+                                        permute_quadrature_quadrilateral(
+                                            quadrature_rule.points, ref, rot
+                                        ),
+                                        cell,
+                                        integral_type,
+                                        element,
+                                        avg,
+                                        entity_type,
+                                        local_derivatives,
+                                        flat_component,
+                                        codim,
+                                    )
                                 )
+                        t = new_table[0]
+                        t["array"] = np.vstack([td["array"] for td in new_table])
+            elif entity_type == "ridge":
+                if tdim > 2:
+                    new_table = []
+                    for ref in range(2):
+                        new_table.append(
+                            get_ffcx_table_values(
+                                permute_quadrature_interval(quadrature_rule.points, ref),
+                                cell,
+                                integral_type,
+                                element,
+                                avg,
+                                entity_type,
+                                local_derivatives,
+                                flat_component,
+                                codim,
                             )
+                        )
                     t = new_table[0]
                     t["array"] = np.vstack([td["array"] for td in new_table])
+                else:
+                    # If ridge integral over vertex no permutation is needed
+                    t = get_ffcx_table_values(
+                        quadrature_rule.points,
+                        cell,
+                        integral_type,
+                        element,
+                        avg,
+                        entity_type,
+                        local_derivatives,
+                        flat_component,
+                        codim,
+                    )
         else:
             t = get_ffcx_table_values(
                 quadrature_rule.points,