Mimetic Spectral Elements

requirements-git.txt

@@ -1,5 +1,5 @@
 git+https://github.com/coneoproject/COFFEE.git#egg=coffee
-git+https://github.com/firedrakeproject/ufl.git#egg=ufl
-git+https://github.com/firedrakeproject/fiat.git#egg=fiat
-git+https://github.com/FInAT/FInAT.git#egg=finat
+git+https://github.com/firedrakeproject/ufl.git@mse#egg=ufl
+git+https://github.com/firedrakeproject/fiat.git@mse#egg=fiat
+git+https://github.com/FInAT/FInAT.git@mse#egg=finat
 git+https://github.com/firedrakeproject/loopy.git@firedrake#egg=loopy

tests/test_create_fiat_element.py

@@ -29,7 +29,7 @@ def test_triangle_basic(ufl_element):
     assert isinstance(element, supported_elements[ufl_element.family()])
 
 
-@pytest.fixture(params=["CG", "DG"], scope="module")
+@pytest.fixture(params=["CG", "DG", "DG L2"], scope="module")
 def tensor_name(request):
     return request.param
 
@@ -62,7 +62,8 @@ def test_tensor_prod_simple(ufl_A, ufl_B):
 
 @pytest.mark.parametrize(('family', 'expected_cls'),
                          [('P', FIAT.Lagrange),
-                          ('DP', FIAT_DiscontinuousLagrange)])
+                          ('DP', FIAT_DiscontinuousLagrange),
+                          ('DP L2', FIAT_DiscontinuousLagrange)])
 def test_interval_variant_default(family, expected_cls):
     ufl_element = ufl.FiniteElement(family, ufl.interval, 3)
     assert isinstance(create_element(ufl_element), expected_cls)
@@ -71,8 +72,16 @@ def test_interval_variant_default(family, expected_cls):
 @pytest.mark.parametrize(('family', 'variant', 'expected_cls'),
                          [('P', 'equispaced', FIAT.Lagrange),
                           ('P', 'spectral', FIAT.GaussLobattoLegendre),
+                          ('P', 'mse', FIAT.GaussLobattoLegendre),
+                          ('P', 'dualmse', FIAT.ExtendedGaussLegendre),
                           ('DP', 'equispaced', FIAT_DiscontinuousLagrange),
-                          ('DP', 'spectral', FIAT.GaussLegendre)])
+                          ('DP', 'spectral', FIAT.GaussLegendre),
+                          ('DP', 'mse', FIAT.EdgeGaussLobattoLegendre),
+                          ('DP', 'dualmse', FIAT.EdgeExtendedGaussLegendre),
+                          ('DP L2', 'equispaced', FIAT_DiscontinuousLagrange),
+                          ('DP L2', 'spectral', FIAT.GaussLegendre),
+                          ('DP L2', 'mse', FIAT.EdgeGaussLobattoLegendre),
+                          ('DP L2', 'dualmse', FIAT.EdgeExtendedGaussLegendre)])
 def test_interval_variant(family, variant, expected_cls):
     ufl_element = ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
     assert isinstance(create_element(ufl_element), expected_cls)
@@ -84,18 +93,42 @@ def test_triangle_variant_spectral_fail():
         create_element(ufl_element)
 
 
+def test_triangle_variant_spectral_fail_l2():
+    ufl_element = ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
+    with pytest.raises(ValueError):
+        create_element(ufl_element)
+
+
 def test_quadrilateral_variant_spectral_q():
     element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
     assert isinstance(element.element.A, FIAT.GaussLobattoLegendre)
     assert isinstance(element.element.B, FIAT.GaussLobattoLegendre)
 
 
+def test_quadrilateral_variant_mse_q():
+    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='mse'))
+    assert isinstance(element.element.A, FIAT.GaussLobattoLegendre)
+    assert isinstance(element.element.B, FIAT.GaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_dualmse_q():
+    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='dualmse'))
+    assert isinstance(element.element.A, FIAT.ExtendedGaussLegendre)
+    assert isinstance(element.element.B, FIAT.ExtendedGaussLegendre)
+
+
 def test_quadrilateral_variant_spectral_dq():
     element = create_element(ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='spectral'))
     assert isinstance(element.element.A, FIAT.GaussLegendre)
     assert isinstance(element.element.B, FIAT.GaussLegendre)
 
 
+def test_quadrilateral_variant_spectral_dq_l2():
+    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
+    assert isinstance(element.element.A, FIAT.GaussLegendre)
+    assert isinstance(element.element.B, FIAT.GaussLegendre)
+
+
 def test_quadrilateral_variant_spectral_rtcf():
     element = create_element(ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='spectral'))
     assert isinstance(element.element._elements[0].A, FIAT.GaussLobattoLegendre)
