3.9. Finite Element Approximation

In Discretization, we discretized the governing equations in space using Galerkin’s method. In this chapter, we discuss two important, unresolved issues. The first issue is the method for constructing the basis functions $\symWeightFunc$ , which appear in the weak statements given by the stationary weak description and the transient weak description. The second issue is the calculation of the associated integrals over the domain and its boundary. The finite element method solves both of these issues, and is concisely described by Becker, et al [22]:

The finite element method provides a general and systematic technique for constructing basis functions for Galerkin approximations of boundary value problems. The main idea is that the basis functions […] can be defined piecewise over subregions of the domain called finite elements and that over any subdomain the [basis functions] can be chosen to be very simple functions such as polynomials of low degree.

If the finite elements are also chosen to have convenient shapes (like triangles and quadrilaterals) for calculating integrals, then the integration over the entire domain is also facilitated.

To elaborate, consider Fig. 3.10, which shows a one-dimensional domain partitioned into four finite elements. In Fig. 3.11, the global, piecewise linear basis function for node i, $\symWeightFunc_i$ , is shown as the standard hat function. This function has a value of unity at node i and zero in elements not connected to node i. The finite element approach to constructing $\symWeightFunc_i$ is shown in Fig. 3.11. On each element, two linear functions $\symWeightFunc^e_i$ are defined, $e\in\{1,2,3,4\}$ , $i\in\{1,2\}$ . If we add $\symWeightFunc^2_2$ from the left element to ${\symWeightFunc}^3_1$ from the right element, then we see that the global basis function $\symWeightFunc$ is constructed.

one-dimensional basis function — Fig. 3.10 A linear finite element global basis function at a node is formed from the nodal shape functions of two neighboring elements. The global basis function for node $i$ has support over the two elements $e$ and $e+1$ that share the node, and takes on the value of zero everywhere else.

Fig. 3.10 A linear finite element global basis function at a node is formed from the nodal shape functions of two neighboring elements. The global basis function for node $i$ has support over the two elements $e$ and $e+1$ that share the node, and takes on the value of zero everywhere else.

one-dimensional shape function — Fig. 3.11 The associated element shape functions, $\symWeightFunc^e_2$ (on local node 2 of element $e$ ), and $\symWeightFunc^{e+1}_1$ (on local node 1 of element $e+1$ ), are combined to define the global basis function for node $i$ .

Fig. 3.11 The associated element shape functions, $\symWeightFunc^e_2$ (on local node 2 of element $e$ ), and $\symWeightFunc^{e+1}_1$ (on local node 1 of element $e+1$ ), are combined to define the global basis function for node $i$ .

Locally, the element’s shape functions can be used to interpolate values defined at the nodes to arbitrary location within the element. Specifically, for some quantity $q$ defined at the $N$ (local) nodes of the element $e$ , the value of $q$ within the element is

$q(\vector{x},t) = \sum_{i = 1}{N} q_i(t) \symWeightFunc^e_i(\vector{x})$

Now, consider the global stiffness matrix $\bfK$ , with row $i$ and column $j$ , defined in the previous section. We now rewrite this integral as the sum over each finite element, namely

(3.168) $K_{ij} = \sum\limits_{e=1}^E \int_{\Omega^e} {\bf \nabla} \symWeightFunc_i \cdot {\bf k} {\bf \nabla} \symWeightFunc_j \,d\Omega$

In (3.168), the global basis function $\symWeightFunc_i$ is nonzero only on those elements which support node $i$ . Hence we can consider the contributions to the global stiffness matrix from element $e$ , and define

(3.169) $K^e_{ij} = \int_{\Omega^e} {\bf \nabla} {\symWeightFunc}^e_i \cdot {\bf k} {\bf \nabla} {\symWeightFunc}^e_j \,d\Omega$

