diff --git a/firedrake/bcs.py b/firedrake/bcs.py index 9f340e4181..ae20ab0e2c 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -134,11 +134,19 @@ def nodes(self): raise NotImplementedError("Strong BCs not implemented for element %r, use Nitsche-type methods until we figure this out" % V.finat_element) def hermite_stride(bcnodes): - if isinstance(self._function_space.finat_element, finat.Hermite) and \ - self._function_space.mesh().topological_dimension() == 1: - return bcnodes[::2] # every second dof is the vertex value - else: - return bcnodes + fe = self._function_space.finat_element + tdim = self._function_space.mesh().topological_dimension() + if isinstance(fe, finat.Hermite) and tdim == 1: + bcnodes = bcnodes[::2] # every second dof is the vertex value + elif fe.complex.is_macrocell() and self._function_space.ufl_element().sobolev_space == ufl.H1: + # Skip derivative nodes for supersmooth H1 functions + nodes = fe.fiat_equivalent.dual_basis() + deriv_nodes = [i for i, node in enumerate(nodes) + if len(node.deriv_dict) != 0] + if len(deriv_nodes) > 0: + deriv_ids = self._function_space.cell_node_list[:, deriv_nodes] + bcnodes = np.setdiff1d(bcnodes, deriv_ids) + return bcnodes sub_d = (self.sub_domain, ) if isinstance(self.sub_domain, str) else as_tuple(self.sub_domain) sub_d = [s if isinstance(s, str) else as_tuple(s) for s in sub_d] diff --git a/firedrake/embedding.py b/firedrake/embedding.py index cfb3b3c771..404673641b 100644 --- a/firedrake/embedding.py +++ b/firedrake/embedding.py @@ -4,7 +4,7 @@ import ufl -def get_embedding_dg_element(element): +def get_embedding_dg_element(element, broken_cg=False): cell = element.cell degree = element.degree() family = lambda c: "DG" if c.is_simplex() else "DQ" @@ -16,6 +16,8 @@ def get_embedding_dg_element(element): for (c, d) in zip(cell.sub_cells(), degree))) else: scalar_element = finat.ufl.FiniteElement(family(cell), cell=cell, degree=degree) + if broken_cg: + scalar_element = finat.ufl.BrokenElement(scalar_element.reconstruct(family="Lagrange")) shape = element.value_shape if len(shape) == 0: DG = scalar_element diff --git a/firedrake/mg/embedded.py b/firedrake/mg/embedded.py index 3a7fe74f12..1ea27a6663 100644 --- a/firedrake/mg/embedded.py +++ b/firedrake/mg/embedded.py @@ -12,15 +12,22 @@ __all__ = ("TransferManager", ) -native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ"]) -alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced Hsieh-Clough-Tocher", "Johnson-Mercier"]) -non_native_variants = frozenset(["integral", "fdm"]) +native_families = frozenset(["Lagrange", "Discontinuous Lagrange", "Real", "Q", "DQ", "BrokenElement"]) +alfeld_families = frozenset(["Hsieh-Clough-Tocher", "Reduced-Hsieh-Clough-Tocher", "Johnson-Mercier", + "Alfeld-Sorokina", "Arnold-Qin", "Reduced-Arnold-Qin", "Christiansen-Hu", + "Guzman-Neilan", "Guzman-Neilan Bubble"]) +non_native_variants = frozenset(["integral", "fdm", "alfeld"]) def get_embedding_element(element): - dg_element = get_embedding_dg_element(element) - if element.family() in alfeld_families: - dg_element = dg_element.reconstruct(variant="alfeld") + broken_cg = element.sobolev_space in {ufl.H1, ufl.H2} + dg_element = get_embedding_dg_element(element, broken_cg=broken_cg) + variant = element.variant() or "default" + family = element.family() + # Elements on Alfeld splits are embedded onto DG Powell-Sabin. + # This yields supermesh projection + if (family in alfeld_families) or ("alfeld" in variant.lower() and family != "Discontinuous Lagrange"): + dg_element = dg_element.reconstruct(variant="powell-sabin") return dg_element @@ -67,7 +74,7 @@ def is_native(self, element): return True if isinstance(element.cell, ufl.TensorProductCell) and len(element.sub_elements) > 0: return reduce(and_, map(self.is_native, element.sub_elements)) - return element.family() in native_families and not element.variant() in non_native_variants + return (element.family() in native_families) and not (element.variant() in non_native_variants) def _native_transfer(self, element, op): try: diff --git a/tests/macro/test_macro_interp_project.py b/tests/macro/test_macro_interp_project.py index d71be2a673..12d175fdb2 100644 --- a/tests/macro/test_macro_interp_project.py +++ b/tests/macro/test_macro_interp_project.py @@ -1,34 +1,41 @@ import pytest +import numpy from firedrake import * def interp(u, f): u.interpolate(f) + return assemble(inner(u-f, u-f)*dx)**0.5 -def proj(u, f): - u.project(f) +def proj(u, f, bcs=None): + u.project(f, bcs=bcs) + return assemble(inner(u-f, u-f)*dx)**0.5 def proj_bc(u, f): - u.project(f, bcs=DirichletBC(u.function_space(), f, "on_boundary")) + return proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary")) def h1_proj(u, f, bcs=None): # compute h1 projection of f into u's function # space, store the result in u. v = TestFunction(u.function_space()) - F = (inner(grad(u-f), grad(v)) * dx + d = {H2: grad, H1: grad, HCurl: curl, HDiv: div, HDivDiv: div}[u.ufl_element().sobolev_space] + F = (inner(d(u-f), d(v)) * dx + inner(u-f, v) * dx) + fcp = {"mode": "vanilla"} solve(F == 0, u, bcs=bcs, solver_parameters={"snes_type": "ksponly", "ksp_type": "preonly", - "pc_type": "cholesky"}) + "pc_type": "cholesky"}, + form_compiler_parameters=fcp) + return assemble(F(u-f), form_compiler_parameters=fcp)**0.5 def h1_proj_bc(u, f): - h1_proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary")) + return h1_proj(u, f, bcs=DirichletBC(u.function_space(), f, "on_boundary")) @pytest.fixture(params=("square", "cube")) @@ -53,21 +60,120 @@ def test_projection_scalar_monomial(op, mesh, degree, variant): u = Function(V) x = SpatialCoordinate(mesh) f = sum(x) ** degree - op(u, f) - error = sqrt(assemble(inner(u - f, u - f) * dx)) + error = op(u, f) assert error < 1E-7 -@pytest.mark.parametrize(('el', 'accorder'), - [('HCT', 3), - ('HCT-red', 2)]) +@pytest.fixture +def hierarchy(request): + refine = 1 + mh2 = MeshHierarchy(UnitSquareMesh(5, 5), refine) + mh3 = MeshHierarchy(UnitCubeMesh(3, 3, 3), refine) + return {2: mh2, 3: mh3} + + +def mesh_sizes(mh): + mesh_size = [] + for msh in mh: + DG0 = FunctionSpace(msh, "DG", 0) + h = Function(DG0).interpolate(CellDiameter(msh)) + with h.dat.vec as hvec: + _, maxh = hvec.max() + mesh_size.append(maxh) + return mesh_size + + +def conv_rates(x, h): + x = numpy.asarray(x) + h = numpy.asarray(h) + return numpy.log2(x[:-1] / x[1:]) / numpy.log2(h[:-1] / h[1:]) + + +def run_convergence(mh, el, deg, convrate, op): + errors = [] + for msh in mh: + V = FunctionSpace(msh, el, deg) + u = Function(V) + x = SpatialCoordinate(msh) + f = sum(x) ** (deg+1) + if u.ufl_shape != (): + f = f * Constant(numpy.ones(u.ufl_shape)) + errors.append(op(u, f)) + + conv = conv_rates(errors, mesh_sizes(mh)) + assert numpy.all(conv > convrate - 0.25) + + +# Test L2/H1 convergence on C1 elements +@pytest.mark.parametrize(('dim', 'el', 'deg', 'convrate'), + [(2, 'PS6', 2, 2), + (2, 'PS12', 2, 2), + (2, 'HCT', 3, 3), + (2, 'HCT-red', 3, 2), + ]) @pytest.mark.parametrize('op', (proj, h1_proj)) -def test_projection_hct(el, accorder, op): - msh = UnitSquareMesh(1, 1) - V = FunctionSpace(msh, el, 3) +def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op): + if op == proj: + convrate += 1 + run_convergence(hierarchy[dim], el, deg, convrate, op) + + +# Test L2/H1 convergence on Stokes elements +@pytest.mark.parametrize(('dim', 'el', 'deg', 'convrate'), + [(2, 'Alfeld-Sorokina', 2, 2), + (3, 'Alfeld-Sorokina', 2, 2), + (2, 'Reduced-Arnold-Qin', 2, 1), + (3, 'Guzman-Neilan', 3, 1), + (3, 'Christiansen-Hu', 1, 1), + (2, 'Bernardi-Raugel', 2, 1), + (3, 'Bernardi-Raugel', 3, 1), + (2, 'Johnson-Mercier', 1, 1), + (3, 'Johnson-Mercier', 1, 1), + ]) +@pytest.mark.parametrize('op', (proj, h1_proj)) +def test_piola_convergence(hierarchy, dim, el, deg, convrate, op): + if op == proj: + convrate += 1 + run_convergence(hierarchy[dim], el, deg, convrate, op) + + +# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions +@pytest.mark.parametrize('enriched', (False, True)) +def test_supersmooth_bcs(mesh, enriched): + tdim = mesh.topological_dimension() + if enriched: + cell = mesh.ufl_cell() + element = EnrichedElement(FiniteElement("AS", cell=cell, degree=2), + FiniteElement("GNB", cell=cell, degree=tdim)) + V = FunctionSpace(mesh, element) + else: + V = FunctionSpace(mesh, "Alfeld-Sorokina", 2) + + # check that V in H1 + assert V.ufl_element().sobolev_space == H1 + + # check that V is supersmooth + nodes = V.finat_element.fiat_equivalent.dual.nodes + deriv_nodes = [i for i, node in enumerate(nodes) if len(node.deriv_dict)] + assert len(deriv_nodes) == tdim + 1 + + deriv_ids = V.cell_node_list[:, deriv_nodes] u = Function(V) - x = SpatialCoordinate(msh) - f = sum(x) ** accorder - op(u, f) - error = sqrt(assemble(inner(u - f, u - f) * dx)) - assert error < 1E-9 + + CG = FunctionSpace(mesh, "Lagrange", 2) + RT = FunctionSpace(mesh, "RT", 1) + for sub in [1, (1, 2), "on_boundary"]: + bc = DirichletBC(V, 0, sub) + + # check that we have the expected number of bc nodes + nnodes = len(bc.nodes) + expected = tdim * len(DirichletBC(CG, 0, sub).nodes) + if enriched: + expected += len(DirichletBC(RT, 0, sub).nodes) + assert nnodes == expected + + # check that the bc does not set the derivative nodes + u.assign(111) + u.dat.data_wo[deriv_ids] = 42 + bc.zero(u) + assert numpy.allclose(u.dat.data_ro[deriv_ids], 42) diff --git a/tests/macro/test_macro_multigrid.py b/tests/macro/test_macro_multigrid.py index 5239cd3515..0c91119383 100644 --- a/tests/macro/test_macro_multigrid.py +++ b/tests/macro/test_macro_multigrid.py @@ -152,7 +152,12 @@ def test_macro_multigrid_poisson(hierarchy, degree, variant): uh = Function(V) problem = LinearVariationalProblem(a, L, uh, bcs=bcs) solver = LinearVariationalSolver(problem, solver_parameters=mg_params) - solver.solve() + if complex_mode and variant == "alfeld": + with pytest.raises(NotImplementedError): + solver.solve() + else: + solver.solve() + expected = 10 if mesh.geometric_dimension() == 3 and variant == "alfeld": expected = 14 diff --git a/tests/macro/test_macro_quadrature.py b/tests/macro/test_macro_quadrature.py index 59c9a30e68..7330965d4f 100644 --- a/tests/macro/test_macro_quadrature.py +++ b/tests/macro/test_macro_quadrature.py @@ -7,7 +7,7 @@ def alfeld_split(msh): dim = msh.geometric_dimension() coords = msh.coordinates.dat.data.reshape((-1, dim)) - coords = numpy.row_stack((coords, numpy.average(coords, 0))) + coords = numpy.vstack((coords, numpy.average(coords, 0))) cells = [list(map(lambda i: dim+1 if i == j else i, range(dim+1))) for j in range(dim+1)] plex = plex_from_cell_list(dim, cells, coords, msh.comm) return Mesh(plex) @@ -68,8 +68,8 @@ def test_macro_quadrature_piecewise(degree, variant, meshes): c = Constant(numpy.arange(1, gdim+1)) expr = abs(dot(c, x - x0)) ** degree elif variant == "iso": - vecs = list(map(Constant, numpy.row_stack([numpy.eye(gdim), - numpy.ones((max(degree-gdim, 0), gdim))]))) + vecs = list(map(Constant, numpy.vstack([numpy.eye(gdim), + numpy.ones((max(degree-gdim, 0), gdim))]))) expr = numpy.prod([abs(dot(c, x)) for c in vecs[:degree]]) else: raise ValueError("Unexpected variant") diff --git a/tests/macro/test_macro_solve.py b/tests/macro/test_macro_solve.py index 4e8b75a401..8a58807d5b 100644 --- a/tests/macro/test_macro_solve.py +++ b/tests/macro/test_macro_solve.py @@ -33,23 +33,35 @@ def mixed_element(mh, variant): return Vel, Pel -def conv_rates(x): +def mesh_sizes(mh): + mesh_size = [] + for msh in mh: + DG0 = FunctionSpace(msh, "DG", 0) + h = Function(DG0).interpolate(CellDiameter(msh)) + with h.dat.vec as hvec: + _, maxh = hvec.max() + mesh_size.append(maxh) + return mesh_size + + +def conv_rates(x, h): x = np.asarray(x) - return np.log2(x[:-1] / x[1:]) + h = np.asarray(h) + return np.log2(x[:-1] / x[1:]) / np.log2(h[:-1] / h[1:]) @pytest.fixture def convergence_test(variant): if variant == "iso": - def check(uerr, perr): - return (conv_rates(uerr)[-1] >= 1.9 + def check(uerr, perr, h): + return (conv_rates(uerr, h)[-1] >= 1.9 and np.allclose(perr, 0, atol=1.e-8)) elif variant == "alfeld": - def check(uerr, perr): + def check(uerr, perr, h): return (np.allclose(uerr, 0, atol=5.e-9) and np.allclose(perr, 0, atol=5.e-7)) elif variant == "th": - def check(uerr, perr): + def check(uerr, perr, h): return (np.allclose(uerr, 0, atol=1.e-10) and np.allclose(perr, 0, atol=1.e-8)) return check @@ -100,7 +112,7 @@ def test_riesz(mh, variant, mixed_element, convergence_test): u_err.append(errornorm(as_vector(zexact[:dim]), uh)) p_err.append(errornorm(zexact[-1], ph)) - assert convergence_test(u_err, p_err) + assert convergence_test(u_err, p_err, mesh_sizes(mh)) def stokes_mms(Z, zexact): @@ -128,7 +140,6 @@ def test_stokes(mh, variant, mixed_element, convergence_test): dim = mh[0].geometric_dimension() if variant == "iso" and dim == 3: pytest.xfail("P2:P1 iso x P1 is not inf-sup stable in 3D") - u_err = [] p_err = [] el1, el2 = mixed_element @@ -155,7 +166,7 @@ def test_stokes(mh, variant, mixed_element, convergence_test): u_err.append(errornorm(as_vector(zexact[:dim]), uh)) p_err.append(errornormL2_0(zexact[-1], ph)) - assert convergence_test(u_err, p_err) + assert convergence_test(u_err, p_err, mesh_sizes(mh)) def test_div_free(mh, variant, mixed_element, div_test): diff --git a/tests/regression/test_helmholtz_zany.py b/tests/regression/test_helmholtz_zany.py index c0b3bec3d6..32906353c0 100644 --- a/tests/regression/test_helmholtz_zany.py +++ b/tests/regression/test_helmholtz_zany.py @@ -54,11 +54,7 @@ def helmholtz(n, el_type, degree, perturb): # It is, somewhat oddly, a suitable C^1 nonconforming element but # not a suitable C^0 nonconforming one. @pytest.mark.parametrize(('el', 'deg', 'convrate'), - [('PS6', 2, 2), - ('PS12', 2, 2), - ('Hermite', 3, 3), - ('HCT', 3, 3), - ('HCT', 4, 4), + [('Hermite', 3, 3), ('Bell', 5, 4), ('Argyris', 5, 5), ('Argyris', 6, 6)])