Fix ONPolynomialSet normalization

FIAT/expansions.py

@@ -10,7 +10,6 @@
 import numpy
 import math
 from FIAT import reference_element
-from FIAT import jacobi
 
 
 def morton_index2(p, q=0):
@@ -52,7 +51,7 @@ def jacobi_factors(x, y, z, dx, dy, dz):
     return fa, fb, fc, dfa, dfb, dfc
 
 
-def dubiner_recurrence(dim, n, order, ref_pts, jacobian):
+def dubiner_recurrence(dim, n, order, ref_pts, Jinv):
     """Dubiner recurrence from (Kirby 2010)"""
     if order > 2:
         raise ValueError("Higher order derivatives not supported")
@@ -65,8 +64,11 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian):
     sym_outer = lambda x, y: outer(x, y) + outer(y, x)
 
     pad_dim = dim + 2
-    dX = pad_jacobian(jacobian, pad_dim)
-    phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 1.)
+    dX = pad_jacobian(Jinv, pad_dim)
+    ref_vol = 2 ** dim / math.factorial(dim)
+    scale = math.sqrt(abs(numpy.linalg.det(Jinv)) / ref_vol)
+    phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), scale)
+
     if dphi is not None:
         dphi[0] = (phi[0] - phi[0]) * dX[0]
     if ddphi is not None:
@@ -115,9 +117,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian):
                                 c * (fc * ddphi[iprev] + sym_outer(dphi[iprev], dfc) + phi[iprev] * ddfc))
 
         # normalize
-        for alpha in reference_element.lattice_iter(0, n+1, codim+1):
-            icur = idx(*alpha)
-            scale = math.sqrt(sum(alpha) + 0.5 * len(alpha))
+        d = codim + 1
+        for index in reference_element.lattice_iter(0, n+1, d):
+            icur = idx(*index)
+            scale = math.sqrt((2.0 * sum(index) + d) / d)
             for result in results:
                 result[icur] *= scale
     return results
@@ -284,24 +287,6 @@ def __init__(self, ref_el):
             raise Exception("Must have a line")
         super(LineExpansionSet, self).__init__(ref_el)
 
-    def _tabulate(self, n, pts, order=0):
-        """Returns a tuple of (vals, derivs) such that
-        vals[i,j] = phi_i(pts[j]), derivs[i,j] = D vals[i,j]."""
-        xs = self._mapping(pts).T
-        results = []
-        scale = numpy.sqrt(0.5 + numpy.arange(n+1))
-        for k in range(order+1):
-            v = numpy.zeros((n + 1, len(xs)), xs.dtype)
-            if n >= k:
-                v[k:] = jacobi.eval_jacobi_batch(k, k, n-k, xs)
-            for p in range(n + 1):
-                v[p] *= scale[p]
-                scale[p] *= 0.5 * (p + k + 1) * self.A[0, 0]
-            shape = v.shape
-            shape = shape[:1] + (1,) * k + shape[1:]
-            results.append(v.reshape(shape))
-        return tuple(results)
-
 
 class TriangleExpansionSet(ExpansionSet):
     """Evaluates the orthonormal Dubiner basis on a triangular

test/unit/test_fiat.py

@@ -530,16 +530,21 @@ 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, quadrature
-    U = expansions.ExpansionSet(cell)
-    degree = 10
-    rule = quadrature.make_quadrature(cell, degree + 1)
-    phi = U.tabulate(degree, rule.pts)
-    w = rule.get_weights()
-    scale = 0.5 ** -cell.get_spatial_dimension()
-    results = scale * np.dot(phi, w[:, None] * phi.T)
+@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
+    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]))