There is a subtle difference in the interpretation of the indices in (3.168) as opposed to (3.169). In (3.168), the indices $i$ and $j$ are global, and span all nodes in the finite element mesh. In (3.169), the indices $i$ and $j$ are local, and only span the number of element shape functions. Aria calculates the element stiffness matrix contributions defined in (3.169) and assembles the results into the correct locations in the global stiffness matrix using a mapping from the local node numbers of a given element to the associated global node numbers. In fact, each of the integrals in the weak statements derived in Discretization involves a basis function and is calculated in the same way.

In the following section, Element Integration, we discuss how the integrals are computed on a master element. Next, in Linear Master Elements, we discuss the master element interpolation functions for the elements used in Aria.

3.9.1. Element Integration

Element integrals of the form given by (3.169) are defined on the element coordinates in physical space. When calculating these quantities numerically, it is convenient to transform these integrals from physical space to a reference element, also called a master element, in parametric space. This transformation is performed locally, for each finite element in the mesh. In this section we discuss the details of this process for two-dimensional physical space in Cartesian coordinates. The approach is completely general, and extensible to three-dimensional space and other coordinate systems. This situation is illustrated in Fig. 3.12, showing such a transformation from the quadrilateral element, which spans the region $[-1,1]\times[-1,1]$ , in parametric coordinates $(\xi, \eta)$ to an arbitrary region in physical coordinates $(x,y)$ .

quad mapping in master element — Fig. 3.12 Mapping from the master element to physical coordinates: The quadrilateral master element in $(\xi,\eta)$ on the region $[-1,1]\times[-1,1]$ is mapped to its image in physical coordinates $(x,y)$ .

Fig. 3.12 Mapping from the master element to physical coordinates: The quadrilateral master element in $(\xi,\eta)$ on the region $[-1,1]\times[-1,1]$ is mapped to its image in physical coordinates $(x,y)$ .

Returning to (3.169), given some function $f(x,y)$ , we wish to calculate the integral

(3.170) $I = \int_{\Omega^e} f(x,y) \, dx \, dy$

on the master element in $(\xi, \eta)$ space. The change of variables theorem [28] from vector calculus states that (3.170) is equivalent to

(3.171) $I = \int_{\hat{\Omega}^e} \hat{f}(\xi, \eta) |{\bf J}| \, d\xi \, d\eta$

where $\hat{\Omega}^e$ is the master element and $|{\bf J}|$ is the determinant of the non-singular Jacobian matrix of the transformation

(3.172) $\left ( \begin{array}{c} x \\ y \end{array} \right ) = \left ( \begin{array}{c} x(\xi,\eta) \\ y(\xi, \eta) \end{array}\right)$

Recall the Jacobian matrix of (3.172)

(3.173) $\mathbf{J} = % \begin{pmatrix} \dfrac{\partial x}{\partial \xi} & \dfrac{\partial x}{\partial \eta}\\[12pt] % \dfrac{\partial y}{\partial \xi} & \dfrac{\partial y}{\partial \eta} % \end{pmatrix}$

and its determinant

(3.174) $|{\bf J}| = \frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \eta} - \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \xi}$

Using the local finite element shape functions to construct isoparametric transformations for (3.172), we have

(3.175) $x(\xi,\eta) & = & \sum \limits_{i=1}^{n} x_i \hat{\symWeightFunc}^e_i \\ y(\xi,\eta) & = & \sum \limits_{i=1}^{n} y_i \hat{\symWeightFunc}^e_i ,$

where $n$ is the number of shape functions on element $e$ , and the coordinates of each node of the element is given by the pair $(x_i, y_i)$ .

Since we seek a numerical solution, it is convenient to compute the integrals numerically; we use Gaussian quadrature rules of the form

(3.176) $I = \int_{\hat{\Omega}^e} \hat{f}(\xi, \eta) |{\bf J}| \, d\xi \, d\eta \approx \sum\limits_{g=1}^{G} w_g\hat{f}(\xi_g, \eta_g) |{\bf J}| _{(\xi_g, \eta_g)},$

where the number of Gauss points is denoted $G$ , $(\xi_g,\eta_g)$ are the Gauss point coordinates (abscissas), $|{\bf J}|_ {(\xi_g, \eta_g)}$ indicates that the Jacobian determinant is to be evaluated at the Gauss point coordinates, and $w_g$ is the associated weight. For master element shapes in common use, such as triangles, hexahedra, quadrilaterals, etc., the quadrature points and weights are well-known and tabulated, e.g., see [22].

Table 3.1 Weights and coordinates for the Gauss quadrature rules on elements with linear shape functions.
element topology	degree $G$	weights $w_g$	coordinates $\xi_g$
1D Edge	$2$	$1, 1$	$-\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}}$
2D Triangle	$3$	$\frac{1}{6}, \frac{1}{6}, \frac{1}{6}$	$(\frac{1}{2}, \frac{1}{2}), (0, \frac{1}{2}), (\frac{1}{2}, 0)$
2D Quadrilateral	$4$	$\frac{1}{4}, \frac{1}{4}, \frac{1}{4}, \frac{1}{4}$	$(-\gamma, -\gamma), ( \gamma, -\gamma),$ $(\gamma, \gamma), (-\gamma, \gamma)$
3D Tetrahedron	$4$	$\frac{1}{24}, \frac{1}{24}, \frac{1}{24}, \frac{1}{24}$	$(\alpha, \beta, \beta), ( \beta, \alpha, \beta),$ $( \beta, \beta, \alpha), ( \beta, \beta, \beta)$
3D Hexahedron	$8$	$\frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}, \frac{1}{8}$	$(-\gamma,-\gamma,-\gamma), ( \gamma,-\gamma,-\gamma),$ $( \gamma, \gamma,-\gamma), (-\gamma, \gamma,-\gamma),$ $(-\gamma,-\gamma, \gamma), ( \gamma,-\gamma, \gamma),$ $( \gamma, \gamma, \gamma), (-\gamma, \gamma, \gamma)$

where $\gamma=\sqrt{3}/6$ , $\alpha=0.58541020$ and $\beta=0.13819660$ . For linear elements, Aria uses the rules shown in Table 3.1, which correspond to the master element topologies given in Linear Master Elements. Because the order of quadrature is fixed for a given element (and the associated polynomial degree of element shape functions), the user of Aria should not forget to account for possible errors in the quadrature. Error in quadrature could be significant when user subroutines for boundary conditions, source terms, etc., attempt to represent highly oscillatory, or unsmooth data.

3.9.2. Linear Master Elements

In this section, we present the shape functions $\hat{\psi}^e_i$ for the linear master elements currently supported in Aria. Presently, all elements use the Lagrange basis functions. For the sake of brevity, we drop the superscript $e$ on the element shape function. Hence, the notation $\hat{\psi}_i$ is equivalent to $\hat{\psi}^e_i$ .

3.9.2.1. 2D Triangular Element

The master element for the two-dimensional triangle with nodes numbered locally as shown. The Gauss quadrature rule given in Table 3.1 for the 2D Triangle is used to integrate this element. The shape functions are

(3.177) $\hat{\psi}_1 & = & {1 - \xi - \eta } \hat{\psi}_2 & = & \xi \hat{\psi}_3 & = & \eta$

3.9.2.2. 2D Quadrilateral Element

The master element for the two-dimensional quadrilateral is shown in Fig. 3.12. The Gauss quadrature rule for the 2D Quadrilateral given in Table 3.1 is used to integrate this element. The shape functions are

(3.178) $\hat{\psi}_1 & = & \frac{1}{4}(1 - \xi )(1 - \eta ) \hat{\psi}_2 & = & \frac{1}{4}(1 + \xi )(1 - \eta ) \hat{\psi}_3 & = & \frac{1}{4}(1 + \xi )(1 + \eta ) \hat{\psi}_4 & = & \frac{1}{4}(1 - \xi )(1 + \eta )$

