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

                                    O  F  E  L  I

                            Object  Finite  Element  Library

  ==============================================================================

   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

  ==============================================================================

             An example of a Finite Element Code using OFELI

            Solution of a 2-D Transient Heat Equation using :
                  - P1 Finite elements for space discretization
                  - Implicit Euler scheme for time discretization

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


#include "OFELI.h"
#include "Therm.h"
#include "User.h"
using namespace OFELI;

int main(int argc, char *argv[])
{
   Element *el;
   Side    *sd;

// Read and output mesh data
   if (argc <= 1) {
      cout << " Usage : ex4 <parameter_file>" << endl;
      exit(1);
   }
   IPF data(argv[1]);
   double max_time = data.MaxTime();
   double deltat = data.TimeStep();
   Mesh ms(data.MeshFile(1));
   User ud(ms);

// Declare problem data (matrix, rhs, boundary conditions, body forces)
   SkSMatrix<double> a(ms);
   Vect<double> b(ms.NbDOF()), u(ms.NbDOF()), bc(ms.NbDOF());

// Read in initial solution
   ud.SetInitialData(u);

   int nb_step = int(max_time/deltat);
   double time = 0;

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

   for (int step=1; step<=nb_step; step++) {
      time += deltat;
      b = 0;

//    Read in boundary conditions
      ud.Time(time);
      ud.SetDBC(bc);

      for (ms.TopElement(); (el=ms.GetElement());) {

//       Declare an instance of class DC2DT3
         DC2DT3 eq(el,u,time);

//       Capacity contribution to matrix and to RHS
         eq.LCapacityToLHS(1./deltat);
         eq.LCapacityToRHS(1./deltat);

//       Diffusion contribution to matrix
         eq.Diffusion();

//       Assemble element matrix and RHS
         if (step==1)
           a.Assembly(el,eq.A());
         b.Assembly(el,eq.b());
      }

//    Loop over edges
      for (ms.TopSide(); (sd=ms.GetSide());) {
         DC2DT3 eq(sd,u,time);
         eq.BoundaryRHS(ud);
         b.Assembly(sd,eq.b());
      }

//    Impose boundary conditions on matrix and RHS
      a.Prescribe(ms,b,bc,step-1);

//    Solve the linear system of equations
      if (step == 1) 
        a.Factor();
      a.Solve(b);
      u = b;

//    Output solution
      cout << "\nSolution for time : " << time << endl << u;
   }

#ifdef WITH_PAUSE
   system("PAUSE");
#endif
   return 0;
}