**2. Mesh generation**

After several steps, numerical solution of Eq. (1) is obtained by solving a linear

where A and *b* are called stiffness matrix and right hand side vector, respectively. In practice each term of diffusion, advection or reaction is accumulated

Several challenges make finding the numerical solution of (1) difficult. 3 of the most important are 1) complex geometry of Ω, 2) heterogeneity and high contrast in coefficient *κ* that produce very ill-conditioned linear systems and 3) large scale domains. A very common numerical method to solve (1) is FEM that uses a mesh representing the complex geometries accurately and overcomes the first issue. We can employ powerful linear solvers such as multigrid or multiscale methods to resolve the effect of heterogeneity and high contrast in *κ* in the linear system. Finally a careful implementation in parallel machines reduces computational time of large scale simulations considerably. However, lots of effort and research are still needed to propose a method obtaining a reliable solution in a reasonable time for real problems. Other issues such as uncertainties in the data or nonlinear coupled

Considering a positive definite *κ* we explain (and implement) how to obtain a

1. In Matlab we generate arbitrary 2D and simple 3D domains. The main domain is divided into a collection of subdomains, called elements where each element includes some nodes. Numbered nodes and elements generate a mesh. Actually by a mesh we mean 2 data structures, Cells (integer valued) that stores identifier or index of nodes in each element and Nodes (real valued) that stores coordinates of each node. Many software such as Matlab, GID and Gmesh generate reliable meshes. Corresponding to each element *e*, a positive definite matrix (a constant or a 2 � 2 matrix in 2D or a 3 � 3 matrix in 3D), named *κe*, will be set to form *κ*. In each element *e*, we can also set *μ<sup>e</sup>* and *γ<sup>e</sup>* to

form *μ* and *γ* in Eq. (1), respectively. See Section 2 for details.

considered for any type of elements exist in the mesh.

refer to [1] for further discussions and references.

accuracy and efficiency.

**98**

2.A brief introduction to FEM is presented that covers weak formulation (Section 3), shape functions and reference elements (Section 4).

4.Efficient implementation of boundary conditions, particularly Dirichlet boundary condition is presented in Section 5.2. We show that how a correct but careless implementation of Dirichlet boundary conditions decreases both

5.After assembly of the linear system (2), we use Matlab functions to solve the linear system and plot the result. Linear solvers are the core of a numerical simulator and their efficiency has a direct impact on the overall simulation. We

3.Evaluating of element integrals and preparing table of calculated integrals are discussed in Section 4.1. Then with help of affine mappings (Section 4.2) we use a linear combination of element integrals to accumulate into the stiffness matrix which is explained in Section 5. Note that different types of elements (for example triangles and rectangles in 2D) might exist in the mesh and we can have an unstructured mesh, generally. Hence element integrals should be

A*x* ¼ *b* (2)

system

separately in the stiffness matrix.

*Recent Advances in Numerical Simulations*

system of equations should be addressed, as well.

finite element solution of (1) through the following steps.

Implementation of FEM is began with mesh generation which is dividing the main domain into several subdomains or elements. Each element is finite (finite

**Figure 1.** *Creating a geometry and generating its mesh shown in Figure 2.*

ð Ω

*DOI: http://dx.doi.org/10.5772/intechopen.98291*

*<sup>κ</sup>*∇*<sup>u</sup>* � <sup>∇</sup>*<sup>v</sup>* � ∮ <sup>Γ</sup>

and our goal is to find unknown coefficients *u <sup>j</sup>*

*Nn*

*j*¼1 *u j* ð Ω

formula we obtain

*φj* n o *j*

write

X *Nn*

*j*¼1 *u j* ð Ω

form.

**101**

ð Ω

*<sup>κ</sup>*∇*<sup>φ</sup> <sup>j</sup>* � <sup>∇</sup>*φ<sup>i</sup>* <sup>þ</sup><sup>X</sup>

elements and then sum them up Ð

on the boundary of the elements.

**4. Shape functions**

∇ � �ð Þ *κ*∇*u v* ¼

ð Ω

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction…*

*ν* � *κ*∇*uv* þ

*<sup>u</sup>* <sup>¼</sup> <sup>X</sup> *Nn*

*j*¼1

setting *v* ¼ *φi*, *i* ¼ 1, ⋯, *Nn*, which is called standard Galerkin method we obtain

*Nn*

*j*¼1 *u j* ð Ω *μφ<sup>j</sup>*

Eq. (6) consists of evaluation of (from left to right) diffusion, advection and reaction terms, boundary integral and right hand side that we calculate them in

its gradient in boundary integral are known, hence we do not write its expansion

In practice definite integrals of basis functions and their derivatives in (6) are evaluated in an element, hence restriction of basis functions over elements, called shape functions, contribute in computations. With help of affine mapping we map a typical element into a fixed element, called reference element and evaluate integrals in it. So only shape functions in reference element determine values of integrals in (6). We use Lagrange elements which means corresponding to each node of an element a shape function is defined such that it is a polynomial in the element, takes value 1 at that specific node and 0 at other nodes of the element. Therefore a basis function in (5) becomes a continues function with value 1 at a specific node and 0 at neighboring nodes with trivial (zero) extension in the rest of the domain, as plotted in **Figure 3**. Note that the derivatives of Lagrange basis functions are not continues

Some examples of reference elements and their shape functions are presented in **Figure 4**. We see that in linear triangle element where nodes are the vertices of the triangle, the first shape function *S*1ð Þ¼ *r*, *s* 1 � *r* � *s* corresponding to the first node *n*<sup>1</sup> ¼ ð Þ 0, 0 has value 1 at *n*<sup>1</sup> and 0 at *n*<sup>2</sup> ¼ ð Þ 1, 0 and *n*<sup>3</sup> ¼ ð Þ 0, 1 . It is called a linear element since shape functions are combination of first order polynomials 1, f g *r*, *s* .

*<sup>γ</sup>* � <sup>∇</sup>*<sup>φ</sup> <sup>j</sup>φ<sup>i</sup>* <sup>þ</sup><sup>X</sup>

<sup>Ω</sup> <sup>¼</sup> <sup>P</sup> *e* Ð *e*

where *ν* is the (outward) normal vector. Regularity means that *u* and *v* have piecewise continuous (at least) first order partial derivatives to make the above integrals meaningful. So multiplying both side of (1) by *v* and applying Green's

> ð Ω

Approximation of functions in FEM is done by means of basis functions. Assume

, *j* ¼ 1, … , *Nn* is a set of basis functions, introduced in Section 4, that we can

*u <sup>j</sup>φ<sup>j</sup>*

� � *j*

*<sup>φ</sup><sup>i</sup>* <sup>þ</sup> ∮ <sup>Γ</sup>

� � to assemble the linear system (2). *u* or

*<sup>κ</sup>*∇*<sup>u</sup>* � <sup>∇</sup>*<sup>v</sup>* � ∮ <sup>Γ</sup>

*γ* � ∇*uv* þ

ð Ω *μuv* ¼

*ν* � *κ*∇*uv* (3)

ð Ω

, (5)

. Substituting (5) in (4) and

*ν* � �ð Þ *κ*∇*u φ<sup>i</sup>* ¼

ð Ω *f φ<sup>i</sup>*

(6)

*f v:* (4)

**Figure 2.**

*A 3D domain consists of 3 stacked cylinders (left) and the generated mesh (middle) and the solution of Eq. (21) in right.*

elements) and has a regular shape such as triangle or cuboid. Each element has a specific type, such as linear or bilinear, which is defined by shape functions and explained in Section 4. We prefer to use Lagrange triangles and parallelograms in 2D and tetrahedron and cuboid in 3D space as elements. The reason of such selection is to avoid numerical integration to evaluate definite integrals (of Lagrange elements) arise in FEM.

The main output of the mesh generation is two sets, nodes and elements that form the mesh. Nodes are not necessarily the vertices of the elements, for example in linear Crouzeix-Raviart element nodes are placed at edge midpoints or in quadratic Lagrange (triangle) elements midpoints and vertices together form the nodes.

