Loading

## 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 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

 ``` 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.

 ` 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. 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. Note that the vector constructor uses here the mesh instance. This simplifies the use of algebraic expressions for vector assignments.

 ``` TrMatrix A(NbN); Vect 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). We can compute the right-hand side of the linear system by multiplying it by the mesh size.

 ``` double h = L/double(N); b *= h;```

• 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.

 ``` for (int i=2; i

• 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. Note that the exact solution vector use also the parser.

 ``` cout << "\nSolution :\n"<< b; Vect sol(ms); sol.set("sin(4*pi*x)"); cout << "Error = " << (b-sol).getNormMax() << endl;```

• We now end the program.

 ``` return 0; }```

• 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

Solution :
1       0.00000000e+00
2       1.08674760e+00
3       6.71646957e-01
4      -6.71646957e-01
5      -1.08674760e+00
6      -1.62832710e-15
7       1.08674760e+00
8       6.71646957e-01
9      -6.71646957e-01
10      -1.08674760e+00
11       0.00000000e+00

Error = 1.35691088e-01
```

Copyright © 1998-2016 Rachid Touzani