Wu-Xu H3 non-conforming triangle

FIAT/__init__.py

@@ -62,6 +62,7 @@
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.hdiv_trace import HDivTrace
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen
+from FIAT.wuxu import WuXuH3NC, WuXuRobustH3NC
 from FIAT.walkington import Walkington
 from FIAT.histopolation import Histopolation
 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
@@ -88,6 +89,8 @@
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,
+                      "Nonconforming Wu-Xu": WuXuH3NC,
+                      "Nonconforming Robust Wu-Xu": WuXuRobustH3NC,
                       "Hsieh-Clough-Tocher": HsiehCloughTocher,
                       "QuadraticPowellSabin6": QuadraticPowellSabin6,
                       "QuadraticPowellSabin12": QuadraticPowellSabin12,

FIAT/wuxu.py

@@ -0,0 +1,192 @@
+# Copyright (C) 2022 Robert C. Kirby (Baylor University)
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+
+from FIAT import expansions, polynomial_set, dual_set, finite_element
+from FIAT.functional import (IntegralMomentOfDerivative,
+                             IntegralMomentOfBidirectionalDerivative,
+                             PointDerivative, PointEvaluation)
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.polynomial_set import mis
+from FIAT.bubble import Bubble
+from FIAT.lagrange import Lagrange
+import numpy
+
+polydim = expansions.polynomial_dimension
+
+
+def WuXuRobustH3NCSpace(ref_el):
+    """Constructs a basis for the the Wu Xu H^3 nonconforming space
+    P^{(3,2)}(T) = P_3(T) + b_T P_1(T) + b_T^2 P_1(T),
+    where b_T is the standard cubic bubble."""
+
+    sd = ref_el.get_spatial_dimension()
+    assert sd == 2
+
+    em_deg = 7
+
+    # Unfortunately,  b_T^2 P_1 has degree 7 (cubic squared times a linear)
+    # so we need a high embedded degree!
+    p7 = polynomial_set.ONPolynomialSet(ref_el, 7)
+
+    dimp1 = polydim(ref_el, 1)
+    dimp3 = polydim(ref_el, 3)
+    dimp7 = polydim(ref_el, 7)
+
+    # Here's the first bit we'll work with.  It's already expressed in terms
+    # of the ON basis for P7, so we're golden.
+    p3fromp7 = p7.take(list(range(dimp3)))
+
+    # Rather than creating the barycentric coordinates ourself, let's
+    # reuse the existing bubble functionality
+    bT = Bubble(ref_el, 3)
+    p1 = Lagrange(ref_el, 1)
+
+    # next, we'll have to project b_T P1 and b_T^2 P1 onto P^7
+    Q = create_quadrature(ref_el, 14)
+    Qpts = numpy.array(Q.get_points())
+    Qwts = numpy.array(Q.get_weights())
+
+    zero_index = tuple([0 for i in range(sd)])
+
+    # it's just one bubble function: let's get a 1d array!
+    bT_at_qpts = bT.tabulate(0, Qpts)[zero_index][0, :]
+    p1_at_qpts = p1.tabulate(0, Qpts)[zero_index]
+
+    # Note: difference in signature because bT, p1 are FE and p7 is a
+    # polynomial set
+    p7_at_qpts = p7.tabulate(Qpts)[zero_index]
+
+    bubble_coeffs = numpy.zeros((6, dimp7), "d")
+
+    # first three: bT P1, last three will be bT^2 P1
+    foo = bT_at_qpts * p1_at_qpts * Qwts
+    bubble_coeffs[:dimp1, :] = numpy.dot(foo, p7_at_qpts.T)
+
+    foo = bT_at_qpts * foo
+    bubble_coeffs[dimp1:2*dimp1, :] = numpy.dot(foo, p7_at_qpts.T)
+
+    bubbles = polynomial_set.PolynomialSet(ref_el, 3, em_deg,
+                                           p7.get_expansion_set(),
+                                           bubble_coeffs)
+
+    return polynomial_set.polynomial_set_union_normalized(p3fromp7, bubbles)
+
+
+def WuXuH3NCSpace(ref_el):
+    """Constructs a basis for the the Wu Xu H^3 nonconforming space
+    P(T) = P_3(T) + b_T P_1(T),
+    where b_T is the standard cubic bubble."""
+
+    assert ref_el.get_spatial_dimension() == 2
+
+    em_deg = 4
+    p4 = polynomial_set.ONPolynomialSet(ref_el, em_deg)
+
+    # Here's the first bit we'll work with.  It's already expressed in terms
+    # of the ON basis for P4, so we're golden.
+    dimp3 = polydim(ref_el, 3)
+    p3fromp4 = p4.take(list(range(dimp3)))
+
+    # Rather than creating the barycentric coordinates ourself, let's
+    # reuse the existing bubble functionality
+    bT = Bubble(ref_el, 4)
+
+    return polynomial_set.polynomial_set_union_normalized(p3fromp4, bT.get_nodal_basis())
+
+
+class WuXuRobustH3NCDualSet(dual_set.DualSet):
+    """Dual basis for WuXu H3 nonconforming element consisting of
+    vertex values and gradients and first and second normals at edge midpoints."""
+
+    def __init__(self, ref_el, degree):
+        sd = ref_el.get_spatial_dimension()
+        assert sd == 2
+        top = ref_el.get_topology()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
+
+        # jet at each vertex
+        verts = ref_el.get_vertices()
+        for v in sorted(top[0]):
+            cur = len(nodes)
+            nodes.append(PointEvaluation(ref_el, verts[v]))
+            nodes.extend(PointDerivative(ref_el, verts[v], alpha) for alpha in mis(sd, 1))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # average of first and second normal derivative along each edge
+        Q_ref = create_quadrature(ref_el.construct_subelement(1), degree-1)
+        for e in sorted(top[1]):
+            n = ref_el.compute_normal(e)
+            Q = FacetQuadratureRule(ref_el, 1, e, Q_ref)
+            f = numpy.full(Q.get_weights().shape, 1/Q.jacobian_determinant())
+            cur = len(nodes)
+            nodes.append(IntegralMomentOfDerivative(ref_el, n, Q, f))
+            nodes.append(IntegralMomentOfBidirectionalDerivative(ref_el, Q, f, n, n))
+            entity_ids[1][e].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class WuXuH3NCDualSet(dual_set.DualSet):
+    """Dual basis for WuXu H3 nonconforming element consisting of
+    vertex values and gradients and second normals at edge midpoints."""
+
+    def __init__(self, ref_el, degree):
+        sd = ref_el.get_spatial_dimension()
+        assert sd == 2
+        top = ref_el.get_topology()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
+
+        # jet at each vertex
+        verts = ref_el.get_vertices()
+        for v in sorted(top[0]):
+            cur = len(nodes)
+            nodes.append(PointEvaluation(ref_el, verts[v]))
+            nodes.extend(PointDerivative(ref_el, verts[v], alpha) for alpha in mis(sd, 1))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # average of second normal derivative along each edge
+        Q_ref = create_quadrature(ref_el.construct_subelement(1), degree-2)
+        for e in sorted(top[1]):
+            n = ref_el.compute_normal(e)
+            Q = FacetQuadratureRule(ref_el, 1, e, Q_ref)
+            f = numpy.full(Q.get_weights().shape, 1/Q.jacobian_determinant())
+            cur = len(nodes)
+            nodes.append(IntegralMomentOfBidirectionalDerivative(ref_el, Q, f, n, n))
+            entity_ids[1][e].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class WuXuRobustH3NC(finite_element.CiarletElement):
+    """The Wu-Xu robust H3 nonconforming finite element"""
+    def __init__(self, ref_el, degree=7):
+        poly_set = WuXuRobustH3NCSpace(ref_el)
+        assert degree == poly_set.degree
+        dual = WuXuRobustH3NCDualSet(ref_el, degree)
+        super().__init__(poly_set, dual, degree)
+
+
+class WuXuH3NC(finite_element.CiarletElement):
+    """The Wu-Xu H3 nonconforming finite element"""
+    def __init__(self, ref_el, degree=4):
+        poly_set = WuXuH3NCSpace(ref_el)
+        assert degree == poly_set.degree
+        dual = WuXuH3NCDualSet(ref_el, degree)
+        super().__init__(poly_set, dual, degree)