In presence of complex geometries mesh generation can be a challenge for software mesh generators. For not very complicated geometry, Matlab, with a good quality, can generate triangular and tetrahedron mesh for 2D and 3D domains, respectively. We give an example, as shown in Part 1 of **Figure 1**.

First with createpde (1) we create a raw model considering one partial differential equation. Then we import or create geometry for this model. For our purposes 3 stacked cylinders with radius 2 and heights 1, 3 and 2 would be fine as shown in **Figure 2**. pdegplot shows how Matlab has specified regions (cylinders) and faces by numbering them. Then generateMesh generates a mesh with linear (or quadratic by default) tetrahedron elements and can be viewed by pdeplot3D. Smaller value for Hmax gives a finer mesh. Setting *Nn* and *Ne* for number of nodes and elements, respectively, Nodes that stores 3 Cartesian coordinates of each node is a matrix of size 3 *Nn*. Since all elements are linear tetrahedron, Cells that stores index of nodes of each element, would be a 4 *Ne* matrix. For example the first element of the mesh is formed by 4 nodes, stored in Cells(:,1) and Nodes(:,Cells(:,1)) returns 3D coordinates of that 4 nodes. So we can simply traverse over nodes and elements by their index. Moreover, we can extract element indices of a region with findElements or node indices of a region or a face by findNodes, where are necessary to set boundary conditions. For example we find node indices of bottom and top faces of the domain in Lines 21–22, since we will set boundary conditions on them.

## **3. Weak form**

In a bounded domain Ω with Lipschitz boundary Γ, Green's formula says for two regular functions *u* and *v* we have

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction… DOI: http://dx.doi.org/10.5772/intechopen.98291*

$$\int\_{\Omega} \nabla \cdot (-\kappa \nabla u) v = \int\_{\Omega} \kappa \nabla u \cdot \nabla v - \oint\_{\Gamma} \nu \cdot \kappa \nabla u v \tag{3}$$

where *ν* is the (outward) normal vector. Regularity means that *u* and *v* have piecewise continuous (at least) first order partial derivatives to make the above integrals meaningful. So multiplying both side of (1) by *v* and applying Green's formula we obtain

$$\int\_{\Omega} \kappa \nabla u \cdot \nabla v - \oint\_{\Gamma} \nu \cdot \kappa \nabla u v + \int\_{\Omega} \gamma \cdot \nabla u v + \int\_{\Omega} \mu u v = \int\_{\Omega} f v. \tag{4}$$

Approximation of functions in FEM is done by means of basis functions. Assume *φj* n o *j* , *j* ¼ 1, … , *Nn* is a set of basis functions, introduced in Section 4, that we can write

$$
\mu = \sum\_{j=1}^{N\_u} \mu\_j \rho\_j,\tag{5}
$$

and our goal is to find unknown coefficients *u <sup>j</sup>* � � *j* . Substituting (5) in (4) and setting *v* ¼ *φi*, *i* ¼ 1, ⋯, *Nn*, which is called standard Galerkin method we obtain

$$\sum\_{j=1}^{N\_n} u\_j \left[ \begin{matrix} \kappa \nabla \rho\_j \cdot \nabla \rho\_i + \sum\_{j=1}^{N\_n} u\_j \end{matrix} \right]\_{\Omega} \gamma \cdot \nabla \rho\_j \rho\_i + \sum\_{j=1}^{N\_n} u\_j \left[ \begin{matrix} \mu \rho\_j \rho\_i + \oint\_{\Gamma} \nu \cdot (-\kappa \nabla u) \, \rho\_i = \int\_{\Omega} f \, \rho\_i \right] \right] \tag{6}$$

Eq. (6) consists of evaluation of (from left to right) diffusion, advection and reaction terms, boundary integral and right hand side that we calculate them in elements and then sum them up Ð <sup>Ω</sup> <sup>¼</sup> <sup>P</sup> *e* Ð *e* � � to assemble the linear system (2). *u* or its gradient in boundary integral are known, hence we do not write its expansion form.

### **4. Shape functions**

elements) and has a regular shape such as triangle or cuboid. Each element has a specific type, such as linear or bilinear, which is defined by shape functions and explained in Section 4. We prefer to use Lagrange triangles and parallelograms in 2D and tetrahedron and cuboid in 3D space as elements. The reason of such selection is to avoid numerical integration to evaluate definite integrals (of Lagrange

*A 3D domain consists of 3 stacked cylinders (left) and the generated mesh (middle) and the solution of Eq. (21)*

The main output of the mesh generation is two sets, nodes and elements that form the mesh. Nodes are not necessarily the vertices of the elements, for example in linear Crouzeix-Raviart element nodes are placed at edge midpoints or in quadratic Lagrange (triangle) elements midpoints and vertices together form the nodes. In presence of complex geometries mesh generation can be a challenge for software mesh generators. For not very complicated geometry, Matlab, with a good quality, can generate triangular and tetrahedron mesh for 2D and 3D domains,

First with createpde (1) we create a raw model considering one partial differential equation. Then we import or create geometry for this model. For our purposes 3 stacked cylinders with radius 2 and heights 1, 3 and 2 would be fine as shown in **Figure 2**. pdegplot shows how Matlab has specified regions (cylinders) and faces by numbering them. Then generateMesh generates a mesh with linear (or quadratic by default) tetrahedron elements and can be viewed by pdeplot3D. Smaller value for Hmax gives a finer mesh. Setting *Nn* and *Ne* for number of nodes and elements, respectively, Nodes that stores 3 Cartesian coordinates of each node is a matrix of size 3 *Nn*. Since all elements are linear tetrahedron, Cells that stores index of nodes of each element, would be a 4 *Ne* matrix. For example the first element of the mesh is formed by 4 nodes, stored in Cells(:,1) and Nodes(:,Cells(:,1)) returns 3D coordinates of that 4 nodes. So we can simply traverse over nodes and elements by their index. Moreover, we can extract element indices of a region with findElements or node indices of a region or a face by findNodes, where are necessary to set boundary conditions. For example we find node indices of bottom and top faces of the domain

In a bounded domain Ω with Lipschitz boundary Γ, Green's formula says for two

respectively. We give an example, as shown in Part 1 of **Figure 1**.

in Lines 21–22, since we will set boundary conditions on them.

elements) arise in FEM.

*Recent Advances in Numerical Simulations*

**Figure 2.**

*in right.*

**3. Weak form**

**100**

regular functions *u* and *v* we have

In practice definite integrals of basis functions and their derivatives in (6) are evaluated in an element, hence restriction of basis functions over elements, called shape functions, contribute in computations. With help of affine mapping we map a typical element into a fixed element, called reference element and evaluate integrals in it. So only shape functions in reference element determine values of integrals in (6).

We use Lagrange elements which means corresponding to each node of an element a shape function is defined such that it is a polynomial in the element, takes value 1 at that specific node and 0 at other nodes of the element. Therefore a basis function in (5) becomes a continues function with value 1 at a specific node and 0 at neighboring nodes with trivial (zero) extension in the rest of the domain, as plotted in **Figure 3**. Note that the derivatives of Lagrange basis functions are not continues on the boundary of the elements.

Some examples of reference elements and their shape functions are presented in **Figure 4**. We see that in linear triangle element where nodes are the vertices of the triangle, the first shape function *S*1ð Þ¼ *r*, *s* 1 � *r* � *s* corresponding to the first node *n*<sup>1</sup> ¼ ð Þ 0, 0 has value 1 at *n*<sup>1</sup> and 0 at *n*<sup>2</sup> ¼ ð Þ 1, 0 and *n*<sup>3</sup> ¼ ð Þ 0, 1 . It is called a linear element since shape functions are combination of first order polynomials 1, f g *r*, *s* .

passing *n*<sup>2</sup> and *n*<sup>6</sup> which is 0*:*5 � *r* � *s* and *n*3, *n*<sup>4</sup> and *n*<sup>5</sup> which is 1 � *r* � *s*. So *S*<sup>1</sup> has to be of the form *α*ð Þ 0*:*5 � *r* � *s* ð Þ 1 � *r* � *s* . Substituting *n*<sup>1</sup> in equation of *S*1, we can find *α* such that it gives 1 for *n*1. This element is called quadratic since its shape

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction…*