@@ -104,6 +137,14 @@ def test_quadrilateral_variant_spectral_rtcf():
     assert isinstance(element.element._elements[1].B, FIAT.GaussLobattoLegendre)
 
 
+def test_quadrilateral_variant_spectral_rtce():
+    element = create_element(ufl.FiniteElement('RTCE', ufl.quadrilateral, 2, variant='spectral'))
+    assert isinstance(element.element._elements[0].A, FIAT.GaussLobattoLegendre)
+    assert isinstance(element.element._elements[0].B, FIAT.GaussLegendre)
+    assert isinstance(element.element._elements[1].A, FIAT.GaussLegendre)
+    assert isinstance(element.element._elements[1].B, FIAT.GaussLobattoLegendre)
+
+
 def test_cache_hit(ufl_element):
     A = create_element(ufl_element)
     B = create_element(ufl_element)

tests/test_create_finat_element.py

@@ -39,7 +39,7 @@ def test_triangle_vector(ufl_element, ufl_vector_element):
     assert scalar == vector.base_element
 
 
-@pytest.fixture(params=["CG", "DG"])
+@pytest.fixture(params=["CG", "DG", "DG L2"])
 def tensor_name(request):
     return request.param
 
@@ -70,7 +70,8 @@ def test_tensor_prod_simple(ufl_A, ufl_B):
 
 @pytest.mark.parametrize(('family', 'expected_cls'),
                          [('P', finat.Lagrange),
-                          ('DP', finat.DiscontinuousLagrange)])
+                          ('DP', finat.DiscontinuousLagrange),
+                          ('DP L2', finat.DiscontinuousLagrange)])
 def test_interval_variant_default(family, expected_cls):
     ufl_element = ufl.FiniteElement(family, ufl.interval, 3)
     assert isinstance(create_element(ufl_element), expected_cls)
@@ -79,8 +80,16 @@ def test_interval_variant_default(family, expected_cls):
 @pytest.mark.parametrize(('family', 'variant', 'expected_cls'),
                          [('P', 'equispaced', finat.Lagrange),
                           ('P', 'spectral', finat.GaussLobattoLegendre),
+                          ('P', 'mse', finat.GaussLobattoLegendre),
+                          ('P', 'dualmse', finat.ExtendedGaussLegendre),
                           ('DP', 'equispaced', finat.DiscontinuousLagrange),
-                          ('DP', 'spectral', finat.GaussLegendre)])
+                          ('DP', 'spectral', finat.GaussLegendre),
+                          ('DP', 'mse', finat.EdgeGaussLobattoLegendre),
+                          ('DP', 'dualmse', finat.EdgeExtendedGaussLegendre),
+                          ('DP L2', 'equispaced', finat.DiscontinuousLagrange),
+                          ('DP L2', 'spectral', finat.GaussLegendre),
+                          ('DP L2', 'mse', finat.EdgeGaussLobattoLegendre),
+                          ('DP L2', 'dualmse', finat.EdgeExtendedGaussLegendre)])
 def test_interval_variant(family, variant, expected_cls):
     ufl_element = ufl.FiniteElement(family, ufl.interval, 3, variant=variant)
     assert isinstance(create_element(ufl_element), expected_cls)
@@ -92,6 +101,12 @@ def test_triangle_variant_spectral_fail():
         create_element(ufl_element)
 
 
+def test_triangle_variant_spectral_fail_l2():
+    ufl_element = ufl.FiniteElement('DP L2', ufl.triangle, 2, variant='spectral')
+    with pytest.raises(ValueError):
+        create_element(ufl_element)
+
+
 def test_quadrilateral_variant_spectral_q():
     element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='spectral'))
     assert isinstance(element.product.factors[0], finat.GaussLobattoLegendre)
@@ -104,6 +119,80 @@ def test_quadrilateral_variant_spectral_dq():
     assert isinstance(element.product.factors[1], finat.GaussLegendre)
 
 
