interpol.h

/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/

// ------------------------------------------------------------
//
// variable.h
//
// ------------------------------------------------------------

#ifndef INTERPOL_H
#define INTERPOL_H

/////////////////////////////////////////////////////////////
// 1. Interpolate from  blockgrid to rectangular blockgrid
// 2. Interpolate from  blockgrid  to  blockgrid
// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid
// 4. Interpolate from blockgrid direct to any point
/////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////
// 1. Interpolate from  blockgrid to rectangular blockgrid
/////////////////////////////////////////////////////////////

/** \addtogroup InterpolationOperators **/
/* @{ */  
/**
 * Interpolation 3D:
data structure on structured grid
 \verbatim
ny    *  *  *  *  *  *
ny-1  *  +  +  +  +  *
3     *  +  +  +  +  *
2     *  +  +  +  +  *
1     *  +  +  +  +  *
0     *  *  *  *  *  *
      0  1  2  3 nx-1 nx
  \endverbatim      
 * interpolates data from blockgrid_ to rectangular block [pWSD, pENT]
**/
class Interpolate_direct;
class Interpolate_on_structured_grid {
 public:
/**
 * preparation for interpolation
 @param  nx_ number of structured grid points x-direction
 @param  ny_ number of structured grid points y-direction
 @param  nz_ number of structured grid points z-direction
 @param  pWSD Corner WSD of structured grid
 @param  pENT Corner WSD of structured grid
 @param  blockgrid_ of unstructured grid
**/
  Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
				 D3vector pWSD, D3vector pENT,
                 Blockgrid& blockgrid_, bool trilinearInterpolationFlag_ = false);
  Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
                 Blockgrid& blockgrid_, bool trilinearInterpolationFlag_ = false);

  ~Interpolate_on_structured_grid();
/**
 * interpolates from Variable u(of unstructured blockgrid) to structured data
 @param  u:  Variable on unstructured blockgrid
 @param  data: pointer to structured grid data i+nx*(j+ny*k)
**/
  template <class DTyp>
    void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);

    void update_Interpolate_on_structured_grid(Blockgrid& blockgrid_, bool onlyOnSurfaceZ = false);
    void setInterpolateDirect(Interpolate_direct *id_){id = id_;}

    int getNx(){return nx;}
    int getNy(){return ny;}
    int getNz(){return nz;}

    double getHx(){return hx;}
    double getHy(){return hy;}
    double getHz(){return hz;}
public:

    D3vector trilinarInterpolation(D3vector X, int id_Hex, int i,  int j,  int k );
    bool trilinearInterpolationFlag{true};
  int nx, ny, nz;
  D3vector pENT,pWSD;
  private:
  int* ids_hex;
  int *ids_i, *ids_j, *ids_k;
  int *typ_tet;

  D3vector *lambda;

  
  double hx, hy, hz;

  Blockgrid*        blockgrid;
  Unstructured_grid*       ug;

  Interpolate_direct* id{NULL};
};
/* @} */


template <class DTyp>
class Give_corner_data_of_cube {
 public:
  Give_corner_data_of_cube(Variable<DTyp>& u_, int hex_, int i_, int j_, int k_);
  
  DTyp WSD() { return u->template Give_data<hexahedronEl>(id, i,   j,   k,   Nx, Ny); }
  DTyp ESD() { return u->template Give_data<hexahedronEl>(id, i+1, j,   k,   Nx, Ny); }
  DTyp WND() { return u->template Give_data<hexahedronEl>(id, i,   j+1, k,   Nx, Ny); }
  DTyp END() { return u->template Give_data<hexahedronEl>(id, i+1, j+1, k,   Nx, Ny); }
  DTyp WST() { return u->template Give_data<hexahedronEl>(id, i,   j,   k+1, Nx, Ny); }
  DTyp EST() { return u->template Give_data<hexahedronEl>(id, i+1, j,   k+1, Nx, Ny); }
  DTyp WNT() { return u->template Give_data<hexahedronEl>(id, i,   j+1, k+1, Nx, Ny); }
  DTyp ENT() { return u->template Give_data<hexahedronEl>(id, i+1, j+1, k+1, Nx, Ny); }

 private:
  int id, i, j, k;
  int Nx, Ny;
  Variable<DTyp>* u;
};

template <class DTyp>
Give_corner_data_of_cube<DTyp>::Give_corner_data_of_cube(Variable<DTyp>& u_,
							 int hex_, int i_, int j_, int k_) {
  Blockgrid*        blockgrid;

  u = &u_;
  blockgrid = u->Give_blockgrid();
  id = hex_; i = i_;  j = j_;  k = k_;
  
  Nx = blockgrid->Give_Nx_hexahedron(id);
  Ny = blockgrid->Give_Ny_hexahedron(id);
}


