Reduced Johnson-Mercier

FIAT/__init__.py

@@ -27,7 +27,7 @@
 from FIAT.arnold_qin import ArnoldQin
 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div
 from FIAT.christiansen_hu import ChristiansenHu
-from FIAT.johnson_mercier import JohnsonMercier
+from FIAT.johnson_mercier import JohnsonMercier, ReducedJohnsonMercier
 from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
 from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace
 from FIAT.SminusDiv import TrimmedSerendipityDiv
@@ -104,6 +104,7 @@
                       "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1,
                       "Guzman-Neilan H1(div)": GuzmanNeilanH1div,
                       "Johnson-Mercier": JohnsonMercier,
+                      "Reduced Johnson-Mercier": ReducedJohnsonMercier,
                       "Lagrange": Lagrange,
                       "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen,
                       "Gauss-Lobatto-Legendre": GaussLobattoLegendre,

FIAT/bernardi_raugel.py

@@ -116,6 +116,9 @@ def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None
                         Q = Qt[f]
                         phi = ft_at_qpts
                         udir = thats[f][i-1]
+                        if sd == 3:
+                            udir = numpy.cross(perp(*thats[f]), udir)
+
                     nodes.append(FrobeniusIntegralMoment(ref_el, Q, numpy.outer(udir, phi)))
                     entity_ids[sd-1][f].extend(range(cur, len(nodes)))
         super().__init__(nodes, ref_el, entity_ids)

FIAT/functional.py

@@ -459,6 +459,7 @@ def __init__(self, ref_el, Q, f_at_qpts):
         self.Q = Q
 
         sd = ref_el.get_spatial_dimension()
+        shp = f_at_qpts.shape[1:-1] + (sd,)
 
         points = Q.get_points()
         self.dpts = points
@@ -468,7 +469,7 @@ def __init__(self, ref_el, Q, f_at_qpts):
         dpt_dict = {tuple(pt): [(wt, alphas[i], (i,)) for i in range(sd)]
                     for pt, wt in zip(points, weights)}
 
-        super().__init__(ref_el, tuple(), {}, dpt_dict,
+        super().__init__(ref_el, shp, {}, dpt_dict,
                          "IntegralMomentOfDivergence")
 
 
@@ -492,7 +493,7 @@ def __init__(self, ref_el, Q, f_at_qpts):
         dpt_dict = {tuple(pt): [(wt[i], alphas[j], (i, j)) for i, j in numpy.ndindex(shp)]
                     for pt, wt in zip(points, weights)}
 
-        super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence")
+        super().__init__(ref_el, shp, {}, dpt_dict, "IntegralMomentOfDivergence")
 
 
 class PointNormalEvaluation(Functional):

FIAT/johnson_mercier.py

@@ -1,7 +1,9 @@
 from FIAT import finite_element, dual_set, macro, polynomial_set
 from FIAT.check_format_variant import parse_quadrature_scheme
-from FIAT.functional import TensorBidirectionalIntegralMoment
+from FIAT.functional import (FrobeniusIntegralMoment, IntegralMomentOfTensorDivergence,
+                             TensorBidirectionalIntegralMoment)
 from FIAT.quadrature import FacetQuadratureRule
+from FIAT.nedelec_second_kind import NedelecSecondKind
 import numpy
 
 
@@ -28,7 +30,12 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             cur = len(nodes)
             Q = FacetQuadratureRule(ref_el, dim, f, Qref, avg=True)
             thats = ref_el.compute_tangents(dim, f)
-            nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats)
+            if sd == 2:
+                nhat = numpy.dot(R, *thats)
+            else:
+                nhat = numpy.cross(*thats)
+                thats = numpy.cross(nhat[None, :], thats, axis=1)
+
             nodes.extend(TensorBidirectionalIntegralMoment(ref_el, nhat, comp, Q, phi)
                          for phi in phis for comp in (nhat, *thats))
             entity_ids[dim][f].extend(range(cur, len(nodes)))
@@ -57,3 +64,124 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None):
         formdegree = ref_el.get_spatial_dimension() - 1
         mapping = "double contravariant piola"
         super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
