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

                             *************************
                             *     L E L A S 2 D     *
                             *************************

                A Finite Element Code for Linearized Elastostatics in
                                  Two Dimensions

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

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

int main(int argc, char *argv[])
{
   Mesh     ms;
   Element  *el;
   Side     *sd;
   ifstream mf, bcf, bodyf, boundf;

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

   IPF data("lelas2d - 1.0",argv[1]);
   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 2 D                     *\n";
     cout << "    *     Two-Dimensional Linearized Elastostatics        *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "         A Finite Element Code for Linearized Elastostatics\n";
     cout << "                             in 2-D Geometries\n\n";
     cout << "            LELAS2D 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 (output_flag > 1)
     cout << "Reading mesh data ...\n";
   ms.Get(data.MeshFile(1));
   int nb_dof = ms.Dim();
   if (output_flag > 1)
     cout << ms;

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

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

   Vect<double> body_f(ms.NbDOF());
   if (body_flag) {
     if (output_flag > 1)
       cout << "Reading Body Forces ...\n";
     FDF bf_file(data.BFFile(),"r");
     NodeVect<double> ui(ms,ms.NbDOF());
     bf_file.Get(ui);
     body_f = Vect<double>(ui);
   }

   if (output_flag > 1)
     cout << "Reading Boundary Tractions ...\n";
   Vect<double> bound_f(ms.NbDOF());
   if (bound_flag) {
     FDF bn_file(data.SFFile(),"r");
     NodeVect<double> ui(ms,ms.NbDOF());
     bn_file.Get(ui);
     bound_f = Vect<double>(ui);
   }

   NodeVect<double> uf(ms,nb_dof);

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

   if (output_flag > 1)
     cout << "Looping over elements ...\n";
   for (ms.TopElement(); (el=ms.GetElement());) {
      Elas2DQ4 eq(el);
      eq.Deviator();
      eq.Dilatation();
      if (body_flag)
        eq.BodyRHS(Vect<double>(el,body_f));
      else
        eq.BodyRHS(ud);
      a.Assembly(el,eq.A());
      b.Assembly(el,eq.b());
   }

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

   if (output_flag > 1)
     cout << "Looping over sides ...\n";
   for (ms.TopSide(); (sd=ms.GetSide());) {
      Elas2DQ4 eq(sd);
      if (bound_flag)
        eq.BoundaryRHS(Vect<double>(sd,bound_f));
      else
        eq.BoundaryRHS(ud);
      b.Assembly(sd,eq.b());
   }

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

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

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

   if (save_flag) {
     FDF pl_file(data.PlotFile(),FDF_WRITE);
     pl_file.Put(uf);
   }

// Calculate principal and Von-Mises stresses
// ------------------------------------------

   if (output_flag > 1)
     cout << "Calculating stresses ...\n";
   ElementVect<double> st(ms, 3, "Principal Stress");
   ElementVect<double> vm(ms, 1, "Von-Mises Stress");
   for (ms.TopElement(); (el=ms.GetElement());) {
      Vect<double> ste(3);
      Elas2DQ4 eq(el,b);
      eq.Stress(ste,vm(el->Label(),1));
      for (int k=1; k<=3; k++)
         st(el->Label(),k) = ste(k);
   }
   return 0;
}