Fix ONPolynomialSet normalization

FIAT/brezzi_douglas_marini.py

@@ -30,13 +30,12 @@ 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=1)
             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
+                n = ref_el.compute_normal(f)
                 phis = n[None, :, None] * Pq_at_qpts[:, None, :]
                 nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
                              for phi in phis)

FIAT/expansions.py

@@ -252,8 +252,8 @@ def __init__(self, ref_el, scale=None, variant=None):
         self.mapping = lambda x: numpy.dot(self.A, x) + self.b
         self._dmats_cache = {}
         if scale is None:
-            scale = math.sqrt(1.0 / self.base_ref_el.volume())
-        elif isinstance(scale, str):
+            scale = "orthonormal"
+        if isinstance(scale, str):
             scale = scale.lower()
             if scale == "orthonormal":
                 scale = math.sqrt(1.0 / ref_el.volume())

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

@@ -121,8 +121,8 @@ class ONPolynomialSet(PolynomialSet):
 
     """
 
-    def __init__(self, ref_el, degree, shape=tuple(), scale=None, variant=None):
 
+    def __init__(self, ref_el, degree, shape=tuple(), scale=None, variant=None):
         if shape == tuple():
             num_components = 1
         else:

FIAT/raviart_thomas.py

@@ -75,13 +75,12 @@ 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=1)
             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
+                n = ref_el.compute_normal(f)
                 phis = n[None, :, None] * Pq_at_qpts[:, None, :]
                 nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
                              for phi in phis)
@@ -91,7 +90,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=1)
                 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)

test/unit/test_fiat.py

@@ -530,20 +530,22 @@ def test_error_point_high_order(element):
         eval(element)
 
 
-@pytest.mark.parametrize('cell', [I, T, S])
-def test_expansion_orthonormality(cell):
-    from FIAT import expansions
+@pytest.mark.parametrize('degree', [0, 10])
+@pytest.mark.parametrize('dim', range(1, 4))
+@pytest.mark.parametrize('cell', ["default", "ufc"])
+def test_expansion_orthonormality(cell, dim, degree):
+    from FIAT.expansions import ExpansionSet
     from FIAT.quadrature_schemes import create_quadrature
-    U = expansions.ExpansionSet(cell)
-    degree = 10
-    rule = create_quadrature(cell, 2*degree)
-    phi = U.tabulate(degree, rule.pts)
-    qwts = rule.get_weights()
-    scale = 2 ** cell.get_spatial_dimension()
-    results = scale * np.dot(np.multiply(phi, qwts), phi.T)
-
-    assert np.allclose(results, np.diag(np.diag(results)))
-    assert np.allclose(np.diag(results), 1.0)
+    from FIAT import reference_element
+    make_cell = {"default": reference_element.default_simplex,
+                 "ufc": reference_element.ufc_simplex}[cell]
+    ref_el = make_cell(dim)
+    U = ExpansionSet(ref_el)
+    Q = create_quadrature(ref_el, 2 * degree)
+    qpts, qwts = Q.get_points(), Q.get_weights()
+    phi = U.tabulate(degree, qpts)
+    results = np.dot(np.multiply(phi, qwts), phi.T)
+    assert np.allclose(results, np.eye(results.shape[0]))
 
 
 @pytest.mark.parametrize('dim', range(1, 4))