## 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 output the banner and get the program argument

 ``` banner(); if (argc <= 1) { cout << " Usage: lesson2 " << endl; exit(1); }```

• We construct an instance of class Mesh by using the file that contains mesh data in file. Note that to use fist a default constructor then invoke the member function get to read mesh data in the mesh file.

 ` Mesh ms(argv[1]);`

• The problem matrix is symmetric and will be stored in skyline format (thus using class SkSMatrix<double>)

 ` SkSMatrix A(ms);`

• Vectors b and bc (as instances of class Vect<double>) will contain the right-hand side and the solution, and imposed boundary conditions at nodes:

 ` Vect b(ms), bc(ms);`

• We assign imposed boundary conditions to vector bc by using member function setNodeBC which allows using an interpreted function of node coordinates: We prescribe the function y to nodes with code 1.

 ` bc.setNodeBC(1,"y");`

• 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 theElement enables access to current element data. The library defines the macro MeshElements(ms) that stands for a shorthand for the line

`      for (ms.topElement(); (theElement=ms.getElement());)`

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

 ``` MeshElements(ms) { DC2DT3 eq(theElement); eq.Diffusion(); eq.ElementAssembly(A); eq.ElementAssembly(b); }```

Of course, in the present case, assembling the right-hand side is 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(b,bc);`

• Solution is obtained by factorizing and backsubstituting:

 ``` A.solve(b); ```

Vector b contains now the solution.

• We finally output the solution and end the program.

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

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 test.m. Note, in this file, that a code equal to 1 is associated to nodes with y=0 and y=1.

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