Lesson 1

Lesson 1 : A 1-D problem

This lesson concerns a simple one-dimensional two-point boundary value problem. The whole program can be found in Example 1 in the examples directory. We will examine it here below line by line.

We start by including the header file OFELI.h that itself includes all kernel class definitions.
```
   #include "OFELI.h" 
```
The OFELI library is embedded in namespace called OFELI.
```
   using namespace OFELI;
```
Our program has arguments that will described later.
```
   int main(int argc, char *argv[])
   {
```
Lmin and Lmax are the ends 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
```
      double Lmin=0, Lmax=1;
      int N=10;
```
The functions f and error return respectively the right-hand side function and the error vector. Their implementation is given hereafter.
```
      double f(double x);
      double error(const Mesh &ms, const Vect<double> &u);
```
The OFELI 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.
```
   Mesh ms(Lmin,Lmax,N);
```
We denote by NbN the number of unknowns, the solution being prescribed at x=0 and then we print out the mesh. Note that member function setVerbose with a large value for argument enables printing the whole mesh information.
```
   int NbN = N+1;
   ms.setVerbose(10);
   cout << ms;
```
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.
```
   TrMatrix<double> a(NbN);
   Vect<double> b(NbN);
```
h is the mesh size (length of an element).
```
   double h = (Lmax-Lmin)/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. Here x is the i-th node coordinate. The function f is invoked to calculate the right-hand side.
```
   for (int i=2; i<NbN; i++) {
      double x = ms.getPtrNode(i)->getCoord(1);
      a(i,i) = 2./h;
      a(i,i+1) = -1./h;
      a(i,i-1) = -1./h;
      b(i) = f(x)*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 output the solution, calculate the error at each node and output it.

   cout << "\nSolution :\n" << b;
   cout << "Error = " << error(ms,b) << endl;

We can now end the program.
```
   return 0;
}
```
We give now an example of right-hand side. Namely if we set f(x)=2*a*(1-2*a*x*x)*exp(-a*x*x) then the solution u that satisfies u(0)=u(1)=0 is given by u(x) = exp(-a*x*x) + (1-exp(-a))*x - 1
```
double f(double x)
{
   double a = 10;
   return 2*a*(1-2*a*x*x)*exp(-a*x*x);
}
```

Knowing the exact solution, we can compute the error via an appropriate function


double error(const Mesh &ms, const Vect<double> &u)
{
   double err=0;
   double a = 10;
   for (int i=1; i<=u.getSize(); i++) {
      double x = ms.getPtrNode(i)->getCoord(1);
      err = max(err,fabs(u(i) - exp(-a*x*x) - x*(1-exp(-a)) + 1));
   }
   return err;
}

If you execute the program without any argument, you will obtain as output :


M E S H     D A T A
===================

Space Dimension        :      1
Number of nodes        :     11
Number of elements     :     10
Number of sides        :      0


LIST OF ELEMENTS :

 Element :     1
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  1 2

 Element :     2
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  2 3

 Element :     3
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  3 4

 Element :     4
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  4 5

 Element :     5
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  5 6

 Element :     6
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  6 7

 Element :     7
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  7 8

 Element :     8
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  8 9

 Element :     9
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  9 10

 Element :    10
 Shape : line ** Material : GenericMaterial ** Code :     1 ** Nodes :  10 11



LIST OF NODES :

 Node :        1
 Coordinates : 0.000000e+000  0.000000e+000  0.000000e+000
  1 d.o.f. :        1
 Codes :          0

 Node :        2
 Coordinates : 1.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        2
 Codes :          0

 Node :        3
 Coordinates : 2.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        3
 Codes :          0

 Node :        4
 Coordinates : 3.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        4
 Codes :          0

 Node :        5
 Coordinates : 4.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        5
 Codes :          0

 Node :        6
 Coordinates : 5.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        6
 Codes :          0

 Node :        7
 Coordinates : 6.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        7
 Codes :          0

 Node :        8
 Coordinates : 7.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        8
 Codes :          0

 Node :        9
 Coordinates : 8.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :        9
 Codes :          0

 Node :       10
 Coordinates : 9.000000e-001  0.000000e+000  0.000000e+000
  1 d.o.f. :       10
 Codes :          0

 Node :       11
 Coordinates : 1.000000e+000  0.000000e+000  0.000000e+000
  1 d.o.f. :       11
 Codes :          0



Solution :
     1     0.00000000e+000
     2     1.71545445e-003
     3    -1.41343078e-001
     4    -3.11214412e-001
     5    -4.16034601e-001
     6    -4.32020322e-001
     7    -3.82338044e-001
     8    -2.98774350e-001
     9    -2.02104670e-001
    10    -1.01513714e-001
    11     0.00000000e+000
Error = 1.7912959e-002
Removing Mesh instance ...