main.cpp

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

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

                             *************************
                             *     C O N T A C T     *
                             *************************


                            A Signorini Contact Problem
                      for Plane Strain Linearized Elasticity

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

   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[])
{
   Element  *el;
   Side     *sd;
   ifstream mf, bcf, bodyf, boundf;
   ofstream dmf;
   double   penal = 1.e07;
   double   err = 1.;

   cout << "\n\n";
   cout << "contact, version 1.1, Copyright (c) 1998 - 2004 by Rachid Touzani\n";
   cout << "contact comes with ABSOLUTELY NO WARRANTY.\n";
   cout << "This is free software, and your are allowed to redistribute it\n";
   cout << "under certain conditions. Details are distributed with the software." << endl;

   if (argc < 2) {
     cout << "\nUsage:  contact  <parameter_file>\n";
     return 0;
   }
   IPF data("contact - 1.0",argv[1]);
   int verbose = data.Verbose();

   if (verbose) {
     cout << endl << endl;
     cout << "    *********************************************************\n";
     cout << "    *                     C o n t a c t                     *\n";
     cout << "    * Two-Dimensional Linearized Elastostatics with contact *\n";
     cout << "    *********************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "         A Finite Element Code for Linearized Elastostatics\n";
     cout << "                      in 2-D Geometries with contact\n\n";
     cout << "          Contact 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 parameters and mesh data
   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();
   int it_max = data.NbIter();
   double tolerance = data.Tolerance();
   Mesh ms(data.MeshFile(1));

   VDF vdf(ms,data.DataFile());
   if (verbose > 1)
     cout << "Reading mesh data ...\n";
   int nb_dof = ms.Dim();
   int nn = ms.NbDOF();
   if (verbose > 1)
     cout << 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> u(nn), v(nn);

// Read boundary conditions, body and boundary forces
   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);
   NodeVect<double> uf(ms,u,2,1,"Displacement");

//---------------
// Iteration Loop
//---------------

   int it = 0;
   double oldn = 1;
   while (it < it_max && err > tolerance) {
        it++;
        v = 0.; a = 0.;
        for (ms.TopElement(); (el=ms.GetElement());) {
           Elas2DT3 eq(el);
           eq.Deviator();
           eq.Dilatation();
           eq.BodyRHS(body_f);
           a.Assembly(el,eq.A());
           v.Assembly(el,eq.b());
        }
        for (ms.TopSide(); (sd=ms.GetSide());) {
           Elas2DT3 eq(sd,u);
           eq.BoundaryRHS(bound_f);
           eq.SignoriniContact(bound_f,penal);
           a.Assembly(sd,eq.A());
           v.Assembly(sd,eq.b());
        }

        a.Prescribe(ms,v,bc);
        a.Factor();
        a.Solve(v);

        err = 0.;
        for (int i=1; i<=nn; i++)
           err += (v(i)-u(i))*(v(i)-u(i));
        err = sqrt(err)/oldn;
        if (verbose > 0)
          cout << "Iteration : " << it << ". Discrepancy : " << err << endl;
        u = v;
        oldn = u.Norm2();
   }

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

   cout << "\nResults after : " << it << " iterations.\n";
   cout << "Error = " << err << endl << endl;
   if (save_flag) {
     FDF pl_file(data.AuxFile(1),"w");
     pl_file.Put(uf);
   }

   return 0;
}