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

   Copyright (C) 1998 - 2002 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 "Mesh.h"
#include "Vect.h"
#include "LocalVect.h"
#include "LocalMatrix.h"
#include "Assembly.h"
#include "CalcElem.h"
#include "Precond.h"
#include "CG.h"

#include "SpMatrix.h"

using namespace OFELI;

int CalPres(double deltat, const Mesh &ms, const SpMatrix<> &A, 
            const ILUPrec< double,SpMatrix<> > &P, const Vect<> &u, 
            const Vect<> &ub, Vect<> &p, Vect<> &q)
//------------------------------------------------------------------------------
//                           Calculate Pressure
//------------------------------------------------------------------------------
{
   real_t toler = 1.e-12;
   LocalVect<double,3> dx, dy;
   const Element *e;
   Vect<> b(p.Size());
   for (ms.TopElement(); (e=ms.GetElement());) {
      size_t n = e->Label();
      CalcElem(*e,dx,dy);
      LocalVect<double,6> ue(e,u);
      double d = OFELI_SIXTH*(dx[0]*ue[0] + dy[0]*ue[1] +
                              dx[1]*ue[2] + dy[1]*ue[3] +
                              dx[2]*ue[4] + dy[2]*ue[5]);
      b(e->NodeLabel(1)) += 0.225*(dx[0]*ub(2*n-1) + dy[0]*ub(2*n)) - d;
      b(e->NodeLabel(2)) += 0.225*(dx[1]*ub(2*n-1) + dy[1]*ub(2*n)) - d;
      b(e->NodeLabel(3)) += 0.225*(dx[2]*ub(2*n-1) + dy[2]*ub(2*n)) - d;
   }
   int nb_it = CG(A,P,b,q,1000,toler,0);
   p.MultAdd(q,2);
   return nb_it;
}