merge conflict

FIAT/P0.py

@@ -40,7 +40,7 @@ def __init__(self, ref_el):
                 entity_ids[dim][entity] = [entity] if dim == sd else []
                 entity_permutations[dim][entity] = perms
 
-        super(P0Dual, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
+        super().__init__(nodes, ref_el, entity_ids, entity_permutations)
 
 
 class P0(finite_element.CiarletElement):
@@ -49,4 +49,4 @@ def __init__(self, ref_el):
         dual = P0Dual(ref_el)
         degree = 0
         formdegree = ref_el.get_spatial_dimension()  # n-form
-        super(P0, self).__init__(poly_set, dual, degree, formdegree)
+        super().__init__(poly_set, dual, degree, formdegree)

FIAT/Sminus.py

@@ -497,7 +497,7 @@ def __init__(self, ref_el, degree):
             self.basis = {(0, 0): Array(Sminus_list)}
         else:
             self.basis = {(0, 0, 0): Array(Sminus_list)}
-        super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola")
+        super().__init__(ref_el=ref_el, degree=degree, mapping="covariant piola")
 
 
 class TrimmedSerendipityFace(TrimmedSerendipity):
@@ -526,4 +526,4 @@ def __init__(self, ref_el, degree):
         Sminus_list = EL + FL
         Sminus_list = [[-a[1], a[0]] for a in Sminus_list]
         self.basis = {(0, 0): Array(Sminus_list)}
-        super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
+        super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")

FIAT/SminusCurl.py

@@ -88,11 +88,11 @@ def __init__(self, ref_el, degree, mapping):
 
         entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids)
 
-        super(TrimmedSerendipity, self).__init__(ref_el=ref_el,
-                                                 dual=None,
-                                                 order=degree,
-                                                 formdegree=formdegree,
-                                                 mapping=mapping)
+        super().__init__(ref_el=ref_el,
+                         dual=None,
+                         order=degree,
+                         formdegree=formdegree,
+                         mapping=mapping)
 
         topology = ref_el.get_topology()
         unflattening_map = compute_unflattening_map(topology)
@@ -231,7 +231,7 @@ def __init__(self, ref_el, degree):
 
             Sminus_list = EL + FL
             self.basis = {(0, 0): Array(Sminus_list)}
-            super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
+            super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
 
 
 def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid):

FIAT/SminusDiv.py

@@ -78,11 +78,11 @@ def __init__(self, ref_el, degree, mapping):
 
         entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids)
 
-        super(TrimmedSerendipity, self).__init__(ref_el=ref_el,
-                                                 dual=None,
-                                                 order=degree,
-                                                 formdegree=formdegree,
-                                                 mapping=mapping)
+        super().__init__(ref_el=ref_el,
+                         dual=None,
+                         order=degree,
+                         formdegree=formdegree,
+                         mapping=mapping)
 
         topology = ref_el.get_topology()
         unflattening_map = compute_unflattening_map(topology)
@@ -198,7 +198,7 @@ def __init__(self, ref_el, degree):
                 IL = ()
             Sminus_list = FL + IL
             self.basis = {(0, 0, 0): Array(Sminus_list)}
-            super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
+            super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
 
         else:
             # Put all 2 dimensional stuff here.
@@ -213,7 +213,7 @@ def __init__(self, ref_el, degree):
             Sminus_list = EL + FL
             Sminus_list = [[-a[1], a[0]] for a in Sminus_list]
             self.basis = {(0, 0): Array(Sminus_list)}
-            super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
+            super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")
 
 
 def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid):

FIAT/__init__.py

@@ -7,9 +7,14 @@
 # Import finite element classes
 from FIAT.finite_element import FiniteElement, CiarletElement  # noqa: F401
 from FIAT.argyris import Argyris
+from FIAT.bernardi_raugel import BernardiRaugel
 from FIAT.bernstein import Bernstein
 from FIAT.bell import Bell
 from FIAT.hct import HsiehCloughTocher
+from FIAT.alfeld_sorokina import AlfeldSorokina
+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.brezzi_douglas_marini import BrezziDouglasMarini
 from FIAT.Sminus import TrimmedSerendipityEdge  # noqa: F401
@@ -31,8 +36,7 @@
 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.powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12
 from FIAT.hierarchical import Legendre, IntegratedLegendre
 from FIAT.P0 import P0
 from FIAT.raviart_thomas import RaviartThomas
