Rescale integral moment dofs to match lowest order point and integral variants

Author: pbrubeck

FIAT/brezzi_douglas_marini.py

@@ -30,7 +30,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             # 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 = polynomial_set.ONPolynomialSet(facet, degree, scale="L2 piola")
             Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
                 cur = len(nodes)

FIAT/expansions.py

@@ -141,31 +141,36 @@ def xi_tetrahedron(eta):
 
 
 class ExpansionSet(object):
-    def __new__(cls, ref_el, *args, **kwargs):
+    def __new__(cls, *args, **kwargs):
         """Returns an ExpansionSet instance appopriate for the given
         reference element."""
         if cls is not ExpansionSet:
             return super(ExpansionSet, cls).__new__(cls)
+        ref_el = args[0]
         if ref_el.get_shape() == reference_element.POINT:
-            return PointExpansionSet(ref_el)
+            return PointExpansionSet(*args, **kwargs)
         elif ref_el.get_shape() == reference_element.LINE:
-            return LineExpansionSet(ref_el)
+            return LineExpansionSet(*args, **kwargs)
         elif ref_el.get_shape() == reference_element.TRIANGLE:
-            return TriangleExpansionSet(ref_el)
+            return TriangleExpansionSet(*args, **kwargs)
         elif ref_el.get_shape() == reference_element.TETRAHEDRON:
-            return TetrahedronExpansionSet(ref_el)
+            return TetrahedronExpansionSet(*args, **kwargs)
         else:
             raise ValueError("Invalid reference element type.")
 
-    def __init__(self, ref_el):
+    def __init__(self, ref_el, scale=None):
+        if isinstance(scale, str) and scale.lower() == "l2 piola":
+            scale = 1.0 / ref_el.volume()
+        elif scale is None:
+            scale = math.sqrt(1.0 / ref_el.volume())
         self.ref_el = ref_el
+        self.scale = scale
         dim = ref_el.get_spatial_dimension()
         self.base_ref_el = reference_element.default_simplex(dim)
         v1 = ref_el.get_vertices()
         v2 = self.base_ref_el.get_vertices()
         self.A, self.b = reference_element.make_affine_mapping(v1, v2)
         self.mapping = lambda x: numpy.dot(self.A, x) + self.b
-        self.scale = numpy.sqrt(1.0 / ref_el.volume())
         self._dmats_cache = {}
 
     def get_num_members(self, n):
@@ -266,10 +271,10 @@ def tabulate_jet(self, n, pts, order=1):
 
 class PointExpansionSet(ExpansionSet):
     """Evaluates the point basis on a point reference element."""
-    def __init__(self, ref_el):
+    def __init__(self, ref_el, **kwargs):
         if ref_el.get_spatial_dimension() != 0:
             raise ValueError("Must have a point")
-        super(PointExpansionSet, self).__init__(ref_el)
+        super(PointExpansionSet, self).__init__(ref_el, **kwargs)
 
     def tabulate(self, n, pts):
         """Returns a numpy array A[i,j] = phi_i(pts[j]) = 1.0."""
@@ -279,10 +284,10 @@ def tabulate(self, n, pts):
 
 class LineExpansionSet(ExpansionSet):
     """Evaluates the Legendre basis on a line reference element."""
-    def __init__(self, ref_el):
+    def __init__(self, ref_el, **kwargs):
         if ref_el.get_spatial_dimension() != 1:
             raise Exception("Must have a line")