+def test_quadrilateral_variant_spectral_dq_l2():
+    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='spectral'))
+    assert isinstance(element.product.factors[0], finat.GaussLegendre)
+    assert isinstance(element.product.factors[1], finat.GaussLegendre)
+
+
+def test_quadrilateral_variant_mse_q():
+    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='mse'))
+    assert isinstance(element.product.factors[0], finat.GaussLobattoLegendre)
+    assert isinstance(element.product.factors[1], finat.GaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_dualmse_q():
+    element = create_element(ufl.FiniteElement('Q', ufl.quadrilateral, 3, variant='dualmse'))
+    assert isinstance(element.product.factors[0], finat.ExtendedGaussLegendre)
+    assert isinstance(element.product.factors[1], finat.ExtendedGaussLegendre)
+
+
+def test_quadrilateral_variant_mse_dq():
+    element = create_element(ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='mse'))
+    assert isinstance(element.product.factors[0], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.factors[1], finat.EdgeGaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_mse_dq_l2():
+    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='mse'))
+    assert isinstance(element.product.factors[0], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.factors[1], finat.EdgeGaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_dualmse_dq():
+    element = create_element(ufl.FiniteElement('DQ', ufl.quadrilateral, 1, variant='dualmse'))
+    assert isinstance(element.product.factors[0], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.factors[1], finat.EdgeExtendedGaussLegendre)
+
+
+def test_quadrilateral_variant_dualmse_dq_l2():
+    element = create_element(ufl.FiniteElement('DQ L2', ufl.quadrilateral, 1, variant='dualmse'))
+    assert isinstance(element.product.factors[0], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.factors[1], finat.EdgeExtendedGaussLegendre)
+
+
+def test_quadrilateral_variant_mse_rtcf():
+    element = create_element(ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='mse'))
+    assert isinstance(element.product.elements[0].wrappee.factors[0], finat.GaussLobattoLegendre)
+    assert isinstance(element.product.elements[0].wrappee.factors[1], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[0], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[1], finat.GaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_mse_rtce():
+    element = create_element(ufl.FiniteElement('RTCE', ufl.quadrilateral, 2, variant='mse'))
+    assert isinstance(element.product.elements[0].wrappee.factors[0], finat.GaussLobattoLegendre)
+    assert isinstance(element.product.elements[0].wrappee.factors[1], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[0], finat.EdgeGaussLobattoLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[1], finat.GaussLobattoLegendre)
+
+
+def test_quadrilateral_variant_dualmse_rtcf():
+    element = create_element(ufl.FiniteElement('RTCF', ufl.quadrilateral, 2, variant='dualmse'))
+    assert isinstance(element.product.elements[0].wrappee.factors[0], finat.ExtendedGaussLegendre)
+    assert isinstance(element.product.elements[0].wrappee.factors[1], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[0], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[1], finat.ExtendedGaussLegendre)
+
+
+def test_quadrilateral_variant_dualmse_rtce():
+    element = create_element(ufl.FiniteElement('RTCE', ufl.quadrilateral, 2, variant='dualmse'))
+    assert isinstance(element.product.elements[0].wrappee.factors[0], finat.ExtendedGaussLegendre)
+    assert isinstance(element.product.elements[0].wrappee.factors[1], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[0], finat.EdgeExtendedGaussLegendre)
+    assert isinstance(element.product.elements[1].wrappee.factors[1], finat.ExtendedGaussLegendre)
+
+
 def test_cache_hit(ufl_element):
     A = create_element(ufl_element)
     B = create_element(ufl_element)

tests/test_underintegration.py

@@ -10,9 +10,9 @@
                  action, interval, quadrilateral, dot, grad)
 
 from FIAT import ufc_cell
-from FIAT.quadrature import GaussLobattoLegendreQuadratureLineRule, GaussLegendreQuadratureLineRule
+from FIAT.quadrature import GaussLobattoLegendreQuadratureLineRule, GaussLegendreQuadratureLineRule, ExtendedGaussLegendreQuadratureLineRule
 
-from finat.point_set import GaussLobattoLegendrePointSet, GaussLegendrePointSet
+from finat.point_set import GaussLobattoLegendrePointSet, GaussLegendrePointSet, ExtendedGaussLegendrePointSet
 from finat.quadrature import QuadratureRule, TensorProductQuadratureRule
 
 from tsfc import compile_form
@@ -38,6 +38,16 @@ def gl_quadrature_rule(cell, elem_deg):
     return finat_rule
 
 
+def egl_quadrature_rule(cell, elem_deg):
+    fiat_cell = ufc_cell("interval")
+    fiat_rule = ExtendedGaussLegendreQuadratureLineRule(fiat_cell, elem_deg + 1)
+    line_rules = [QuadratureRule(ExtendedGaussLegendrePointSet(fiat_rule.get_points()),
+                                 fiat_rule.get_weights())
+                  for _ in range(cell.topological_dimension())]
+    finat_rule = reduce(lambda a, b: TensorProductQuadratureRule([a, b]), line_rules)
+    return finat_rule
+
+
 def mass_cg(cell, degree):
     m = Mesh(VectorElement('Q', cell, 1))
     V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='spectral'))
@@ -46,6 +56,14 @@ def mass_cg(cell, degree):
     return u*v*dx(scheme=gll_quadrature_rule(cell, degree))
 
 
+def mass_cg_egl(cell, degree):
+    m = Mesh(VectorElement('Q', cell, 1))
+    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='dualmse'))
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    return u*v*dx(scheme=egl_quadrature_rule(cell, degree))
+
+
 def mass_dg(cell, degree):
     m = Mesh(VectorElement('Q', cell, 1))
     V = FunctionSpace(m, FiniteElement('DQ', cell, degree, variant='spectral'))
