Enable macroelement variants for most elements (#152)

* Macroize all the non-zany elements

* SimplicialComplex.get_facet_element

Author: pbrubeck

FIAT/argyris.py

@@ -53,13 +53,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             # interior dofs
             q = degree - 6
             if q >= 0:
-                Q = create_quadrature(ref_el, interpolant_deg + q)
-                Pq = polynomial_set.ONPolynomialSet(ref_el, q, scale=1)
-                phis = Pq.tabulate(Q.get_points())[(0,) * sd]
-                scale = ref_el.volume()
-                cur = len(nodes)
-                nodes.extend(IntegralMoment(ref_el, Q, phi/scale) for phi in phis)
-                entity_ids[sd][0] = list(range(cur, len(nodes)))
+                cell = ref_el.construct_subelement(sd)
+                Q_ref = create_quadrature(cell, interpolant_deg + q)
+                Pq = polynomial_set.ONPolynomialSet(cell, q, scale=1)
+                Phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd]
+                for entity in sorted(top[sd]):
+                    Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
+                    scale = 1 / Q.jacobian_determinant()
+                    phis = scale * Phis
+                    cur = len(nodes)
+                    nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis)
+                    entity_ids[sd][entity] = list(range(cur, len(nodes)))
 
         elif variant == "point":
             # edge dofs
@@ -77,9 +81,10 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             # interior dofs
             if degree > 5:
                 cur = len(nodes)
-                internalpts = ref_el.make_points(2, 0, degree - 3)
-                nodes.extend(PointEvaluation(ref_el, pt) for pt in internalpts)
-                entity_ids[2][0] = list(range(cur, len(nodes)))
+                for entity in sorted(top[sd]):
+                    internalpts = ref_el.make_points(sd, entity, degree - 3)
+                    nodes.extend(PointEvaluation(ref_el, pt) for pt in internalpts)
+                    entity_ids[sd][entity] = list(range(cur, len(nodes)))
         else:
             raise ValueError("Invalid variant for Argyris")
         super().__init__(nodes, ref_el, entity_ids)
@@ -103,7 +108,9 @@ class Argyris(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree=5, variant=None):
 
-        variant, interpolant_deg = check_format_variant(variant, degree)
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
+        if splitting is not None:
+            raise NotImplementedError(f"{type(self).__name__} is not implemented as a macroelement.")
 
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble")
         dual = ArgyrisDualSet(ref_el, degree, variant, interpolant_deg)

FIAT/brezzi_douglas_marini.py

@@ -10,6 +10,7 @@
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
+import numpy
 
 
 class BDMDualSet(dual_set.DualSet):
@@ -26,7 +27,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 entity_ids[dim][entity] = []
 
         if variant == "integral":
-            facet = ref_el.get_facet_element()
+            facet = ref_el.construct_subelement(sd-1)
             # Facet nodes are \int_F v\cdot n p ds where p \in P_{q}
             # degree is q
             Q_ref = create_quadrature(facet, interpolant_deg + degree)
@@ -56,14 +57,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
         if degree > 1:
             if interpolant_deg is None:
                 interpolant_deg = degree
-            cur = len(nodes)
-            Q = create_quadrature(ref_el, interpolant_deg + degree - 1)
-            Nedel = nedelec.Nedelec(ref_el, degree - 1, variant)
-            Nedfs = Nedel.get_nodal_basis()
-            Ned_at_qpts = Nedfs.tabulate(Q.get_points())[(0,) * sd]
-            nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
-                         for phi in Ned_at_qpts)
-            entity_ids[sd][0] = list(range(cur, len(nodes)))
+
+            cell = ref_el.construct_subelement(sd)
+            Q_ref = create_quadrature(cell, interpolant_deg + degree - 1)
+            Nedel = nedelec.Nedelec(cell, degree - 1, variant)
+            Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd]
+            for entity in top[sd]:
+                Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
+                Jinv = numpy.linalg.inv(Q.jacobian())
+                phis = numpy.tensordot(Jinv.T, Ned_at_qpts, (1, 1)).transpose((1, 0, 2))
+                cur = len(nodes)
+                nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis)
+                entity_ids[sd][entity] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -90,7 +95,9 @@ class BrezziDouglasMarini(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree, variant=None):
 
-        variant, interpolant_deg = check_format_variant(variant, degree)
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
+        if splitting is not None:
+            ref_el = splitting(ref_el)
 
         if degree < 1:
             raise Exception("BDM_k elements only valid for k >= 1")

FIAT/check_format_variant.py

