6  Finite Element Approximation in Multiple Dimensions

Chapter 5 developed the finite element method on an interval. The essential ideas were local even though the final algebraic problem was global: divide the domain into elements, construct polynomial basis functions on those elements, join neighboring pieces through shared degrees of freedom, evaluate the variational problem element by element, and assemble the local contributions into one sparse system.

The same structure carries over to two and three dimensions. What changes is the geometry. An interval element has two endpoints and one spatial coordinate. A triangle or tetrahedron has several facets and requires more than one coordinate; a quadrilateral or hexahedron introduces tensor-product basis functions and generally nonconstant geometry mappings. Gradients must be transformed between physical and reference coordinates, and natural boundary terms are integrated over element facets rather than over endpoints.

This chapter develops those changes without restarting the finite element method from the beginning. We first construct a linear triangular element directly on a physical triangle and then build the corresponding global space. The support of the global basis functions explains the sparsity of the assembled matrix. We next introduce the reference triangle, coordinate map, Jacobian, and quadrature rules before applying these constructions to one scalar-diffusion problem from element equations through a FEniCSx computation. The remaining sections extend the construction to higher-order and other element types and finally to vector-valued elasticity. The presentation follows the physical-element and local-to-global constructions commonly used in finite element texts and lecture notes (Jog 2011; Larson and Bengzon 2010; Arnold 2018).

6.1 From intervals to multidimensional meshes

Let \(\Om\subset\R^d\), with \(d=2\) or \(3\), be the computational domain. In one dimension, the mesh consisted of intervals. In multiple dimensions, the domain is partitioned into elements. Typical two-dimensional elements are triangles and quadrilaterals; typical three-dimensional elements are tetrahedra and hexahedra. For the straight-sided elements developed first, \(\Om\) is taken to be polygonal in two dimensions or polyhedral in three dimensions. Curved boundaries require an appropriate curved or isoparametric boundary representation, which is not developed in this chapter.

We write

\[ \mathcal T_h = \{\mathcal K_1,\mathcal K_2,\ldots,\mathcal K_{N_e}\} \tag{6.1}\]

for a mesh containing \(N_e\) elements. This chapter assumes a conforming, boundary-fitted mesh. The interiors of distinct elements do not overlap, the closures of the elements cover the computational domain,

\[ \overline{\Om} = \bigcup_{e=1}^{N_e}\overline{\mathcal K_e}. \]

In addition, the intersection of two distinct closed elements is either empty or a complete geometric entity shared by both elements. In two dimensions, the intersection may be a vertex or a complete edge. In three dimensions, it may be a vertex, a complete edge, or a complete face. Thus one element is not allowed to place a vertex in the interior of a neighboring element edge. Such a vertex is called a hanging node; meshes containing hanging nodes require additional constraints and are not considered in this chapter.

The lower-dimensional parts of an element also become important. In two dimensions, the boundary of an element consists of facets, which are edges. In three dimensions, the facets are two-dimensional faces. A triangle has three facets, a quadrilateral has four, a tetrahedron has four, and a hexahedron has six. Vertices and edges are lower-dimensional entities used to describe element geometry and finite element degrees of freedom.

Figure 6.1: Common element geometries in two and three dimensions.

Figure 6.2: Elements, vertices, shared facets, and boundary facets in a two-dimensional triangular mesh.

Two neighboring elements meet along a shared facet. A facet contained in \(\partial\Om\) is a boundary facet. This distinction matters because domain integrals are evaluated over elements, whereas natural boundary conditions produce integrals over boundary facets.

We also assume that the boundary mesh resolves the boundary-condition partition: every boundary facet belongs entirely to one prescribed boundary portion, such as \(\Gam_D\) or \(\Gam_N\). A change from one boundary type to another occurs only at a boundary vertex in two dimensions or along a boundary edge in three dimensions.

For later use, the size of an element is measured by its diameter,

\[ h_e = \operatorname{diam}(\mathcal K_e) = \sup_{\bx,\bm y\in\mathcal K_e}\norm{\bx-\bm y}_2, \qquad h=\max_e h_e. \tag{6.2}\]

The number \(h\) records the largest element size. Chapter 7 will add conditions on element shape when it studies approximation error; those conditions are not needed for the constructions in this chapter.

For a triangular mesh, let the global mesh vertices be

\[ \bx_1,\bx_2,\ldots,\bx_{N_v}\in\R^2. \]

Each triangle has three local vertices. The connectivity map

\[ I^e:\{1,2,3\}\longrightarrow\{1,\ldots,N_v\} \tag{6.3}\]

assigns each local vertex of element \(\mathcal K_e\) to its global vertex. If

\[ I^e=(4,7,9), \]

then local vertices 1, 2, and 3 of element \(e\) are global vertices 4, 7, and 9. The physical coordinates of the element are therefore

\[ \bx_1^e=\bx_4, \qquad \bx_2^e=\bx_7, \qquad \bx_3^e=\bx_9. \]

The geometry and connectivity play the same roles as in Chapter 5. The coordinates determine the physical shape of each element; the connectivity tells us which global degrees of freedom are shared by neighboring elements.

It is useful to distinguish a vertex, which is a geometric entity of the mesh, from an interpolation node, at which a nodal degree of freedom is defined. For the linear triangular and bilinear quadrilateral elements used in this chapter, the nodes are the element vertices. For higher-order elements, additional nodes may lie on edges, faces, or in the element interior.

A mesh may contain many more entities than vertices and elements. In practical finite element software, facets and edges are often stored or generated because boundary conditions, interface terms, and adjacency queries require them. For the constructions below, however, the key starting data are the element coordinates and connectivity.

Exercise 6.1 (Reading a triangular mesh) Suppose the global vertices of two triangles are numbered \(1,2,3,4\), with

\[ \bx_1=(0,0),\quad \bx_2=(1,0),\quad \bx_3=(0,1),\quad \bx_4=(1,1), \]

and the element connectivities are

\[ I^1=(1,2,3), \qquad I^2=(2,4,3). \]

1. Identify the facet shared by the two elements.
2. List the boundary facets of the two-element mesh.
3. State which global vertices receive contributions from both elements.
4. Explain why the shared facet is not part of the exterior boundary.
5. Check whether the two elements form a conforming mesh according to the definition above.

6.2 The \(P_1\) triangular finite element space

We begin with the simplest conforming two-dimensional element: a scalar linear triangular element. The element

\[ \mathcal K_e = \operatorname{conv} \{\bx_1^e,\bx_2^e,\bx_3^e\} \]

is the convex hull of the three noncollinear vertices \(\bx_a^e=(x_a^e,y_a^e)\), \(a=1,2,3\). In this setting, the convex hull is simply the filled triangle consisting of the vertices, edges, and interior.

The local polynomial space is

\[ \mathbb P_1(\mathcal K_e) = \operatorname{span}\{1,x,y\}. \tag{6.4}\]

It has dimension three. We use the three vertex evaluations

\[ N_a(w)=w(\bx_a^e), \qquad a=1,2,3, \]

as the local degrees of freedom. In the terminology of Definition 5.1,

\[ \left( \mathcal K_e, \mathbb P_1(\mathcal K_e), \{N_1,N_2,N_3\} \right) \]

is a finite element.

As in the one-dimensional construction, we want one local basis function for each nodal degree of freedom. The basis function \(\phi_a^e\) must be linear on the triangle and satisfy

\[ \phi_a^e(\bx_b^e)=\delta_{ab}, \qquad a,b\in\{1,2,3\}. \tag{6.5}\]

A linear polynomial in two variables has the form

\[ \phi_a^e(x,y) = \alpha_a+\beta_a x+\gamma_a y. \tag{6.6}\]

The three nodal conditions determine the three coefficients. For each local basis function,

\[ \begin{bmatrix} 1 & x_1^e & y_1^e\\ 1 & x_2^e & y_2^e\\ 1 & x_3^e & y_3^e \end{bmatrix} \begin{bmatrix} \alpha_a\\ \beta_a\\ \gamma_a \end{bmatrix} = \begin{bmatrix} \delta_{a1}\\ \delta_{a2}\\ \delta_{a3} \end{bmatrix}. \tag{6.7}\]

The determinant of the matrix on the left has magnitude \(2|\mathcal K_e|\). It is therefore nonsingular when the three vertices are not collinear. The vertex values uniquely determine a function in \(\mathbb P_1(\mathcal K_e)\), so the three degrees of freedom are unisolvent.

Because \(\phi_a^e\) is linear, its gradient is constant on the element:

\[ \grad\phi_a^e = \begin{bmatrix} \beta_a\\ \gamma_a \end{bmatrix}. \tag{6.8}\]

This is the two-dimensional analogue of the constant derivative of a one-dimensional linear shape function.

The three nodal basis functions are also called the barycentric coordinate functions of the triangle. The function \(\phi_a^e\) equals one at vertex \(a\), vanishes on the edge opposite that vertex, and varies linearly between these values. Inside the triangle,

\[ 0\le \phi_a^e\le 1, \qquad \phi_1^e+\phi_2^e+\phi_3^e=1. \tag{6.9}\]

Consequently, every point of the triangle can be written as

\[ \bx = \sum_{a=1}^{3}\phi_a^e(\bx)\bx_a^e. \tag{6.10}\]

The values \(\phi_a^e(\bx)\) are the barycentric coordinates of \(\bx\) relative to the three vertices.

6.2.1 A concrete physical triangle

Take the triangle with vertices

\[ \bx_1=(0,0), \qquad \bx_2=(2,0), \qquad \bx_3=(0,1). \]

Solving Equation 6.7 gives

\[ \phi_1(x,y) = 1-\frac{x}{2}-y, \qquad \phi_2(x,y) = \frac{x}{2}, \qquad \phi_3(x,y) = y. \tag{6.11}\]

Figure 6.3: Geometry, polynomial space, and nodal property of a physical \(P_1\) triangular element.

Figure 6.4: The three nodal basis functions over the physical triangle. Each basis function is one at its associated vertex and zero along the opposite edge.

Their gradients are

\[ \grad\phi_1 = \begin{bmatrix} -1/2\\ -1 \end{bmatrix}, \qquad \grad\phi_2 = \begin{bmatrix} 1/2\\ 0 \end{bmatrix}, \qquad \grad\phi_3 = \begin{bmatrix} 0\\ 1 \end{bmatrix}. \tag{6.12}\]

For this triangle, the general partition-of-unity identity in Equation 6.9 becomes

\[ \phi_1+\phi_2+\phi_3=1. \tag{6.13}\]

For nodal values \(U_1^e,U_2^e,U_3^e\), the local finite element field is

\[ u_h|_{\mathcal K_e} = U_1^e\phi_1^e + U_2^e\phi_2^e + U_3^e\phi_3^e. \tag{6.14}\]

The nodal interpolation property gives \(u_h(\bx_a^e)=U_a^e\). Thus the interpretation of the coefficients is exactly the same as in Chapter 5.

Example 6.1 (Interpolation on the physical triangle) For the triangle in Figure 6.3, let

\[ U_1=1,\qquad U_2=3,\qquad U_3=2. \]

Using Equation 6.11,

\[ u_h(x,y) = 1\left(1-\frac{x}{2}-y\right) +3\left(\frac{x}{2}\right) +2y = 1+x+y. \]

The gradient is therefore

\[ \grad u_h = \begin{bmatrix} 1\\ 1 \end{bmatrix}. \]

The local \(P_1\) approximation is a plane over the triangle, just as the one-dimensional \(P_1\) approximation is a straight line over an interval.

Exercise 6.2 (Basis functions on a general physical triangle) For the triangle with vertices \((1,0)\), \((3,0)\), and \((2,2)\):

1. solve Equation 6.7 for the three basis functions;
2. verify the nodal interpolation property;
3. verify that the basis functions sum to one;
4. compute their gradients;
5. identify the edge on which each basis function vanishes; and
6. use Equation 6.10 to recover the point \((2,1)\) from its barycentric coordinates.

Exercise 6.3 (Reproduction of an affine function) Let

\[ p(x,y)=c_0+c_1x+c_2y \]

be any affine function, and let \(p_a=p(\bx_a^e)\) be its values at the vertices of a nondegenerate triangle. Use Equation 6.10 and the partition-of-unity property to show that

\[ p(\bx)=\sum_{a=1}^{3}p_a\phi_a^e(\bx) \]

at every point \(\bx\in\mathcal K_e\). Thus a \(P_1\) triangle reproduces every affine scalar field exactly.

6.2.2 From local triangles to a global conforming space

A finite element function is constructed element by element, but the final function must be globally admissible. For the scalar diffusion problem, the natural energy space is \(H^1(\Om)\), so the finite element field must be continuous across shared facets.

For a triangular \(P_1\) mesh, define

\[ \mathcal V_h = \left\{ v_h\in C(\overline{\Om}): v_h|_{\mathcal K_e}\in\mathbb P_1(\mathcal K_e) \text{ for every }\mathcal K_e\in\mathcal T_h \right\}. \tag{6.15}\]

One global basis function \(\phi_i\) is associated with each global mesh vertex \(\bx_i\). It satisfies

\[ \phi_i(\bx_j)=\delta_{ij} \]

and is linear on every element. If \(\mathcal E(i)\) is the set of elements incident to vertex \(\bx_i\), then

\[ \operatorname{supp}\phi_i = \bigcup_{e\in\mathcal E(i)}\overline{\mathcal K_e}. \tag{6.16}\]

Thus \(\phi_i\) is zero outside the patch of elements that contain vertex \(i\).

The local and global basis functions are related exactly as in Chapter 5:

\[ \phi_{I^e(a)}|_{\mathcal K_e} = \phi_a^e. \tag{6.17}\]

This relation explains continuity across a shared facet. Consider two triangles that share an edge with endpoints \(\bx_i\) and \(\bx_j\). The restriction of a \(P_1\) function to that edge is a one-dimensional linear polynomial. Both elements use the same nodal values at \(\bx_i\) and \(\bx_j\), so both restrictions are the same linear function along the entire shared edge. The elementwise pieces therefore join continuously.

Figure 6.5: Continuity of two \(P_1\) pieces along a shared edge follows from their common endpoint values.

Every \(u_h\in\mathcal V_h\) has the global representation

\[ u_h(\bx) = \sum_{j=1}^{N_v} U_j\phi_j(\bx), \tag{6.18}\]

and the nodal property again gives

\[ U_j=u_h(\bx_j). \]

The difference from one dimension is geometric rather than algebraic. A global basis function may now be supported on several triangles surrounding a vertex, and the pattern of neighboring interactions depends on the mesh connectivity.

Figure 6.6: A global nodal basis function is nonzero only on the elements incident to its node.

This local support determines the sparsity of the stiffness matrix. For the diffusion problem introduced later in this chapter,

\[ K_{ij} = \int_\Om k\,\grad\phi_j\cdot\grad\phi_i\,\dx. \]

If the supports of \(\phi_i\) and \(\phi_j\) do not overlap on an element, the integrand is zero throughout the domain and \(K_{ij}=0\). A basis function therefore interacts only with basis functions whose node patches overlap its own patch.

Figure 6.7: Nonzero matrix entries follow the overlap of basis-function supports.

Exercise 6.4 (Continuity across a shared edge) Two \(P_1\) triangles share the edge connecting global vertices \(i\) and \(j\).

1. Restrict the finite element function from each triangle to the shared edge.
2. Show that each restriction is a linear function of one edge coordinate.
3. Explain why equality of the values at the two shared vertices implies equality along the complete edge.
4. State why the conclusion would fail if one triangle placed an additional, independent degree of freedom in the interior of the shared edge.

6.3 Reference triangle, mappings, and quadrature

The physical-triangle construction shows exactly what the basis functions mean, but an implementation should not solve a new interpolation system on every element. As in Chapter 5, we define basis functions once on a reference element and map them to the physical elements.

For triangles, use the reference element

\[ \widehat{\mathcal K} = \left\{ (\xi,\eta): \xi\ge0,\ \eta\ge0,\ \xi+\eta\le1 \right\}, \tag{6.19}\]

with reference vertices

\[ \widehat{\bx}_1=(0,0), \qquad \widehat{\bx}_2=(1,0), \qquad \widehat{\bx}_3=(0,1). \]

Figure 6.8: The reference \(P_1\) triangle, its local vertex numbering, and its reference coordinates.

We collect the reference coordinates as

\[ \bm{\xi} = \begin{bmatrix}\xi&\eta\end{bmatrix}^T. \]

The reference \(P_1\) basis functions are

\[ \widehat\phi_1(\xi,\eta) = 1-\xi-\eta, \qquad \widehat\phi_2(\xi,\eta) = \xi, \qquad \widehat\phi_3(\xi,\eta) = \eta. \tag{6.20}\]

The reference basis also constructs the geometry map. For a physical triangle with vertices \(\bx_1^e,\bx_2^e,\bx_3^e\), define

\[ \bx^e:\widehat{\mathcal K}\longrightarrow\mathcal K_e \]

by

\[ \bx = \bx^e(\bm\xi) = \sum_{a=1}^{3} \widehat\phi_a(\bm\xi)\bx_a^e. \tag{6.21}\]

The nodal property gives \(\bx^e(\widehat{\bx}_a)=\bx_a^e\), so each reference vertex is sent to the corresponding physical vertex. Substituting the three reference basis functions and collecting terms gives

\[ \bx = \bx^e(\bm{\xi}) = \bx_1^e + \bm{J}_e \begin{bmatrix} \xi\\ \eta \end{bmatrix}, \tag{6.22}\]

where

\[ \bm{J}_e = \begin{bmatrix} x_2^e-x_1^e & x_3^e-x_1^e\\ y_2^e-y_1^e & y_3^e-y_1^e \end{bmatrix}. \tag{6.23}\]

Thus

\[ \bm J_e = \frac{\partial\bx}{\partial\bm\xi}, \]

and its two columns are the physical edge vectors from local vertex 1 to local vertices 2 and 3.

The finite element field is interpolated on the same reference triangle:

\[ u_h(\bx^e(\bm\xi)) = \sum_{a=1}^{3} U_a^e\widehat\phi_a(\bm\xi). \tag{6.24}\]

The geometry and field therefore use the same \(P_1\) interpolation. According to the classification introduced in Section 5.6.2, this is an isoparametric description. The map is affine because the \(P_1\) basis functions are affine.

Figure 6.9: Affine mapping from the reference triangle to a physical triangle.

The determinant measures the change in area:

\[ \dd x = |\det\bm J_e|\, \dd\bm{\xi}, \tag{6.25}\]

where \(\dd\bm{\xi}=\dd\xi\,\dd\eta\). Since the reference triangle has area \(1/2\),

\[ |\mathcal K_e| = \frac12|\det\bm J_e|. \tag{6.26}\]

A valid affine triangle map requires \(\det\bm J_e\ne0\). In a consistently oriented mesh, the local vertex ordering is usually chosen so that \(\det\bm J_e>0\). In two dimensions this is obtained by listing the three local vertices counterclockwise.

6.3.1 Transformation of gradients

The physical basis function is

\[ \phi_a^e(\bx) = \widehat\phi_a((\bx^e)^{-1}(\bx)). \]

The gradient transformation follows from the chain rule. Since

\[ \widehat\phi_a(\bm{\xi}) = \phi_a^e(\bx^e(\bm{\xi})), \]

differentiation with respect to the reference coordinate gives

\[ \grad_{\bm{\xi}}\widehat\phi_a = \bm J_e^T \grad_{\bx}\phi_a^e. \]

Therefore,

