Lesson 2

Lesson 2 : A 2-D steady state diffusion code

We consider here a 2-D steady state boundary value problem. We solve a Poisson equation (Diffusion) with "simple" data. Concerning boundary conditions, we impose a Dirichlet (essential) boundary condition on a portion of the domain and a homogeneous Neumann (natural) condition on the remaining boundary. Note that owing to the variational formulation the Neumann condition is implicit (We have nothing to do for it).

The Finite Element Code

Here is a description of the source code. The source file can be found Example 2 in the examples directory

As usual we start by including the principal header file and the file Therm.h that includes all classes related to heat transfer problems.
```
#include "OFELI.h"
#include "Therm.h"
using namespace OFELI;
```
Our program will have as argument the mesh file name.
```
int main(int argc, char *argv[])
{
```
We declare an instance of class Mesh and a pointer to class Element that will be invoked later
```
   Mesh ms;
   Element *el;
```

We output the OFELI banner and get the program argument

   banner();
   if (argc <= 1) {
      cout << " Usage : ex2 <mesh_file>" << endl;
      exit(1);
   }

We use the member function Get to read mesh data in file. Note that it is possible to directly invoke the mesh file in the constructor in order to avoid calling member function Get. This will be done in forthcoming lessons.
```
   ms.Get(argv[1]);
```
The problem matrix is symmetric and will be stored in skyline format (thus using class SkSMatrix<double>)
```
   SkSMatrix<double> a(ms);
```
The vector b (as instance of class Vect<double>) will contain the right-hand side and the solution :
```
   Vect<double> b(ms.getNbDOF());
```
The vector bc contains imposed boundary conditions at nodes. Here, we present a simple method to fill it : We define a function BC for this end.
```
   Vect<double> bc(ms.getNbDOF());
   void BC(const Mesh &ms, Vect<double> &bc);
   BC(ms,bc);
```
We now start building the linear system. For this, we have to implement a loop over all element meshes by using member functions topElement() and getElement(). The returned pointer el enables access to current element data. For each element, we construct an instance of class DC2DT3 for diffusion-convection problems in 2-D using 3-node triangles. The member function Diffusion calculates the contribution to element matrix of diffusion term. We then assemble matrix and right-hand side (which is 0 here). Note that element matrix and right-hand side are accessible via the functions A() and b() respectively.
```
   for (ms.topElement(); (el=ms.getElement());) {
      DC2DT3 eq(el);
      eq.Diffusion();
      a.Assembly(el,eq.A());
      b.Assembly(el,eq.b());
   }
```
Of course, in the present case, assembling the right-hand side is actually useless.
Once the linear system is assembled, Dirichlet boundary conditions are imposed by a penalty technique. This is implemented via the function Prescribe, member of all matrix classes.
```
   a.Prescribe(ms,b,bc);
```
Solution is obtained by factorizing and backsubstituting :
```
   a.Factor();
   a.Solve(b);
```
Vector b contains now the solution.

We finally output the solution and end the program.

   cout << "\nSolution :\n" << b;
   return 0;
}

Let us give now the implementation of function BC. A few things can be said about this: we simply say that nodes with code =1 have as imposed boundary condition the value 1.

   void BC(const Mesh &ms, Vect<double> &bc)
   {
      const Node *nd;
      for (ms.topNode(); (nd=ms.getNode());) {
         if (nd->getCode(1)==2)
            bc(nd->getLabel()) = 1.;
      }
   }

A finite element mesh

To test this program we use a finite element mesh of a rectangle [0,3]x[0,1]. The imposed boundary conditions are

   u(x,0)=0, u(x,1)=1, 0<x<3

Homogeneous Neumann boundary conditions are "imposed" on the portions x=0 and x=3. The solution is then

   u(x,y)=y

The mesh file is ex2.m. Note, in this file, that a code equal to 0 is associated to nodes with y=0 and 1 is associated to nodes with y=1. The line starting with BC gives the associated value to the code 1, the case of value 0 is by default.

You can now execute the code to obtain the exact solution.