@@ -65,6 +69,7 @@
 # List of supported elements and mapping to element classes
 supported_elements = {"Argyris": Argyris,
                       "Bell": Bell,
+                      "Bernardi-Raugel": BernardiRaugel,
                       "Bernstein": Bernstein,
                       "Brezzi-Douglas-Marini": BrezziDouglasMarini,
                       "Brezzi-Douglas-Fortin-Marini": BrezziDouglasFortinMarini,
@@ -79,7 +84,13 @@
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,
-                      "HsiehCloughTocher": HsiehCloughTocher,
+                      "Hsieh-Clough-Tocher": HsiehCloughTocher,
+                      "Alfeld-Sorokina": AlfeldSorokina,
+                      "Arnold-Qin": ArnoldQin,
+                      "Christiansen-Hu": ChristiansenHu,
+                      "Guzman-Neilan 1st kind H1": GuzmanNeilanFirstKindH1,
+                      "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1,
+                      "Guzman-Neilan H1(div)": GuzmanNeilanH1div,
                       "Johnson-Mercier": JohnsonMercier,
                       "Lagrange": Lagrange,
                       "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen,

FIAT/alfeld_sorokina.py

@@ -0,0 +1,84 @@
+# Copyright (C) 2024 Pablo D. Brubeck
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Pablo D. Brubeck ([email protected]), 2024
+
+from FIAT import finite_element, dual_set, polynomial_set
+from FIAT.functional import ComponentPointEvaluation, PointDivergence
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.macro import CkPolynomialSet, AlfeldSplit
+
+import numpy
+
+
+def AlfeldSorokinaSpace(ref_el, degree):
+    """Return a vector-valued C0 PolynomialSet on an Alfeld split with C0
+    divergence. This works on any simplex and for all polynomial degrees."""
+    ref_complex = AlfeldSplit(ref_el)
+    sd = ref_complex.get_spatial_dimension()
+    C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), variant="bubble")
+    expansion_set = C0.get_expansion_set()
+    num_members = C0.get_num_members()
+    coeffs = C0.get_coeffs()
+
+    facet_el = ref_complex.construct_subelement(sd-1)
+    phi = polynomial_set.ONPolynomialSet(facet_el, 0 if sd == 1 else degree-1)
+    Q = create_quadrature(facet_el, 2 * phi.degree)
+    qpts, qwts = Q.get_points(), Q.get_weights()
+    phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)]
+    weights = numpy.multiply(phi_at_qpts, qwts)
+
+    rows = []
+    for facet in ref_complex.get_interior_facets(sd-1):
+        n = ref_complex.compute_normal(facet)
+        jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=1)
+        div_jump = n[:, None, None] * jumps[1][None, ...]
+        r = numpy.tensordot(div_jump, weights, axes=(-1, -1))
+        rows.append(r.reshape(num_members, -1).T)
+
+    if len(rows) > 0:
+        dual_mat = numpy.vstack(rows)
+        nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True)
+        coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0))
+    return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs)
+
+
+class AlfeldSorokinaDualSet(dual_set.DualSet):
+    def __init__(self, ref_el, degree):
+        if degree != 2:
+            raise NotImplementedError(f"{type(self).__name__} only defined for degree = 2")
+
+        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 = []
+        for dim in sorted(top):
+            for entity in sorted(top[dim]):
+                cur = len(nodes)
+
+                dpts = ref_el.make_points(dim, entity, degree-1)
+                nodes.extend(PointDivergence(ref_el, pt) for pt in dpts)
+
+                pts = ref_el.make_points(dim, entity, degree)
+                nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt)
+                             for pt in pts for k in range(sd))
+                entity_ids[dim][entity].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class AlfeldSorokina(finite_element.CiarletElement):
+    """The Alfeld-Sorokina C0 quadratic macroelement with C0 divergence.
+
+    This element belongs to a Stokes complex, and is paired with CG1(Alfeld).
+    """
+    def __init__(self, ref_el, degree=2):
+        dual = AlfeldSorokinaDualSet(ref_el, degree)
+        poly_set = AlfeldSorokinaSpace(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension() - 1  # (n-1)-form
+        super().__init__(poly_set, dual, degree, formdegree,
+                         mapping="contravariant piola")

FIAT/argyris.py

@@ -82,7 +82,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 entity_ids[2][0] = list(range(cur, len(nodes)))
         else:
             raise ValueError("Invalid variant for Argyris")
-        super(ArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids)
+        super().__init__(nodes, ref_el, entity_ids)
 
 
 class Argyris(finite_element.CiarletElement):
@@ -107,4 +107,4 @@ def __init__(self, ref_el, degree=5, variant=None):
 
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble")
         dual = ArgyrisDualSet(ref_el, degree, variant, interpolant_deg)
-        super(Argyris, self).__init__(poly_set, dual, degree)
+        super().__init__(poly_set, dual, degree)

FIAT/arnold_qin.py

@@ -0,0 +1,71 @@
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Pablo D. Brubeck ([email protected]), 2024
+
+from FIAT import finite_element, polynomial_set
+from FIAT.bernardi_raugel import BernardiRaugelDualSet
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.reference_element import TRIANGLE
+from FIAT.macro import CkPolynomialSet
+from FIAT.hct import HsiehCloughTocher
+
+import numpy
+
+
+def ArnoldQinSpace(ref_el, degree, reduced=False):
+    """Return a basis for the Arnold-Qin space
+    curl(HCT-red) + P0 x if reduced = True, and
+    curl(HCT) + P0 x if reduced = False."""
+    if ref_el.get_shape() != TRIANGLE:
+        raise ValueError("Arnold-Qin only defined on triangles")
+    if degree != 2:
+        raise ValueError("Arnold-Qin only defined for degree = 2")
+    sd = ref_el.get_spatial_dimension()
+    HCT = HsiehCloughTocher(ref_el, degree+1, reduced=True)
+    ref_complex = HCT.get_reference_complex()
+    Q = create_quadrature(ref_complex, 2 * degree)
+    qpts, qwts = Q.get_points(), Q.get_weights()
+
+    x = qpts.T
+    bary = numpy.asarray(ref_el.make_points(sd, 0, sd+1))
+    P0x_at_qpts = x[None, :, :] - bary[:, :, None]
+
+    tab = HCT.tabulate(1, qpts)
+    curl_at_qpts = numpy.stack([tab[(0, 1)], -tab[(1, 0)]], axis=1)
+    if reduced:
+        curl_at_qpts = curl_at_qpts[:9]
+
+    C0 = CkPolynomialSet(ref_complex, degree, order=0, scale=1, variant="bubble")
+    C0_at_qpts = C0.tabulate(qpts)[(0,) * sd]
+    duals = numpy.multiply(C0_at_qpts, qwts)
+    M = numpy.dot(duals, C0_at_qpts.T)
+    duals = numpy.linalg.solve(M, duals)
+
+    # Remove the constant nullspace
+    ids = [0, 3, 6]
+    A = numpy.asarray([[1, 1, 1], [1, -1, 0], [0, -1, 1]])
+    phis = curl_at_qpts
+    phis[ids] = numpy.tensordot(A, phis[ids], axes=(-1, 0))
+    # Replace the constant nullspace with P_0 x
+    phis[0] = P0x_at_qpts
+    coeffs = numpy.tensordot(phis, duals, axes=(-1, -1))
+    return polynomial_set.PolynomialSet(ref_complex, degree, degree,
+                                        C0.get_expansion_set(), coeffs)
+
+
+class ArnoldQin(finite_element.CiarletElement):
+    """The Arnold-Qin C0(Alfeld) quadratic macroelement with divergence in P0.
+    This element belongs to a Stokes complex, and is paired with unsplit DG0."""
+    def __init__(self, ref_el, degree=2, reduced=False):
+        poly_set = ArnoldQinSpace(ref_el, degree)
+        if reduced:
+            order = 1
+            mapping = "contravariant piola"
+        else:
+            order = degree
+            mapping = "affine"
+        dual = BernardiRaugelDualSet(ref_el, order, degree=degree)
+        formdegree = ref_el.get_spatial_dimension() - 1  # (n-1)-form
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)