\[ \boxed{ \grad_{\bx}\phi_a^e = \bm J_e^{-T} \grad_{\bm{\xi}}\widehat\phi_a. } \tag{6.27}\]

This formula is the multidimensional counterpart of the one-dimensional factor \(1/J_e\) in Chapter 5. The difference is that the Jacobian is now a matrix and the gradient is a vector.

For a two-dimensional triangle, \(\bm J_e\in\R^{2\times2}\), \(\grad_{\bm{\xi}}\widehat\phi_a\in\R^2\), and \(\grad_{\bx}\phi_a^e\in\R^2\). The inverse transpose follows directly from solving the chain-rule relation \(\grad_{\bm{\xi}}\widehat\phi_a=\bm J_e^T\grad_{\bx}\phi_a^e\) for the physical gradient.

Example 6.2 (Mapping the concrete triangle) For the physical triangle with vertices \((0,0)\), \((2,0)\), and \((0,1)\),

\[ \bm J_e = \begin{bmatrix} 2&0\\ 0&1 \end{bmatrix}, \qquad \det\bm J_e=2. \]

Hence \(|\mathcal K_e|=\tfrac12(2)=1\), as used previously. Moreover,

\[ \grad_{\bx}\phi_a^e = \begin{bmatrix} 1/2&0\\ 0&1 \end{bmatrix} \grad_{\bm{\xi}}\widehat\phi_a, \]

which recovers the physical gradients in Equation 6.12.

Exercise 6.5 (Checking the gradient transformation) For a physical triangle with vertices \((1,1)\), \((4,1)\), and \((2,3)\):

1. form \(\bm J_e\) and compute its determinant;
2. calculate the physical area;
3. compute \(\bm J_e^{-T}\);
4. transform the three reference gradients; and
5. construct the physical basis functions \(\phi_a^e=\widehat\phi_a\circ (\bx^e)^{-1}\) and verify that their directly computed gradients agree with the transformed gradients from Part 4.

6.3.2 Numerical integration on triangles and edges

After transformation to a reference element, a typical element integral has the form

\[ \int_{\mathcal K_e} g(\bx)\,\dx = \int_{\widehat{\mathcal K}} g(\bx^e(\bm{\xi})) |\det\bm J_e(\bm{\xi})| \,\dd\bm{\xi}. \tag{6.28}\]

For affine triangles, \(\bm J_e\) is constant. For more general isoparametric elements, the Jacobian may vary across the element.

A quadrature rule on the reference triangle approximates

\[ \int_{\widehat{\mathcal K}} g(\xi,\eta)\,\dd\xi\,\dd\eta \approx \sum_{r=1}^{n_q} w_r g(\xi_r,\eta_r). \tag{6.29}\]

For the reference triangle in Equation 6.19, useful low-order rules include:

rule	points	weights	polynomial exactness
one point	\((1/3,1/3)\)	\(1/2\)	degree 1
three point	\((1/6,1/6)\), \((2/3,1/6)\), \((1/6,2/3)\)	\(1/6\) each	degree 2

The weights sum to \(1/2\), the area of the reference triangle.

Figure 6.10: One-point and three-point quadrature rules on the reference triangle.

The transformed scalar-diffusion stiffness entry is

\[ K_{ab}^e = \int_{\widehat{\mathcal K}} k(\bx^e(\bm{\xi})) \left( \bm J_e^{-T} \grad_{\bm{\xi}}\widehat\phi_b \right) \cdot \left( \bm J_e^{-T} \grad_{\bm{\xi}}\widehat\phi_a \right) |\det\bm J_e| \,\dd\bm{\xi}. \tag{6.30}\]

For a \(P_1\) affine triangle with constant \(k\), the physical gradients and the Jacobian are constant, so the stiffness integrand is constant. A one-point rule is therefore sufficient. More complicated coefficients, higher-order basis functions, or nonaffine geometry may require higher-order quadrature.

For example, the element load entry is approximated by

\[ F_a^e \approx \sum_{r=1}^{n_q} w_r f(\bx^e(\bm\xi_r)) \widehat\phi_a(\bm\xi_r) |\det\bm J_e(\bm\xi_r)|. \tag{6.31}\]

For constant \(f\) on an affine \(P_1\) triangle, the integrand is linear. The one-point rule gives

\[ F_a^e = \frac12 f\frac13|\det\bm J_e| = \frac{f|\mathcal K_e|}{3}, \]

which agrees with Equation 6.44.

Boundary-facet integrals are transformed separately. In two dimensions a facet is a line segment, so the one-dimensional quadrature rules from Chapter 5 can be reused on a reference edge. In three dimensions a facet is a two-dimensional surface and requires a surface quadrature rule.

Exercise 6.6 (Triangle quadrature)  

1. Verify that the one-point rule integrates the three reference basis functions exactly.
2. Verify that the three-point rule integrates \(\xi^2\), \(\xi\eta\), and \(\eta^2\) exactly.
3. Explain why the one-point rule is sufficient for the stiffness matrix of a \(P_1\) affine triangle with constant diffusivity.
4. Determine whether the one-point rule exactly evaluates the load vector when \(f(x,y)=x+y\) and the geometry map is affine.
5. For the \(P_2\) triangle of Section 6.5.1 on an affine element with constant diffusivity, determine the polynomial degree of the stiffness integrand in Equation 6.30 and decide whether the one-point or the three-point rule integrates it exactly.

6.4 Scalar diffusion: from the variational problem to computation

We now apply the multidimensional space to the scalar diffusion problem from Section 3.2.3. Let \(\Gam_D\) and \(\Gam_N\) partition the boundary, and consider

\[ -\divg(k\grad u)=f \qquad\text{in }\Om, \]

with

\[ u=g\quad\text{on }\Gam_D, \qquad \bq\cdot\bn=q \quad\text{on }\Gam_N. \]

Here

\[ \bq=k\grad u, \qquad q=\bq\cdot\bn. \]

The scalar \(q=\bq\cdot\bn\) is the prescribed normal component. In heat conduction, physical heat flow is \(-\bq\), so \(q>0\) means heat enters the body. The symbol \(h\) remains reserved for the mesh size.

The variational problem is

\[ \text{find }u\in\mathcal U_g \text{ such that}\qquad a(u,v)=\ell(v) \qquad\forall v\in\mathcal V_0, \]

where

\[ a(w,v) = \int_\Om k\,\grad w\cdot\grad v\,\dx, \]

and

\[ \ell(v) = \int_\Om f v\,\dx + \int_{\Gam_N}q v\,\ds. \tag{6.32}\]

The homogeneous finite element test space is

\[ \mathcal V_{h,0} = \mathcal V_h\cap\mathcal V_0 = \left\{ v_h\in\mathcal V_h: v_h=0\text{ on }\Gam_D \right\}. \tag{6.33}\]

To describe nonzero prescribed data, suppose \(g\) is continuous on \(\Gam_D\) and let \(g_h\) denote its nodal finite element interpolant on that boundary. The finite element trial space is

\[ \mathcal U_{h,g_h} = \left\{ w_h\in\mathcal V_h: w_h=g_h\text{ on }\Gam_D \right\}. \tag{6.34}\]

In general, \(g_h\) is an approximation of \(g\) rather than the exact boundary function. If \(g\) is itself piecewise linear on the boundary mesh, then \(g_h=g\). The discrete variational problem is

\[ \text{find }u_h\in\mathcal U_{h,g_h} \text{ such that}\qquad a(u_h,v_h)=\ell(v_h) \qquad \forall v_h\in\mathcal V_{h,0}. \tag{6.35}\]

With the global basis, \(u_h=\sum_j U_j\phi_j\). This sum includes both prescribed and unknown nodal coefficients. Testing with a basis function \(\phi_i\) associated with a free degree of freedom gives

\[ \sum_j K_{ij}U_j = F_i, \]

where

\[ K_{ij} = \int_\Om k\,\grad\phi_j\cdot\grad\phi_i\,\dx, \]

and

\[ F_i = \int_\Om f\phi_i\,\dx + \int_{\Gam_N}q\phi_i\,\ds. \tag{6.36}\]

After assembly, the terms containing prescribed coefficients are moved to the right-hand side, exactly as in Equation 5.46.

The global domain integral is decomposed element by element:

\[ K_{ij} = \sum_{e=1}^{N_e} \int_{\mathcal K_e} k\,\grad\phi_j\cdot\grad\phi_i\,\dx. \]

Using Equation 6.17, define the local stiffness entries

\[ K_{ab}^e = \int_{\mathcal K_e} k\, \grad\phi_b^e\cdot\grad\phi_a^e \,\dx, \tag{6.37}\]

and the local domain-load entries

\[ F_a^e = \int_{\mathcal K_e} f\,\phi_a^e\,\dx. \tag{6.38}\]

Natural boundary data contribute through boundary facets. If a facet \(\mathcal F\subset\partial\mathcal K_e\cap\Gam_N\), its contribution is

\[ F_a^{\mathcal F} = +\int_{\mathcal F} q\,\phi_a^e\,\ds. \tag{6.39}\]

To evaluate this integral, map the reference interval \(\widehat{\mathcal F}=[-1,1]\) to the straight physical edge with endpoints \(\bx_i\) and \(\bx_j\):

\[ \bx(\zeta) = \frac{1-\zeta}{2}\bx_i + \frac{1+\zeta}{2}\bx_j. \tag{6.40}\]

The edge Jacobian is the length of the tangent vector,

\[ J_{\mathcal F} = \norm{\frac{\dd\bx}{\dd\zeta}}_2 = \frac{\norm{\bx_j-\bx_i}_2}{2} = \frac{L_{\mathcal F}}{2}. \tag{6.41}\]

Therefore,

\[ F_a^{\mathcal F} = +\int_{-1}^{1} q(\bx(\zeta)) \phi_a^e(\bx(\zeta)) J_{\mathcal F} \,\dd\zeta. \tag{6.42}\]