@@ -27,6 +27,8 @@
 
 
 def check_format_variant(variant, degree):
+    splitting, variant = parse_lagrange_variant(variant, integral=True)
+
     if variant is None:
         variant = "integral"
 
@@ -44,7 +46,7 @@ def check_format_variant(variant, degree):
         raise ValueError('Choose either variant="point" or variant="integral"'
                          'or variant="integral(q)"')
 
-    return variant, interpolant_degree
+    return splitting, variant, interpolant_degree
 
 
 def parse_lagrange_variant(variant, discontinuous=False, integral=False):
@@ -61,7 +63,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
 
     default = "integral" if integral else "spectral"
     if integral:
-        supported_point_variants = {"integral": None}
+        supported_point_variants = {"integral": None, "point": "point"}
     elif discontinuous:
         supported_point_variants = supported_dg_variants
     else:
@@ -81,6 +83,8 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
             k, = match.groups()
             call_split = IsoSplit
             splitting_args = (int(k),)
+        elif opt.startswith("integral"):
+            point_variant = opt
         elif opt in supported_point_variants:
             point_variant = supported_point_variants[opt]
         else:

FIAT/crouzeix_raviart.py

@@ -81,12 +81,13 @@ class CrouzeixRaviart(finite_element.CiarletElement):
     """
 
     def __init__(self, ref_el, degree, variant=None):
-
-        variant, interpolant_deg = check_format_variant(variant, degree)
-
         if degree % 2 != 1:
             raise ValueError("Crouzeix-Raviart only defined for odd degree")
 
-        space = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble")
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
+        if splitting is not None:
+            ref_el = splitting(ref_el)
+
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
         dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg)
-        super().__init__(space, dual, degree)
+        super().__init__(poly_set, dual, degree)

FIAT/expansions.py

@@ -274,11 +274,18 @@ def __init__(self, ref_el, scale=None, variant=None):
         if scale is None:
             scale = math.sqrt(1.0 / base_ref_el.volume())
         self.scale = scale
+        self.variant = variant
         self.continuity = "C0" if variant == "bubble" else None
         self.recurrence_order = 2
         self._dmats_cache = {}
         self._cell_node_map_cache = {}
 
+    def reconstruct(self, ref_el=None, scale=None, variant=None):
+        """Reconstructs this ExpansionSet with modified arguments."""
+        return ExpansionSet(ref_el or self.ref_el,
+                            scale=scale or self.scale,
+                            variant=variant or self.variant)
+
     def get_scale(self, n, cell=0):
         scale = self.scale
         sd = self.ref_el.get_spatial_dimension()
@@ -371,7 +378,7 @@ def _tabulate(self, n, pts, order=0):
                     phi[alpha] /= mult[None, ipts]
 
         # Insert subcell tabulations into the corresponding submatrices
-        idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args)
+        idx = lambda *args: args if args[-1] is Ellipsis else numpy.ix_(*args)
         num_phis = self.get_num_members(n)
         cell_node_map = self.get_cell_node_map(n)
         result = {}
@@ -599,7 +606,10 @@ def polynomial_dimension(ref_el, n, continuity=None):
             raise ValueError("Only degree zero polynomials supported on point elements.")
         return 1
     top = ref_el.get_topology()
-    if continuity == "C0":
+
+    if isinstance(continuity, dict):
+        space_dimension = sum(len(continuity[dim][0]) * len(top[dim]) for dim in top)
+    elif continuity == "C0":
         space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top)
     else:
         dim = ref_el.get_spatial_dimension()
@@ -620,7 +630,9 @@ def polynomial_entity_ids(ref_el, n, continuity=None):
     entity_ids = {}
     cur = 0
     for dim in sorted(top):
-        if continuity == "C0":
+        if isinstance(continuity, dict):
+            dofs, = set(len(continuity[dim][entity]) for entity in continuity[dim])
+        elif continuity == "C0":
             dofs = math.comb(n - 1, dim)
         else:
             # DG numbering

FIAT/finite_element.py

@@ -15,6 +15,7 @@
 from FIAT.dual_set import DualSet
 from FIAT.polynomial_set import PolynomialSet
 from FIAT.quadrature_schemes import create_quadrature
+from FIAT.macro import MacroPolynomialSet
 
 
 class FiniteElement(object):
@@ -134,6 +135,11 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref
         ref_complex = ref_complex or poly_set.get_reference_element()
         super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex)
 
+        # Tile the poly_set
+        if ref_complex.is_macrocell() and len(poly_set) > len(dual):
+            base_element = type(self)(ref_el, order)
+            poly_set = MacroPolynomialSet(ref_complex, base_element)
+
         if len(poly_set) != len(dual):
             raise ValueError(f"Dimension of function space is {len(poly_set)}, but got {len(dual)} nodes.")
 

FIAT/gopalakrishnan_lederer_schoberl.py

@@ -1,4 +1,5 @@
 from FIAT import finite_element, dual_set, polynomial_set, expansions
+from FIAT.check_format_variant import check_format_variant
 from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
@@ -12,34 +13,29 @@ def __init__(self, ref_el, degree):
         nodes = []
         entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
 
-        # Face dofs: moments of normal-tangential components against a basis for Pk
-        ref_facet = ref_el.construct_subelement(sd-1)
-        Qref = create_quadrature(ref_facet, 2*degree)
-        P = polynomial_set.ONPolynomialSet(ref_facet, degree)
-        phis = P.tabulate(Qref.get_points())[(0,) * (sd-1)]
-        for f in sorted(top[sd-1]):
-            cur = len(nodes)
-            Q = FacetQuadratureRule(ref_el, sd-1, f, Qref)
-            Jdet = Q.jacobian_determinant()
-            normal = ref_el.compute_scaled_normal(f)
-            tangents = ref_el.compute_tangents(sd-1, f)
-            n = normal / Jdet
-            nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
-                         for phi in phis for t in tangents)
-            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
-
+        # Face dofs: moments of normal-tangential components against a basis for P_k
         # Interior dofs: moments of normal-tangential components against a basis for P_{k-1}
-        if degree > 0:
-            cur = len(nodes)
-            Q = create_quadrature(ref_el, 2*degree-1)
-            P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
-            phis = P.tabulate(Q.get_points())[(0,) * sd]
-            for f in sorted(top[sd-1]):
-                n = ref_el.compute_scaled_normal(f)
-                tangents = ref_el.compute_tangents(sd-1, f)
-                nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
-                             for phi in phis for t in tangents)
-            entity_ids[sd][0].extend(range(cur, len(nodes)))
+        for dim in (sd-1, sd):
+            q = degree + sd-1 - dim
+            if q < 0:
+                continue
+
+            ref_facet = ref_el.construct_subelement(dim)
+            Q_ref = create_quadrature(ref_facet, degree + q)
+            P = polynomial_set.ONPolynomialSet(ref_facet, q, scale=1)
+            phis = P.tabulate(Q_ref.get_points())[(0,) * dim]
+
+            for entity in sorted(top[dim]):
+                cur = len(nodes)
+                Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
+                Jdet = Q.jacobian_determinant()
+                for f in ref_el.get_connectivity()[(dim, sd-1)][entity]:
+                    normal = ref_el.compute_scaled_normal(f)
+                    tangents = ref_el.compute_tangents(sd-1, f)
+                    n = normal / Jdet
+                    nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                                 for phi in phis for t in tangents)
+                entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -59,7 +55,14 @@ class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement):
     the weak symmetry constraint.
 
     """
