Implement a traceless H(curl div) element (#69)

* TracelessTensorPolynomialSet

* Add GLS H(curl div) element

FIAT/__init__.py

@@ -48,6 +48,8 @@
 from FIAT.raviart_thomas import RaviartThomas
 from FIAT.crouzeix_raviart import CrouzeixRaviart
 from FIAT.regge import Regge
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
 from FIAT.arnold_winther import ArnoldWinther
 from FIAT.arnold_winther import ArnoldWintherNC
@@ -108,6 +110,8 @@
                       "BrokenElement": DiscontinuousElement,
                       "HDiv Trace": HDivTrace,
                       "Hellan-Herrmann-Johnson": HellanHerrmannJohnson,
+                      "Gopalakrishnan-Lederer-Schoberl 1st kind": GopalakrishnanLedererSchoberlFirstKind,
+                      "Gopalakrishnan-Lederer-Schoberl 2nd kind": GopalakrishnanLedererSchoberlSecondKind,
                       "Conforming Arnold-Winther": ArnoldWinther,
                       "Nonconforming Arnold-Winther": ArnoldWintherNC,
                       "Hu-Zhang": HuZhang,

FIAT/functional.py

@@ -631,13 +631,13 @@ class TensorBidirectionalIntegralMoment(FrobeniusIntegralMoment):
     field, f a function tabulated at points, and v,w be vectors. This implements the evaluation
     \int v^T u(x) w f(x).
     Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed
-    from the Frobenius inner product of u with wv^T. This gives the
+    from the Frobenius inner product of u with vw^T. This gives the
     correct weights.
     """
 
     def __init__(self, ref_el, v, w, Q, f_at_qpts):
-        wvT = numpy.outer(w, v)
-        F_at_qpts = numpy.multiply(wvT[..., None], f_at_qpts)
+        vwT = numpy.outer(v, w)
+        F_at_qpts = numpy.multiply(vwT[..., None], f_at_qpts)
         super().__init__(ref_el, Q, F_at_qpts, "TensorBidirectionalMomentInnerProductEvaluation")
 
 

FIAT/gopalakrishnan_lederer_schoberl.py

@@ -0,0 +1,88 @@
+from FIAT import finite_element, dual_set, polynomial_set, expansions
+from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.restricted import RestrictedElement
+
+
+class GLSDual(dual_set.DualSet):
+    def __init__(self, ref_el, degree):
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
+        nodes = []
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+        # Face dofs: moments of normal-tangential components against a basis for Pk
+        ref_facet = ref_el.construct_subelement(sd-1)
+        Qref = create_quadrature(ref_facet, 2*degree)
+        P = polynomial_set.ONPolynomialSet(ref_facet, degree)
+        phis = P.tabulate(Qref.get_points())[(0,) * (sd-1)]
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            Q = FacetQuadratureRule(ref_el, sd-1, f, Qref)
+            Jdet = Q.jacobian_determinant()
+            normal = ref_el.compute_scaled_normal(f)
+            tangents = ref_el.compute_tangents(sd-1, f)
+            n = normal / Jdet
+            nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                         for phi in phis for t in tangents)
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        # Interior dofs: moments of normal-tangential components against a basis for P_{k-1}
+        if degree > 0:
+            cur = len(nodes)
+            Q = create_quadrature(ref_el, 2*degree-1)
+            P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
+            phis = P.tabulate(Q.get_points())[(0,) * sd]
+            for f in sorted(top[sd-1]):
+                n = ref_el.compute_scaled_normal(f)
+                tangents = ref_el.compute_tangents(sd-1, f)
+                nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                             for phi in phis for t in tangents)
+            entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement):
+    """The GLS element used for the Mass-Conserving mixed Stress (MCS)
+    formulation for Stokes flow with weakly imposed stress symmetry.
+
+    GLS^2(k) is the space of trace-free polynomials of degree k with
+    continuous normal-tangential components.
+
+    Reference: https://doi.org/10.1137/19M1248960
+
+    Notes
+    -----
+    This element does not include the bubbles required for inf-sup stability of
+    the weak symmetry constraint.
+
+    """
+    def __init__(self, ref_el, degree):
+        poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree)
+        dual = GLSDual(ref_el, degree)
+        sd = ref_el.get_spatial_dimension()
+        formdegree = (1, sd-1)
+        mapping = "covariant contravariant piola"
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
+
+
+def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree):
+    """The GLS element used for the Mass-Conserving mixed Stress (MCS)
+    formulation for Stokes flow.
+
+    GLS^1(k) is the space of trace-free polynomials of degree k with
+    continuous normal-tangential components of degree k-1.
+
+    Reference: https://doi.org/10.1093/imanum/drz022
+    """
+    fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree)
+    entity_dofs = fe.entity_dofs()
+    sd = ref_el.get_spatial_dimension()
+    dimPkm1 = (sd-1)*expansions.polynomial_dimension(ref_el.construct_subelement(sd-1), degree-1)
+    indices = []
+    for f in entity_dofs[sd-1]:
+        indices.extend(entity_dofs[sd-1][f][:dimPkm1])
+    indices.extend(entity_dofs[sd][0])
+    return RestrictedElement(fe, indices=indices)

FIAT/guzman_neilan.py

@@ -109,15 +109,12 @@ def GuzmanNeilanH1div(ref_el, degree=2, reduced=False):
     """
     order = 0
     AS = AlfeldSorokina(ref_el, 2)
-    if reduced:
+    if reduced or ref_el.get_spatial_dimension() <= 2:
         order = 1
         # Only extract the div bubbles
         div_nodes = [i for i, node in enumerate(AS.dual_basis())
                      if len(node.deriv_dict) > 0]
         AS = RestrictedElement(AS, indices=div_nodes)
-    elif ref_el.get_spatial_dimension() <= 2:
-        # Quadratic bubbles are already included in 2D
-        return AS
     GN = GuzmanNeilanH1(ref_el, order=order)
     return NodalEnrichedElement(AS, GN)
 

FIAT/polynomial_set.py

@@ -125,16 +125,13 @@ def __init__(self, ref_el, degree, shape=tuple(), **kwargs):
         if shape == tuple():
             coeffs = numpy.eye(num_members)
         else:
-            coeffs_shape = (num_members, *shape, num_exp_functions)
-            coeffs = numpy.zeros(coeffs_shape, "d")
-            # use functional's index_iterator function
-            cur_bf = 0
+            coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+            cur = 0
+            exp_bf = range(num_exp_functions)
             for idx in index_iterator(shape):
-                n = expansion_set.get_num_members(embedded_degree)
-                for exp_bf in range(n):
-                    cur_idx = (cur_bf, *idx, exp_bf)
-                    coeffs[cur_idx] = 1.0
-                    cur_bf += 1
+                cur_bf = range(cur, cur+num_exp_functions)
+                coeffs[(cur_bf, *idx, exp_bf)] = 1.0
+                cur += num_exp_functions
 
         super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
 
@@ -209,19 +206,49 @@ def __init__(self, ref_el, degree, size=None, **kwargs):
         embedded_degree = degree
 
         # set up coefficients for symmetric tensors
-        coeffs_shape = (num_members, *shape, num_exp_functions)
-        coeffs = numpy.zeros(coeffs_shape, "d")
-        cur_bf = 0
+        coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+        cur = 0
+        exp_bf = range(num_exp_functions)
         for i, j in index_iterator(shape):
+            if i > j:
+                continue
+            cur_bf = range(cur, cur+num_exp_functions)
+            coeffs[cur_bf, i, j, exp_bf] = 1.0
+            coeffs[cur_bf, j, i, exp_bf] = 1.0
+            cur += num_exp_functions
+
+        super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
+
+
+class TracelessTensorPolynomialSet(PolynomialSet):
+    """Constructs an orthonormal basis for traceless-tensor-valued
+    polynomials on a reference element.
+    """
+    def __init__(self, ref_el, degree, size=None, **kwargs):
+        expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
+
+        sd = ref_el.get_spatial_dimension()
+        if size is None:
+            size = sd
+
+        shape = (size, size)
+        num_exp_functions = expansion_set.get_num_members(degree)
+        num_components = size * size - 1
+        num_members = num_components * num_exp_functions
+        embedded_degree = degree
+
+        # set up coefficients for traceless tensors
+        coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+        cur = 0
+        exp_bf = range(num_exp_functions)
+        for i, j in index_iterator(shape):
+            if i == size-1 and j == size-1:
+                continue
+            cur_bf = range(cur, cur+num_exp_functions)
+            coeffs[cur_bf, i, j, exp_bf] = 1.0
             if i == j:
-                for exp_bf in range(num_exp_functions):
-                    coeffs[cur_bf, i, j, exp_bf] = 1.0
-                    cur_bf += 1
-            elif i < j:
-                for exp_bf in range(num_exp_functions):
-                    coeffs[cur_bf, i, j, exp_bf] = 1.0
-                    coeffs[cur_bf, j, i, exp_bf] = 1.0
-                    cur_bf += 1
+                coeffs[cur_bf, -1, -1, exp_bf] = -1.0
+            cur += num_exp_functions
 
         super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
 

test/unit/test_fiat.py

@@ -35,6 +35,8 @@
 from FIAT.regge import Regge                                    # noqa: F401
 from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet   # noqa: F401
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson  # noqa: F401
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind  # noqa: F401
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind  # noqa: F401
 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini  # noqa: F401
 from FIAT.gauss_legendre import GaussLegendre                   # noqa: F401
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre    # noqa: F401
@@ -273,6 +275,18 @@ def __init__(self, a, b):
     "HellanHerrmannJohnson(S, 2)",
     "HellanHerrmannJohnson(T, 1, variant='point')",
     "HellanHerrmannJohnson(S, 1, variant='point')",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 1)",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 2)",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 3)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 1)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 2)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 3)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 0)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 1)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 2)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 0)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 1)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 2)",
     "BrezziDouglasFortinMarini(T, 2)",
     "GaussLegendre(I, 0)",
     "GaussLegendre(I, 1)",