This is an ordinary one-dimensional integral on the reference edge.

Figure 6.11: A natural boundary condition contributes through an integral over a boundary facet.

The global load vector therefore receives contributions from both elements and boundary facets.

Example 6.3 (Element matrix for the concrete triangle) For the triangle in Figure 6.3, the area is

\[ |\mathcal K_e|=1. \]

If \(k\) is constant on the element, the basis gradients in Equation 6.12 are constant, so

\[ K_{ab}^e = k|\mathcal K_e| \grad\phi_b^e\cdot\grad\phi_a^e. \]

Hence

\[ \bK^e = k \begin{bmatrix} 5/4 & -1/4 & -1\\ -1/4 & 1/4 & 0\\ -1 & 0 & 1 \end{bmatrix}. \tag{6.43}\]

If the source \(f\) is constant, each linear basis function has average value \(1/3\) over the triangle. Therefore

\[ \bF^e = \frac{f|\mathcal K_e|}{3} \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix} = \frac{f}{3} \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}. \tag{6.44}\]

As in one dimension,

\[ \bK^e \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix} = \bm{0}, \]

because a constant nodal field has zero gradient.

Example 6.4 (A constant flux on one edge) Let a boundary edge of a \(P_1\) triangle have length \(L_{\mathcal F}\) and let the prescribed normal datum be the constant \(q_0\). Only the two vertex basis functions associated with the edge are nonzero on that facet. Their restrictions are \((1-\zeta)/2\) and \((1+\zeta)/2\). For the endpoint whose restriction is \((1-\zeta)/2\), Equation 6.42 gives

\[ \int_{\mathcal F}q_0\phi_a^e\,\ds = \frac{q_0L_{\mathcal F}}{2} \int_{-1}^{1}\frac{1-\zeta}{2}\,\dd\zeta = \frac{q_0L_{\mathcal F}}{2}, \]

and the same value is obtained for the other endpoint. Hence the facet load vector is

\[ \bF^{\mathcal F} = \frac{q_0L_{\mathcal F}}{2} \begin{bmatrix} 1\\ 1 \end{bmatrix}. \]

These two entries are assembled into the global load vector at the two global vertices of the boundary edge.

Exercise 6.7 (A linearly varying boundary flux) On a straight edge of length \(L_{\mathcal F}\), let the prescribed flux vary linearly from \(q_i\) at endpoint \(i\) to \(q_j\) at endpoint \(j\).

1. Express \(q(\zeta)\) using the two linear reference-edge basis functions.
2. Use Equation 6.42 to derive the two-entry consistent facet load vector.
3. Verify that the sum of the entries equals \(L_{\mathcal F}(q_i+q_j)/2\). For heat conduction, verify that the physical outward heat-flow integral is \(-L_{\mathcal F}(q_i+q_j)/2\).

Exercise 6.8 (A \(P_1\) diffusion element) For the triangle with vertices \((0,0)\), \((1,0)\), and \((0,1)\):

1. derive the three basis functions;
2. compute their gradients;
3. derive the element stiffness matrix for constant \(k\);
4. derive the element load vector for constant \(f\); and
5. verify that the row sums of the stiffness matrix are zero.

Figure 6.12: A two-triangle mesh used to reproduce the linear solution \(u(x,y)=x\).

Example 6.5 (Assembly and solution on two triangles) Consider the unit square with vertices

\[ \bx_1=(0,0),\quad \bx_2=(1,0),\quad \bx_3=(0,1),\quad \bx_4=(1,1), \]

and triangle connectivities

\[ I^1=(1,2,3), \qquad I^2=(2,4,3). \]

Let \(k=1\) and \(f=0\). Prescribe \(u=0\) on the left boundary, \(q=1\) on the right boundary, and \(q=0\) on the top and bottom boundaries. Because \(\bq=\grad u\), these are the boundary data satisfied by \(u(x,y)=x\).

Direct calculation gives

\[ \bK^1 = \frac12 \begin{bmatrix} 2&-1&-1\\ -1&1&0\\ -1&0&1 \end{bmatrix}, \qquad \bK^2 = \frac12 \begin{bmatrix} 1&-1&0\\ -1&2&-1\\ 0&-1&1 \end{bmatrix}. \]

Scatter-adding these matrices with the two connectivity maps produces

\[ \bK = \begin{bmatrix} 1&-1/2&-1/2&0\\ -1/2&1&0&-1/2\\ -1/2&0&1&-1/2\\ 0&-1/2&-1/2&1 \end{bmatrix}. \tag{6.45}\]

The right boundary edge has length one and \(q=1\). The positive sign in the variational boundary term gives two positive load entries. The constant-flux calculation in Example 6.4 gives

\[ \bF = \begin{bmatrix} 0&1/2&0&1/2 \end{bmatrix}^T. \]

The entries are ordered by global vertices \(1,2,3,4\).

The prescribed coefficients are \(U_1=U_3=0\). The equations for the two free coefficients are

\[ \begin{bmatrix} 1&-1/2\\ -1/2&1 \end{bmatrix} \begin{bmatrix} U_2\\U_4 \end{bmatrix} = \begin{bmatrix} 1/2\\1/2 \end{bmatrix}, \]

which gives \(U_2=U_4=1\). Therefore the assembled finite element field is \(u_h(x,y)=x\) on both triangles. This linear field satisfies the differential equation and all four boundary conditions, so the two-element mesh reproduces the exact solution.

Exercise 6.9 (Reusing the two-triangle mesh) On the mesh in Example 6.5, set \(k=1\), \(f=1\), prescribe \(u=0\) on the left boundary, and prescribe zero normal flux elsewhere.

1. Derive the element load vector on each triangle.
2. Assemble the global load vector.
3. Solve the two-by-two reduced system for \(U_2\) and \(U_4\).
4. Write the resulting piecewise-linear function separately on each triangle.

6.4.1 Assembly and boundary conditions

For scalar \(P_1\) diffusion, each triangular element has three local degrees of freedom. For two-dimensional displacement elasticity, the same triangle has six. The local-to-global assembly rule itself is unchanged.

If local scalar indices \(a,b\) correspond to global indices \(I^e(a),I^e(b)\),

\[ K_{I^e(a),I^e(b)} \mathrel{+}= K_{ab}^e, \qquad F_{I^e(a)} \mathrel{+}= F_a^e. \tag{6.46}\]

Figure 6.13: Scatter-add assembly of a three-by-three element matrix into the rows and columns selected by its connectivity.

The gather form of the same rule also carries over unchanged. Let \(\mathcal E(i)\) now denote the set of elements incident on global vertex \(i\). For \(e\in\mathcal E(i)\), the inverse \((I^e)^{-1}(i)\) gives the local number of that vertex in element \(e\). The argument that led to Equation 5.32 gives

\[ K_{ij} = \sum_{e\in\mathcal E(i)\cap\mathcal E(j)} K_{(I^e)^{-1}(i),(I^e)^{-1}(j)}^e, \qquad F_i = \sum_{e\in\mathcal E(i)}F_{(I^e)^{-1}(i)}^e \tag{6.47}\]

for the element contributions, with boundary-facet contributions added to \(F_i\) in the same way. The formulas are identical to those of Chapter 5; what changes is the size of the incident sets. In one dimension \(\mathcal E(i)\) contains at most two elements, whereas around an interior vertex of a triangular mesh it contains every triangle of the surrounding patch.

For the two-triangle mesh of Example 6.5, vertex 2 is local vertex 2 of \(\mathcal K_1\) and local vertex 1 of \(\mathcal K_2\), and vertex 3 is local vertex 3 of both elements. Hence

\[ K_{22}=K_{22}^1+K_{11}^2=\tfrac12+\tfrac12=1, \qquad K_{23}=K_{23}^1+K_{13}^2=0+0=0, \]

in agreement with Equation 6.45. The second entry is zero even though vertices 2 and 3 share an edge. On each of these right triangles the gradients of the two corresponding basis functions are orthogonal. Mesh adjacency determines which entries can be nonzero, but an individual entry within that pattern may still vanish because of the element geometry.

For vector-valued elements, the element degree-of-freedom map replaces the vertex map in this rule. A vertex shared by several elements contributes the same global displacement components to each incident element.

The global stiffness matrix remains sparse because a basis function interacts only with basis functions supported on elements that overlap its support. In two and three dimensions, the sparsity pattern is no longer a simple narrow band under arbitrary numbering, but it still reflects mesh adjacency.

Essential boundary conditions are imposed on global degrees of freedom associated with the prescribed boundary. The same partitioning from Chapter 5 applies:

\[ \bK_{ff}\bU_f = \bF_f-\bK_{fp}\bU_p. \tag{6.48}\]

The geometric task is more involved because a boundary condition is usually specified on a boundary portion rather than at a single endpoint. The mesh must therefore identify the boundary facets belonging to \(\Gam_D\) or \(\Gam_u\), and the finite element software must identify the degrees of freedom attached to those facets.

Natural boundary conditions require the same boundary-facet identification, but their contribution is assembled into the load vector through facet integrals rather than through prescribed coefficient values.

After the solve, the prescribed values are inserted into \(\bU\) and the residual vector \(\bR=\bF-\bK\bU\) is formed exactly as in Equation 5.48. Its free entries vanish up to rounding error, and its negative prescribed entries are the discrete reactions, \(\bm r_p=-\bR_p\). For elasticity they are the nodal support forces. For scalar diffusion, the entry associated with a prescribed nodal degree of freedom \(i\) is \(r_i=a(u_h,\phi_i)-\ell(\phi_i)\). If \(u_h\) is replaced by a sufficiently smooth exact solution, integration by parts shows that this quantity equals \(\int_{\Gam_D}q\,\phi_i\,\ds\). The reactions are discrete counterparts of these weighted boundary fluxes; they need not equal direct boundary integrals of the computed gradient. Their sum is also a useful check. The nodal basis functions sum to one, the free equations give \(a(u_h,\phi_i)=\ell(\phi_i)\) for every free nodal degree of freedom, and \(a(u_h,1)=0\). Together these facts give

