Demo: A 1-D elliptic problem

This very simple code illustrates the numerical solution of a 1-D elliptic problem using P1 finite element method using the library .

  • All C++ programs using the library must include the header file OFELI.h that itself includes all kernel class definitions.

    #include "OFELI.h"

  • The library is embedded in namespace called OFELI.

    using namespace OFELI;

  • Our program has arguments that will described later.

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

  • L is the length of the interval in which the problem is defined. Here we have fixed their respective values at 0 and 1. N is the number of finite elements. Its default value is 10

       const double L = 1;
   int N = 10;

  • The function banner outputs the official banner of the library.

       banner();

  • N is the program argument (if this one is present) 

       if (argc > 1)
      N = atoi(argv[1]);

  • We now declare an instance of class Mesh with the appropriate constructor. Note this is not necessary in our example, but we use it here to present some associated functionalities like the assignment of vector components using algebraic expressions. 

       Mesh ms(L,N);

  • We denote by NbN the number of unknowns, the solution being prescribed at x=0 and then we print out the mesh. 

       int NbN = N+1;

  • In our demo programs, we have made use of exception handling. This is presented here for this first demo code. For this, the remaining of this program is contained in try block. 

       try {

  • We declare an instance of class TrMatrix<double> for the tridiagonal matrix, with size NbN and an instance of class Vect<double> for the right-hand side and the solution. Note that the vector constructor uses here the mesh instance. This simplifies the use of algebraic expressions for vector assignments. 

          TrMatrix<double> A(NbN);
      Vect<double> b(ms);

  • In order to test the code, we choose an exact solution: We set for right-hand side the function f(x) = 16π2sin(4πx) and the boundary conditions u(0) = u(1) = 0. This yields the solution u(x) = sin(4πx). To implement this without implementing additional functions we resort to the 's parser. For this, the member function set of class Vect<double> enables assigning regular expressions to nodes in function of their coordinates. The variables are x, y and z. 

          b.set("16*pi*pi*sin(4*pi*x)");

  • h is the mesh size (length of an element). 

          double h = L/double(N);

  • We now build up the matrix and the right-hand side. Note that we skip, for the moment, the first and the last lines for boundary condition treatment. The right-hand side vector is constructed by multiplying by the mesh size. 

          for (int i=2; i<NbN; i++) {
         A(i,i  ) =  2./h;
         A(i,i+1) = -1./h;
         A(i,i-1) = -1./h;
      }
      b *= h;

  • We modify the first and last equation in order to take account for boundary conditions. 

          A(1,1) = 1.;
      A(1,2) = 0.;
      b(1) = 0;
      A(NbN,NbN) = 1.;
      A(NbN-1,NbN) = 0.;
      b(NbN) = 0;

  • The linear system is solved. 

          A.solve(b);

  • We finally calculate the maximum-norm of the error and output it. The exact solution vector is defined here. 

          Vect<double> sol(ms);
      sol.set("sin(4*pi*x)");
      cout << "Max-Norm error = " << (b-sol).Norm(NORM_MAX) << endl;

  • We then close the block defining the exception region by invoking the macro that catches the exception. 

       } CATCH_EXCEPTION

  • We end the program. 

       return 0;
}