Loading

We introduce, in this lesson, new aspects of OFELI 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

- 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 <parameter_file>" << endl; exit(1) } IPF data(argv[1]);

- Time integration parameters are retrieved as IPF class members. These data are stored in global variables of OFELI
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<double> A(ms); Vect<double> 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 OFELI: 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.

- 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<double> {

- 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<double>(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<double> &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<double> &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.

Copyright © 1998-2016 Rachid Touzani