TwoCEulerEvaluator.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_TWOCOMP_EULER_EVALUATOR_H
#define AMROC_TWOCOMP_EULER_EVALUATOR_H

#include <iostream>
#include <cstdio>
#include <cfloat>

#include "Evaluator.h"

template <class VectorType>
class TwoCompEulerEvaluator : public Evaluator<VectorType> {
  typedef Evaluator<VectorType> base;
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef typename base::vector_block_type vector_block_type;
  typedef typename base::data_block_type data_block_type;

  TwoCompEulerEvaluator() : evaluator_base_type() {
    std::sprintf(title,"3D-Euler equations - 2 Components");
  }

  virtual void SetUp(VisualGridBase *gr) {
    evaluator_base_type::SetUp(gr);
    const int base=evaluator_base_type::NKeys();
    std::sprintf(&(tkeys[   base *LEN_TKEYS]), "Pressure                ");
    std::sprintf(&(tkeys[(base+1)*LEN_TKEYS]), "Gamma                   ");
    std::sprintf(&(tkeys[(base+2)*LEN_TKEYS]), "Schl.-Plot Pressure     ");
    std::sprintf(&(tkeys[(base+3)*LEN_TKEYS]), "Schl.-Plot Gamma        ");
    base::fkeys[base  ] = base::Grid->FKeys();
    base::fkeys[base+1] = base::Grid->FKeys();
    base::fkeys[base+2] = base::Grid->FKeys();
    base::fkeys[base+3] = base::Grid->FKeys();
    base::ikeys[base  ] = 112; 
    base::ikeys[base+1] = 103; 
    base::ikeys[base+2] = 80; 
    base::ikeys[base+3] = 71; 
  }

  virtual int NKeys() const { return evaluator_base_type::NKeys()+4; }

  virtual void SetScalCells(int key,float v[],int& offset, 
                            vector_block_type& DB, float* dx, bool* CompRead) {
    int idx=0;
    const int base=evaluator_base_type::NKeys();
    BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
    if (key == base+1) {          
      if (!CompRead[4] || !CompRead[0] || !CompRead[5])
        return;

      float rho, u1, u2, u3, E, g, pinf;
      for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
        rho = FastIndex3(dat,i,j,k)(0);               
        u1  = FastIndex3(dat,i,j,k)(1);
        u2  = FastIndex3(dat,i,j,k)(2);
        u3  = FastIndex3(dat,i,j,k)(3);
        E   = FastIndex3(dat,i,j,k)(4); 
        g   = FastIndex3(dat,i,j,k)(5);
        pinf = FastIndex3(dat,i,j,k)(6);
        v[idx+offset]=(E-0.5*rho*(u1*u1+u2*u2+u3*u3)-pinf)/g;
        mvals[2*base]   = Min(mvals[2*base],   v[idx+offset]) ;
        mvals[2*base+1] = Max(mvals[2*base+1], v[idx+offset]) ;
        idx++;        
      end_for
    }
    else if (key == base+2) {
      if (!CompRead[5])
        return;

      float g;
      for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
        g = FastIndex3(dat,i,j,k)(5);
        v[idx+offset]=1.0/g+1.0;
        mvals[2*(base+1)]   = Min(mvals[2*(base+1)],   v[idx+offset]) ;
        mvals[2*(base+1)+1] = Max(mvals[2*(base+1)+1], v[idx+offset]) ;
        idx++;        
      end_for
    }
    else
      evaluator_base_type::SetScalCells(key,v,offset,DB,dx,CompRead);

    EndFastIndex3(dat); 
    offset+=idx;
  }
    


