main.cpp

来自「OFELI is an object oriented library of C」· C++ 代码 · 共 229 行

/*==============================================================================

                                *******************
                                *     T T D 2     *
                                *******************


                        A Finite Element Code for Transient
                       Analysis of Thermal Diffusion Problems
                                  in 2-D Geometries

  ------------------------------------------------------------------------------

   Copyright (C) 1998 - 2004 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free 
   Software Foundation; Version 2 of the License.

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
   details.

   You should have received a copy of the GNU General Public License 
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================*/


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

int main(int argc, char *argv[])
{
   Element  *el;
   Side     *sd;
   double   time = 0;
   int      step, nb_step;
   ifstream mf, inf, bcf, bodyf, boundf;
   FDF      *pl_file=NULL;

// Expand arguments
   if (argc < 2) {
     cout << "\nUsage:  ttd2  <parameter_file>\n";
     return 0;
   }

   IPF data("ttd2 - 1.0",argv[1]);
   double max_time = data.MaxTime();
   double deltat = data.TimeStep();
   int output_flag = data.Output();
   int save_flag = data.Save();
   int init_flag = data.Init();
   int bc_flag = data.BC();
   int body_flag = data.BF();
   int bound_flag = data.SF();
   int verbose = data.Verbose();

   if (save_flag)
     pl_file = new FDF(data.PlotFile(),"w");

   if (verbose) {
     cout << endl << endl;
     cout << "    *******************************************************\n";
     cout << "    *                        T  T  D  2                   *\n";
     cout << "    *    Transient Thermal Diffusion in 2-D Geometries    *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "                   A Finite Element Code for Transient\n";
     cout << "       Analysis of Thermal Diffusion Problems in 2-D Geometries\n\n";
     cout << "            ttd2 uses OFELI Library of Finite Element Classes\n\n";
     cout << "                           V E R S I O N   1.0\n\n";
     cout << "                     Copyright R. Touzani, 1998\n\n";
     cout << "=====================================================================\n\n";
   }

//-----------
// Read data
//-----------

// Read Mesh data
   if (verbose > 1)
     cout << "Reading mesh data ...\n";
   Mesh ms(data.MeshFile());
   int nb_dof = 1;

   if (verbose > 1)
     cout << ms;
   User ud(ms);

// Declare problem data (matrix, rhs, boundary conditions, body forces)
   if (verbose > 1)
     cout << "Allocating memory for matrix and R.H.S. ...\n";
   SkSMatrix<double> a(ms);
   Vect<double> b(ms.NbDOF());

   Vect<double> u(ms.NbDOF());
   NodeVect<double> ui(ms,nb_dof);
   if (!init_flag)
     ud.SetInitialData(u);
   else {
     FDF ff(data.InitFile(),"r");
     ff.Get(ui);
     u = ui;
   }

   Vect<double> bc(ms.NbDOF()), body_f(ms.NbDOF()), bound_f(ms.NbDOF());
   NodeVect<double> uf(ms,nb_dof);
   nb_step = int(max_time/deltat);

// ------------------------
// Loop over time steps
// ------------------------

   for (step=1; step<=nb_step; step++) {

      if (verbose > 0)
        cout << "Performing time step " << step << endl;

      time += deltat;
      ud.Time(time);
      b = 0;

//    Read boundary temperature
      if (!bc_flag)
        ud.SetDBC(bc);
      else {
        FDF ff(data.BCFile(),"r");
        ff.Get(ui);
        bc = ui;
      }

//   Read boundary fluxes
     if (bound_flag) {
       FDF ff(data.SFFile(),"r");
       ff.Get(ui);
       bound_f = ui;
     }

// Read body sources
   if (body_flag) {
     FDF ff(data.BFFile(),"r");
     ff.Get(ui);
     body_f = ui;
   }

//    Loop over elements
//    ------------------

      if (verbose > 1)
        cout << "Looping over elements ...\n";

      for (ms.TopElement(); (el=ms.GetElement());) {

         DC2DT3 eq(el,u,time); 
         eq.LCapacity((double)(1./deltat));
         eq.Diffusion();
         if (body_flag)
           eq.BodyRHS(body_f);
         else
           eq.BodyRHS(ud);
         if (step==1)
           a.Assembly(el,eq.A());
         b.Assembly(el,eq.b());
      }

//    Loop over sides
//    ---------------

      if (verbose > 1)
        cout << "Looping over sides ...\n";

      for (ms.TopSide(); (sd=ms.GetSide());) {
         DC2DT3 eq(sd,u,time);
         if (bound_flag)
           eq.BoundaryRHS(bound_f);
         else
           eq.BoundaryRHS(ud);
         b.Assembly(sd,eq.b());
      }

//    Take account for boundary conditions and solve system
//    -----------------------------------------------------

      if (verbose > 1)
        cout << "Imposing boundary conditions ...\n";
      a.Prescribe(ms,b,bc,step-1);

      if (verbose > 1)
        cout << "Solving linear system ...\n";
      if (step == 1)
        a.Factor();
      a.Solve(b);
      u = b;

      uf.FromVect(u,1,"Temperature",time);
      if (verbose > 1)
        cout << "\nSolution at time : " << time << endl << uf;
      if (step%save_flag == 0)
        pl_file->Put(uf);
   }

   if (save_flag)
     delete pl_file;

// Store mesh in a 'Gnuplot' file
   char file[40];
   strcpy(file,data.Project());
   MDF2Gnuplot(strcat(file,"-gpl1.dat"),ms);

// Store solution in a 'Gnuplot' file
   strcpy(file,data.Project());
   FDF gf(data.PlotFile(),"r");
   FDF2Gnuplot(strcat(file,"-gpl2.dat"), ms, gf, uf);

// Calculate error for the demo
   void error(double time, const Mesh &ms, const Vect<double> &u);
   error(time,ms,u);

   return 0;
}