tl;dr: defelement link
There is an equivalence between all polynomial basis which span the same space so in theory you can use the standard monomial basis to fit the various parameters using a Vandermonde matrix.
The Bell elements have 6 DOF per vertex, and then an additional 1 DOF per edge for a total of 21 DOF.
The vertex DOF map:
- The vertex value $v_i$
- The 5 derivatives of order 2 or less $\partial_x v_i$, $\partial_y v_i$, $\partial_{xx} v_i$, $\partial_{xy} v_i$, and $\partial_{yy} v_i$
The face DOF map $\int_{e_i} \frac{\partial v}{\partial \hat{n}_i}$
Let's now take the standard degree 5 monomial basis. Numerically this might not behave so well, but mathematically it's equivalent. The Vandermonde matrix is the evaluation of each basis function at the various DOF:
$$
V_{ij} = \phi_j|_{\mathrm{DOF}_i}
$$
So for example, the monomial $x^3 y^2$ evaluated for the DOF $v_2 = v(0,1) = 0^3 \cdot 1^2 = 0$, and for the DOF $\partial_{xy} v_1 = \partial_{xy} v(0,1) = 3 \cdot 0^2 \cdot 2 \cdot 1 = 0$.
The remaining trick is evaluating the monomial basis for the edge DOFs.
These need to compute the directional derivative with respect to the edge normal and then integrate this along the edge.
For edge 0 which runs between $(1,0)$ and $(0,1)$ with $\hat{n}_0 = \left(\frac{-1}{\sqrt{2}},\frac{-1}{\sqrt{2}}\right)$ we get
$$
-\int_{e_0} \frac{1}{\sqrt{2}} \partial_x \phi + \frac{1}{\sqrt{2}} \partial_y \phi ds
$$
For edge 1 which runs between $(0,0)$ and $(0,1)$ with $\hat{n}_1 = (-1,0)$ we get
$$
-\int_{e_1} \partial_x \phi ds
$$
For edge 2 which runs between $(0,0)$ and $(1,0)$ with $\hat{n}_2 = (0,1)$ we get
$$
\int_{e_2} \partial_y \phi ds
$$
Once we know the Vandermonde matrix, we simply need to invert it to get the inverse mapping and consequently the coefficients of each monomial.
The tl;dr link at the top has these listed out for you, and here is a Python code using SymPy which does this analytically:
import sympy
x, y, t = sympy.symbols(["x", "y", "t"])
# Monomial basis
basis = [
1+0*x,
x,
x**2,
x**3,
x**4,
x**5,
y,
x * y,
x**2 * y,
x**3 * y,
x**4 * y,
y**2,
x * y**2,
x**2 * y**2,
x**3 * y**2,
y**3,
x * y**3,
x**2 * y**3,
y**4,
x * y**4,
y**5,
]
# various derivatives of the Monomial basis
dx_basis = [sympy.diff(b, x) for b in basis]
dy_basis = [sympy.diff(b, y) for b in basis]
dxx_basis = [sympy.diff(b, (x, 2)) for b in basis]
dxy_basis = [sympy.diff(b, x, y) for b in basis]
dyy_basis = [sympy.diff(b, (y, 2)) for b in basis]
# evaluate the basis on the edges
edge0_vals = [
-sympy.integrate((dx + dy).subs({x: t, y: 1 - t}) / sympy.sqrt(2), (t, 0, 1))
for dx, dy in zip(dx_basis, dy_basis)
]
edge1_vals = [-sympy.integrate((dx).subs({x: 0, y: t}), (t, 0, 1)) for dx in dx_basis]
edge2_vals = [sympy.integrate((dy).subs({x: t, y: 0}), (t, 0, 1)) for dy in dy_basis]
# now build the Vandermonde matrix
V = sympy.Matrix([[0]*len(basis) for b in basis])
for i in range(len(basis)):
V[i,0] = basis[i].subs({x:0,y:0})
V[i,1] = dx_basis[i].subs({x:0,y:0})
V[i,2] = dy_basis[i].subs({x:0,y:0})
V[i,3] = dxx_basis[i].subs({x:0,y:0})
V[i,4] = dxy_basis[i].subs({x:0,y:0})
V[i,5] = dyy_basis[i].subs({x:0,y:0})
V[i,6] = basis[i].subs({x:1,y:0})
V[i,7] = dx_basis[i].subs({x:1,y:0})
V[i,8] = dy_basis[i].subs({x:1,y:0})
V[i,9] = dxx_basis[i].subs({x:1,y:0})
V[i,10] = dxy_basis[i].subs({x:1,y:0})
V[i,11] = dyy_basis[i].subs({x:1,y:0})
V[i,12] = basis[i].subs({x:0,y:1})
V[i,13] = dx_basis[i].subs({x:0,y:1})
V[i,14] = dy_basis[i].subs({x:0,y:1})
V[i,15] = dxx_basis[i].subs({x:0,y:1})
V[i,16] = dxy_basis[i].subs({x:0,y:1})
V[i,17] = dyy_basis[i].subs({x:0,y:1})
V[i,18] = edge0_vals[i]
V[i,19] = edge1_vals[i]
V[i,20] = edge2_vals[i]
Vinv = V.inv()
# now expand to the polynomials
pbasis = [sympy.simplify(sum([Vinv[i,j] * basis[j] for j in range(len(basis))])) for i in range(len(basis))]