# 10. Discretization

## 10.1. Weak Form for Large Deformation Problems

We begin by reviewing the field equations to be considered. The reference for this chapter is [[1]]. The problem to be solved is shown schematically in Fig. 10.1, in which we want to the compute finite deformation response of a body \(\Omega\) in its reference configuration.

Fig. 10.1 Large deformation initial/boundary value problem

Assuming that this time dependent configuration mapping is denoted by \(\varphi_t\), the following problem is solved for each time, \(t\), in the time interval of interest:

(10.1)\[\nabla \cdot \mathbf{T} + \mathbf{f} = \rho \mathbf{a} \text{ on } \varphi_t (\Omega),\]
(10.2)\[\varphi_t = \bar{\varphi_t} \text{ on } \varphi_t (\Gamma_u),\]

and

(10.3)\[\mathbf{t} = \bar{\mathbf{t}} \text{ on } \varphi_t (\Gamma_{\sigma}),\]

where all notations are as discussed in Section 4. In particular, \(\mathbf{a}\) is the material acceleration expressed in spatial coordinates, \(\mathbf{f}\) is the body force per unit (spatial) volume, and \(\mathbf{T}\) is the Cauchy stress tensor. The vector \(\mathbf{t}\) is the Cauchy traction vector, obtained via \(\mathbf{t} = \mathbf{T}\mathbf{n}\), where \(\mathbf{n}\) is the outward unit normal to the spatial surface \(\varphi_t (\Gamma_{\sigma})\).

The problem is also subject to initial conditions of the form

(10.4)\[\varphi(\mathbf{X},0) = \varphi_0 (\mathbf{X}) \text{ on } \Omega,\]

and

(10.5)\[\frac{\partial \varphi}{\partial t} (\mathbf{X},0) = \mathbf{X}_0 (\mathbf{X}) \text{ on } \Omega.\]

Recall that (10.1) through (10.3) are written in the so-called spatial configuration, but we still consider ourselves working in a Lagrangian framework where all quantities are ultimately indexed to material points through the mapping \(\mathbf{x} = \varphi_t(\mathbf{X})\) (see Lagrangian and Eulerian Descriptions in Section 4).

A prerequisite of the finite element method is that a weak, or variational, form of the above field equations be available for discretization. This can be obtained following the general procedure outlined for linear problems in Section 3 by considering weighting functions \(\varphi^*\) defined over \(\Omega\) which satisfy the following condition:

(10.6)\[\varphi^* = \mathbf{0}\text{ on } \Gamma_u,\]

where we also assume that all \(\varphi^*\) are sufficiently smooth so that any desired partial derivatives can be computed. In treating large deformation problems, it is useful to consider spatial forms of the functions \(\varphi^*\) obtained by composition with the (unknown) mapping \(\varphi^{-1}_t\). We denote these spatial variations by \(\mathbf{w}\) and note that they may be obtained via

(10.7)\[\mathbf{w}(\mathbf{x}) = \varphi^* \left( \varphi^{-1}_t(\mathbf{x}) \right)\]

for any \(\mathbf{x} \in \varphi_t(\mathbf{X})\). (10.6) means

(10.8)\[\mathbf{w} = \mathbf{0} \text{ on } \varphi_t (\Gamma_u ).\]

Assuming the configuration mapping \(\varphi_t\) is smooth, all required partial derivatives of \(\mathbf{w}\) can be computed.

With these definitions, the development in Section 3 can be reproduced in the current context to provide the following spatial representation of the variational form for large deformations:

Given the boundary conditions \(\bar{\mathbf{t}}\) on \(\varphi_t(\Gamma_{\sigma})\), \(\bar{\varphi_t}\) on \(\varphi_t(\Gamma_{u})\), the initial conditions \(\varphi_0\) and \(\mathbf{V}_0\) on \(\Omega\), and the distributed body for \(\mathbf{f}\) on \(\varphi_t (\Omega)\), find \(\varphi_t \in S_t\) for each time \(t \in (0,T)\) such that:

(10.9)\[\int_{\varphi_t (\Omega)} \rho \mathbf{w} \cdot \mathbf{a} \mathrm{d}v + \int_{\varphi_t (\Omega)} \nabla \mathbf{w}:\mathbf{T} \mathrm{d}v = \int_{\varphi_t (\Omega)} \mathbf{w} \cdot \mathbf{f} \mathrm{d}v + \int_{\varphi_t (\Gamma_{\sigma})} \mathbf{w} \cdot \bar{\mathbf{t}} \mathrm{d}a\]

