Functional factory

FIAT/brezzi_douglas_marini.py

@@ -9,38 +9,27 @@
                   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):
     def __init__(self, ref_el, degree, variant, interpolant_deg):
-        nodes = []
         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] = []
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
         if variant == "integral":
             facet = ref_el.get_facet_element()
-            # Facet nodes are \int_F v\cdot n p ds where p \in P_{q}
+            # Facet nodes are \int_F v.n p ds where p \in P_{q}
             # degree is q
-            Q_ref = create_quadrature(facet, interpolant_deg + degree)
+            Q = create_quadrature(facet, interpolant_deg + degree)
             Pq = polynomial_set.ONPolynomialSet(facet, degree)
-            Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
+            ells = functional.NormalMoments(ref_el, Q, Pq)
             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.extend(ells.generate(sd-1, f))
+                entity_ids[sd - 1][f].extend(range(cur, len(nodes)))
 
         elif variant == "point":
             # Define each functional for the dual set
@@ -50,7 +39,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 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)))
+                entity_ids[sd - 1][f].extend(range(cur, len(nodes)))
 
         # internal nodes
         if degree > 1:
@@ -60,10 +49,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             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)))
+            ells = functional.FrobeniusIntegralMoments(ref_el, Q, Nedfs)
+            nodes.extend(ells.generate(sd, 0))
+            entity_ids[sd][0].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 

FIAT/functional.py

@@ -16,7 +16,7 @@
 import numpy
 import sympy
 
-from FIAT import polynomial_set, jacobi
+from FIAT import polynomial_set, quadrature, jacobi
 from FIAT.quadrature import GaussLegendreQuadratureLineRule
 from FIAT.reference_element import UFCInterval as interval
 
@@ -754,3 +754,64 @@ def __init__(self, ref_el, Q, P_at_qpts, facet):
         pt_dict = {tuple(pt): [(wt*nnT[idx], idx) for idx in index_iterator(shp)]
                    for pt, wt in zip(points, weights)}
         super().__init__(ref_el, shp, pt_dict, {}, "IntegralMomentOfNormalNormalEvaluation")
+
+
+class FunctionalGenerator():
+    def __init__(self, ref_el):
+        self.ref_el = ref_el
+
+    def generate(self, dim, entity):
+        raise NotImplementedError
+
+
+class IntegralMoments(FunctionalGenerator):
+    def __init__(self, ref_el, Q, P, comp=tuple(), shp=None):
+        if shp is None:
+            shp = P.get_shape()
+        sd = P.ref_el.get_spatial_dimension()
+        self.phis = P.tabulate(Q.get_points())[(0,)*sd]
+        self.Q = Q
+        self.comp = comp
+        self.shape = shp
+        super().__init__(ref_el)
+
+    def generate(self, dim, entity):
+        if dim == self.ref_el.get_spatial_dimension():
+            assert entity == 0
+            for phi in self.phis:
+                yield IntegralMoment(self.ref_el, self.Q, phi, comp=self.comp, shp=self.shape)
+
+
+class FrobeniusIntegralMoments(IntegralMoments):
+    def __init__(self, ref_el, Q, P):
+        super().__init__(ref_el, Q, P, comp=slice(None))
+
+    def generate(self, dim, entity):
+        if dim == self.ref_el.get_spatial_dimension():
+            assert entity == 0
+            for phi in self.phis:
+                yield FrobeniusIntegralMoment(self.ref_el, self.Q, phi)
+
+
+class FacetIntegralMoments(FrobeniusIntegralMoments):
+    def get_entity_transform(self, dim, entity):
+        return lambda f: f
+
+    def generate(self, dim, entity):
+        transform = self.get_entity_transform(dim, entity)
+        Q_mapped = quadrature.FacetQuadratureRule(self.ref_el, dim, entity, self.Q)
+        Jdet = Q_mapped.jacobian_determinant()
+        for phi in self.phis:
+            yield FrobeniusIntegralMoment(self.ref_el, Q_mapped, transform(phi)/Jdet)
+
+
+class NormalMoments(FacetIntegralMoments):
+    def get_entity_transform(self, dim, entity):
+        n = self.ref_el.compute_scaled_normal(entity)
+        return lambda f: numpy.multiply(n[..., None], f)
+
+
+class TangentialMoments(FacetIntegralMoments):
+    def get_entity_transform(self, dim, entity):
+        ts = self.ref_el.compute_tangents(dim, entity)
+        return lambda f: numpy.dot(ts.T, f)

FIAT/nedelec.py

