Working towards 3D-1D

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.CellEdgeJacobian: self.cell_edge_jacobian
             ufl.geometry.ReferenceCellVolume: self.reference_cell_volume,
             ufl.geometry.ReferenceFacetVolume: self.reference_facet_volume,
             ufl.geometry.ReferenceCellEdgeVectors: self.reference_cell_edge_vectors,

ffcx/element_interface.py

@@ -52,3 +52,19 @@ 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", "edge": "_E"}[entity_type]
     name += f"_Q{quadrature_rule.id()}"
     return name
 
@@ -357,29 +356,91 @@ 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:
+            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,
+                    )
+                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,
+                            )
+                        )
+
+                    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":
+                        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 == "edge":
                 new_table = []
                 for ref in range(2):
                     new_table.append(
@@ -398,48 +459,6 @@ def build_optimized_tables(
 
                 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":
-                    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])
         else:
             t = get_ffcx_table_values(
                 quadrature_rule.points,

ffcx/ir/integral.py

@@ -60,7 +60,6 @@ def compute_integral_ir(cell, integral_type, entity_type, integrands, argument_s
 
     for quadrature_rule, integrand in integrands.items():
         expression = integrand
-
         # Rebalance order of nested terminal modifiers
         expression = balance_modifiers(expression)
 
@@ -87,7 +86,7 @@ def compute_integral_ir(cell, integral_type, entity_type, integrands, argument_s
         for domain in ufl.domain.extract_domains(integrand):
             if domain.topological_dimension() != cell.topological_dimension():
                 is_mixed_dim = True
-
+        breakpoint()
         mt_table_reference = build_optimized_tables(
             quadrature_rule,
             cell,

ffcx/ir/representation.py

@@ -196,6 +196,7 @@ def _compute_integral_ir(
         "interior_facet": "facet",
         "vertex": "vertex",
         "custom": "cell",
+        "edge": "edge"
     }
 
     # Iterate over groups of integrals
@@ -267,10 +268,15 @@ def _compute_integral_ir(
                 # prescribed in certain cases.
 
                 degree = md["quadrature_degree"]
-                if integral_type != "cell":
+                if integral_type == "facet":
                     facet_types = cell.facet_types()
                     assert len(facet_types) == 1
                     cellname = facet_types[0].cellname()
+                elif integral_type == "edge":
+                    edge_types = cell.edge_types()
+                    assert len(edge_types) == 1
+                    cellname = edge_types[0].cellname()
+
                 if degree > 1:
                     warnings.warn(
                         "Explicitly selected vertex quadrature (degree 1), "
@@ -463,7 +469,7 @@ def _compute_form_ir(
     # it has to know their names for codegen phase
     ir["integral_names"] = {}
     ir["subdomain_ids"] = {}
-    ufcx_integral_types = ("cell", "exterior_facet", "interior_facet")
+    ufcx_integral_types = ("cell", "exterior_facet", "interior_facet", "edge")
     ir["subdomain_ids"] = {itg_type: [] for itg_type in ufcx_integral_types}
     ir["integral_names"] = {itg_type: [] for itg_type in ufcx_integral_types}
     for itg_index, itg_data in enumerate(form_data.integral_data):

ffcx/ir/representationutils.py

@@ -12,7 +12,7 @@
 import numpy as np
 import ufl
 
-from ffcx.element_interface import create_quadrature, map_facet_points, reference_cell_vertices
+from ffcx.element_interface import create_quadrature, map_facet_points, reference_cell_vertices, map_edge_points
 
 logger = logging.getLogger("ffcx")
 
@@ -56,7 +56,6 @@ def create_quadrature_points_and_weights(
     pts = None
     wts = None
     tensor_factors = None
-
     if integral_type == "cell":
         if cell.cellname() in ["quadrilateral", "hexahedron"] and use_tensor_product:
             if cell.cellname() == "quadrilateral":
@@ -78,7 +77,13 @@ def create_quadrature_points_and_weights(
         # Raise exception for cells with more than one facet type e.g. prisms
         if len(facet_types) > 1:
             raise Exception(f"Cell type {cell} not supported for integral type {integral_type}.")
-        pts, wts = create_quadrature(facet_types[0].cellname(), degree, rule, elements)
+        pts, wts = create_quadrature(facet_types[0].cellname(), degree, rule, elements) 
+    elif integral_type in ufl.measure.edge_integral_types:
+        edge_types = cell.edge_types()
+        if len(edge_types) > 1:
+            raise Exception(f"Cell type {cell} not supported for integral type {integral_type}.")
+        pts, wts = create_quadrature(edge_types[0].cellname(), degree, rule, elements)
+        print(pts, wts, edge_types)
     elif integral_type in ufl.measure.point_integral_types:
         pts, wts = create_quadrature("vertex", degree, rule, elements)
     elif integral_type == "expression":
@@ -95,6 +100,8 @@ def integral_type_to_entity_dim(integral_type, tdim):
         entity_dim = tdim
     elif integral_type in ufl.measure.facet_integral_types:
         entity_dim = tdim - 1
+    elif integral_type in ufl.measure.edge_integral_types:
+        entity_dim = tdim - 2
     elif integral_type in ufl.measure.point_integral_types:
         entity_dim = 0
     elif integral_type in ufl.custom_integral_types:
@@ -117,6 +124,9 @@ def map_integral_points(points, integral_type, cell, entity):
     elif entity_dim == tdim - 1:
         assert points.shape[1] == tdim - 1
         return np.asarray(map_facet_points(points, entity, cell.cellname()))
+    elif entity_dim == tdim - 2:
+        assert points.shape[1] == tdim - 2
+        return np.asarray(map_edge_points(points, entity, cell.cellname()))
     elif entity_dim == 0:
         return np.asarray([reference_cell_vertices(cell.cellname())[entity]])
     else: