momentum.cpp

OFELI is an object oriented library of C++ classes for development of finite element codes. Its main

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

   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 "SpMatrix.h"
#include "Precond.h"
#include "CG.h"
#include "LocalMatrix.h"
#include "Vect.h"
#include "LocalVect.h"
#include "UserData.h"
#include "CalcElem.h"
#include "VDF.h"
#include "Line2.h"
#include "util.h"
#include "Assembly.h"

using namespace OFELI;


void Momentum(int step, double time, double deltat, Mesh &ms, VDF &vdf,
              Vect<> &u, Vect<> &v, Vect<> &ub, Vect<> &c, const SpMatrix<> &A,
              const ILUPrec< double,SpMatrix<> > &P, const Vect<> &p, double dens,
              double visc, double &cfl, int verbose)
//------------------------------------------------------------------------------
//                    Construct and Solve Momentum Equations
//------------------------------------------------------------------------------
{
   double toler = 1.e-12;
   Element *e;
   size_t i, j;
   double cc, ax, ay, dxbp, dybp, bm, bv, dudx, dudy, dvdx, dvdy;
   LocalVect<double,3> dx, dy;
   LocalVect<double,6> b, ce, us;
   Vect<> f(2*ms.NbNodes()), ubc(2*ms.NbNodes()), sf(2*ms.NbNodes());

   Scale(0.5,c,v);
   c = 0;
   vdf.Get(BODY_FORCE,f,time);
   vdf.Get(TRACTION,sf,time);

//===================================================================
//                         LOOP OVER ELEMENTS
//===================================================================

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

      size_t n = e->Label();
      LocalVect<double,3> pe(e,p,1);
      LocalVect<double,6> ue(e,u);
      us[0] = 0.5*(ue[2] + ue[4]); us[1] = 0.5*(ue[3] + ue[5]);
      us[2] = 0.5*(ue[4] + ue[0]); us[3] = 0.5*(ue[5] + ue[1]);
      us[4] = 0.5*(ue[0] + ue[2]); us[5] = 0.5*(ue[1] + ue[3]);
      double det = CalcElem(*e,dx,dy);

//-------------------------------------------------------------------
//    cfl in element
//-------------------------------------------------------------------

      double d1 = ue[0] + ue[2] + ue[4];
      double d2 = ue[1] + ue[3] + ue[5];
      cfl = max(cfl,deltat*sqrt(2*(d1*d1+d2*d2)/(9*det)));

//-------------------------------------------------------------------
//    Mass
//-------------------------------------------------------------------

      cc = OFELI_SIXTH*dens*det/deltat;
      Scale(cc,ue,b);
      cc *= 0.5;
      b[0] = cc*(ue[0]+us[0]); b[1] = cc*(ue[1]+us[1]);
      b[2] = cc*(ue[2]+us[2]); b[3] = cc*(ue[3]+us[3]);
      b[4] = cc*(ue[4]+us[4]); b[5] = cc*(ue[5]+us[5]);

//-------------------------------------------------------------------
//    Viscous Term (R.H.S.)
//-------------------------------------------------------------------

      cc = 0.25*visc/det;
      for (i=0; i<3; i++) {
         ax = cc*dx[i]; ay = cc*dy[i];
         for (j=0; j<3; j++) {
            b[2*i  ] -= (2*ax*dx[j]+ay*dy[j])*ue[2*j] + ay*dx[j]*ue[2*j+1];
            b[2*i+1] -= ax*dy[j]*ue[2*j] + (2*ay*dy[j]+ax*dx[j])*ue[2*j+1];
         }
      }

//-------------------------------------------------------------------
//    Convection
//-------------------------------------------------------------------

      cc = 0.5*OFELI_TWELVETH*dens;
      dudx = dx[0]*ue[0] + dx[1]*ue[2] + dx[3]*ue[4];
      dudy = dy[0]*ue[0] + dy[1]*ue[2] + dy[3]*ue[4];
      dvdx = dx[0]*ue[1] + dx[1]*ue[3] + dx[3]*ue[5];
      dvdy = dy[0]*ue[1] + dy[1]*ue[3] + dy[3]*ue[5];
      for (i=0; i<3; i++) {
         ce[2*i  ] = cc*((d1 + ue[2*i])*dudx + (d2 + ue[2*i+1])*dudy);
         ce[2*i+1] = cc*((d1 + ue[2*i])*dvdx + (d2 + ue[2*i+1])*dvdy);
      }
      Axpy(-1.5,ce,b);
      Assembly(e,ce,c);

//-------------------------------------------------------------------
//    Pressure Gradient
//-------------------------------------------------------------------

      cc = OFELI_SIXTH*(pe(1) + pe(2) + pe(3));
      for (i=0; i<3; i++) {
         b[2*i  ] += cc*dx[i];
         b[2*i+1] += cc*dy[i];
      }

//-------------------------------------------------------------------
//    Body Force
//-------------------------------------------------------------------

      for (i=1; i<=3; i++) {
         b(2*i-1) += f(2*e->NodeLabel(i)-1)*det*OFELI_SIXTH;
         b(2*i  ) += f(2*e->NodeLabel(i)  )*det*OFELI_SIXTH;
      }

//-------------------------------------------------------------------
//    Bubble velocities
//-------------------------------------------------------------------

      dxbp = -0.225*Dot(dx,pe); 
      dybp = -0.225*Dot(dy,pe);
      bm = 0.144642857*det*dens/deltat;
      bv = 1.0125*visc/det*(Dot(dx,dx) + Dot(dy,dy));
      ub(2*n-1) = (dxbp + ub(2*n-1)*(bm-bv))/(bm+bv);
      ub(2*n  ) = (dybp + ub(2*n  )*(bm-bv))/(bm+bv);

//-------------------------------------------------------------------
//    Assembly of R.H.S.
//-------------------------------------------------------------------

      Assembly(e,b,v);
   }

//-------------------------------------------------------------------
// Tractions
//-------------------------------------------------------------------

   Side *s;
   LocalVect<double,4> ss;
   for (ms.TopSide(); (s=ms.GetSide());) {
      Line2 ln(s);
      Point<> x[2];
      x[0] = s->PtrNode(1)->Coord();
      x[1] = s->PtrNode(2)->Coord();
      cc = 0.5*ln.Length();
      for (i=0; i<2; i++) {
         ss[2*i  ] = cc*sf(2*s->NodeLabel(i+1)-1);
         ss[2*i+1] = cc*sf(2*s->NodeLabel(i+1)  );
      }
      Assembly(s,ss,v);
   }

//-------------------------------------------------------------------
// Boundary Conditions
//-------------------------------------------------------------------

   vdf.Get(BOUNDARY_CONDITION,ubc,time);
   for (i=1; i<=ms.NbNodes(); i++) {
      j = ms.NodeLabel(i);
      if (ms.PtrNode(j)->Code(1))
        v(2*j-1) = ubc(2*j-1)*A(2*j-1,2*j-1);
      if (ms.PtrNode(j)->Code(2))
        v(2*j  ) = ubc(2*j  )*A(2*j  ,2*j  );
   }

//-------------------------------------------------------------------
// Solve the linear system
//-------------------------------------------------------------------

   int it = CG(A,P,v,u,1000,toler,0);
   if (verbose > 1)
     cout << "Nb. of CG iterations for Momentum : " << it << endl;
}