Go back

Assembling the linear system

In this section, the user can learn how to assemble a Poisson problem by means of the Expression Template Assembler tool. It provides an easy way to build the linear system by simply writing an expression similar to the weak formulation of the problem.

First of all, we have to include:

1. the header file Integrate.hpp that contains the definitions of the classes and functions that actually allow to perform the integration;

2. the header file MatrixEpetra.hpp that defines the algebraic matrix.

   #!C++
   // ...
   
   // Header for the integrate expressions
   #include <lifev/eta/expression/Integrate.hpp>
   // Header for the matrices
   #include <lifev/core/array/MatrixEpetra.hpp>
   
   // ...
   

   We can also define the function that describes the source term of the Poisson problem:

   #!C++
   // ...
   
   // Raw function that describes the source term
   Real fSourceTerm ( Real /*t*/, Real /*x*/, Real /*y*/, Real /*z*/ , ID /*i*/ )
   {
       return 1.;
   }
   
   // ...
   

   Inside the main function, we introduce some aliases of the type we will use below.

   #!C++
   // ...
   int main ( int argc, char** argv )
   {
       // ...
   
       typedef MapEpetra map_Type;
       typedef std::shared_ptr<map_Type> mapPtr_Type;
       typedef MatrixEpetra<Real> matrix_Type;
       typedef std::shared_ptr<matrix_Type> matrixPtr_Type;
       typedef VectorEpetra vector_Type;
       typedef FESpace<mesh_Type, map_Type>::function_Type function_Type;
   
       // ...
   

   To assemble the linear matrix associated with the bilinear form of the Laplace operator, we have to:

3. create a matrix by specifying the map of the degrees of freedom;

4. initialize the matrix with zeros;
5. call the integrate() function by indicating i) the elements over which we want to integrate, ii) the quadrature rule, iii) the FE space of the test functions, iv) the FE space of the trial functions and v) the expression to integrate;
6. close the matrix.

   #!C++
       // ...
   
       // Creating the matrix for the system based on the Dofs' map
       matrixPtr_Type systemMatrix ( new matrix_Type ( uETFESpace->map() ) );
   
       // Initializing the matrix with zeros
       systemMatrix->zero();
   
       Real matrixNorm( systemMatrix->normInf() );
       if ( Comm->MyPID() == 0 )
       {
           // Matrix norm is set to -1 if the matrix it is not filled
           std::cout << "\nMatrix inf norm: " << matrixNorm << "\n\n";
       }
   
       // Assembling the stiffness matrix
       {
           using namespace ExpressionAssembly;
   
           // a( u_h, v_h ) = ( grad (u_h) , grad (v_h) )
           integrate ( elements ( uETFESpace->mesh() ),
                       uFESpace->qr(), // quadRuleTetra4pt,
                       uETFESpace,
                       uETFESpace,
                       dot ( grad (phi_j) , grad (phi_i)  )
                     ) >> systemMatrix;
       }
   
       // Closing the matrix: collect the repeated data and store in a unique way.
       systemMatrix->globalAssemble();
   
       matrixNorm = systemMatrix->normInf();
       if ( Comm->MyPID() == 0 )
       {
           std::cout << "\nMatrix inf norm: " << matrixNorm << "\n\n";
       }
   
       // ...
   

   Before assembling the right hand side, we have to interpolate the source function defined above into the discrete FE space. To do this, it is sufficient to create a vector using the map of the degrees of freedom and by calling the interpolate() method of the FESpace class, as follows:

   #!C++
       // ...
   
       // Creating the function for the Right Hand Side from the raw function
       function_Type fRhs (fSourceTerm);
   
       // Creating the vector that will contain the interpolated function f
       // on the Finite Element Space
       vector_Type fInterpolated (uETFESpace->map(), Repeated);
   
       // Initializing the vector with zeros
       fInterpolated.zero();
   
       // Interpolating the function f on the Finite Element Space
       uFESpace->interpolate (fRhs, fInterpolated, 0.0);
   
       // ...
   

Now we are ready to assemble the right hand side of the linear system. By following the same procedure used to build the matrix, we have to:

1. create a vector by specifying the map of the degrees of freedom;
2. initialize the vector with zeros;
3. call the integrate() function by indicating i) the elements over which we want to integrate, ii) the quadrature rule, iii) the FE space of the test functions and iv) the expression to integrate;
4. close the vector.

   #!C++
       // ...
   
       // Creating the vector for the Right Hand Side based on the Dofs' map
       vector_Type rhs (uETFESpace->map(), Repeated);
   
       // Initializing the vector with zeros
       rhs.zero();
   
       Real vectorNorm( rhs.normInf() );
       if ( Comm->MyPID() == 0 )
       {
           std::cout << "\nVector inf norm: " << vectorNorm << "\n\n";
       }
   
       // Assembling the RHS vector
       {
           using namespace ExpressionAssembly;
   
           // F (v_h) = ( f_h , v_h )
           integrate ( elements ( uETFESpace->mesh() ),
                       uFESpace->qr() ,
                       uETFESpace,
                       value (uETFESpace, fInterpolated) * phi_i
                     ) >> rhs;
       }
   
       // Closing the matrix: collect the repeated data and store in a unique way.
       rhs.globalAssemble();
   
       vectorNorm = rhs.normInf();
       if ( Comm->MyPID() == 0 )
       {
           std::cout << "\nVector inf norm: " << vectorNorm << "\n\n";
       }
   
       // ...
   }
   

Here you can download the complete main file mainAssembling.cpp, the mesh file cube.mesh and dataAssembling.in file.