Mesh (triangulation) ${\mathcal T}$ consisting of
Idea: the material consists of small triangles (segments, tetrahedra) that act like basic building blocks of a material
Functions that are piecewise linear (more precisely, piecewise affine) with respect to $\mathcal T$:
$$
u_h(x) = \sum_{i=1}^N c_i v_i(x)
$$
Such space of function is denoted as ${\mathcal P}^1(\mathcal T)$
Often we'll use $$ \mathcal P^1_0(\mathcal T) := \{u\in\mathcal P^1_0(\mathcal T) : u|_{\Omega} = 0 \} $$
Discretize the equations...
First we need the variational (weak) form
Start with the equation (now I like to have minus in front of $\Delta$): $$ \begin{align*} -\Delta u &= f \\ u|_{\Gamma} &= 0 \end{align*} $$
Multiply the equation by a function $v$ (called the test function) such that $v|_{\Gamma}=0$ and integrate: $$ -\int_{\Omega} (\Delta u) v =- \int_{\Gamma} \nabla u v \cdot n + \int_{\Omega} \nabla u \cdot \nabla v =\int_{\Omega} f v \qquad\forall v $$
because $v=0$ on the boundary, we have $$ \int_{\Omega} \nabla u \cdot \nabla v =\int_{\Omega} f v \qquad\forall v $$
Recall from Lecture 2 that the energy of interaction of a two-dimensional network of springs was something like (in 1D): $$ E^{\rm 1d}_h(u) = \sum_x k (u(x+h) - u(x))^2/2 - f(x) u(x), $$ where $f(x)$ is the external force acting on the mass $x$.
if we let $k=1/h$, scale $f$ accordingly, and let $h\to 0$ then we have that $$ E(u) = \int_\Omega |\nabla u|^2/2 - \int_\Omega f u $$ (this formula works in any dimension)
Now, we can ask a question: $$ \text{find $u$ such that $E(u)$ is minimal} $$
To find the equation for $u$, one can take a functional derivative (equivalent to the principle of virtual displacements of mechanics): $$ \langle \partial E(u), v \rangle = \lim_{\epsilon\to 0}\frac{E(u+\epsilon v)}{\epsilon} = \int_\Omega \nabla u \cdot \nabla v - \int_\Omega f v = 0, \qquad\forall v $$ so we derived the variational form without the strong form
In general, the variational form is obtained directly from the energy minimum principle.
(BC = boundary condition)
Suppose we have a problem $$ \begin{align*} -\Delta u &= f\qquad\text{on $\Omega$} \\ u &= 0\qquad\text{on $\Gamma_1$} \\ u_n &= 0\qquad\text{on $\Gamma_2$} \end{align*} $$ where $u_n$ is the normal derivative and $\Gamma_1 \cup\Gamma_2 = \partial\Omega$.
The trick is then to introduce the space $X = \{\text{function }u : u|_{\Gamma_1}=0 \}$. Then we have, for $v\in X$ $$ -\int_{\Omega} \Delta u v = - \int_{\Gamma_1} u_n v - \int_{\Gamma_2} u_n v + \int_{\Omega} \nabla u\cdot \nabla v $$
Magically, the first integral =0 because $v=0$ on $\Gamma_1$ and the second integral =0 because $u_n=0$ on $\Gamma_2$.
Same variational formulation, only spaces are different
Denote $A(u,v) = \int_{\Omega} \nabla u\cdot\nabla v$ and $F(v) := \int_\Omega f v$.
Discrete equations: $$ A(u_h, v_h) = F(v_h) \qquad \forall v_h\in \mathcal P^1_0(\mathcal T) $$ (That's it!)
But to implement, we do $$ A(u_h, v_\ell) = F(v_\ell), \qquad \ell=1,...,N-1 $$ (i.e., it is enough to test only with basis functions)
Then we substitute $u_h = \sum_{k=1}^{N-1} c_k v_k$: $$ \sum_{k=1}^{N-1} c_k A(v_k, v_\ell) = F(v_\ell) \qquad \ell=1,\ldots,N-1, $$
Thus, $A(v_k, v_\ell)$ are the elements of the matrix (often called the stiffness matrix), and $F(v_\ell)$ are the components of the right-hand side
It remains to solve the linear system...
The mesh consists of triangles (tetrahedra, quadrilaterals, ...)
The mesh is defined by, typically, 3 arrays:
nodes:
$x_1$, $y_1$
$x_2$, $y_2$
...
$x_N$, $y_N$
(i.e., the first node has coords $(x_1, y_1)$, etc.)
triangles:
$n_1$, $n_2$, $n_3$
$m_1$, $m_2$, $m_3$
...
(i.e., the first triangle has nodes $n_1$, $n_2$, $n_3$ as vertices, etc.)
faces:
$n_1$, $n_2$
...
(i.e., the first face on the boundary has nodes $n_1$, $n_2$, etc. )
Basis functions look like this:
Classical way:
Notice that $\nabla v_k$ are piecewise constant: $$ A_{k,\ell} = \int_{\Omega} \nabla v_k \nabla v_\ell = \sum_{T\in\mathcal T} |T| (\nabla v_k|_T) (\nabla v_\ell|_T), $$ where $|T|$ is the volume of $T$
The algorithm would the look something like this:
for k = 1 to n
for l = 1 to n
for all T
if(k and l are nodes of T)
// otherwise the integral is zero
A(k,l) += |T|
*(\nabla v_k|_T)
*(\nabla v_l|_T)The loops in the algorithm are often reversed:
for all T
for k, vertices of T
for l, vertices of T
A(k,l) += <<as before>>Advantage: can now loop only over 3 vertices of each triangle, not waisting time
We have $$ 2 |T| = \det \left( \begin{array}{cc} x _2 - x _1 & y _2 - y _1 \\ x _3 - x _1 & y _3 - y _1 \\ \end{array} \right) $$ Let $T$'s nodes being $n_1$, $n_2$, $n_3$. Denote $\eta_i := v_{n_i}$. Then it can be shown $$ \nabla \eta _j = \frac{1}{2|T|} \left( \begin{array}{cc} y _{j+1} - y _{j+2} \\ x _{j+2} - x _{j+1} \\ \end{array} \right) $$ Here we mean $(x_4,y_4) = (x_1,y_1)$, $(x_5,y_5) = (x_2,y_2)$
One can do calculations and derive that the matrix $M_{j,k} = |T| (\nabla \eta_j)\cdot(\nabla \eta_k)$ can be evaluated as $$ M = \frac{|T|}{2} G G^T \qquad\text{where}\qquad G = \left( \begin{array}{ccc} 1 & 1 & 1 \\ x _1 & x _2 & x _3 \\ y _1 & y _2 & y _3 \end{array} \right)^{-1} \left( \begin{array}{ccc} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{array} \right) $$
the pseudo-code would then look like this:
for all T
calculate M
for k=1..3
for l=1..3
A(triangles(k),triangles(l)) += M(k,l)More details in the Remarks around 50 lines of Matlab: short finite element implementation
for all T
calculate f(xS, yS)
for k=1..3
f(triangles(k)) += 1/3 * area(T) * f(xS, yS)A = zero matrix for free nodes
for k = free nodes
for l = free nodes
for all T
<<SAME>>This fills only the needed rows & columns of the matrix
The loops in the algorithm are often reversed:
A = zero matrix for all nodes
for all T
<<SAME>>
remove rows and columns from A corr. to non-free nodesThe free nodes are those that are on $\Gamma_1$ in $$ \begin{align*} -\Delta u &= 0\qquad\text{on $\Omega$} \\ u &= 0\qquad\text{on $\Gamma_1$} \\ u_n &= 0\qquad\text{on $\Gamma_2$} \end{align*} $$
Alternatively, one can replace the corresponding row and column by $$ \begin{pmatrix} \cdot & 0 & \cdot & \cdot \\ 0 & 1 & 0 & 0\\ \cdot & 0 & \cdot & \cdot \\ \cdot & 0 & \cdot & \cdot \\ \end{pmatrix} $$
Exercise: think why
If $u$ is the exact solution and $u_h$ is the approximate solution then one can prove that $$ \|\nabla u_h - \nabla u\|_{L^2} \leq h \|\nabla^2 u\|_{L^2}, $$ where $h$ is the maximal size of the triangle and $\nabla^2 u$ is the Hessian matrix.
Comes from a truly beautiful argument:
In [6]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/alex.css", "r").read()
return HTML(styles)
css_styling()
Out[6]: