Basis function in a tetraedron for finite elements contex

In the finite element method we need to know a base for the fem spaces. For example, a base for the space

$$P_1(\hat{K})=<\{1-x-y,x,y,z\}>$$

is a typical base for the polynomials of degree less than or equal to 1 on a tetrahedron with vertex (0,0,0), (1,0,0), (0,1,0) and (0,0,1).

But, which is the base for $$P_2(\hat{K})$$ and $$P_3(\hat{K})$$? I know that the dimension is 10 and 20, respectively, but I do not can see how to derive it.

• Note that for the P1 element, the first basis function also includes $z$. The correct function is $1-x-y-z$. Aug 21 at 9:02

You can follow the concept of Pascal's pyramid (see the picture) for identifying the basis functions for elements of different orders.

The picture is taken from what-when-how.com. Refer to the linked site for the details on computing the coefficients.

I derived P2 before. P3 should be similar.

You already know the dimension, which is good. We also need to know the position of nodes. For P2, we have 4 at vertices, and 6 at edges. Then the rest is to determine the coefficients of an interpolating polynomial.

Here is the general idea. The 10 monomials $$m_i$$ are 1, x, y, z, x*y, ... etc. For each term we want to determine its coefficient $$c_i$$, So that interpolating polynomial $$q = \sum_i m_i c_i$$ evaluates exact node position given x, y, and z. Solving the system will give us $$c_i$$, and substitute them back to $$q$$, regrouping the terms will produce the basis in terms of x,y,and z.

Here is some maple code to help with the derivation:

restart;
with(LinearAlgebra):
# the 10 monomials
m:=[1,x,y,z,x*y,y*z,z*x,x^2,y^2,z^2];
# We are going to determine ci
q:=+(seq(c||i * m[i],i=1..10));

# generate node positions
node:=[<0,0,0>,<1,0,0>,<0,1,0>,<0,0,1>];
edge:=[seq(seq((node[i]+node[j])/2,j=i+1..4),i=1..4)];
all:=[op(node),op(edge)];

# Evaluate at node coordinates, and require them to be equal to node positions qi.
E:=[seq(subs(x=all[i][1],y=all[i][2],z=all[i][3],q=q||(i-1)),i=1..10)];
# Solve for coefficients ci, and substitute back to q
q:=subs(eliminate(E,[seq(c||i,i=1..10)])[1],q);
# Find the coefficient of each qi (the node position)
N:=[seq(diff(q,q||i),i=0..9)];


Then the P2 basis from above code:

> lprint(N);
[2*x^2+4*x*y+4*x*z+2*y^2+4*y*z+2*z^2-3*x-3*y-3*z+1, 2*x^2-x, 2*y^2-y, 2*z^2-z, -4*x^2-4*x*y-4*x*z+4*x, -4*x
*y-4*y^2-4*y*z+4*y, -4*x*z-4*y*z-4*z^2+4*z, 4*x*y, 4*z*x, 4*y*z]