HEqSolver

HEqSolver is a Finite Element Solver written in R for solving the diffusion equation. HEqSolver employes piecewise linear basis functions on uniform triangulation. It solves the diffusion equation on a two dimensional rectangular domain with homogeneous boundary conditions of Dirichlet type. We illustrate HEqSolver through an example of the diffusion equation posed on a rectangular domain in which we present both the relevant theory of the finite element method a long with the time-stepping scheme. All this is illustrated through the presentation of the mathematical theory of the finite element method on an example, for which the solution is obtained by HEqSolver.

Finite element method overview

Consider $\Omega=(0\,\,\,1)\times(0\,\,\,1)\subset\mathbb{R}^2$ , a rectangular open bounded domain with $\partial\Omega\subset\mathbb{R}^2$ denoting its boundary. We define the initial conditions on each point in space $(x,y)\in\Omega$ at time $t=0$ and hence consider the diffusion equation in two space and one time dimension for $0<t\leq\,T_{max}$ , where $T_{max}$ denotes the final time. Let $u(x,y,t)$ satisfy the probelm $\frac{\partial\,u}{\partial\,t}=D\Delta\,u(x,y)$ for $\!(x,y)\in\Omega,\,\,\,t\in(0\,\,\,T_{max}\]$ and $u(x,y,t)=0$ for $\quad(x,y)\in\partial\Omega$ and $t>0”/>. Here, <img src=$ 0"/> denotes the diffusion coefficient. Initial conditions for the problem are prescribed such that $u(x,y,0)=u_0(x,y)$ , where $u_0(x,y)\in\mathbb{L}_2(\Omega)$ , which means that it is a square integrable function of $x$ and $y$ . The approximate solution of time-dependent problems through numerical methods can in general be achieved through two different types of approaches. The two approaches mainly depend on whether the scheme is formulated to discretise time first followed by the space discretisation or vice vera. We adopt the approach in which we apply Galerkin finite element method to implement the space discretisation first and then we proceed with time-stepping using the implicit Euler scheme. Consider the diffusion problem to be posed on a space-time cylinder $\Omega\times(0\,\,\,T_{max}\]$ , in the form of a sequence of Poisson’s equations, that are all well-posed on $\Omega$ , i.e. satisfying all the requirements stated by the Lax-Milgram Theorem at each time point in the interval $(0\,\,\,T_{max}\]$ .
This allows to introduce at each time point $t$ , the solution and test spaces of functions for the corresponding sequence of Poisson equations to be denoted by $H^1_{\Omega}=\{w\in\,H^1:\|w(x,y)\|^2<\infty,\,\|\nabla\,w(x,y)\|^2<\infty\}$ and $H^1_{\Omega_0}=\{w\in\,H^1:\,w(x,y)=0,\;\,x,y\in\partial\Omega\}$ respectively, where $H^1_{\Omega}$ and $H^1_{\Omega_0}$ are a special class of Sobelov spaces. Assuming that the diffusion problem is well-posed on $\Omega$ for all $0<t\leq\,T$ </t\leq\,T”>, then we multiply both sides of the diffusion equation by a test function $w(x,y)\in\,H^1_{\Omega_0}$ and integrate by parts with application of Green’s Theorem. Using the prescribed boundary conditions lead to vanish the boundary terms in the application of integration by parts to the right hand-side of the equation, which results in the weak formulation of the diffusion problem and that is to find some $u\in\,H^1_{\Omega}$ at each time point $t$ such that $\int_{\Omega}\frac{\partial\,u}{\partial\,t}w\,dxdy=-D\int_{\Omega}\nabla\,u\,\cdot\nabla\,w\,dxdy$ for all $w\in\,H^1_{\Omega_0}$ . Given that the problem is well-posed in the continuous space of functions of which both $H^1_{\Omega}$ and $H^1_{\Omega_0}$ are subsets, we define the finite dimensional solution and test spaces namely $V^h\subset\,H^1_{\Omega}$ and $V_0^h\subset\,H^1_{\Omega_0}$ respectively for each time point $t$ . Let $\Omega^h$ be a triangulated approximation of $\Omega$ , containing a discrete integer number of nodes and $K$ triangles, then using the finite dimentional solution and test spaces, the finite element formulation of the problem is to find $u^h\in\,V^h$ such that $\int_{\Omega^h}\frac{\partial\,u^h}{\partial\,t}w^h\,dxdy=-D\int_{\Omega^h}\nabla\,u^h\cdot\nabla\,w^h\,dxdy$ for all $\,w^h\in\,V_0^h$ . We expand $u^h$ and $w^h$ in terms of its finite i.e. $N$ dimensional basis functions in the form $u^h=\sum_{i=1}^N\,U_i\phi_i$ and $w^h=\sum_{j=1}^N\,\phi_j$ . Substituting these in the finite element formulation we obtain a discrete set of continuous oridnary differential equations in the form $\frac{d}{dt}\Big(\int_{\Omega}\sum_{i=1}^NU_i(t)\phi_i\sum_{j=1}^N\phi_j\,dxdy\Big)=-D\int_{\Omega}\sum_{i=1}^NU_i(t)\nabla\,u\cdot\sum_{j=1}^N\nabla\,\phi_j\,dxdy$ . We use the commutative property of summation and integration to write the semi-discrete system of differential equations in matrix notation in the form $M\frac{dU}{dt}=-DSU(t)$ , where $S$ and $M$ respectively denote the stiffness and mass matrices. The vector $U$ contains the unknown approximate solution values at each node in the domain $\Omega^h$ at time $t$ . The entries of $S$ and $M$ are given by $\[S\]_{i,j}=\int_{\Omega^h}\,\nabla\phi_i\cdot\nabla\phi_j\,dxdy$ and $\[M\]_{i,j}=\int_{\Omega^h}\phi_i\phi_j\,dxdy$ respectively. The entries of $U$ are the finite element approximate solution values to be found from the system of ordinary differential equations.
In order to solve the semi-discrete problem we introduce a uniform mesh of $T$ discrete number of time points in the form $(0,\Delta\,t,2\Delta\,t,3\Delta\,t...,T\Delta\,t=T_{max})$ , which is equally spaced by time step-size $\Delta\,t$ . We use the Implicit Euler Method, to fully discretise the semi-discrete system of equations. Let $U^m$ denote the finite element approximate solution at the time-step $m\Delta\,t$ , then full discretisation of the problem is achieved by writing $M\frac{U^{m+1}-U^m}{\Delta\,t}=-DSU^{m+1}$ , which can be rearranged in the form of a system of $T$ linear equations for each time-step written as $U^{m+1}=\Big(S+D\Delta\,t\,M\Big)^{-1}MU^{m}$ .
Let $\Omega$ be discretised by $2N^2$ uniform triangles connected through $(N+1)^2$ nodes. Let $\mathcal{T}$ denote the uniform triangulation of $\Omega$ , consisting of all triangles. The algorithm is programmed to compute locally on each triangle $K\in\mathcal{T}$ , the entries of the mass and stiffness matrices namely $M$ and $S$ . The entries of the global mass and stiffness matrices are computed from contributions by the local $3\times3$ mass and stiffness matrices namely $M$ and $S$ on each triangle. The global mass and stiffness matrices are therefore, given by the formulea $M_{i,j}=\sum_{K\in\mathcal{T}}\int_{K}\phi_i\phi_jdxdy$ and $S_{i,j}=\sum_{K\in\mathcal{T}}\int_{K}\nabla\phi_i\cdot\nabla\phi_jdxdy$ respectively.