+
+
+class DG0Alfeld(finite_element.CiarletElement):
+
+    def __init__(self, ref_complex, P1=None, quad_scheme=None):
+        ref_el = ref_complex.get_parent()
+        sd = ref_el.get_spatial_dimension()
+        if P1 is None:
+            if sd == 1:
+                P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(sd,))
+            else:
+                P1 = NedelecSecondKind(ref_el, 1)
+
+        Q = parse_quadrature_scheme(ref_complex, P1.degree(), quad_scheme)
+        phis = P1.tabulate(0, Q.get_points())[(0,)*sd]
+        nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis]
+        dual = dual_set.DualSet(nodes, ref_el, P1.entity_dofs())
+
+        poly_set = polynomial_set.ONPolynomialSet(ref_complex, 0, shape=(sd,))
+        formdegree = sd
+        super().__init__(poly_set, dual, poly_set.degree, formdegree)
+
+
+def rbm_complement(ref_complex):
+    """Constructs a basis for the complement of the rigid body motions over P0(Alfeld)."""
+    P0 = DG0Alfeld(ref_complex)
+    P0 = DG0Alfeld(ref_complex, P1=P0)
+
+    entity_ids = P0.entity_dofs()
+    ids = []
+    for entity in entity_ids[1]:
+        ids.extend(entity_ids[1][entity][1:])
+    return P0.get_nodal_basis().take(ids)
+
+
+class ReducedJohnsonMercierDualSet(dual_set.DualSet):
+    def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
+        if degree != 1:
+            raise ValueError("Reduced Johnson-Mercier only defined for degree=1")
+        if variant is not None:
+            raise ValueError(f"Reduced Johnson-Mercier does not have the {variant} variant")
+        ref_el = ref_complex.get_parent()
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+        nodes = []
+
+        ref_facet = ref_el.get_facet_element()
+        Q = parse_quadrature_scheme(ref_facet, degree+1, quad_scheme)
+        P1 = polynomial_set.ONPolynomialSet(ref_facet, degree, scale="orthonormal")
+        P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)]
+        dimP1 = len(P1)*(sd-1)
+        dimNed1 = dimP1 // 2
+        if sd == 3:
+            # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)]
+            RT_at_qpts = numpy.zeros((dimP1, sd-1, P1_at_qpts.shape[-1]))
+            RT_at_qpts[0, 0] = P1_at_qpts[0]
+            RT_at_qpts[1, 1] = P1_at_qpts[0]
+            RT_at_qpts[2, 0] = P1_at_qpts[1]
+            RT_at_qpts[2, 1] = P1_at_qpts[2]
+            # Basis for the complement of RT [(y, x), (x, -y), (y, -x)]
+            RT_at_qpts[3, 0] = P1_at_qpts[2]
+            RT_at_qpts[3, 1] = P1_at_qpts[1]
+            RT_at_qpts[4, 0] = P1_at_qpts[1]
+            RT_at_qpts[4, 1] = -P1_at_qpts[2]
+            RT_at_qpts[5, 0] = P1_at_qpts[2]
+            RT_at_qpts[5, 1] = -P1_at_qpts[1]
+        else:
+            RT_at_qpts = P1_at_qpts[:, None, :]
+
+        # Facet degrees of freedom
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            n = ref_el.compute_scaled_normal(f)
+            Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True)
+            Jf = Qf.jacobian()
+            # Normal-normal moments against P1
+            nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n, n, Qf, phi) for phi in P1_at_qpts)
+            # Normal-tangential moments against n x RT
+            phis = numpy.tensordot(Jf, RT_at_qpts[:dimNed1].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2)
+            if sd == 3:
+                phis = numpy.cross(n[None, :, None], phis, axis=1)
+            phis = phis[:, :, None, :] * n[None, None, :, None]
+            nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis)
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        # Facet constraints
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            n = ref_el.compute_scaled_normal(f)
+            Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True)
+            Jf = Qf.jacobian()
+            # Normal-tangential moments against n x (P1 \ RT)
+            phis = numpy.tensordot(Jf, RT_at_qpts[dimNed1:].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2)
+            if sd == 3:
+                phis = numpy.cross(n[None, :, None], phis, axis=1)
+            phis = phis[:, :, None, :] * n[None, None, :, None]
+            nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis)
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        # Interior constraints: moments of divergence against (P1 \ Nedelec)
+        ref_complex = macro.AlfeldSplit(ref_el)
+        Q = parse_quadrature_scheme(ref_complex, 2*(degree-1))
+        comp_space = rbm_complement(ref_complex)
+        phis = comp_space.tabulate(Q.get_points())[(0,)*sd]
+        cur = len(nodes)
+        nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis)
+        entity_ids[sd][0].extend(range(cur, len(nodes)))
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class ReducedJohnsonMercier(finite_element.CiarletElement):
+    """The Reduced Johnson-Mercier finite element."""
+
+    def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None):
+        ref_complex = macro.AlfeldSplit(ref_el)
+        poly_set = macro.HDivSymPolynomialSet(ref_complex, degree)
+        dual = ReducedJohnsonMercierDualSet(ref_complex, degree, variant=variant, quad_scheme=quad_scheme)
+        formdegree = ref_el.get_spatial_dimension() - 1
+        mapping = "double contravariant piola"
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)