for all admissible \(\mathbf{w}\), where \(S_t\) is defined as

(10.10)\[S_t = \left\{ \varphi_t \middle| \varphi_t = \bar{\varphi}(t) \text{ on } \Gamma_u \text{ , } \varphi_t \text{ is smooth} \right\}\]

and where admissible \(\mathbf{w}\) are related in a one-to-one manner via (10.7) to the material variations \(\varphi^* \in W\) with the definition of \(W\) being

(10.11)\[W = \left\{ \varphi^* \middle| \varphi^* = \mathbf{0} \text{ on } \Gamma_u \text{ , } \varphi^* \text{ is smooth} \right\} .\]

Note that in contrast to the previous development, the constitutive relation governing \(\mathbf{T}\) is left unspecified, but it can in general be subject to both geometric and material nonlinearities. Furthermore, it should be implied that geometric nonlinearities include consideration of large deformation kinematics discussed in Section 5, Section 6, and Section 9. The notation \(\mathbf{a}\) for the acceleration is to be understood as the material acceleration as defined by (6.4).

In addition, the solution \(\varphi\) is subject to the following conditions at \(t=0\):

(10.12)\[\int_{\Omega} \varphi^* \cdot (\varphi |_{t=0} - \varphi_0) \mathrm{d} \Omega = 0\]

and

(10.13)\[\int_{\Omega} \varphi^* \cdot \left( \left.\frac{\partial \varphi}{\partial t}\right|_{t=0} - \mathbf{V}_0 \right) \mathrm{d} \Omega = 0 ,\]

both of which must hold for all \(\varphi^* \in W\).

## 10.2. Finite Element Discretization

The process of numerically approximating a continuous problem is generically called discretization. In the finite element method, the entity discretized is a weak form (alternatively, variational equation). The variational form to be considered here is that just summarized in the previous section. We now refer to Fig. 10.2 which gives the general notation to be used in the description of the discretization process.

Fig. 10.2 General notation for finite element discretization of the reference domain.

Referring to Fig. 10.2, the reference domain \(\Omega\) is subdivided into a number of element subdomains \(\Omega^e\). The superscript \(e\) is an index to a specific element, running between \(1\) and the total number of elements in the discretization, \(n_{el}\), of the domain \(\Omega\). We assume in the figure and throughout the ensuing discussion that \(\Omega\) is a subset of \(\mathbb{R}^3\).

Note that a number of nodal points are indicated by the dots in Fig. 10.2. We assume that all degrees of freedom in the discrete system to be proposed will be associated with these nodes. These nodes may lay at corners, edges, and in interiors of the elements with which they are associated. A key feature of the finite element method will be that a specific element can be completely characterized by the coordinates and degrees-of-freedom associated with the nodes attached to it. In the following we will index the nodes with uppercase letters \(A\), \(B\), etc. having values running between \(1\) and \(n_{np}\), the total number of nodal points in the problem discretization.

## 10.3. Galerkin Finite Element Methods

The essence of any finite element method lies in the discretization of the variational form. This discretization process involves approximation of a typical member of both the solution space \(S_t\) and the weighting space \(W\). These approximations are typically expressed as an expansion in terms of prescribed shape or interpolation functions, usually associated with specific nodal points in the mesh. Since the number of nodal points is obviously finite, the expansion is likewise finite, giving rise to the concept of a finite-dimensional approximation of the space.

Roughly speaking, the idea of discretization is as follows. We know from earlier chapters that if the variational equation is enforced considering all \(\varphi_t \in S_t\) and \(\varphi^* \in W\) as mandated by its definition, then the solution of the weak form is completely equivalent to that of the strong form (i.e., the governing partial differential equation with boundary/initial conditions). This fact results because of the arbitrary nature of \(\varphi^*\) and the very general definitions for \(S_t\) and \(W\). If we restrict our attention only to some subset of the above spaces, we inherently incur some error with the solution of our approximated weak form in that it no longer is identical to the solution of the strong form. If our choice for the type of shape functions to be used is reasonable, however, we can represent the full solution and weighting spaces with arbitrary closeness by increasing the number of nodal points and/or the degree of polynomial approximation utilized in the interpolation functions. In the limit of such refinement, we should expect recovery of the exact solution (i.e., convergence).

We represent the shape function associated with node \(A\) as \(N_A\) and assume it to be as follows:

(10.14)\[N_A : \bar{\Omega} \rightarrow \mathbb{R} .\]

Given a time, \(t\), the finite dimensional counterpart of \(\varphi_t\) will be denoted as \(\varphi^h_t\) and is expressed in terms of the shape functions as

(10.15)\[\varphi^h_t = \sum_{B=1}^{n_{np}} N_B \mathbf{d}_B (t) ,\]

where \(\mathbf{d}_B (t)\) is a 3-vector containing the unknown displacements of nodal point \(B\) at time \(t\). Given a prescribed set of nodal shape functions \({N_B}, \text{ } B=1,\dots,n_{np}\), the finite dimensional solution space \(S^h_t\) is defined as the collection of all such \(\varphi^h_t\):

(10.16)\[S^h_t = \left\{ \varphi^h_t = \sum_{B=1}^{n_{np}} N_B \mathbf{d}_B (t) \middle| \varphi^h_t \approx \tilde{\varphi}_t (\mathbf{X} \text{ for all } \mathbf{X} \in \Gamma_u \right\} .\]

In other words, we require members of the discrete solution space to (approximately) satisfy the displacement boundary condition on \(\Gamma_u\). The approximation comes about because, in general, we only force \(\varphi^h_t\) to interpolate the nodal values of \(\bar{\varphi}_t\) on \(\Gamma_u\) with the \(N_B\) serving as the interpolation functions. Note that \(\Gamma_u\) itself is typically geometrically approximated by the finite element discretization, also contributing to the approximation.

This notationally defines the discretization procedure for \(\varphi^h_t\). It still remains, however, to approximate the weighting space. The (Bubnov-) Galerkin finite element method is characterized by utilizing the same shape functions to approximate \(W\) as were used to approximate \(S_t\). Accordingly, we define a member of this space, \(\left(\varphi^{*}\right)^h\), via

(10.17)\[\left(\varphi^{*}\right)^h = \sum_{A=1}^{n_{np}} N_A \mathbf{c}_A ,\]

where \(\mathbf{c}_A\) are 3-vectors of nodal constants. We can then express the finite dimensional weighting space \(W^h\) via

(10.18)\[W^h = \left\{ \left(\varphi^{*}\right)^h = \sum_{A=1}^{n_{np}} N_A \mathbf{c}_A \middle| \left(\varphi^{*}\right)^h = 0 \text{ for all } \mathbf{X} \in \Gamma_u \right\} .\]

Analogous to the situation for \(S^h_t\), (10.18) features a discrete version of the boundary condition on \(\Gamma_u\). In other words, \(W^h\) consists of all functions of the form (10.17) resulting in satisfaction of this condition. Note that the only restriction on \(\mathbf{c}_A\) is that they result in satisfaction of the homogeneous boundary condition on \(\Gamma_u\).

With these ideas in hand, the approximate Galerkin solution to the initial/boundary value problem takes the form described below.

Given the boundary conditions \(\bar{\mathbf{t}}\) on \(\varphi^h_t (\Gamma_{\sigma})\), \(\bar{\varphi}_t\) on \(\varphi^h_t(\Gamma_u)\), the initial conditions \(\varphi_0\) and \(\mathbf{V}_0\) on \(\Omega\), and the distributed body force \(\mathbf{f}\) on \(\varphi^h_t(\Omega)\), find \(\varphi^h_t \in S^h_t\) for each time \(t \in (0,T)\) such that:

(10.19)\[\int_{\varphi^h_t(\Omega)} \rho \mathbf{w}^h \cdot \mathbf{a}^h \mathrm{d}v + \int_{\varphi^h_t(\Omega)} \nabla \mathbf{w}^h:\mathbf{T}^h \mathrm{d}v = \int_{\varphi^h_t(\Omega)} \mathbf{w}^h \cdot \mathbf{f} \mathrm{d}v + \int_{\varphi^h_t(\Gamma_{\sigma})} \mathbf{w}^h \cdot \bar{\mathbf{t}} \mathrm{d}a\]

for all admissible \(\mathbf{w}^h\), where \(S_t\) is defined in (10.16) and where admissible \(\mathbf{w}^h\) are related to the material variations \(\left(\varphi^{*}\right)^h \in W^h\) via

(10.20)\[\mathbf{w}^h (\mathbf{x}) = \left(\varphi^{*}\right)^h \in \left(\varphi^h_t\right)^{-1} (\mathbf{x}).\]

In (10.19), \(\mathbf{T}^h\) refers to the Cauchy stress field computed from the discrete mapping \(\varphi^h_t\) through the constitutive relations, whereas \(\mathbf{a}^h\) is the discrete material acceleration.

The initial conditions are ordinarily simplified in the discrete case to read

(10.21)\[\mathbf{d}_B (0) = \bar{\varphi}_0 (\mathbf{X}_B)\]

and

(10.22)\[\dot{\mathbf{d}}_B (0) = \mathbf{V}_0 (\mathbf{X}_B),\]

both of which must hold for all nodes \(B=1,\dots,n_{np}\), where \(\mathbf{X}_B\) are the reference coordinates of the node in question.

section{Notation for Discrete Problem}

In preparation for generating vector/matrix equations for the discrete system, it will be helpful to be explicit with our notation. We therefore express the nodal vectors \(\mathbf{c}_A\) and \(\mathbf{d}_B\) in terms of their components via

(10.23)\[\mathbf{c}_A = \left\{ c_{iA} \right\}, \text{ } i=1,2,3\]

and

(10.24)\[\mathbf{d}_B = \left\{ d_{jB} \right\}, \text{ } j=1,2,3.\]

Note that indices \(i\) and \(j\) are spatial indices, in general. It is useful in generating matrix equations to have indices referring not to nodes \(A\) and \(B\) or spatial directions \(i\) and \(j\), but rather to degree of freedom numbers in the problem. Thus, we define for notational convenience the concept of an \(ID\) array that is set up as follows:

(10.25)\[ID(i,A) = P \text{ (global degree of freedom number)}.\]

In other words, the \(ID\) array takes the spatial direction index and nodal point number as arguments and assigns a global degree of freedom number to the corresponding unknown. For three-dimensional deformation problems, the number of degrees of freedom \(n_{dof}\) is

(10.26)\[n_{dof} = 3 \times n_{np} .\]

With this notation, the equation numbers \(P\) and \(Q\) corresponding to the degrees of freedom are defined as

(10.27)\[P = ID(i,A)\]

and

(10.28)\[Q = ID(j,B) .\]

## 10.4. Discrete Equations

We now generate the discrete equations by substitution of (10.15) and (10.17) into (10.19), causing the variational equation to read

(10.29)\[\begin{split}\begin{aligned} \int_{\varphi^h_t(\Omega)} \rho \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) c_A \right) \cdot \left( \sum_{B=1}^{n_{np}} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \ddot{\mathbf{d}}_B (t) \right) \mathrm{d}v \\ + \int_{\varphi^h_t(\Omega)} \left( \sum_{A=1}^{n_{np}} \nabla N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) \otimes c_A \right) : \mathbf{T}^h \mathrm{d}v \\ + \int_{\varphi^h_t(\Omega)} \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) c_A \right) \cdot \mathbf{f} \mathrm{d}v + \int_{\varphi^h_t(\Gamma_{\sigma})} \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) c_A \right) \cdot \bar{\mathbf{t}} \mathrm{d}a \end{aligned}\end{split}\]

