Add Johnson-Mercier elements (#67)

* HDivSymPolynomialSet add the Johnson-Mercier macro element

---------

Co-authored-by: Rob Kirby <robert.c.kirby@gmail.com>

Author: pbrubeck

FIAT/__init__.py

@@ -6,11 +6,11 @@
 
 # Import finite element classes
 from FIAT.finite_element import FiniteElement, CiarletElement  # noqa: F401
-from FIAT.argyris import Argyris
-from FIAT.hct import HsiehCloughTocher
+from FIAT.argyris import Argyris, QuinticArgyris
 from FIAT.bernstein import Bernstein
 from FIAT.bell import Bell
-from FIAT.argyris import QuinticArgyris
+from FIAT.hct import HsiehCloughTocher
+from FIAT.johnson_mercier import JohnsonMercier
 from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
 from FIAT.Sminus import TrimmedSerendipityEdge  # noqa: F401
 from FIAT.Sminus import TrimmedSerendipityFace  # noqa: F401
@@ -62,7 +62,6 @@
 
 # List of supported elements and mapping to element classes
 supported_elements = {"Argyris": Argyris,
-                      "HsiehCloughTocher": HsiehCloughTocher,
                       "Bell": Bell,
                       "Bernstein": Bernstein,
                       "Brezzi-Douglas-Marini": BrezziDouglasMarini,
@@ -78,6 +77,8 @@
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,
+                      "HsiehCloughTocher": HsiehCloughTocher,
+                      "Johnson-Mercier": JohnsonMercier,
                       "Lagrange": Lagrange,
                       "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen,
                       "Gauss-Lobatto-Legendre": GaussLobattoLegendre,

FIAT/expansions.py

@@ -84,9 +84,6 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
 
     if variant == "bubble":
         scale = -scale
-    if n == 0:
-        # Always return 1 for n=0 to make regression tests pass
-        scale = 1.0
 
     num_members = math.comb(n + dim, dim)
     results = tuple([None] * num_members for i in range(order+1))
@@ -286,18 +283,24 @@ def __init__(self, ref_el, scale=None, variant=None):
                                 base_verts) for cell in top[sd]]
         if scale is None:
             scale = math.sqrt(1.0 / self.base_ref_el.volume())
-        elif isinstance(scale, str):
-            scale = scale.lower()
-            if scale == "orthonormal":
-                scale = math.sqrt(1.0 / ref_el.volume())
-            elif scale == "l2 piola":
-                scale = 1.0 / ref_el.volume()
         self.scale = scale
         self.continuity = "C0" if variant == "bubble" else None
         self.recurrence_order = 2
         self._dmats_cache = {}
         self._cell_node_map_cache = {}
 
+    def get_scale(self, cell=0):
+        scale = self.scale
+        if isinstance(scale, str):
+            sd = self.ref_el.get_spatial_dimension()
+            vol = self.ref_el.volume_of_subcomplex(sd, cell)
+            scale = scale.lower()
+            if scale == "orthonormal":
+                scale = math.sqrt(1.0 / vol)
+            elif scale == "l2 piola":
+                scale = 1.0 / vol
+        return scale
+
     def get_num_members(self, n):
         return polynomial_dimension(self.ref_el, n, self.continuity)
 
@@ -317,8 +320,11 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         ref_pts = apply_mapping(A, b, numpy.transpose(pts))
         Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
         sd = self.ref_el.get_spatial_dimension()
+
+        # Always return 1 for n=0 to make regression tests pass
+        scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell)
         phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
-                                 self.scale, variant=self.variant)
+                                 scale, variant=self.variant)
         if self.continuity == "C0":
             phi = C0_basis(sd, n, phi)
 
@@ -494,7 +500,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         Jinv = A[0, 0] if direction is None else numpy.dot(A, direction)
         xs = apply_mapping(A, b, numpy.transpose(pts)).T
         results = {}
-        scale = self.scale * numpy.sqrt(2 * numpy.arange(n+1) + 1)
+        scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
         for k in range(order+1):
             v = numpy.zeros((n + 1, len(xs)), "d")
             if n >= k:

FIAT/functional.py

@@ -457,18 +457,17 @@ def __init__(self, ref_el, Q, f_at_qpts):
         nqp = len(qpts)
         dpts = qpts
         self.dpts = dpts
+        sd = ref_el.get_spatial_dimension()
 
         assert len(f_at_qpts.shape) == 2
-        assert f_at_qpts.shape[0] == 2
+        assert f_at_qpts.shape[0] == sd
         assert f_at_qpts.shape[1] == nqp
 
-        sd = ref_el.get_spatial_dimension()
-
         dpt_dict = OrderedDict()
 
         alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)]
         for q, pt in enumerate(dpts):
-            dpt_dict[tuple(pt)] = [(qwts[q]*f_at_qpts[i, q], alphas[j], (i, j)) for i in range(2) for j in range(2)]
+            dpt_dict[tuple(pt)] = [(qwts[q]*f_at_qpts[i, q], alphas[j], (i, j)) for i in range(sd) for j in range(sd)]
 
         super().__init__(ref_el, tuple(), {}, dpt_dict,
                          "IntegralMomentOfDivergence")

