It is very common to use Lagrange polynomials for basis functions
PDEs are normally not solved using tensor products in multiple dimensions
Solving PDEs with FEM
The Galerkin method is as always: Find \(u_N \in V_N\) such that \[
(\mathcal{R}_N, v) = 0, \quad \forall \, v \in V_N,
\] where \(\mathcal{R}_N = \mathcal{L}(u_N)-f\) is a residual and \(\mathcal{L}(u)\) is some mathematical operator.
Choose \(V_N\) and everything follows.
Insert for \(u_N = \sum_{j=0}^N \hat{u}_j \psi_j\) and \(v=\psi_i\) to obtain the linear algebra problem.
All the complexity lies in the implementation, using elementwise integration and mapping to a reference domain.
It is very common (and some times necessary) to use integration by parts or Green’s first identity in order to manipulate the variational equation and to incorporate natural Neumann boundary conditions.
Boundary conditions are not incorporated into the basis functions.
We will only use (local) Lagrange polynomials as basis functions, but there are many other options.
Example - Poisson’s equation with inhomogeneous Dirichlet boundary conditions
\[
u'' = f, \quad x \in (0, L), \, u(0) = a, u(L) = b,
\]
First choose \(V_N\). We can here use piecewise linear polynomials or even higher order
\[
V_N = \text{span}\{\psi_j(x)\}_{j=0}^N
\]
Remember
\[
\psi_j(x_i) = \delta_{ij}
\]
so each basis function is one in exactly one node and zero in all the others. Quadratic:
Use integration by parts
The variational Galerkin form is now to find \(u_N \in V_N\) such that
Since \(\psi_i(0)=1\) for \(i=0\) and \(\psi_i(L)=1\) for \(i=N\) and zero for all \(i=1, 2, \ldots, N-1\) the last two terms are only relevant for \(i=0\) and \(i=N\).
Dirichlet boundary conditions
We have \(N+1\) equations for \(N+1\) unknowns \(\{\hat{u}_i\}_{i=0}^N\)\[
\sum_{j=0}^N(\psi'_j, \psi'_i) \hat{u}_j = -(f, \psi_i) + u'_N(L)\psi_i(L) - u'_N(0)\psi_i(0)
\]
The boundary term \(u'_N(L)\psi_i(L) - u'_N(0)\psi_i(0)\) can be used to specify Neumann boundary conditions. For Dirichlet we have
is the element stiffness matrix of shape \((N+1) \times (N+1)\). We also use the same mapping from local to global degrees of freedom, and define a dense element matrix as
for finite elements of order \(d\) and \(\mathcal{I}_d=\{0, 1, \ldots, d\}\). Remember that \(d=1\) for linear elements, and then higher order elements (\(d>1\)) simply use more nodes within each element.
Mapping to reference domain \(\tilde{\Omega} = [-1, 1]\)
The finite element stiffness matrix is implemented using a mapping to the reference domain \(X\in[-1, 1]\), and the reference basis functions
\[
\psi_{q(e, r)}(x) = \ell_r(X)
\]
where the element is \(\Omega^{(e)} = [x_L, x_R]\) and \(h(e)=x_R-x_L\), such that
where \(\ell^{(n)}_r = \frac{d^{n} \ell_r}{dX^n}\).
Note
This represents all the element matrices you will ever need, and it is easily implemented using a few lines of Sympy code.
Note
Higher order of the derivatives \((m, n)\) requires higher order elements. Remember, piecewise linear elements can only be differentiated once.
Sympy implementation of generic element matrix
Using uniform nodes
\[
X_j = -1 + \frac{2j}{d}, \quad j=0,1,\ldots, d
\]
Xj =lambda d: np.array([sp.Rational(2*i, d) for i in np.arange(d+1)])-1qe =lambda l, r, s, m, n: sp.integrate(l[r].diff(x, m)*l[s].diff(x, n), (x, -1, 1))def Qe(d=1, m=0, n=0): A = np.zeros((d+1, d+1), dtype=object) l = Lagrangebasis(Xj(d))for r inrange(d+1):for s inrange(d+1): A[r, s] = qe(l, r, s, m, n)return (h/2)**(1-m-n)*sp.Matrix(A)
Assume manufactured solution \(u(x)=J_0(x)\), where \(J_0(x)\) is the 0 order Bessel function of the first kind. Use Dirichlet boundary conditions, \(u(0)=J_0(0)\) and \(u(10)=J_0(10)\), and thus neglect \([u'_Nv]_{x=0}^{x=10}\).
N =20xj = np.linspace(0, 10, N+1)d =1S = assemble_q(xj, d=d, m=1, n=1) M = assemble_q(xj, d=d, m=0, n=0)alpha =10
The boundary term \(J'_0(10)\psi_i(10)\) is only nonzero for \(i=N\) and \(J'_0(0) \psi_i(0)\) is only nonzero for \(i=0\). This follows since \(\psi_i(x_j) = \delta_{ij}\) and \(x_0=0\) and \(x_N=10\).
Note
The term \((J'_0(10)\psi_i(10) - J'_0(0) \psi_i(0))\) goes into \(\boldsymbol{b}\), because it is a linear form containing the test function \(\psi_i\) but not the trial function \(\psi_j\).
Solve
The coefficient matrix \(-S+\alpha M\) is unmodified (no identing first/last rows) and the right hand side \(\boldsymbol{b}\) takes care of the Neumann boundary conditions
But naturally the method is harder to implement. We need to use position vectors in real and reference space
\[
\boldsymbol{x} = x \boldsymbol{i} + y \boldsymbol{j} + z \boldsymbol{k} \quad \text{and} \quad \boldsymbol{X} = X \boldsymbol{i} + Y \boldsymbol{j} + Z \boldsymbol{k}
\]
and for 2D the piecewise linear basis function now look like “tents”
Reference triangle and reference basis functions
The reference triangle has three nodes \(\boldsymbol{X}_0= (0, 0), \boldsymbol{X}_1= (1, 0)\) and \(\boldsymbol{X}_2 = (0, 1)\), and the three nonzero, linear reference basis functions on this element are
In order to use piecewise linear (P1) elements for Poisson’s equation we need to use the “multidimensional integration by parts”
\[
\int_{\Omega} \nabla^2 u \,v \,d\Omega = -\int_{\Omega} \nabla u \cdot \nabla v \, d \Omega + \oint_{\partial \Omega} \nabla u \cdot \boldsymbol{n} \, v \, d \sigma,
\]
where \(\partial \Omega\) is the enclosing boundary of the domain \(\Omega\), \(\boldsymbol{n}\) is an outward pointing unit normal vector and \(d\sigma\) is a line element in 2D and a surface element in 3D. With \(L^2\) inner product notation:
\[
(\nabla^2 u, v)_{L^2(\Omega)} = - (\nabla u, \nabla v)_{L^2(\Omega)} + (\nabla u \cdot \boldsymbol{n}, v)_{L^2(\partial \Omega)}.
\]
Note
Piecewise linear basis functions have only one nonzero derivative in each direction. So for P1 elements \((\nabla^2u, v)=0\), but \((\nabla u, \nabla v)\) may be nonzero.
We here use the Jacobian \(J\) of the transformation, defined such that \[ \small
(J)_{ij} = \frac{\partial x_j}{\partial X_i}, \quad (J^{-1})_{ij} = \frac{\partial X_j}{\partial x_i}
\]
For each element collect the true coordinates of the nodes \(\boldsymbol{x}_{q(e)}\)
Use a precomputed matrix \(\nabla_X \tilde{\psi}\) and compute the Jacobian \(J\)
Compute the inverse of the Jacobian \(J^{-1}\)
Assemble element matrix, for example stiffness \[
\int_{\Omega^{(e)}} \nabla \psi_{q(e, s)}(\boldsymbol{x}) \cdot \nabla \psi_{q(e, r)}(\boldsymbol{x}) \, d\Omega = \int_{\tilde{\Omega}} \left(J^{-1} \cdot \nabla_X \tilde{\psi}_s(\boldsymbol{X})\right) \cdot \left(J^{-1} \cdot \nabla_X \tilde{\psi}_r(\boldsymbol{X})\right) \, \text{det}J\, d\tilde{\Omega}
\] For piecewise linear elements the integrand is a constant and you only need to compute the area of each cell. For higher order elements the integral requires some work!
Add element matrix to global matrix (finite element assembly). For all \(e\) and \(r\) and \(s\) do \[
a_{q(e, r), q(e, s)} \mathrel += \tilde{a}_{rs}
\]
A better implementation would use vectorization for the inner \(r, s\) loops.
Summary
The finite element method is easily adapted to complex meshes consisting of triangles in 2D or tetrahedra in 3D.
The implementation is mapping the true cells to reference triangles or tetrahedra. The mapping involves a Jacobian that is easy to compute from the linear Lagrange basis functions and cell nodes.
Higher order elements use internal points on each element.
Green’s first identity, which is similar to an integration by substitution in multiple dimensions, is used to transform weak variational forms.
Neumann boundary conditions are implemented directly in the weak form.
Dirichlet boundary conditions are implemented as always by identing rows in the coefficient matrix belonging to a Dirichlet boundary. The procedure eliminates boundary terms from the weak form.
There are several good open source finite element solvers out there! (e.g., FEniCS, Elmer, FreeFEM)
