Poisson equation by the finite element method

(If the reader is interested in the details, she/he is advised to read an excellent and cheap book Numerical Solution of Partial Differential Equations by the Finite Element Method by Claes Johnson as a reference.)

Let us consider the boundary value problem: Find $u \in H^1_0(\Omega)$ such that
$$
\begin{aligned}
-\Delta u &= F,\quad\text{in $\Omega$},\\
u &= 0,\quad\text{on $\partial \Omega$},\end{aligned}$$ where $\Omega = (0,1)^2$. Analytical solution exists for some $F$ but if $\Omega$ is arbitrary, the typical approach is to use the finite element method: The problem is casted into the weak formulation and then solved by choosing an appropriate discrete subspace $V_h \subset H^1_0(\Omega)$ and by testing with all the basis functions of $V_h$ to get a matrix equation. The weak formulation reads: Find $u \in V:=H^1_0(\Omega)$ such that
$$(\nabla u, \nabla v) = (F,v),$$for every $v \in V$ where we use the notation$$(a,b)_\Omega := \int_\Omega a\cdot b \, \mathrm{d}x = (a,b).$$After we have defined a basis $\{\phi_j\}, j=1,\dots,N$ for the discrete space $V_h \subset V$ (of course $V_h = \text{span}(\phi_1, \dots \phi_N)$) we get the discrete weak formulation: Find $u_h \in V_h$ such that
\begin{equation}
\label{eq:discreteweak}
(\nabla u_h, \nabla v_h) = (F,v_h),
\end{equation} for every $v_h \in V_h$. Céa's lemma guarantees that the solution $u_h$ is the best in the space $V_h$ as long as $V_h \subset V$. Furthermore, if we define the basis $\phi_j, j=1,\dots,N$ by using a triangulation $\mathcal{T}_h = \{K\}$ ($h$ refers to the maximum circumcircle diameter of $K$'s) of the domain $\Omega = \cup_{K \in \mathcal{T}_h} \overline{K}$ such that
$$\phi_j(N_i) = \delta_{ij}$$ where $\delta_{ij}$ is the Kronecker delta symbol and such that $\phi_j \in \mathcal{P}_1(K)$ for every $j$ then by using something known as the Bramble-Hilbert lemma and some scaling argument we can prove that
$$\|u-u_h\|_{H^1(\Omega)} \leq C h |u|_{H^2(\Omega)},$$ if the domain $\Omega$ contains no "evil singularities" and the load $F$ is smooth enough. Such a priori estimate gives some sort of guarantee that the numerical solution $u_h$ will eventually start to be twice as close to the real solution if we refine the mesh such that $h$ gets halved. This is something we next verify numerically by building our own solver.

By substituting $u_h = \sum_{j=1}^N a_j \phi_j$ to the equation \eqref{eq:discreteweak} and by testing with just the basis functions we get
$$ \sum_{j=1}^N a_j (\nabla \phi_j, \nabla \phi_i) = (F,\phi_i),\quad i=1,\dots,N$$which can be written by using matrix notation as
$$ Sx = f,$$where $S_{ij} = (\nabla \phi_i, \nabla \phi_j)$, $f_i = (F,\phi_i)$ and $x_i = a_i$. Calculation of the entries $S_{ij}$ is rather cumbersome in the global coordinate frame and therefore we will do the numerical integration on a reference triangle $\hat{K}$ with vertices $(0,0), (1,0), (0,1)$.

Exercise: Suppose we have a global triangle $K$ with vertices $n_1,n_2$ and $n_3$. Then the affine mapping $\xi: \hat{K} \rightarrow K$ can be written explicitly as
\begin{equation}
\label{eq:affmap}
\xi(\hat{x}) = A\hat{x}+b,
\end{equation} where
$$A = \begin{bmatrix}n_2 -n_1 &n_3-n_1\end{bmatrix},\quad b=n_1,$$ and $n_j$ for all $j$ are considered as column vectors. $\blacksquare$

To make the calculation more efficient, we note that in the summation
\begin{equation}
\label{eq:summand}
S_{ij} =  (\nabla \phi_i, \nabla \phi_j)_\Omega = \sum_{K \in \mathcal{T}_h} (\nabla \phi_i, \nabla \phi_j)_K,
\end{equation} most of the summands integrate to zero. Thus, we will end up with a more efficient procedure if we're not even trying to evaluate such terms in the first place. It is quite straightforward to see that if we just calculate the matrices
\begin{equation} \label{eq:localstiffness} \hat{S}_K = \begin{bmatrix}
(\nabla \phi_{K_1},\nabla \phi_{K_1})_K & (\nabla \phi_{K_1},\nabla \phi_{K_2})_K & (\nabla \phi_{K_1},\nabla \phi_{K_3})_K \\
(\nabla \phi_{K_2},\nabla \phi_{K_1})_K & (\nabla \phi_{K_2},\nabla \phi_{K_2})_K & (\nabla \phi_{K_2},\nabla \phi_{K_3})_K \\
(\nabla \phi_{K_3},\nabla \phi_{K_1})_K & (\nabla \phi_{K_3},\nabla \phi_{K_2})_K & (\nabla \phi_{K_3},\nabla \phi_{K_3})_K
\end{bmatrix},\end{equation} for all $K\in \mathcal{T}_h$, we end up having all the terms required to perform the summation \eqref{eq:summand}. ($K_j$ is the $j$'th node of triangle $K$) However, we require some sort of knowledge how these so-called "local stiffness matrices" should be summed into the "global stiffness matrix" $S$. For this procedure (often referred to as "assembly of the global stiffness matrix") we define a operator $\mathcal{A}_K: \mathbb{R}^{3\times 3} \rightarrow \mathbb{R}^{N \times N}$ such that
$$(\mathcal{A}_K(X))_{ij} = \begin{cases}
X_{kl}, &\text{if $K_k =i$ and $K_l = j$ for $k,l = 1,2,3$}\\
0,&\text{otherwise.}
\end{cases}$$By using this "assembly operator" we may write
$$S = \sum_{K \in \mathcal{T}_h} \mathcal{A}_K(\hat{S}_K).$$

Let us start coding. Today's language of choice is Mathematica. Unfortunately, the author is not a very experienced Mathematica hacker and therefore the quality of code should be questioned. Any improvements I hear gladly in the commenting section.

We store our mesh node points in a N-by-2 matrix where the first column contains the x-coordinates and the second column the y-coordinates. Thus, let us generate nodes uniformly on the unit square.

mesh[p] =
  Flatten[Outer[{#1,#2}&, Range[0,1,0.1], Range[0,1,0.1]], 1];

Next we need to represent our triangles somehow. For that we use a N-by-3 matrix where each row $j$ contains indices of the $j$'th triangle corresponding to the row indices of the node point matrix. I could have chosen to use the transpose of such matrix instead but I think this is more convenient in the Mathematica-context.

mesh[t] = DelaunayFast[mesh[p]];

The function DelaunayFast generates a Delaunay triangulation of the points in mesh[p]. The function uses the built-in function ListDensityPlot and the implementation I found from StackOverflow reads

Needs["ComputationalGeometry`"];
DelaunayFast[points_] := 
 Module[{pl},
  pl = ListDensityPlot[ArrayPad[points,{{0,0},{0,1}}]];
  Cases[pl, Polygon[a_] :> a, Infinity][[1]]]

In order to visualize the mesh I made an implementation of the function TriMesh which is more or less familiar from the MATLAB environment.

LeaveIfContains[x_, y_] := Select[x, MemberQ[#, y] &];
Tri2Graph[t_] := {#,Union[Flatten[LeaveIfContains[t, #]]]}& /@
   Union[Flatten[t]];
TriMesh[mesh_] := PlanarGraphPlot[mesh[p], Tri2Graph[mesh[t]]];

Now by running

TriMesh[mesh]

we get a picture like

A triangular mesh of the unit square.

For convenience (and to fight Mathematica's annoying [[]]-syntax), we define the following shorthands

mesh[p, k_] := mesh[p][[k]]
mesh[t, k_, l_] := mesh[t][[k, l]]

Then the code to generate all the affine mappings \eqref{eq:affmap} is

A = Table[{mesh[p,mesh[t,k,2]] - mesh[p,mesh[t,k,1]], 
     mesh[p,mesh[t,k,3]] - mesh[p,mesh[t,k,1]]}//Transpose,
     {k, 1, Length[mesh[t]]}];

It is generally a good idea to try to vectorize any interpreted code for performance. Even though this is the first time I'm implementing a solver for Mathematica, I decide to follow the same route I'd use in MATLAB or Numpy: I will write a function which evaluates the elements (i,j) of the local stiffness matrices for all elements at once. To actually evaluate the elements of $\hat{S}_K$ we will transform the integration into the local reference element and use a numerical quadrature. On the reference element the linear basis functions corresponding to the nodes $(0,0), (1,0)$ and $(0,1)$ are $$\begin{aligned} \hat{\phi}_1 &= 1-\hat{x}-\hat{y},\\ \hat{\phi}_2 &= \hat{x}, \\ \hat{\phi}_3 &= \hat{y}. \end{aligned}$$ Exercise: Show that the transformation formula for the gradients is $$\nabla \phi_{K_i} = A_K^{-T} \hat{\nabla} \hat{\phi}_i,$$where $\hat{\nabla}$ denotes the gradient in the local coordinates and $i=1,2,3$. $\blacksquare$

Because the basis functions are linear and defined on a triangle, we may just use a very simple numerical quadrature (which is exact for linear $f$): $$\int_{\hat{K}} f(\hat{x},\hat{y})\,\mathrm{d}x\,\mathrm{d}y \approx \frac 12 f(\frac 13, \frac 13).$$ Thus, the element of the local stiffness matrix is $$(\hat{S}_K)_{ij} = (A_K^{-T} \hat{\nabla} \hat{\phi}_i,A_K^{-T} \hat{\nabla} \hat{\phi}_j)_{\hat{K}} |\det A| = \frac 12\langle A_K^{-T} \hat{\nabla} \hat{\phi}_i,A_K^{-T} \hat{\nabla} \hat{\phi}_j\rangle |\det A|,$$ where $\langle.\rangle$ denotes the standard vector dot product. Now we just may define our local basis as

Phi = {1 - #1 - #2 &, #1 &, #2 &};
DPhi = {{-1., -1.}, {1., 0.}, {0., 1.}};

and then the code to build the elements $(i,j)$ of the local stiffness matrices (for all elements at once) reads

Shat[i_, j_] := 
  0.5 Dot[Inverse[Transpose[#]].DPhi[[i]], 
      Inverse[Transpose[#]].DPhi[[j]]]*Abs[Det[#]] & /@ A;

Similarly for the local load vector, i.e. $$(\hat{f}_K)_i = \frac 12 F(A_K \frac 13 + b_K) \hat{\phi}_i(\frac 13, \frac 13) | \det A_k |,$$ we have

F[x_,y_] := 1.0;
b = Table[mesh[p,mesh[t,k,1]], {k, 1, Length[mesh[t]]}];
fhat[i_] := Phi[[i]][1/3, 1/3]*
    Thread[X[A,b]]/.{X->(0.5 F@@(#1.{1/3,1/3}+#2)*Abs[Det[#1]]&)};

(Sorry for the confusing substitution including X. For some reason Thread didn't work as I expected based on the manual, so I had to come up with a work-around.)

After these function definitions, the global stiffness matrix and load vector assembly is as simple as this:

S = Sum[SparseArray[
    Table[{mesh[t,k,i],mesh[t,k,j]}, {k,1,Length[mesh[t]]}] -> 
    Shat[i,j],{Length[mesh[p]],Length[mesh[p]]}],{i,1,3},{j,1,3}];
f = Sum[SparseArray[
    Table[{mesh[t,k,i],1},{k,1,Length[mesh[t]]}] -> 
    fhat[i],{Length[mesh[p]],1}],{i,1,3}];

Note! This of course assumes that we have set

SetSystemOptions["SparseArrayOptions" ->
    "TreatRepeatedEntries"->Total];

since we get double index pair entries into the SparseArray function call.

We're left with finding the Dirichlet nodes (the boundary nodes) and then solving the system in the inner nodes. To find the boundary nodes I wrote the following short query.

ID = 
  Union[Flatten[Position[mesh[p], #]]] & /@ {{0.,_}, {1.,_},
   {_,0.}, {_, 1.}} // Flatten // Union;

The inner nodes are just a complement of this set:

II = 
  Complement[Range[1, Length[mesh[p]]], ID];

Now we can solve the system by

u = ConstantArray[0., Length[mesh[p]]];
u[[II]] = LinearSolve[S[[II, II]],f[[II]] // Flatten;

and we can visualize the solution by using an implementation of MATLAB's TriSurf I wrote for the occasion:

DrawTriangle[st_, p_, u_] := 
  Graphics3D[Polygon[Append[p[[#]], u[[#]]] & /@ st]];
TriSurf[mesh_, u_] := 
  Show[DrawTriangle[#, mesh[p], u] & /@ mesh[t],
   PlotRange -> All, Axes -> True, BoxRatios -> {1,1,.5}];

By running

TriSurf[mesh,u]

we get the beautiful picture

The solution $u_h$ of the Poisson equation with unit load.

To be continued ...