FIAT/mardal_tai_winther.py

@@ -8,26 +8,52 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 
+import numpy
 from FIAT import dual_set, expansions, finite_element, polynomial_set
-from FIAT.functional import (IntegralMomentOfNormalEvaluation,
-                             IntegralMomentOfTangentialEvaluation,
-                             IntegralLegendreNormalMoment,
-                             IntegralMomentOfDivergence)
-
+from FIAT.functional import FrobeniusIntegralMoment
+from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
 
 
-def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg):
+def curl(table_u):
+    grad_u = {alpha.index(1): table_u[alpha] for alpha in table_u if sum(alpha) == 1}
+    d = len(grad_u)
+    if d == 2:
+        curl_u = [grad_u[1], -grad_u[0]]
+        return numpy.concatenate(curl_u, axis=1)
+    else:
+        indices = ((i, j) for i in reversed(range(d)) for j in reversed(range(i+1, d)))
+        curl_u = [((-1)**k) * (grad_u[j][:, i, :] - grad_u[i][:, j, :]) for k, (i, j) in enumerate(indices)]
+        return numpy.transpose(curl_u, (1, 0, 2))
+
+
+def MardalTaiWintherSpace(ref_el):
+    """Construct the MTW space [P1]^d + curl(B [P1]^d)."""
     sd = ref_el.get_spatial_dimension()
-    P = polynomial_set.ONPolynomialSet(ref_el, stop_deg)
-    Q = create_quadrature(ref_el, comp_deg + stop_deg)
+    # Polynomials of degree sd+1
+    degree = sd + 1
+    Pk = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), scale="orthonormal")
+
+    # Grab [P1]^d from [Pk]^d
+    dimP1 = expansions.polynomial_dimension(ref_el, 1)
+    dimPk = expansions.polynomial_dimension(ref_el, degree)
+    ids = [i+dimPk*j for i in range(dimP1) for j in range(sd)]
+    P1 = Pk.take(ids)
+    # Project curl(B [P1]^d) into [Pk]^d
+    BP1 = polynomial_set.make_bubbles(ref_el, degree+1, shape=((sd*(sd-1))//2,))
 
-    dim0 = expansions.polynomial_dimension(ref_el, start_deg-1)
-    dim1 = expansions.polynomial_dimension(ref_el, stop_deg)
-    indices = list(range(dim0, dim1))
-    phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd]
-    for phi in phis:
-        yield IntegralMomentOfDivergence(ref_el, Q, phi)
+    Q = create_quadrature(ref_el, 2*degree)
+    qpts = Q.get_points()
+    qwts = Q.get_weights()
+    Pk_at_qpts = Pk.tabulate(qpts)
+    BP1_at_qpts = BP1.tabulate(qpts, 1)
+
+    inner = lambda u, v, qwts: numpy.tensordot(u, numpy.multiply(v, qwts), axes=(range(1, u.ndim),)*2)
+    C = inner(curl(BP1_at_qpts), Pk_at_qpts[(0,)*sd], qwts)
+    coeffs = numpy.tensordot(C, Pk.get_coeffs(), axes=(1, 0))
+    curlBP1 = polynomial_set.PolynomialSet(ref_el, degree, degree, Pk.get_expansion_set(), coeffs)
+
+    return polynomial_set.polynomial_set_union_normalized(P1, curlBP1)
 
 
 class MardalTaiWintherDual(dual_set.DualSet):
@@ -36,70 +62,55 @@ def __init__(self, ref_el, degree):
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
 
-        if sd != 2:
-            raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2.")
+        if sd not in (2, 3):
+            raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2 and 3.")
 
-        if degree != 3:
-            raise ValueError("Mardal-Tai-Winther elements are only defined for degree 3.")
+        if degree != sd+1:
+            raise ValueError("Mardal-Tai-Winther elements are only defined for degree = dim+1.")
 
         entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
         nodes = []
 
-        # no vertex dofs
-
         # On each facet, let n be its normal.  We need to integrate
-        # u.n and u.t against the first Legendre polynomial (constant)
-        # and u.n against the second (linear).
-        facet = ref_el.get_facet_element()
-        # Facet nodes are \int_F v.n p ds where p \in P_{q-1}
-        # degree is q - 1
-        Q = create_quadrature(facet, degree+1)
-        Pq = polynomial_set.ONPolynomialSet(facet, 1)
-        phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)]
-        for f in sorted(top[sd-1]):
-            cur = len(nodes)
-            nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[0], f))
-            nodes.append(IntegralMomentOfTangentialEvaluation(ref_el, Q, phis[0], f))
-            nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[1], f))
-            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+        # u.n against a Dubiner basis for P1
+        # and u x n against a basis for lowest-order RT.
+        ref_facet = ref_el.get_facet_element()
+        Q = create_quadrature(ref_facet, degree+1)
+        P1 = polynomial_set.ONPolynomialSet(ref_facet, 1)
+        P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)]
+        if sd == 2:
+            # For 2D just take the constant
+            RT_at_qpts = P1_at_qpts[:1, None, :]
+        else:
+            # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)]
+            RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1]))
+            RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :]
+            RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :]
+            RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :]
+            RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :]
 