template <class DTyp>
void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation) {
  int i,j,k, id_hex, typ;
  int ind_global;



  for(int id=0;id<ug->Give_number_hexahedra();++id)
    u.template Update<hexahedronEl>(id);

  for(int ks = 0; ks < nz;++ks)
    for(int js = 0; js < ny;++js)
      for(int is = 0; is < nx;++is) {
    ind_global = is+nx*(js+ny*ks);
    i = ids_i[ind_global];
    j = ids_j[ind_global];
    k = ids_k[ind_global];
    id_hex = ids_hex[ind_global];

    /*
    // test : das geht nicht::
    if(id_hex < 0) {
      ind_global = is+nx*((ny-1-js)+ny*ks);
      i = ids_i[ind_global];
      j = ids_j[ind_global];
      k = ids_k[ind_global];
      id_hex = ids_hex[ind_global];
    }
*/
	if(id_hex < 0) data[ind_global] = defaultInterpolation;
	else {
	  Give_corner_data_of_cube<DTyp> du(u, id_hex, i,j,k);

	  typ = typ_tet[ind_global];

	
	  /*
	  if( ind_global == 3)
	    cout << " \n ind_global: " << is + nx * (js + ny * ks)
	         << " typ: " << typ
	         << " id_hex " << id_hex
	         << " i: " << i
	         << " j: " << j
	         << " k: " << k
	         << endl;
  	  */

      if(typ==0) data[ind_global] = interpolate_in_tet(lambda[ind_global],
							   du.WND(),du.WNT(),du.WST(),du.EST());
      if(typ==1) data[ind_global] = interpolate_in_tet(lambda[ind_global],
							   du.EST(),du.WND(),du.WST(),du.ESD());
      if(typ==2) data[ind_global] = interpolate_in_tet(lambda[ind_global],
							   du.WND(),du.WSD(),du.WST(),du.ESD());
      if(typ==3) data[ind_global] = interpolate_in_tet(lambda[ind_global],
                               du.EST(),du.WND(),du.ESD(),du.END());
      if(typ==4) data[ind_global] = interpolate_in_tet(lambda[ind_global],
							   du.ENT(),du.WNT(),du.EST(),du.END());
      if(typ==5) data[ind_global] = interpolate_in_tet(lambda[ind_global],
							   du.WNT(),du.WND(),du.EST(),du.END());
      if(typ==6) data[ind_global] = interpolate_in_tet_trilinear(lambda[ind_global],du.END(),du.ESD(),du.WSD(),du.WND(),
                                                                 du.WNT(),du.WST(),du.EST(),du.ENT());
	}
      }
}

/////////////////////////////////////////////////////////////
// 2. Interpolate from  blockgrid  to  blockgrid
/////////////////////////////////////////////////////////////

/** \addtogroup InterpolationOperators **/
/* @{ */  
/**
 * interpolates data from blockgrid_from to blockgrid_to_ using an internal rectangular grid if size nx,ny,nz
**/
class PointInterpolator;

class Interpolate_on_block_grid_from_pointinterpolator {
public:
/**
* If a pointinterpolator is accesible, it can be used directly to interpolate on blockgrid_to
@param  blockgrid_ of unstructured grid
**/

 Interpolate_on_block_grid_from_pointinterpolator(PointInterpolator *interp, Blockgrid* blockgrid_to_ );
 //~Interpolate_on_block_grid_from_pointinterpolator();
 void interpolate(Variable<double>* U_to, Boundary_Marker *marker = NULL);

 double evaluate(double coord_x, double coord_y, double coord_z);


private:
 int nx, ny, nz;
 double hx, hy, hz;
 PointInterpolator *pointInterpolator;
// double* data;
 Blockgrid* blockgrid_to;

 D3vector pWSD;
 D3vector pENT;

};

class Interpolate_on_block_grid {
 public:
/**
 * preparation for interpolation
 @param  nx_ number of structured grid points x-direction
 @param  ny_ number of structured grid points y-direction
 @param  nz_ number of structured grid points z-direction
 @param  blockgrid_ of unstructured grid
**/
  Interpolate_on_block_grid(int nx_, int ny_, int nz_,
                    Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_);
  Interpolate_on_block_grid(PointInterpolator *interp, Blockgrid* blockgrid_to_ );
  ~Interpolate_on_block_grid();
  
