Implementation of finite element method (FEM) needs special cares, particularly for essential boundary conditions that have an important effect on symmetry and number of unknowns in the linear systems. Moreover, avoiding numerical integration and using (off-line) calculated element integrals decrease the computational cost significantly. In this chapter we briefly present theoretical topics of FEM. Instead we focus on what is important (and how) to carefully implement FEM for equations that can be the core of a numerical simulator for a diffusion–advection-reaction problem. We consider general 2D and 3D domains, high contrast and heterogeneous diffusion coefficients and generalize the method to nonlinear parabolic equations. Although we use Matlab codes to simplify the explanation of the proposed method, we have implemented it in C++ to reveal the efficiency and examples are presented to admit it.
- finite element method
- diffusion–advection-reaction equation
- boundary condition
In a domain in consider a partial differential equation of the form of
where the unknown function , reaction coefficient and the given function are scalar functions . Diffusion coefficient and the (convergence free) advection coefficient are defined in and give (normally piecewise constant) values in and , respectively. Eq. (1) is called a diffusion–advection-reaction equation and has immense applications in science and engineering. Transport of heat, momentum and energy, solid mechanics, CO2 sequestration, computational fluid dynamics and biophysics are a few number of research fields that as a part of the solution strategy need to solve (1) numerically.
After several steps, numerical solution of Eq. (1) is obtained by solving a linear system
where and are called stiffness matrix and right hand side vector, respectively. In practice each term of diffusion, advection or reaction is accumulated separately in the stiffness matrix.
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 system of equations should be addressed, as well.
Considering a positive definite we explain (and implement) how to obtain a finite element solution of (1) through the following steps.
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 , a positive definite matrix (a constant or a matrix in 2D or a matrix in 3D), named , will be set to form . In each element , we can also set and to form and in Eq. (1), respectively. See Section 2 for details.
A brief introduction to FEM is presented that covers weak formulation (Section 3), shape functions and reference elements (Section 4).
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 considered for any type of elements exist in the mesh.
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 accuracy and efficiency.
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 refer to  for further discussions and references.
Generalization to more complex equations such as Eq. (25) now becomes straightforward and solving strategy is given in Section 6.
We finish this chapter by introducing the necessary topics to have an independent and efficient simulator in Section 7. Most of the topics presented here are fully discussed in [1, 2, 3] which are of great value for further reading.
2. Mesh generation
Implementation of FEM is began with mesh generation which is dividing the main domain into several subdomains or elements. Each element is finite (finite 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.
3. Weak form
In a bounded domain with Lipschitz boundary , Green’s formula says for two regular functions and we have
where is the (outward) normal vector. Regularity means that and have piecewise continuous (at least) first order partial derivatives to make the above integrals meaningful. So multiplying both side of (1) by and applying Green’s formula we obtain
Approximation of functions in FEM is done by means of basis functions. Assume , is a set of basis functions, introduced in Section 4, that we can write
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 to assemble the linear system (2). or its gradient in boundary integral are known, hence we do not write its expansion form.
4. Shape functions
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 corresponding to the first node has value 1 at and 0 at and . It is called a linear element since shape functions are combination of first order polynomials . Starting from 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 for bilinear quadrilateral element where nodes are vertices of the square, it suffices to consider the line passing and which is multiplied by the line passing and which is . So has to be of the form which gives 0 for all nodes except . Substituting in equation of , we can find such that it gives 1 for . This element is called bilinear since its shape functions are linear combination of .
For quadratic triangle element where vertices and midpoints form nodes, for shape function corresponding to , we need the multiplication of two lines passing and which is and and which is . So has to be of the form . Substituting in equation of , we can find such that it gives 1 for . This element is called quadratic since its shape functions are linear combination of . Note that we first numbered nodes at vertices and then at edge midpoints.
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 which has nodes and partial derivatives of shape functions by which means we define the matrix (for diffusion term) such that its th entry is
Similarly we define and . In 3D space we also have to define and . In Figure 5 we show how to evaluate for linear tetrahedron element. For example and for linear tetrahedron are
Other element integrals such as (for advection term) and (for reaction or right hand side terms) can be defined where their th entries are
Note that and are symmetric matrices while is not. Moreover, since 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
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 , in 2D space with vertices , and the reference linear triangle, . Set
and define the affine (composition of a linear mapping and a translation) mapping with the rule
We see that the affine mapping (10) maps the nodes of reference linear triangle, , to the vertices of the triangle,
maps the reference bilinear quadrilateral to an arbitrary parallelogram. Remember in a parallelogram where and are opposite vertices, we have and .
maps the reference linear and quadratic tetrahedron to an arbitrary tetrahedron.
4.3 Change of variables in integrals
Consider a scalar function (a typical basis function) , where in 2D space. Gradient of is , where denotes the transpose of the vector. The chain rule and applying the affine mapping (10) on give
where and is the transpose of the inverse of . Note that if is the restriction of a basis function in element , then would be the corresponding shape function in . Recalling , we can write
Setting matrix and dropping subscript as well as and from , then we have
2. for a vector , defined for each element in 3D space
5. Assembly of the linear system
5.1 Diffusion term
suggests how to accumulate into the stiffness matrix:
traverse elements and in each element map into the reference element
with help of element integrals and Eq. (17) evaluate the desired term
map the local matrix 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 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 can store of all elements. We also assume that reaction coefficient is constant, otherwise a vector of size 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
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.
where and include indices of unknown and known values of , respectively. If with permutation vector we reorder rows and columns of the linear system in (2), then we can write
The first line gives a smaller 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 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 and 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 -axis to simulate a fracture in the domain. We solve where in fracture is 100 times larger than rest of the domain. Setting and as Dirichlet boundary condition on left and right faces (along -axis) of the domain, respectively, a correct solution should be started from and monotonically reached to . 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 .
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 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.
5.3 Right hand side term
Function 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 and therefore right hand side term becomes with possibly too many . We prefer nodal value approximation of , since in applications that need to solve parabolic equation, has a value in all elements; see 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
where and are nonlinear functions of and boundary condition and initial value are provided. A fully implicit time discretization of (25) gives
which simplifies to
In Eq. (26), the unknown -dependent variables and should be evaluated at the current time , while 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 , by solving the following equation
Having bounded derivatives and a good starting point (which is normally chosen ), the resulted sequence in Eq. (27) converges quadratically to , 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.
7. Conclusion and future works
In this introductory chapter we presented the framework of FEM in brief (but effective) and implemented a numerical solver for diffusion–advection-reaction equation. Accumulation of different terms and setting boundary conditions correctly as well as evaluating of definite integrals without numerical integration were explained and their codes were also given. Finally we showed that nonlinear parabolic equations can be solved by combining of Newton method and diffusion–reaction equation where the latter was explained in this chapter.
However, to have a reliable and efficient simulator several topics and challenges should be addressed. Assuming correct mathematical model, mesh generation and assigning (measured or suggested) values to the coefficients of the problem are very important to make close the numerical model and its solution to the real problem [4, 5]. In a real problem, 2D and 3D elements appear together. For example 2D elements are employed to model faults in a 3D reservoir. Moreover, the generated mesh is normally unstructured and different types of elements are included in the mesh. Hence to assembly the linear system, traversing elements and nodes should have been implemented very efficiently.
Solving linear systems is normally the most time consuming part of the simulation and efficient implementation of linear solvers particularly in parallel machines, are of great interest . Multigrid and multiscale methods, particularly their algebraic form, have attracted interest to solve large scale linear systems since they have shown good scalability in addition to resolving low frequency parts of the error in solving linear systems [7, 8]. Standard Galerkin method that we used in this chapter is not a conservative scheme hence modifications or other methods such as discontinues Galerkin method would be necessary to solve problems in computational fluid dynamics [9, 10]. Proposing and implementing of advanced numerical algorithms to linearize and solve a system of nonlinear coupled initial-boundary value problem in a large scale domain become necessary. At last comparison with real data (if available) and quantifying uncertainties might force us to revisit all steps of the simulation again.
Hani Akbari would like to thank Prof. Lars K. Nielsen, scientific director at Novo Nordisk Foundation Center for Bio-Sustainability. Hani Akbari is grateful to DTU-Biosustain since this work was completed when he was a postdoc at DTU-Biosustain.