main.cpp

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

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

                             *************************
                             *      L E L A S 3 D    *
                             *************************

                A Finite Element Code for 3-D Linearized Elastostatics in

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

   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 "Solid.h"
using namespace OFELI;

int main(int argc, char *argv[])
{
   double toler = 1.e-8;
   Element  *el;
   ifstream mf, bcf, bodyf, boundf;

   if (argc < 2) {
     cout << "\nUsage:  lelas3d  <parameter_file>\n";
     return 0;
   }

   IPF data("lelas3d - 1.0",argv[1]);
   int verbose = data.Verbose();
   int output_flag = data.Output();
   int save_flag = data.Save();
   int bc_flag = data.BC();
   int body_flag = data.BF();
   int bound_flag = data.SF();

   if (output_flag) {
     cout << endl << endl;
     cout << "    *******************************************************\n";
     cout << "    *                      L E L A S 3 D                  *\n";
     cout << "    *     Three-Dimensional Linearized Elastostatics      *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "         A Finite Element Code for Linearized Elastostatics\n";
     cout << "                             in 3-D Geometries\n\n";
     cout << "            lelas3d 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());
   ms.EliminateImposedDOF();
   ms.NumberEquations();
   int nb_dof = ms.Dim();
   if (output_flag > 1)
     cout << ms;
   VDF vdf(ms,data.DataFile());

// 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> b(ms.NbEq()), u(ms.NbEq());

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

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

   if (output_flag > 1)
     cout << "Looping over elements ...\n";
   for (ms.TopElement(); (el=ms.GetElement());) {
      Elas3DT4 eq(el);
      eq.Deviator();
      eq.Dilatation();
      eq.BodyRHS(body_f);
      eq.UpdateBC(BCVect<double>(bc));
      A.Assembly(el,eq.A());
      b.Assembly(el,eq.b());
   }

// Solve system
// ------------

   DiagPrec<double,SpMatrix<double> > p(A);
   int nb_it = CG(A,p,b,u,1000,toler,0);
   cout << "Nb. of iterations : " << nb_it << endl;

   NodeVect<double> uf(ms,u,BCVect<double>(bc),"Displacement");
   if (output_flag > 0)
     cout << uf;

   if (save_flag) {
     FDF pl_file(data.PlotFile(),FDF_WRITE);
     pl_file.Put(uf);
     DeformMesh(ms,uf,1);
     ms.Put(strcat(data.Project(),".1.m"));
   }
   return 0;
}