vmatrix.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 <stdio.h>

#include "Mesh.h"
#include "SpMatrix.h"
#include "Vect.h"
#include "CalcElem.h"
#include "LocalMatrix.h"
#include "LocalVect.h"
#include "Line2.h"
#include "Assembly.h"

using namespace OFELI;

void VMatrix(double deltat, double dens, double visc, const Mesh &ms,
             SpMatrix<> &A)
{
   void BCMat(const Element &e, LocalMatrix<double,6,6> &a);
   const Element *e;
   const Side *s;
   double aa, ax, ay;
   LocalVect<double,3> dx, dy;
   LocalMatrix<double,6,6> a;
   const double penal = 1.e20;

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

   for (ms.TopElement(); (e=ms.GetElement());) {
      double det = CalcElem(*e,dx,dy);
      double c = OFELI_SIXTH*dens*det/deltat;
      double aa = 0.25*visc/det;
      for (size_t i=0; i<3; i++) {
         ax = aa*dx[i]; ay = aa*dy[i];
         for (size_t j=i; j<3; j++) {
            a(2*i+1,2*j+1) = 2*ax*dx[j] + ay*dy[j];
            a(2*i+1,2*j+2) = ay*dx[j];
            a(2*i+2,2*j+1) = ax*dy[j];
            a(2*i+2,2*j+2) = 2*ay*dy[j] + ax*dx[j];
         }
      }
      a(1,1) += c; a(2,2) += c; a(3,3) += c;
      a(4,4) += c; a(5,5) += c; a(6,6) += c;
      a.Symmetrize();

//    Boundary conditions
      BCMat(*e,a);

//    Assembly
      Assembly<double,6,Element>(e,a,A);
   }

// Periodic Boundary condition
   for (ms.TopSide(); (s=ms.GetSide());) {
      Line2 ln(s);
      for (size_t i=1; i<=2; i++) {
         double c = 0.5*ln.Length()*penal;
         if (s->Code(1) == PERIODIC_A)
           a(2*i-1,2*i-1) += c;
         else if (s->Code(1) == PERIODIC_B)
           a(2*i-1,2*i-1) -= c;
         if (s->Code(2) == PERIODIC_A)
           a(2*i  ,2*i  ) += c;
         else if (s->Code(2) == PERIODIC_B)
           a(2*i  ,2*i  ) -= c;
      }
      Assembly<double,6,Side>(s,a,A);
   }
}


void BCMat(const Element &e, LocalMatrix<double,6,6> &a)
{
   for (size_t i=1; i<=3; i++) {
      for (size_t j=1; j<=3; j++)
         if (i!=j) {
           if (e.PtrNode(i)->Code(1))
             a(2*i-1,2*j-1) = 0;
           if (e.PtrNode(i)->Code(2))
             a(2*i  ,2*j  ) = 0;
         }
   }
}