Add GLS H(curl div) element

FIAT/__init__.py

@@ -36,6 +36,7 @@
 from FIAT.raviart_thomas import RaviartThomas
 from FIAT.crouzeix_raviart import CrouzeixRaviart
 from FIAT.regge import Regge
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberl
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
 from FIAT.arnold_winther import ArnoldWinther
 from FIAT.arnold_winther import ArnoldWintherNC

FIAT/gopalakrishnan_lederer_schoberl.py

@@ -0,0 +1,131 @@
+from FIAT import finite_element, polynomial_set, dual_set, functional
+import numpy
+
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.expansions import polynomial_entity_ids
+
+
+def mask_facet_ids(ref_el, dim, constrained_ids, mask):
+    closure_ids = dual_set.make_entity_closure_ids(ref_el, constrained_ids)
+    mask.fill(1)
+    for facet in closure_ids[dim]:
+        mask[facet][..., closure_ids[dim][facet]] = 0
+    indices = numpy.flatnonzero(mask)
+    return indices
+
+
+def make_polynomial_sets(ref_el, degree):
+    sd = ref_el.get_spatial_dimension()
+    phi = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree, variant="bubble")
+    expansion_set = phi.get_expansion_set()
+    entity_ids = polynomial_entity_ids(ref_el, degree, continuity=expansion_set.continuity)
+    mask = numpy.ones((phi.get_num_members(),), int).reshape(sd+1, sd-1, -1)
+
+    # extract bubbles
+    bubble_ids = mask_facet_ids(ref_el, sd-1, entity_ids, mask)
+    P_bubble = phi.take(bubble_ids)
+
+    # build constrained space with normal-tangential component in P_{k-1}
+    constrained_ids = {}
+    for dim in entity_ids:
+        constrained_ids[dim] = {}
+        if dim == 0 or dim == sd:
+            for entity in entity_ids[dim]:
+                constrained_ids[dim][entity] = []
+        else:
+            dimPkm1 = len(ref_el.make_points(dim, 0, degree-1))
+            for entity in entity_ids[dim]:
+                constrained_ids[dim][entity] = entity_ids[dim][entity][dimPkm1:]
+
+    indices = mask_facet_ids(ref_el, sd-1, constrained_ids, mask)
+    Sigma = phi.take(indices)
+    return P_bubble, Sigma
+
+
+class GLSDualSet(dual_set.DualSet):
+
+    def __init__(self, ref_el, degree, bubbles):
+        FIM = functional.FrobeniusIntegralMoment
+        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: bidirectional nt Legendre moments
+        dim = sd - 1
+        ref_facet = ref_el.construct_subelement(dim)
+        Qref = create_quadrature(ref_facet, 2*degree-1)
+        P = polynomial_set.ONPolynomialSet(ref_facet, degree-1)
+        phis = P.tabulate(Qref.get_points())[(0,) * dim]
+
+        for facet in sorted(top[dim]):
+            cur = len(nodes)
+            Q = FacetQuadratureRule(ref_el, dim, facet, Qref)
+            Jdet = Q.jacobian_determinant()
+            tangents = ref_el.compute_normalized_tangents(dim, facet)
+            normal = ref_el.compute_normal(facet)
+            normal /= numpy.linalg.norm(normal)
+            scaled_normal = normal / Jdet
+            comps = [numpy.outer(scaled_normal, that) for that in tangents]
+            nodes.extend(FIM(ref_el, Q, comp[:, :, None] * phi[None, None, :])
+                         for phi in phis for comp in comps)
+            entity_ids[dim][facet].extend(range(cur, len(nodes)))
+
+        # Interior dofs: moments against nt bubbles
+        Q = create_quadrature(ref_el, 2*degree)
+        phis = bubbles.tabulate(Q.get_points())[(0,) * sd]
+        cur = len(nodes)
+        nodes.extend(FIM(ref_el, Q, phi) for phi in phis)
+        entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super(GLSDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+
+class GopalakrishnanLedererSchoberl(finite_element.CiarletElement):
+
+    def __init__(self, ref_el, degree):
+        bubbles, poly_set = make_polynomial_sets(ref_el, degree)
+        dual = GLSDualSet(ref_el, degree, bubbles)
+        mapping = "covariant contravariant piola"
+        super(GopalakrishnanLedererSchoberl, self).__init__(poly_set, dual, degree, mapping=mapping)
+
+
+if __name__ == "__main__":
+    from FIAT import ufc_simplex
+    from FIAT.expansions import polynomial_dimension
+
+    degree = 4
+    ref_el = ufc_simplex(3)
+    GLS(ref_el, degree)
+
+    bubbles, poly_set = make_polynomial_sets(ref_el, degree)
+    expansion_set = poly_set.get_expansion_set()
+
+    # test dimension of constrained space
+    sd = ref_el.get_spatial_dimension()
+    facet_el = ref_el.construct_subelement(sd-1)
+    dimPkm1 = polynomial_dimension(facet_el, degree-1)
+    dimPk = polynomial_dimension(facet_el, degree)
+    expected = (sd**2-1)*(expansion_set.get_num_members(degree) - (dimPk - dimPkm1))
+    assert poly_set.get_num_members() == expected
+
+    # test normal-tangential bubble support
+    facet_el = ref_el.construct_subelement(sd-1)
+    Qref = create_quadrature(facet_el, 2*degree)
+    norms = numpy.zeros((bubbles.get_num_members(),))
+    top = ref_el.get_topology()
+    for facet in top[sd-1]:
+        Q = FacetQuadratureRule(ref_el, sd-1, facet, Qref)
+        qpts, qwts = Q.get_points(), Q.get_weights()
+        phi_at_pts = bubbles.tabulate(qpts)[(0,) * sd]
+        n = ref_el.compute_normal(facet)
+        rts = ref_el.compute_normalized_tangents(sd-1, facet)
+        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)
+    expected = (sd**2-1)*expansion_set.get_num_members(degree-1)
+    assert bubbles.get_num_members() == expected