3.9.2.3. 2D Face Element

When applying flux boundary conditions, it is often necessary to integrate on the surfaces of the two-dimensional elements presented in 2D Triangular Element and 2D Quadrilateral Element. In both cases, the surface is topologically a one-dimensional edge. The two shape functions that define this element are

(3.179) $\hat{\psi}_1 & = & \frac{1}{2}(1 - \xi ) \hat{\psi}_2 & = & \frac{1}{2}(1 + \xi )$

The Gauss quadrature rule for the 1D Edge topology that is given in Table 3.1 is used to integrate this element.

Note that the location along this edge is given by the single parametric coordinate, $\xi$ , although there is a pair of physical Cartesian coordinates, $(x, y)$ associated with each value of $\xi$ . This fact makes the construction of the transformation given as (3.172) a little less obvious. In this case, the determinant of the Jacobian matrix is the ratio of the edge length in physical coordinates to the edge length in computational coordinates,

$|{\bf J}| = \sqrt{ \left(\frac{\partial x}{\partial \xi}\right)^2 + \left(\frac{\partial y}{\partial \xi}\right)^2}$

Since we are using isoparametric mappings, it is easy to show that, using the linear shape functions given as (3.179), this determinant is given by

$|{\bf J}| = \frac{1}{2}\sqrt{ \left(x_2 - x_1\right)^2 + \left(y_2 - y_1\right)^2}$

3.9.2.4. 3D Tetrahedron Element

The topology of the four-node, linear tetrahedron element is shown in Fig. 3.13.

Fig. 3.13 The master tetrahedron element is mapped to its image in physical coordinates $(x,y,z)$ .

The Gauss quadrature rule for the 3D Tetrahedron given in Table 3.1 is used to integrate this element. The shape functions are

(3.180) $\hat{\psi}_1 & = & {1 - \xi - \eta -\zeta } \hat{\psi}_2 & = & \xi \hat{\psi}_3 & = & \eta \hat{\psi}_4 & = & \zeta$

3.9.2.5. 3D Triangular Face Element

The tetrahedral element presented in 3D Tetrahedron Element has three-node, triangular faces, which are topologically equivalent to the two-dimensional triangular element that was presented in 2D Triangular Element. Accordingly, the shape functions are identical (3.177). The Gauss quadrature rule for the 2D Triangle that is given in Table 3.1 is used to integrate this element.

Since the physical coordinates are three-dimensional, but there are only two parameters on the surface, $(\xi, \eta)$ , the transformation from computational coordinates to physical coordinates requires a little elaboration. Vector calculus provides a formula for the integral of a scalar function over a parameterized surface [28]. We are interested in the surface integral

(3.181) $I = \int_{\Gamma^e} f(x,y,z) \, dx\, dy\, dz$

where $\Gamma^e$ is the three dimensional surface of our face element, and is given by the isoparametric mapping

(3.182) $\left ( \begin{array}{c} x(\xi,\eta) \\ y(\xi,\eta) \\ z(\xi,\eta) \end{array} \right ) = \left ( \begin{array}{c} \sum \limits_{i=1}^3 x_i\hat{\psi}_i(\xi,\eta)\\ \sum \limits_{i=1}^3 y_i\hat{\psi}_i(\xi,\eta)\\ \sum \limits_{i=1}^3 z_i\hat{\psi}_i(\xi,\eta) \end{array}\right)$

It can be shown that (3.181) may be computed as

(3.183) $I = \int_{\hat{\Gamma}^e} \hat{f}(\xi, \eta) |{\bf J}| \, d\xi \, d\eta,$

where $\hat{\Gamma}^e$ is the master element, $\hat{f}(\xi, \eta) = f(x(\xi,\eta), y(\xi, \eta), z(\xi, \eta))$ , and the determinant of the Jacobian matrix may be written as

