Some Zany H(div) elements (#56)

Implement Arnold-Winther and Mardal-Tai-Winther elements

Co-authored-by: Patrick Farrell <[email protected]>
Co-authored-by: Francis Aznaran <[email protected]>

FIAT/__init__.py

@@ -30,6 +30,9 @@
 from FIAT.crouzeix_raviart import CrouzeixRaviart
 from FIAT.regge import Regge
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
+from FIAT.arnold_winther import ArnoldWinther
+from FIAT.arnold_winther import ArnoldWintherNC
+from FIAT.mardal_tai_winther import MardalTaiWinther
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.tensor_product import TensorProductElement
 from FIAT.enriched import EnrichedElement
@@ -77,7 +80,10 @@
                       "TensorProductElement": TensorProductElement,
                       "BrokenElement": DiscontinuousElement,
                       "HDiv Trace": HDivTrace,
-                      "Hellan-Herrmann-Johnson": HellanHerrmannJohnson}
+                      "Hellan-Herrmann-Johnson": HellanHerrmannJohnson,
+                      "Conforming Arnold-Winther": ArnoldWinther,
+                      "Nonconforming Arnold-Winther": ArnoldWintherNC,
+                      "Mardal-Tai-Winther": MardalTaiWinther}
 
 # List of extra elements
 extra_elements = {"P0": P0,

FIAT/arnold_winther.py

@@ -0,0 +1,166 @@
+# -*- coding: utf-8 -*-
+"""Implementation of the Arnold-Winther finite elements."""
+
+# Copyright (C) 2020 by Robert C. Kirby (Baylor University)
+# Modified by Francis Aznaran (Oxford University)
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+
+from FIAT.finite_element import CiarletElement
+from FIAT.dual_set import DualSet
+from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet
+from FIAT.functional import (
+    PointwiseInnerProductEvaluation as InnerProduct,
+    FrobeniusIntegralMoment as FIM,
+    IntegralMomentOfTensorDivergence,
+    IntegralLegendreNormalNormalMoment,
+    IntegralLegendreNormalTangentialMoment)
+
+from FIAT.quadrature import make_quadrature
+
+import numpy
+
+
+class ArnoldWintherNCDual(DualSet):
+    def __init__(self, cell, degree):
+        if not degree == 2:
+            raise ValueError("Nonconforming Arnold-Winther elements are"
+                             "only defined for degree 2.")
+        dofs = []
+        dof_ids = {}
+        dof_ids[0] = {0: [], 1: [], 2: []}
+        dof_ids[1] = {0: [], 1: [], 2: []}
+        dof_ids[2] = {0: []}
+
+        dof_cur = 0
+        # no vertex dofs
+        # proper edge dofs now (not the contraints)
+        # moments of normal . sigma against constants and linears.
+        for entity_id in range(3):                  # a triangle has 3 edges
+            for order in (0, 1):
+                dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
+                         IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
+            dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
+            dof_cur += 4
+
+        # internal dofs: constant moments of three unique components
+        Q = make_quadrature(cell, 2)
+
+        e1 = numpy.array([1.0, 0.0])              # euclidean basis 1
+        e2 = numpy.array([0.0, 1.0])              # euclidean basis 2
+        basis = [(e1, e1), (e1, e2), (e2, e2)]    # basis for symmetric matrices
+        for (v1, v2) in basis:
+            v1v2t = numpy.outer(v1, v2)
+            fatqp = numpy.zeros((2, 2, len(Q.pts)))
+            for i, y in enumerate(v1v2t):
+                for j, x in enumerate(y):
+                    for k in range(len(Q.pts)):
+                        fatqp[i, j, k] = x
+            dofs.append(FIM(cell, Q, fatqp))
+        dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
+        dof_cur += 3
+
+        # put the constraint dofs last.
+        for entity_id in range(3):
+            dof = IntegralLegendreNormalNormalMoment(cell, entity_id, 2, 6)
+            dofs.append(dof)
+            dof_ids[1][entity_id].append(dof_cur)
+            dof_cur += 1
+
+        super(ArnoldWintherNCDual, self).__init__(dofs, cell, dof_ids)
+
+
+class ArnoldWintherNC(CiarletElement):
+    """The definition of the nonconforming Arnold-Winther element.
+    """
+    def __init__(self, cell, degree):
+        assert degree == 2, "Only defined for degree 2"
+        Ps = ONSymTensorPolynomialSet(cell, degree)
+        Ls = ArnoldWintherNCDual(cell, degree)
+        mapping = "double contravariant piola"
+
+        super(ArnoldWintherNC, self).__init__(Ps, Ls, degree,
+                                              mapping=mapping)
+
+
+class ArnoldWintherDual(DualSet):
+    def __init__(self, cell, degree):
+        if not degree == 3:
+            raise ValueError("Arnold-Winther elements are"
+                             "only defined for degree 3.")
+        dofs = []
+        dof_ids = {}
+        dof_ids[0] = {0: [], 1: [], 2: []}
+        dof_ids[1] = {0: [], 1: [], 2: []}
+        dof_ids[2] = {0: []}
+
+        dof_cur = 0
+
+        # vertex dofs
+        vs = cell.get_vertices()
+        e1 = numpy.array([1.0, 0.0])
+        e2 = numpy.array([0.0, 1.0])
+        basis = [(e1, e1), (e1, e2), (e2, e2)]
+
+        dof_cur = 0
+
+        for entity_id in range(3):
+            node = tuple(vs[entity_id])
+            for (v1, v2) in basis:
+                dofs.append(InnerProduct(cell, v1, v2, node))
+            dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3))
+            dof_cur += 3
+
+        # edge dofs now
+        # moments of normal . sigma against constants and linears.
+        for entity_id in range(3):
+            for order in (0, 1):
+                dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
+                         IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
+            dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
+            dof_cur += 4
+
+        # internal dofs: constant moments of three unique components
+        Q = make_quadrature(cell, 3)
+
+        e1 = numpy.array([1.0, 0.0])              # euclidean basis 1
+        e2 = numpy.array([0.0, 1.0])              # euclidean basis 2
+        basis = [(e1, e1), (e1, e2), (e2, e2)]    # basis for symmetric matrices
+        for (v1, v2) in basis:
+            v1v2t = numpy.outer(v1, v2)
+            fatqp = numpy.zeros((2, 2, len(Q.pts)))
+            for k in range(len(Q.pts)):
+                fatqp[:, :, k] = v1v2t
+            dofs.append(FIM(cell, Q, fatqp))
+        dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
+        dof_cur += 3
+
+        # Constraint dofs
+
+        Q = make_quadrature(cell, 5)
+
+        onp = ONPolynomialSet(cell, 2, (2,))
+        pts = Q.get_points()
+        onpvals = onp.tabulate(pts)[0, 0]
+
+        for i in list(range(3, 6)) + list(range(9, 12)):
+            dofs.append(IntegralMomentOfTensorDivergence(cell, Q,
+                                                         onpvals[i, :, :]))
+
+        dof_ids[2][0] += list(range(dof_cur, dof_cur+6))
+
+        super(ArnoldWintherDual, self).__init__(dofs, cell, dof_ids)
+
+
+class ArnoldWinther(CiarletElement):
+    """The definition of the conforming Arnold-Winther element.
+    """
+    def __init__(self, cell, degree):
+        assert degree == 3, "Only defined for degree 3"
+        Ps = ONSymTensorPolynomialSet(cell, degree)
+        Ls = ArnoldWintherDual(cell, degree)
+        mapping = "double contravariant piola"
+        super(ArnoldWinther, self).__init__(Ps, Ls, degree, mapping=mapping)