\[ \sum_{i\text{ prescribed}}r_i = -\ell(1) = -\int_\Om f\,\dx-\int_{\Gam_N}q\,\ds. \]

The discrete reactions therefore balance the source and the prescribed boundary flux exactly. This is the multidimensional counterpart of the equilibrium check in Example 5.4.

6.4.2 Complete multidimensional finite element procedure

The multidimensional finite element method can now be organized as the same local-to-global procedure used in one dimension, with additional geometric transformations.

read the mesh coordinates, element connectivity, and boundary tags
choose the finite element and quadrature rules
initialize the global matrix K and vector F

for each element e:
    read the element coordinates and global degree-of-freedom map
    initialize the element matrix K^e and domain-load vector F^e

    for each element quadrature point:
        evaluate reference basis functions and reference gradients
        map the quadrature point to the physical element
        evaluate J_e, det(J_e), and J_e^{-T}
        transform the basis gradients to physical coordinates
        evaluate coefficients and source or body-force data
        add the quadrature contribution to K^e and F^e

    scatter-add K^e and F^e into K and F

for each boundary facet carrying natural data:
    read the facet degree-of-freedom map and geometry
    initialize the facet load vector F^facet

    for each facet quadrature point:
        evaluate the restricted basis functions and facet Jacobian
        evaluate the prescribed flux or traction
        add the positive natural-boundary load contribution to F^facet

    scatter-add F^facet into F

identify prescribed and free degrees of freedom
solve K_ff U_f = F_f - K_fp U_p
insert prescribed values into the global coefficient vector
postprocess flux, strain, stress, and reactions as required

The outer algorithm is largely independent of whether the element is a triangle, quadrilateral, tetrahedron, or hexahedron. What changes are the reference basis, geometry map, Jacobian, quadrature rule, and number of local degrees of freedom.

Exercise 6.10 (Trace one triangle through the complete algorithm) Use element \(\mathcal K_1\) with connectivity \(I^1=(1,2,3)\) in Example 6.5.

1. Write its coordinate array, reference basis, and reference gradients.
2. Compute \(\bm J_1\), \(\det\bm J_1\), \(\bm J_1^{-T}\), and the physical basis gradients.
3. Use one-point quadrature to recover \(\bK^1\).
4. List all nine global matrix entries that receive contributions from \(\bK^1\).
5. Identify the boundary facet of \(\mathcal K_2\) that produces the right-edge flux vector and list the two global vector entries that receive it.
6. Starting from the assembled system, reproduce the reduced two-by-two system and solution reported in the example.

6.4.3 FEniCSx implementation and mesh refinement

The mathematical construction above corresponds closely to the objects used by FEniCSx. The software automates reference-element evaluation, coordinate mapping, quadrature, degree-of-freedom maps, local element kernels, and sparse assembly, but the objects being manipulated are the same ones developed in this chapter. The theory in Chapters 5 and 6 uses one-based numbering. Python arrays and FEniCSx degree-of-freedom maps use zero-based indices, so mathematical index \(i\) is represented by Python index i - 1 when these objects are stored in code. Accordingly, x[0] in the example below denotes the first spatial coordinate.

Mathematical object	FEniCSx counterpart
mesh \(\mathcal T_h\)	dolfinx.mesh.Mesh
element type and geometry	mesh topology and geometry
finite element space \(\mathcal V_h\)	fem.functionspace(...)
trial and test functions	ufl.TrialFunction, ufl.TestFunction
gradient and inner products	UFL grad, inner
domain integral	ufl.dx
boundary-facet integral	ufl.ds
Dirichlet degrees of freedom	fem.locate_dofs_*
essential boundary condition	fem.dirichletbc
assembled linear problem	dolfinx.fem.petsc.LinearProblem

A manufactured solution gives a direct check of the implementation. On the unit square, choose

\[ u_\star(x,y)=1+x+x(1-x)\sin(\pi y). \]

Prescribe \(u=g=1+x\) on the vertical sides and \(q=\grad u\cdot\bn\) on the horizontal sides. Substitution gives

\[ f(x,y)=\left[2+\pi^2x(1-x)\right]\sin(\pi y), \qquad q(x,y)=-\pi x(1-x). \tag{6.49}\]

The same expression for \(q\) holds on the top and bottom edges because the outward unit normal changes direction between the two edges. The domain and boundary partition are shown in Figure 6.14.

Figure 6.14: Mixed boundary conditions for the manufactured diffusion problem.

The central part of the FEniCSx implementation is:

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh
from dolfinx.fem.petsc import LinearProblem

domain = mesh.create_unit_square(
    MPI.COMM_WORLD, 16, 16, cell_type=mesh.CellType.triangle)

V = fem.functionspace(domain, ("Lagrange", 1))

x = ufl.SpatialCoordinate(domain)
u_exact = 1 + x[0] + x[0] * (1 - x[0]) * ufl.sin(np.pi * x[1])
f = (2 + np.pi**2 * x[0] * (1 - x[0])) * ufl.sin(np.pi * x[1])
q = -np.pi * x[0] * (1 - x[0])

fdim = domain.topology.dim - 1
vertical_facets = mesh.locate_entities_boundary(
    domain, fdim,
    lambda X: np.isclose(X[0], 0.0) | np.isclose(X[0], 1.0)
)
horizontal_facets = mesh.locate_entities_boundary(
    domain, fdim,
    lambda X: np.isclose(X[1], 0.0) | np.isclose(X[1], 1.0)
)
facet_tags = mesh.meshtags(
    domain, fdim, horizontal_facets,
    np.ones(len(horizontal_facets), dtype=np.int32)
)

prescribed_dofs = fem.locate_dofs_topological(
    V, fdim, vertical_facets)
g = fem.Function(V)
g.interpolate(lambda X: 1.0 + X[0])
bc = fem.dirichletbc(g, prescribed_dofs)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx + q * v * ds(1)

problem = LinearProblem(
    a,
    L,
    bcs=[bc],
    petsc_options_prefix="diffusion_",
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
)

u_h = problem.solve()
error_form = fem.form((u_h - u_exact) ** 2 * ufl.dx)
error_local = fem.assemble_scalar(error_form)
error_L2 = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))

if domain.comm.rank == 0:
    print(f"L2 error: {error_L2:.6e}")

The boundary-facet term in L has a plus sign because \(q=\bq\cdot\bn\) and \(\bq=\grad u\). The notebook diffusion_2d_demo.ipynb contains the complete calculation, including boundary tags, plots, error evaluation, and a refinement loop.

The UFL expression

ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx

represents the same bilinear form as

\[ \int_\Om k\,\grad u_h\cdot\grad v_h\,\dx. \]

FEniCSx evaluates this form element by element, transforms basis gradients from reference to physical coordinates, applies quadrature, and assembles the local contributions into a global PETSc matrix. The code does not replace the derivation; it implements it.

Figure 6.15 shows the finite element solution on the coarsest mesh. Refining the mesh from 4 to 16 and then 64 triangles gives the results in Figure 6.16.

Figure 6.15: Finite element solution on a four-triangle mesh.

Figure 6.16: Finite element solutions on three successively refined triangular meshes.

number of triangles	\(h\)	\(L^2\) error	observed rate
4	1.118	\(1.291\times10^{-1}\)	–
16	0.559	\(3.870\times10^{-2}\)	1.74
64	0.280	\(1.079\times10^{-2}\)	1.84

The error decreases as the mesh is refined. The rates approach two, which is the \(L^2\) rate predicted for sufficiently smooth solutions approximated by \(P_1\) elements. Chapter 7 develops this estimate and its assumptions.

Figure 6.17: The \(L^2\) error decreases under mesh refinement.

Exercise 6.11 (Verifying the FEniCSx solution) Using the accompanying notebook:

1. Identify the mesh, finite element space, trial function, test function, bilinear form, linear functional, and prescribed boundary condition in the code.
2. Verify Equation 6.49 by differentiating \(u_\star\) and evaluating \(q=\grad u_\star\cdot\bn\) on both horizontal edges.
3. Run the calculation with \((n_x,n_y)=(1,2)\), \((2,4)\), and \((4,8)\) and reproduce the error table above.
4. Explain which steps of the element loop in Section 6.4.2 are carried out internally by FEniCSx.

6.5 Extensions to other elements and polynomial degrees

6.5.1 Quadratic triangular elements

Chapter 5 raised the polynomial degree on an interval by adding an interior node. The corresponding construction on a triangle uses the six-dimensional space

\[ \mathbb P_2(\widehat{\mathcal K}) = \operatorname{span}\{1,\xi,\eta,\xi^2,\xi\eta,\eta^2\}. \tag{6.50}\]

Six degrees of freedom are therefore needed. The quadratic Lagrange triangle uses point evaluation at the three vertices and at the three edge midpoints (Larson and Bengzon 2010, 152–54). We number the local nodes as

\[ \widehat{\bx}_1=(0,0),\quad \widehat{\bx}_2=(1,0),\quad \widehat{\bx}_3=(0,1),\quad \widehat{\bx}_4=\left(\tfrac12,0\right),\quad \widehat{\bx}_5=\left(\tfrac12,\tfrac12\right),\quad \widehat{\bx}_6=\left(0,\tfrac12\right), \]

so that nodes 4, 5, and 6 are the midpoints of the edges joining vertices \(1\)–\(2\), \(2\)–\(3\), and \(3\)–\(1\), respectively. Software libraries may order the edge nodes differently; as with vector components, every element array must use one consistent ordering.

Figure 6.18: Local nodes of the quadratic reference triangle: three vertices and three edge midpoints.

The basis is most easily written with the reference \(P_1\) functions Equation 6.20, which are the barycentric coordinates of the reference triangle:

\[ \lambda_1=1-\xi-\eta, \qquad \lambda_2=\xi, \qquad \lambda_3=\eta. \]

The six quadratic nodal basis functions are

\[ \begin{aligned} \widehat\phi_a&=\lambda_a(2\lambda_a-1),\qquad a=1,2,3,\\ \widehat\phi_4&=4\lambda_1\lambda_2, \qquad \widehat\phi_5=4\lambda_2\lambda_3, \qquad \widehat\phi_6=4\lambda_3\lambda_1. \end{aligned} \tag{6.51}\]

Each is a product of two linear functions and so belongs to \(\mathbb P_2(\widehat{\mathcal K})\). The nodal property \(\widehat\phi_a(\widehat{\bx}_b)=\delta_{ab}\) can be checked directly from the barycentric values: at a vertex one \(\lambda_a\) equals one and the others are zero, and at an edge midpoint the two barycentric coordinates of that edge equal \(1/2\) while the third is zero. For example, \(\widehat\phi_1\) vanishes at vertices 2 and 3 because \(\lambda_1=0\) there, and at the midpoints of the two edges through vertex 1 because \(2\lambda_1-1=0\) there. Six functions satisfying the nodal property are linearly independent, so the six degrees of freedom are unisolvent on \(\mathbb P_2(\widehat{\mathcal K})\).

Continuity between neighboring elements follows from the same edge argument used for \(P_1\) triangles. On the edge joining vertices 1 and 2, \(\lambda_3=0\), so only \(\widehat\phi_1\), \(\widehat\phi_2\), and \(\widehat\phi_4\) are nonzero, and the restriction of the field to that edge is a quadratic polynomial in one edge coordinate. A quadratic polynomial in one variable is determined by its values at three distinct points. Two elements sharing an edge use the same two vertex degrees of freedom and the same midpoint degree of freedom on that edge, so their restrictions agree along the entire edge. Consequently, a global \(P_2\) space contains one degree of freedom for each mesh vertex and one for each mesh edge. A vertex degree of freedom is shared by every element incident on that vertex, whereas an edge degree of freedom is shared only by the elements that contain the edge.

There are two equivalent ways to describe a straight-sided quadratic triangle. The field may use the six quadratic basis functions while the geometry retains the affine \(P_1\) map Equation 6.22; this is a subparametric description. Alternatively, an isoparametric quadratic geometry map is

\[ \bx^{e,(2)}(\bm\xi) = \sum_{a=1}^{6} \widehat\phi_a(\bm\xi)\bx_a^e. \tag{6.52}\]

When the three physical edge nodes are placed at the edge midpoints, the quadratic terms in this map cancel and \(\bx^{e,(2)}\) reduces exactly to the affine map Equation 6.22. Moving an edge geometry node away from the midpoint produces a curved edge and a nonconstant Jacobian. Curved quadratic geometry is not developed further in this chapter.

For either straight-sided description, the physical basis functions are \(\phi_a^e=\widehat\phi_a\circ (\bx^e)^{-1}\), their gradients follow from Equation 6.27, and the physical interpolation nodes are the images of the six reference nodes. Here \(\bx^e\) denotes the geometry map selected for the element.

Unlike the \(P_1\) case, the gradients of the quadratic basis functions are linear rather than constant on the element. The quadrature requirement for a \(P_2\) element matrix is therefore higher, as Exercise 6.6 asks you to verify.

Exercise 6.12 (The six-node triangle)  

1. Verify the nodal property of Equation 6.51 at all six reference nodes.
2. Using \(\lambda_1+\lambda_2+\lambda_3=1\), show that the six functions sum to one.
3. Restrict Equation 6.51 to the edge \(\eta=0\) and write the three nonzero restrictions as functions of \(\xi\). Compare them with the quadratic reference basis of Chapter 5 after mapping \([0,1]\) to \([-1,1]\).
4. On the two-triangle mesh in Exercise 6.1, count the global vertex and edge degrees of freedom of a scalar \(P_2\) space, and state which of them are shared by both elements.

6.5.2 Quadrilateral elements

An affine map sends the reference triangle to any physical triangle. A general four-sided element requires a different construction because an affine map can produce only a parallelogram. The standard four-node quadrilateral begins with the reference square

\[ \widehat{\mathcal K}=[-1,1]^2 \]

and the tensor-product polynomial space

\[ \mathbb Q_1(\widehat{\mathcal K}) = \operatorname{span}\{1,\xi,\eta,\xi\eta\}. \tag{6.53}\]

The notation \(\mathbb Q_1\) means that the polynomial degree is at most one in each coordinate separately. The product \(\xi\eta\) is therefore included even though its total degree is two. The four local nodes are ordered

\[ (-1,-1),\qquad (1,-1),\qquad (1,1),\qquad (-1,1), \]

and the corresponding nodal basis functions are

\[ \begin{aligned} \widehat\phi_1(\xi,\eta) &=\frac14(1-\xi)(1-\eta),& \widehat\phi_2(\xi,\eta) &=\frac14(1+\xi)(1-\eta),\\ \widehat\phi_3(\xi,\eta) &=\frac14(1+\xi)(1+\eta),& \widehat\phi_4(\xi,\eta) &=\frac14(1-\xi)(1+\eta). \end{aligned} \tag{6.54}\]

They satisfy the nodal interpolation property and sum to one. Each restriction to a reference edge is a one-dimensional linear function.

6.5.2.1 Isoparametric geometry and field interpolation

The geometry map is defined by

\[ \bx = \bx^e(\xi,\eta) = \sum_{a=1}^{4} \widehat\phi_a(\xi,\eta)\bx_a^e. \tag{6.55}\]

The same reference functions are used to interpolate the finite element field,

\[ u_h(\bx^e(\xi,\eta)) = \sum_{a=1}^{4} U_a^e\widehat\phi_a(\xi,\eta). \tag{6.56}\]

Using the same basis for the geometry and the field is the isoparametric construction introduced in Chapter 5.

Figure 6.19: Isoparametric mapping from the reference square to a four-node physical quadrilateral.

The physical basis is again obtained by composition,

\[ \phi_a^e(\bx) = \widehat\phi_a((\bx^e)^{-1}(\bx)). \]

For a general quadrilateral, the Jacobian matrix

\[ \bm J_e(\xi,\eta) = \begin{bmatrix} \partial x/\partial\xi&\partial x/\partial\eta\\ \partial y/\partial\xi&\partial y/\partial\eta \end{bmatrix} \tag{6.57}\]

varies across the element. Consequently, \(\bm J_e^{-T}\) and \(|\det\bm J_e|\) must be evaluated at every quadrature point. A valid, orientation-preserving element must satisfy

\[ \det\bm J_e(\xi,\eta)>0 \qquad \text{throughout }\widehat{\mathcal K}. \tag{6.58}\]

If the physical element is a parallelogram, the map reduces to an affine map and \(\bm J_e\) is constant. For a non-parallelogram, the map is bilinear and the mapped basis functions are generally not polynomials in the physical coordinates \((x,y)\).

Continuity across a shared edge follows by the same reasoning used for triangles, but it now involves the geometry map as well as the field. On the reference edge \(\eta=-1\), for example, \(\widehat\phi_3\) and \(\widehat\phi_4\) vanish, while \(\widehat\phi_1=(1-\xi)/2\) and \(\widehat\phi_2=(1+\xi)/2\). Substituting these restrictions into Equation 6.55 and Equation 6.56 shows that the map sends this reference edge linearly onto the straight segment from \(\bx_1^e\) to \(\bx_2^e\), and that the field varies linearly along that segment between \(U_1^e\) and \(U_2^e\). The same holds on every edge. If two conforming \(Q_1\) elements share a physical edge, both therefore describe the field on that edge as the linear function determined by the same two endpoint degrees of freedom, and the elementwise pieces join continuously. The finite element functions belong to \(H^1(\Om)\) for the same reason as in the \(P_1\) case.

Example 6.6 (Jacobian of a nonrectangular quadrilateral) Consider the physical vertices

\[ \bx_1=(0,0),\quad \bx_2=(2,0),\quad \bx_3=(5/2,1),\quad \bx_4=(0,1). \]

Substitution into Equation 6.55 gives

\[ x(\xi,\eta) = (1+\xi)\left(\frac98+\frac18\eta\right), \qquad y(\xi,\eta)=\frac{1+\eta}{2}. \]

Therefore,

\[ \bm J_e(\xi,\eta) = \begin{bmatrix} 9/8+\eta/8&(1+\xi)/8\\ 0&1/2 \end{bmatrix}, \qquad \det\bm J_e = \frac{9+\eta}{16}. \tag{6.59}\]

The determinant varies between \(1/2\) and \(5/8\) on the reference square and is positive everywhere. The mapping is valid, but geometric factors must be recomputed at each quadrature point.

A tensor-product Gauss rule is natural on the reference square. A \(2\times2\) rule uses the four combinations

\[ (\xi_r,\eta_r) \in \left\{-\frac1{\sqrt3},\frac1{\sqrt3}\right\} \times \left\{-\frac1{\sqrt3},\frac1{\sqrt3}\right\}, \tag{6.60}\]

with unit weights. At each point, the implementation evaluates the reference basis, \(\bm J_e\), \(\det\bm J_e\), the transformed gradients, and the integrand.

Exercise 6.13 (Properties and integration of a \(Q_1\) element)  

1. Verify the nodal interpolation and partition-of-unity properties of
2. Show that the restriction of the four basis functions to each reference edge is linear.
3. For the element in Example 6.6, evaluate \(\bx^e\), \(\bm J_e\), and \(\det\bm J_e\) at the four points in Equation 6.60.
4. Map the four Gauss points to the physical element.
5. Explain why one constant Jacobian cannot be used for this element.

Exercise 6.14 (An invalid bilinear map) Consider the four vertices