FIAT/polynomial_set.py

@@ -246,40 +246,39 @@ def __init__(self, ref_el, degree, size=None, **kwargs):
         embedded_degree = degree
 
         # set up coefficients for traceless tensors
-        normalize = lambda u: u / numpy.linalg.norm(u)
         top = ref_el.get_topology()
         verts = ref_el.get_vertices()
         v0 = numpy.array(verts[0])
         rts = [numpy.array(v1) - v0 for v1 in verts[1:]]
         rts.insert(0, -sum(rts))
 
-        #rts = [ref_el.compute_tangents(sd-1, e) for e in top[sd-1]]
-        #rts = list(map(normalize, rts))
+        normalize = lambda u: u / numpy.linalg.norm(u)
+        rts = list(map(normalize, rts))
 
         dev = lambda S: S - (numpy.trace(S) / S.shape[0]) * numpy.eye(*S.shape)
-        basis = numpy.zeros((num_components, *shape), "d")
+        basis = numpy.zeros((len(top[sd-1]), sd-1, *shape), "d")
         if size == 2:
             R = numpy.array([[0, 1], [-1, 0]])
-            for i in range(num_components):
+            for i in top[sd-1]:
                 i1 = (i + 1) % (sd+1)
                 i2 = (i + 2) % (sd+1)
-                basis[i] = dev(numpy.outer(rts[i1], R @ rts[i2]))
+                basis[i, 0] = dev(numpy.outer(rts[i1], R @ rts[i2]))
         elif size == 3:
-            for i in range(num_components):
+            for i in top[sd-1]:
                 i1 = (i + 1) % (sd+1)
                 i2 = (i + 2) % (sd+1)
                 i3 = (i + 3) % (sd+1)
-                if i > sd:
-                    i1, i2, i3 = i2, i3, i1
-                basis[i] = dev(numpy.outer(rts[i1], numpy.cross(rts[i2], rts[i3])))
+                for j in range(sd-1):
+                    if j > 0:
+                        i1, i2, i3 = i2, i3, i1
+                    basis[i, j] = dev(numpy.outer(rts[i1], numpy.cross(rts[i2], rts[i3])))
         else:
             raise NotImplementedError("TODO")
 
-        print(basis)
         coeffs_shape = (num_members, *shape, num_exp_functions)
         coeffs = numpy.zeros(coeffs_shape, "d")
         cur_bf = 0
-        for S in basis:
+        for S in basis.reshape(-1, *shape):
             for exp_bf in range(num_exp_functions):
                 coeffs[cur_bf, :, :, exp_bf] = S
                 cur_bf += 1
@@ -311,51 +310,3 @@ def make_bubbles(ref_el, degree, codim=0, shape=()):
         indices = list((numpy.array(indices)[:, None] + dimPk * numpy.arange(ncomp)[None, :]).flat)
     poly_set = poly_set.take(indices)
     return poly_set
-
-
-if __name__ == "__main__":
-
-    from FIAT import ufc_simplex
-    from FIAT.quadrature_schemes import create_quadrature
-    from FIAT.quadrature import FacetQuadratureRule
-    from FIAT.dual_set import make_entity_closure_ids
-    from FIAT.expansions import polynomial_entity_ids
-
-    ref_el = ufc_simplex(3)
-    sd = ref_el.get_spatial_dimension()
-    degree = 4
-    phi = TracelessTensorPolynomialSet(ref_el, degree, variant="bubble")
-    expansion_set = phi.get_expansion_set()
-
-    ncomp = sd**2 - 1
-    entity_ids = polynomial_entity_ids(ref_el, degree, continuity=expansion_set.continuity)
-    closure_ids = make_entity_closure_ids(ref_el, entity_ids)
-    mask = numpy.ones((phi.get_num_members(),), int).reshape((ncomp, -1))
-    for component in range(ncomp):
-        facet = component % len(closure_ids[sd-1])
-        mask[component][closure_ids[sd-1][facet]] = 0
-    bubble_ids = numpy.flatnonzero(mask)
-    print(bubble_ids)
-
-    facet_el = ref_el.construct_subelement(sd-1)
-    Qref = create_quadrature(facet_el, 2*degree)
-    top = ref_el.get_topology()
-    norms = numpy.zeros((phi.get_num_members(),))
-    for facet in top[sd-1]:
-        Q = FacetQuadratureRule(ref_el, sd-1, facet, Qref)
-        qpts, qwts = Q.get_points(), Q.get_weights()
-        phi_at_pts = phi.tabulate(qpts)[(0,) * sd]
-        n = ref_el.compute_normal(facet)
-        rts = ref_el.compute_normalized_tangents(sd-1, facet)
-        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)
-
-    bubbles = numpy.flatnonzero(norms < 1E-12)
-    expected = (3*degree*(degree+1))//2 if sd == 2 else (8*degree*(degree+1)*(degree+2))//6
-    assert len(bubbles) == expected
-    assert numpy.allclose(bubbles, bubble_ids)
-
-    print(len(bubbles), phi.get_num_members())
-    print(bubbles)