@@ -62,12 +80,20 @@ def laplace(cell, degree):
     return dot(grad(u), grad(v))*dx(scheme=gll_quadrature_rule(cell, degree))
 
 
+def laplace_dualmse(cell, degree):
+    m = Mesh(VectorElement('Q', cell, 1))
+    V = FunctionSpace(m, FiniteElement('Q', cell, degree, variant='dualmse'))
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    return dot(grad(u), grad(v))*dx(scheme=egl_quadrature_rule(cell, degree))
+
+
 def count_flops(form):
     kernel, = compile_form(form, parameters=dict(mode='spectral'))
     return EstimateFlops().visit(kernel.ast)
 
 
-@pytest.mark.parametrize('form', [mass_cg, mass_dg])
+@pytest.mark.parametrize('form', [mass_cg, mass_dg, mass_cg_egl])
 @pytest.mark.parametrize(('cell', 'order'),
                          [(quadrilateral, 2),
                           (TensorProductCell(interval, interval), 2),
@@ -79,7 +105,7 @@ def test_mass(form, cell, order):
     assert (rates < order).all()
 
 
-@pytest.mark.parametrize('form', [mass_cg, mass_dg])
+@pytest.mark.parametrize('form', [mass_cg, mass_dg, mass_cg_egl])
 @pytest.mark.parametrize(('cell', 'order'),
                          [(quadrilateral, 2),
                           (TensorProductCell(interval, interval), 2),
@@ -102,6 +128,17 @@ def test_laplace(cell, order):
     assert (rates < order).all()
 
 
+@pytest.mark.parametrize(('cell', 'order'),
+                         [(quadrilateral, 4),
+                          (TensorProductCell(interval, interval), 4),
+                          (TensorProductCell(quadrilateral, interval), 5)])
+def test_laplace_dualmse(cell, order):
+    degrees = numpy.arange(4, 10)
+    flops = [count_flops(laplace_dualmse(cell, int(degree))) for degree in degrees]
+    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees + 1))
+    assert (rates < order).all()
+
+
 if __name__ == "__main__":
     import os
     import sys

tsfc/fiatinterface.py

@@ -61,7 +61,16 @@
     "NCE": None,
     "NCF": None,
     "DPC": FIAT.DPC,
-    "S": FIAT.Serendipity
+    "S": FIAT.Serendipity,
+    "DPC L2": FIAT.DPC,
+    "Discontinuous Lagrange L2": FIAT.DiscontinuousLagrange,
+    "Gauss-Legendre L2": FIAT.GaussLegendre,
+    "DQ L2": None,
+    "Extended-Gauss-Legendre": FIAT.ExtendedGaussLegendre,
+    "Gauss-Lobatto-Legendre Edge": FIAT.EdgeGaussLobattoLegendre,
+    "Gauss-Lobatto-Legendre Edge L2": FIAT.EdgeGaussLobattoLegendre,
+    "Extended-Gauss-Legendre Edge": FIAT.EdgeExtendedGaussLegendre,
+    "Extended-Gauss-Legendre Edge L2": FIAT.EdgeExtendedGaussLegendre
 }
 """A :class:`.dict` mapping UFL element family names to their
 FIAT-equivalent constructors.  If the value is ``None``, the UFL
@@ -131,16 +140,24 @@ def convert_finiteelement(element, vector_is_mixed):
             lmbda = FIAT.Lagrange
         elif kind == 'spectral' and element.cell().cellname() == 'interval':
             lmbda = FIAT.GaussLobattoLegendre
+        elif kind == 'mse' and element.cell().cellname() == 'interval':
+            lmbda = FIAT.GaussLobattoLegendre
+        elif kind == 'dualmse' and element.cell().cellname() == 'interval':
+            lmbda = FIAT.ExtendedGaussLegendre
         else:
             raise ValueError("Variant %r not supported on %s" % (kind, element.cell()))
-    elif element.family() == "Discontinuous Lagrange":
+    elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]:
         if kind == 'equispaced':
             lmbda = FIAT.DiscontinuousLagrange
         elif kind == 'spectral' and element.cell().cellname() == 'interval':
             lmbda = FIAT.GaussLegendre
+        elif kind == 'mse' and element.cell().cellname() == 'interval':
+            lmbda = FIAT.EdgeGaussLobattoLegendre
+        elif kind == 'dualmse' and element.cell().cellname() == 'interval':
+            lmbda = FIAT.EdgeExtendedGaussLegendre
         else:
             raise ValueError("Variant %r not supported on %s" % (kind, element.cell()))
-    elif element.family() == "DPC":
+    elif element.family() in ["DPC", "DPC L2"]:
         if element.cell().geometric_dimension() == 2:
             element = element.reconstruct(cell=ufl.hypercube(2))
         elif element.cell.geometric_dimension() == 3: