Rckirby/ps (#78)

* Add Powell-Sabin elements

* Add test for PS akin to HCT test

---------

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

Author: rckirby

FIAT/__init__.py

@@ -31,6 +31,8 @@
 from FIAT.morley import Morley
 from FIAT.nedelec import Nedelec
 from FIAT.nedelec_second_kind import NedelecSecondKind
+from FIAT.powell_sabin import QuadraticPowellSabin6
+from FIAT.powell_sabin import QuadraticPowellSabin12
 from FIAT.hierarchical import Legendre, IntegratedLegendre
 from FIAT.P0 import P0
 from FIAT.raviart_thomas import RaviartThomas
@@ -93,6 +95,8 @@
                       "Regge": Regge,
                       "EnrichedElement": EnrichedElement,
                       "NodalEnrichedElement": NodalEnrichedElement,
+                      "QuadraticPowellSabin6": QuadraticPowellSabin6,
+                      "QuadraticPowellSabin12": QuadraticPowellSabin12,
                       "TensorProductElement": TensorProductElement,
                       "BrokenElement": DiscontinuousElement,
                       "HDiv Trace": HDivTrace,

FIAT/macro.py

@@ -3,9 +3,10 @@
 
 import numpy
 
-from FIAT.quadrature import FacetQuadratureRule, QuadratureRule
-from FIAT.reference_element import SimplicialComplex, lattice_iter, make_lattice
 from FIAT import expansions, polynomial_set
+from FIAT.quadrature import FacetQuadratureRule, QuadratureRule
+from FIAT.reference_element import (TRIANGLE, SimplicialComplex, lattice_iter,
+                                    make_lattice)
 
 
 def bary_to_xy(verts, bary, result=None):
@@ -320,6 +321,47 @@ def construct_subcomplex(self, dimension):
             return PowellSabinSplit(ref_el)
 
 
+class PowellSabin12Split(SplitSimplicialComplex):
+    """Splits a triangle (only!) by connecting each vertex to the opposite
+    edge midpoint and edge midpoints to each other.
+    """
+    def __init__(self, ref_el):
+        assert ref_el.get_shape() == TRIANGLE
+        verts = ref_el.get_vertices()
+        new_verts = list(verts)
+        new_verts.extend(
+            map(tuple, bary_to_xy(verts,
+                [(1/3, 1/3, 1/3),
+                 (1/2, 1/2, 0),
+                 (1/2, 0, 1/2),
+                 (0, 1/2, 1/2),
+                 (1/2, 1/4, 1/4),
+                 (1/4, 1/2, 1/4),
+                 (1/4, 1/4, 1/2)])))
+
+        edges = [(0, 4), (0, 7), (0, 5),
+                 (1, 4), (1, 8), (1, 6),
+                 (2, 5), (2, 9), (2, 6),
+                 (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9),
+                 (4, 7), (4, 8), (5, 7), (5, 9), (6, 8), (6, 9)]
+
+        super(PowellSabin12Split, self).__init__(
+            ref_el, tuple(new_verts), make_topology(2, len(new_verts), edges))
+
+    def construct_subcomplex(self, dimension):
+        """Constructs the reference subcomplex of the parent cell subentity
+        specified by subcomplex dimension.
+        """
+        if dimension == 2:
+            return self
+        elif dimension == 1:
+            return AlfeldSplit(self.construct_subelement(1))
+        elif dimension == 0:
+            return self.construct_subelement(0)
+        else:
+            raise ValueError("Illegal dimension")
+
+
 class MacroQuadratureRule(QuadratureRule):
     """Composite quadrature rule on parent facets that respects the splitting.
 

FIAT/powell_sabin.py

@@ -0,0 +1,108 @@
+# Copyright (C) 2024 Robert C. Kirby
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Robert C. Kirby (robert.c.kirby@gmail.com), 2024
+
+from FIAT import dual_set, finite_element, macro, polynomial_set
+from FIAT.functional import (IntegralMomentOfNormalDerivative, PointDerivative,
+                             PointEvaluation)
+from FIAT.jacobi import eval_jacobi_batch
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.reference_element import TRIANGLE, ufc_simplex
+
+
+class QuadraticPowellSabin6DualSet(dual_set.DualSet):
+    def __init__(self, ref_complex, degree=2):
+        if degree != 2:
+            raise ValueError("PS6 only defined for degree = 2")
+        ref_el = ref_complex.get_parent()
+        if ref_el.get_shape() != TRIANGLE:
+            raise ValueError("PS6 only defined on triangles")
+        top = ref_el.get_topology()
+        verts = ref_el.get_vertices()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+        # get first order jet at each vertex
+        alphas = polynomial_set.mis(sd, 1)
+        nodes = []
+
+        for v in sorted(top[0]):
+            pt = verts[v]
+            cur = len(nodes)
+            nodes.append(PointEvaluation(ref_el, pt))
+            nodes.extend(PointDerivative(ref_el, pt, alpha) for alpha in alphas)
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        super(QuadraticPowellSabin6DualSet, self).__init__(
+            nodes, ref_el, entity_ids)
+
+
+class QuadraticPowellSabin6(finite_element.CiarletElement):
+    """The PS6 macroelement is a C^1 quadratic macroelement defined
+    on the 6-way Powell-Sabin split of a triangle.
+    """
+    def __init__(self, ref_el, degree=2):
+        if degree != 2:
+            raise ValueError("PS6 only defined for degree = 2")
+        ref_complex = macro.PowellSabinSplit(ref_el)
+        dual = QuadraticPowellSabin6DualSet(ref_complex, degree)
+        poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1)
+
+        super(QuadraticPowellSabin6, self).__init__(poly_set, dual, degree)
+
+
+class QuadraticPowellSabin12DualSet(dual_set.DualSet):
+    def __init__(self, ref_complex, degree=2):
+        if degree != 2:
+            raise ValueError("PS12 only defined for degree = 2")
+        ref_el = ref_complex.get_parent()
+        if ref_el.get_shape() != TRIANGLE:
+            raise ValueError("PS12 only defined on triangles")
+        top = ref_el.get_topology()
+        verts = ref_el.get_vertices()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+        # get first order jet at each vertex
+        alphas = polynomial_set.mis(sd, 1)
+        nodes = []
+
+        for v in sorted(top[0]):
+            pt = verts[v]
+            cur = len(nodes)
+            nodes.append(PointEvaluation(ref_el, pt))
+            nodes.extend(PointDerivative(ref_el, pt, alpha) for alpha in alphas)
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # integral moment of normal derivatives
+        rline = macro.AlfeldSplit(ufc_simplex(1))
+        Q = create_quadrature(rline, degree-1)
+        qpts = Q.get_points()
+
+        x = 2.0*qpts - 1
+        phis = eval_jacobi_batch(1, 1, 0, x)
+        for e in sorted(top[1]):
+            cur = len(nodes)
+            nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis)
+            entity_ids[1][e].extend(range(cur, len(nodes)))
+
+        super(QuadraticPowellSabin12DualSet, self).__init__(
+            nodes, ref_el, entity_ids)
+
+
+class QuadraticPowellSabin12(finite_element.CiarletElement):
+    """The PS12 macroelement is a C^1 quadratic macroelement defined
+    on the 12-way Powell-Sabin split of a triangle.
+    """
+    def __init__(self, ref_el, degree=2):
+        if degree != 2:
+            raise ValueError("PS12 only defined for degree = 2")
+        ref_complex = macro.PowellSabin12Split(ref_el)
+        dual = QuadraticPowellSabin12DualSet(ref_complex, degree)
+        poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1)
+
+        super(QuadraticPowellSabin12, self).__init__(poly_set, dual, degree)

test/unit/test_powell_sabin.py

@@ -0,0 +1,32 @@
+import numpy
+import pytest
+from FIAT import QuadraticPowellSabin6 as PS6
+from FIAT import QuadraticPowellSabin12 as PS12
+from FIAT.functional import PointEvaluation
+from FIAT.reference_element import make_lattice, ufc_simplex
+
+
+@pytest.fixture
+def cell():
+    return ufc_simplex(2)
+
+
+@pytest.mark.parametrize("el", (PS6, PS12))
+def test_powell_sabin_constant(cell, el):
+    # Test that bfs associated with point evaluation sum up to 1
+    fe = el(cell)
+
+    pts = make_lattice(cell.get_vertices(), 3)
+    tab = fe.tabulate(2, pts)
+
+    coeffs = numpy.zeros((fe.space_dimension(),))
+    nodes = fe.dual_basis()
+    entity_dofs = fe.entity_dofs()
+    for v in entity_dofs[0]:
+        for k in entity_dofs[0][v]:
+            if isinstance(nodes[k], PointEvaluation):
+                coeffs[k] = 1.0
+
+    for alpha in tab:
+        expected = 1 if sum(alpha) == 0 else 0
+        assert numpy.allclose(coeffs @ tab[alpha], expected)