diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 43d6184e..3049d44f 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -81,14 +81,12 @@ 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
 
     num_members = math.comb(n + dim, dim)
 
     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)
@@ -98,7 +96,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
                for k in range(order+1)]
 
     phi, dphi, ddphi = results + [None] * (2-order)
-    phi[0] += scale
+    phi[0] = scale
     if dim == 0 or n == 0:
         return results
     if dim > 3 or dim < 0:
@@ -127,29 +125,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] = phi[icur] * dfcur
+                dphi[inext] += fcur * dphi[icur]
                 if ddphi is not None:
-                    ddphi[inext] = factor * ddphi[icur] + sym_outer(dphi[icur], dfactor)
+                    ddphi[inext] = outer(dphi[icur], dfcur)
+                    ddphi[inext] += outer(dfcur, dphi[icur])
+                    ddphi[inext] += fcur * ddphi[icur]
 
             # 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] = phi[icur] * dfcur
+                dphi[inext] += phi[iprev] * dfprev
+                dphi[inext] += fcur * dphi[icur]
+                dphi[inext] += fprev * dphi[iprev]
                 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] = phi[iprev] * ddfprev
+                ddphi[inext] += outer(dphi[icur], dfcur)
+                ddphi[inext] += outer(dfcur, dphi[icur])
+                ddphi[inext] += outer(dphi[iprev], dfprev)
+                ddphi[inext] += outer(dfprev, dphi[iprev])
+                ddphi[inext] += fcur * ddphi[icur]
+                ddphi[inext] += fprev * ddphi[iprev]
 
         # normalize
         d = codim + 1
diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py
index eead0cfa..cb605ec8 100644
--- a/FIAT/lagrange.py
+++ b/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)
diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py
index 1b258e6d..75ff407c 100644
--- a/FIAT/nodal_enriched.py
+++ b/FIAT/nodal_enriched.py
@@ -12,6 +12,7 @@
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import CiarletElement
 from FIAT.barycentric_interpolation import LagrangeLineExpansionSet
+from FIAT.quadrature_schemes import create_quadrature
 
 __all__ = ['NodalEnrichedElement']
 
@@ -50,6 +51,7 @@ def __init__(self, *elements):
         expansion_set = elem.get_nodal_basis().get_expansion_set()
         mapping = elem.mapping()[0]
         value_shape = elem.value_shape()
+        sd = ref_el.get_spatial_dimension()
 
         # Sanity check
         assert all(e.get_reference_complex() == ref_el for e in elements)
@@ -61,11 +63,21 @@ def __init__(self, *elements):
             # Obtain coefficients via interpolation
             points = expansion_set.get_points()
             coeffs = np.vstack([e.tabulate(0, points)[(0,)] for e in elements])
-        else:
-            assert all(e.get_nodal_basis().get_expansion_set() == expansion_set
-                       for e in elements)
+        elif all(e.get_nodal_basis().get_expansion_set() == expansion_set for e in elements):
+            # All elements have the same ExpansionSet, just merge coefficients
             coeffs = [e.get_coeffs() for e in elements]
             coeffs = _merge_coeffs(coeffs, ref_el, embedded_degrees, expansion_set.continuity)
+        else:
+            # Obtain coefficients via projection
+            Q = create_quadrature(ref_el, 2*embedded_degree)
+            qpts, qwts = Q.get_points(), Q.get_weights()
+            phis = expansion_set._tabulate(embedded_degree, qpts, 0)[(0,)*sd]
+            PhiW = phis * qwts[None, :]
+            M = np.tensordot(PhiW, phis, (-1, -1))
+            MinvPhiW = np.linalg.solve(M, PhiW)
+            coeffs = [np.tensordot(e.tabulate(0, qpts)[(0,)*sd], MinvPhiW, (-1, -1)) for e in elements]
+            coeffs = np.concatenate(coeffs, axis=0)
+
         poly_set = PolynomialSet(ref_el,
                                  degree,
                                  embedded_degree,