@@ -11,7 +11,6 @@
 import numpy
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
-from FIAT.quadrature import FacetQuadratureRule
 
 
 def NedelecSpace2D(ref_el, degree):
@@ -122,21 +121,13 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 phi_deg = degree - dim
                 if phi_deg >= 0:
                     facet = ref_el.construct_subelement(dim)
-                    Q_ref = create_quadrature(facet, interpolant_deg + phi_deg)
+                    Q = create_quadrature(facet, interpolant_deg + phi_deg)
                     Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,))
-                    Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim]
-                    Phis = numpy.transpose(Phis, (0, 2, 1))
-
+                    ells = functional.TangentialMoments(ref_el, Q, Pqmd)
                     for entity in top[dim]:
                         cur = len(nodes)
-                        Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                        Jdet = Q.jacobian_determinant()
-                        R = numpy.array(ref_el.compute_tangents(dim, entity))
-                        phis = numpy.dot(Phis, R / Jdet)
-                        phis = numpy.transpose(phis, (0, 2, 1))
-                        nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
-                                     for phi in phis)
-                        entity_ids[dim][entity] = list(range(cur, len(nodes)))
+                        nodes.extend(ells.generate(dim, entity))
+                        entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
         elif variant == "point":
             for i in top[1]:
@@ -145,7 +136,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 pts_cur = ref_el.make_points(1, i, degree + 1)
                 nodes.extend(functional.PointEdgeTangentEvaluation(ref_el, i, pt)
                              for pt in pts_cur)
-                entity_ids[1][i] = list(range(cur, len(nodes)))
+                entity_ids[1][i].extend(range(cur, len(nodes)))
 
             if sd > 2 and degree > 1:  # face tangents
                 for i in top[2]:  # loop over faces
@@ -155,21 +146,19 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                                  for k in range(2)  # loop over tangents
                                  for pt in pts_cur  # loop over points
                                  )
-                    entity_ids[2][i] = list(range(cur, len(nodes)))
+                    entity_ids[2][i].extend(range(cur, len(nodes)))
 
         # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-d}^3(T)
-        dim = sd
-        phi_deg = degree - dim
+        phi_deg = degree - sd
         if phi_deg >= 0:
             if interpolant_deg is None:
                 interpolant_deg = degree
             cur = len(nodes)
             Q = create_quadrature(ref_el, interpolant_deg + phi_deg)
-            Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg)
-            Phis = Pqmd.tabulate(Q.get_points())[(0,) * dim]
-            nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (dim,))
-                         for d in range(dim) for phi in Phis)
-            entity_ids[dim][0] = list(range(cur, len(nodes)))
+            Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, shape=(sd,))
+            ells = functional.FrobeniusIntegralMoments(ref_el, Q, Pqmd)
+            nodes.extend(ells.generate(sd, 0))
+            entity_ids[sd][0].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 

FIAT/nedelec_second_kind.py

@@ -5,15 +5,12 @@
 #
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
-import numpy
-
 from FIAT.finite_element import CiarletElement
 from FIAT.dual_set import DualSet
 from FIAT.polynomial_set import ONPolynomialSet
 from FIAT.functional import PointEdgeTangentEvaluation as Tangent
-from FIAT.functional import FrobeniusIntegralMoment as IntegralMoment
+from FIAT.functional import TangentialMoments
 from FIAT.raviart_thomas import RaviartThomas
-from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.check_format_variant import check_format_variant
 
@@ -109,61 +106,42 @@ def _generate_edge_dofs(self, cell, degree, offset, variant, interpolant_deg):
 
         return (dofs, ids)
 
-    def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant_deg):
+    def _generate_facet_dofs(self, dim, cell, degree, offset, variant, interpolant_deg):
         """Generate degrees of freedom (dofs) for facets."""
 
         # Initialize empty dofs and identifiers (ids)
-        num_facets = len(cell.get_topology()[codim])
+        top = cell.get_topology()
         dofs = []
-        ids = {i: [] for i in range(num_facets)}
+        ids = {i: [] for i in top[dim]}
 
         # Return empty info if not applicable
-        rt_degree = degree - codim + 1
+        rt_degree = degree - dim + 1
         if rt_degree < 1:
             return (dofs, ids)
 
         if interpolant_deg is None:
             interpolant_deg = degree
 
         # Construct quadrature scheme for the reference facet
-        ref_facet = cell.construct_subelement(codim)
-        Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree)
-        if codim == 1:
-            Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,))
+        facet = cell.construct_subelement(dim)
+        Q = create_quadrature(facet, interpolant_deg + rt_degree)
+        if dim == 1:
+            P = ONPolynomialSet(facet, rt_degree, (dim,))
         else:
             # Construct Raviart-Thomas on the reference facet
