Heat Transfer Demo 3: A 3-D steady state diffusion code

This code is developed to solve 3-D steady-state (stationary) heat transfer problems. The code is very similar to previous demo codes: stdc2 and ttd2.

  • We construct an instance of class IPF to manage project data. We extract useful information from this instance: An output flag to deal with output level and a save flag to say if the solution is to be saved in file. 

       IPF data("std3 - 1.0",argv[1]);
   int output_flag = data.getOutput();
   int save_flag = data.getSave();

  • After outputting some information about the code, we declare an instance of class Mesh and remove imposed boundary conditions: In fact, in view of using an iterative procedure to solve the linear system. Penalty techniques are not efficient in this case. 

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

  • In order to read various prescriptions (boundary conditions, heat sources, fluxes), we use the class Prescription. The name of the file containing these prescriptions is given in data.getDataFile() 

       Prescription p(ms,data.getDataFile());

  • The problem matrix is stored in sparse format in order to use an iterative method to solve the linear system. We declare an instance of class SpMatrix. We also declare vectors that will contain the solution and the right-hand side. 

       SpMatrix<double> A(ms);
   Vect<double> u(ms.getNbEq()), b(ms.getNbEq());

  • Now we can read various prescriptions: we first instantiate vectors bc and body_f that will contain respectively Dirichlet boundary conditions and sources. Neumann boundary conditions (imposed fluxes) are stored in a SideVect instance. Then we use the member function get with various arguments to fill properly these vectors. 

       Vect<double> bc(ms.getNbNodes()), body_f(ms.getNbNodes());
   SideVect<double> bound_f(ms,1,"bound_f");
   p.get(BOUNDARY_CONDITION,bc);
   p.get(BODY_FORCE,body_f);
   p.get(BOUNDARY_FORCE,bound_f,0);

  • We define a loop over elements to compute and assemble element arrays into the global matrix. Note the use of the macro MeshElements that updates the pointer theElement to current element at each step. For each element, we instantiate class DC3DT4 for thermal analysis in 3-D using 4-node tetrahedra, we compute the diffusion contribution, the source term and update the right-hand side with Dirichlet boundary conditions. Then we assemble element left-hand and right-hand sides. 

       MeshElements(ms) {
      DC3DT4 eq(theElement);
      eq.Diffusion();
      eq.BodyRHS(body_f);
      eq.updateBC(bc);
      eq.ElementAssembly(A);
      b.ElementAssembly(b);
   }

  • In this demo, we illustrate the use of side integrals that take account for Neumann (flux) boundary conditions. This works exactly like for elements. We use for this the macro MeshSides that updates the Side pointer theSide. 

       MeshSides(ms) {
      DC3DT4 eq(theSide);
      eq.BoundaryRHS(bound_f);
      eq.SideAssembly(b);
   }

  • We can now solve the linear system of equations using the Conjugate Gradient method preconditioned with the incomplete factorization preconditioner. The parameter toler is the tolerance value. On output it gives the actual discrepancy between the two last iterations. 

       double toler = 1.e-8;
   int nb_it = CG(A,Prec<double>(A,ILU_PREC),b,u,1000,toler,0);
   cout << "Nb. of iterations: " << nb_it << endl;

  • We want now to store the solution vector. Noting that the vector b does contain imposed boundary condition, we construct a Vect instance and put in it the obtained solution as well as imposed degrees of freedom by using the member function insertBC. 

       Vect<double> v(ms.getNbDOF());
   v.insertBC(ms,u,bc);

  • We next save the obtained vector, after transforming it to a NodeVect instance in the file data.getSaveFile() by using the class IOField. 

       if (save_flag) {
      IOField pf(data.getSaveFile(),OUT);
      pf.put(NodeVect<double>(ms,v,"Temperature"));
   }
      

    An example

    Let us run this program with the data presented in the following project file: 

    <?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
<info>
    <title></title>
    <date></date>
    <author></author>
</info>
<Project name="beam">
   <mesh_file>beam.m</mesh_file>
   <verbose>1</verbose>
   <output>0</output>
   <save>1</save>
</Project>
<Prescription>
   <BoundaryCondition code="1">1.0</BoundaryCondition>
</Prescription>
</OFELI_File>

    In summary, this file looks mainly like the one in the previous example.

     

    Copyright © 1998-2010 Rachid Touzani