FunctionalBlock

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

@@ -9,7 +9,6 @@
                   polynomial_set, nedelec)
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
-from FIAT.quadrature import FacetQuadratureRule
 
 
 class BDMDualSet(dual_set.DualSet):
@@ -18,54 +17,34 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
 
-        entity_ids = {}
-        # set to empty
-        for dim in top:
-            entity_ids[dim] = {}
-            for entity in top[dim]:
-                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)
             Pq = polynomial_set.ONPolynomialSet(facet, degree)
-            Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
-                cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref)
-                Jdet = Q.jacobian_determinant()
-                n = ref_el.compute_scaled_normal(f) / Jdet
-                phis = n[None, :, None] * Pq_at_qpts[:, None, :]
-                nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
-                             for phi in phis)
-                entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
+                nodes.append(functional.FacetNormalIntegralMomentBlock(ref_el, sd-1, f, Q_ref, Pq))
 
         elif variant == "point":
-            # Define each functional for the dual set
-            # codimension 1 facets
+            # Normal component evaluation at lattice points on facets
             for f in top[sd - 1]:
-                cur = len(nodes)
-                pts_cur = ref_el.make_points(sd - 1, f, sd + degree)
-                nodes.extend(functional.PointScaledNormalEvaluation(ref_el, f, pt)
-                             for pt in pts_cur)
-                entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
+                nodes.append(functional.PointDirectionalEvaluationBlock(ref_el, sd-1, f, direction="normal", degree=degree+sd))
 
         # internal nodes
         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)))
 
-        super().__init__(nodes, ref_el, entity_ids)
+            cell = ref_el.construct_subelement(sd)
+            Q_ref = create_quadrature(cell, interpolant_deg + degree - 1)
+            Nedel = nedelec.Nedelec(cell, degree - 1, variant)
+            mapping, = set(Nedel.mapping())
+            Phi = Nedel.get_nodal_basis()
+            for entity in top[sd]:
+                nodes.append(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Phi, mapping=mapping))
+
+        super().__init__(nodes, ref_el)
 
 
 class BrezziDouglasMarini(finite_element.CiarletElement):
@@ -90,7 +69,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/dual_set.py

@@ -12,7 +12,20 @@
 
 
 class DualSet(object):
-    def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None):
+    def __init__(self, nodes, ref_el, entity_ids=None, entity_permutations=None):
+        if entity_ids is None:
+            # flatten nodes
+            top = ref_el.get_topology()
+            entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+            blocks = nodes
+            nodes = []
+            for ell in blocks:
+                cur = len(nodes)
+                nodes.extend(ell.nodes)
+                dim = ell.entity_dim
+                entity = ell.entity_id
+                entity_ids[dim][entity].extend(range(cur, len(nodes)))
+
         nodes, ref_el, entity_ids, entity_permutations = merge_entities(nodes, ref_el, entity_ids, entity_permutations)
         self.nodes = nodes
         self.ref_el = ref_el

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.")