  /**
   * interpoliert Daten. Falls an einem Punkt nicht interpoliert werden kann
   * da U_from keine Zelle hat, dann wird
   *     defaultInterpolation
   * verwendet.
   **/ 
  void interpolate(Variable<double>* U_from, Variable<double>* U_to, double defaultInterpolation = 0.0);

  double evaluate(double coord_x, double coord_y, double coord_z);


private:
  int nx, ny, nz;
  double hx, hy, hz;
  Interpolate_on_structured_grid* interpolatorStructured;
  PointInterpolator *pointInterpolator;
  double* data;
  Blockgrid* blockgrid_to;
  
  D3vector pWSD;
  D3vector pENT;
  
};
/* @} */


/////////////////////////////////////////////////////////////
// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid  
/////////////////////////////////////////////////////////////

/*
class Intermadiate_grid_for_PointInterpolator : public Interpolate_on_structured_grid
{
public:
  Intermadiate_grid_for_PointInterpolator(int nx_, int ny_, int nz_, Variable<double>* U_from);
  int nx, ny, nz;
  Interpolate_on_structured_grid* interpolatorStructured(int nx_, int ny_, int nz_,D3vector pWSD, D3vector pENT,
		                                          Variable<double>* U_from);  
  
  D3vector pWSD;
  D3vector pENT;
};
*/

/** \addtogroup InterpolationOperators **/
/* @{ */  
/**
 * interpolates data from blockgrid_from to blockgrid_to_ using an internal rectangular grid if size nx,ny,nz
***/
class PointInterpolator {
  friend Interpolate_on_structured_grid;
 public:
/**
 * preparation for interpolation
 @param  nx_ number of structured grid points x-direction
 @param  ny_ number of structured grid points y-direction
 @param  nz_ number of structured grid points z-direction
 @param  blockgrid_ of unstructured grid
**/
  PointInterpolator(int nx_, int ny_, int nz_,
            Variable<double>* U_from, double defaultInterpolation_ = 0.0, bool trilinearInterpolation_ = false);
  
  PointInterpolator(int nx_, int ny_, int nz_,
	            D3vector pWSD, D3vector pENT,
		    Variable<double>* U_from, double defaultInterpolation_ = 0.0);  
  
  PointInterpolator(Interpolate_on_structured_grid* intermediateGrid,
		    Variable<double>* U_from, double defaultInterpolation_ = 0.0);  
  
  PointInterpolator(Interpolate_on_structured_grid* intermediateGrid, double defaultInterpolation_ = -1.0, bool counter = false);



  /**
   * Calculates an intermediate Grid for above constructor
   */
  Interpolate_on_structured_grid* intermediateGrid;  

  
  ~PointInterpolator();
  
  
  
  /**
   * interpoliert Daten. Falls an einem Punkt nicht interpoliert werden kann
   * da U_from keine Zelle hat, dann wird
   *     defaultInterpolation
   * verwendet.
   * @param coord_x, coord_y, coord_z, Koordinaten des Punktes
   * @return interpolierter Wert
   **/

  void updateVariable(Variable<double>* U_from);
  void updatePointInterpolator(Interpolate_on_structured_grid* intermediateGrid);
  double evaluate(double coord_x, double coord_y, double coord_z);
  std::vector<double> evaluateGradient(double coord_x, double coord_y, double coord_z);
  void smoothGrid();
  void resetInterpolator();

  void normToNumberOfWritings();
  void writeOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value);
  void subtractOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value);
  void shiftInterpolatedGrid(double coord_x, double coord_y, double coord_z);
  void scaleInterpolatedData(double scale, double zeroVal = 0.0);

  void QPrint_VTK(QString DateiName);

  D3vector pWSD;
  D3vector pENT;
  int nx, ny, nz;
  double hx, hy, hz;
private:
  double shiftx{0}, shifty{0}, shiftz{0};
  double rotationx{0}, rotationy{0}, rotationz{0};

  Interpolate_on_structured_grid* interpolatorStructured;
  double* data;
  bool dataCounterFlag{false};
  int* dataCounter;

  
  double defaultInterpolation;
  
};
/* @} */


/////////////////////////////////////////////////////////////
// 4. Interpolate from blockgrid direct to any point
/////////////////////////////////////////////////////////////


/** \addtogroup InterpolationOperators **/
/* @{ */
/**
 * Interpolate from blockgrid direct to any point
***/
class Interpolate_direct {
 public:
    Interpolate_direct (Blockgrid* bg):blockgrid(bg){}

    /**
     * preparation for interpolation : calculates the neighbour index for each cell.
     * also  calculates the neighbour cell, if in other blockgrid.
    **/
    void init();
    void updateBoundaryBoxes();

