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):
@@ -26,20 +25,15 @@ 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)
             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)
+                nodes.extend(functional.FacetNormalIntegralMomentBlock(ref_el, f, Q_ref, Pq_at_qpts))
                 entity_ids[sd - 1][f] = list(range(cur, len(nodes)))
 
         elif variant == "point":
@@ -56,14 +50,16 @@ 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)
+            mapping, = set(Nedel.mapping())
+            Ned_at_qpts = Nedel.tabulate(0, Q_ref.get_points())[(0,) * sd]
+            for entity in top[sd]:
+                cur = len(nodes)
+                nodes.extend(functional.FacetIntegralMomentBlock(ref_el, sd, entity, Q_ref, Ned_at_qpts, mapping=mapping))
+                entity_ids[sd][entity] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -90,7 +86,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/functional.py

@@ -15,7 +15,7 @@
 from itertools import chain
 import numpy
 
-from FIAT import polynomial_set, jacobi, quadrature_schemes
+from FIAT import polynomial_set, jacobi, quadrature, quadrature_schemes
 
 
 def index_iterator(shp):
@@ -25,7 +25,127 @@ def index_iterator(shp):
     return numpy.ndindex(shp)
 
 
-class Functional(object):
+def pullback(mapping, J, detJ, phi):
+    try:
+        formdegree = {
+            "identity": (0,),
+            "L2 piola": (3,),
+            "covariant piola": (1,),
+            "contravariant piola": (2,),
+            "double covariant piola": (1, 1),
+            "double contravariant piola": (2, 2),
+            "covariant contravariant piola": (1, 2),
+            "contravariant covariant piola": (2, 1), }[mapping]
+    except KeyError:
+        raise ValueError(f"Unrecognized mapping {mapping}")
+
+    F1 = numpy.linalg.pinv(J).T
+    F2 = J / detJ
+
+    perm = [*range(1, phi.ndim-1), 0, -1]
+    phi = phi.transpose(perm)
+
+    for i, k in enumerate(formdegree):
+        if k == 0:
+            continue
+        elif k == 3:
+            phi = phi / detJ
+        else:
+            F = F1 if k == 1 else F2
+            phi = numpy.tensordot(F, phi, (1, i))
+
+    perm = [-2, *reversed(range(phi.ndim-2)), -1]
+    phi = phi.transpose(perm)
+    return phi
+
+
+class FunctionalBlock:
+    def __init__(self, nodes):
+        self.nodes = nodes
+
+    def __iter__(self):
+        for ell in self.nodes:
+            yield ell
+
+    def __len__(self):
+        return len(self.nodes)
+
+    def completion(self):
+        return self
+
+
+class FunctionalBlockView(FunctionalBlock):
+    def __init__(self, block, index):
+        self.block = block
+        nodes = [block.nodes[index]]
+        super().__init__(nodes)
+
+    def completion(self):
+        return self.block
+
+
+class PointFunctionalBlock(FunctionalBlock):
+    def __init__(self, ref_el, x, order):
+        from FIAT.polynomial_set import mis
+        sd = ref_el.get_spatial_dimension()
+        nodes = [PointDerivative(ref_el, x, alpha) for alpha in mis(sd, order)]
+        super().__init__(nodes)
+
+
+class PointGradient(PointFunctionalBlock):
+    def __init__(self, ref_el, x):
+        super().__init__(ref_el, x, 1)
+
+
+class PointHessian(PointFunctionalBlock):
+    def __init__(self, ref_el, x):
+        super().__init__(ref_el, x, 2)
+
+
+class PointDirectionalDerivativeBlock(PointFunctionalBlock):
+    def __init__(self, ref_el, x, basis):
+        nodes = [PointDirectionalDerivative(ref_el, x, e) for e in basis]
+        super().__init__(nodes)
+
+
+class PointNormalTangentialDerivativeBlock(PointDirectionalDerivativeBlock):
+    def __init__(self, ref_el, x, facet_id):
+        sd = ref_el.get_spatial_dimension()
+        basis = numpy.zeros((sd, sd))
+        basis[0] = ref_el.compute_scaled_normal(facet_id)
+        basis[1:] = ref_el.compute_tangents(sd-1, facet_id)
+        super().__init__(self, x, basis)
+
+
+class PointNormalDerivativeView(FunctionalBlockView):
+    def __init__(self, ref_el, x, facet_id):
+        block = PointNormalTangentialDerivativeBlock(ref_el, x, facet_id)
+        super().__init__(block, 0)
+
+
+class FacetIntegralMomentBlock(FunctionalBlock):
+    def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, mapping="L2 piola"):
+        Q = quadrature.FacetQuadratureRule(ref_el, entity_dim, entity_id, Q_ref)
+        phis = pullback(mapping, Q.jacobian(), Q.jacobian_determinant(), Phis)
+        nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis]
+        super().__init__(nodes)
+
+
+class FacetDirectionalIntegralMomentBlock(FacetIntegralMomentBlock):
+    def __init__(self, ref_el, entity_dim, entity_id, Q_ref, Phis, direction):
+        direction = numpy.asarray(direction)
+        Phis = Phis[(slice(None), *(None for _ in range(direction.ndim)), slice(None))] * direction[(*(None for _ in range(Phis.ndim-1)), Ellipsis, None)]
+        super().__init__(ref_el, entity_dim, entity_id, Q_ref, Phis)
+
+
+class FacetNormalIntegralMomentBlock(FacetDirectionalIntegralMomentBlock):
+    def __init__(self, ref_el, facet_id, Q_ref, Phis):
+        direction = ref_el.compute_scaled_normal(facet_id)
+        sd = ref_el.get_spatial_dimension()
+        super().__init__(ref_el, sd-1, facet_id, Q_ref, Phis, direction)
+
+
+class Functional(FunctionalBlock):
     r"""Abstract class representing a linear functional.
     All FIAT functionals are discrete in the sense that
     they are written as a weighted sum of (derivatives of components of) their