refactor bubbles, add a test

FIAT/expansions.py

@@ -24,7 +24,7 @@ def morton_index3(p, q=0, r=0):
 def jrc(a, b, n):
     """Jacobi recurrence coefficients"""
     an = (2*n+1+a+b)*(2*n+2+a+b) / (2*(n+1)*(n+1+a+b))
-    bn = (a*a-b*b) * (2*n+1+a+b) / (2*(n+1)*(2*n+a+b)*(n+1+a+b))
+    bn = (a+b)*(a-b)*(2*n+1+a+b) / (2*(n+1)*(n+1+a+b)*(2*n+a+b))
     cn = (n+a)*(n+b)*(2*n+2+a+b) / ((n+1)*(n+1+a+b)*(2*n+a+b))
     return an, bn, cn
 
@@ -36,9 +36,10 @@ def bubble_jrc(a, b, n):
         bn = (a - 3*b - 2) / 4
         cn = 0.0
     else:
+        # an, bn, cn = jrc(a-1, b+1, n-1)
         # TODO check dependence on b
-        an = (2*n-1+a+b) * (2*n+a+b) / (2*(n+1)*(n+a+b))
-        bn = (a+b)*(a-b-2) * (2*n-1+a+b) / (2*(n+1)*(2*n-2+a+b)*(n+a+b))
+        an = (2*n-1+a+b)*(2*n+a+b) / (2*(n+1)*(n+a+b))
+        bn = (a+b)*(a-b-2)*(2*n-1+a+b) / (2*(n+1)*(n+a+b)*(2*n-2+a+b))
         cn = (n+a-2)*(n+b-1)*(2*n+a+b) / ((n+1)*(n+a+b)*(2*n-2+a+b))
     return an, bn, cn
 
@@ -98,13 +99,14 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
         fa, fb, fc, dfa, dfb, dfc = jacobi_factors(*X[codim:codim+3], *dX[codim:codim+3])
         ddfc = 2 * outer(dfb, dfb)
         for sub_index in reference_element.lattice_iter(0, n, codim):
-            alpha = 2 * sum(sub_index) + len(sub_index)
             # handle i = 1
             icur = idx(*sub_index, 0)
             inext = idx(*sub_index, 1)
             if bubble:
+                alpha = 2 * sum(sub_index)
                 a = b = 1.0
             else:
+                alpha = 2 * sum(sub_index) + len(sub_index)
                 a = 0.5 * alpha + 1.0
                 b = 0.5 * alpha
             factor = a * fa - b * fb
@@ -317,7 +319,6 @@ def _tabulate(self, n, pts, order=0):
             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]

FIAT/hierarchical.py

@@ -14,8 +14,7 @@
 from FIAT.barycentric_interpolation import make_dmat
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
-from FIAT.polynomial_set import ONPolynomialSet, mis
-from FIAT.expansions import morton_index2, morton_index3, polynomial_dimension
+from FIAT.polynomial_set import ONPolynomialSet, make_bubbles
 
 
 class LegendreDual(dual_set.DualSet):
@@ -59,12 +58,12 @@ def __init__(self, ref_el, degree, poly_set):
         entity_permutations = {}
 
         # vertex dofs
-        top = ref_el.get_topology()
-        nodes = [functional.PointEvaluation(ref_el, pt) for pt in ref_el.vertices]
-        nvertices = len(top[0])
-        entity_ids[0] = {k: [k] for k in range(nvertices)}
-        entity_permutations[0] = {k: {0: [0]} for k in range(nvertices)}
+        vertices = ref_el.get_vertices()
+        nodes = [functional.PointEvaluation(ref_el, pt) for pt in vertices]
+        entity_ids[0] = {k: [k] for k in range(len(vertices))}
+        entity_permutations[0] = {k: {0: [0]} for k in range(len(vertices))}
 
+        top = ref_el.get_topology()
         for dim in range(1, len(top)):
             entity_ids[dim] = {}
             entity_permutations[dim] = {}
@@ -73,7 +72,8 @@ def __init__(self, ref_el, degree, poly_set):
             quad_degree = 2 * degree
             ref_facet = ufc_simplex(dim)
             Q_ref = create_quadrature(ref_facet, quad_degree)
-            phis = self._tabulate_bubbles(ref_facet, degree, Q_ref, P=poly_set if dim == len(top)-1 else None)
+            P = poly_set if dim == len(top)-1 else None
+            phis = self._tabulate_L2_duals(ref_facet, degree, Q_ref, poly_set=P)
 
             for entity in range(len(top[dim])):
                 cur = len(nodes)
