GridGenerator.H

Go to the documentation of this file.

#include <vector>
#include <cmath>

namespace gridgen {

  template<int numDim,size_t numX> 
  int GenerateRegularGrid(std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double scale[] = {numX};
    double xMax = 1.0;
    double xMin = 0.0;
    double xRange = xMax - xMin;
    std::vector<double> dX(numDim,0.0);
    for(int iDim = 0;iDim < numDim;iDim++){
      dX[iDim] = xRange/(scale[iDim]-1.0);
      xyz[iDim] = ijk[iDim]*dX[iDim] + xMin;
    }
    return(0);
  };

  template<int numDim,size_t numX,size_t numY> 
  int GenerateRegularGrid(std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double scale[] = {numX,numY};
    double xMax = 1.0;
    double xMin = 0.0;
    double xRange = xMax - xMin;
    std::vector<double> dX(numDim,0.0);
    for(int iDim = 0;iDim < numDim;iDim++){
      dX[iDim] = xRange/(scale[iDim]-1.0);
      xyz[iDim] = ijk[iDim]*dX[iDim] + xMin;
    }
    return(0);
  };

  template<int numDim,size_t numX,size_t numY,size_t numZ> 
  int GenerateRegularGrid(std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double scale[] = {numX,numY,numZ};
    double xMax = 1.0;
    double xMin = 0.0;
    double xRange = xMax - xMin;
    std::vector<double> dX(numDim,0.0);
    for(int iDim = 0;iDim < numDim;iDim++){
      dX[iDim] = xRange/(scale[iDim]-1.0);
      xyz[iDim] = ijk[iDim]*dX[iDim] + xMin;
    }
    return(0);
  };

  template<int numDim,size_t numX,size_t numY,size_t numZ>
  int RectilinearGrid1(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double iVar = ((double)(ijk[0]+1)/(double(numX)));
    xyz[0] = std::pow(iVar,2.0);
    iVar = ((double)(ijk[1]+1)/(double(numY)));
    xyz[1] = 2.0*std::pow(iVar,2.0);
    iVar = ((double)(ijk[2]+1)/(double(numZ)));
    xyz[2] = 3.0*std::pow(iVar,2.0);
    return(0);
  };

  template<int numDim,size_t numX,size_t numY>
  int RectilinearGrid1(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double iVar = ((double)(ijk[0]+1)/(double(numX)));
    xyz[0] = std::pow(iVar,2.0);
    iVar = ((double)(ijk[1]+1)/(double(numY)));
    xyz[1] = 2.0*std::pow(iVar,2.0);
    return(0);
  }

  template<int numDim,size_t numX>
  int RectilinearGrid1(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double iVar = ((double)(ijk[0]+1)/(double(numX)));
    xyz[0] = std::pow(iVar,2.0);
    return(0);
  }
  
  template<size_t numX,size_t numY,size_t numZ,int nWave>
  int VGWavy(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double aX = 1.0;
    double aY = 1.0;
    double aZ = 1.0;
    double lX = 4.0;
    double lY = 4.0;
    double lZ = 4.0;

    double pi = 3.141592653589793238;
    double pix2 = 2.0*pi;
    double pixWave = pi*nWave;
    double dX = lX/(numX);
    double dY = lY/(numY);
     double dZ = lZ/(numZ);
    
    xyz[0] = -lX/2.0 + dX*(ijk[0]+aX*std::sin(pixWave*ijk[1]*dY/lY)*
                           std::sin(pixWave*ijk[2]*dZ/lZ));
    
    xyz[1] = -lY/2.0 + dY*(ijk[1]+aY*std::sin(pixWave*ijk[0]*dX/lX)*
                           std::sin(pixWave*ijk[2]*dZ/lZ));
    
    xyz[2] = -lZ/2.0 + dZ*(ijk[2]+aZ*std::sin(pixWave*ijk[0]*dX/lX)*
                           std::sin(pixWave*ijk[1]*dY/lY));
    
    return(0);
  };
  