where we note in particular that \(\mathbf{T}^h\) is a function of $varphi^h_t = sum_{B=1}^{n_{np}} N_B mathbf{d}_B (t)$ through the strain-displacement relations (nonlinear, in general) and the constitutive law (as yet unspecified and perhaps likewise nonlinear).

The inertial term in (10.29) can be expanded as

(10.30)\[\begin{split}\begin{aligned} \int_{\varphi^h_t(\Omega)} \rho \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) c_A \right) \cdot \left( \sum_{B=1}^{n_{np}} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \ddot{\mathbf{d}}_B (t) \right) \mathrm{d}v \\ = \sum_{A=1}^{n_{np}} \sum_{i=1}^{3} c_{iA} \int_{\varphi^h_t(\Omega)} \left( \rho N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) c_{iA} \left( \sum_{B=1}^{n_{np}} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \ddot{d}_{iB} \right) \mathrm{d}v \right) \\ = \sum_{A=1}^{n_{np}} \sum_{i=1}^{3} c_{iA} \left[ \sum_{B=1}^{n_{np}} \sum_{i=1}^{3} \int_{\varphi^h_t(\Omega)} \rho N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v \ddot{d}_{iB} \right] \\ = \sum_{P=1}^{n_{dof}} c_P \left( \sum_{Q=1}^{n_{dof}} M_{PQ} \ddot{d}_Q \right) \end{aligned}\end{split}\]

