2.5 Stiffness Matrix and Right Hand Side

When a regular triangulation $\mathcal{T}$ has been generated for the domain $\Omega$, the space $S$ of the numerical solution $U$ has to be defined. A common choice of basis functions for the spline spaces $S$ and $S_D$ are ``hat functions'' (Fig. 2.2), which are defined for every node $(x_j,y_j,z_j)$ of the finite element mesh as

$$\eta_j(x_k,y_k,z_k)=\delta_{jk} \quad (j,k=1, \ldots, N) \quad.$$ (2.19)

We are using isoparametric elements, i.e. we use the same polynomials (linear basis functions) for the approximation of the geometry and the solution. If we define

$$S_D:=\mathrm{span}\{\eta_j\vert(x_j,y_j,z_j) \not\in \Gamma_D\}$$ (2.20)

then

$$U_D \in S_D$$

and we can calculate the stiffness matrix $A_{jk}$ (Eq. (2.15)) and the right hand side $b_j$ (Eq. (2.16)) as a sum over all elements $T$ and surface triangles $E$ on $\Gamma_N$

$$A_{jk}=\sum_{T \in \mathcal{T}} \int_T \nabla \eta_j \cdot \nabla \eta_k d{v}$$ (2.21)

and

$$b_j= \sum_{T \in \mathcal{T}} \int_T f \cdot \eta_j d{v}... ... \int_\Omega \nabla \eta_j \cdot \nabla \eta_k d{v} \quad.$$ (2.22)

Figure 2.2: Hat function (linear basis function) for a triangulation in 2D.

It is most convenient to calculate the stiffness matrix on an element by element basis (local or element matrices) and finally assemble the contributions from the local matrices to the global stiffness matrix.

If we assume that the four vertices of a tetrahedral element $T$ are given by $(x_j,y_j,z_j)$ with $j=1, \ldots, 4$, then the volume $\vert T\vert$ of the element is given by

$$\vert T\vert= \frac{1}{6} \mathrm{det} \left( \begin{ar... ...j+2} \\ 1 & x_{j+3} & y_{j+3} & z_{j+3} \end{array} \right)$$ (2.23)

where the local numbering $j=1, \ldots, 4$ is chosen in such a way, that the right hand side of Eq. (2.23) is positive.

The corresponding basis functions are given by

$$\eta_j(x_k,y_k,z_k)=\delta_{jk} \quad, \quad j,k=1, \ldots, 4 \quad.$$

Thus, $\eta_j$ can also be written as

$$\eta_j(x,y,z)= \mathrm{det} \left( \begin{array}{cccc} ... ...j+2} \\ 1 & x_{j+3} & y_{j+3} & z_{j+3} \end{array} \right)$$ (2.24)

and

$${ \nabla \eta_j(x,y,z)= \frac{1}{6\vert T\vert} \cdot }$$

$$\left( \begin{array}{c} y_{j+2} z_{j+1}-y_{j+3} z_{j+1}-y_{j+1} z... ...j+1} y_{j+2}+x_{j+3} y_{j+2}+x_{j+1} y_{j+3}-x_{j+2} y_{j+3} \end{array}\right)$$

where all indices are understood modulo 4.

As a result we can easily calculate the stiffness matrix entries

$$A_{jk}=\sum_{T \in \mathcal{T}}\int_T \nabla \eta_j (\nabla \eta_k)^T d{v} \quad.$$ (2.25)

For the right hand side of Eq. (2.16) we need to evaluate $\int_T f \cdot \eta_j d{v}$. If we use the value of $f$ in the center of gravity $(x_S,y_S,z_S)$ of $T$, we can make the approximation

$$\int_T f \cdot \eta_j d{v} \approx \frac{\vert T\vert}{4} f(x_S,y_S,z_S) \quad.$$ (2.26)

The second term of the right hand side (Eq. (2.16)) can be evaluated in a similar way. However, Neumann boundary conditions will not be required in the following.

Finally, the Dirichlet boundary conditions have to be incorporated. One straight forward and easy to implement method is to replace all rows of the stiffness matrix $A_{jk}$, which correspond to Dirichlet boundary nodes, with zero and a single one in the main diagonal. On the right hand side, the entries of the Dirichlet nodes are replaced with their boundary values.