Implementation of Finite Element Methods

Per-Olof Persson
persson@berkeley.edu

Math 228B Numerical Solutions of Differential Equations

The Poisson Problem in 2-D

  • Consider the problem \[\begin{aligned} -\nabla^2 u &= f \text{ in }\Omega \\ n\cdot\nabla u &= g \text{ on } \Gamma \end{aligned}\] for a domain \(\Omega\) with boundary \(\Gamma\)

  • Seek solution \(\hat{u}\in\hat{X}\), multiply by a test function \(v\in\hat{X}\), and integrate: \[\begin{aligned} \int_\Omega -\nabla^2 \hat{u} v\,d\Omega = \int_\Omega fv\,d\Omega \end{aligned}\]

  • Apply the divergence theorem and use the Neumann condition, to get the Galerkin form \[\begin{aligned} \int_\Omega \nabla \hat{u} \nabla v\,d\Omega = \int_\Omega fv\,d\Omega + \oint_\Gamma g v\,ds \end{aligned}\]

Finite Element Formulation

  • Expand in basis \(\hat{u}=\sum_i \hat{u}_i \phi_i(x)\), insert into the Galerkin form, and set \(v = \phi_i\), \(i=1,\ldots,n\): \[\begin{aligned} \int_\Omega \left[\sum_{j=1}^n \hat{u}_j\nabla\phi_j\right]\cdot\nabla \phi_i\,d\Omega = \int_\Omega f\phi_i\,d\Omega + \oint_\Gamma g\phi_i\,ds \end{aligned}\] Switch order of integration and summation to get the finite element formulation: \[\begin{aligned} \sum_{j=1}^n A_{ij} \hat{u}_j = b_i,\quad\text{or}\quad A\hat{u} = b \end{aligned}\] where \[\begin{aligned} A_{ij} = \int_\Omega \nabla\phi_i\cdot\nabla\phi_j\,d\Omega, \quad b_i = \int_\Omega f \phi_i + \oint_\Gamma g\phi_i\,ds \end{aligned}\]

Discretization

  • Find a tringulation of the domain \(\Omega\) into triangular elements \(T^k\), \(k=1,\ldots,K\) and nodes \(x_i\), \(i=1,\ldots,n\)

  • Consider the space \(\hat{X}\) of continuous functions that are linear within each element

  • Use a nodal basis \(\hat{X}=\mathrm{span}\{\phi_1,\ldots,\phi_n\}\) defined by \[\begin{aligned} \phi_i\in \hat{X},\quad \phi_i(x_j)=\delta_{ij},\quad 1\le i,j\le n \end{aligned}\]

  • A function \(v\in\hat{X}\) can then be written \[\begin{aligned} v=\sum_{i=1}^n v_i\phi_i(x) \end{aligned}\] with the nodal interpretation \[\begin{aligned} v(x_j)=\sum_{i=1}^n v_i\phi_i(x) = \sum_{i=1}^n v_i\delta_{ij} = v_j \end{aligned}\]

image

Local Basis Functions

  • Consider a triangular element \(T^k\) with local nodes \(x_1^k, x_2^k, x_3^k\)

  • The local basis functions \(\mathcal{H}_1^k,\mathcal{H}_2^k,\mathcal{H}_3^k\) are linear functions: \[\begin{aligned} \mathcal{H}_\alpha^k = c_\alpha^k + c_{x,\alpha}^k x + c_{y,\alpha}^k y,\quad \alpha=1,2,3 \end{aligned}\] with the property that \(\mathcal{H}_\alpha^k(x_\beta) = \delta_{\alpha\beta}\), \(\beta=1,2,3\)

  • This leads to linear systems of equations for the coefficients: \[\begin{aligned} \begin{pmatrix} 1 & x_1^k & y_1^k \\ 1 & x_2^k & y_2^k \\ 1 & x_3^k & y_3^k \end{pmatrix} \begin{pmatrix} c_\alpha^k \\ c_{x,\alpha}^k \\ x_{y,\alpha}^k \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \text{ or } \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \text{ or } \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \end{aligned}\] or \(C = V^{-1}\) with coefficient matrix \(C\) and Vandermonde matrix \(V\)

Elementary Matrices and Loads

  • The elementary matrix for an element \(T^k\) becomes \[\begin{aligned} A^k_{\alpha\beta} &= \int_{T^k} \frac{\partial \mathcal{H}_\alpha^k}{\partial x} \frac{\partial \mathcal{H}_\beta^k}{\partial x} + \frac{\partial \mathcal{H}_\alpha^k}{\partial y} \frac{\partial \mathcal{H}_\beta^k}{\partial y}\,d\Omega \\ &= \mathrm{Area}^k (c_{x,\alpha}^k c_{x,\beta}^k + c_{y,\alpha}^k c_{y,\beta}^k),\quad \alpha,\beta=1,2,3 \end{aligned}\]

  • The elementary load becomes \[\begin{aligned} b^k_\alpha &= \int_{T^k} f\,\mathcal{H}_\alpha^k\,d\Omega \\ &= \text{(if $f$ constant)} \\ &= \frac{\mathrm{Area}^k}{3}f,\quad \alpha=1,2,3 \end{aligned}\]

image

Assembly, The Stamping Method

  • Assume a local-to-global mapping \(t(k,\alpha)\), giving the global node number for local node number \(\alpha\) in element \(k\)

  • The global linear system is then obtained from the elementary matrices and loads by the stamping method:

    \(A=0\), \(b=0\)
    for \(k=1,\ldots,K\)
    \(A(t(k,:),t(k,:)) = A(t(k,:),t(k,:)) + A^k\)
    \(b(t(k,:)) = b(t(k,:)) + b^k\)

image

Dirichlet Conditions

  • Suppose Dirichlet conditions \(u=u_D\) are imposed on part of the boundary \(\Gamma_D\)

  • Enforce \(\hat{u}_i = u_D\) for all nodes \(i\) on \(\Gamma_D\) directly in the linear system of equations: \[\begin{aligned} & \qquad\qquad\qquad\,i \\ \begin{array}{c} \\ \\ i \\ \\ \\ \end{array} &\begin{pmatrix} \\ \\ 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\ \\ \\ \end{pmatrix} \hat{u} = \begin{pmatrix} \\ \\ u_D \\ \\ \\ \end{pmatrix} \end{aligned}\]