moments against L2 duals fo bubbles

FIAT/expansions.py

@@ -29,18 +29,17 @@ def jrc(a, b, n):
     return an, bn, cn
 
 
-def bubble_jrc(a, b, n):
+def integrated_jrc(a, b, n):
     """Integrated Jacobi recurrence coefficients"""
     if n == 1:
         an = (a + b + 2) / 4
         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)*(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))
+        an, bn, cn = jrc(a-1, b+1, n-1)
+        an *= n / (n+1)
+        bn *= n / (n+1)
+        cn *= (n-1) / (n+1)
     return an, bn, cn
 
 
@@ -67,7 +66,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, bubble=False):
+def dubiner_recurrence(dim, n, order, ref_pts, jacobian, variant=None):
     """Dubiner recurrence from (Kirby 2010)"""
     if order > 2:
         raise ValueError("Higher order derivatives not supported")
@@ -91,7 +90,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
     if dim > 3 or dim < 0:
         raise ValueError("Invalid number of spatial dimensions")
 
-    coefficients = bubble_jrc if bubble else jrc
+    coefficients = integrated_jrc if variant == "integral" else jrc
     X = pad_coordinates(ref_pts, pad_dim)
     idx = (lambda p: p, morton_index2, morton_index3)[dim-1]
     for codim in range(dim):
@@ -102,13 +101,15 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
             # handle i = 1
             icur = idx(*sub_index, 0)
             inext = idx(*sub_index, 1)
-            if bubble:
+
+            beta = 0
+            if variant == "integral":
                 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
+                a = 0.5 * (alpha + beta) + 1.0
+                b = 0.5 * (alpha - beta)
             factor = a * fa - b * fb
             phi[inext] = factor * phi[icur]
             if dphi is not None:
@@ -120,7 +121,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
             # general i by recurrence
             for i in range(1, n - sum(sub_index)):
                 iprev, icur, inext = icur, inext, idx(*sub_index, i + 1)
-                a, b, c = coefficients(alpha, 0, i)
+                a, b, c = coefficients(alpha, beta, i)
                 factor = a * fa - b * fb
                 phi[inext] = factor * phi[icur] - c * (fc * phi[iprev])
                 if dphi is None:
@@ -135,8 +136,9 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
 
         # normalize
         for index in reference_element.lattice_iter(0, n+1, codim+1):
+            alpha = 2 * sum(index) + len(index)
+            scale = math.sqrt(0.5 * alpha)
             icur = idx(*index)
-            scale = math.sqrt(sum(index) + 0.5 * len(index))
             for result in results:
                 result[icur] *= scale
     return results
@@ -177,9 +179,9 @@ def __new__(cls, *args, **kwargs):
         else:
             raise ValueError("Invalid reference element type.")
 
-    def __init__(self, ref_el, bubble=False):
+    def __init__(self, ref_el, variant=None):
         self.ref_el = ref_el
-        self.bubble = bubble
+        self.variant = variant
         dim = ref_el.get_spatial_dimension()
         self.base_ref_el = reference_element.default_simplex(dim)
         v1 = ref_el.get_vertices()
@@ -206,7 +208,7 @@ def _tabulate(self, n, pts, order=0):
         """A version of tabulate() that also works for a single point.
         """
         D = self.ref_el.get_spatial_dimension()
-        return dubiner_recurrence(D, n, order, self._mapping(pts), self.A, bubble=self.bubble)
+        return dubiner_recurrence(D, n, order, self._mapping(pts), self.A, variant=self.variant)
 
     def get_dmats(self, degree):
         """Returns a numpy array with the expansion coefficients dmat[k, j, i]
@@ -309,7 +311,7 @@ def __init__(self, ref_el, **kwargs):
     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]."""
-        if self.bubble:
+        if self.variant is not None:
             return super(LineExpansionSet, self)._tabulate(n, pts, order=order)
 
         xs = self._mapping(pts).T

FIAT/hierarchical.py

@@ -69,11 +69,13 @@ def __init__(self, ref_el, degree, poly_set):
             entity_permutations[dim] = {}
             perms = make_entity_permutations_simplex(dim, degree - dim)
 
-            quad_degree = 2 * degree
             ref_facet = ufc_simplex(dim)
-            Q_ref = create_quadrature(ref_facet, quad_degree)
-            P = poly_set if dim == len(top)-1 else None
-            phis = self._tabulate_L2_duals(ref_facet, degree, Q_ref, poly_set=P)
+            if dim == 1:
+                Q_ref = create_quadrature(ref_facet, 2 * degree)
+                phis = self._tabulate_H1_duals(ref_facet, degree, Q_ref)
+            else:
+                Q_ref = create_quadrature(ref_facet, max(0, 2 * degree - dim - 1))
+                phis = self._tabulate_L2_duals(ref_facet, degree, Q_ref)
 
             for entity in range(len(top[dim])):
                 cur = len(nodes)
