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

   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 "Mesh.h"
#include "LocalVect.h"
#include "CalcElem.h"

using namespace OFELI;

void UpdateV(const Mesh &ms, double dt, const Vect<> &M, const Vect<> &q,
             Vect<> &u, Vect<> &ub)
//------------------------------------------------------------------------------
//                             Update velocity field
//------------------------------------------------------------------------------
{
   const Element *e;
   LocalVect<double,3> dx, dy, b;
   static double z = 126*dt/81.0;
   for (ms.TopElement(); (e=ms.GetElement());) {
      double det = CalcElem(*e,dx,dy);
      LocalVect<double,3> qe(e,q,1);
      ub(2*e->Label()-1) -= z*Dot(dx,qe)/det;
      ub(2*e->Label()  ) -= z*Dot(dy,qe)/det;
      double d = OFELI_SIXTH*(qe[0] + qe[1] + qe[2]);
      for (size_t i=0; i<3; i++) {
         size_t n = e->NodeLabel(i+1);
         u(2*n-1) += d*dx[i]/M(n);
         u(2*n  ) += d*dy[i]/M(n);
      }
   }
}