-    def __init__(self, ref_el, degree):
+    def __init__(self, ref_el, degree, variant=None):
+
+        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
+        assert variant == "integral"
+
+        if splitting is not None:
+            ref_el = splitting(ref_el)
+
         poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree)
         dual = GLSDual(ref_el, degree)
         sd = ref_el.get_spatial_dimension()
@@ -68,7 +71,7 @@ def __init__(self, ref_el, degree):
         super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
 
 
-def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree):
+def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree, variant=None):
     """The GLS element used for the Mass-Conserving mixed Stress (MCS)
     formulation for Stokes flow.
 
@@ -77,12 +80,16 @@ def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree):
 
     Reference: https://doi.org/10.1093/imanum/drz022
     """
-    fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree)
+    fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree, variant=variant)
     entity_dofs = fe.entity_dofs()
     sd = ref_el.get_spatial_dimension()
-    dimPkm1 = (sd-1)*expansions.polynomial_dimension(ref_el.construct_subelement(sd-1), degree-1)
+    facet = ref_el.construct_subelement(sd-1)
+    dimPkm1 = (sd-1)*expansions.polynomial_dimension(facet, degree-1)
+
     indices = []
-    for f in entity_dofs[sd-1]:
+    for f in sorted(entity_dofs[sd-1]):
         indices.extend(entity_dofs[sd-1][f][:dimPkm1])
-    indices.extend(entity_dofs[sd][0])
+    for cell in sorted(entity_dofs[sd]):
+        indices.extend(entity_dofs[sd][cell])
+
     return RestrictedElement(fe, indices=indices)