Boundary Quadrature elements (#3973)

* Extract quad degree from Quadrature elements

* Projection: support for Boundary Quadrature element

* Projector: refactoring

Author: pbrubeck

firedrake/output/vtk_output.py

@@ -81,7 +81,7 @@ def is_dg(V):
 
     :arg V: A FunctionSpace.
     """
-    return V.finat_element.entity_dofs() == V.finat_element.entity_closure_dofs()
+    return V.finat_element.is_dg()
 
 
 def is_linear(V):

firedrake/preconditioners/fdm.py

@@ -720,10 +720,7 @@ def _element_kernels(self):
 
     @cached_property
     def insert_mode(self):
-        is_dg = {}
-        for Vsub in self.V:
-            element = Vsub.finat_element
-            is_dg[Vsub] = element.entity_dofs() == element.entity_closure_dofs()
+        is_dg = {Vsub: Vsub.finat_element.is_dg() for Vsub in self.V}
 
         insert_mode = {}
         for Vrow, Vcol in product(self.V, self.V):
@@ -1935,9 +1932,7 @@ def assemble_reference_tensor(self, V):
 
         degree = max(e.degree() for e in line_elements)
         eta = float(self.appctx.get("eta", degree*(degree+1)))
-        element = V.finat_element
-        is_dg = element.entity_dofs() == element.entity_closure_dofs()
-
+        is_dg = V.finat_element.is_dg()
         Afdm = []  # sparse interval mass and stiffness matrices for each direction
         Dfdm = []  # tabulation of normal derivatives at the boundary for each direction
         bdof = []  # indices of point evaluation dofs for each direction

firedrake/projection.py

@@ -12,7 +12,6 @@
 from firedrake import functionspaceimpl
 from firedrake import function
 from firedrake.adjoint_utils import annotate_project
-from finat import HDivTrace
 
 
 __all__ = ['project', 'Projector']
@@ -164,24 +163,25 @@ def __init__(
         self.form_compiler_parameters = form_compiler_parameters
         self.bcs = bcs
         self.constant_jacobian = constant_jacobian
-        try:
-            element = self.target.function_space().finat_element
-            is_dg = element.entity_dofs() == element.entity_closure_dofs()
-            is_variable_layers = self.target.function_space().mesh().variable_layers
-        except AttributeError:
-            # Mixed space
-            is_dg = False
-            is_variable_layers = True
-        self.use_slate_for_inverse = (use_slate_for_inverse and is_dg and not is_variable_layers
-                                      and (not complex_mode or SLATE_SUPPORTS_COMPLEX))
+
+        F = self.target.function_space()
+        if type(F.ufl_element()) is finat.ufl.MixedElement:
+            slate_supported = False
+            needs_trace = False
+        else:
+            slate_supported = (F.finat_element.is_dg() and not F.mesh().variable_layers
+                               and (not complex_mode or SLATE_SUPPORTS_COMPLEX))
+            needs_trace = F.ufl_element().family() in {"HDiv Trace", "Boundary Quadrature"}
+
+        self.use_slate_for_inverse = use_slate_for_inverse and slate_supported
+        self.needs_trace = needs_trace
 
     @cached_property
     def A(self):
-        u = firedrake.TrialFunction(self.target.function_space())
-        v = firedrake.TestFunction(self.target.function_space())
         F = self.target.function_space()
-        mixed = isinstance(F.ufl_element(), finat.ufl.MixedElement)
-        if not mixed and isinstance(F.finat_element, HDivTrace):
+        u = firedrake.TrialFunction(F)
+        v = firedrake.TestFunction(F)
+        if self.needs_trace:
             if F.extruded:
                 a = (
                     firedrake.inner(u, v)*firedrake.ds_t
@@ -247,12 +247,11 @@ class BasicProjector(ProjectorBase):
     def rhs_form(self):
         v = firedrake.TestFunction(self.target.function_space())
         F = self.target.function_space()
-        mixed = isinstance(F.ufl_element(), finat.ufl.MixedElement)
-        if not mixed and isinstance(F.finat_element, HDivTrace):
+        if self.needs_trace:
             # Project onto a trace space by supplying the respective form on the facets.
             # The measures on the facets differ between extruded and non-extruded mesh.
             # FIXME The restrictions of cg onto the facets is also a trace space,
-            # but we only cover DGT.
+            # but we only cover DGT and BQ.
             if F.extruded:
                 form = (
                     firedrake.inner(self.source, v)*firedrake.ds_t

tests/firedrake/regression/test_projection.py

@@ -391,3 +391,32 @@ def run_extr_trace_projection(x, degree=1, family='DGT'):
                          + area * (w - ref) * (w - ref) * ds_b
                          + area * (w('+') - ref('+')) * (w('+') - ref('+')) * dS_h
                          + area * (w('+') - ref('+')) * (w('+') - ref('+')) * dS_v))
+
+
+def test_project_scalar_boundary_quadrature():
+    msh = UnitSquareMesh(4, 4)
+    x = SpatialCoordinate(msh)
+    f = x[0]*(2-x[0])*x[1]*(2-x[1])
+
+    T = FunctionSpace(msh, "BQ", 8)
+    w = Function(T)
+    w.project(f, solver_parameters={'mat_type': 'matfree', 'ksp_type': 'preonly', 'pc_type': 'jacobi'})
+
+    err = w - f
+    err2 = inner(err, err)
+    error = sqrt(assemble(err2 * ds + err2('+') * dS))
+    assert error < 1E-8
+
+
+def test_project_vector_boundary_quadrature():
+    msh = UnitCubeMesh(4, 4, 4)
+    n = FacetNormal(msh)
+
+    T = VectorFunctionSpace(msh, "BQ", 0)
+    w = Function(T)
+    w.project(n, solver_parameters={'mat_type': 'matfree', 'ksp_type': 'preonly', 'pc_type': 'jacobi'})
+
+    err = w - n
+    err2 = inner(err, err)
+    error = sqrt(assemble(err2 * ds))
+    assert error < 1E-8

tests/firedrake/regression/test_projection_symmetric_tensor.py

@@ -32,7 +32,7 @@ def test_projection_symmetric_tensor(mesh, degree, family, tdim):
         sp = {"mat_type": "matfree",
               "ksp_type": "preonly",
               "pc_type": "jacobi"}
-        fcp = {"quadrature_degree": Nq}
+        fcp = {}
     else:
         Q = TensorFunctionSpace(mesh, family, degree=degree, shape=shape, symmetry=None)
         Qs = TensorFunctionSpace(mesh, family, degree=degree, shape=shape, symmetry=True)

tests/firedrake/regression/test_quadrature.py

@@ -1,8 +1,13 @@
 from firedrake import *
+import pytest
 
 
-def test_hand_specified_quadrature():
-    mesh = UnitSquareMesh(5, 5)
+@pytest.fixture
+def mesh(request):
+    return UnitSquareMesh(5, 5)
+
+
+def test_hand_specified_quadrature(mesh):
     V = FunctionSpace(mesh, 'CG', 2)
     v = TestFunction(V)
 
@@ -12,3 +17,18 @@ def test_hand_specified_quadrature():
     a_q2 = assemble(a, form_compiler_parameters={'quadrature_degree': 2})
 
     assert not np.allclose(a_q0.dat.data, a_q2.dat.data)
+
+
+@pytest.mark.parametrize("diagonal", [False, True])
+@pytest.mark.parametrize("mat_type", ["matfree", "aij"])
+@pytest.mark.parametrize("family", ["Quadrature", "Boundary Quadrature"])
+def test_quadrature_element(mesh, family, mat_type, diagonal):
+    V = FunctionSpace(mesh, family, 2)
+    v = TestFunction(V)
+    u = TrialFunction(V)
+    if family == "Boundary Quadrature":
+        a = inner(u, v) * ds + inner(u('+'), v('+')) * dS
+    else:
+        a = inner(u, v) * dx
+
+    assemble(a, mat_type=mat_type, diagonal=diagonal)

tsfc/driver.py

@@ -95,10 +95,10 @@ def compile_integral(integral_data, form_data, prefix, parameters, interface, *,
         interface = firedrake_interface_loopy.KernelBuilder
     scalar_type = parameters["scalar_type"]
     integral_type = integral_data.integral_type
-    if integral_type.startswith("interior_facet") and diagonal:
-        raise NotImplementedError("Sorry, we can't assemble the diagonal of a form for interior facet integrals")
     mesh = integral_data.domain
     arguments = form_data.preprocessed_form.arguments()
+    if integral_type.startswith("interior_facet") and diagonal and any(a.function_space().finat_element.is_dg() for a in arguments):
+        raise NotImplementedError("Sorry, we can't assemble the diagonal of a form for interior facet integrals")
     kernel_name = f"{prefix}_{integral_type}_integral"
     # Dict mapping domains to index in original_form.ufl_domains()
     domain_numbering = form_data.original_form.domain_numbering()

tsfc/fem.py

@@ -641,15 +641,13 @@ def fiat_to_ufl(fiat_dict, order):
 
 @translate.register(Argument)
 def translate_argument(terminal, mt, ctx):
-    argument_multiindex = ctx.argument_multiindices[terminal.number()]
-    sigma = tuple(gem.Index(extent=d) for d in mt.expr.ufl_shape)
     element = ctx.create_element(terminal.ufl_element(), restriction=mt.restriction)
 
     def callback(entity_id):
         finat_dict = ctx.basis_evaluation(element, mt, entity_id)
         # Filter out irrelevant derivatives
-        filtered_dict = {alpha: table
-                         for alpha, table in finat_dict.items()
+        filtered_dict = {alpha: finat_dict[alpha]
+                         for alpha in finat_dict
                          if sum(alpha) == mt.local_derivatives}
 
         # Change from FIAT to UFL arrangement
@@ -664,7 +662,8 @@ def callback(entity_id):
         quad_multiindex_permuted = _make_quad_multiindex_permuted(mt, ctx)
         mapper = gem.node.MemoizerArg(gem.optimise.filtered_replace_indices)
         table = mapper(table, tuple(zip(quad_multiindex, quad_multiindex_permuted)))
-    return gem.ComponentTensor(gem.Indexed(table, argument_multiindex + sigma), sigma)
+    argument_multiindex = ctx.argument_multiindices[terminal.number()]
+    return gem.partial_indexed(table, argument_multiindex)
 
 
 @translate.register(TSFCConstantMixin)

tsfc/kernel_interface/common.py

@@ -295,21 +295,28 @@ def create_context(self):
 
 
 def set_quad_rule(params, cell, integral_type, functions):
-    # Check if the integral has a quad degree attached, otherwise use
-    # the estimated polynomial degree attached by compute_form_data
+    # Check if the integral has a quad degree or quad element attached,
+    # otherwise use the estimated polynomial degree attached by compute_form_data
+    quad_rule = params.get("quadrature_rule", "default")
     try:
         quadrature_degree = params["quadrature_degree"]
     except KeyError:
-        quadrature_degree = params["estimated_polynomial_degree"]
-        function_degrees = [f.ufl_function_space().ufl_element().degree()
-                            for f in functions]
-        if all((asarray(quadrature_degree) > 10 * asarray(degree)).all()
-               for degree in function_degrees):
-            logger.warning("Estimated quadrature degree %s more "
-                           "than tenfold greater than any "
-                           "argument/coefficient degree (max %s)",
-                           quadrature_degree, max_degree(function_degrees))
-    quad_rule = params.get("quadrature_rule", "default")
+        elements = [f.ufl_function_space().ufl_element() for f in functions]
+        quad_data = set((e.degree(), e.quadrature_scheme() or "default") for e in elements
+                        if e.family() in {"Quadrature", "Boundary Quadrature"})
+        if len(quad_data) == 0:
+            quadrature_degree = params["estimated_polynomial_degree"]
+            if all((asarray(quadrature_degree) > 10 * asarray(e.degree())).all() for e in elements):
+                logger.warning("Estimated quadrature degree %s more "
+                               "than tenfold greater than any "
+                               "argument/coefficient degree (max %s)",
+                               quadrature_degree, max_degree([e.degree() for e in elements]))
+        else:
+            try:
+                (quadrature_degree, quad_rule), = quad_data
+            except ValueError:
+                raise ValueError("The quadrature rule cannot be inferred from multiple Quadrature elements")
+
     if isinstance(quad_rule, str):
         scheme = quad_rule
         fiat_cell = as_fiat_cell(cell)