\[ \bx_1=(0,0),\qquad \bx_2=(2,0),\qquad \bx_3=(1/2,1/2),\qquad \bx_4=(0,2), \]

listed in the local order used in 1. The quadrilateral they define is not convex: its interior angle at \(\bx_3\) exceeds \(180^\circ\).

1. Show that at the reference vertex \((1,1)\), \(\partial\bx/\partial\xi=(\bx_3-\bx_4)/2\) and \(\partial\bx/\partial\eta=(\bx_3-\bx_2)/2\). Derive the corresponding expressions at the other three reference vertices.
2. Evaluate \(\det\bm J_e\) at the four reference vertices.
3. Identify the vertex at which the map ceases to be orientation-preserving, and explain why Equation 6.58 fails for this element.
4. Repeat Part 2 for the unit square with vertices listed in the crossed order \((0,0)\), \((1,0)\), \((0,1)\), \((1,1)\). Explain why this ordering is also invalid even though the four points form a square.

6.5.3 Tetrahedral elements

The three-dimensional analogue of a triangle is a tetrahedron. The reference tetrahedron is

\[ \widehat{\mathcal K} = \left\{(\xi,\eta,\zeta): \xi,\eta,\zeta\ge0, \ \xi+\eta+\zeta\le1 \right\}. \tag{6.61}\]

Its four vertices are the origin and the three Cartesian unit vectors. The \(P_1\) nodal basis is

\[ \widehat\phi_1=1-\xi-\eta-\zeta, \qquad \widehat\phi_2=\xi, \qquad \widehat\phi_3=\eta, \qquad \widehat\phi_4=\zeta. \tag{6.62}\]

In this subsection the reference-coordinate vector is

\[ \bm\xi = \begin{bmatrix}\xi&\eta&\zeta\end{bmatrix}^T. \]

For physical vertices \(\bx_1^e,\ldots,\bx_4^e\), the affine map is

\[ \bx^e(\bm\xi) = \bx_1^e + \bm J_e \begin{bmatrix}\xi\\\eta\\\zeta\end{bmatrix}, \qquad \bm J_e = \begin{bmatrix} \bx_2^e-\bx_1^e& \bx_3^e-\bx_1^e& \bx_4^e-\bx_1^e \end{bmatrix}. \tag{6.63}\]

The reference tetrahedron has volume \(1/6\), so

\[ |\mathcal K_e| = \frac{|\det\bm J_e|}{6}. \tag{6.64}\]

With a consistent orientation, the local vertex ordering is chosen so that \(\det\bm J_e>0\).

The volume and gradient transformations have exactly the same form as in two dimensions:

\[ \dx=|\det\bm J_e|\,\dd\bm\xi, \qquad \grad_{\bx}\phi_a^e = \bm J_e^{-T}\grad_{\bm\xi}\widehat\phi_a. \]

The centroid rule evaluates the integrand at

\[ \bm\xi_1 = \begin{bmatrix}1/4&1/4&1/4\end{bmatrix}^T \]

with weight \(w_1=1/6\). It integrates every polynomial of total degree at most one exactly on the reference tetrahedron.

Exercise 6.15 (A physical tetrahedron) Consider the tetrahedron with vertices

\[ (0,0,0),\qquad(2,0,0),\qquad(0,1,0),\qquad(0,0,3). \]

1. Form \(\bm J_e\) and compute the physical volume.
2. Compute the four physical basis gradients.
3. Verify that the gradients sum to zero.
4. For constant \(k\), write the \(4\times4\) element stiffness matrix in the form \(K_{ab}^e=k|\mathcal K_e|\grad\phi_b^e\cdot\grad\phi_a^e\).

Exercise 6.16 (A flux on a tetrahedral face) For the tetrahedron in Exercise 6.15, consider the face \(\mathcal F\) with vertices \((2,0,0)\), \((0,1,0)\), and \((0,0,3)\).

1. Parameterize \(\mathcal F\) affinely from the reference triangle, starting at \((2,0,0)\), and compute \(\bm g_1\), \(\bm g_2\), and \(\bm g_1\times\bm g_2\).
2. Use Equation 6.66 to compute the area of \(\mathcal F\).
3. Determine the outward unit normal and verify that it points away from the fourth vertex \((0,0,0)\).
4. For a constant normal datum \(q_0\) on this face, compute the three facet load entries associated with the \(P_1\) basis functions of the face vertices, and verify that they sum to \(q_0|\mathcal F|\). For heat conduction, verify that the physical outward heat-flow integral is \(-q_0|\mathcal F|\).

6.5.4 Hexahedral elements

A hexahedron is the three-dimensional tensor-product analogue of a quadrilateral. On the reference cube \([-1,1]^3\), associate each vertex with a sign triple \((\xi_a,\eta_a,\zeta_a)\in\{-1,1\}^3\). The eight trilinear basis functions can be written compactly as

\[ \widehat\phi_a(\xi,\eta,\zeta) = \frac18 (1+\xi_a\xi) (1+\eta_a\eta) (1+\zeta_a\zeta), \qquad a=1,\ldots,8. \tag{6.65}\]

One consistent numbering is \(1=(-1,-1,-1)\), \(2=(1,-1,-1)\), \(3=(1,1,-1)\), \(4=(-1,1,-1)\), \(5=(-1,-1,1)\), \(6=(1,-1,1)\), \(7=(1,1,1)\), and \(8=(-1,1,1)\). The first four nodes form the lower face and the last four the corresponding upper face. This is the numbering in Figure 6.20; a software library may use another order.

Figure 6.20: Local reference-node numbering for the tetrahedron and hexahedron. The tetrahedral nodes follow the four basis functions above; the cube follows the sign triples specified in the text.

Their span is

\[ \mathbb Q_1 = \operatorname{span} \{1,\xi,\eta,\zeta,\xi\eta,\xi\zeta,\eta\zeta,\xi\eta\zeta\}. \]

The geometry map \(\bx=\sum_a\widehat\phi_a\bx_a^e\) is trilinear, and a general physical hexahedron has a Jacobian that varies within the element. Volume integrals use

\[ \int_{\mathcal K_e}g(\bx)\,\dx = \int_{[-1,1]^3} g(\bx^e(\bm\xi))|\det\bm J_e(\bm\xi)|\,\dd\bm\xi. \]

A valid orientation-preserving hexahedral map must satisfy \(\det\bm J_e(\bm\xi)>0\) throughout the reference cube.

Tensor-product Gauss quadrature forms the three-dimensional quadrature points from one-dimensional Gauss points. The same sequence—map points, evaluate the Jacobian, transform gradients, and accumulate local arrays—therefore applies to triangles, quadrilaterals, tetrahedra, and hexahedra.

Exercise 6.17 (The reference \(Q_1\) hexahedron)  

1. Verify the nodal property of Equation 6.65 at the eight sign triples.
2. Verify that the eight functions sum to one.
3. How many points are produced by a \(2\times2\times2\) tensor-product Gauss rule?
4. Explain why the Jacobian is constant for an affine parallelepiped, including a rectangular box, but may vary for a general distorted hexahedron.

6.5.5 Boundary-face integration in three dimensions

The three-dimensional facet transformation is the analogue of the edge Jacobian Equation 6.41. Let a facet \(\mathcal F\) of a tetrahedral or hexahedral element be parameterized by a map \(\bx(\zeta_1,\zeta_2)\) from a reference facet \(\widehat{\mathcal F}\), which is the reference triangle for a tetrahedral face and the square \([-1,1]^2\) for a hexahedral face. The two vectors

\[ \bm g_1=\frac{\partial\bx}{\partial\zeta_1}, \qquad \bm g_2=\frac{\partial\bx}{\partial\zeta_2} \]

are tangent to the facet. The parallelogram they span has area \(\norm{\bm g_1\times\bm g_2}_2\), so the surface measure transforms as

\[ \ds = \norm{\bm g_1\times\bm g_2}_2\,\dd\zeta_1\dd\zeta_2 \tag{6.66}\]

(Belytschko et al. 2000, 4–44). A facet load entry then becomes

\[ \int_{\mathcal F}q\,\phi_a^e\,\ds = \int_{\widehat{\mathcal F}} q(\bx(\bm\zeta))\, \phi_a^e(\bx(\bm\zeta))\, \norm{\bm g_1\times\bm g_2}_2 \,\dd\zeta_1\dd\zeta_2, \qquad \bm\zeta=(\zeta_1,\zeta_2), \]

which is evaluated with a quadrature rule on \(\widehat{\mathcal F}\). The unit normal is \(\bn=\pm(\bm g_1\times\bm g_2)/\norm{\bm g_1\times\bm g_2}_2\), with the sign chosen so that \(\bn\) points out of the element. For a planar triangular face with vertices \(\bx_i\), \(\bx_j\), and \(\bx_k\), the affine parameterization \(\bx=\bx_i+\zeta_1(\bx_j-\bx_i)+\zeta_2(\bx_k-\bx_i)\) gives constant tangent vectors, and \(\norm{\bm g_1\times\bm g_2}_2=2|\mathcal F|\). For a face of a distorted hexahedron, the tangent vectors vary over the face and the facet Jacobian must be evaluated at every facet quadrature point.

6.6 Vector-valued approximation and linear elasticity

The scalar construction extends naturally to vector-valued fields. In two-dimensional elasticity, the displacement has two components,

\[ \bu(\bx) = \begin{bmatrix} u_1(\bx)\\ u_2(\bx) \end{bmatrix}, \]

so each nodal vertex of a displacement-based element carries two displacement degrees of freedom.

Let \(\mathcal V_h\) be the scalar finite element space. Its vector-valued counterpart is

\[ \mathcal V_h(\Om;\R^2) = \left\{ \bm w_h:\Om\to\R^2: w_{h1},w_{h2}\in\mathcal V_h \right\}. \tag{6.67}\]

Using the displacement boundary data from Equation 4.34, define

