main.cpp

OFELI is an object oriented library of C++ classes for development of finite element codes. Its main

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

                                *******************
                                *    S T D C 2    *
                                *******************


                        A Finite Element Code for Steady-State
                Analysis of Thermal Diffusion - Convection 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"

using namespace OFELI;


int main(int argc, char *argv[])
{
   Mesh     ms;
   Element  *el;
   Side     *sd;
   int      ret;
   FDF      *pl_file=NULL;
   ifstream mf, bcf, bodyf, boundf, vf;
   Point<double> vv(1,0,0);

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

   IPF data("stdc2 - 1.0",argv[1]);
   int output_flag = data.Output();
   int verbose = data.Verbose();
   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  C  2                  *\n";
     cout << "    *     Steady State Thermal Diffusion - Convection     *\n";
     cout << "    *                 in 2-D Geometries                   *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "               A Finite Element Code for Steady-State\n";
     cout << "           Analysis of Thermal Diffusion - Convection Problems\n";
     cout << "                             in 2-D Geometries\n\n";
     cout << "            STDC2 uses OFELI Library of Finite Element Classes\n\n";
     cout << "                           V E R S I O N   1.1\n\n";
     cout << "                     Copyright R. Touzani, 2003\n\n";
     cout << "=====================================================================\n\n";
   }

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

// Read Mesh data
   if (output_flag > 1)
     cout << "Reading mesh data ...\n";
   ms.Get(data.MeshFile());
   VDF vdf(ms,data.DataFile());
   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";
   SkMatrix<double> a(ms);
   Vect<double> u(ms.NbDOF());

// Read boundary conditions, body and boundary forces
/*   if (output_flag > 1)
     cout << "Reading boundary conditions ...\n";
   Vect<double> bc(ms.NbDOF());
   int bc_flag = data.BC();
   if (!bc_flag)
     ud.SetDBC(bc);
   else {
     FDF ff(data.BCFile(),"r");
     NodeVect<double> ui(ms,nb_dof);
     ff.Get(ui);
     bc = ui;
   }*/

// Read boundary conditions, body and boundary forces
   if (verbose > 1)
     cout << "Reading boundary conditions, body and boundary forces ...\n";
   Vect<double> bc(ms.NbNodes());
   Vect<double> 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);
   NodeVect<double> uf(ms,nb_dof);
   NodeVect<double> v(ms,ms.Dim());
   int vel_flag = data.IntPar(1);
   if (vel_flag) {
     if (output_flag > 1)
       cout << "Reading Velocity in file ...\n";
     FDF ff(data.AuxFile(1),"r");
     ff.Get(v);
   }

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

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

   for (ms.TopElement(); (el=ms.GetElement());) {
      DC2DT3 eq(el); 
      eq.Diffusion();

      if (vel_flag)
        eq.Convection(Vect<double>(v));
      else
        eq.Convection(vv);
         
      eq.BodyRHS(body_f);

      a.Assembly(el,eq.A());
      u.Assembly(el,eq.b());
   }

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

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

   for (ms.TopSide(); (sd=ms.GetSide());) {
      DC2DT3 eq(sd);
      eq.BoundaryRHS(Vect<double>(bound_f));
      u.Assembly(sd,eq.b());
   }

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

   if (output_flag > 1)
     cout << "Imposing boundary conditions ...\n";
   a.Prescribe(ms,u,bc);

   if (output_flag > 1)
     cout << "Solving linear system ...\n";
   ret = a.Factor();
   if (ret) {
     cout << "Error in matrix factorization : " << ret << "-th pivot";
     cout << " is too small = " << a(ret,ret) << endl;
     exit(1);
   }
   a.Solve(u);

   uf.FromVect(u,1,"Temperature",0);
   if (output_flag > 0)
     cout << uf;

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