-        # Generate constraint nodes on the cell and facets
-        # * div(v) must be constant on the cell.  Since v is a cubic and
-        #   div(v) is quadratic, we need the integral of div(v) against the
-        #   linear and quadratic Dubiner polynomials to vanish.
-        #   There are two linear and three quadratics, so these are five
-        #   constraints
-        # * v.n must be linear on each facet.  Since v.n is cubic, we need
-        #   the integral of v.n against the cubic and quadratic Legendre
-        #   polynomial to vanish on each facet.
-
-        # So we introduce functionals whose kernel describes this property,
-        # as described in the FIAT paper.
-        start_order = 2
-        stop_order = 3
-        qdegree = degree + stop_order
         for f in sorted(top[sd-1]):
             cur = len(nodes)
-            nodes.extend(IntegralLegendreNormalMoment(ref_el, f, order, qdegree)
-                         for order in range(start_order, stop_order+1))
+            n = ref_el.compute_scaled_normal(f)
+            Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True)
+            # Normal moments against P1
+            nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts)
+            # Map the RT basis into the facet
+            Jf = Qf.jacobian()
+            phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2)
+            if sd == 3:
+                # Moments against cross(n, RT)
+                phis = numpy.cross(n[None, :, None], phis, axis=1)
+            nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis)
             entity_ids[sd-1][f].extend(range(cur, len(nodes)))
-
-        cur = len(nodes)
-        nodes.extend(DivergenceDubinerMoments(ref_el, start_order-1, stop_order-1, degree))
-        entity_ids[sd][0].extend(range(cur, len(nodes)))
-
         super().__init__(nodes, ref_el, entity_ids)
 
 
 class MardalTaiWinther(finite_element.CiarletElement):
     """The definition of the Mardal-Tai-Winther element.
     """
     def __init__(self, ref_el, degree=3):
-        sd = ref_el.get_spatial_dimension()
-        assert degree == 3, "Only defined for degree 3"
-        assert sd == 2, "Only defined for dimension 2"
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd,))
         dual = MardalTaiWintherDual(ref_el, degree)
-        formdegree = sd-1
-        mapping = "contravariant piola"
-        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
+        poly_set = MardalTaiWintherSpace(ref_el)
+        formdegree = ref_el.get_spatial_dimension() - 1
+        super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

finat/__init__.py

@@ -28,7 +28,7 @@
 from .guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanBubble, GuzmanNeilanH1div  # noqa: F401
 from .powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12  # noqa: F401
 from .hermite import Hermite                                       # noqa: F401
-from .johnson_mercier import JohnsonMercier                        # noqa: F401
+from .johnson_mercier import JohnsonMercier, ReducedJohnsonMercier  # noqa: F401
 from .mtw import MardalTaiWinther                                  # noqa: F401
 from .morley import Morley                                         # noqa: F401
 from .walkington import Walkington                                 # noqa: F401