Define P1 with trivial Vandermonde matrix

Author: pbrubeck

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
 
@@ -94,9 +93,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
 
     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 +127,32 @@ 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] + phi[icur] * dfcur
                 if ddphi is not None:
-                    ddphi[inext] = factor * ddphi[icur] + sym_outer(dphi[icur], dfactor)
+                    ddphi[inext] = fcur * ddphi[icur] + sym_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] + 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] + phi[icur] * dfcur +
+                               fprev * dphi[iprev] + 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))
+                ddfcur = -c * ddfc
+                ddphi[inext] = (fcur * ddphi[icur] + sym_outer(dphi[icur], dfcur) +
+                                fprev * ddphi[iprev] + sym_outer(dphi[iprev], dfcur) + phi[iprev] * ddfcur)
 
         # normalize
         d = codim + 1
@@ -159,9 +162,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
             if variant is not None:
                 p = index[-1] + shift
                 alpha = 2 * (sum(index[:-1]) + d * shift) - 1
-                norm2 = (0.5 + d) / d
                 if p > 0 and p + alpha > 0:
-                    norm2 *= (p + alpha) * (2*p + alpha) / p
+                    norm2 = (2*d+1) * (p + alpha) * (2*p + alpha) / (2*d*p)
+                else:
+                    norm2 = 1.0
             else:
                 norm2 = (2*sum(index) + d) / d
             scale = math.sqrt(norm2)
@@ -183,9 +187,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)