Boundary Quadrature element

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,15 +60,18 @@ 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
         sub_el = ref_el.construct_subelement(dimension)
         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)

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)
 
 

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)