C^1 continuous element for a triangle?

Question

I am looking for an element for FEM that is piecewise $C^1$ continuous across triangles (i.e. $C^1$ continuous on the edge separating 2 triangles of the mesh).

I have heard about the Bell element: https://arxiv.org/pdf/1706.09017.pdf https://defelement.com/elements/bell.html

But I am struggling finding a simple definition of the form (or similar):

Given the known points $(x_0, y_0) ... (x_n,y_n)$ and the known partial derivatives, or gradients, or whatever is necessary e.g. ($\frac{\partial p}{\partial x} (x_0, y_0)$).

The bell element is the quintic polynomial with coefficients $k_0 = \text{some function of } x_0,y_0, k_1= ...$

i.e. a simplified formulation that tells you what your inputs are and how to turn them into the coefficients for the polynomial.

Answer 1

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:

1. The vertex value $v_i$
2. 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))]

Answer 2

The $C^1$ elements are all very challenging to implement, which is why you don't see them get used very often. If you just want to see what the basis functions are on the reference (or other) triangles, you can try out symfem, which will spit out the basis functions as sympy expressions. This is what I got with a few minutes of playing around with it:

In [1]: import symfem

In [2]: help(symfem)


In [3]: bell = symfem.create_element("triangle", "Bell", 5)

In [4]: bell
Out[4]: <symfem.elements.bell.Bell at 0x7f699e217190>

In [5]: bell.get_polynomial_basis()
Out[5]: 
[1,
 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]

Symfem includes the Bell, Argyris, and Hsieh-Clough-Tocher elements if you want to try those out too, although HCT is defined on a macro-element.

While $C^1$ elements are difficult to implement, this paper, showed that they're substantially better for biharmonic problems than the more common $C^0$ interior penalty method. So it's possible that surmounting the implementation challenges is really worth it, especially if you can build off of someone else's work.