pmatrix.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 "Vect.h"
#include "CalcElem.h"
#include "LocalVect.h"
#include "Assembly.h"

using namespace OFELI;

void PMatrix(double dt, const Mesh &ms, SpMatrix<> &A, Vect<> &M)
//------------------------------------------------------------------------------
//                           Calculate Pressure Matrix
//------------------------------------------------------------------------------
{
   size_t ii, jj;
   SpMatrix<> Dx(1,ms,0), Dy(1,ms,0);
   LocalVect<double,3> dx, dy;
   const Element *e;
   for (ms.TopElement(); (e=ms.GetElement());) {
      double det = CalcElem(*e,dx,dy);
      double z = OFELI_SIXTH*det/dt;
      double c = 0.35*dt/det;
      for (size_t i=1; i<=3; i++) {
         M(e->NodeLabel(i)) += z;
         for (size_t j=1; j<=3; j++) {
            ii = e->NodeLabel(i);
            jj = e->NodeLabel(j);
            A.Add(ii,jj,c*(dx(i)*dx(j)+dy(i)*dy(j)));
            Dx.Add(ii,jj,dx(j));
            Dy.Add(ii,jj,dy(j));
         }
      }
   }

   for (size_t i=1; i<=A.Size(); i++) {
      for (size_t jj=A.RowPtr(i); jj<A.RowPtr(i+1); jj++) {
         int j = A.ColInd(jj);
         for (size_t kk=Dx.RowPtr(i); kk<Dx.RowPtr(i+1); kk++) {
            int k = Dx.ColInd(kk);
            double d = OFELI_SIXTH*OFELI_SIXTH/M(k);
            for (size_t ll=Dx.RowPtr(j); ll<Dx.RowPtr(j+1); ll++) {
               if (k == Dx.ColInd(ll)) {
                 A.Add(i,j,d*(Dx(i,k)*Dx(j,k) + Dy(i,k)*Dy(j,k)));
                 break;
               }
            }
         }
      }
   }
}