22. Finite element system assembling

As a first step, we assume there are no Dirichlet boundary conditions. The finite element problem is

\[\begin{equation*} \mbox{Find } u_h \in V_{\cal T} \mbox{ such that } : A (u_h, v_h) = f(v_h) \qquad \forall \, v_h \in V_{\cal T} \end{equation*}\]

The nodal basis and the dual functionals provides the one to one relation between \({\mathbb R}^N\) and \(V_{\cal T}\):

\[ {\mathbb R}^N \ni \underline u \leftrightarrow u_h \in V_{\cal T} \qquad \mbox{with} \qquad u_h = \sum_{i=1}^N \varphi_i \underline u_i \qquad \mbox{and} \qquad \underline u_i = \psi_i (u_h). \]

Using the nodal basis expansion of \(u_h\), and testing only with the set of basis functions, one has

\[ A( \sum_{i=1}^N u_i \varphi_i, \varphi_j) = f(\varphi_j) \qquad \forall \, j = 1 \ldots N \]

With

\[ A_{ji} = A(\varphi_i, \varphi_j) \qquad \mbox{and} \qquad \underline f_j = f(\varphi_j), \]

one obtains the linear system of equations

\[ A \underline u = \underline f \]

The preferred way to compute the matrix \(A\) and vector \(f\) is a sum over element contributions. The restrictions of the bilinear and linear form to the elements are

\[ A_T(u,v) = \int_T \lambda \nabla u \cdot \nabla v \, dx + \int_{\partial \Omega \cap T} \alpha u v \, ds \]

and

\[ f_T(v) = \int_T f v \, dx + \int_{\partial \Omega \cap T} g v \, ds \]

Then

\[ A(u,v) = \sum_{T \in \cal T} A_T(u,v) \qquad f(v) = \sum_{T \in \cal T} f_T(v) \]

On each element, one defines the \(N_T \times N_T\) element matrix and element vector in terms of the local basis on \(T\):

\[ A_{T,\alpha \beta} = A_T (\varphi^T_\beta, \varphi^T_\alpha) \qquad \qquad \underline{f}_{T,\alpha} = \underline{f}_T (\varphi^T_\alpha) \]

Then, the global matrix and the global vector are

\[ A = \sum_{T \in \cal T} C_T A_T C_T^t \]

and

\[ \underline f = \sum_{T \in \cal T} C_T \underline f_T \]

Namely,

\[\begin{align*} \underline f_i & = f(\varphi_i) = \sum_{T \in \cal T} f_T (\varphi_i|_T) = \sum_{T \in \cal T} f_T (\sum_\alpha C_{T,i\alpha} \varphi_T^\alpha) \\ & = \sum_{T \in \cal T} \sum_{\alpha} C_{T,i\alpha} f_T (\varphi_T^\alpha) = \sum_{T \in \cal T} \sum_{\alpha} C_{T,i\alpha} \underline f_\alpha \end{align*}\]

and

\[\begin{align*} A_{ji} & = \sum_{T \in \cal T} A(\varphi_i|_T, \varphi_j|_T) = \sum_{T \in \cal T} A (\sum_\alpha C_{T,i\alpha} \varphi_T^\alpha, \sum_\beta C_{T,j\beta} \varphi_T^\beta) \\ & = \sum_{T \in \cal T} \sum_\alpha \sum_\beta C_{T,i\alpha} A_{T,\alpha \beta} C_{T,j\beta} \end{align*}\]

22.1. Computation of element-vectors and element-matrices

The element-vectors and element-matrices are usually computed by the pull-back to a reference element \(\widehat T\). We assume to have equivalent elements on \(T\) and \(\widehat T\). The mapping, the Jacobi-matrix \(F\) and the Jacobi-determinant \(J\) are:

\[ T = \Phi(\widehat T) \qquad F = \Phi^\prime \qquad J = \operatorname{det} F \]

To compute the entries of the element-vector we use the multi-dimensional substitution rule to map the integral over \(T\) to an integral over the reference element.

\[\begin{align*} f_{T,\alpha} & = \int_T f(x) \varphi_{T,\alpha}(x) \, dx \\ &= \int_{\widehat T} f (\Phi(\hat x)) \hat \varphi_{\widehat T, \alpha} (\hat x) \, |J(\hat x)| \, d \hat x \end{align*}\]

The integrands are usually smooth functions. Thus, numerical integration rules with integration points \(\{ x_1, \ldots x_m \}\) and integration weights \(\{ \omega_1, \ldots \omega_m \}\) can be applied:

\[ \int_{\widehat T} g(\hat x) d \hat x \approx \sum_{k=1}^m \omega_k g(x_k) \]

The order of an integration rule is the largest number such that all polynomials up to this order are integrated exact.

We write \(\widehat \varphi_{\widehat T}\) for the vector of \(N_T\) basis functions. Then, up to accuracy of numerical integration, we obtain the element vector \(\underline f_T \in {\mathbb R}^{N_T}\)

\[ \underline f_T = \sum_{k=1}^m \omega_k | J(x_k) | (f \circ \Phi_T) (x_k) \; \widehat \varphi_{\widehat T} (x_k) \]

For the pull-back of gradients we use that

\[ \varphi_\alpha (\Phi_T (\hat x)) = \hat \varphi_\alpha (\hat x), \]

apply the chain rule

\[ \varphi_\alpha^\prime \, \Phi_T^\prime = \hat \varphi_\alpha^\prime, \]

multiply with the inverse of the Jacobi-matrix, write the derivative as a column vector \(\nabla \varphi\) to obtain

\[ \nabla \varphi_\alpha (x) = F^{-t} \nabla \widehat \varphi_\alpha (\hat x) \]

for \(x = \Phi_T(\hat x)\), and the inverse of the transpose \(F^{-t}\).

The transformation of gradients allows to compute the element-matrix by integration on the reference element:

\[\begin{align*} A_{T,\alpha \beta} &= \int_T \lambda(x) \nabla \varphi_\alpha(x) \nabla \varphi_\beta(x) \, dx \\ &= \int_{\widehat T} \lambda(\Phi(\hat x)) \big[ F^{-t} \nabla \widehat \varphi_\alpha (\hat x) \big] \big[ F^{-t} \nabla \widehat \varphi_\beta (\hat x) \big] \, |J| \, d \hat x \end{align*}\]

We define the \(d \times N_T\) matrix-function \(B\) with columns

\[ B_{:\alpha}(\hat x) = F^{-t} \nabla \widehat \varphi_\alpha (\hat x), \qquad 1 \leq \alpha \leq N_T \]

With numerical integration we obtain the element matrix

\[ A_T = \sum_{k=1}^m \omega_k \lambda(\Phi(\hat x_k)) | J(\hat x_k) | \, B^t(\hat x) B(\hat x) \]