Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
I am trying to learn about 2D FEM methods. I am trying to understand the generalization of Lagrange polynomial basis from 1D into 2 variable polynomials over triangle domains.
The most basic element if I understand correctly is piece Wise linearity, obtained through barycentric coordinates.
What would the quadratic element be? Or am I mistaken about the linear element?
You can rather easily write Lagrange bases over any dimension simplices. Define a $k$-simplex as the convex combinations of $k+1$ points, e.g. a triangle is a $2$-simplex.
As a triangle $K$ of vertices $P_1$, $P_2$, $P_3$ is the convex combination of these three, it is convenient to talk about the set of barycentric coordinates:
$\hat K = \{\xi \in [0,1]^{k+1}\ / \ |\xi|_1 = 1 \}$
Then $K$ is mapped to from $\hat K$ by a mapping $F_K : \xi \in \hat K \mapsto \xi_1 P_1 + \xi_2 P_2 + \xi_3 P_3$. You can also define non-linear simplices this way.
Next, we define interpolation nodes. These are arbitrary, which is not to say their choice is without consequences. The most basic set of interpolation nodes is the set of regularly spaced nodes. For this, and to index the Lagrange basis functions, we define the set of multi-indices for degree $d$ $k$-simplices:
$\hat K ^d_k = \{ \alpha \in [\![0,d]\!]\ / \ |\alpha|_1 = d \}$
For instance, for triangles,
$\hat K^d_2 = \{ 0 \leq i,j,k \leq d, i+j+k = d\}$
These multi-indices are in the same number as there are elements in a basis of degree $d$ polynomials of $k$ variables or homogeneous degree $d+1$ polynomials of $k+1$ variables. You can see them as the monomials
$\alpha \to \xi_1^{\alpha_1}\xi_2^{\alpha_2} ... \xi_{k+1}^{\alpha_{k+1}}$
which are none other than the Bernstein polynomials up to scaling. The former represents polynomials over the reference element (which I haven't talked about), whereas the latter represents polynomials over the barycentric domain $\hat K$. You go from $\xi \in \hat K$ to $x$ in the reference element by, for instance, setting $x = \xi_2, y = \xi_3$; reciprocically, by setting $\xi = (1-x-y,x,y)$. This may be counter-intuitive but it corresponds to the most widely used reference element with $\hat P_1 = 0$, $\hat P_2 = (1,0)$ and $\hat P_3 = (0,1)$.
Using these multi-indices, we define the regularly spaced nodes as the barycentric coordinates $\{\xi^\alpha\}_{\alpha \in \hat K^d_k}$ such that, for all $\alpha \in K^d_k$,
$\xi^\alpha = \alpha / d$.
Next, we define the Lagrange basis as the set of polynomials $\{\phi_\alpha\}_{\alpha \in \hat K^d_k}$ such that, for all $\alpha,\beta \in K^d_k$,
$\phi_{\alpha}(\xi^\beta) = \delta_{\alpha\beta}$
where $\delta$ is the Kronecker symbol, equal to 0, or 1 iff $\alpha = \beta$.
You can take the above formula and power through it. If writing in a basis of monomials similar to the above, that involves solving a $(d+1)(d+2)/2$ size system for triangles per polynomial, or size $(d+1)(d+2)(d+3)/6$ for tetrahedra. It quickly becomes unmanageable!
I have found a rather nice geometric construction method, which I believe is equivalent to existing formulas. It has the merit to be very easy to reconstruct with bad memory.
First, we draw the regularly spaced nodes, as follows:
Recall these nodes are ordered by a multi-index. The node $300$ is at the first vertex ($\xi = (3,0,0) / 3$), the node $210$ is $2/3$ of the first, $1/3$ of the second, and so on. Now we want a Lagrange function, say $\phi_{300}$. We know it's a homogeneous degree $3$ function of the $k+1$ barycentric coordinates. Therefore, we're going to write it as a product of $k+1$ degree $1$ monomials, of the form $(\xi_i - c_i)$. We also know it has to vanish at all other nodes than the node $\xi^{300}$. Therefore, we're going to choose the monomials to offer roots at those other nodes. These are represented by segments crossing out nodes on the element. The "puzzle" reduces thus to choosing three segments such that all nodes but the one marked $u$ are crossed out:
Those lines have the equations $\{ u = 0 \}$, $\{ u = 1/3 \}$ and $\{ u = 2/3 \}$. As such, our polynomial is
$\phi_{300}(u,v,w) = C(u-0)(u-1/3)(u-2/3)$
which guarantees that it vanishes on all those nodes crossed out. However, we still need to fit the constant to have $\phi_{300}(1,0,0) = 1$. That's simply a matter of evaluating at the node; we get $C = 9/2$. Then, it's clear to see that:
$\phi_{210}(u,v,w) = C(u-0)(v-0)(u-1/3)$
and
$\phi_{111}(u,v,w) = C(u-0)(v-0)(w-0)$
The others are obtained by permutations. If you want those in physical coordinates, you simply replace $(u,v,w)$ by $(1-x-y,x,y)$. If you ever have to derive this basis for degree $5$, $6$, etc... you'll be glad not to have to invert a size $21$, $28$... system for each polynomial!
As a final note, I work here in barycentric coordinates because that is how these bases are defined in Finite Elements. When working on non-linear elements, you do not know the bases explicitly in the physical frame, as those are composed by $F_K^{-1}$ (the inverse of a two-variate polynomial or another mapping). The fact these are taught on linear elements most of the time can make this fact slightly confusing.
Carlos Felippa, in Chapter 24 of the class notes for his introductory FEM class, has an excellent discussion of the formulation of the quadratic, isoparametric triangle element. These notes do not seem to be available on his web site any longer but you can find copies various places on the web, Chapter 24, for example.
To answer your specific question, the shape functions for the quadratic triangle in terms of the three barycentric coordinates (equation 24.12 in his notes) are
$$\begin{equation} {\bf N} = \left\{ \begin{array}{l} \zeta_1(2\zeta_1 - 1) \\ \zeta_2(2\zeta_2 - 1) \\ \zeta_3(2\zeta_3 - 1) \\ 4\zeta_1\zeta_2 \\ 4\zeta_2\zeta_3 \\ 4\zeta_3\zeta_1 \end{array} \right\} \end{equation} $$
