Heat Transfer Demo 2: A 2-D transient diffusion code

This code is developed to solve 2-D transient (time dependent) heat transfer problems. We use here an alternative to the presentation of the previous demo code. More specifically, we show how to perform all important phases of a time dependent finite element code.

• As usual, the program starts by including the header file OFELI.h that itself includes all kernel class definitions, then we include thermal dedicated classes described in Therm.h. We also declare the namespace to simplify the code writing:
  

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

• We construct an instance of class IPF to manage project data. We extract useful information from this instance: The maximal time, the time step, a flag for saving results and another for controlling output.

  IPF data("ttd2 - 1.1",argv[1]); double max_time = data.getMaxTime(); double deltat = data.getTimeStep(); int save_flag = data.getSave(); int verbose = data.getVerbose();

• We instantiate class IOField to store the solution for plotting:

  IOField pf(data.getPlotFile(),OUT);

• We now define the Mesh instance and the matrix (class SkSMatrix<double> handles skyline symmetric storage). We next instantiate vectors b and u that will contain respectively the right-hand side and the solution.

  Mesh ms(data.getMeshFile()); SkSMatrix<double> A(ms); Vect<double> b(ms.getNbDOF()), u(ms.getNbDOF());

• The initial solution is given through the class Prescription.

  Prescription pr(ms,data.getDataFile()); pr.get(INITIAL_FIELD,u);

• Boundary conditions and source vector will also be defined in the same way, but since all these properties are time dependent, we will do this each time step. For now, we just instantiate these vectors and define a NodeVect instance. In addition, we calculate the number of time steps:

  Vect<double> bc(ms.getNbDOF()), body_f(ms.getNbDOF()), bound_f(ms.getNbDOF()); int nb_step = int(max_time/deltat);

• We start the time step loop after initializing time value at 0:

  double time = 0; for (int step=1; step<=nb_step; step++) { time += deltat; b = 0;

• For the current time step, we get (Dirichlet) boundary conditions, sources and fluxes (Neumann boundary conditions):

  pr.get(BOUNDARY_CONDITION,bc,time); pr.get(SOURCE,body_f,time); pr.get(FLUX,bound_f,time);

• We define a loop over mesh elements to fill element arrays. We use here the macro MeshElementLoop that stands for a shorthand of the element loop. It uses the mesh instance and gives a value to the pointer el to class Element:

  Element *el; MeshElementLoop(ms,el) {

• For each element, we first define an instance of class DC2DT3 which handles diffusion-convection problems in 2-D domains using 3-node triangular finite elements (P₁ elements). We make use of the constructor for time dependent problems, this one involving the solution at previous time step and the current value of time. We then add the (lumped) capacity contribution to the element matrix and right-hand side. Note that this contribution is multiplied by the inverse of the time step. We also add a contribution of the diffusion matrix. All this corresponds to an implementation of the backward Euler scheme. We finally add the body right-hand side. Here the instance of Vect<double> corresponds contains the local element vector (local numbering of the nodes):

  DC2DT3 eq(el,u,time); eq.LCapacity(double(1./deltat)); eq.Diffusion(); eq.BodyRHS(Vect<double>(el,body_f));

• Once element arrays are computed, we assemble them to the global matrix and right-hand side. Note that since the matrix is computed and factorized only once, we avoid reassembling it.

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

• When Neumann boundary conditions are concerned (fluxes), we have to add a loop over sides. Such a loop is very similar to element loops.

  Side *sd; MeshSideLoop(ms,sd) { DC2DT3 eq(sd,u,time); eq.BoundaryRHS(Vect<double>(el,bound_f)); b.Assembly(sd,eq.b()); }

• Once the linear system is constructed, we impose (Dirichlet) boundary conditions by a penalty technique. This is done by the member function Prescribe that belongs to all matrix classes. We then solve the linear system by factorization and backsubstitution. The obtained solution is contained in b and transferred to vector u.

  A.Prescribe(ms,b,bc,step-1); if (step == 1) A.Factor(); A.Solve(b); u = b;

• It remains now to store the solution in a file if the parameter save_flag enables it. The procedure it to store the solution each multiple of save_flag step.

  NodeVect<double> U(ms,u,"Temperature",time); if (step%save_flag == 0) pf.put(U); }

• To end the program, we want to save the finite element mesh in a file to be treated by the software gnuplot. The utility function saveGnuplot can be used for this.

  saveGnuplot(data.getProject()+"-gpl1.dat",ms);

An example