finat/__init__.py

@@ -18,6 +18,7 @@
 from .aw import ArnoldWinther, ArnoldWintherNC                     # noqa: F401
 from .hz import HuZhang                                            # noqa: F401
 from .bell import Bell                                             # noqa: F401
+from .wuxu import WuXuH3NC, WuXuRobustH3NC                         # noqa: F401
 from .bernardi_raugel import BernardiRaugel, BernardiRaugelBubble  # noqa: F401
 from .hct import HsiehCloughTocher, ReducedHsiehCloughTocher       # noqa: F401
 from .arnold_qin import ArnoldQin, ReducedArnoldQin                # noqa: F401

finat/citations.py

@@ -260,6 +260,7 @@ def cite(*args, **kwargs):
 }
 """)
     petsctools.add_citation("Walkington2010", """
+@article{Walkington2010,
 author = {Walkington, Noel J.},
 title = {{A \\$C^1\\$ Tetrahedral Finite Element without Edge Degrees of Freedom}},
 journal = {SIAM Journal on Numerical Analysis},
@@ -268,4 +269,14 @@ def cite(*args, **kwargs):
 pages = {330-342},
 year = {2014},
 doi = {10.1137/130912013},
+}""")
+    petsctools.add_citation("WuXu2019", """
+@article{WuXu2019,
+title={Nonconforming finite element spaces for $2m^{\\mathrm{th}}$ order partial differential equations on $\\mathbb{R}^n$ simplicial grids when $\\mathbb{m} = \\mathbb{n} + 1$},
+author={Wu, Shuonan and Xu, Jinchao},
+journal={Mathematics of Computation},
+volume={88},
+number={316},
+pages={531--551},
+year={2019}
 }""")

finat/element_factory.py

@@ -91,6 +91,8 @@
                       "Hu-Zhang": finat.HuZhang,
                       "Mardal-Tai-Winther": finat.MardalTaiWinther,
                       "Walkington": finat.Walkington,
+                      "Nonconforming Wu-Xu": finat.WuXuH3NC,
+                      "Nonconforming Robust Wu-Xu": finat.WuXuRobustH3NC,
                       # These require special treatment
                       "Q": None,
                       "DQ": None,

finat/ufl/elementlist.py

@@ -22,7 +22,7 @@
 from numpy import asarray
 
 from ufl.cell import Cell, TensorProductCell
-from ufl.sobolevspace import H1, H2, L2, HCurl, HDiv, HDivDiv, HCurlDiv, HEin, HInf
+from ufl.sobolevspace import H1, H2, H3, L2, HCurl, HDiv, HDivDiv, HCurlDiv, HEin, HInf
 from ufl.utils.formatting import istr
 
 # List of valid elements
@@ -116,6 +116,8 @@ def show_elements():
 register_element("Argyris", "ARG", 0, H2, "custom", (5, None), ("triangle",))
 register_element("Bell", "BELL", 0, H2, "custom", (5, 5), ("triangle",))
 register_element("Morley", "MOR", 0, H2, "custom", (2, 2), simplices[1:])
+register_element("Nonconforming Wu-Xu", "WXnc", 0, H3, "custom", (4, 4), ("triangle",))
+register_element("Nonconforming Robust Wu-Xu", "WXncr", 0, H3, "custom", (7, 7), ("triangle",))
 
 # Macro elements
 register_element("QuadraticPowellSabin6", "PS6", 0, H2, "custom", (2, 2), ("triangle",))

finat/wuxu.py

@@ -0,0 +1,101 @@
+from finat.fiat_elements import ScalarFiatElement
+from finat.physically_mapped import identity, PhysicallyMappedElement
+from finat.argyris import _vertex_transform
+from finat.citations import cite
+from gem import ListTensor
+
+import FIAT
+import numpy
+
+
+class WuXuRobustH3NC(PhysicallyMappedElement, ScalarFiatElement):
+    def __init__(self, cell, degree=7):
+        if degree != 7:
+            raise ValueError("Degree must be 7 for robust Wu-Xu element")
+        cite("WuXu2019")
+        super().__init__(FIAT.WuXuRobustH3NC(cell))
+
+    def basis_transformation(self, coordinate_mapping):
+        return wuxu_transformation(self, coordinate_mapping)
+
+
+class WuXuH3NC(PhysicallyMappedElement, ScalarFiatElement):
+    def __init__(self, cell, degree=4):
+        if degree != 4:
+            raise ValueError("Degree must be 4 for the Wu-Xu element")
+        cite("WuXu2019")
+        super().__init__(FIAT.WuXuH3NC(cell))
+
+    def basis_transformation(self, coordinate_mapping):
+        return wuxu_transformation(self, coordinate_mapping)
+
+
+def hessian_transform(J):
+    return numpy.array(
+        [[J[0, 0] * J[0, 0], J[0, 0] * J[1, 0] + J[0, 0] * J[1, 0], J[1, 0] * J[1, 0]],
+         [J[0, 1] * J[0, 0], J[0, 1] * J[1, 0] + J[0, 0] * J[1, 1], J[1, 0] * J[1, 1]],
+         [J[0, 1] * J[0, 1], J[0, 1] * J[1, 1] + J[0, 1] * J[1, 1], J[1, 1] * J[1, 1]]],
+        dtype=object)
+
+
+def wuxu_transformation(self, coordinate_mapping):
+    top = self.cell.topology
+    sd = self.cell.get_spatial_dimension()
+    entity_ids = self._element.entity_dofs()
+
+    V = identity(self.space_dimension())
+    _vertex_transform(V, 1, self.cell, coordinate_mapping)
+
+    bary, = self.cell.make_points(sd, 0, sd+1)
+    J = coordinate_mapping.jacobian_at(bary)
+    Thetainv = hessian_transform(J)
+    J = numpy.array([[J[i, j] for j in range(sd)] for i in range(sd)])
+
+    ns = coordinate_mapping.physical_normals()
+    ts = coordinate_mapping.physical_tangents()
+    lens = coordinate_mapping.physical_edge_lengths()
+    nhats = coordinate_mapping.reference_normals()
+    thats = coordinate_mapping.normalized_reference_edge_tangents()
+
+    for e in top[1]:
+        v0, v1 = top[1][e]
+        vid0 = entity_ids[0][v0]
+        vid1 = entity_ids[0][v1]
+
+        G = numpy.array([[u[e, j] for j in range(sd)] for u in (ns, ts)])
+        Ghat = numpy.array([[u[e, j] for j in range(sd)] for u in (nhats, thats)])
+
+        if len(entity_ids[1][e]) > 1:
+            # first derivative moments
+            eid = entity_ids[1][e][0]
+            B1 = (Ghat @ J.T) @ G.T
+            alpha = B1[0, 1] / lens[e]
+            V[eid, eid] = B1[0, 0]
+            V[eid, vid0[0]] = -1*alpha
+            V[eid, vid1[0]] = alpha
+
+        # second derivative moments
+        eid = entity_ids[1][e][-1]
+        Gamma = hessian_transform(G)
+        Gammainvhat = hessian_transform(Ghat.T)
+        B2 = (Gammainvhat @ Thetainv) @ Gamma
+        beta = B2[0, 1:] @ G / lens[e]
+        V[eid, eid] = B2[0, 0]
+        V[eid, vid0[1:]] = -1*beta
+        V[eid, vid1[1:]] = beta
+
+    # Now let's fix the scaling.
+    h = coordinate_mapping.cell_size()
+
+    # This gets the vertex gradients
+    for v in top[0]:
+        vids = entity_ids[0][v][1:]
+        V[:, vids] *= 1 / h[v]
+
+    # this scales second derivative moments.  First should be ok.
+    for e in top[1]:
+        eid = entity_ids[1][e][-1]
+        he = (1/len(top[1][e])) * sum(h[v] for v in top[1][e])
+        V[:, eid] *= 1 / (he * he)
+
+    return ListTensor(V.T)

test/FIAT/unit/test_fiat.py

@@ -64,6 +64,7 @@
 from FIAT.nodal_enriched import NodalEnrichedElement
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen    # noqa: F401
 from FIAT.histopolation import Histopolation                    # noqa: F401
+from FIAT.wuxu import WuXuH3NC, WuXuRobustH3NC                  # noqa: F401
 
 P = Point()
 I = UFCInterval()  # noqa: E741
@@ -346,6 +347,8 @@ def __init__(self, a, b):
     "Argyris(T, 5, 'point')",
     "Argyris(T, 5, 'integral')",
     "Argyris(T, 6, 'integral')",
+    "WuXuH3NC(T, 4)",
+    "WuXuRobustH3NC(T, 7)",
     "CubicHermite(I)",
     "CubicHermite(T)",
     "CubicHermite(S)",

test/finat/test_mass_conditioning.py

@@ -15,6 +15,8 @@
     (2, finat.Argyris, 5, "point"),
     (2, finat.Argyris, 5, None),
     (2, finat.Argyris, 6, None),
+    (2, finat.WuXuH3NC, 4, None),
+    (2, finat.WuXuRobustH3NC, 7, None),
     (3, finat.Walkington, 5, None),
 ])
 def test_mass_scaling(scaled_ref_to_phys, sd, element, degree, variant):