add bubbles from Beuchler and Schoberl

FIAT/expansions.py

@@ -29,6 +29,20 @@ def jrc(a, b, n):
     return an, bn, cn
 
 
+def bubble_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:
+        # 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))
+        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
+
+
 def pad_coordinates(ref_pts, embedded_dim):
     """Pad reference coordinates by appending -1.0."""
     return tuple(ref_pts) + (-1.0, )*(embedded_dim - len(ref_pts))
@@ -52,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, alpha=0, beta=0):
+def dubiner_recurrence(dim, n, order, ref_pts, jacobian, bubble=False):
     """Dubiner recurrence from (Kirby 2010)"""
     if order > 2:
         raise ValueError("Higher order derivatives not supported")
@@ -76,24 +90,23 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, alpha=0, beta=0):
     if dim > 3 or dim < 0:
         raise ValueError("Invalid number of spatial dimensions")
 
+    coefficients = bubble_jrc if bubble else jrc
     X = pad_coordinates(ref_pts, pad_dim)
     idx = (lambda p: p, morton_index2, morton_index3)[dim-1]
-    base_alpha = alpha
     for codim in range(dim):
         # Extend the basis from codim to codim + 1
         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 = base_alpha + 2 * sum(sub_index) + len(sub_index)
+            alpha = 2 * sum(sub_index) + len(sub_index)
             # handle i = 1
             icur = idx(*sub_index, 0)
             inext = idx(*sub_index, 1)
-            apb = alpha + beta
-            if apb == 0 or apb == -1:
-                a = 0.5 * (alpha + beta) + 1.0
-                b = 0.5 * (alpha - beta)
+            if bubble:
+                a = b = 1.0
             else:
-                a, b, c = jrc(alpha, beta, 0)
+                a = 0.5 * alpha + 1.0
+                b = 0.5 * alpha
             factor = a * fa - b * fb
             phi[inext] = factor * phi[icur]
             if dphi is not None:
@@ -105,7 +118,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian, alpha=0, beta=0):
             # 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 = jrc(alpha, beta, i)
+                a, b, c = coefficients(alpha, 0, i)
                 factor = a * fa - b * fb
                 phi[inext] = factor * phi[icur] - c * (fc * phi[iprev])
                 if dphi is None:
@@ -162,17 +175,17 @@ def __new__(cls, *args, **kwargs):
         else:
             raise ValueError("Invalid reference element type.")
 
-    def __init__(self, ref_el, alpha=0, beta=0):
+    def __init__(self, ref_el, bubble=False):
         self.ref_el = ref_el
-        self.alpha = alpha
-        self.beta = beta
+        self.bubble = bubble
         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(numpy.linalg.det(self.A))
+        detA = numpy.sqrt(numpy.linalg.det(numpy.dot(self.A.T, self.A)))
+        self.scale = numpy.sqrt(detA)
         self._dmats_cache = {}
 
     def get_num_members(self, n):
@@ -191,7 +204,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, alpha=self.alpha, beta=self.beta)
+        return dubiner_recurrence(D, n, order, self._mapping(pts), self.A, bubble=self.bubble)
 
     def get_dmats(self, degree):
         """Returns a numpy array with the expansion coefficients dmat[k, j, i]
@@ -294,13 +307,17 @@ 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:
+            return super(LineExpansionSet, self)._tabulate(n, pts, order=order)
+
         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(self.alpha+k, self.beta+k, n-k, xs)
+                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,7 +14,8 @@
 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
+from FIAT.polynomial_set import ONPolynomialSet, mis
+from FIAT.expansions import morton_index2, morton_index3, polynomial_dimension
 
 
 class LegendreDual(dual_set.DualSet):
@@ -72,7 +73,7 @@ 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)
+            phis = self._tabulate_bubbles(ref_facet, degree, Q_ref, P=poly_set if dim == len(top)-1 else None)
 
             for entity in range(len(top[dim])):
                 cur = len(nodes)
@@ -83,18 +84,28 @@ 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):
+    def _tabulate_bubbles(self, ref_el, degree, Q, P=None):
         qpts = Q.get_points()
         dim = ref_el.get_spatial_dimension()
         k = degree - dim - 1
         if k < 0:
             return numpy.zeros((0, len(qpts)))
-        P = ONPolynomialSet(ref_el, k, alpha=1, beta=1)
-        P_at_qpts = P.tabulate(qpts)[(0,) * dim]
 
-        x = qpts.T
-        bubble = (1.0 - numpy.sum(x, axis=0)) * numpy.prod(x, axis=0)
-        bubbles = P_at_qpts * bubble
+        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()
         if dim == 1:
@@ -103,15 +114,15 @@ def _tabulate_bubbles(self, ref_el, degree, Q):
         else:
             # Get ON basis
             P = ONPolynomialSet(ref_el, degree)
-            tab = P.tabulate(qpts, 1)
+            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(tab[alpha]) for alpha in tab if sum(alpha) == 1)
+            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(tab[(0, ) * dim], W)
+            v = numpy.multiply(P_table[(0, ) * dim], W)
             K = numpy.dot(numpy.dot(v.T, K), v)
 
-        phis = numpy.multiply(numpy.dot(bubbles, K), 1/W)
+        phis = numpy.multiply(numpy.dot(bubbles_table, K), 1/W)
         phis = numpy.concatenate([phis[1::2], phis[0::2]])
         return phis
 
@@ -122,7 +133,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)
+        poly_set = ONPolynomialSet(ref_el, degree, bubble=True)
         dual = IntegratedLegendreDual(ref_el, degree, poly_set)
         formdegree = 0  # 0-form
         super(IntegratedLegendre, self).__init__(poly_set, dual, degree, formdegree)

FIAT/polynomial_set.py

@@ -120,7 +120,7 @@ class ONPolynomialSet(PolynomialSet):
 
     """
 
-    def __init__(self, ref_el, degree, shape=tuple(), alpha=0, beta=0):
+    def __init__(self, ref_el, degree, shape=tuple(), bubble=False):
 
         if shape == tuple():
             num_components = 1
@@ -130,7 +130,7 @@ def __init__(self, ref_el, degree, shape=tuple(), alpha=0, beta=0):
         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, alpha=alpha, beta=beta)
+        expansion_set = expansions.ExpansionSet(ref_el, bubble=bubble)
 
         # set up coefficients
         if shape == tuple():