## Lesson 4: A time-dependent problem

We introduce, in this lesson, new aspects of programmation:

• We consider a time-dependent heat transfer problem that we solve by Backward Euler time stepping scheme.
• We consider the case of Neumann Boundary conditions.
• Data (problem parameters) are introduced by a data file using IPF format.

The Finite Element Code

#### The main program

Here is a description of the source code. The source file can be found in Example 4 in the examples directory.
• We start, as usual, by including required headers and naming the appropriate namespace.

 ```#include "OFELI.h" #include "Therm.h" using namespace OFELI;```

• The program will have as argument the name of the parameter data file:

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

• We expand program arguments and declare an instance of class IPF for parameter file:

 ``` if (argc <= 1) { cout << " Usage: lesson4 " << endl; exit(1) } IPF data(argv[1]);```

• Time integration parameters are retrieved as IPF class members. These data are stored in global variables of which are theFinalTime and theTimeStep.

 ``` theFinalTime = data.getMaxTime(); theTimeStep = data.getTimeStep();```

• The mesh instance is constructed by giving the mesh file.

 ` Mesh ms(data.getMeshFile());`

• In the present example, we introduce boundary conditions through a user defined class. This may be optional for Dirichlet conditions but necessary for Neumann ones.

 ` User ud(ms);`

Implementation of class User will be given later.

• We declare matrix and vector data: first, the matrix A is declared as instance of class SkSMatrix. The vectors b, u and bc will contain respectively, alternatively the right-hand side and the current solution, the previous solution and prescribed Dirichlet boundary conditions.

 ``` SkSMatrix A(ms); Vect b(ms.getNbDOF()), u(ms.getNbDOF()), bc(ms.getNbDOF()); ```

• Since, we are dealing with a transient problem, we need initial data. This is retrieved from class member setInitialData of class User :

 ` ud.setInitialData(u);`

• We start a loop over time steps. For this we use a macro defined in : TimeLoop stands for a loop over time steps, that updates both time value (global variable theTime) and the step index (global variable theStep):

 ` TimeLoop {`

• The first thing to do here is to update time value and initialize the right-hand side to zero since this one will be assembled.

 ` b = 0;`

• We next transmit the user data class instance ud the time value:

 ` ud.setTime(theTime);`
• In order to deal with a problem with time-dependent boundary condition we re-fill vector bc at this level.

 ` ud.setDBC(bc);`

• We write a loop over finite elements as in the previous lessons:

 ` MeshElements(ms) {`

• We use here class DC2DT3 with the constructor that involves time. Instance eq will then be used to build matrix and right-hand side.

 ` DC2DT3 eq(theElement,u,theTime);`

• The element matrix is constructed with capacity term (chosen here to be lumped) and diffusion term:

 ``` eq.LCapacity(1./theTimeStep); eq.Diffusion();```

Note that capacity matrix is multiplied by the inverse of time step. This is necessary to implement the backward Euler scheme.

• We assemble matrix and right-hand side (useless for the present example). Note that, since the matrix does not depend on time, it is assembled once and factorized once.

 ``` if (theStep==1) eq.ElementAssembly(A); eq.ElementAssembly(b); }```

The loop on elements is closed.

• To deal with Neumann boundary conditions (involving boundary integrals), we have to loop over given sides. The loop looks like the one over elements:

 ` MeshSides(ms) {`

• For each side (pointed by eq) we invoke a constructor that involves sides.

 ` DC2DT3 eq(theSide,u,theTime);`

• We fill the side vector using the instance ud of class User. The function BoundaryRHS calculates the side integral.

 ` eq.BoundaryRHS(ud);`

• We assemble side vectors just like for elements and close the loop.

 ``` eq.SideAssembly(b); }```

• Once the linear system is assembled, we impose Dirichlet boundary conditions by a penalty techniques implemented in member function Prescribe:

 ` A.Prescribe(b,bc,theStep-1);`

• As said before, factorization is carried out at the first time step only. Obviously, solution is called each time step.

 ` A.solve(b);`

• Now, vector b contains the solution. We copy it to u to store it as a previous solution.

 ` u = b;`
• We may want to output the solution each time step:

 ``` cout << "\nSolution for time: " << time << endl << u; } return 0; }```

and then close the time stepping loop and the program.

#### A user defined class

We have now to implement class User that defines boundary conditions, initial conditions. The class is defined in file User.h in the Example 4 in the tutorial directory.
• Of course, we start by including file OFELI.h and invoking the namespace:

 ```#include "OFELI.h" using namespace OFELI;```

• Class User derives from abstract class UserData.

 `class User : public UserData {`

• This class has only public members and no attributes.

 `public :`

• We have a constructor that provides the mesh to the class : nothing to do, the parent class does the job for you.

 ` User(Mesh &mesh) : UserData(mesh) { }`

• We define a member class to give a value to prescribe for boundary condition in function of node code, node coordinates, time value and degree of freedom: Here, we impose that a code 2 imposes the value 1.0. Any other code will impose the default value 0.0.

 ``` double BoundaryCondition(const Point &x, int code, double time=0., size_t dof=1) { double ret = 0.0; if (code == 2) ret = 1.0; return ret; }```

• The same scheme works for Neumann boundary condition:

 ``` double SurfaceForce(const Point &x, int code, double time, size_t dof) { if (code) return 1.0; else return 0.0; } ```

• The class definition ends here

 `};`

• Let us finally note that since no implementation is given for initial condition, the default one is 0.0 for each degree of freedom.

A test

We use here exactly the same mesh file as in the previous lesson. Of course, we obtain the same solution. The convergence is obtained after 3 iterations.