CartesianGrid.H

Go to the documentation of this file.

 
#ifndef _CARTESIANGRID_H_
#define _CARTESIANGRID_H_
 
#include <array>
#include <vector>
#include "Gas.H"
#include "../Core/Geometry.H"
#include "../IO/ParmParser.H"
#include "../Utils/ErrCodes.H"
#include "../Utils/RealTensor2.H"
#include "../Utils/Vec3.H"
 
namespace criptic {
  namespace gas {
 
    class CartesianGrid : public Gas {
 
    public:
 
      typedef enum {
    denIdx = 0,
    ionDenIdx,
    vxIdx,
    vyIdx,
    vzIdx,
    BxIdx,
    ByIdx,
    BzIdx,
    xH0Idx,
    xHpIdx,
    xHe0Idx,
    xHepIdx,
    xHep2Idx,
    xeIdx,
    ZIdx,
    nf
      } gasIdxType;
      
      CartesianGrid(const ParmParser& pp,
            const Geometry& geom_,
            const IdxType nBB_,
            const int nGhost_ = 1);
 
      CartesianGrid(const IdxVec& nGrid_,
            const Geometry& geom_,
            const IdxType nBB_,
            const int nGhost_ = 1);
      
      virtual void frame(const RealVec &x,
             const Real t,
             RealVec& v,
             TNBBasis& tnb) const override;
 
      virtual GasData gasData(const RealVec &x,
                  const Real t) const override;
      
      virtual Real dxGhost() const override { return nGhost * dx.min(); }
 
      void applyBC();
 
      const IdxVec& gridSize() const { return nGrid; }
 
      IdxType nComponents() const { return nComp; }
 
      const RealBox& xBox() const { return geom.xBox(); }
 
      const RealVec& xLo() const { return geom.xLo(); }
 
      const RealVec& xHi() const { return geom.xHi(); }
 
      const RealVec& xSize() const { return geom.xSize(); }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      IdxType idx(const Vec3<T>& iv, IdxType n) const {
    IdxType i = n + nComp *
      (iv[0] + nGhost + (nGrid[0] + 2*nGhost) *
       (iv[1] + nGhost + (nGrid[1] + 2*nGhost) *
        (iv[2] + nGhost)));
    return i;
      }
 
      void posToIdx(const RealVec& x,
            SIdxVec& iv,
            RealVec& f) const {
    RealVec ix = (x - geom.xLo()) / dx - 0.5;
    iv = floor(ix);
    f = ix - iv;
      }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      RealVec idxToPos(const Vec3<T>& iv,
               const RealVec& f = zeroIdxVec) const {
    return geom.xLo() + (iv + 0.5) * dx;
      }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      Real& operator()(const Vec3<T>& iv, const IdxType n) {
    return data[idx(iv, n)];
      }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      const Real& operator()(const Vec3<T>& iv, const IdxType n) const {
    return data[idx(iv, n)];
      }
      
      Real& operator()(const SignedIdxType i,
               const SignedIdxType j,
               const SignedIdxType k,
               const IdxType n) {
    SIdxVec iv(i,j,k);
    return data[idx(iv, n)];
      }
    