  template<int numDim,size_t numX,size_t numY>
  int CurvilinearGrid1(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    double r0   = 0.25;
    double rmax = 1.25;
    double rt   = 2.0;
    double dr   = (rmax - r0)/(double(numX-1));
    double d0   = (2.0*M_PI)/(double(numY));
    double r    = r0 + ijk[0]*dr;
    double a0   = ijk[1]*d0;
    
    
    xyz[0] = r*std::cos(a0);
    xyz[1] = r*std::sin(a0);  

    return(0);
  };

  int Cylinder2D(const std::vector<int>    &inOptions,
                 const std::vector<double> &inParameters,
                 const std::vector<size_t> &gridSizes,
                 const std::vector<size_t> &ij,
                 std::vector<double>       &xy)
  {
    double r0   = inParameters[0];
    double rmax = inParameters[1];
    double dr   = (rmax - r0)/(double(gridSizes[0]-1));
    double d0   = (2.0*M_PI)/(double(gridSizes[1]));
    double r    = r0 + ij[0]*dr;
    double a0   = ij[1]*d0;

    xy[0] = r*std::cos(a0);
    xy[1] = r*std::sin(a0);  
    
    return(0);
  };
  
  template<int numDim,size_t numX,size_t numY,size_t numZ>
  int CurvilinearGrid1(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    //  double scale[] = {CLSIZE1,CLSIZE2,CLSIZE3};
    double r0 = 0.25;
    double rmax = 1.25;
    double rt   = 2.0;
    double dr   = (rmax - r0)/(double(numX-1));
    double d0   = (2.0*M_PI)/(double(numY));
    double dp   = (2.0*M_PI)/(double(numZ));
    double r  = r0 + ijk[0]*dr;
    double a0 = ijk[1]*d0;
    
    
    double phi = ijk[2]*dp;
    xyz[0] = (rt + r*std::cos(a0))*std::sin(phi);
    xyz[1] = (rt + r*std::cos(a0))*std::cos(phi);
    xyz[2] = r*std::sin(a0);
    
    return(0);
  };

  template<int numDim,size_t numX,size_t numY>
  int CurvilinearGrid2(const std::vector<size_t> &ijk,std::vector<double> &xyz)
  {
    //  double scale[] = {CLSIZE1,CLSIZE2,CLSIZE3};
    
    double I = (ijk[0] + 1)/numX;
    double J = (ijk[1] + 1)/numY;
    
    if(numDim == 2){
      xyz[0] = I*I + J;
      xyz[1] = 2.0*J*J + I;
    } else if (numDim == 3){
      double K = ijk[2] + 1;
      xyz[0] = I*I + 2.0*J*J + K;
      xyz[1] = J*J + 2.0*K*K + I;
      xyz[2] = 2.0*I*I + 3.0*K*K + J;
    }
    return(0);
  };

  

}
//   template<int numDim,std::vector<double> scale>
//   int RectilinearMetric1(const std::vector<size_t> &ijk,
//                          std::vector<double> &xyz,double &jacm1)
//   {
//     //    double scale[] = {101,51,21};
//     double dxDXi  = 2.0*(ijk[0]+1)/(scale[0]*scale[0]);
//     double dyDEta = 4.0*(ijk[1]+1)/(scale[1]*scale[1]);
//     xyz[0] = dyDEta;
//     xyz[1] = dxDXi;
//     jacm1 = xyz[0]*xyz[1];
//     if (numDim == 3){
//       double dzDZeta = 6.0*(ijk[2]+1)/(scale[2]*scale[2]);
//       xyz[0] = dyDEta*dzDZeta;
//       xyz[1] = dxDXi*dzDZeta;
//       xyz[2] = dxDXi*dyDEta;
//       jacm1 = xyz[0]*xyz[1]*xyz[2];
//     }
//     return(0);
//   };