    /**
     * Iterates throuh all cells : if D3vector v is inside of a cell, the weighting of the 8 cell points are calculated (D3vector lambda).
     * Uses previously cell indexes and its neighbours, since it's assumed that two consecutive find_cell(v) calls are close to each other.
     * Saves : idNow, idPrev, ifPrevPrev, idPrevPrevPrev
     *
    **/
    void find_cell(D3vector v);
    void find_surface_cell(D3vector v);
    /**
     * Print a vtk file with a box, surrounding the cell.
    **/
    void writeBox(D3vector vWSD, D3vector vENT, std::string str);
    /**
     * Print a vtk file with a point.
    **/
    static void writePoint(D3vector v, std::string str, int Counter);
    /**
     * Writes interpolation efficiency: how many direct hits, second try, iterate through all, ect...
    **/
   void writeInterpolationEfficiency();
   //   void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
    template <class DTyp>
    DTyp evaluate(Variable<DTyp>& u);

    std::vector<std::vector<std::vector<int> > > getBoundaryBox(){return array_box_boundary;}
    std::vector<std::vector<std::vector<D3vector> > > * getArrayBox() {return &arrayBoxWSDENT;}
    void setArrayBox(std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT_){arrayBoxWSDENT = arrayBoxWSDENT_;}



    int counterFast{}, counterFastest{}, counterSamePoint{}, counterSecondTry{}, counterThirdTry{}, counterSlow{}, counterHexa{}, counterCorner{}, counterEdge{};
    int checkCounter{};
    int boxCounter{};
    int typCounter0{},typCounter1{},typCounter2{},typCounter3{},typCounter4{},typCounter5{};
    bool debugTest;
    bool badLambdaFound;
    bool trilinearInterpolationFlag{false};
    D3vector lambda;
    D3vector vNow{1e10,1e10,1e10}, vPrev, vPrevPrev, vPrevPrevPrev;

    int getI(){return iNow;}
    int getJ(){return jNow;}
    int getK(){return kNow;}
    D3vector getLam(){return lambda;}
    int getTet(){return typ_tet;}
    int getHex(){return idHexNow;}

    D3vector bilinearInterpolation(D3vector v,D3vector a,D3vector b,D3vector c,D3vector d);
    D3vector trilinarInterpolation(D3vector v, int id_Hex,int i,  int j,  int k );
        bool vectorInBox(D3vector vWSD, D3vector vENT, D3vector v, double eps = 1e-10);
        int checkBox(int idHex, int i, int j, int k, D3vector v);
        int checkBoxSurface(int idHex, int i, int j, int k, D3vector v);

    private:
      int idHexPrev{-1}, iPrev{-1}, jPrev{-1}, kPrev{-1};
      int idHexPrevPrevPrev{-1}, iPrevPrevPrev{-1}, jPrevPrevPrev{-1}, kPrevPrevPrev{-1};
      int idHexPrevPrev{-1}, iPrevPrev{-1}, jPrevPrev{-1}, kPrevPrev{-1};
      int idHexNow{-1}, iNow{-1}, jNow{-1}, kNow{-1};

      double lamLowerLimit{-0.1};
      double lamUpperLimit{1.1};
      int typ_tet;

      D3vector cWSD, cESD;
      D3vector cWND, cEND;

      D3vector cWST, cEST;
      D3vector cWNT, cENT;
      Blockgrid    *blockgrid;
      std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
      std::vector<std::vector<bool > > isOuterBoundaryFlag;

      std::vector<std::vector<std::vector<int> > > array_box_boundary;
      std::vector<std::vector<std::vector<int> > > array_point_boundary;
private:
      /**
       * Necessary to calculate the neighbour relation.
      **/
      std::vector<std::vector<int> > calculateNeighbourIndexRelation(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer);
      std::vector<int> calculateNeighbourIndex(std::vector<std::vector<int> > relation, int id_hex_outside,int id_hex_inside, int i, int j, int k, int Nx, int Ny, int Nz);
      std::vector<std::vector<int> > filterCorrectNeighbours(std::vector<std::vector<int> > outer);
      bool compareIndicies(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer, int notCheck);
      std::vector<std::vector<int> > switchIJ(std::vector<std::vector<int> > v);
      std::vector<std::vector<int> > switchIK(std::vector<std::vector<int> > v);
      std::vector<std::vector<int> > switchJK(std::vector<std::vector<int> > v);
      std::vector<std::vector<int> > invertI (std::vector<std::vector<int> > v);
      std::vector<std::vector<int> > invertJ (std::vector<std::vector<int> > v);
      std::vector<std::vector<int> > invertK (std::vector<std::vector<int> > v);