@@ -84,15 +86,16 @@ def __init__(self, ref_el, degree, poly_set):
 
         super(IntegratedLegendreDual, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
 
-    def _tabulate_L2_duals(self, ref_el, degree, Q, poly_set=None):
+    def _tabulate_H1_duals(self, ref_el, degree, Q):
         qpts = Q.get_points()
         qwts = Q.get_weights()
-        dim = ref_el.get_spatial_dimension()
         moments = lambda v: numpy.dot(numpy.multiply(v, qwts), v.T)
+
+        dim = ref_el.get_spatial_dimension()
         if dim == 1:
             # Assemble a stiffness matrix in the Lagrange basis
             dmat, _ = make_dmat(qpts.flatten())
-            K = numpy.dot(numpy.multiply(dmat, qwts), dmat.T)
+            K = moments(dmat)
         else:
             # Get ON basis
             P = ONPolynomialSet(ref_el, degree)
@@ -103,9 +106,17 @@ def _tabulate_L2_duals(self, ref_el, degree, Q, poly_set=None):
             v = numpy.multiply(P_table[(0, ) * dim], qwts)
             K = numpy.dot(numpy.dot(v.T, K), v)
 
-        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)
+        B = make_bubbles(ref_el, degree)
+        phis = B.tabulate(qpts)[(0,) * dim]
+        phis = numpy.multiply(numpy.dot(phis, K), 1/qwts)
+        return phis
+
+    def _tabulate_L2_duals(self, ref_el, degree, Q):
+        dim = ref_el.get_spatial_dimension()
+        B = make_bubbles(ref_el, degree)
+        phis = B.tabulate(Q.pts)[(0,) * dim]
+        if len(phis) > 0:
+            phis = phis / phis[0]
         return phis
 
 
@@ -115,7 +126,7 @@ class IntegratedLegendre(finite_element.CiarletElement):
     def __init__(self, ref_el, degree):
         if ref_el.shape not in {POINT, LINE, TRIANGLE, TETRAHEDRON}:
             raise ValueError("%s is only defined on simplices." % type(self))
-        poly_set = ONPolynomialSet(ref_el, degree, bubble=True)
+        poly_set = ONPolynomialSet(ref_el, degree, variant="integral")
         dual = IntegratedLegendreDual(ref_el, degree, poly_set)
         formdegree = 0  # 0-form
         super(IntegratedLegendre, self).__init__(poly_set, dual, degree, formdegree)

FIAT/polynomial_set.py

@@ -16,6 +16,7 @@
 # an entire set of polynomials)
 
 import numpy
+from itertools import chain
 from FIAT import expansions
 from FIAT.functional import index_iterator
 
@@ -120,7 +121,7 @@ class ONPolynomialSet(PolynomialSet):
 
     """
 
-    def __init__(self, ref_el, degree, shape=tuple(), bubble=False):
+    def __init__(self, ref_el, degree, shape=tuple(), variant=None):
 
         if shape == tuple():
             num_components = 1
@@ -130,7 +131,7 @@ def __init__(self, ref_el, degree, shape=tuple(), bubble=False):
         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, bubble=bubble)
+        expansion_set = expansions.ExpansionSet(ref_el, variant=variant)
 
         # set up coefficients
         if shape == tuple():
@@ -243,14 +244,14 @@ def __init__(self, ref_el, degree, size=None):
                                                        expansion_set, coeffs)
 
 
-def make_bubbles(ref_el, degree, shape=(), poly_set=None):
+def make_bubbles(ref_el, degree, shape=()):
     """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))
+    poly_set = ONPolynomialSet(ref_el, degree, shape=shape, variant="integral")
+    degrees = chain(range(dim + 1, degree+1, 2), range(dim + 2, degree+1, 2))
+
     if dim == 1:
         indices = list(degrees)
     else:
@@ -260,9 +261,5 @@ def make_bubbles(ref_el, degree, shape=(), poly_set=None):
             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
+    poly_set = poly_set.take(indices)
+    return poly_set

test/unit/test_fiat.py

@@ -603,12 +603,15 @@ def eval_basis(f, pt):
 @pytest.mark.parametrize('cell', [I, T, S])
 def test_make_bubbles(cell):
     from FIAT.reference_element import make_lattice
+    from FIAT.quadrature_schemes import create_quadrature
     from FIAT.expansions import polynomial_dimension
-    from FIAT.polynomial_set import make_bubbles, PolynomialSet
+    from FIAT.polynomial_set import make_bubbles, PolynomialSet, ONPolynomialSet
 
     degree = 10
-    sd = cell.get_spatial_dimension()
     B = make_bubbles(cell, degree)
+
+    # basic tests
+    sd = cell.get_spatial_dimension()
     assert isinstance(B, PolynomialSet)
     assert B.degree == degree
     assert B.get_num_members() == polynomial_dimension(cell, degree - sd - 1)
@@ -629,6 +632,17 @@ def test_make_bubbles(cell):
     assert values.shape == (m, m)
     assert np.linalg.matrix_rank(values.T) == m
 
+    # test that B intersected with span(P_{degree+2} \ P_{degree}) is empty
+    P = ONPolynomialSet(cell, degree + 2)
+    P = P.take(list(range(polynomial_dimension(cell, degree),
+                          P.get_num_members())))
+
+    Q = create_quadrature(cell, P.degree + B.degree)
+    qpts = Q.get_points()
+    qwts = Q.get_weights()
+    P_at_qpts = P.tabulate(qpts)[(0,) * sd]
+    B_at_qpts = B.tabulate(qpts)[(0,) * sd]
+    assert np.allclose(np.dot(np.multiply(P_at_qpts, qwts), B_at_qpts.T), 0.0)
 
 
 if __name__ == '__main__':