**Exercise 1** Find shape functions of linear tetrahedron as presented in **Figure 4**. **Exercise 2** Quadratic tetrahedron element has nodes at midpoint of edges of linear tetrahedron element, so it has 6 other nodes (numbered 5 to 10) in addition to vertices (numbered 1 to 4 similar to linear tetrahedron). Find the 10 shape functions. For example for node at origin we have *S*<sup>1</sup> ¼ 2 1ð Þ � *r* � *s* � *t* ð Þ 0*:*5 � *r* � *s* � *t*

Definite integrals of shape functions and their derivatives in reference elements,

Similarly we define *Drs* and *Dss*. In 3D space we also have to define *Drt*, *Dst* and *Dtt*. In **Figure 5** we show how to evaluate *Drs* for linear tetrahedron element. For

CCCA, *Dst* <sup>¼</sup> <sup>1</sup>

3010 �1 �4 �1101 1 0 �1 0 �30 1 30 �1 00000 0 0 000 00000 0 0 000 �40 0 0 4 4 0 �40 0 �1 0 �30 4 4 4 �4 0 �4 �10 1 0 4 0 8 �4 0 �8 1030 �4 �4 �4404 1 0 �1 0 �4 0 �8408 00000 0 0 000

Other element integrals such as *Dr* (for advection term) and *D* (for reaction or

ð ^*e* 6

0

BBB@

1

called element integrals, play an important role in assembly of stiffness matrix. Denoting the reference element by ^*e* which has *N*^*<sup>e</sup>* nodes and partial derivatives of

<sup>2</sup> � �. Note that we first numbered

*<sup>∂</sup><sup>r</sup>* we define the matrix *Drr* (for diffusion term)

1 00 �1 0 00 0 �100 1 0 00 0

*Si S <sup>j</sup>*, *i*, *j* ¼ 1, ⋯, *N*^*e:* (8)

1

CCCA*:*

1

CCCCCCCCCCCCCCCCCCA

*<sup>∂</sup>rSi <sup>∂</sup>rS <sup>j</sup>*, *<sup>i</sup>*, *<sup>j</sup>* <sup>¼</sup> 1, <sup>⋯</sup>, *<sup>N</sup>*^*e:* (7)

functions are linear combination of 1,*r*, *s*,*rs*,*r*2, *s*

or for node at 0, 0, 1 ð Þ we have *S*<sup>4</sup> ¼ 2*t t*ð Þ � 0*:*5 .

*Dij rr* ¼ ð ^*e*

1 0 �1 0 �10 1 0 0000 0000

**Exercise 3** For quadratic tetrahedron element find *Drs* as

right hand side terms) can be defined where their *ij*th entries are

*<sup>∂</sup>rSi <sup>S</sup> <sup>j</sup>*, *<sup>D</sup>ij* <sup>¼</sup>

example *Drs* and *Dst* for linear tetrahedron are

0

BBB@

0

BBBBBBBBBBBBBBBBBB@

*Dij <sup>r</sup>* ¼ ð ^*e*

**103**

shape functions by *∂rS* which means *<sup>∂</sup><sup>S</sup>*

**4.1 Element integrals**

such that its *ij*th entry is

*Drs* <sup>¼</sup> <sup>1</sup> 6

*Drs* <sup>¼</sup> <sup>1</sup> 30

nodes at vertices and then at edge midpoints.

*DOI: http://dx.doi.org/10.5772/intechopen.98291*

**Figure 3.**

*Left) a shape function of reference bilinear quadrilateral. Right) 2 basis functions where their restrictions in a triangle is linear and in a quadrilateral is bilinear.*

**Figure 4.** *Some of reference Lagrange elements and their shape functions.*

Starting from *n*<sup>1</sup> nodes are numbered counterclockwise and we consider it for all Lagrange elements.

For Lagrange elements it is easy to find shape functions. For example to find *S*<sup>1</sup> for bilinear quadrilateral element where nodes are vertices of the square, it suffices to consider the line passing *n*<sup>2</sup> and *n*<sup>3</sup> which is 1 � *r* multiplied by the line passing *n*<sup>3</sup> and *n*<sup>4</sup> which is 1 � *s*. So *S*<sup>1</sup> has to be of the form *α*ð Þ 1 � *r* ð Þ 1 � *s* which gives 0 for all nodes except *n*<sup>1</sup> ¼ �ð Þ 1, �1 . Substituting *n*<sup>1</sup> in equation of *S*1, we can find *α* such that it gives 1 for *n*1. This element is called bilinear since its shape functions are linear combination of 1, f g *r*, *s*,*rs* .

For quadratic triangle element where vertices and midpoints form nodes, for shape function *S*<sup>1</sup> corresponding to *n*1, we need the multiplication of two lines

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction… DOI: http://dx.doi.org/10.5772/intechopen.98291*

passing *n*<sup>2</sup> and *n*<sup>6</sup> which is 0*:*5 � *r* � *s* and *n*3, *n*<sup>4</sup> and *n*<sup>5</sup> which is 1 � *r* � *s*. So *S*<sup>1</sup> has to be of the form *α*ð Þ 0*:*5 � *r* � *s* ð Þ 1 � *r* � *s* . Substituting *n*<sup>1</sup> in equation of *S*1, we can find *α* such that it gives 1 for *n*1. This element is called quadratic since its shape functions are linear combination of 1,*r*, *s*,*rs*,*r*2, *s* <sup>2</sup> � �. Note that we first numbered nodes at vertices and then at edge midpoints.

**Exercise 1** Find shape functions of linear tetrahedron as presented in **Figure 4**.

**Exercise 2** Quadratic tetrahedron element has nodes at midpoint of edges of linear tetrahedron element, so it has 6 other nodes (numbered 5 to 10) in addition to vertices (numbered 1 to 4 similar to linear tetrahedron). Find the 10 shape functions. For example for node at origin we have *S*<sup>1</sup> ¼ 2 1ð Þ � *r* � *s* � *t* ð Þ 0*:*5 � *r* � *s* � *t* or for node at 0, 0, 1 ð Þ we have *S*<sup>4</sup> ¼ 2*t t*ð Þ � 0*:*5 .

## **4.1 Element integrals**

Definite integrals of shape functions and their derivatives in reference elements, called element integrals, play an important role in assembly of stiffness matrix. Denoting the reference element by ^*e* which has *N*^*<sup>e</sup>* nodes and partial derivatives of shape functions by *∂rS* which means *<sup>∂</sup><sup>S</sup> <sup>∂</sup><sup>r</sup>* we define the matrix *Drr* (for diffusion term) such that its *ij*th entry is

$$D\_{rr}^{\vec{\mu}} = \int\_{\hat{\imath}} \partial\_r \mathbb{S}\_i \, \, \partial\_r \mathbb{S}\_j, \quad i, j = 1, \cdots, N\_{\hat{\imath}}.\tag{7}$$

Similarly we define *Drs* and *Dss*. In 3D space we also have to define *Drt*, *Dst* and *Dtt*. In **Figure 5** we show how to evaluate *Drs* for linear tetrahedron element. For example *Drs* and *Dst* for linear tetrahedron are

$$D\_{\pi} = \frac{1}{6} \begin{pmatrix} 1 & 0 & -1 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad D\_{\text{if}} = \frac{1}{6} \begin{pmatrix} 1 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{pmatrix}.$$

**Exercise 3** For quadratic tetrahedron element find *Drs* as

$$D\_{\mathcal{H}} = \frac{1}{30} \begin{pmatrix} 3 & 0 & 1 & 0 & -1 & -4 & -1 & 1 & 0 & 1 \\ 1 & 0 & -1 & 0 & -3 & 0 & 1 & 3 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -4 & 0 & 0 & 0 & 4 & 4 & 0 & -4 & 0 & 0 \\ -1 & 0 & -3 & 0 & 4 & 4 & 4 & -4 & 0 & -4 \\ -1 & 0 & 1 & 0 & 4 & 0 & 8 & -4 & 0 & -8 \\ 1 & 0 & 3 & 0 & -4 & -4 & -4 & 4 & 0 & 4 \\ 1 & 0 & -1 & 0 & -4 & 0 & -8 & 4 & 0 & 8 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}$$

