Add optional variant kwarg to FunctionSpace() (#3295)

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
 

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)

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

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

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))

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')

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)