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

                              **********************
                              *    T  I  F  F 2    *
                              **********************

                 A Finite Element Code for Transient Incompressible
                       Fluid Flow Simulations 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 "Fluid.h"
#include "User.h"
using namespace OFELI;

int main(int argc, char *argv[])
{
   Mesh     ms;
   Element  *el;
   Side     *sd;
   int      i, step;
   ifstream mf, bcf, bodyf, boundf, inf;
   FDF      *p_file=NULL, *v_file=NULL;
   double   time = 0;

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

   IPF data("tiff2 - 1.0",argv[1]);
   int output_flag = data.Output();
   int save_flag = data.Save();
   int init_flag = data.Init();
   int bc_flag = data.BC();
   int body_flag = data.BF();
   int bound_flag = data.SF();
   int pres_flag = data.IntPar(1);
   double deltat = data.TimeStep();
   double max_time = data.MaxTime();

   if (output_flag) {
     cout << "=====================================================================\n\n";
     cout << "                             T  I  F  F  2\n\n";
     cout << "         A Finite Element Code for Transient Incompressible\n";
     cout << "                 Fluid Flow Simulation in 2-D Geometries\n\n\n";
     cout << "              tiff2 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));
   User ud(ms);
   int nb_dof = ms.Dim();

// Print Mesh data
   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";
   SkSMatrix<double> a(ms);
   Vect<double> b(ms.NbDOF());

// Read initial condition, boundary conditions, body and boundary forces
   Vect<double> u(ms.NbDOF());
   if (output_flag > 1)
     cout << "Reading initial condition ...\n";
   if (!init_flag)
     ud.SetInitialData(u);
   else {
     FDF in_file(data.InitFile(),"r");
     NodeVect<double> ui(ms,nb_dof);
     u = in_file.Get(ui);
   }

   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,nb_dof);
     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,nb_dof);
     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 sf_file(data.SFFile(),"r");
     NodeVect<double> ui(ms,nb_dof);
     sf_file.Get(ui);
     bound_f = Vect<double>(ui);
   }

   int transient = 0;
   int nb_step = 1;
   if (deltat <= max_time) {
     transient = 1;
     nb_step = (int)(max_time/deltat);
   }
   NodeVect<double> uf(ms,nb_dof);
   if (save_flag) {
     v_file = new FDF(data.AuxFile(1),"w");
     p_file = new FDF(data.AuxFile(2),"w");
   }

// Loop over time steps
// --------------------

   for (step=1; step<=nb_step; step++) {

      if (output_flag > 1 && transient)
        cout << "Performing time step " << step << " ..." << endl;
      time += deltat;
      cout << "Step : " << step << ", Time : " << time << endl;
      b = 0;

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

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

      for (ms.TopElement(); (el=ms.GetElement());) {
         NSP2DQ41 eq(el,u,time);
         eq.LMass(1./deltat);
         eq.Penal(1.e07);
         eq.Viscous(0.001);
         eq.RHS_Convection();
         eq.BodyRHS(ud);
         if (step==1) 
           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());) {
         NSP2DQ41 eq(sd);
         eq.BoundaryRHS(ud);
         b.Assembly(sd,eq.b());
      }

//    Impose Boundary Conditions and Solve the linear system
//    ------------------------------------------------------

      a.Prescribe(ms,b,bc,step-1);
      if (step == 1)
        a.Factor();
      a.Solve(b);
      u = b;

//    Output and/or Store solution
//    ----------------------------
 
      uf.FromVect(u,1,"Velocity",time);
      if (output_flag > 0)
        cout << uf;
      if (save_flag && step==nb_step)
        v_file->Put(uf);

//    Calculate and smooth pressure
//    -----------------------------

      if (pres_flag) {
        double pres = 0.;
        Vect<double> pm(ms.NbNodes());
        NodeVect<double> p(ms,1,"Pressure",time);

        for (ms.TopElement(); (el=ms.GetElement());) {
           NSP2DQ41 eq(el,u,time);
           pres += eq.Pressure(1.e07);
           eq.PresMat();
           for (i=1; i<=eq.NbNodes(); i++) {
              int n = el->PtrNode(i)->Label();
              pm(n) += eq.PM()(i);
              p(n,1) += pres;
           }
        }
        for (int n=1; n<=ms.NbNodes(); n++)
           p(n,1) /= pm(n);
        if (output_flag > 0)
          cout << p;
        if (save_flag  && step==nb_step)
          p_file->Put(p);
      }
   }

   if (save_flag > 0) {
     delete p_file;
     delete v_file;
   }
   return 0;
}