Other element integrals such as *Dr* (for advection term) and *D* (for reaction or right hand side terms) can be defined where their *ij*th entries are

$$D^{\vec{j}}\_r = \int\_{\vec{\varepsilon}} \partial\_r \mathbb{S}\_i \, \, \mathbb{S}\_j, \quad D^{\vec{j}} = \int\_{\vec{\varepsilon}} \mathbb{S}\_i \, \, \mathbb{S}\_j, \quad i, j = \mathbb{1}, \cdots, N\_{\vec{\varepsilon}}.\tag{8}$$

Starting from *n*<sup>1</sup> nodes are numbered counterclockwise and we consider it for all

*Left) a shape function of reference bilinear quadrilateral. Right) 2 basis functions where their restrictions in a*

For Lagrange elements it is easy to find shape functions. For example to find *S*<sup>1</sup> for bilinear quadrilateral element where nodes are vertices of the square, it suffices to consider the line passing *n*<sup>2</sup> and *n*<sup>3</sup> which is 1 � *r* multiplied by the line passing *n*<sup>3</sup> and *n*<sup>4</sup> which is 1 � *s*. So *S*<sup>1</sup> has to be of the form *α*ð Þ 1 � *r* ð Þ 1 � *s* which gives 0 for all nodes except *n*<sup>1</sup> ¼ �ð Þ 1, �1 . Substituting *n*<sup>1</sup> in equation of *S*1, we can find *α* such that it gives 1 for *n*1. This element is called bilinear since its shape functions are

For quadratic triangle element where vertices and midpoints form nodes, for shape function *S*<sup>1</sup> corresponding to *n*1, we need the multiplication of two lines

Lagrange elements.

**Figure 4.**

**102**

**Figure 3.**

*triangle is linear and in a quadrilateral is bilinear.*

*Recent Advances in Numerical Simulations*

linear combination of 1, f g *r*, *s*,*rs* .

*Some of reference Lagrange elements and their shape functions.*

#### **Figure 5.**

*Element integrals* D*rs for linear tetrahedron.*

Note that *Drr* and *D* are symmetric matrices while *Dr* is not. Moreover, since *Drs* ¼ *Dsr* we calculate only one of them. All element integrals are evaluated once and saved for future use. In **Figure 1** we load element integrals of linear tetrahedron that we have calculated and saved in LinearTetraElementIntegral.mat.

#### **4.2 Affine mapping on reference elements**

So far we introduced shape functions and calculated element integrals, all in reference element. Now we explain how to map an arbitrary element in physical space into the reference element and evaluate integrals in (6) only in terms of element integrals and without any numerical integration. Clearly change of variable is necessary to evaluate integrals when we use a mapping which is explained in next subsection.

Consider an arbitrary triangle, say *e*, in 2D space with vertices *vi* ¼ *xi*, *yi* , *<sup>i</sup>* <sup>¼</sup> 1, 2, 3, and the reference linear triangle, ^*e*. Set

$$T = \begin{pmatrix} \mathbf{x}\_2 - \mathbf{x}\_1 & \mathbf{x}\_3 - \mathbf{x}\_1 \\ \mathbf{y}\_2 - \mathbf{y}\_1 & \mathbf{y}\_3 - \mathbf{y}\_1 \end{pmatrix}, \quad \eta = \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix}, \quad \xi = \begin{pmatrix} r \\ \mathbf{s} \end{pmatrix}, \quad \eta\_o = \begin{pmatrix} \mathbf{x}\_1 \\ \mathbf{y}\_1 \end{pmatrix} \tag{9}$$

and define the affine (composition of a linear mapping and a translation) mapping ^*e* ! *e* with the rule

$$
\xi \to \eta = T\xi + \eta\_o. \tag{10}
$$

*ξ* ¼

ð *e*

*κe*∇*ηφ<sup>i</sup>* � ∇*ηφ<sup>j</sup> dη* ¼

from ∇, then we have

ð *e*

> ð *e*

1*:*

**105**

Setting matrix *Me* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup>

*<sup>κ</sup>*∇*φ<sup>i</sup>* � <sup>∇</sup>*φ<sup>j</sup> <sup>d</sup><sup>η</sup>* <sup>¼</sup> *<sup>M</sup>*11<sup>ð</sup>

**Exercise 7** Show that in 3D space,

where *<sup>M</sup>* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup> *<sup>κ</sup>tT*�*<sup>t</sup>*

where vector *<sup>M</sup>* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup>

*r s t* 1

CA ! *<sup>η</sup>* <sup>¼</sup>

*x y z*

1

CA <sup>¼</sup>

0

B@

*x*<sup>2</sup> � *x*<sup>1</sup> *x*<sup>3</sup> � *x*<sup>1</sup> *x*<sup>4</sup> � *x*<sup>1</sup> *y*<sup>2</sup> � *y*<sup>1</sup> *y*<sup>3</sup> � *y*<sup>1</sup> *y*<sup>4</sup> � *y*<sup>1</sup> *z*<sup>2</sup> � *z*<sup>1</sup> *z*<sup>3</sup> � *z*<sup>1</sup> *z*<sup>4</sup> � *z*<sup>1</sup>

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction…*

maps the reference linear and quadratic tetrahedron to an arbitrary tetrahedron.

Consider a scalar function (a typical basis function) *φ* ¼ *φ η*ð Þ, where *η* ¼ ð Þ *x*, *y* in

where *φ ξ* ^ð Þ¼ *φ η*ð Þ and *<sup>T</sup>*�*<sup>t</sup>* is the transpose of the inverse of *<sup>T</sup>*. Note that if *<sup>φ</sup>* is the restriction of a basis function in element *e*, then *φ*^ would be the corresponding

> ð ^*e* ð Þ <sup>∇</sup>*<sup>ξ</sup> Si <sup>t</sup>*

*<sup>∂</sup>rSi <sup>∂</sup>sS <sup>j</sup> <sup>d</sup><sup>ξ</sup>* <sup>þ</sup> *<sup>M</sup>*22<sup>ð</sup>

∇*ηφ<sup>j</sup> dη* ¼

*∂rSi ∂rS <sup>j</sup> dξ*

^*e*

¼ *M*11*Drr* þ ð Þ *M*<sup>12</sup> þ *M*<sup>21</sup> *Drs* þ *M*22*Dss*

*κ*∇*φ<sup>i</sup>* � ∇*φ <sup>j</sup> dη* ¼ *M*11*Drr* þ ð Þ *M*<sup>12</sup> þ *M*<sup>21</sup> *Drs* þ ð Þ *M*<sup>13</sup> þ *M*<sup>31</sup> *Drt*

þ*M*22*Dss* þ ð Þ *M*<sup>23</sup> þ *M*<sup>32</sup> *Dst*

ð *e*

2. for a vector *γ*, defined for each element *e* in 3D space

*γ:*

ð *e*

vector. The chain rule and applying the affine mapping (10) on *φ* give

1

*r s t*

, where *t* denotes the transpose of the

<sup>∇</sup>*ηφ* <sup>¼</sup> *<sup>T</sup>*�*<sup>t</sup>* <sup>∇</sup>*<sup>ξ</sup> <sup>φ</sup>*^ (14)

*T*�<sup>1</sup> *κt eT*�*<sup>t</sup>*

*eT*�*<sup>t</sup>* and dropping subscript *e* as well as *η* and *ξ*

^*e*

*φ<sup>i</sup> φ <sup>j</sup> dη* ¼ ∣*T*∣*D* (18)

*γ* � ∇*φ<sup>i</sup> φ<sup>j</sup> dη* ¼ *M*1*Dr* þ *M*2*Ds* þ *M*3*Dt* (19)

*∂sSi ∂sS <sup>j</sup> dξ*

1

CA <sup>þ</sup>

*x*1 *y*1 *z*1

∇*<sup>ξ</sup> S <sup>j</sup>* ∣*T*∣*dξ:* (15)

(16)

(17)

1

CA (13)

0

B@

0

B@

CA

0

*DOI: http://dx.doi.org/10.5772/intechopen.98291*