      /**
       * Necessary to evaluate the neighbour cell relation.
      **/
      void setPrevIndex();
      void setPrevPrevIndex();
      void setPrevPrevPrevIndex();
      int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
      int checkForHexaNeighboursSurface(int idHex, int i, int j, int k, D3vector v);
      int checkBoxSurrounding(int idHex, int i, int j, int k, D3vector v);
      int checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v, bool neglectBadLambda = false);
      int checkBoxSurroundingSurface(int idHex, int i, int j, int k, D3vector v, bool neglectBadLambda = false);
      int checkCorner(int idHex, int i, int j, int k, D3vector v);
      int checkCornerSurface(int idHex, int i, int j, int k, D3vector v);
      int checkEdge(int idHex, int i, int j, int k, D3vector v);
      int checkEdgeSurface(int idHex, int i, int j, int k, D3vector v);
      bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
      double calculateOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);




};

/* @} */
template<class DTyp>
DTyp Interpolate_direct::evaluate(Variable<DTyp> &u)
{
 int i,j,k, id_hex, typ;
 int ind_global;

 DTyp returnVal;

// for(int id=0;id<5;++id)
//   u.template Update<hexahedronEl>(id);


   if(idHexNow < 0)
   {returnVal = -1;}
   else {
     Give_corner_data_of_cube<DTyp> du(u, idHexNow , iNow,jNow,kNow);

     int typ = typ_tet;

     if (lambda == D3vector(0,0,0))
     {
         cout << "no point found " << endl;
     }
     //cout << "lambda : " ; lambda.Print();cout<<endl;
     if(typ==0) returnVal = interpolate_in_tet(lambda,
                              du.WND(),du.WNT(),du.WST(),du.EST());
     if(typ==1) returnVal = interpolate_in_tet(lambda,
                              du.EST(),du.WND(),du.WST(),du.ESD());
     if(typ==2) returnVal = interpolate_in_tet(lambda,
                              du.WND(),du.WSD(),du.WST(),du.ESD());
     if(typ==3) returnVal = interpolate_in_tet(lambda,
                              du.EST(),du.WND(),du.ESD(),du.END());
     if(typ==4) returnVal = interpolate_in_tet(lambda,
                              du.ENT(),du.WNT(),du.EST(),du.END());
     if(typ==5) returnVal = interpolate_in_tet(lambda,
                              du.WNT(),du.WND(),du.EST(),du.END());
     if(typ==6) returnVal = interpolate_in_tet_trilinear(lambda,du.END(),du.ESD(),du.WSD(),du.WND(),
                                                           du.WNT(),du.WST(),du.EST(),du.ENT());
     if(typ==7) returnVal = interpolate_in_tet(lambda,du.EST(),du.ESD(),du.ENT(),0.0);
     if(typ==8) returnVal = interpolate_in_tet(lambda,du.END(),du.ESD(),du.ENT(),0.0);
     if(typ==9) returnVal = interpolate_in_tet(lambda,du.WND(),du.END(),du.WSD(),0.0);
     if(typ==10) returnVal = interpolate_in_tet(lambda,du.ESD(),du.END(),du.WSD(),0.0);
     if(typ==11) returnVal = interpolate_in_tet(lambda,du.ENT(),du.WNT(),du.EST(),0.0);
     if(typ==12) returnVal = interpolate_in_tet(lambda,du.WST(),du.WNT(),du.EST(),0.0);
     if(typ==13) returnVal = interpolate_in_tet_bilinear(lambda,du.WND(),du.END(),du.WSD(),du.ESD());
     if(typ==14) returnVal = interpolate_in_tet_bilinear(lambda,du.ENT(),du.WNT(),du.WST(),du.EST());





   }
   //    D3vector cWSD = {1,1,0};      // 1,1,0 : x2
   //    D3vector cESD = {2,0,-1};     // 1,0,0 : x1
   //    D3vector cWND = {0,1,0};      // 0,1,0 : x3
   //    D3vector cEND = {-1,0,-1.0};  // 0,0,0 : x0
   //    D3vector cWST = {1,1,1};      // 1,1,1 : x5
   //    D3vector cEST = {2,0,2};      // 1,0,1 : x6
   //    D3vector cWNT = {0,1,1};      // 0,1,1 : x4
   //    D3vector cENT = {-1,0,2};     // 0,0,1 : x7
    return returnVal;
};


#endif // INTERPOL_H