(3.184) $|{\bf J}| = \sqrt{ \left| \frac{\partial (x,y)}{\partial (\xi, \eta)}\right|^2 + \left| \frac{\partial (y,z)}{\partial (\xi, \eta)}\right|^2 + \left| \frac{\partial (x,z)}{\partial (\xi, \eta)}\right|^2 }$

Equation (3.184) involves the calculation of the determinants of three sub-matrices. For example,

$\left| \frac{\partial (x,y)}{\partial (\xi, \eta)}\right| = \frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \eta} - \frac{\partial x}{\partial \eta} \frac{\partial y}{\partial \xi}$

3.9.2.6. 3D Hexahedron Element

Fig. 3.14 shows the topology of the eight-node, linear hexahedral master element.

Fig. 3.14 The master hexahedron element in $(\xi,\eta,\zeta)$ on the region $[-1,1]\times[-1,1]\times[-1,1]$ is mapped to its image in physical coordinates $(x,y,z)$ .

The Gauss quadrature rule for the 3D Hexahedron that is given in Table 3.1 is used to integrate this element. The shape functions are

$\hat{\psi}_1 & = & {\frac{1}{8}(1 - \xi )(1 - \eta )(1 - \zeta)} \hat{\psi}_2 & = & {\frac{1}{8}(1 + \xi )(1 - \eta )(1 - \zeta)} \hat{\psi}_3 & = & {\frac{1}{8}(1 + \xi )(1 + \eta )(1 - \zeta)} \hat{\psi}_4 & = & {\frac{1}{8}(1 - \xi )(1 + \eta )(1 - \zeta)} \hat{\psi}_5 & = & {\frac{1}{8}(1 - \xi )(1 - \eta )(1 + \zeta)} \hat{\psi}_6 & = & {\frac{1}{8}(1 + \xi )(1 - \eta )(1 + \zeta)} \hat{\psi}_7 & = & {\frac{1}{8}(1 + \xi )(1 + \eta )(1 + \zeta)} \hat{\psi}_8 & = & {\frac{1}{8}(1 - \xi )(1 + \eta )(1 + \zeta)}$

3.9.2.7. 3D Quadrilateral Face Element

The hexahedron element described in 3D Hexahedron Element has quadrilateral faces. The shape functions that are used for interpolating on one of these faces are identical to the shape functions used for the 2D Quadrilateral Element in (3.178). As discussed in 3D Triangular Face Element, integrals on this surface master element may be computed according to (3.183), with the determinant given by (3.184). The Gauss quadrature rule for the 2D Quadrilateral that is given in Table 3.1 is used to integrate this element.

3.9.2.8. 3D Triangular Shell Element

Beta Capability

Shell support in Aria is currently a beta feature. This capability is a work-in-progress and is not thoroughly tested and should be used with caution.

Sometimes, a subdomain $\Omega_i\subset \Omega$ may be very thin, and may also possibly consist of a highly conductive material like aluminum. In this case, it is reasonable to neglect the thermal gradient through the thickness of $\Omega_i$ . Shell elements are specialized volume elements that are often used to model this situation. The topology of the linear triangular shell is illustrated in Fig. 3.15, which has three nodes, and two faces. Topologically, the shell has no thickness that can be derived solely from its node coordinates; instead, the thickness is an attribute specified as input data. The shape functions are identical to (3.177) in 2D Triangular Element.

The transformation from computational coordinates to physical coordinates for the three-dimensional shell element depends on the thickness attribute. In this case, we wish to calculate the volume integral

(3.185) $I = \int_{\Omega^e} f(x,y,z) \, dx \, dy \, dz$

The parametric mapping from computational coordinates to physical coordinates may be regarded as the surface mapping, plus a component that comes up out of the shell mid-plane. With reference to Fig. 3.15, we may write this mapping as