B@

**4.3 Change of variables in integrals**

2D space. Gradient of *<sup>φ</sup>* is <sup>∇</sup>*ηφ* <sup>¼</sup> *<sup>∂</sup>xφ*, *<sup>∂</sup>y<sup>φ</sup>* � �*<sup>t</sup>*

shape function in ^*e*. Recalling *dη* ¼ ∣*T*∣*dξ*, we can write

*κe*∇*ηφ<sup>i</sup>* � �*<sup>t</sup>*

*κt*

^*e*

þð Þ *<sup>M</sup>*<sup>12</sup> <sup>þ</sup> *<sup>M</sup>*<sup>21</sup> <sup>ð</sup>

þ*M*33*Dtt*,

. **Exercise 8** Referring to notation (8) show that

ð *e*

0

B@

We see that the affine mapping (10) maps the nodes of reference linear triangle, f g *ni <sup>i</sup>*¼1,2,3, to the vertices of the triangle,

$$v\_i = Tn\_i + v\_1, \quad i = 1, 2, 3. \tag{11}$$

**Exercise 4** Show that the affine mapping in (9) and (10) also maps the reference quadratic triangle to an arbitrary triangle with 6 nodes (3 at vertices and 3 at edges midpoint).

**Exercise 5** Show that the affine mapping,

$$\xi = \begin{pmatrix} r \\ s \end{pmatrix} \rightarrow \eta = \begin{pmatrix} x \\ y \end{pmatrix} = \frac{1}{2} \begin{pmatrix} \mathbf{x}\_2 - \mathbf{x}\_1 & \mathbf{x}\_4 - \mathbf{x}\_1 \\ \mathbf{y}\_2 - \mathbf{y}\_1 & \mathbf{y}\_4 - \mathbf{y}\_1 \end{pmatrix} \begin{pmatrix} r \\ s \end{pmatrix} + \frac{1}{2} \begin{pmatrix} \mathbf{x}\_1 + \mathbf{x}\_3 \\ \mathbf{y}\_1 + \mathbf{y}\_3 \end{pmatrix} \tag{12}$$

maps the reference bilinear quadrilateral to an arbitrary parallelogram. Remember in a parallelogram where *v*<sup>1</sup> and *v*<sup>3</sup> are opposite vertices, we have *x*<sup>1</sup> þ *x*<sup>3</sup> ¼ *x*<sup>2</sup> þ *x*<sup>4</sup> and *y*<sup>1</sup> þ *y*<sup>3</sup> ¼ *y*<sup>2</sup> þ *y*4.

**Exercise 6** Show that the affine mapping,

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction… DOI: http://dx.doi.org/10.5772/intechopen.98291*

$$\xi = \begin{pmatrix} r \\ s \\ t \end{pmatrix} \rightarrow \eta = \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \\ z \end{pmatrix} = \begin{pmatrix} \mathbf{x}\_2 - \mathbf{x}\_1 & \mathbf{x}\_3 - \mathbf{x}\_1 & \mathbf{x}\_4 - \mathbf{x}\_1 \\ \mathbf{y}\_2 - \mathbf{y}\_1 & \mathbf{y}\_3 - \mathbf{y}\_1 & \mathbf{y}\_4 - \mathbf{y}\_1 \\ \mathbf{z}\_2 - \mathbf{z}\_1 & \mathbf{z}\_3 - \mathbf{z}\_1 & \mathbf{z}\_4 - \mathbf{z}\_1 \end{pmatrix} \begin{pmatrix} r \\ s \\ t \end{pmatrix} + \begin{pmatrix} \mathbf{x}\_1 \\ \mathbf{y}\_1 \\ \mathbf{z}\_1 \end{pmatrix} \tag{13}$$

maps the reference linear and quadratic tetrahedron to an arbitrary tetrahedron.

#### **4.3 Change of variables in integrals**

Note that *Drr* and *D* are symmetric matrices while *Dr* is not. Moreover, since *Drs* ¼ *Dsr* we calculate only one of them. All element integrals are evaluated once and saved for future use. In **Figure 1** we load element integrals of linear tetrahedron

So far we introduced shape functions and calculated element integrals, all in reference element. Now we explain how to map an arbitrary element in physical space into the reference element and evaluate integrals in (6) only in terms of element integrals and without any numerical integration. Clearly change of variable is necessary to evaluate integrals when we use a mapping which is explained in next

Consider an arbitrary triangle, say *e*, in 2D space with vertices *vi* ¼ *xi*, *yi*

*y* 

and define the affine (composition of a linear mapping and a translation) map-

We see that the affine mapping (10) maps the nodes of reference linear triangle,

**Exercise 4** Show that the affine mapping in (9) and (10) also maps the reference quadratic triangle to an arbitrary triangle with 6 nodes (3 at vertices and 3 at edges

> *x*<sup>2</sup> � *x*<sup>1</sup> *x*<sup>4</sup> � *x*<sup>1</sup> *y*<sup>2</sup> � *y*<sup>1</sup> *y*<sup>4</sup> � *y*<sup>1</sup> *r*

maps the reference bilinear quadrilateral to an arbitrary parallelogram. Remem-

ber in a parallelogram where *v*<sup>1</sup> and *v*<sup>3</sup> are opposite vertices, we have *x*<sup>1</sup> þ *x*<sup>3</sup> ¼

, *<sup>ξ</sup>* <sup>¼</sup> *<sup>r</sup>*

*s* 

*ξ* ! *η* ¼ *Tξ* þ *ηo:* (10)

*vi* ¼ *T ni* þ *v*1, *i* ¼ 1, 2, 3*:* (11)

*s* 

þ 1 2

*x*<sup>1</sup> þ *x*<sup>3</sup> *y*<sup>1</sup> þ *y*<sup>3</sup> 

, *<sup>η</sup><sup>o</sup>* <sup>¼</sup> *<sup>x</sup>*<sup>1</sup>

*y*1 

, *<sup>η</sup>* <sup>¼</sup> *<sup>x</sup>*

, *<sup>i</sup>* <sup>¼</sup>

(9)

(12)

that we have calculated and saved in LinearTetraElementIntegral.mat.

**4.2 Affine mapping on reference elements**

*Element integrals* D*rs for linear tetrahedron.*

*Recent Advances in Numerical Simulations*

1, 2, 3, and the reference linear triangle, ^*e*. Set

*<sup>T</sup>* <sup>¼</sup> *<sup>x</sup>*<sup>2</sup> � *<sup>x</sup>*<sup>1</sup> *<sup>x</sup>*<sup>3</sup> � *<sup>x</sup>*<sup>1</sup> *y*<sup>2</sup> � *y*<sup>1</sup> *y*<sup>3</sup> � *y*<sup>1</sup> 

f g *ni <sup>i</sup>*¼1,2,3, to the vertices of the triangle,

! *<sup>η</sup>* <sup>¼</sup> *<sup>x</sup>*

*x*<sup>2</sup> þ *x*<sup>4</sup> and *y*<sup>1</sup> þ *y*<sup>3</sup> ¼ *y*<sup>2</sup> þ *y*4.

**Exercise 5** Show that the affine mapping,

*y* 

**Exercise 6** Show that the affine mapping,

¼ 1 2

ping ^*e* ! *e* with the rule

subsection.

**Figure 5.**

midpoint).

**104**

*<sup>ξ</sup>* <sup>¼</sup> *<sup>r</sup> s* 

Consider a scalar function (a typical basis function) *φ* ¼ *φ η*ð Þ, where *η* ¼ ð Þ *x*, *y* in 2D space. Gradient of *<sup>φ</sup>* is <sup>∇</sup>*ηφ* <sup>¼</sup> *<sup>∂</sup>xφ*, *<sup>∂</sup>y<sup>φ</sup>* � �*<sup>t</sup>* , where *t* denotes the transpose of the vector. The chain rule and applying the affine mapping (10) on *φ* give

$$
\nabla\_{\eta} \varrho = T^{-t} \,\, \nabla\_{\xi} \hat{\varrho} \,\, \tag{14}
$$

