Finite elements in 2D and 3D

Finite element approximation is particularly powerful in 2D and 3D because the method can handle a geometrically complex domain $ \Omega $ with ease. The principal idea is, as in 1D, to divide the domain into cells and use polynomials for approximating a function over a cell. Two popular cell shapes are triangles and the quadrilaterals. Figures 34, 35, and 36 provide examples. P1 elements means linear functions ($ a_0 + a_1x + a_2y $) over triangles, while Q1 elements have bilinear functions ($ a_0 + a_1x + a_2y + a_3xy $) over rectangular cells. Higher-order elements can easily be defined.

Figure 34: Examples on 2D P1 elements.

Figure 35: Examples on 2D P1 elements in a deformed geometry.

Figure 36: Examples on 2D Q1 elements.

Basis functions over triangles in the physical domain

Cells with triangular shape will be in main focus here. With the P1 triangular element, $ u $ is a linear function over each cell, as depicted in Figure 37, with discontinuous derivatives at the cell boundaries.

Figure 37: Example on piecewise linear 2D functions defined on triangles.

We give the vertices of the cells global and local numbers as in 1D. The degrees of freedom in the P1 element are the function values at a set of nodes, which are the three vertices. The basis function $ \basphi_i(x,y) $ is then 1 at the vertex with global vertex number $ i $ and zero at all other vertices. On an element, the three degrees of freedom uniquely determine the linear basis functions in that element, as usual. The global $ \basphi_i(x,y) $ function is then a combination of the linear functions (planar surfaces) over all the neighboring cells that have vertex number $ i $ in common. Figure 38 tries to illustrate the shape of such a "pyramid"-like function.

Figure 38: Example on a piecewise linear 2D basis function over a patch of triangles.

Element matrices and vectors

As in 1D, we split the integral over $ \Omega $ into a sum of integrals over cells. Also as in 1D, $ \basphi_i $ overlaps $ \basphi_j $ (i.e., $ \basphi_i\basphi_j\neq 0 $) if and only if $ i $ and $ j $ are vertices in the same cell. Therefore, the integral of $ \basphi_i\basphi_j $ over an element is nonzero only when $ i $ and $ j $ run over the vertex numbers in the element. These nonzero contributions to the coefficient matrix are, as in 1D, collected in an element matrix. The size of the element matrix becomes $ 3\times 3 $ since there are three degrees of freedom that $ i $ and $ j $ run over. Again, as in 1D, we number the local vertices in a cell, starting at 0, and add the entries in the element matrix into the global system matrix, exactly as in 1D. All details and code appear below.

Basis functions over triangles in the reference cell

As in 1D, we can define the basis functions and the degrees of freedom in a reference cell and then use a mapping from the reference coordinate system to the physical coordinate system. We also have a mapping of local degrees of freedom numbers to global degrees of freedom numbers.

The reference cell in an $ (X,Y) $ coordinate system has vertices $ (0,0) $, $ (1,0) $, and $ (0,1) $, corresponding to local vertex numbers 0, 1, and 2, respectively. The P1 element has linear functions $ \refphi_r(X,Y) $ as basis functions, $ r=0,1,2 $. Since a linear function $ \refphi_r(X,Y) $ in 2D is on the form $ C_{r,0} + C_{r,1}X + C_{r,2}Y $, and hence has three parameters $ C_{r,0} $, $ C_{r,1} $, and $ C_{r,2} $, we need three degrees of freedom. These are in general taken as the function values at a set of nodes. For the P1 element the set of nodes is the three vertices. Figure 39 displays the geometry of the element and the location of the nodes.

Figure 39: 2D P1 element.

Requiring $ \refphi_r=1 $ at node number $ r $ and $ \refphi_r=0 $ at the two other nodes, gives three linear equations to determine $ C_{r,0} $, $ C_{r,1} $, and $ C_{r,2} $. The result is $$ \begin{align} \refphi_0(X,Y) &= 1 - X - Y,\\ \refphi_1(X,Y) &= X,\\ \refphi_2(X,Y) &= Y \end{align} $$

Higher-order approximations are obtained by increasing the polynomial order, adding additional nodes, and letting the degrees of freedom be function values at the nodes. Figure 40 shows the location of the six nodes in the P2 element.

Figure 40: 2D P2 element.

A polynomial of degree $ p $ in $ X $ and $ Y $ has $ n_p=(p+1)(p+2)/2 $ terms and hence needs $ n_p $ nodes. The values at the nodes constitute $ n_p $ degrees of freedom. The location of the nodes for $ \refphi_r $ up to degree 6 is displayed in Figure 41.

Figure 41: 2D P1, P2, P3, P4, P5, and P6 elements.