      const Real& operator()(const SignedIdxType i,
                 const SignedIdxType j,
                 const SignedIdxType k,
                 const IdxType n) const {
    SIdxVec iv(i,j,k);
    return data[idx(iv, n)];
      }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      std::vector<Real> interp(const Vec3<T>& iv,
                   const RealVec& f,
                   const IdxType sComp = 0,
                   const IdxType eComp = maxIdx) const {
    IdxType eComp_ = eComp > nComp ? nComp : eComp;
    assert(sComp < eComp_);
    std::vector<Real> q(eComp_ - sComp);
    SIdxVec ivTmp;
    for (ivTmp[2] = iv[2]; ivTmp[2] <= iv[2]+1; ivTmp[2]++) {
      Real fz = (ivTmp[2] == iv[2]) ? 1.0 - f[2] : f[2];
      for (ivTmp[1] = iv[1]; ivTmp[1] <= iv[1]+1; ivTmp[1]++) {
        Real fy = (ivTmp[1] == iv[1]) ? 1.0 - f[1] : f[1];
        for (ivTmp[0] = iv[0]; ivTmp[0] <= iv[0]+1; ivTmp[0]++) {
          Real fx = (ivTmp[0] == iv[0]) ? 1.0 - f[0] : f[0];
          for (IdxType c = sComp; c < eComp_; c++) {
        q[c - sComp] += fx * fy * fz * (*this)(ivTmp, c);
          }
        }
      }
    }
    return q;
      }
      
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      RealVec interpVec(const Vec3<T>& iv,
            const RealVec& f,
            const IdxType c) const {
    RealVec q;
    SIdxVec ivTmp;
    for (ivTmp[2] = iv[2]; ivTmp[2] <= iv[2]+1; ivTmp[2]++) {
      Real fz = (ivTmp[2] == iv[2]) ? 1.0 - f[2] : f[2];
      for (ivTmp[1] = iv[1]; ivTmp[1] <= iv[1]+1; ivTmp[1]++) {
        Real fy = (ivTmp[1] == iv[1]) ? 1.0 - f[1] : f[1];
        for (ivTmp[0] = iv[0]; ivTmp[0] <= iv[0]+1; ivTmp[0]++) {
          Real fx = (ivTmp[0] == iv[0]) ? 1.0 - f[0] : f[0];
          for (IdxType c_ = 0; c_ < 3; c_++) {
        q[c_] += fx * fy * fz * (*this)(ivTmp, c + c_);
          }
        }
      }
    }
    return q;
      }
      
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      RealVec interpGradScalar(const Vec3<T>& iv,
                   const RealVec& f,
                   const IdxType c) const {
    RealVec q;
    for (int d = 0; d < 3; d++) {
      SIdxVec ivTmp;
      for (ivTmp[2] = iv[2]; ivTmp[2] <= iv[2]+1; ivTmp[2]++) {
        Real fz;
        if (d != 2) {
          fz = (ivTmp[2] == iv[2]) ? 1.0 - f[2] : f[2];
        } else {
          fz = (ivTmp[2] == iv[2]) ? -1.0 : 1.0;
        }
        for (ivTmp[1] = iv[1]; ivTmp[1] <= iv[1]+1; ivTmp[1]++) {
          Real fy;
          if (d != 1) {
        fy = (ivTmp[1] == iv[1]) ? 1.0 - f[1] : f[1];
          } else {
        fy = (ivTmp[1] == iv[1]) ? -1.0 : 1.0;
          }
          for (ivTmp[0] = iv[0]; ivTmp[0] <= iv[0]+1; ivTmp[0]++) {
        Real fx;
        if (d != 0) {
          fx = (ivTmp[0] == iv[0]) ? 1.0 - f[0] : f[0];
        } else {
          fx = (ivTmp[0] == iv[0]) ? -1.0 : 1.0;
        }
        q[d] += fx * fy * fz * (*this)(ivTmp, c);
          }
        }
      }
    }
    return q/dx;
      }
 
      template<class T,
           std::enable_if_t<std::is_integral<T>::value, bool> = true >
      RealTensor2 interpGradVector(const Vec3<T>& iv,
                   const RealVec& f,
                   const IdxType c) const {
    RealTensor2 q;
    for (int d = 0; d < 3; d++) {
      SIdxVec ivTmp;
      for (ivTmp[2] = iv[2]; ivTmp[2] <= iv[2]+1; ivTmp[2]++) {
        Real fz;
        if (d != 2) {
          fz = (ivTmp[2] == iv[2]) ? 1.0 - f[2] : f[2];
        } else {
          fz = (ivTmp[2] == iv[2]) ? -1.0 : 1.0;
        }
        for (ivTmp[1] = iv[1]; ivTmp[1] <= iv[1]+1; ivTmp[1]++) {
          Real fy;
          if (d != 1) {
        fy = (ivTmp[1] == iv[1]) ? 1.0 - f[1] : f[1];
          } else {
        fy = (ivTmp[1] == iv[1]) ? -1.0 : 1.0;
          }
          for (ivTmp[0] = iv[0]; ivTmp[0] <= iv[0]+1; ivTmp[0]++) {
        Real fx;
        if (d != 0) {
          fx = (ivTmp[0] == iv[0]) ? 1.0 - f[0] : f[0];
        } else {
          fx = (ivTmp[0] == iv[0]) ? -1.0 : 1.0;
        }
        for (int c_ = 0; c_ < 3; c_++) {
          q[3*c_ + d] += fx * fy * fz * (*this)(ivTmp, c + c_);
        }
          }
        }
      }
    }
    for (int c_ = 0; c_ < 3; c_++) {
      for (int d = 0; d < 3; d++) {
        q[3*c_ + d] /= dx[d];
      }
    }
    return q;
      }
 
    private:
 
      // Private data
      const Geometry& geom;   
      IdxVec nGrid;           
      IdxType nComp;          
      IdxType nBB;            
      IdxType nGhost;         
      RealVec dx;             
      std::vector<Real> data; 
    };    
  }
}
 
#endif
// _CARTESIANGRID_H_