where *φ ξ* ^ð Þ¼ *φ η*ð Þ and *<sup>T</sup>*�*<sup>t</sup>* is the transpose of the inverse of *<sup>T</sup>*. Note that if *<sup>φ</sup>* is the restriction of a basis function in element *e*, then *φ*^ would be the corresponding shape function in ^*e*. Recalling *dη* ¼ ∣*T*∣*dξ*, we can write

$$\int\_{\boldsymbol{\varepsilon}} \boldsymbol{\kappa}\_{\boldsymbol{\varepsilon}} \nabla\_{\eta} \boldsymbol{\rho}\_{i} \cdot \nabla\_{\eta} \boldsymbol{\varrho}\_{j} \, d\boldsymbol{\eta} = \int\_{\boldsymbol{\varepsilon}} \left( \boldsymbol{\kappa}\_{\boldsymbol{\varepsilon}} \nabla\_{\eta} \boldsymbol{\varrho}\_{i} \right)^{t} \nabla\_{\eta} \boldsymbol{\varrho}\_{j} \, d\boldsymbol{\eta} = \int\_{\boldsymbol{\delta}} (\nabla\_{\boldsymbol{\xi}} \mathbb{S}\_{i})^{t} \boldsymbol{T}^{-1} \boldsymbol{\kappa}\_{\boldsymbol{\varepsilon}}^{t} \boldsymbol{T}^{-t} \nabla\_{\boldsymbol{\xi}} \mathbb{S}\_{j} \, | \, |\boldsymbol{T}| d\boldsymbol{\xi}. \tag{15}$$

Setting matrix *Me* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup> *κt eT*�*<sup>t</sup>* and dropping subscript *e* as well as *η* and *ξ* from ∇, then we have

$$\begin{aligned} \int\_{\epsilon} \kappa \nabla \boldsymbol{\rho}\_{i} \cdot \nabla \boldsymbol{\rho}\_{j} \, d\boldsymbol{\eta} &= \boldsymbol{M}\_{11} \int\_{\hat{\epsilon}} \boldsymbol{\partial\_{r}} \mathbf{S}\_{i} \cdot \boldsymbol{\partial\_{r}} \mathbf{S}\_{j} \, d\boldsymbol{\xi} \\ &+ (\boldsymbol{M}\_{12} + \boldsymbol{M}\_{21}) \int\_{\hat{\epsilon}} \boldsymbol{\partial\_{r}} \mathbf{S}\_{i} \cdot \boldsymbol{\partial\_{s}} \mathbf{S}\_{j} \, d\boldsymbol{\xi} + \boldsymbol{M}\_{22} \int\_{\hat{\epsilon}} \boldsymbol{\partial\_{s}} \mathbf{S}\_{i} \cdot \boldsymbol{\partial\_{s}} \mathbf{S}\_{j} \, d\boldsymbol{\xi} \\ &= \boldsymbol{M}\_{11} \boldsymbol{D}\_{rr} + (\boldsymbol{M}\_{12} + \boldsymbol{M}\_{21}) \boldsymbol{D\_{rs}} + \boldsymbol{M}\_{22} \boldsymbol{D\_{ss}} \end{aligned} \tag{16}$$

**Exercise 7** Show that in 3D space,

$$\begin{aligned} \int\_{\epsilon} \kappa \nabla \rho\_i \cdot \nabla \rho\_j \, d\eta &= M\_{11} D\_{rr} + (M\_{12} + M\_{21}) D\_{rt} + (M\_{13} + M\_{31}) D\_{rt} \\ &+ M\_{22} D\_{tt} + (M\_{23} + M\_{32}) D\_{st} \\ &+ M\_{33} D\_{tt}, \end{aligned} \tag{17}$$

where *<sup>M</sup>* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup> *<sup>κ</sup>tT*�*<sup>t</sup>* . **Exercise 8** Referring to notation (8) show that