where \(M_{PQ}\) is defined as

(10.31)\[M_{PQ} = \int_{\varphi^h_t(\Omega)} \rho N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v.\]

The second term of (10.29) can be simplified via

(10.32)\[\begin{split}\begin{aligned} \int_{\varphi^h(\Omega)} \left( \sum_{A=1}^{n_{np}} \nabla N_A \left( \varphi^{-1}_t(\mathbf{x}) \right) \otimes \mathbf{c}_A \right) : \mathbf{T}^h \mathrm{d}v \\ = \int_{\varphi^h(\Omega)} \left( \sum_{A=1}^{n_{np}} \sum_{i=1}^{3} \sum_{j=1}^{3} N_{A,j} \left( \varphi^{-1}_t (\mathbf{x}) \right) c_{iA} T^{h}_{ij} \right) \mathrm{d}v = \sum_{P=1}^{n_{dof}} c_P F^{\mathrm{int}}_P \end{aligned}\end{split}\]

where

(10.33)\[F^{\mathrm{int}}_{P} = \int_{\varphi^h_t(\Omega)} \left[ \sum_{j=1}^3 N_{A,j} \left( \varphi^{-1}_t (\mathbf{x}) \right) T_{ij} N_B \right] \mathrm{d}v\]

Finally, the last two terms of (10.29) can be treated as

(10.34)\[\int_{\varphi^h_t(\Omega)} \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathbf{c}_A \right) \cdot \mathbf{f} \mathrm{d}v + \int_{\varphi^h_t(\Gamma_{\sigma})} \left( \sum_{A=1}^{n_{np}} N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathbf{c}_A \right) \cdot \bar{\mathbf{t}} \mathrm{d}a = \sum_{P=1}^{n_{dof}} c_P F^{\mathrm{ext}}_P\]

where

(10.35)\[F^{\mathrm{ext}}_{P} = \int_{\varphi^h_t(\Omega)} N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) f_i \mathrm{d}v + \int_{\varphi^h_t(\Gamma_{\sigma})} N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \cdot \bar{t}_i \mathrm{d}a .\]

## 10.5. Generation of Vector/Matrix Equations

We now define the following vectors and matrices of global variables, all with dimension \(n_{dof}\):

(10.36)\[\begin{split}\begin{aligned} \mathbf{c} &= \left\{ c_P \right\} \\ \mathbf{d}(t) &= \left\{ d_Q(t) \right\} \\ \mathbf{F}^{\mathrm{int}} (\mathbf{d}(t)) &= \left\{ F^{\mathrm{int}}_P \right\} \\ \mathbf{F}^{\mathrm{ext}} &= \left\{ F^{\mathrm{ext}}_P \right\} \\ \mathbf{M} &= \left[ M_{PQ} \right] \\ \end{aligned}\end{split}\]

The results of (10.30)–(10.35) can now be summarized as

(10.37)\[\mathbf{c}^T \left[ \mathbf{M} \ddot{\mathbf{d}}(t) + \mathbf{F}^{\mathrm{int}}(\mathbf{d}(t)) - \mathbf{F}^{\mathrm{ext}} \right] = 0 ,\]

which must hold for all \(n_{dof}\)-vectors \(\mathbf{c}\) that result in satisfaction of the homogeneous boundary condition imposed on \(W\) (i.e., (10.18)).

Finally we observe that not all of the members of \(\mathbf{d}(t)\) are unknown; for nodes lying on \(\Gamma_u\), these degrees of freedom are prescribed. Furthermore, the corresponding entries of \(\mathbf{c}\) at these nodes are typically taken to be zero, so that the aforementioned condition on \(W^h\) is obeyed. Since the remainder of the vector \(\mathbf{c}\) is arbitrary, it must be the case that the elements of the bracketed term in (10.37) corresponding to free degrees of freedom must be identically zero, so that (10.37) will hold for arbitrary combinations of the \(c_P\). Thus we can write the nonlinear equation that expresses the discrete equations of motion:

(10.38)\[\mathbf{M} \ddot{\mathbf{d}}(t) + \mathbf{F}^{\mathrm{int}}(\mathbf{d}(t)) = \mathbf{F}^{\mathrm{ext}} .\]

Here we employ a slight abuse of notation because we have asserted in (10.36) that all vectors and matrices have dimension \(n_{dof}\), yet we only enforce (10.38) for free degrees of freedom. Denoting the number of free degrees of freedom as \(n_{eq}\), on can account for this difference in practice by calculating the vector and matrix entries for all degrees of freedom and then merely disregarding the \(n_{dof}-n_{eq}\) equations corresponding to the prescribed degrees of freedom. The members of \(\mathbf{d}(t)\) that are prescribed do need to be retained in its definition, however, since they enter into both terms on the left-hand side of (10.38). It should simply be remembered that only \(n_{eq}\) members of \(\mathbf{d}(t)\) are, in fact, unknown. We will have an opportunity to visit the general topic of constraint enforcement in greater detail when discussing solutions to these nonlinear equations (see Section 13).

## 10.6. Localization and Assembly

The description to this point is mostly a matter of mathematical manipulation with little insight gained into the character of the interpolation functions, \(N_A\). In fact, the basic nature of these interpolation functions distinguishes the finite element method from other variational solution techniques.

The detail of shape function construction will be discussed in Section 14 in the context of element programming. However it is useful to discuss here the basic character of finite element approximation functions to give general insight into the structure of the method. We refer to Fig. 10.3 which depicts a node \(A\) in \(\Omega\), along with the elements attached to it. A basic starting point for the development of a finite element method is as follows: the shape function associated with Node \(A\), \(N_A\), is only nonzero in that sub-portion of \(\Omega\) encompassed by the elements associated with Node \(A\) and is zero everywhere else in \(\Omega\).

This property of the shape functions is crucial to the modular character of the finite element method. Shape functions \(N_A\) having this property are said to possess local support.

Fig. 10.3 Local support of finite element interpolation functions. The region of support for \(N_A\) shown as shaded.

To gain insight into the effect of this property, we examine the expression given in (10.31) for an element of the mass matrix \(M_{PQ}\). We note in particular that the integrand of (10.31) will be nonzero if both nodes \(A\) and \(B\) share a common element in the mesh. Otherwise \(M_{PQ}\) must be zero. If we fix our attention on a given Node \(A\) in the mesh, we can conclude that very few Nodes \(B\) will produce nonzero column entries in \(\mathbf{M}\). This matrix is therefore sparse, and it would be a tremendous waste of time to compute \(\mathbf{M}\) by looping over all the possible combinations of node numbers and spatial indices without regard to elements and the node numbers attached to them.

Instead the global matrices and vectors needed in the solution of (10.38) are more typically computed using two important concepts: localization and assembly. Still considering the matrix \(\mathbf{M}\) as an example, we note that by the elementary properties of integration, we can write

(10.39)\[\begin{split}\begin{aligned} M_{PQ} & = \int_{\varphi^h_t(\Omega)} \rho N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v \\ & = \sum_{e=1}^{n_{el}} \int_{\varphi^h_t(\Omega^e)} \rho N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v \\ & = \sum_{e=1}^{n_{el}} M^e_{PQ} , \end{aligned}\end{split}\]

where

(10.40)\[M^e_{PQ} = \int_{\varphi^h_t(\Omega^e)} \rho N_A \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_B \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v.\]

Thus the global mass matrix can be computed as the sum of a number of element mass matrices. This fact in itself is not especially useful because each of the \(\mathbf{M}^e\) is extremely sparse, even more so than \(\mathbf{M}\). In fact, the only entries of \(\mathbf{M}^e\) that will be nonzero will be those for which both \(P\) and \(Q\) are degrees of freedom associated with element \(e\).