FIAT/arnold_winther.py

@@ -25,7 +25,7 @@
 
 
 class ArnoldWintherNCDual(DualSet):
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree=2):
         if not degree == 2:
             raise ValueError("Nonconforming Arnold-Winther elements are"
                              "only defined for degree 2.")
@@ -70,24 +70,23 @@ def __init__(self, cell, degree):
             dof_ids[1][entity_id].append(dof_cur)
             dof_cur += 1
 
-        super(ArnoldWintherNCDual, self).__init__(dofs, cell, dof_ids)
+        super().__init__(dofs, cell, dof_ids)
 
 
 class ArnoldWintherNC(CiarletElement):
     """The definition of the nonconforming Arnold-Winther element.
     """
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree=2):
         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)
+        super().__init__(Ps, Ls, degree, mapping=mapping)
 
 
 class ArnoldWintherDual(DualSet):
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree=3):
         if not degree == 3:
             raise ValueError("Arnold-Winther elements are"
                              "only defined for degree 3.")
@@ -152,15 +151,15 @@ def __init__(self, cell, degree):
 
         dof_ids[2][0] += list(range(dof_cur, dof_cur+6))
 
-        super(ArnoldWintherDual, self).__init__(dofs, cell, dof_ids)
+        super().__init__(dofs, cell, dof_ids)
 
 
 class ArnoldWinther(CiarletElement):
     """The definition of the conforming Arnold-Winther element.
     """
-    def __init__(self, cell, degree):
+    def __init__(self, cell, degree=3):
         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)
+        super().__init__(Ps, Ls, degree, mapping=mapping)

FIAT/barycentric_interpolation.py

@@ -73,7 +73,7 @@ def __init__(self, ref_el, pts):
 
         self.degree = max(len(wts) for wts in self.weights)-1
         self.recurrence_order = self.degree + 1
-        super(LagrangeLineExpansionSet, self).__init__(ref_el)
+        super().__init__(ref_el)
         self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0"
 
     def get_num_members(self, n):
@@ -123,5 +123,4 @@ def __init__(self, ref_el, pts, shape=tuple()):
                     coeffs[cur_idx] = 1.0
                     cur_bf += 1
 
-        super(LagrangePolynomialSet, self).__init__(ref_el, degree, embedded_degree,
-                                                    expansion_set, coeffs)
+        super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)