(3.186) $\left ( \begin{array}{c} x(\xi,\eta, \zeta) \\ y(\xi,\eta, \zeta) \\ z(\xi,\eta, \zeta) \end{array} \right ) = \left ( \begin{array}{c} \sum \limits_{i=1}^I x_i\hat{\psi}_i(\xi,\eta) + \zeta \theta_x \sum \limits_{i=1}^I \tau_i\hat{\psi}_i(\xi,\eta) \\ \sum \limits_{i=1}^I y_i\hat{\psi}_i(\xi,\eta) + \zeta \theta_y \sum \limits_{i=1}^I \tau_i\hat{\psi}_i(\xi,\eta)\\ \sum \limits_{i=1}^I z_i\hat{\psi}_i(\xi,\eta) + \zeta \theta_z \sum \limits_{i=1}^I \tau_i\hat{\psi}_i(\xi,\eta) \end{array}\right)$

where the shape functions $\hat{\psi}_i$ are given by (3.177), $I=3$ , $\tau_i$ denotes the thickness of the shell at node $i$ , $\zeta \in [-1,1]$ is the parametric coordinate in the direction locally orthogonal to the shell mid-plane, and we have defined the inner products

(3.187) $\begin{array}{lll} \theta_x & =& \hat{\bf n}_{\zeta} \cdot \hat{\bf i} \\ \theta_y & = &\hat{\bf n}_{\zeta} \cdot \hat{\bf j}\\ \theta_z & = &\hat{\bf n}_{\zeta} \cdot \hat{\bf k} \end{array}$

Triangular shell element — Fig. 3.15 The midplane of the master triangular shell element is mapped to its image in physical coordinates $(x,y,z)$ , and the variation in the shell thickness is represented with the element shape functions.

Fig. 3.15 The midplane of the master triangular shell element is mapped to its image in physical coordinates $(x,y,z)$ , and the variation in the shell thickness is represented with the element shape functions.

In (3.187), $\hat{\bf n}_{\zeta}$ is the unit vector in the direction locally orthogonal to the element mid-plane, and $\hat{\bf i}, \hat{\bf j}$ and $\hat{\bf k}$ are the unit vectors in the $x,y$ and $z$ coordinate directions, respectively. The volume integral in (3.185) on the master element is

(3.188) $I = \int_{\hat{\Omega}^e} \hat{f}(\xi, \eta, \zeta) | {\bf J} | \; d\xi \, d\eta \, d\zeta$

where $|{\bf J}|$ is the determinant of the Jacobian matrix obtained from the transformation given by (3.186).

Aria currently only supports constant thickness shells. In this case, $I$ reduces to

(3.189) $I = \tau \int_{\hat{\Gamma}^e} \hat{f}(\xi, \eta, 0) | {\bf J} | \; d\xi \, d\eta$

where $\tau$ is the constant shell thickness, $\hat{\Gamma}^e$ denotes the triangular surface of the shell mid-plane, $|{\bf J}|$ is the determinant of the transformation given by (3.182), and the Gauss quadrature rule for the 2D Triangle that is given in Table 3.1 is used to perform the numerical integration.

3.9.2.9. 3D Quadrilateral Shell Element

Beta Capability

Shell support in Aria is currently a beta feature. This capability is a work-in-progress and is not thoroughly tested and should be used with caution.

The three-dimensional, linear, quadrilateral shell element is, in principle, no different from the triangular shell element which we discussed in 3D Triangular Shell Element. As illustrated in Fig. 3.16, the quadrilateral shell has four nodes, instead of three, with the basis functions given in (3.178). The transformation from computational coordinates to physical coordinates is again given by (3.186), but with $I=4$ , and the $\hat{\psi}_i$ are defined by (3.178).

Fig. 3.16 The midplane of the master quadrilateral shell element is mapped to its image in physical coordinates $(x,y,z)$ , and the variation in the shell thickness is represented with the element shape functions.

Currently, Aria only supports constant thickness shells, and the simplifications discussed in 3D Triangular Shell Element also apply in this case. The volume integral therefore reduces to the product of the shell thickness and the integral over the mid-plane surface, with the transformation from computational to physical coordinates identical to that which is described in 2D Quadrilateral Element.