$$\begin{aligned} \text{(1. } \qquad \qquad \qquad \qquad \qquad \int\_{\epsilon} \rho\_i \, \, \rho\_j \, d\eta = |T|D \tag{18} \end{aligned} \tag{18}$$

2. for a vector *γ*, defined for each element *e* in 3D space

$$\int\_{\varepsilon} \chi \cdot \nabla \rho\_i \, \, \rho\_j \, d\eta = M\_1 D\_r + M\_2 D\_s + M\_3 D\_t \tag{19}$$

where vector *<sup>M</sup>* <sup>¼</sup> <sup>∣</sup>*T*∣*T*�<sup>1</sup> *γ:*

## **5. Assembly of the linear system**

With help of element integrals, affine mapping and Eqs. (17)–(19) we can accumulate each term of (6) into the stiffness matrix of the linear system (2).

#### **5.1 Diffusion term**

Equation

$$\int\_{\Omega} \kappa \nabla \rho\_j \cdot \nabla \rho\_i = \sum\_{\epsilon} \int\_{\epsilon} \kappa\_{\epsilon} \nabla \rho\_j \cdot \nabla \rho\_i \tag{20}$$

suggests how to accumulate Ð <sup>Ω</sup>*κ*∇*u* � ∇*v* into the stiffness matrix:


**Exercise 9** Follow the above steps to accumulate advection and reaction terms into the stiffness matrix.

In **Figure 6** we implemented accumulation of diffusion and reaction terms for linear or quadratic tetrahedron. Note that other elements only have a different affine mapping *T* and dimension, if we had elements in 2D space. We assume that *κ* is a diagonal matrix for each element, so a matrix storage of size 3 � *Nc* can store *κ* of all elements. We also assume that reaction coefficient *μ* is constant, otherwise a vector of size *Nc* can store *μ* for all elements and should be an input argument. The evaluated integrals are saved such that the indices of row and column and value of the integral are set in IM, JM and DDvec (for diffusion term), respectively. So Line 33 of **Figure 1** is to assembly the stiffness matrix of the equation

$$\int\_{\mathfrak{Q}} \kappa \nabla u \cdot \nabla v + \mathbf{0}.\mathbf{1} \Big|\_{\mathfrak{Q}} u v = \mathbf{0} \tag{21}$$

reaction terms and right hand side function means that we have considered a pure Neumann problem and this is exactly what we do in practice. After that we add requested boundary conditions to the linear system of (2). For example if *g* is nonzero on part of the boundary, say <sup>Γ</sup>*N*, then *<sup>g</sup>*∮ *<sup>φ</sup><sup>i</sup>* is known for indices that lie on

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction…*

*DOI: http://dx.doi.org/10.5772/intechopen.98291*

Γ*<sup>N</sup>* and should be added to the *i*th entry of *b* in (2).

*Assembly of diffusion–reaction terms for (linear or quadratic) tetrahedron.*

**Figure 6.**

**107**

where *κ* is set in Lines 25–26. We also set boundary conditions at bottom and top faces of the boundary (with values 1 and 10, respectively) and explain how to implement it in next subsection to obtain the solution of the problem, plotted in **Figure 2**.

#### **5.2 Boundary integral**

Probably correct implementation of boundary integrals is the most important part of the FEM, since boundary conditions mainly determine the situation of physical problem. Neumann, Robin and Dirichlet are 3 types of boundary conditions that might be considered on different parts of the boundary. Although boundary conditions are set on the boundary, they only affect on boundary nodes.

**Neumann boundary condition.** *ν* � �ð Þ¼ *κ*∇*u g* is a Neumann boundary condition, so boundary integral in (6) becomes *<sup>g</sup>*∮ *<sup>φ</sup>i*. If *<sup>g</sup>* is zero, then nothing has to be done for the boundary integrals. In fact incorporating only diffusion–advection-

#### **Figure 6.**

**5. Assembly of the linear system**

*Recent Advances in Numerical Simulations*

suggests how to accumulate Ð

ð Ω

1. traverse elements and in each element map Ð

3.map the local matrix into the stiffness matrix.

33 of **Figure 1** is to assembly the stiffness matrix of the equation

*κ*∇*u* � ∇*v* þ 0*:*1

faces of the boundary (with values 1 and 10, respectively) and explain how to implement it in next subsection to obtain the solution of the problem, plotted in

ð Ω

where *κ* is set in Lines 25–26. We also set boundary conditions at bottom and top

Probably correct implementation of boundary integrals is the most important part of the FEM, since boundary conditions mainly determine the situation of physical problem. Neumann, Robin and Dirichlet are 3 types of boundary conditions that might be considered on different parts of the boundary. Although boundary conditions are set on the boundary, they only affect on boundary nodes. **Neumann boundary condition.** *ν* � �ð Þ¼ *κ*∇*u g* is a Neumann boundary condition, so boundary integral in (6) becomes *<sup>g</sup>*∮ *<sup>φ</sup>i*. If *<sup>g</sup>* is zero, then nothing has to be done for the boundary integrals. In fact incorporating only diffusion–advection-

ð Ω

**5.1 Diffusion term**

Equation

element

into the stiffness matrix.

**Figure 2**.

**106**

**5.2 Boundary integral**

With help of element integrals, affine mapping and Eqs. (17)–(19) we can accumulate each term of (6) into the stiffness matrix of the linear system (2).

*<sup>κ</sup>*∇*<sup>φ</sup> <sup>j</sup>* � <sup>∇</sup>*φ<sup>i</sup>* <sup>¼</sup> <sup>X</sup>

2.with help of element integrals and Eq. (17) evaluate the desired term

**Exercise 9** Follow the above steps to accumulate advection and reaction terms

In **Figure 6** we implemented accumulation of diffusion and reaction terms for linear or quadratic tetrahedron. Note that other elements only have a different affine mapping *T* and dimension, if we had elements in 2D space. We assume that *κ* is a diagonal matrix for each element, so a matrix storage of size 3 � *Nc* can store *κ* of all elements. We also assume that reaction coefficient *μ* is constant, otherwise a vector of size *Nc* can store *μ* for all elements and should be an input argument. The evaluated integrals are saved such that the indices of row and column and value of the integral are set in IM, JM and DDvec (for diffusion term), respectively. So Line

*e*

ð *e*

<sup>Ω</sup>*κ*∇*u* � ∇*v* into the stiffness matrix:

*e*

*κe*∇*φ<sup>j</sup>* � ∇*φ<sup>i</sup>* (20)

*κe*∇*φ<sup>i</sup>* � ∇*φ<sup>j</sup>* into the reference

*uv* ¼ 0 (21)

*Assembly of diffusion–reaction terms for (linear or quadratic) tetrahedron.*

reaction terms and right hand side function means that we have considered a pure Neumann problem and this is exactly what we do in practice. After that we add requested boundary conditions to the linear system of (2). For example if *g* is nonzero on part of the boundary, say <sup>Γ</sup>*N*, then *<sup>g</sup>*∮ *<sup>φ</sup><sup>i</sup>* is known for indices that lie on Γ*<sup>N</sup>* and should be added to the *i*th entry of *b* in (2).

**Dirichlet boundary condition.** *u* ¼ *g* is a Dirichlet boundary condition, which means value of *u* at some nodes is known. So if we decompose node indices into 2 sets, say *U* and *K*, then we can write

$$u = \sum\_{j \in U} u\_j \rho\_j + \sum\_{j \in K} u\_j \rho\_j,\tag{22}$$

On the other hand reduction in size of the linear system by using (24), has great advantage in applications that large number of nodes have Dirichlet (more precisely essential) boundary condition. For example, in pore scale modeling and to approximate the permeability of a sandstone model, we had to set Dirichlet boundary condition for 75% of the nodes. Since the model included almost 10 million elements, solving a linear system with 2.5 million unknowns gave a significant speed up in computational time, in addition to significant improvement in accuracy. **Exercise 10** Robin boundary condition is a combination of Neumann and Dirichlet conditions which states *u* � *βν* � *κ*∇*u* ¼ *g* on Γ*R*. *β* is a non-zero function (normally a constant number) on Γ*R*. Considering �*ν* � *κ*∇*u* ¼ ð Þ *g* � *u =β*, explain

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction…*

Function *f* in (1) is the source or sink term and hence normally is defined on a very small part of the domain. It can be defined as a constant number over an element or can be approximated by its nodal values in the form *<sup>f</sup>* <sup>¼</sup> <sup>P</sup> *<sup>f</sup> <sup>j</sup>*

*f <sup>j</sup>* ¼ 0. We prefer nodal value approximation of *f*, since in applications that need to solve parabolic equation, *f* has a value in all elements; see *q*½ � *<sup>k</sup>* in Eqs. (27) and (28).

In applications such as simulation of compressible flows we need to solve the

where *ρ* and *f* are nonlinear functions of *u* and boundary condition and initial

*<sup>ρ</sup>* <sup>þ</sup> <sup>∇</sup> � �ð Þ� *ρκ*∇*<sup>u</sup> f u*ð Þ� <sup>1</sup>

In Eq. (26), the unknown *u*-dependent variables *ρ* and *f* should be evaluated at

*<sup>n</sup>*þ1, while *ρ<sup>n</sup>* is known from the previous time step.

Newton–Raphson linearization is the most common method to solve the nonlinear Eq. (26) which generates a sequence *<sup>u</sup>*½ � *<sup>k</sup>* � �, *<sup>k</sup>* <sup>¼</sup> 0, 1, 2, … by solving the

<sup>∇</sup> � �*ρ*½ � *<sup>k</sup> <sup>κ</sup>*∇*<sup>u</sup>*

*<sup>κ</sup>*∇*u<sup>n</sup>*þ<sup>1</sup> � � <sup>¼</sup> *<sup>f</sup>*

<sup>þ</sup> <sup>∇</sup> � �*ρ<sup>n</sup>*þ<sup>1</sup>

value are provided. A fully implicit time discretization of (25) gives

þ ∇ � �ð Þ¼ *ρ*ð Þ *u κ*∇*u f u*ð Þ, in Ω (25)

*n*þ1

*<sup>ρ</sup><sup>n</sup>* <sup>¼</sup> <sup>0</sup>*:* (26)

Δ*t*

� � <sup>þ</sup> *<sup>μ</sup>*½ � *<sup>k</sup> <sup>u</sup>* <sup>¼</sup> *<sup>q</sup>*½ � *<sup>k</sup>* , (27)

*<sup>j</sup>*¼<sup>1</sup> *<sup>f</sup> <sup>j</sup>* Ð *φ<sup>j</sup>* and

*φ <sup>j</sup>φ<sup>i</sup>* with possibly too many

how to implement it.

**5.3 Right hand side term**

time dependent equation

which simplifies to

the current time *t*

following equation

where

**109**

therefore right hand side term becomes P*Nn*

*DOI: http://dx.doi.org/10.5772/intechopen.98291*

**6. Generalization to parabolic equation**

*<sup>∂</sup>ρ*ð Þ *<sup>u</sup> ∂t*

> *<sup>ρ</sup><sup>n</sup>*þ<sup>1</sup> � *<sup>ρ</sup><sup>n</sup>* Δ*t*

> > Δ*t*

*G u*ð Þ¼ <sup>1</sup>

where *U* and *K* include indices of unknown and known values of *u*, respectively. If with permutation vector ½ � *U*, *K* we reorder rows and columns of the linear system in (2), then we can write

$$\mathcal{A}u = \begin{pmatrix} \mathcal{A}\_{11} & \mathcal{A}\_{12} \\ \mathcal{A}\_{21} & \mathcal{A}\_{22} \end{pmatrix} \begin{pmatrix} u\_U \\ u\_K \end{pmatrix} = \begin{pmatrix} b\_U \\ b\_K \end{pmatrix} \tag{23}$$

The first line gives a smaller linear system

$$\mathcal{A}\_{11}\mathfrak{u}\_U = \mathfrak{b}\_U - \mathcal{A}\_{12}\mathfrak{u}\_K \tag{24}$$

which preserves the symmetry or positive definiteness of the original problem. In Lines 35–42 of **Figure 1** we showed how to implement it in Matlab. Setting values 1 and 10 for bottom and top faces, we expect a smooth solution starting from 1 at the bottom and reaching to 10 at the top as shown in **Figure 2**.

The second approach to set Dirichlet condition that does not use permutation of the linear system is setting A<sup>21</sup> and A<sup>22</sup> in (23) to zero and identity matrices, receptively and solving the modified linear system. Although this approach gives the correct solution theoretically, it is very error-prone numerically as shown in a test case in **Figure 7** and explained as follows.

A domain with cuboid elements is refined in the middle and along the *x*-axis to simulate a fracture in the domain. We solve ∇ � �ð Þ¼ *κ*∇*u* 0 where *κ* in fracture is 100 times larger than rest of the domain. Setting 10<sup>6</sup> and 5 � <sup>10</sup><sup>6</sup> as Dirichlet boundary condition on left and right faces (along *x*-axis) of the domain, respectively, a correct solution should be started from 10<sup>6</sup> and monotonically reached to <sup>5</sup> � <sup>10</sup>6. The linear system is symmetric positive definite and preconditioned conjugate gradient method is employed to solve it. The correct solution shown in left image of **Figure 7** is obtained by (24). However, the second approach gives the nonphysical (wrong) solution along the fracture as shown in right. Moreover, the correct solution is obtained 3 times faster than the wrong one, mainly because the second approach violates the symmetry (but preserves the positive definiteness) of the linear system, due to setting *A*<sup>21</sup> ¼ 0.

#### **Figure 7.**

*Left) Result of correct implementation of Dirchlet BC by using Eq. (24). Right) nonphysical (wrong) result due to using the modified linear system (the second approach to implement Dirchlet BC).*

