diff --git a/docs/source/variational-problems.rst b/docs/source/variational-problems.rst index 9a8a7645e2..632e14346a 100644 --- a/docs/source/variational-problems.rst +++ b/docs/source/variational-problems.rst @@ -288,9 +288,9 @@ For discontinuous elements these are the Gauss-Legendre points, and for continuous elements these are the Gauss-Lobatto-Legendre points. For CG and DG spaces on simplices, Firedrake offers both equispaced points and the better conditioned recursive Legendre points from :cite:`Isaac2020` via the -`recursivenodes`_ module. These are selected by passing -`variant="equispaced"` or `variant="spectral"` to the -:py:class:`~ufl.classes.FiniteElement` constructor. For example: +`recursivenodes`_ module. These are selected by passing `variant="equispaced"` +or `variant="spectral"` to the :py:class:`~ufl.classes.FiniteElement` or +:py:func:`~.FunctionSpace` constructors. For example: .. code-block:: python3 diff --git a/firedrake/functionspace.py b/firedrake/functionspace.py index e79c0c2275..8345a0eb7e 100644 --- a/firedrake/functionspace.py +++ b/firedrake/functionspace.py @@ -20,7 +20,7 @@ @PETSc.Log.EventDecorator() -def make_scalar_element(mesh, family, degree, vfamily, vdegree): +def make_scalar_element(mesh, family, degree, vfamily, vdegree, variant): """Build a scalar :class:`finat.ufl.finiteelement.FiniteElement`. Parameters @@ -31,6 +31,8 @@ def make_scalar_element(mesh, family, degree, vfamily, vdegree): The finite element family. degree : The degree of the finite element. + variant : + The variant of the finite element. vfamily : The finite element in the vertical dimension (extruded meshes only). @@ -57,20 +59,20 @@ def make_scalar_element(mesh, family, degree, vfamily, vdegree): and vfamily is not None and vdegree is not None: la = finat.ufl.FiniteElement(family, cell=cell.sub_cells()[0], - degree=degree) + degree=degree, variant=variant) # If second element was passed in, use it lb = finat.ufl.FiniteElement(vfamily, cell=ufl.interval, - degree=vdegree) + degree=vdegree, variant=variant) # Now make the TensorProductElement return finat.ufl.TensorProductElement(la, lb) else: - return finat.ufl.FiniteElement(family, cell=cell, degree=degree) + return finat.ufl.FiniteElement(family, cell=cell, degree=degree, variant=variant) @PETSc.Log.EventDecorator("CreateFunctionSpace") -def FunctionSpace(mesh, family, degree=None, name=None, vfamily=None, - vdegree=None): +def FunctionSpace(mesh, family, degree=None, name=None, + vfamily=None, vdegree=None, variant=None): """Create a :class:`.FunctionSpace`. Parameters @@ -89,6 +91,8 @@ def FunctionSpace(mesh, family, degree=None, name=None, vfamily=None, vdegree : The degree of the element in the vertical dimension (extruded meshes only). + variant : + The variant of the finite element. Notes ----- @@ -97,13 +101,13 @@ def FunctionSpace(mesh, family, degree=None, name=None, vfamily=None, are ignored and the appropriate :class:`.FunctionSpace` is returned. """ - element = make_scalar_element(mesh, family, degree, vfamily, vdegree) + element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) return impl.WithGeometry.make_function_space(mesh, element, name=name) @PETSc.Log.EventDecorator() -def DualSpace(mesh, family, degree=None, name=None, vfamily=None, - vdegree=None): +def DualSpace(mesh, family, degree=None, name=None, + vfamily=None, vdegree=None, variant=None): """Create a :class:`.FunctionSpace`. Parameters @@ -122,6 +126,8 @@ def DualSpace(mesh, family, degree=None, name=None, vfamily=None, vdegree : The degree of the element in the vertical dimension (extruded meshes only). + variant : + The variant of the finite element. Notes ----- @@ -130,13 +136,13 @@ def DualSpace(mesh, family, degree=None, name=None, vfamily=None, other arguments are ignored and the appropriate :class:`.FunctionSpace` is returned. """ - element = make_scalar_element(mesh, family, degree, vfamily, vdegree) + element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) return impl.FiredrakeDualSpace.make_function_space(mesh, element, name=name) @PETSc.Log.EventDecorator() def VectorFunctionSpace(mesh, family, degree=None, dim=None, - name=None, vfamily=None, vdegree=None): + name=None, vfamily=None, vdegree=None, variant=None): """Create a rank-1 :class:`.FunctionSpace`. Parameters @@ -158,6 +164,8 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, vdegree : The degree of the element in the vertical dimension (extruded meshes only). + variant : + The variant of the finite element. Notes ----- @@ -174,7 +182,7 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, pass it to :class:`.FunctionSpace` directly instead. """ - sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree) + sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) if dim is None: dim = mesh.ufl_cell().geometric_dimension() if not isinstance(dim, numbers.Integral) and dim > 0: @@ -186,7 +194,7 @@ def VectorFunctionSpace(mesh, family, degree=None, dim=None, @PETSc.Log.EventDecorator() def TensorFunctionSpace(mesh, family, degree=None, shape=None, symmetry=None, name=None, vfamily=None, - vdegree=None): + vdegree=None, variant=None): """Create a rank-2 FunctionSpace. Parameters @@ -211,6 +219,8 @@ def TensorFunctionSpace(mesh, family, degree=None, shape=None, vdegree : The degree of the element in the vertical dimension (extruded meshes only). + variant : + The variant of the finite element. Notes ----- @@ -226,7 +236,7 @@ def TensorFunctionSpace(mesh, family, degree=None, shape=None, `FunctionSpace` directly instead. """ - sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree) + sub_element = make_scalar_element(mesh, family, degree, vfamily, vdegree, variant) shape = shape or (mesh.ufl_cell().geometric_dimension(),) * 2 element = finat.ufl.TensorElement(sub_element, shape=shape, symmetry=symmetry) return FunctionSpace(mesh, element, name=name) diff --git a/tests/regression/test_fdm.py b/tests/regression/test_fdm.py index 1fd6112000..58fc7782de 100644 --- a/tests/regression/test_fdm.py +++ b/tests/regression/test_fdm.py @@ -137,8 +137,7 @@ def test_p_independence_hgrad(mesh, variant): expected = [16, 12] if mesh.topological_dimension() == 3 else [9, 7] solvers = [fdmstar] if variant is None else [fdmstar, facetstar] for degree in range(3, 6): - element = FiniteElement(family, cell=mesh.ufl_cell(), degree=degree, variant=variant) - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, family, degree, variant=variant) problem = build_riesz_map(V, grad) for sp, expected_it in zip(solvers, expected): assert solve_riesz_map(problem, sp) <= expected_it @@ -150,8 +149,7 @@ def test_p_independence_hcurl(mesh): expected = [13, 10] if mesh.topological_dimension() == 3 else [6, 6] solvers = [fdmstar, facetstar] for degree in range(3, 6): - element = FiniteElement(family, cell=mesh.ufl_cell(), degree=degree, variant="fdm") - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, family, degree, variant="fdm") problem = build_riesz_map(V, curl) for sp, expected_it in zip(solvers, expected): assert solve_riesz_map(problem, sp) <= expected_it @@ -163,8 +161,7 @@ def test_p_independence_hdiv(mesh): expected = [6, 6] solvers = [fdmstar, facetstar] for degree in range(3, 6): - element = FiniteElement(family, cell=mesh.ufl_cell(), degree=degree, variant="fdm") - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, family, degree, variant="fdm") problem = build_riesz_map(V, div) for sp, expected_it in zip(solvers, expected): assert solve_riesz_map(problem, sp) <= expected_it @@ -207,18 +204,17 @@ def test_variable_coefficient(mesh): def fs(request, mesh): degree = 3 tdim = mesh.topological_dimension() - cell = mesh.ufl_cell() element = request.param variant = "fdm_ipdg" if element == "rt": family = "RTCF" if tdim == 2 else "NCF" - return FunctionSpace(mesh, FiniteElement(family, cell, degree=degree, variant=variant)) + return FunctionSpace(mesh, family, degree, variant=variant) else: if tdim == 1: family = "DG" if element == "dg" else "CG" else: family = "DQ" if element == "dg" else "Q" - return VectorFunctionSpace(mesh, FiniteElement(family, cell, degree=degree, variant=variant), dim=5-tdim) + return VectorFunctionSpace(mesh, family, degree, dim=5-tdim, variant=variant) @pytest.mark.skipcomplex diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 79859f5768..220c02bfde 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -83,6 +83,11 @@ def test_function_space_different_family_differ(mesh): assert FunctionSpace(mesh, "CG", 1) != FunctionSpace(mesh, "DG", 1) +def test_function_space_different_variant_differ(mesh): + "FunctionSpaces defined with different element variants differ." + assert FunctionSpace(mesh, "CG", 3, variant="equispaced") != FunctionSpace(mesh, "CG", 3) + + def test_function_space_vector_function_space_differ(mesh): """A FunctionSpace and a VectorFunctionSpace defined with the same family and degree differ.""" @@ -106,6 +111,14 @@ def test_function_space_collapse(cg1): assert cg1 == cg1.collapse() +@pytest.mark.parametrize("space", + [FunctionSpace, VectorFunctionSpace, + TensorFunctionSpace]) +def test_function_space_variant(mesh, space): + element = FiniteElement("DG", degree=1, variant="equispaced") + assert space(mesh, element) == space(mesh, "DG", 1, variant="equispaced") + + @pytest.mark.parametrize("modifier", [BrokenElement, HDivElement, HCurlElement]) @@ -219,8 +232,8 @@ def test_reconstruct_mesh_degree(family, dual): @pytest.mark.parametrize("family", ["CG", "DG"]) def test_reconstruct_variant(family, dual): m1 = UnitIntervalMesh(1) - V1 = FunctionSpace(m1, FiniteElement(family, cell=m1.ufl_cell(), degree=4, variant="spectral")) - V2 = FunctionSpace(m1, FiniteElement(family, cell=m1.ufl_cell(), degree=4, variant="equispaced")) + V1 = FunctionSpace(m1, family, degree=4, variant="spectral") + V2 = FunctionSpace(m1, family, degree=4, variant="equispaced") if dual: V1, V2 = V1.dual(), V2.dual() assert V1.reconstruct(variant="equispaced") == V2 diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index 5a5f93b9df..89dacae95f 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -290,45 +290,12 @@ def test_interpolator_Pk(degree): @pytest.mark.parametrize("degree", range(1, 4)) -def test_interpolator_spectral(degree): - mesh = UnitSquareMesh(10, 10, quadrilateral=True) +@pytest.mark.parametrize("variant", ("spectral", "equispaced")) +@pytest.mark.parametrize("quads", (True, False), ids=("quads", "triangles")) +def test_interpolator(quads, degree, variant): + mesh = UnitSquareMesh(10, 10, quadrilateral=quads) x = SpatialCoordinate(mesh) - fe1 = FiniteElement("CG", mesh.ufl_cell(), degree, variant="spectral") - P1 = FunctionSpace(mesh, fe1) - P2 = FunctionSpace(mesh, "CG", degree + 1) - - expr = x[0]**degree + x[1]**degree - x_P1 = assemble(interpolate(expr, P1)) - interpolator = Interpolator(TestFunction(P1), P2) - x_P2 = assemble(interpolator.interpolate(x_P1)) - x_P2_direct = assemble(interpolate(expr, P2)) - - assert np.allclose(x_P2.dat.data, x_P2_direct.dat.data) - - -@pytest.mark.parametrize("degree", range(1, 4)) -def test_interpolator_equiquads(degree): - mesh = UnitSquareMesh(10, 10, quadrilateral=True) - x = SpatialCoordinate(mesh) - fe1 = FiniteElement("CG", mesh.ufl_cell(), degree, variant="equispaced") - P1 = FunctionSpace(mesh, fe1) - P2 = FunctionSpace(mesh, "CG", degree + 1) - - expr = x[0]**degree + x[1]**degree - x_P1 = assemble(interpolate(expr, P1)) - interpolator = Interpolator(TestFunction(P1), P2) - x_P2 = assemble(interpolator.interpolate(x_P1)) - x_P2_direct = assemble(interpolate(expr, P2)) - - assert np.allclose(x_P2.dat.data, x_P2_direct.dat.data) - - -@pytest.mark.parametrize("degree", range(1, 4)) -def test_interpolator_equitris(degree): - mesh = UnitSquareMesh(10, 10, quadrilateral=False) - x = SpatialCoordinate(mesh) - fe1 = FiniteElement("CG", mesh.ufl_cell(), degree, variant="equispaced") - P1 = FunctionSpace(mesh, fe1) + P1 = FunctionSpace(mesh, "CG", degree, variant=variant) P2 = FunctionSpace(mesh, "CG", degree + 1) expr = x[0]**degree + x[1]**degree @@ -400,8 +367,7 @@ def test_adjoint_Pk(degree): def test_adjoint_quads(): mesh = UnitSquareMesh(10, 10) - fe1 = FiniteElement("CG", mesh.ufl_cell(), 1, variant="equispaced") - P1 = FunctionSpace(mesh, fe1) + P1 = FunctionSpace(mesh, "CG", 1) P2 = FunctionSpace(mesh, "CG", 2) v = conj(TestFunction(P2)) diff --git a/tests/regression/test_interpolation_nodes.py b/tests/regression/test_interpolation_nodes.py index 24ca4de2ca..6c6e338afb 100644 --- a/tests/regression/test_interpolation_nodes.py +++ b/tests/regression/test_interpolation_nodes.py @@ -34,10 +34,9 @@ def degree(request): ("N2div")], ids=lambda x: "%s" % x) def V(request, mesh, degree): - space = request.param + family = request.param over_integration = max(0, 9 - degree) - V_el = FiniteElement(space, mesh.ufl_cell(), degree, variant=f"integral({over_integration})") - return FunctionSpace(mesh, V_el) + return FunctionSpace(mesh, family, degree, variant=f"integral({over_integration})") def test_div_curl_preserving(V): @@ -48,18 +47,18 @@ def test_div_curl_preserving(V): elif dim == 3: x, y, z = SpatialCoordinate(mesh) if dim == 2: - if "Nedelec" in V.ufl_element().family(): + if V.ufl_element().sobolev_space == HCurl: expression = grad(sin(x)*cos(y)) else: expression = as_vector([cos(y), sin(x)]) elif dim == 3: - if "Nedelec" in V.ufl_element().family(): + if V.ufl_element().sobolev_space == HCurl: expression = grad(sin(x)*exp(y)*z) else: expression = as_vector([sin(y)*z, cos(x)*z, exp(x)]) f = assemble(interpolate(expression, V)) - if "Nedelec" in V.ufl_element().family(): + if V.ufl_element().sobolev_space == HCurl: norm_exp = sqrt(assemble(inner(curl(f), curl(f))*dx)) else: norm_exp = sqrt(assemble(inner(div(f), div(f))*dx)) @@ -80,11 +79,10 @@ def compute_interpolation_error(baseMesh, nref, space, degree): expression = as_vector([sin(x)*cos(y), exp(x)*y]) elif dim == 3: expression = as_vector([sin(y)*z*cos(x), cos(x)*z*x, exp(x)*y]) - V_el = FiniteElement(space, mesh.ufl_cell(), degree, variant="integral") - V = FunctionSpace(mesh, V_el) + V = FunctionSpace(mesh, space, degree, variant="integral") f = assemble(interpolate(expression, V)) error_l2 = errornorm(expression, f, 'L2') - if "Nedelec" in V.ufl_element().family(): + if V.ufl_element().sobolev_space == HCurl: error_hD = errornorm(expression, f, 'hcurl') else: error_hD = errornorm(expression, f, 'hdiv') diff --git a/tests/regression/test_variants.py b/tests/regression/test_variants.py index 6b53160aeb..0f6e052ba2 100644 --- a/tests/regression/test_variants.py +++ b/tests/regression/test_variants.py @@ -24,9 +24,7 @@ def mesh(request): @pytest.mark.parametrize('degree', [0, 2, 4]) def test_dg_integral_orthogonality(mesh, degree): - element = FiniteElement("DG", cell=mesh.ufl_cell(), degree=degree, variant="integral") - - V = FunctionSpace(mesh, element) + V = FunctionSpace(mesh, "DG", degree, variant="integral") x = SpatialCoordinate(mesh) x2 = dot(x, x) expr = exp(-x2)