Define P1 with trivial Vandermonde matrix

FIAT/brezzi_douglas_fortin_marini.py

@@ -70,7 +70,7 @@ def BDFMSpace(ref_el, order):
     # Linear vector valued space. Since the embedding degree of this element
     # is 2, this is implemented by taking the quadratic space and selecting
     # the linear polynomials.
-    vec_poly_set = polynomial_set.ONPolynomialSet(ref_el, order, (sd,))
+    vec_poly_set = polynomial_set.ONPolynomialSet(ref_el, order, (sd,), variant="bubble")
     # Linears are the first three polynomials in each dimension.
     vec_poly_set = vec_poly_set.take([0, 1, 2, 6, 7, 8])
 

FIAT/expansions.py

@@ -81,7 +81,6 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
         raise ValueError("Higher order derivatives not supported")
     if variant not in [None, "bubble", "dual"]:
         raise ValueError(f"Invalid variant {variant}")
-
     if variant == "bubble":
         scale = -scale
 
@@ -90,13 +89,13 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
     phi, dphi, ddphi = results + (None,) * (2-order)
 
     outer = lambda x, y: x[:, None, ...] * y[None, ...]
-    sym_outer = lambda x, y: outer(x, y) + outer(y, x)
 
     pad_dim = dim + 2
     dX = pad_jacobian(Jinv, pad_dim)
-    phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), scale)
+    phi[0] = sum((ref_pts[i] * 0 for i in range(dim)), scale)
+
     if dphi is not None:
-        dphi[0] = (phi[0] - phi[0]) * dX[0]
+        dphi[0] = (phi[0] * 0) * dX[0]
     if ddphi is not None:
         ddphi[0] = outer(dphi[0], dX[0])
     if dim == 0 or n == 0:
@@ -127,29 +126,46 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
                 a = 0.5 * (alpha + beta) + 1.0
                 b = 0.5 * (alpha - beta)
 
-            factor = a * fa - b * fb
-            phi[inext] = factor * phi[icur]
+            fcur = a * fa - b * fb
+            phi[inext] = fcur * phi[icur]
             if dphi is not None:
-                dfactor = a * dfa - b * dfb
-                dphi[inext] = factor * dphi[icur] + phi[icur] * dfactor
+                dfcur = a * dfa - b * dfb
+                dphi[inext] = fcur * dphi[icur]
+                dphi[inext] += phi[icur] * dfcur
                 if ddphi is not None:
-                    ddphi[inext] = factor * ddphi[icur] + sym_outer(dphi[icur], dfactor)
+                    ddphi[inext] = fcur * ddphi[icur]
+                    ddphi[inext] += outer(dfcur, dphi[icur])
+                    ddphi[inext] += outer(dphi[icur], dfcur)
 
             # 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, beta, i)
-                factor = a * fa - b * fb
-                phi[inext] = factor * phi[icur] - c * (fc * phi[iprev])
+
+                fcur = a * fa - b * fb
+                fprev = -c * fc
+                phi[inext] = fcur * phi[icur]
+                phi[inext] += fprev * phi[iprev]
                 if dphi is None:
                     continue
-                dfactor = a * dfa - b * dfb
-                dphi[inext] = (factor * dphi[icur] + phi[icur] * dfactor -
-                               c * (fc * dphi[iprev] + phi[iprev] * dfc))
+
+                dfcur = a * dfa - b * dfb
+                dfprev = -c * dfc
+                dphi[inext] = fcur * dphi[icur]
+                dphi[inext] += fprev * dphi[iprev]
+                dphi[inext] += phi[icur] * dfcur
+                dphi[inext] += phi[iprev] * dfprev
                 if ddphi is None:
                     continue
-                ddphi[inext] = (factor * ddphi[icur] + sym_outer(dphi[icur], dfactor) -
-                                c * (fc * ddphi[iprev] + sym_outer(dphi[iprev], dfc) + phi[iprev] * ddfc))
+
+                ddfprev = -c * ddfc
+                ddphi[inext] = fcur * ddphi[icur]
+                ddphi[inext] += fprev * ddphi[iprev]
+                ddphi[inext] += outer(dfcur, dphi[icur])
+                ddphi[inext] += outer(dphi[icur], dfcur)
+                ddphi[inext] += outer(dfprev, dphi[iprev])
+                ddphi[inext] += outer(dphi[iprev], dfprev)
+                ddphi[inext] += phi[iprev] * ddfprev
 
         # normalize
         d = codim + 1
@@ -183,9 +199,9 @@ def C0_basis(dim, n, tabulations):
     # Recover facet bubbles
     for phi in tabulations:
         icur = 0
-        phi[icur] *= -1
         for inext in range(1, dim+1):
-            phi[icur] -= phi[inext]
+            phi[icur] += phi[inext]
+        phi[icur] *= -1
         if dim == 2:
             for i in range(2, n+1):
                 phi[idx(0, i)] -= phi[idx(1, i-1)]

FIAT/lagrange.py

@@ -83,7 +83,6 @@ def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False):
             points = get_lagrange_points(dual)
             poly_set = LagrangePolynomialSet(ref_el, points)
         else:
-            poly_variant = "bubble" if ref_el.is_macrocell() else None
-            poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant=poly_variant)
+            poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble", scale=1)
         formdegree = 0  # 0-form
         super().__init__(poly_set, dual, degree, formdegree)