Update Basix quadrature interface

demo/ExpressionInterpolation.py

@@ -52,7 +52,7 @@
 
 # Find quadrature points for quadrature element
 b_rule = basix.quadrature.string_to_type(q_rule)
-quadrature_points, _ = basix.make_quadrature(b_rule, b_cell, q_degree)
+quadrature_points, _ = basix.quadrature.make_quadrature(b_cell, q_degree, rule=b_rule)
 
 # Get interpolation points for output space
 family = basix.finite_element.string_to_family("Lagrange", cell)

ffcx/codegeneration/basix_custom_element_template.py

@@ -30,7 +30,8 @@
   .discontinuous = {discontinuous},
   .highest_complete_degree = {highest_complete_degree},
   .interpolation_nderivs = {interpolation_nderivs},
-  .highest_degree = {highest_degree}
+  .highest_degree = {highest_degree},
+  .polyset_type = {polyset_type}
 }};
 
 // End of code for custom element {factory_name}

ffcx/codegeneration/finite_element.py

@@ -119,6 +119,7 @@ def generate_custom_element(name, ir):
     d = {}
     d["factory_name"] = name
     d["cell_type"] = int(ir.cell_type)
+    d["polyset_type"] = int(ir.polyset_type)
     d["map_type"] = int(ir.map_type)
     d["sobolev_space"] = int(ir.sobolev_space)
     d["highest_complete_degree"] = ir.highest_complete_degree

ffcx/codegeneration/ufcx.h

@@ -200,6 +200,9 @@ extern "C"
 
     /// The highest degree of a polynomial in the element
     int highest_degree;
+
+    /// The polyset type of the element
+    int polyset_type;
   } ufcx_basix_custom_finite_element;
 
   typedef struct ufcx_dofmap

ffcx/element_interface.py

@@ -71,14 +71,19 @@ def basix_index(indices: typing.Tuple[int]) -> int:
     return basix.index(*indices)
 
 
-def create_quadrature(cellname, degree, rule) -> typing.Tuple[npt.NDArray[np.float64],
-                                                              npt.NDArray[np.float64]]:
+def create_quadrature(
+    cellname: str, degree: int, rule: str, elements: typing.List[basix.ufl._ElementBase]
+) -> typing.Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
     """Create a quadrature rule."""
     if cellname == "vertex":
         return (np.ones((1, 0), dtype=np.float64), np.ones(1, dtype=np.float64))
     else:
-        return basix.make_quadrature(basix.quadrature.string_to_type(rule),
-                                     basix.cell.string_to_type(cellname), degree)
+        celltype = basix.cell.string_to_type(cellname)
+        polyset_type = basix.PolysetType.standard
+        for e in elements:
+            polyset_type = basix.polyset_superset(celltype, polyset_type, e.polyset_type)
+        return basix.make_quadrature(
+            celltype, degree, rule=basix.quadrature.string_to_type(rule), polyset_type=polyset_type)
 
 
 def reference_cell_vertices(cellname: str) -> npt.NDArray[np.float64]:
@@ -111,7 +116,7 @@ def __init__(self, cellname: str, value_shape: typing.Tuple[int, ...], scheme: t
             assert points is None
             assert weights is None
             repr = f"QuadratureElement({cellname}, {scheme}, {degree})"
-            self._points, self._weights = create_quadrature(cellname, degree, scheme)
+            self._points, self._weights = create_quadrature(cellname, degree, scheme, [])
         else:
             assert degree is None
             assert points is not None
@@ -262,6 +267,11 @@ def map_type(self) -> basix.MapType:
         """The Basix map type."""
         return basix.MapType.identity
 
+    @property
+    def polyset_type(self) -> basix.PolysetType:
+        """The polyset type of the element."""
+        raise NotImplementedError()
+
 
 class RealElement(basix.ufl._ElementBase):
     """A real element."""

ffcx/ir/elementtables.py

@@ -109,7 +109,7 @@ def get_ffcx_table_values(points, cell, integral_type, element, avg, entitytype,
         else:
             # Make quadrature rule and get points and weights
             points, weights = create_quadrature_points_and_weights(
-                integral_type, cell, element.highest_degree(), "default")
+                integral_type, cell, element.highest_degree(), "default", [element])
 
     # Tabulate table of basis functions and derivatives in points for each entity
     tdim = cell.topological_dimension()

ffcx/ir/representation.py

@@ -70,6 +70,7 @@ class CustomElementIR(typing.NamedTuple):
     discontinuous: bool
     highest_complete_degree: int
     highest_degree: int
+    polyset_type: basix.PolysetType
 
 
 class ElementIR(typing.NamedTuple):
@@ -274,6 +275,7 @@ def _compute_custom_element_ir(basix_element: basix.finite_element.FiniteElement
     ir["interpolation_nderivs"] = basix_element.interpolation_nderivs
     ir["highest_complete_degree"] = basix_element.highest_complete_degree
     ir["highest_degree"] = basix_element.highest_degree
+    ir["polyset_type"] = basix_element.polyset_type
 
     return CustomElementIR(**ir)
 
@@ -435,7 +437,7 @@ def _compute_integral_ir(form_data, form_index, element_numbers, integral_names,
             else:
                 degree = md["quadrature_degree"]
                 points, weights = create_quadrature_points_and_weights(
-                    integral_type, cell, degree, scheme)
+                    integral_type, cell, degree, scheme, [convert_element(e) for e in form_data.argument_elements])
 
             points = np.asarray(points)
             weights = np.asarray(weights)

ffcx/ir/representationutils.py

@@ -44,18 +44,18 @@ def id(self):
         return self.hash_obj.hexdigest()[-3:]
 
 
-def create_quadrature_points_and_weights(integral_type, cell, degree, rule):
+def create_quadrature_points_and_weights(integral_type, cell, degree, rule, elements):
     """Create quadrature rule and return points and weights."""
     if integral_type == "cell":
-        return create_quadrature(cell.cellname(), degree, rule)
+        return create_quadrature(cell.cellname(), degree, rule, elements)
     elif integral_type in ufl.measure.facet_integral_types:
         facet_types = cell.facet_types()
         # 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}.")
-        return create_quadrature(facet_types[0].cellname(), degree, rule)
+        return create_quadrature(facet_types[0].cellname(), degree, rule, elements)
     elif integral_type in ufl.measure.point_integral_types:
-        return create_quadrature("vertex", degree, rule)
+        return create_quadrature("vertex", degree, rule, elements)
     elif integral_type == "expression":
         return (None, None)