The generalization to 3D is straightforward: the reference element is a tetrahedron with vertices $ (0,0,0) $, $ (1,0,0) $, $ (0,1,0) $, and $ (0,0,1) $ in a $ X,Y,Z $ reference coordinate system. The P1 element has its degrees of freedom as four nodes, which are the four vertices, see Figure 42. The P2 element adds additional nodes along the edges of the cell, yielding a total of 10 nodes and degrees of freedom, see Figure 43.

Figure 42: P1 elements in 1D, 2D, and 3D.

Figure 43: P2 elements in 1D, 2D, and 3D.

The interval in 1D, the triangle in 2D, the tetrahedron in 3D, and its generalizations to higher space dimensions are known as simplex cells (the geometry) or simplex elements (the geometry, basis functions, degrees of freedom, etc.). The plural forms simplices and simplexes are also a much used shorter terms when referring to this type of cells or elements. The side of a simplex is called a face, while the tetrahedron also has edges.

Acknowledgment. Figures 39 to 43 are created by Anders Logg and taken from the FEniCS book: Automated Solution of Differential Equations by the Finite Element Method, edited by A. Logg, K.-A. Mardal, and G. N. Wells, published by Springer, 2012.

Affine mapping of the reference cell

Let $ \refphi_r^{(1)} $ denote the basis functions associated with the P1 element in 1D, 2D, or 3D, and let $ \xdno{q(e,r)} $ be the physical coordinates of local vertex number $ r $ in cell $ e $. Furthermore, let $ \X $ be a point in the reference coordinate system corresponding to the point $ \x $ in the physical coordinate system. The affine mapping of any $ \X $ onto $ \x $ is then defined by $$ \begin{equation} \x = \sum_{r} \refphi_r^{(1)}(\X)\xdno{q(e,r)}, \tag{44} \end{equation} $$ where $ r $ runs over the local vertex numbers in the cell. The affine mapping essentially stretches, translates, and rotates the triangle. Straight or planar faces of the reference cell are therefore mapped onto straight or planar faces in the physical coordinate system. The mapping can be used for both P1 and higher-order elements, but note that the mapping itself always applies the P1 basis functions.

Figure 44: Affine mapping of a P1 element.

Isoparametric mapping of the reference cell

Instead of using the P1 basis functions in the mapping (44), we may use the basis functions of the actual Pd element: $$ \begin{equation} \x = \sum_{r} \refphi_r(\X)\xdno{q(e,r)}, \tag{45} \end{equation} $$ where $ r $ runs over all nodes, i.e., all points associated with the degrees of freedom. This is called an isoparametric mapping. For P1 elements it is identical to the affine mapping (44), but for higher-order elements the mapping of the straight or planar faces of the reference cell will result in a curved face in the physical coordinate system. For example, when we use the basis functions of the triangular P2 element in 2D in (45), the straight faces of the reference triangle are mapped onto curved faces of parabolic shape in the physical coordinate system, see Figure 45.

Figure 45: Isoparametric mapping of a P2 element.

From (44) or (45) it is easy to realize that the vertices are correctly mapped. Consider a vertex with local number $ s $. Then $ \refphi_s=1 $ at this vertex and zero at the others. This means that only one term in the sum is nonzero and $ \x=\xdno{q(e,s)} $, which is the coordinate of this vertex in the global coordinate system.

Computing integrals

Let $ \tilde\Omega^r $ denote the reference cell and $ \Omega^{(e)} $ the cell in the physical coordinate system. The transformation of the integral from the physical to the reference coordinate system reads $$ \begin{align} \int_{\Omega^{(e)}}\basphi_i (\x) \basphi_j (\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) \refphi_j (\X) \det J\, \dX,\\ \int_{\Omega^{(e)}}\basphi_i (\x) f(\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) f(\x(\X)) \det J\, \dX, \end{align} $$ where $ \dx $ means the infinitesimal area element $ dx dy $ in 2D and $ dx dy dz $ in 3D, with a similar definition of $ \dX $. The quantity $ \det J $ is the determinant of the Jacobian of the mapping $ \x(\X) $. In 2D, $$ \begin{equation} J = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial x}{\partial Y}\\ \frac{\partial y}{\partial X} & \frac{\partial y}{\partial Y} \end{array}\right], \quad \det J = \frac{\partial x}{\partial X}\frac{\partial y}{\partial Y} - \frac{\partial x}{\partial Y}\frac{\partial y}{\partial X} \tp \tag{46} \end{equation} $$ With the affine mapping (44), $ \det J=2\Delta $, where $ \Delta $ is the area or volume of the cell in the physical coordinate system.

Remark. Observe that finite elements in 2D and 3D builds on the same ideas and concepts as in 1D, but there is simply much more to compute because the specific mathematical formulas in 2D and 3D are more complicated and the book keeping with dof maps also gets more complicated. The manual work is tedious, lengthy, and error-prone so automation by the computer is a must.