diff --git a/firedrake/output/vtk_output.py b/firedrake/output/vtk_output.py index 23072c5019..7f71a16cbe 100644 --- a/firedrake/output/vtk_output.py +++ b/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): diff --git a/firedrake/preconditioners/fdm.py b/firedrake/preconditioners/fdm.py index d0039a3139..1b6aa9c9aa 100644 --- a/firedrake/preconditioners/fdm.py +++ b/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 diff --git a/firedrake/projection.py b/firedrake/projection.py index a1fdeeb6bd..3ce346f664 100644 --- a/firedrake/projection.py +++ b/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 diff --git a/tests/firedrake/regression/test_projection.py b/tests/firedrake/regression/test_projection.py index 176b56c93f..1424f5ffde 100644 --- a/tests/firedrake/regression/test_projection.py +++ b/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 diff --git a/tests/firedrake/regression/test_projection_symmetric_tensor.py b/tests/firedrake/regression/test_projection_symmetric_tensor.py index bc0e551b73..ebec1c9657 100644 --- a/tests/firedrake/regression/test_projection_symmetric_tensor.py +++ b/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) diff --git a/tests/firedrake/regression/test_quadrature.py b/tests/firedrake/regression/test_quadrature.py index 6ec5f8b855..9c738b888b 100644 --- a/tests/firedrake/regression/test_quadrature.py +++ b/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) diff --git a/tsfc/driver.py b/tsfc/driver.py index 38a780a0ad..9e518292bc 100644 --- a/tsfc/driver.py +++ b/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() diff --git a/tsfc/fem.py b/tsfc/fem.py index 10c10dc97d..d86be550db 100644 --- a/tsfc/fem.py +++ b/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) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py index 3f7daa5383..667dab75a8 100644 --- a/tsfc/kernel_interface/common.py +++ b/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)