test/unit/test_gopalakrishnan_lederer_schoberl.py

@@ -0,0 +1,77 @@
+import pytest
+import numpy
+
+from FIAT import (GopalakrishnanLedererSchoberlFirstKind,
+                  GopalakrishnanLedererSchoberlSecondKind)
+from FIAT.reference_element import ufc_simplex
+from FIAT.expansions import polynomial_dimension
+from FIAT.polynomial_set import ONPolynomialSet
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
+
+
+@pytest.fixture(params=("T", "S"))
+def cell(request):
+    dim = {"I": 1, "T": 2, "S": 3}[request.param]
+    return ufc_simplex(dim)
+
+
+@pytest.mark.parametrize("degree", (1, 2, 3))
+@pytest.mark.parametrize("kind", (1, 2))
+def test_gls_bubbles(kind, cell, degree):
+    if kind == 1:
+        element = GopalakrishnanLedererSchoberlFirstKind
+    else:
+        element = GopalakrishnanLedererSchoberlSecondKind
+    fe = element(cell, degree)
+    sd = cell.get_spatial_dimension()
+    facet_el = cell.construct_subelement(sd-1)
+    poly_set = fe.get_nodal_basis()
+
+    # test dimension of constrained space
+    dimPkm1 = polynomial_dimension(facet_el, degree-1)
+    dimPkp1 = polynomial_dimension(facet_el, degree+1)
+    dimPk = polynomial_dimension(facet_el, degree)
+    if kind == 1:
+        constraints = dimPk - dimPkm1
+    else:
+        constraints = 0
+    expected = (sd**2-1)*(polynomial_dimension(cell, degree) - constraints)
+    assert poly_set.get_num_members() == expected
+
+    # test dimension of the bubbles
+    entity_dofs = fe.entity_dofs()
+    bubbles = poly_set.take(entity_dofs[sd][0])
+    expected = (sd**2-1)*polynomial_dimension(cell, degree-1)
+    assert bubbles.get_num_members() == expected
+
+    top = cell.get_topology()
+    Qref = create_quadrature(facet_el, 2*degree+1)
+    Pk = ONPolynomialSet(facet_el, degree+1)
+    if kind == 1:
+        start, stop = dimPkm1, dimPkp1
+    else:
+        start, stop = dimPk, dimPkp1
+    PkH = Pk.take(list(range(start, stop)))
+    PkH_at_qpts = PkH.tabulate(Qref.get_points())[(0,)*(sd-1)]
+    weights = numpy.transpose(numpy.multiply(PkH_at_qpts, Qref.get_weights()))
+    for facet in top[sd-1]:
+        n = cell.compute_scaled_normal(facet)
+        rts = cell.compute_tangents(sd-1, facet)
+        Q = FacetQuadratureRule(cell, sd-1, facet, Qref)
+        qpts, qwts = Q.get_points(), Q.get_weights()
+
+        # test the degree of normal-tangential components
+        phi_at_pts = fe.tabulate(0, qpts)[(0,) * sd]
+        for t in rts:
+            nt = numpy.outer(t, n)
+            phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2)))
+            assert numpy.allclose(numpy.dot(phi_nt, weights), 0)
+
+        # test the support of the normal-tangential bubble
+        phi_at_pts = bubbles.tabulate(qpts)[(0,) * sd]
+        for t in rts:
+            nt = numpy.outer(t, n)
+            phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2)))
+            norms = numpy.dot(phi_nt**2, qwts)
+            assert numpy.allclose(norms, 0)