\[ \begin{aligned} \mathcal U_{h,\overline{\bu}_h} &=\{\bm w_h\in\mathcal V_h(\Om;\R^2): \bm w_h=\overline{\bu}_h\text{ on }\Gam_u\},\\ \mathcal V_{h,0} &=\{\bv_h\in\mathcal V_h(\Om;\R^2): \bv_h=\bm 0\text{ on }\Gam_u\}. \end{aligned} \tag{6.68}\]

Here \(\overline{\bu}_h\) is the finite element interpolation of the prescribed boundary displacement. The discrete elasticity problem is the restriction of Equation 4.36 to these spaces:

\[ \left. \begin{aligned} &\text{Find }\bu_h\in\mathcal U_{h,\overline{\bu}_h}\text{ such that}\\ &\int_{\Om} \beps(\bv_h):\mathbb C:\beps(\bu_h)\,\dx = \int_{\Om}\rho\bb\cdot\bv_h\,\dx + \int_{\Gam_t}\overline{\bt}\cdot\bv_h\,\ds \quad\forall\bv_h\in\mathcal V_{h,0}. \end{aligned} \right\} \qquad\text{Discrete variational problem} \tag{6.69}\]

Figure 6.21: A two-dimensional triangular displacement element carries two displacement degrees of freedom at each vertex.

Using the scalar basis functions \(\phi_i\), write

\[ \bu_h(\bx) = \sum_i \phi_i(\bx)\bm U_i, \qquad \bm U_i = \begin{bmatrix} U_{i1}\\ U_{i2} \end{bmatrix}. \tag{6.70}\]

On one three-node triangle,

\[ \bu_h|_{\mathcal K_e} = \sum_{a=1}^{3} \phi_a^e\bm U_a^e. \]

Collect the six element coefficients as

\[ \bU^e = \begin{bmatrix} U_{11}^e& U_{12}^e& U_{21}^e& U_{22}^e& U_{31}^e& U_{32}^e \end{bmatrix}^T. \]

The ordering is by vertex: the two components at vertex 1 are followed by the two components at vertex 2 and then those at vertex 3. Other orderings are possible, but the element degree-of-freedom map and every element matrix must use the same ordering.

Then

\[ \bu_h|_{\mathcal K_e} = \bN^e\bU^e, \]

with

\[ \bN^e = \begin{bmatrix} \phi_1^e&0&\phi_2^e&0&\phi_3^e&0\\ 0&\phi_1^e&0&\phi_2^e&0&\phi_3^e \end{bmatrix}. \tag{6.71}\]

For the two-dimensional reduction of Remark 4.4, use the in-plane engineering-strain vector, with \(\gamma_{12}=2\varepsilon_{12}\),

\[ \beps_{\mathrm V} = \begin{bmatrix} \varepsilon_{11}\\ \varepsilon_{22}\\ \gamma_{12} \end{bmatrix}, \]

the element strain is

\[ \beps_{\mathrm V} = \bB^e\bU^e, \]

where

\[ \bB^e = \begin{bmatrix} \phi_{1,x}^e&0& \phi_{2,x}^e&0& \phi_{3,x}^e&0\\ 0&\phi_{1,y}^e& 0&\phi_{2,y}^e& 0&\phi_{3,y}^e\\ \phi_{1,y}^e&\phi_{1,x}^e& \phi_{2,y}^e&\phi_{2,x}^e& \phi_{3,y}^e&\phi_{3,x}^e \end{bmatrix}. \tag{6.72}\]

For a \(P_1\) triangle, the basis gradients are constant, so \(\bB^e\) is constant on the element. This is the familiar constant-strain triangle.

Let \(\bm C_{\mathrm V}\) be the plane-stress or plane-strain constitutive matrix obtained by the appropriate reduction in Remark 4.4. The element stiffness matrix is

\[ \bK^e = \int_{\mathcal K_e} (\bB^e)^T \bm C_{\mathrm V} \bB^e \,\dx. \tag{6.73}\]

The body-force and traction contributions are

\[ \bF^e = \int_{\mathcal K_e} (\bN^e)^T\rho\bb\,\dx + \int_{\partial\mathcal K_e\cap\Gam_t} (\bN^e)^T\overline{\bt}\,\ds. \tag{6.74}\]

These two-dimensional formulas are written per unit out-of-plane thickness. If the model represents a body with constant thickness \(t\), both the domain and edge contributions are multiplied by \(t\).

The mathematical pattern is the same as for scalar diffusion: basis functions define the approximation, derivatives form the element strain or gradient operator, element integrals produce local arrays, and connectivity determines where those arrays enter the global system.

Example 6.7 (Strain in a constant-strain triangle) For the triangle in Figure 6.3, the basis gradients in Equation 6.12 give

\[ \bB^e = \begin{bmatrix} -1/2&0&1/2&0&0&0\\ 0&-1&0&0&0&1\\ -1&-1/2&0&1/2&1&0 \end{bmatrix}. \tag{6.75}\]

Consider the affine displacement field

\[ \bu(x,y) = \begin{bmatrix} 2x+y\\ -x+3y \end{bmatrix}. \]

Its nodal coefficient vector is

\[ \bU^e = \begin{bmatrix} 0&0&4&-2&1&3 \end{bmatrix}^T. \]

Multiplication by Equation 6.75 gives

\[ \beps_{\mathrm V} = \bB^e\bU^e = \begin{bmatrix} 2\\3\\0 \end{bmatrix}. \]

The first two entries are \(\varepsilon_{11}\) and \(\varepsilon_{22}\). The third is the engineering shear strain \(\gamma_{12}=u_{1,2}+u_{2,1}=1-1=0\). Because the element reproduces affine displacements, it reproduces their constant strains exactly.

Exercise 6.18 (Rigid motions and the constant-strain triangle) Use Equation 6.75 to verify the following statements.

1. If every node is given the same translation \((c_1,c_2)^T\), then \(\bB^e\bU^e=\bm 0\).
2. The infinitesimal rigid rotation \(\bu(x,y)=(-\theta y,\theta x)^T\) also gives \(\bB^e\bU^e=\bm 0\).
3. Explain why these results are required of a physically meaningful linear elastic element.

Exercise 6.19 (Stiffness matrix of a constant-strain triangle) Use the triangle of Example 6.7 with unit thickness. Take an isotropic material with Young’s modulus \(E\) and Poisson’s ratio \(\nu=0\). Then \(\lambda=0\) and \(\mu=E/2\), and the plane-stress and plane-strain constitutive matrices coincide:

\[ \bm C_{\mathrm V} = E \begin{bmatrix} 1&0&0\\ 0&1&0\\ 0&0&1/2 \end{bmatrix}. \]

1. Explain why Equation 6.73 reduces to \(\bK^e=|\mathcal K_e|(\bB^e)^T\bm C_{\mathrm V}\bB^e\) for this element.
2. Use Equation 6.75 to compute the \(6\times6\) matrix \(\bK^e\).
3. Verify that \(\bK^e\) is symmetric.
4. Verify that the nodal vectors of the two rigid translations and of the infinitesimal rotation in Exercise 6.18 lie in the nullspace of \(\bK^e\).
5. An unconstrained element matrix must have these three zero-energy modes. Explain why they are physical rigid motions rather than spurious modes, and why they are removed from the assembled problem only by displacement conditions that prevent rigid motion of the whole body.

Exercise 6.20 (Local and global displacement degrees of freedom) A two-dimensional triangular mesh has global vertices \(1,2,3,4\) and elements

\[ I^1=(1,2,3), \qquad I^2=(2,4,3). \]

Assume two displacement components per vertex.

1. Give one consistent global numbering of the eight displacement degrees of freedom.
2. Write the six global degree-of-freedom numbers associated with each element.
3. Identify which displacement degrees of freedom are shared by the two triangles.
4. Explain how the scalar vertex connectivity is expanded into a vector-valued degree-of-freedom map.

6.7 Chapter summary

• A multidimensional mesh partitions the domain into elements. Shared vertices, edges, and facets determine how local polynomial pieces and their degrees of freedom are connected.
• The scalar \(P_1\) triangle has three nodal basis functions with constant physical gradients. Neighboring triangles join continuously because their restrictions to a common facet are determined by the same nodal values.
• A global basis function is supported on the patch of elements incident to its node. This local support produces sparse global matrices and permits assembly from element stiffness matrices, domain load vectors, and boundary-facet load vectors.
• The reference-to-physical map transforms coordinates, measures, and gradients through its Jacobian matrix, determinant, and inverse transpose. Quadrature and boundary-facet integration can therefore be performed on fixed reference entities.
• The six-node \(P_2\) triangle adds edge-midpoint degrees of freedom. Its field interpolation may be combined with either an affine geometry map or a quadratic isoparametric map.
• Tetrahedra extend the simplex construction to three dimensions, while quadrilateral and hexahedral elements use tensor-product basis functions and generally nonconstant isoparametric mappings.
• Vector-valued elasticity interpolates each displacement component with the scalar basis. The strain-displacement matrix collects the required basis derivatives, and assembly and boundary conditions use the same connectivity rules as in the scalar problem.
• FEniCSx automates element evaluation, assembly, boundary-condition treatment, and solution, but its objects correspond directly to the mesh, spaces, variational forms, and algebraic system developed in the chapter.

Arnold, Douglas N. 2018. Lecture Notes on Numerical Analysis of Partial Differential Equations. University of Minnesota.

Belytschko, Ted, Wing Kam Liu, and Brian Moran. 2000. Nonlinear Finite Elements for Continua and Structures. Wiley.

Jog, C. S. 2011. Introduction to the Finite Element Method. Department of Mechanical Engineering, Indian Institute of Science, Bangalore. https://mecheng.iisc.ac.in/csjog/me257Notes/fem.pdf.

Larson, Mats G., and Fredrik Bengzon. 2010. The Finite Element Method: Theory, Implementation, and Practice. Springer.