## 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[]) {```

• We expand the argument of the program:

 ``` if (argc <= 1) { cout << " Usage: lesson3 " << endl; exit(1); }```

• We read mesh data using an option that removes imposed degrees of freedom for the list of unknowns. This is in contrast with the penalty technique used in Lesson 2.

 ` Mesh ms(argv[1],true);`

• We store the matrix in a sparse format, thus using class SpMatrix<double>

 ` SpMatrix 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. The vector bc will store imposed (Dirichlet) boundary conditions, and u is the final vector containing the solution. Note that the vectors bc and u are constructed using the mesh instance. This means that their size is the total number of degrees of freedom, while the other vectors have as size ms.getNbEq() which the number of equations, i.e. the number of degrees of freedom minus the number of prescribed boundary conditions.

 ` Vect b(ms.getNbEq()), x(ms.getNbEq()), bc(ms), u(ms);`

• Let us present a simple method to implement boundary conditions: We use for this the member function setNodeBC of class Vect<double>. This function enables imposing by means of an algebraic expression. Here, we impose on all nodes with code 1 the value y, which is the y-coordinate of the node.

 ` bc.setNodeBC(1,"y");`

• 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. Note that we have used a different member function for assembly than the one in Lesson 2: Here this one is a member function of the equation class. Its use is simpler since it uses as a unique argument either the matrix (for matrix assembly) or the right-hand side.

 ``` MeshElements(ms) { DC2DT3 eq(theElement); eq.Diffusion(); eq.updateBC(bc); eq.ElementAssembly(A); eq.ElementAssembly(b); }```

• In order to use an iterative method to solve the linear system use the class LinearSolver that drives iterative solvers. Here we use the constructor with the maximal number of iterations 1000, the tolerance for convergence criterium 1.e-8 and a verbosity parameter. We solve the linear system using the conjugate gradient method preconditioned with the diagonal ILU preconditioner.

 ``` LinearSolver ls(1000,1.e-8,1); int nb_it = ls.solve(A,b,x,CG_SOLVER,DILU_PREC);```

Note that the member function solve returns the number of performed iterations.

• We print out this number of iterations

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

• We incorporate boundary conditions into the solution vector:

 ` u.insertBC(ms,x,bc);`

• We finally output the solution vector and 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.