  virtual void SetScalNodes(int key,float v[],int& offset,
                            vector_block_type& DB, const BBox& bbox, 
                            float* dx, bool* CompRead) {
    int idx=0;
    const int base=evaluator_base_type::NKeys();
    BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
    const Coords& step = bbox.stepsize();
    if (key == base+1) {          
      if (!CompRead[4] || !CompRead[0] || !CompRead[5])
        return;

      float rho, u1, u2, u3, E, c, g, pinf;
      for_3 (i, j, k, bbox, step)
        BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                     i,j,k,step(0),step(1),step(2));
        AddOver *= DB.bbox();
          
        c=0;
        v[idx+offset] = 0.0;
        for_3 (l, m, n, AddOver, step)
          rho = FastIndex3(dat,l,m,n)(0);                     
          u1  = FastIndex3(dat,l,m,n)(1);
          u2  = FastIndex3(dat,l,m,n)(2);
          u3  = FastIndex3(dat,l,m,n)(3);
          E   = FastIndex3(dat,l,m,n)(4);             
          g   = FastIndex3(dat,l,m,n)(5);
          pinf = FastIndex3(dat,l,m,n)(6);
          if (rho<FLT_MAX && u1<FLT_MAX && u2<FLT_MAX && u3<FLT_MAX && 
              E<FLT_MAX && g<FLT_MAX && pinf<FLT_MAX) {
            v[idx+offset] += (E-0.5*rho*(u1*u1+u2*u2+u3*u3)-pinf)/g;
            c += 1.0;
          }
        end_for
        v[idx+offset] /= c;

        mvals[2*base]   = Min(mvals[2*base],   v[idx+offset]) ;
        mvals[2*base+1] = Max(mvals[2*base+1], v[idx+offset]) ;
        idx++;        
      end_for 
    } 
    else if (key == base+3) {
      if (!CompRead[4] || !CompRead[0] || !CompRead[5])
        return;

      float rho, u1, u2, u3, E, c, g, pinf;
      for_3 (i, j, k, bbox, step)
        BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                     i,j,k,step(0),step(1),step(2));
        AddOver *= DB.bbox();
          
        data_block_type DBHelp(AddOver);
        DBhelp = FLT_MAX;
        BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
        for_3 (l, m, n, AddOver, step)
          rho = FastIndex3(dat,l,m,n)(0);                     
          u1  = FastIndex3(dat,l,m,n)(1);
          u2  = FastIndex3(dat,l,m,n)(2);
          u3  = FastIndex3(dat,l,m,n)(3);
          E   = FastIndex3(dat,l,m,n)(4);             
          g   = FastIndex3(dat,l,m,n)(5);
          pinf = FastIndex3(dat,l,m,n)(6);
          if (rho<FLT_MAX && u1<FLT_MAX && u2<FLT_MAX && u3<FLT_MAX && 
              E<FLT_MAX && g<FLT_MAX && pinf<FLT_MAX) 
            FastIndex3(help,l,m,n) = (E-0.5*rho*(u1*u1+u2*u2+u3*u3)-pinf)/g;
        end_for
        EndFastIndex3(help); 
        v[idx+offset] = SetScalGradNodes(DBHelp,dx);
        
        mvals[2*(base+2)]   = Min(mvals[2*(base+2)],   v[idx+offset]) ;
        mvals[2*(base+2)+1] = Max(mvals[2*(base+2)+1], v[idx+offset]) ;
        idx++;        
      end_for 
    } 
    else if (key == base+2) {
      if (!CompRead[5])
        return;

      float g, c;
      for_3 (i, j, k, bbox, step)
        BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                     i,j,k,step(0),step(1),step(2));
        AddOver *= DB.bbox();
          
        c=0;
        v[idx+offset] = 0.0;
        for_3 (l, m, n, AddOver, step)
          g = FastIndex3(dat,l,m,n)(5);
          if (g<FLT_MAX) {
            v[idx+offset] += 1.0/g+1.0;
            c += 1.0;
          }
        end_for
        v[idx+offset] /= c;

        mvals[2*(base+1)]   = Min(mvals[2*(base+1)],   v[idx+offset]) ;
        mvals[2*(base+1)+1] = Max(mvals[2*(base+1)+1], v[idx+offset]) ;
        idx++;        
      end_for
    }
    else if (key == base+4) {
      if (!CompRead[5])
        return;

      float g, c;
      for_3 (i, j, k, bbox, step)
        BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                     i,j,k,step(0),step(1),step(2));
        AddOver *= DB.bbox();
          
        data_block_type DBHelp(AddOver);
        DBhelp = FLT_MAX;
        BeginFastIndex3(help, DBHelp.bbox(), DBHelp.data(), DataType);
        for_3 (l, m, n, AddOver, step)
          g = FastIndex3(dat,l,m,n)(5);
          if (g<FLT_MAX) 
            FastIndex3(help,l,m,n) = 1.0/g+1.0;
        end_for
        EndFastIndex3(help); 
        v[idx+offset] = SetScalGradNodes(DBHelp,dx);

        mvals[2*(base+3)]   = Min(mvals[2*(base+3)],   v[idx+offset]) ;
        mvals[2*(base+3)+1] = Max(mvals[2*(base+3)+1], v[idx+offset]) ;
        idx++;        
      end_for
    }
    else
      evaluator_base_type::SetScalNodes(key,v,offset,DB,bbox,dx,CompRead);

    EndFastIndex3(dat); 
    offset+=idx;
  }

};

#endif

---