This fact can be exploited by defining another local element matrix \(\mathbf{m}^e\) containing only degrees of freedom associated with that element. We introduce element degrees of freedom indices \(p\) and \(q\), as indicated in Fig. 10.4. Assuming that \(p\) and \(q\) can take on values between \(1\) and \(n_{edof}\), where \(n_{edof}\) is the number of degrees of freedom associated with the element, an \(n_{edof} \times n_{edof}\) matrix \(\mathbf{m}\) is constructed as

(10.41)\[\mathbf{m}^e = \left[ m^e_{pq} \right] .\]

Fig. 10.4 Element (local) degrees of freedom for a sample finite element.

The \(m^e_{pq}\) can be specified by introducing the concept of a local node number \(a\) or \(b\) as shown in Fig. 10.4. With these definitions we can write

(10.42)\[m^e_{pq} = \int_{\varphi^h_t(\Omega^e)} \rho N_a \left( \varphi^{-1}_t (\mathbf{x}) \right) \delta_{ij} N_b \left( \varphi^{-1}_t (\mathbf{x}) \right) \mathrm{d}v\]

where a sample relationship between indices \(i\), \(a\), and \(p\) appropriate for the element at hand might be

(10.43)\[p = (a-1) \times 2 + i\]

(similarly for \(j\), \(b\), and \(q\)). The notation \(N_a\) simply refers to the shape function associated with local Node \(a\). By definition it is the restriction of the global interpolation function \(N_A\) to the element domain.

Calculation of the local element entities, such as \(\mathbf{m}^e\), turns out to be highly modular procedure whose form remains essentially unchanged for any element in a mesh. Detailed discussion of this calculation is deferred until Section 14.

Let us suppose for a moment, however, that we have a procedure in hand for calculating this matrix. We might then propose the following procedure for calculating the global mass matrix \(\mathbf{M}\) and internal force vector \(\mathbf{F}^{\mathrm{int}}\):

- Zero out \(\mathbf{M}\), \(\mathbf{F}^{\mathrm{int}}\).

- For each element \(e\), \(e = 1,\dots, n_{el}\) :

  - Prepare local data necessary for element calculations - e.g., \(\mathbf{X}^e\) (\(n_{edof}\) - vector of element nodal coordinates), \(\mathbf{d}^e\) (\(n_{edof}\)-vector of element nodal configuration mappings), etc.

  - Calculate element internal force vector \(\mathbf{f}^{\mathrm{int},e} = \left\{f^{\mathrm{int},e}_p \right\}\) and element mass matrix \(\mathbf{m}^e = \left[ m^e_{pq} \right]\) via

    (10.44)\[f^{\mathrm{int},e}_p = \int_{\varphi^h_t(\Omega^e)} \left[ \sum_{j=1}^3 N_{a,j} \left( \varphi^{-1}_t(\mathbf{x}) \right) T^h_{ij} \right] \mathrm{d}v\]

    and (10.42).

  - Assemble the element internal force vector and element mass matrix into their global counterparts by performing the following calculations for all local degrees of freedom \(p\) and \(q\):

    (10.45)\[M_{PQ} = M_{PQ} + m^e_{pq}\]

    and

    (10.46)\[F^{\mathrm{int}}_{P} = F^{\mathrm{int}}_{P} + f^{\mathrm{int},e}_p ,\]

    where local degrees of freedom are related to global degrees of freedom via the LM array, defined so that

    (10.47)\[P=LM(p,e)\]

    and

    (10.48)\[Q=LM(q,e) .\]

Step 2a) above is referred to as localization; given a particular element, \(e\), it extracts the local information from the global arrays necessary for element level calculations. Step 2b) consists of element level calculations; these calculations will be discussed in detail in Section 14. Step 2c) is the process known as assembly and takes the data produced by the element level calculations and assembles them in the proper locations of the global arrays.

We can thus now summarize the effect of localization and assembly in a finite element architecture. They act as pre- and post-processors to the element-level calculations, enabling the entities needed for global equilibrium calculations to be computed in a modular manner as summation of element contributions. Of course, the effectiveness of this procedure, as well as the convergence behavior of the numerical method in general, depends crucially on the interpolation functions chosen and their definitions in terms of elements. We defer this topic for now and concentrate in the coming chapters on the classes of problems and global equation-solving strategies to be utilized.

[1]

J.E. Marsden and T.J.R. Hughes. Mathematical Foundations of Elasticity. Dover, 1983.