This repository contains an implementation of the two-dimensional finite element method in the Cartesian coordinate system.
The FEM is a numerical method that solves boundary value problems (i.e., it solves partial differential equations) by discretizing the area. This is achieved by the construction of a mesh of the object: the numerical domain for the solution, which has a finite number of points. The FEM formulation of a boundary value problem finally results in a system of algebraic equations. The simple equations that model these finite elements are then assembled into a larger system of equations that models the entire problem. The method approximates the unknown function over the domain.
We will consider a two-dimensional elliptic equation:
Using the definitions of gradient and divergence (in the case of two variables in the Cartesian coordinate system), the equation can be rewritten as follows:
This equation is set in some area Ω with boundary 
and has boundary conditions:
Our goal is to find u = u(x, y) (i.e., its numerical approximation). The first step is to find an arbitrary trial function which satisfies the given boundary conditions. Substituting it into the equation, we will calculate the residual (the difference between the chosen and theoretical functions) to estimate the accuracy of the approximation. Next, we need to find such a function among some set of functions that approximates the u in the best way, i.e. solve the problem of minimization of a functional.
The best method for minimizing the error in our case is the Galerkin method, which can be described by this equation:
where R0 is an error (residual) and Ψ is some arbitrary function coinciding with the trial function.
Let's take the left side of our elliptic equation as R0 and divide the integral by the sum of integrals:
The FEM proceeds from a different form of the DE, called the weak formulation. In order to get it, we will use integration by parts, and also take the boundary conditions into account. Let's use integration by parts in the first integral:
The integral on the boundary S is divided into three boundaries (S1, S2, S3) on which we have boundary conditions:
And then:
This is where the division into finite elements begins. We can now represent the function u as a decomposition of the basis functions ψ with weights q:
Now the approximation of the function u will be performed by docked local basis functions on finite elements, and the solution of the problem is a vector of weights (q1, q2, ... qn), which can be obtained by solving a system of linear equations.
Let's substitute the decomposition of the function u into the Galerkin equation and obtain the final expression for the linear system with boundary conditions:
We will represent the finite elements as rectangles, and as functions Ψ we will take bilinear functions. The area will be divided into subareas corresponding to each finite element:
Bilinear functions are the product of the following one-dimensional functions on the corresponding axes: 
So, there are four local basis functions (on one finite element):
Each local basis function equals 1 at one node of its finite element Ωk.
At other nodes and finite elements it equals 0. 
In addition to the bilinear functions, we similarly construct biquadratic functions on each finite element, which we will use to decompose the diffusion coefficient λ:
The corresponding one-dimensional functions:
Now we need to proceed to the solution on each finite element from which the solution of the whole problem will be assembled. For this, let's present the equation obtained by Galerkin's method as a sum of integrals over regions Ωk without taking into account boundary conditions, and let's also disclose the gradients:
The first term gives the stiffness matrix G, the second term gives the mass matrix M. The sum of these matrices is called a local matrix A=G+M and has dimension 4x4 (based on the number of nodes). The integral on the right side is called a local vector b.
It's time to build G, M, and b. Let's decompose the diffusion coefficient (i.e., construct an interpolating function) and derive the formula for calculating the elements of the matrix G:
In this expression λ1, λ2, ..., λ9 are the values of λ at the corresponding nodes of the finite element.
The matrix is too large to be shown here because each of its 16 elements has 9 terms.
For the matrix M, we replace the parameter γ by its mean value on each finite element Ωk:
For the vector b we will represent the right side f as an interpolating function:
Now, using the local matrices and the local vector, we need to assemble the global matrix (A') and the global vector of the right side (f). Consider the following area as an example:
This picture shows the finite element numbers and global nodes. To find the solution we need to calculate 5 local matrices and 5 local vectors. The global matrix will be assembled from the local matrices: the assembly takes place in accordance with the local node numbering on the finite elements. For the second finite element, for example, local number 1 corresponds to global number 3, local number 2 corresponds to global number 4, and so on. Let's add the first two local matrices to the global matrix:
The global vector can be found in the same way.
Since each of the elements is related to a limited number of neighboring elements, the system of linear algebraic equations has a sparse form. We will store the global matrix efficiently, in the sparse string format:
- ggl, ggu - arrays for non-zero elements of the lower and upper triangles
- diOrig - array for diagonal elements
- ig - array for indexes of elements from which a string of non-zero elements begins
- jg - array for the column/row numbers in which the lower/upper triangle element is located, respectively
This allows us to store only non-zero elements. Consider the matrix:
Its representation in the sparse string format:
Up to this point, boundary conditions have not been taken into account, so it's time to do it.
To account for natural boundary conditions, we form local matrices and local vectors that will be added to the linear system similarly to the local matrices and local vectors of finite elements. In our (two-dimensional) case, the natural boundary conditions are considered on edges (one-dimensional segments).
Let's denote by Г the edge of length h, and by (i, j) the edge coordinates. We will represent the parameters θ and u as linear basis expansions, the parameter β will be regarded as a constant. Only two basis functions (Ψ1 and Ψ2) are nonzero on each edge, so the corresponding expressions take the following form.
For boundary condition of the second type:
For boundary condition of the third type:
Taking the essential (first type) boundary conditions into account, which is carried out after the global matrix is assembled, is the follows. By going through all the nodes on the edge, we replace the corresponding elements of the diagonal of the global matrix with a number that is much larger than the rest of the matrix elements.
In the global vector, the element with the corresponding number is assigned the product of the large number and the analytic value of the function.
We will decompose the global matrix into lower and upper triangles:
using the formulas:
Lower-upper decomposition (or factorization) is a better way to implement Gaussian elimination, especially for repeated solving a number of equations with the same left-hand side. The original system can now be solved in two steps:
To solve the final system of linear equations, we will use a method called the Local Optimal Scheme (LOS) with incomplete factorization, which is a lesser-known alternative to the conjugate gradient method (CGM).
Incomplete factorization allows to solve the following equivalent system instead of solving LUq=f:
Now let's assume:
Then we start the iterative process:
The process ends when 
becomes small enough. Vector x is the solution (i.e., it contains the coefficients q1, q2, ..., qn).
- main.cpp
The main program module containing function calls, local matrices and local vectors construction, diffusion coefficient decomposition, global matrix assembling and linear system solution. - config.cpp
A module designed to configure the parameters of the area, the source function, the vector of the right side and the LOS parameters. This module also allocates memory, configures pointers, and builds the matrix portrait. - boundaries.cpp
The implementation of a sequential accounting of natural and essential boundary conditions by enumerating the given types of conditions on the boundaries of the domain. - gauss.cpp
The implementation of forward and backward Gaussian elimination, as well as global matrix decomposition. These functions allow the implementation of incomplete factorization for the LOS. - math.cpp
A module containing minor functions: multiplication of global matrix by vector and vector by vector (applied in LOS), sum of products of upper and lower triangles (applied in global matrix decomposition). - io.cpp
A module containing functions to read input data and build a table with the result of program execution on the given data. - common.h
A header file containing preprocessor directives, function declarations and global variables common to all modules. - size.txt
This file specifies the number of nodes on the x-axis and the number of nodes on the y-axis, separated by a space. - nodesX.txt, nodesY.txt
Files containing node coordinates in the X and Y axes, respectively. - edgeTypes.txt
A file containing four space-separated values (1/2/3) describing the type of condition on the boundaries of an area. The first value corresponds to the lower boundary, the second to the right, the third to the top, and the fourth to the left.
Since this is a Windows console application, you can simply clone the repository directly into your Visual Studio environment or use git clone https://github.com/lenferdetroud/finite-element-method.git. The application has no interface.