Code documentation

Recall the problem formulated in terms of $u(x,y,t)$ satisfying the diffusion equation $\frac{\partial\,u}{\partial\,t}=\Delta\,u(x,y)$ on $\Omega=(0\;1)\times(0\;1)$ . It satisfies homogeneous Dirichlet type boundary conditions for all $t\in(0\,\,\,1\]$ , which is formally written as $u(x,y,t)=0,\;\,(x,y)\in\partial\Omega$ . The initial conditions are prescribed such that $u(x,y,0)=1,\;\,(x,y)\in\Omega,\;\,t=0$ .
We present the code in the form of a serie of enumerated documentations for each block of code:

Variable
```
 L = 1
```
stores the value for the side length of the two-dimensional square domain $\Omega$ .
Variables
```
 tm = 1
 dt = 0.001
```
store the values respectively for the final time $T_{max}$ and the time step-size $\Delta\,t$ .
The number of time points $T$ on the time interval $(0\;1\]$ are stored by
```
 M = tm/dt
```
The number of spatial uniform mesh points that discretises L is stored by
```
 N = 50
```
The side length L is discretised by N+1 equally spaced points using
```
 X = seq(0,L,len=N+1)
```
The time interval is also uniformly discretised by M+1 points using
```
 T = seq(0,tm,len=M+1)
```
This step is to compute and store the matrices containing the values for the two dimensional spatial coordinates of such discretisation, which is achieved by
```
 m = length(X); n=length(X);
 x = matrix(rep(X,each=n),nrow=n);
 y = matrix(rep(X,m),nrow=n)
```
Reshape the data structure of the spatial coordinates from matrices into vector by
```
 x = c(x)
 y = c(y)
```
where x and y now store all the values for the x and y coordinates for all the global nodes, therefore, each one is now a vector of $(N+1)^2$ entries.
Store the number of global nodes in the discretised domain by
```
 GNodes = (N+1)^2
```
Use U to store the approximate discrete solution values at each time step and we define the initial state of U to be a vector of constant values of 1 at each note in the domain and this is achieved by
```
U = matrix(1, GNodes, 1)
```
Triangulation of a quadrilateral mesh by drawing a line of slope -1 diangonally through each square leads to the construction of a uniform triangulated domain $\Omega^h$ , which consists of NumTRI triangles, which is defined in terms of N by
```
NumTRI = 2*N^2
```
The connectivity array LocNodes is a matrix that stores the global counting of all the nodes with a local structure in the sense that it must have NumTRI rows and three columns for each vertix of each triangle. The code to achieve this is
```
LocNodes = matrix(0,NumTRI,3)
```
Triangulation of the domain is obtained such that the vertices of each triangle is locally counted in anti-clockwise orientation and such that the local counting of vertices of each triangle fills the global connectivity array LocNodes in the correct order. Note that we must fill in 4+2=6 vertices for each square, because each square is divided into two triangles by a diagonal line, however, each triangle has one distinct but two shared vertices. Therefore, it can be noted that the first three lines inside the nested for loop correspond to the three vertices of all the lower triangles on each square and the latter three lines correspond to all the three vertices of all the upper triangles. The nested loop that distributes such order of counting within the connectivity array is given by
```
for (i in 1:N){
    for (j in 1:N){
            LocNodes[i+2*(j-1)*N,1] = i+(j-1)*(N+1)
        LocNodes[i+2*(j-1)*N,2] = i+j*(N+1)
        LocNodes[i+2*(j-1)*N,3] = (i+1)+(j-1)*(N+1)
        LocNodes[i+N+2*(j-1)*N,1] = i+1+j*(N+1)
        LocNodes[i+N+2*(j-1)*N,2] = (i+1)+(j-1)*(N+1)
        LocNodes[i+N+2*(j-1)*N,3] = i+j*(N+1)
    }
}
```
Introduce the global stiffness and mass matrices filled with zeros each of size GNodes x GNodes
```
SP_Stiff <- matrix(0, GNodes, GNodes)
SP_Mass <- matrix(0, GNodes, GNodes)
```
The next segment of code is a for loop that executes a series of computations on each one of the triangles. It starts with
```
for (n in 1:NumTRI)
    {# This is the start of the for loop on all triangles
```

The first part within the loop defines the local position vectors r1, r2, and r3 in terms of x and y coordinate values through which we refer to different vertices of each triangle.

r1 = matrix(c(x[LocNodes[n,1]],y[LocNodes[n,1]]),nrow=2, byrow=FALSE)
r2 = matrix(c(x[LocNodes[n,2]],y[LocNodes[n,2]]),nrow=2, byrow=FALSE);
r3 = matrix(c(x[LocNodes[n,3]],y[LocNodes[n,3]]),nrow=2, byrow=FALSE);

The second part within the loop defines a jacobian matrix J for the mapping from an arbitrary triangle in to a reference triangle that has vertices in order and .
```
J = matrix(c(r2[1]-r1[1],r2[2]-r1[2]
    ,r3[1]-r1[1],r3[2]-r1[2]), nrow=2, byrow=TRUE)
```
The Jacobian of the mapping namely J serves to reduce the computational cost by a significant amount, particularly due to a property of integration for computing the integral of a function on a reference domain with a given mapping between the arbitrary domain and the reference domain. Further details on this topic can be found on Integral domain transformation.

Here is the animated solution of the diffusion equation .

…work in progress…