FIAT/johnson_mercier.py

@@ -0,0 +1,63 @@
+from FIAT import finite_element, dual_set, macro, polynomial_set
+from FIAT.functional import IntegralMoment, FrobeniusIntegralMoment
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.quadrature_schemes import create_quadrature
+import numpy
+
+
+class JohnsonMercierDualSet(dual_set.DualSet):
+    def __init__(self, ref_complex, degree, variant=None):
+        if degree != 1:
+            raise ValueError("Johnson-Mercier only defined for degree=1")
+        if variant is not None:
+            raise ValueError(f"Johnson-Mercier does not have the {variant} variant")
+        ref_el = ref_complex.get_parent()
+        top = ref_el.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+        nodes = []
+
+        # Face dofs: bidirectional (nn and nt) Legendre moments
+        dim = sd - 1
+        ref_facet = ref_el.construct_subelement(dim)
+        Qref = create_quadrature(ref_facet, 2*degree)
+        P = polynomial_set.ONPolynomialSet(ref_facet, degree)
+        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_tangents(dim, facet)
+            normal = ref_el.compute_normal(facet)
+            normal /= numpy.linalg.norm(normal)
+            scaled_normal = normal * Jdet
+            uvecs = (scaled_normal, *tangents)
+            comps = [numpy.outer(normal, uvec) for uvec in uvecs]
+            nodes.extend(FrobeniusIntegralMoment(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)))
+
+        cur = len(nodes)
+        if variant is None:
+            # Interior dofs: moments for each independent component
+            Q = create_quadrature(ref_complex, 2*degree-1)
+            P = polynomial_set.ONPolynomialSet(ref_el, degree-1)
+            phis = P.tabulate(Q.get_points())[(0,) * sd]
+            nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j))
+                         for j in range(sd) for i in range(j+1) for phi in phis)
+
+        entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super(JohnsonMercierDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+
+class JohnsonMercier(finite_element.CiarletElement):
+    """The Johnson-Mercier finite element."""
+
+    def __init__(self, ref_el, degree=1, variant=None):
+        ref_complex = macro.AlfeldSplit(ref_el)
+        poly_set = macro.HDivSymPolynomialSet(ref_complex, degree)
+        dual = JohnsonMercierDualSet(ref_complex, degree, variant=variant)
+        mapping = "double contravariant piola"
+        super(JohnsonMercier, self).__init__(poly_set, dual, degree,
+                                             mapping=mapping)

test/unit/test_johnson_mercier.py

@@ -0,0 +1,57 @@
+import pytest
+import numpy
+
+from FIAT import JohnsonMercier, Nedelec
+from FIAT.reference_element import ufc_simplex
+from FIAT.quadrature_schemes import create_quadrature
+
+
+@pytest.fixture(params=("T-ref", "T-phys", "S-ref", "S-phys"))
+def cell(request):
+    cell, deform = request.param.split("-")
+    dim = {"T": 2, "S": 3}[cell]
+    K = ufc_simplex(dim)
+    if deform == "phys":
+        if dim == 2:
+            K.vertices = ((0.0, 0.0), (2.0, 0.1), (0.0, 1.0))
+        else:
+            K.vertices = ((0, 0, 0), (1., 0.1, -0.37),
+                          (0.01, 0.987, -.23), (-0.1, -0.2, 1.38))
+    return K
+
+
+def test_johnson_mercier_divergence_rigid_body_motions(cell):
+    # test that the divergence of interior JM basis functions is orthogonal to
+    # the rigid-body motions
+    degree = 1
+    variant = None
+    sd = cell.get_spatial_dimension()
+    JM = JohnsonMercier(cell, degree, variant=variant)
+
+    ref_complex = JM.get_reference_complex()
+    Q = create_quadrature(ref_complex, 2*(degree)-1)
+    qpts, qwts = Q.get_points(), Q.get_weights()
+
+    tab = JM.tabulate(1, qpts)
+    div = sum(tab[alpha][:, alpha.index(1), :, :] for alpha in tab if sum(alpha) == 1)
+
+    # construct rigid body motions
+    N1 = Nedelec(cell, 1)
+    N1_at_qpts = N1.tabulate(0, qpts)
+    rbms = N1_at_qpts[(0,)*sd]
+    ells = rbms * qwts[None, None, :]
+
+    edofs = JM.entity_dofs()
+    idofs = edofs[sd][0]
+    L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2)))
+    assert numpy.allclose(L[idofs], 0)
+
+    if variant == "divergence":
+        edofs = JM.entity_dofs()
+        cdofs = []
+        for entity in edofs[sd-1]:
+            cdofs.extend(edofs[sd-1][entity][:sd])
+        D = L[cdofs]
+        M = numpy.tensordot(rbms, ells, axes=((1, 2), (1, 2)))
+        X = numpy.linalg.solve(M, D.T)
+        assert numpy.allclose(numpy.tensordot(X, rbms, axes=(0, 0)), div[cdofs])