Lesson 3

Lesson 3 : Using an iterative solver

We consider the same example as in Lesson 2 with the following modifications:

We solve the linear system using the Conjugate Gradient method.
In view of an iterative method we prefer to use, to handle boundary conditions, a classical substitution method rather than the penalty formulation.
We use a defined material by giving its name in the mesh file.

The Finite Element Code

The main program

Here is a description of the source code. The source file can be found in Example 3 in the examples directory.

We start like in Lesson 2.

#include "OFELI.h"
#include "Therm.h"
using namespace OFELI;

int main(int argc, char *argv[])
{

As usual, we declare an instance of class Meshspan>.

   Mesh ms(argv[1]);
   Element *el;
   banner();

We expand the argument of the program :

   if (argc <= 1) {
     cout << " Usage : ex3 <mesh_file>" << endl;
     exit(1);
   }

After reading mesh data we note that handling boundary conditions by elimination requires renumbering the equations. For this, we invoke the member class NumberEquations of class Mesh.
```
   ms.Get(argv[1]);
   ms.NumberEquations();
```
We store the matrix in a sparse format, thus using class SpMatrix<double>
```
   SpMatrix<double> a(ms);
```
Vectors b and x will store respectively the right-hand side and the solution. Note that, since imposed degrees of freedom are eliminated from the equations, the vectors have as sizes the actual number of equations.
```
   Vect<double> b(ms.getNbEq()), x(ms.getNbEq());
```
The vector bc, instance of class Vect<double> will store imposed boundary conditions. It is constructed in the same way as in Lesson 2.
```
   Vect<double> bc(ms.getNbDOF());
```
We construct the linear system of equations just as in Lesson 2. The difference here is that element right-hand sides need to be updated to take into account imposed boundary conditions at element level. This is necessary when using an elimination technique for boundary conditions. The member function updateBC is then used before assembly.
```
   for (ms.topElement(); (el=ms.getElement());) {
      DC2DT3 eq(el);
      eq.Diffusion();
      eq.updateBC(bc);
      a.Assembly(el,eq.A());
      b.Assembly(el,eq.b());
   }
```
We then execute the function CG to run the conjugate gradient invoking the ILU (Incomplete LU factorization) preconditioner. We impose a tolerance of 1.e-12.
```
   double toler = 1.e-8;
   int nb_it = CG(a,Prec<double>(a,ILU_PREC),b,x,100,toler,2);
```
Note that CG returns the number of performed iterations.

We print out this number of iterations

   cout << "Nb. of iterations: " << nb_it << endl;

We can incorporate boundary conditions into the solution vector.
```
   Vect<double> u(ms.getNbDOF());
   u.insertBC(ms,x,bc);
```

We finally end the program

   cout << "\nSolution :\n" << u;
   return 0;
}

How to declare a material ?

This is very simple: in the mesh data file, all elements have, in the present example, the code 1. If we say nothing, then a generic material (precisely called GenericMaterial with default properties is used. Otherwise, we can assign, in the mesh file a material, here Aluminium to this code, by the line

   Material  1  Aluminium

This line must be given after all elements lines. Note that the file Aluminium.md must be present in the material's directory. We are now ready to test the package.

A test

We use here a finer mesh than in Lesson 2 except the line defining the material. We obtain the same solution as Lesson 2 after 15 iterations.