@@ -84,51 +84,33 @@ def __init__(self, ref_el, degree, poly_set):
 
         super(IntegratedLegendreDual, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
 
-    def _tabulate_bubbles(self, ref_el, degree, Q, P=None):
+    def _tabulate_L2_duals(self, ref_el, degree, Q, poly_set=None):
         qpts = Q.get_points()
+        qwts = Q.get_weights()
         dim = ref_el.get_spatial_dimension()
-        k = degree - dim - 1
-        if k < 0:
-            return numpy.zeros((0, len(qpts)))
-
-        if P is None:
-            P = ONPolynomialSet(ref_el, degree, bubble=True)
-        if dim == 1:
-            indices = list(range(2, degree+1))
-        else:
-            idx = (morton_index2, morton_index3)[dim-2]
-            indices = []
-            for p in range(1, degree+1):
-                for alpha in mis(dim, p):
-                    if alpha[0] > 1 and min(alpha[1:]) > 0:
-                        indices.append(idx(*alpha))
-
-        assert len(indices) == polynomial_dimension(ref_el, k)
-        bubbles = P.take(indices)
-        bubbles_table = bubbles.tabulate(qpts)[(0,) * dim]
-
-        W = Q.get_weights()
+        moments = lambda v: numpy.dot(numpy.multiply(v, qwts), v.T)
         if dim == 1:
+            # Assemble a stiffness matrix in the Lagrange basis
             dmat, _ = make_dmat(qpts.flatten())
-            K = numpy.dot(numpy.multiply(dmat, W), dmat.T)
+            K = numpy.dot(numpy.multiply(dmat, qwts), dmat.T)
         else:
             # Get ON basis
             P = ONPolynomialSet(ref_el, degree)
             P_table = P.tabulate(qpts, 1)
             # Assemble a stiffness matrix in the ON basis
-            moments = lambda dv: numpy.dot(numpy.multiply(dv, W), dv.T)
             K = sum(moments(P_table[alpha]) for alpha in P_table if sum(alpha) == 1)
             # Change of basis to Lagrange polynomials at the quadrature nodes
-            v = numpy.multiply(P_table[(0, ) * dim], W)
+            v = numpy.multiply(P_table[(0, ) * dim], qwts)
             K = numpy.dot(numpy.dot(v.T, K), v)
 
-        phis = numpy.multiply(numpy.dot(bubbles_table, K), 1/W)
-        phis = numpy.concatenate([phis[1::2], phis[0::2]])
+        B = make_bubbles(ref_el, degree, poly_set=poly_set)
+        B_at_qpts = B.tabulate(qpts)[(0,) * dim]
+        phis = numpy.multiply(numpy.dot(B_at_qpts, K), 1/qwts)
         return phis
 
 
 class IntegratedLegendre(finite_element.CiarletElement):
-    """1D continuous element with integrated Legendre polynomials."""
+    """Simplicial continuous element with integrated Legendre polynomials."""
 
     def __init__(self, ref_el, degree):
         if ref_el.shape not in {POINT, LINE, TRIANGLE, TETRAHEDRON}:

FIAT/polynomial_set.py

@@ -241,3 +241,28 @@ def __init__(self, ref_el, degree, size=None):
 
         super(ONSymTensorPolynomialSet, self).__init__(ref_el, degree, embedded_degree,
                                                        expansion_set, coeffs)
+
+
+def make_bubbles(ref_el, degree, shape=(), poly_set=None):
+    """Construct a polynomial set with bubbles up to the given degree.
+
+    """
+    from itertools import chain
+
+    dim = ref_el.get_spatial_dimension()
+    degrees = chain(range(3, degree+1, 2), range(2, degree+1, 2))
+    if dim == 1:
+        indices = list(degrees)
+    else:
+        idx = (expansions.morton_index2, expansions.morton_index3)[dim-2]
+        indices = []
+        for p in degrees:
+            for alpha in mis(dim, p):
+                if alpha[0] > 1 and min(alpha[1:]) > 0:
+                    indices.append(idx(*alpha))
+
+    assert len(indices) == expansions.polynomial_dimension(ref_el, degree - dim - 1)
+    if poly_set is None:
+        poly_set = ONPolynomialSet(ref_el, degree, shape=shape, bubble=True)
+    bubbles = poly_set.take(indices)
+    return bubbles

test/unit/test_fiat.py

@@ -590,11 +590,7 @@ def basis(dim, p, q=0, r=0):
         return f
 
     def eval_basis(f, pt):
-        fval = f
-        for coord, pval in zip(eta, duffy_coords(pt)):
-            fval = fval.subs(coord, pval)
-        fval = float(fval)
-        return fval
+        return float(f.subs(dict(zip(eta, duffy_coords(pt)))))
 
     for i in range(n + 1):
         for indices in polynomial_set.mis(dim, i):
@@ -604,6 +600,26 @@ def eval_basis(f, pt):
             assert np.allclose(uh, exact, atol=1E-14)
 
 
+@pytest.mark.parametrize('cell', [I, T, S])
+def test_make_bubbles(cell):
+    from FIAT.expansions import polynomial_dimension
+    from FIAT.polynomial_set import make_bubbles, PolynomialSet
+    degree = 10
+
+    top = cell.get_topology()
+    points = []
+    for dim in range(len(top)-1):
+        for entity in range(len(top[dim])):
+            points.extend(cell.make_points(dim, entity, degree, variant="gl"))
+
+    sd = cell.get_spatial_dimension()
+    B = make_bubbles(cell, degree)
+    assert isinstance(B, PolynomialSet)
+    assert B.degree == degree
+    assert B.get_num_members() == polynomial_dimension(cell, degree - sd - 1)
+    assert np.allclose(B.tabulate(points)[(0,)*sd], 0)
+
+
 if __name__ == '__main__':
     import os
     pytest.main(os.path.abspath(__file__))