main.cpp

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

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

                                *******************
                                *     S T D 3     *
                                *******************


                        A Finite Element Code for Steady-State
                        Analysis of Thermal Diffusion Problems
                                  in 3-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[])
{
   FDF           *pl_file=NULL;
   ifstream      mf, bcf, bodyf, boundf, vf;
   Point<double> vv(1,0,0);
   double        toler = 1.e-8;

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

   IPF data("std3 - 1.0",argv[1]);
   int output_flag = data.Output();
   int save_flag = data.Save();
   if (save_flag)
     pl_file = new FDF(data.SaveFile(),"w");

   if (output_flag) {
     cout << endl << endl;
     cout << "    *******************************************************\n";
     cout << "    *                         S  T  D  3                  *\n";
     cout << "    *               Steady State Thermal Diffusion     *\n";
     cout << "    *                 in 3-D Geometries                   *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "               A Finite Element Code for Steady-State\n";
     cout << "                Analysis of Thermal Diffusion Problems\n";
     cout << "                             in 3-D Geometries\n\n";
     cout << "            STD3 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, 2002\n\n";
     cout << "=====================================================================\n\n";
   }

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

// Read Mesh data
   if (output_flag > 1)
     cout << "Reading mesh data ...\n";
   Mesh ms(data.MeshFile());
   VDF vdf(ms,data.DataFile());
   ms.EliminateImposedDOF();
   ms.NumberEquations();
   int nb_dof = 1;
   if (output_flag > 1)
     cout << ms;

// Declare problem data (matrix, rhs, boundary conditions, body forces)
   if (output_flag > 1)
     cout << "Allocating memory for matrix and R.H.S. ...\n";
   SpMatrix<double> A(ms);
   Vect<double> u(ms.NbEq()), b(ms.NbEq());

// Read boundary conditions, body and boundary forces
   if (output_flag > 1)
     cout << "Reading boundary conditions, body and boundary sources ...\n";
   Vect<double> bc(ms.NbNodes()), body_f(ms.NbNodes());
   SideVect<double> bound_f(ms,nb_dof,"bound_f");
   vdf.Get(BOUNDARY_CONDITION,bc);
   vdf.Get(BODY_FORCE,body_f);
   vdf.Get(BOUNDARY_FORCE,bound_f,0);

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

   if (output_flag > 1)
     cout << "Looping over elements ...\n";
   Element *el;
   for (ms.TopElement(); (el=ms.GetElement());) {
      DC3DT4 eq(el);
      eq.Diffusion();
      eq.BodyRHS(body_f);
      eq.UpdateBC(bc);
      A.Assembly(el,eq.A());
      b.Assembly(el,eq.b());
   }

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

   if (output_flag > 1)
     cout << "Looping over sides ...\n";
   Side  *sd;
   for (ms.TopSide(); (sd=ms.GetSide());) {
      DC3DT4 eq(sd);
      eq.BoundaryRHS(bound_f);
      b.Assembly(sd,eq.b());
   }

   if (output_flag > 1)
     cout << "Solving linear system ...\n";
   ILUPrec<double,SpMatrix<double> > p(A);
   int nb_it = CG(A,p,b,u,1000,toler,0);
   cout << "Nb. of iterations : " << nb_it << endl;

// Output solution
   Vect<double> v(ms.NbDOF());
   v.InsertBC(ms,u,bc);

   NodeVect<double> uf(ms,v,"Temperature");
   if (output_flag > 0)
     cout << uf;

   if (save_flag) {
     pl_file->Put(uf);
     delete pl_file;
   }
   return 0;
}