-        super(LineExpansionSet, self).__init__(ref_el)
+        super(LineExpansionSet, self).__init__(ref_el, **kwargs)
 
     def _tabulate(self, n, pts, order=0):
         """Returns a tuple of (vals, derivs) such that
@@ -306,18 +311,18 @@ def _tabulate(self, n, pts, order=0):
 class TriangleExpansionSet(ExpansionSet):
     """Evaluates the orthonormal Dubiner basis on a triangular
     reference element."""
-    def __init__(self, ref_el):
+    def __init__(self, ref_el, **kwargs):
         if ref_el.get_spatial_dimension() != 2:
             raise Exception("Must have a triangle")
-        super(TriangleExpansionSet, self).__init__(ref_el)
+        super(TriangleExpansionSet, self).__init__(ref_el, **kwargs)
 
 
 class TetrahedronExpansionSet(ExpansionSet):
     """Collapsed orthonormal polynomial expansion on a tetrahedron."""
-    def __init__(self, ref_el):
+    def __init__(self, ref_el, **kwargs):
         if ref_el.get_spatial_dimension() != 3:
             raise Exception("Must be a tetrahedron")
-        super(TetrahedronExpansionSet, self).__init__(ref_el)
+        super(TetrahedronExpansionSet, self).__init__(ref_el, **kwargs)
 
 
 def polynomial_dimension(ref_el, degree):

FIAT/nedelec.py

@@ -123,7 +123,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 if phi_deg >= 0:
                     facet = ref_el.construct_subelement(dim)
                     Q_ref = create_quadrature(facet, interpolant_deg + phi_deg)
-                    Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,))
+                    Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,), scale="L2 piola")
                     Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim]
                     Phis = numpy.transpose(Phis, (0, 2, 1))
 
@@ -165,7 +165,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
                 interpolant_deg = degree
             cur = len(nodes)
             Q = create_quadrature(ref_el, interpolant_deg + phi_deg)
-            Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg)
+            Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, scale="L2 piola")
             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)

FIAT/nedelec_second_kind.py

@@ -129,7 +129,7 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant
         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,))
+            Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,), scale="L2 piola")
         else:
             # Construct Raviart-Thomas on the reference facet
             RT = RaviartThomas(ref_facet, rt_degree, variant)
@@ -149,10 +149,10 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant
             # 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()
+            Jdet = Q_facet.jacobian_determinant()
 
             # Map Phis -> phis (reference values to physical values)
-            piola_map = J / detJ
+            piola_map = J / Jdet
             phis = numpy.dot(Phis, piola_map.T)
             phis = numpy.transpose(phis, (0, 2, 1))
 

FIAT/polynomial_set.py

@@ -120,7 +120,7 @@ class ONPolynomialSet(PolynomialSet):
 
     """
 
-    def __init__(self, ref_el, degree, shape=tuple()):
+    def __init__(self, ref_el, degree, shape=tuple(), scale=None):
 
         if shape == tuple():
             num_components = 1
@@ -130,7 +130,7 @@ def __init__(self, ref_el, degree, shape=tuple()):
         num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
         num_members = num_components * num_exp_functions
         embedded_degree = degree
-        expansion_set = expansions.ExpansionSet(ref_el)
+        expansion_set = expansions.ExpansionSet(ref_el, scale=scale)
 
         # set up coefficients
         if shape == tuple():
@@ -211,7 +211,7 @@ class ONSymTensorPolynomialSet(PolynomialSet):
 
     """
 
-    def __init__(self, ref_el, degree, size=None):
+    def __init__(self, ref_el, degree, size=None, scale=None):
 
         sd = ref_el.get_spatial_dimension()
         if size is None:
@@ -222,7 +222,7 @@ def __init__(self, ref_el, degree, size=None):
         num_components = size * (size + 1) // 2
         num_members = num_components * num_exp_functions
         embedded_degree = degree
-        expansion_set = expansions.ExpansionSet(ref_el)
+        expansion_set = expansions.ExpansionSet(ref_el, scale=scale)
 
         # set up coefficients for symmetric tensors
         coeffs_shape = (num_members, *shape, num_exp_functions)

FIAT/raviart_thomas.py

@@ -75,7 +75,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
             # degree is q - 1
             Q_ref = create_quadrature(facet, interpolant_deg + degree - 1)
-            Pq = polynomial_set.ONPolynomialSet(facet, degree - 1)
+            Pq = polynomial_set.ONPolynomialSet(facet, degree - 1, scale="L2 piola")
             Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
                 cur = len(nodes)
@@ -91,7 +91,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg):
             if degree > 1:
                 cur = len(nodes)
                 Q = create_quadrature(ref_el, interpolant_deg + degree - 2)
-                Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2)
+                Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2, scale="L2 piola")
                 Pkm1_at_qpts = Pkm1.tabulate(Q.get_points())[(0,) * sd]
                 nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,))
                              for d in range(sd)