Demo: Using an iterative solver for a PDE

In this example, we solve the Poisson equation (Laplace equation with a right-hand side) by the P1-finite element method and by using the Conjugate Gradient iterative method to solve the linear system.

  • As usual, we start by including the principal header file and the header file dealing with the Laplace equation and defining the namespace . We start the program by reading the argument as the mesh file's name. Note that in the mesh construction phase, the second argument indicates that the nodes with prescribed boundary conditions (here boundary nodes) are eliminates from the list of equations. This is recommended when iterative methods are used. 

    #include "OFELI.h"
#include "Laplace.h"
using namespace OFELI;
                          
int main(int argc, char *argv[])
{
   Mesh ms(argv[1],true);

  • We instantiate the matrix as a sparse matrix (instance of SpMatrix) and the vectors b, x, bc and f. Here the vectors b and x stand for the right-hand side and the solution of the system. They must be dimensioned as the total number of equations (prescribed degrees of freedom are removed). The vectors bc and f stand respectively for the vector containing imposed values at nodes and the one containing the source (right-hand side of the PDE). Theses vectors are constructed using the Mesh instance. To initialize f we use a regular expression. The vector bc is set to zero since we impose homogeneous Dirichlet boundary conditions. 

       SpMatrix<double> A(ms);
   Vect<double> b(ms.getNbEq()), x(ms.getNbEq()), bc(ms), f(ms);
   f = "exp(-20*(x*x+y*y))";
   bc = 0;

  • To construct the finite element equations, we construct a loop over elements, instantiate for each element an instance of class Laplace2DT3 that deals with the 2-D Laplace equation using P1--elements. The member functions LHS and BodyRHS compute the left and right-hand sides of the element equation. The function updateBC adds the contribution of the boundary condition to the right-hand side. Note that this is actually not necessary here since we enforce the value 0. We next assemble left and right-hand sides into the global matrix and right-hand side vector. 

       MeshElements(ms) {
      Laplace2DT3 eq(theElement);
      eq.LHS();
      eq.BodyRHS(f);
      eq.updateBC(bc);
      eq.ElementAssembly(A);
      eq.ElementAssembly(b);
   }

  • To solve the linear system, we use an instance of class LinearSolver where some parameters are provided. Here we assign at maximum 1000 iterations, a tolerance of 1.e-8 as convergence criterium and no verbosity. We next solve the linear system by using the Conjugate Gradient method with the Diagonal preconditioner. This member function returns the number of performed iterations. 

       LinearSolver<double> ls(1000,1.e-8,0);
   int nb_it = ls.solve(A,b,x,CG_SOLVER,DIAG_PREC);
   cout << "Number of iterations: " << nb_it << endl;

  • To end we construct the solution vector by inserting the prescribed values. This is useful for any post-processing purposes. 

       Vect<double> u(ms);
   u.insertBC(ms,x,bc);
   return 0;
}