diff --git a/demo/ExpressionInterpolation.py b/demo/ExpressionInterpolation.py index ebccde1de..f16a4c9b8 100644 --- a/demo/ExpressionInterpolation.py +++ b/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) diff --git a/ffcx/codegeneration/basix_custom_element_template.py b/ffcx/codegeneration/basix_custom_element_template.py index f5fb55ca9..a2370329e 100644 --- a/ffcx/codegeneration/basix_custom_element_template.py +++ b/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} diff --git a/ffcx/codegeneration/finite_element.py b/ffcx/codegeneration/finite_element.py index 1a2fb7114..a5d51c835 100644 --- a/ffcx/codegeneration/finite_element.py +++ b/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 diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index 170961b33..63567e305 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/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 diff --git a/ffcx/element_interface.py b/ffcx/element_interface.py index 87722b798..140a296b5 100644 --- a/ffcx/element_interface.py +++ b/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.""" diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index f658cf8d3..0e9071154 100644 --- a/ffcx/ir/elementtables.py +++ b/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() diff --git a/ffcx/ir/representation.py b/ffcx/ir/representation.py index 3bdf746ff..5513bdc3e 100644 --- a/ffcx/ir/representation.py +++ b/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) diff --git a/ffcx/ir/representationutils.py b/ffcx/ir/representationutils.py index 49a99d6c4..4a2d1e3b7 100644 --- a/ffcx/ir/representationutils.py +++ b/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)