*A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction… DOI: http://dx.doi.org/10.5772/intechopen.98291*

On the other hand reduction in size of the linear system by using (24), has great advantage in applications that large number of nodes have Dirichlet (more precisely essential) boundary condition. For example, in pore scale modeling and to approximate the permeability of a sandstone model, we had to set Dirichlet boundary condition for 75% of the nodes. Since the model included almost 10 million elements, solving a linear system with 2.5 million unknowns gave a significant speed up in computational time, in addition to significant improvement in accuracy.

**Exercise 10** Robin boundary condition is a combination of Neumann and Dirichlet conditions which states *u* � *βν* � *κ*∇*u* ¼ *g* on Γ*R*. *β* is a non-zero function (normally a constant number) on Γ*R*. Considering �*ν* � *κ*∇*u* ¼ ð Þ *g* � *u =β*, explain how to implement it.

## **5.3 Right hand side term**

**Dirichlet boundary condition.** *u* ¼ *g* is a Dirichlet boundary condition, which means value of *u* at some nodes is known. So if we decompose node indices into 2

*<sup>u</sup> <sup>j</sup>φ<sup>j</sup>* <sup>þ</sup> <sup>X</sup>

where *U* and *K* include indices of unknown and known values of *u*, respectively. If with permutation vector ½ � *U*, *K* we reorder rows and columns of the linear system

which preserves the symmetry or positive definiteness of the original problem. In Lines 35–42 of **Figure 1** we showed how to implement it in Matlab. Setting values 1 and 10 for bottom and top faces, we expect a smooth solution starting from 1 at

The second approach to set Dirichlet condition that does not use permutation of

A domain with cuboid elements is refined in the middle and along the *x*-axis to simulate a fracture in the domain. We solve ∇ � �ð Þ¼ *κ*∇*u* 0 where *κ* in fracture is 100 times larger than rest of the domain. Setting 10<sup>6</sup> and 5 � <sup>10</sup><sup>6</sup> as Dirichlet boundary condition on left and right faces (along *x*-axis) of the domain, respectively, a correct solution should be started from 10<sup>6</sup> and monotonically reached to <sup>5</sup> � <sup>10</sup>6. The linear system is symmetric positive definite and preconditioned conjugate gradient method is employed to solve it. The correct solution shown in left image of **Figure 7** is obtained by (24). However, the second approach gives the nonphysical (wrong) solution along the fracture as shown in right. Moreover, the correct solution is obtained 3 times faster than the wrong one, mainly because the second approach violates the symmetry (but preserves the positive definiteness) of

*Left) Result of correct implementation of Dirchlet BC by using Eq. (24). Right) nonphysical (wrong) result due*

*to using the modified linear system (the second approach to implement Dirchlet BC).*

the linear system is setting A<sup>21</sup> and A<sup>22</sup> in (23) to zero and identity matrices, receptively and solving the modified linear system. Although this approach gives the correct solution theoretically, it is very error-prone numerically as shown in a

*j* ∈*K u <sup>j</sup>φ<sup>j</sup>*

*uK* � � <sup>¼</sup> *bU bK* � �

A11*uU* ¼ *bU* � A12*uK* (24)

, (22)

(23)

*<sup>u</sup>* <sup>¼</sup> <sup>X</sup> *j*∈ *U*

<sup>A</sup>*<sup>u</sup>* <sup>¼</sup> <sup>A</sup><sup>11</sup> <sup>A</sup><sup>12</sup>

the bottom and reaching to 10 at the top as shown in **Figure 2**.

The first line gives a smaller linear system

test case in **Figure 7** and explained as follows.

the linear system, due to setting *A*<sup>21</sup> ¼ 0.

**Figure 7.**

**108**

A21 A22

� � *uU*

sets, say *U* and *K*, then we can write

*Recent Advances in Numerical Simulations*

in (2), then we can write

Function *f* in (1) is the source or sink term and hence normally is defined on a very small part of the domain. It can be defined as a constant number over an element or can be approximated by its nodal values in the form *<sup>f</sup>* <sup>¼</sup> <sup>P</sup> *<sup>f</sup> <sup>j</sup> φ<sup>j</sup>* and therefore right hand side term becomes P*Nn <sup>j</sup>*¼<sup>1</sup> *<sup>f</sup> <sup>j</sup>* Ð *φ <sup>j</sup>φ<sup>i</sup>* with possibly too many *f <sup>j</sup>* ¼ 0. We prefer nodal value approximation of *f*, since in applications that need to solve parabolic equation, *f* has a value in all elements; see *q*½ � *<sup>k</sup>* in Eqs. (27) and (28).

### **6. Generalization to parabolic equation**

In applications such as simulation of compressible flows we need to solve the time dependent equation

$$\frac{\partial \rho(u)}{\partial t} + \nabla \cdot (-\rho(u)\kappa \nabla u) = f(u), \quad \text{in } \Omega \tag{25}$$

where *ρ* and *f* are nonlinear functions of *u* and boundary condition and initial value are provided. A fully implicit time discretization of (25) gives

$$\frac{\rho^{n+1} - \rho^n}{\Delta t} + \nabla \cdot \left( -\rho^{n+1} \kappa \nabla u^{n+1} \right) = f^{n+1}$$

which simplifies to

$$G(u) = \frac{1}{\Delta t}\rho + \nabla \cdot (-\rho \kappa \nabla u) - f(u) - \frac{1}{\Delta t}\rho^\eta = 0. \tag{26}$$

In Eq. (26), the unknown *u*-dependent variables *ρ* and *f* should be evaluated at the current time *t <sup>n</sup>*þ1, while *ρ<sup>n</sup>* is known from the previous time step.

Newton–Raphson linearization is the most common method to solve the nonlinear Eq. (26) which generates a sequence *<sup>u</sup>*½ � *<sup>k</sup>* � �, *<sup>k</sup>* <sup>¼</sup> 0, 1, 2, … by solving the following equation

$$\nabla \cdot \left( -\rho^{[k]} \kappa \nabla u \right) + \mu^{[k]} u = q^{[k]},\tag{27}$$

where

$$\mu^{[k]} = \frac{\mathbf{1}}{\Delta t} \frac{d\rho}{du}\Big|\_{u^{[k]}} - \frac{df}{du}\Big|\_{u^{[k]}}, \quad q^{[k]} = \mu^{[k]} u^{[k]} + f^{[k]} - \frac{\mathbf{1}}{\Delta t} \rho^{n}.\tag{28}$$

Having bounded derivatives and a good starting point *u*½ � <sup>0</sup> (which is normally chosen *un*), the resulted sequence *u*½ � *<sup>k</sup>* in Eq. (27) converges quadratically to *un*þ1, the solution of Eq. (26). Note that (26) is a diffusion–reaction equation which is completely explained how to solve and (core) codes were provided, as well.
