From e42f10fcc0b7cd8ce1a11ef0b637082eea58cea6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Dec 2024 13:40:40 +0000 Subject: [PATCH] Boundary Quadrature element --- FIAT/quadrature_schemes.py | 7 +++++-- finat/quadrature.py | 11 +++++++---- finat/quadrature_element.py | 5 +++-- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index d0c03ade6..1dcf21ff2 100644 --- a/FIAT/quadrature_schemes.py +++ b/FIAT/quadrature_schemes.py @@ -42,7 +42,7 @@ from FIAT.macro import MacroQuadratureRule -def create_quadrature(ref_el, degree, scheme="default", entity=None): +def create_quadrature(ref_el, degree, scheme="default", entity=None, trace=False): """ Generate quadrature rule for given reference element that will integrate a polynomial of order 'degree' exactly. @@ -60,6 +60,7 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None): "KMV" -> spectral lumped scheme for low degree (<=5 on triangles, <=3 on tetrahedra). :kwarg entity: A tuple of entity dimension and entity id specifying the integration domain. If not provided, the domain is the entire cell. + :kwarg trace: Are we constructing a quadrature on codimension 1 facets? """ if entity is not None: dimension, entity_id = entity @@ -67,8 +68,10 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None): Q_ref = create_quadrature(sub_el, degree, scheme=scheme) return FacetQuadratureRule(ref_el, dimension, entity_id, Q_ref) - if ref_el.is_macrocell(): + if ref_el.is_macrocell() or trace: dimension = ref_el.get_dimension() + if trace: + dimension = dimension - 1 sub_el = ref_el.construct_subelement(dimension) Q_ref = create_quadrature(sub_el, degree, scheme=scheme) return MacroQuadratureRule(ref_el, Q_ref) diff --git a/finat/quadrature.py b/finat/quadrature.py index ec6def127..796831e02 100644 --- a/finat/quadrature.py +++ b/finat/quadrature.py @@ -11,7 +11,7 @@ from finat.point_set import GaussLegendrePointSet, PointSet, TensorPointSet -def make_quadrature(ref_el, degree, scheme="default"): +def make_quadrature(ref_el, degree, scheme="default", trace=False): """ Generate quadrature rule for given reference element that will integrate an polynomial of order 'degree' exactly. @@ -24,8 +24,11 @@ def make_quadrature(ref_el, degree, scheme="default"): :arg ref_el: The FIAT cell to create the quadrature for. :arg degree: The degree of polynomial that the rule should integrate exactly. + :arg trace: Are we constructing a quadrature on codimension 1 facets? """ if ref_el.get_shape() == TENSORPRODUCT: + if trace: + raise NotImplementedError("Trace quadrature rule not implemented on tensor product cells") try: degree = tuple(degree) except TypeError: @@ -37,12 +40,12 @@ def make_quadrature(ref_el, degree, scheme="default"): return TensorProductQuadratureRule(quad_rules, ref_el=ref_el) if ref_el.get_shape() == QUADRILATERAL: - return make_quadrature(ref_el.product, degree, scheme) + return make_quadrature(ref_el.product, degree, scheme, trace) if degree < 0: raise ValueError("Need positive degree, not %d" % degree) - if ref_el.get_shape() == LINE and not ref_el.is_macrocell(): + if ref_el.get_shape() == LINE and not ref_el.is_macrocell() and not trace: # FIAT uses Gauss-Legendre line quadature, however, since we # symbolically label it as such, we wish not to risk attaching # the wrong label in case FIAT changes. So we explicitly ask @@ -52,7 +55,7 @@ def make_quadrature(ref_el, degree, scheme="default"): point_set = GaussLegendrePointSet(fiat_rule.get_points()) return QuadratureRule(point_set, fiat_rule.get_weights(), ref_el=ref_el, io_ornt_map_tuple=fiat_rule._intrinsic_orientation_permutation_map_tuple) - fiat_rule = fiat_scheme(ref_el, degree, scheme) + fiat_rule = fiat_scheme(ref_el, degree, scheme=scheme, trace=trace) return QuadratureRule(PointSet(fiat_rule.get_points()), fiat_rule.get_weights(), ref_el=ref_el, io_ornt_map_tuple=fiat_rule._intrinsic_orientation_permutation_map_tuple) diff --git a/finat/quadrature_element.py b/finat/quadrature_element.py index 3f17ec399..febfb854d 100644 --- a/finat/quadrature_element.py +++ b/finat/quadrature_element.py @@ -13,7 +13,7 @@ from finat.quadrature import make_quadrature, AbstractQuadratureRule -def make_quadrature_element(fiat_ref_cell, degree, scheme="default"): +def make_quadrature_element(fiat_ref_cell, degree, scheme="default", trace=False): """Construct a :class:`QuadratureElement` from a given a reference element, degree and scheme. @@ -23,9 +23,10 @@ def make_quadrature_element(fiat_ref_cell, degree, scheme="default"): integrate exactly. :param scheme: The quadrature scheme to use - e.g. "default", "canonical" or "KMV". + :param trace: Are we constructing a quadrature on codimension 1 facets? :returns: The appropriate :class:`QuadratureElement` """ - rule = make_quadrature(fiat_ref_cell, degree, scheme=scheme) + rule = make_quadrature(fiat_ref_cell, degree, scheme=scheme, trace=trace) return QuadratureElement(fiat_ref_cell, rule)