streamf.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 "NodeVect.h"
#include "Vect.h"

using OFELI::Mesh;
using OFELI::Element;
using OFELI::Vect;
using OFELI::NodeVect;

void StreamF(const Mesh &ms, const Vect<> &u, NodeVect<> &psi)
//-----------------------------------------------------------------------------
//                               Calculate Stream Function
//-----------------------------------------------------------------------------
{
   double eps = 1.e-5;
   const Element *e;
   static double x[3], y[3], xx[3], yy[3];
   double a, b;
   int n0=0, n1=0, n2=0, n3=0;
   Vect<int> code(ms.NbNodes());
   int nb = 0;
   psi = 0;

   double vv = 0;
   for (ms.TopElement(); (e=ms.GetElement()) && (vv < eps);) {
      n1 = e->NodeLabel(1);
      n2 = e->NodeLabel(2);
      n3 = e->NodeLabel(3);
      vv  = u(2*n1-1)*u(2*n1-1) + u(2*n1)*u(2*n1) + u(2*n2-1)*u(2*n2-1) + 
            u(2*n2)*u(2*n2) + u(2*n3-1)*u(2*n3-1) + u(2*n3)*u(2*n3);
   }
   if (vv < eps) {
     cerr << "Error in StreamF() : No Flow !!\n";
     exit(1);
   }

   vv = 0;
   for (size_t i=1; (i<=3 && vv < eps); i++) {
     n0 = e->NodeLabel(i);
     vv = u(2*n0-1)*u(2*n0-1) + u(2*n0)*u(2*n0);
   }
   code(n0) = 1;

   while (nb<ms.NbElements()) {
      for (ms.TopElement(); (e=ms.GetElement());) {
         for (size_t j=0; j<3; j++) {
            x[j] = e->PtrNode(j+1)->Coord(1);
            y[j] = e->PtrNode(j+1)->Coord(2);
         }
         xx[0] = x[2]-x[1]; xx[1] = x[0]-x[2]; xx[2] = x[1]-x[0];
         yy[0] = y[2]-y[1]; yy[1] = y[0]-y[2]; yy[2] = y[1]-y[0];
         double det = yy[2]*xx[1] - xx[2]*yy[1];
         size_t i1 = e->NodeLabel(1);
         size_t i2 = e->NodeLabel(2);
         size_t i3 = e->NodeLabel(3);
         double uu = OFELI_SIXTH*(u(2*i1-1)+u(2*i2-1)+u(2*i3-1))*det;
         double vv = OFELI_SIXTH*(u(2*i1)+u(2*i2)+u(2*i3))*det;
         int c = code(i1)*100 + code(i2)*10 + code(i3);
         if (c != 0) {
           switch (c) {
             case   1 : n0 = e->NodeLabel(3);
                        n1 = e->NodeLabel(1);
                        n2 = e->NodeLabel(2);
                        a = (2*uu - xx[2]*psi(n0,1))/det;
                        b = (2*vv - yy[2]*psi(n0,1))/det;
                        psi(n1,1) =  yy[1]*a - xx[1]*b;
                        psi(n2,1) = -yy[0]*a + xx[0]*b;
                        break;
             case  10 : n0 = e->NodeLabel(2);
                        n1 = e->NodeLabel(3);
                        n2 = e->NodeLabel(1);
                        a = (2*uu - xx[1]*psi(n0,1))/det;
                        b = (2*vv - yy[1]*psi(n0,1))/det;
                        psi(n1,1) =  yy[0]*a - xx[0]*b;
                        psi(n2,1) = -yy[2]*a + xx[2]*b;                        
                        break;
             case  11 : n0 = e->NodeLabel(1);
                        n1 = e->NodeLabel(3);
                        a = (2*uu - xx[2]*psi(n1,1))/det;
                        b = (2*vv - yy[2]*psi(n1,1))/det;
                        psi(n0,1) = yy[1]*a - xx[1]*b;
                        break;
             case 100 : n0 = e->NodeLabel(1);
                        n1 = e->NodeLabel(2);
                        n2 = e->NodeLabel(3);
                        a = (2*uu - xx[0]*psi(n0,1))/det;
                        b = (2*vv - yy[0]*psi(n0,1))/det;
                        psi(n1,1) =  yy[2]*a - xx[2]*b;
                        psi(n2,1) = -yy[1]*a + xx[1]*b;                        
                        break;
             case 101 : n0 = e->NodeLabel(2);
                        n1 = e->NodeLabel(1);
                        a = (2*uu - xx[0]*psi(n1,1))/det;
                        b = (2*vv - yy[0]*psi(n1,1))/det;
                        psi(n0,1) = yy[2]*a - xx[2]*b;
                        break;
             case 110 : n0 = e->NodeLabel(3);
                        n1 = e->NodeLabel(2);
                        a = (2*uu - xx[1]*psi(n1,1))/det;
                        b = (2*vv - yy[1]*psi(n1,1))/det;
                        psi(n0,1) = yy[0]*a - xx[0]*b;
                        break;
             case 111 : break;
           }
           nb++;
           code(i1) = code(i2) = code(i3) = 1;
         }
      }
   }
}