-            RT = RaviartThomas(ref_facet, rt_degree, variant)
-            Phi = RT.get_nodal_basis()
-
-        # Evaluate basis functions at reference quadrature points
-        Phis = Phi.tabulate(Q_ref.get_points())[(0,) * codim]
-        # Note: Phis has dimensions:
-        # num_basis_functions x num_components x num_quad_points
-        Phis = numpy.transpose(Phis, (0, 2, 1))
-        # Note: Phis has dimensions:
-        # num_basis_functions x num_quad_points x num_components
-
-        # Iterate over the facets
-        cur = offset
-        for facet in range(num_facets):
-            # Get the quadrature and Jacobian on this facet
-            Q_facet = FacetQuadratureRule(cell, codim, facet, Q_ref)
-            J = Q_facet.jacobian()
-            detJ = Q_facet.jacobian_determinant()
-
-            # Map Phis -> phis (reference values to physical values)
-            piola_map = J / detJ
-            phis = numpy.dot(Phis, piola_map.T)
-            phis = numpy.transpose(phis, (0, 2, 1))
+            RT = RaviartThomas(facet, rt_degree, variant)
+            P = RT.get_nodal_basis()
 
+        ells = TangentialMoments(cell, Q, P)
+        for entity in sorted(top[dim]):
+            cur = len(dofs)
             # Construct degrees of freedom as integral moments on this cell,
             # using the face quadrature weighted against the values
             # of the (physical) Raviart--Thomas'es on the face
-            dofs.extend(IntegralMoment(cell, Q_facet, phi) for phi in phis)
+            dofs.extend(ells.generate(dim, entity))
 
             # Assign identifiers (num RTs per face + previous edge dofs)
-            ids[facet].extend(range(cur, cur + len(phis)))
-            cur += len(phis)
+            ids[entity].extend(range(cur, len(dofs)))
 
         return (dofs, ids)
 

FIAT/raviart_thomas.py

@@ -11,7 +11,6 @@
 from itertools import chain
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
-from FIAT.quadrature import FacetQuadratureRule
 
 
 def RTSpace(ref_el, degree):
@@ -59,44 +58,32 @@ class RTDualSet(dual_set.DualSet):
     moments against polynomials"""
 
     def __init__(self, ref_el, degree, variant, interpolant_deg):
-        nodes = []
         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] = []
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
         if variant == "integral":
             facet = ref_el.get_facet_element()
             # Facet nodes are \int_F v\cdot n p ds where p \in P_q
             q = degree - 1
-            Q_ref = create_quadrature(facet, interpolant_deg + q)
+            Q = create_quadrature(facet, interpolant_deg + q)
             Pq = polynomial_set.ONPolynomialSet(facet, q if sd > 1 else 0)
-            Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
+            ells = functional.NormalMoments(ref_el, Q, Pq)
             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.extend(ells.generate(sd-1, f))
+                entity_ids[sd - 1][f].extend(range(cur, len(nodes)))
 
             # internal nodes. These are \int_T v \cdot p dx where p \in P_{q-1}^d
             if q > 0:
                 cur = len(nodes)
                 Q = create_quadrature(ref_el, interpolant_deg + q - 1)
-                Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1)
-                Pqm1_at_qpts = Pqm1.tabulate(Q.get_points())[(0,) * sd]
-                nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
-                             for d in range(sd)
-                             for phi in Pqm1_at_qpts)
-                entity_ids[sd][0] = list(range(cur, len(nodes)))
+                Pqm1 = polynomial_set.ONPolynomialSet(ref_el, q - 1, shape=(sd,))
+                ells = functional.FrobeniusIntegralMoments(ref_el, Q, Pqm1)
+                nodes.extend(ells.generate(sd, 0))
+                entity_ids[sd][0].extend(range(cur, len(nodes)))
 
         elif variant == "point":
             # codimension 1 facets
@@ -105,7 +92,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1)
                 nodes.extend(functional.PointScaledNormalEvaluation(ref_el, i, pt)
                              for pt in pts_cur)
-                entity_ids[sd - 1][i] = list(range(cur, len(nodes)))
+                entity_ids[sd - 1][i].extend(range(cur, len(nodes)))
 
             # internal nodes.  Let's just use points at a lattice
             if degree > 1:
@@ -114,7 +101,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 nodes.extend(functional.ComponentPointEvaluation(ref_el, d, (sd,), pt)
                              for d in range(sd)
                              for pt in pts)
-                entity_ids[sd][0] = list(range(cur, len(nodes)))
+                entity_ids[sd][0].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)