added new interpolator : necessary for efficient interpolation from pointinterpolator to variable

program/source/extemp/variable.h

@@ -1449,7 +1449,7 @@ void Variable<DTyp>::operator= ( const Expr<A>& a ) {
 			for(int j = 1;j < Ny;++j )
 				for (int i = 1;i < Nx;++i ) {
 					data_quadrangles[id][i+ ( Nx+1 ) *j] =
-					    ao.template Give_data<quadrangleEl> ( id, i, j, 0, Nx, 0 );
+                        ao.template Give_data<quadrangleEl> ( id, i, j, 0, Nx, 0 );
 				}
 		}
 		Update<quadrangleEl> ( id );

program/source/grid/blockgrid.cc

@@ -766,6 +766,15 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
 
 D3vector Blockgrid::Give_coord_quadrangle ( int id, int i, int j ) const
 {
+    if(bg_coord != NULL) if(bg_coord->blockgrid_quad_coordinates_calculated)
+        {
+            int Nx = Give_Nx_quadrangle(  id ) +1;
+            //cout << "id_hex " << id_hex << " i " << i << " j " << j << " k " << k << endl;
+            //cout << "index " << i + Nx* (j  + k * (Ny))<< endl;
+            return bg_coord->blockgrid_quad_coordinates.at(id).at(i + Nx* (j  )  );
+        }
+
+
   double eta, xi;
   D3vector PSW, PSE, PNW, PNE, PW, PE, PN, PS;
   D3vector PTransFace, PWithoutTrans;
@@ -968,6 +977,10 @@ D3vector Blockgrid::Give_coord_edge ( int id, int i ) const
 
 D3vector Blockgrid::Give_coord_point ( int id_point ) const
 {
+    if(bg_coord != NULL) if(bg_coord->blockgrid_point_coordinates_calculated)
+        {
+            return bg_coord->blockgrid_point_coordinates.at(id_point);
+        }
   return ug->Give_point ( id_point )->Give_coordinate();
 }
 
@@ -1028,6 +1041,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
 {
     blockgrid_edge_coordinates_calculated = false;
     blockgrid_hexa_coordinates_calculated = false;
+    blockgrid_quad_coordinates_calculated = false;
+    blockgrid_point_coordinates_calculated = false;
 
     blockgrid_hexa_boundaries_calculated = false;
     blockgrid_hexa_coordinates.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
@@ -1067,8 +1082,40 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
             }
     }
 
+    blockgrid_quad_coordinates.resize(bg->Give_unstructured_grid()->Give_number_quadrangles());
+    for (int id_quad = 0 ; id_quad < bg->Give_unstructured_grid()->Give_number_quadrangles() ; id_quad++)
+    {
+        int Nx = bg->Give_Nx_quadrangle(id_quad)+1;
+        int Ny = bg->Give_Ny_quadrangle(id_quad)+1;
+        blockgrid_quad_coordinates.at(id_quad).resize(Nx*Ny);
+    }
+    for (int id_quad = 0 ; id_quad < bg->Give_unstructured_grid()->Give_number_quadrangles() ; id_quad++)
+    {
+        int Nx = bg->Give_Nx_quadrangle(id_quad)+1;
+        int Ny = bg->Give_Ny_quadrangle(id_quad)+1;
+
+        for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j)
+        {
+                blockgrid_quad_coordinates.at(id_quad).at(i + Nx * j) = bg->Give_coord_quadrangle(id_quad,i,j);
+        }
+    }
+
+    blockgrid_point_coordinates.resize(bg->Give_unstructured_grid()->Give_number_points());
+    for (int id_point = 0 ; id_point < bg->Give_unstructured_grid()->Give_number_points() ; id_point++)
+    {
+        int N = bg->Give_unstructured_grid()->Give_number_points();
+
+        for(int i=0;i<N;++i)
+        {
+                blockgrid_point_coordinates.at(i) = bg->Give_coord_point(i);
+        }
+    }
+
+
     blockgrid_edge_coordinates_calculated = true;
     blockgrid_hexa_coordinates_calculated = true;
+    blockgrid_quad_coordinates_calculated = true;
+    blockgrid_point_coordinates_calculated = true;
 }
 
 void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()

program/source/grid/blockgrid.h

@@ -56,18 +56,28 @@ class Blockgrid_coordinates {
 
     std::vector< std::vector< D3vector > > * getBlockgridHexaCoordinates(){return &blockgrid_hexa_coordinates;}
     std::vector< std::vector< D3vector > > * getBlockgridEdgeCoordinates(){return &blockgrid_edge_coordinates;}
+    std::vector< D3vector > * getBlockgridPointCoordinates(){return &blockgrid_point_coordinates;}
+    std::vector< std::vector< D3vector > > * getBlockgridQuadCoordinates(){return &blockgrid_quad_coordinates;}
+
+
     void setBlockgridHexaCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_hexa_coordinates = coords;}
     void setBlockgridEdgeCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_edge_coordinates = coords;}
+    void setBlockgridPointCoordinates(std::vector< D3vector >  coords)             {blockgrid_point_coordinates = coords;}
+    void setBlockgridQuadCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_quad_coordinates = coords;}
 
 
     std::vector< std::vector< D3vector > > blockgrid_hexa_coordinates;
     std::vector< std::vector< D3vector > > blockgrid_edge_coordinates;
+    std::vector< std::vector< D3vector > > blockgrid_quad_coordinates;
+    std::vector< D3vector >  blockgrid_point_coordinates;
 
 
     Blockgrid *bg;
 
     bool blockgrid_edge_coordinates_calculated;
     bool blockgrid_hexa_coordinates_calculated;
+    bool blockgrid_quad_coordinates_calculated;
+    bool blockgrid_point_coordinates_calculated;
     bool blockgrid_hexa_boundaries_calculated;
     std::vector< std::vector< bool > > blockgrid_edge_coordinates_set;
 

program/source/interpol/interpol.cc

@@ -682,6 +682,7 @@ Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
 }
   
 
+
   
 void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<double>* U_to,
 					    double defaultInterpolation) {
@@ -706,6 +707,8 @@ void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<d
     
     (*U_to) = myFunctor(Xc,Yc,Zc);
 }
+
+
  
 double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, double coord_z) {  
   if(coord_x > pENT.x) return 0.0;
@@ -723,8 +726,7 @@ double Interpolate_on_block_grid::evaluate(double coord_x, double coord_y, doubl
   if(i>=nx-1) i=nx-2;  if(j>=ny-1) j=ny-2;  if(k>=nz-1) k=nz-2;
   
   //cout << "i: " << i << " j: " << j << " k: " << k << endl;
-    
-    
+
   double uWSD = data[i    +nx*(j    +ny* k)];
   double uESD = data[(i+1)+nx*(j    +ny* k)];
   double uWND = data[i    +nx*((j+1)+ny* k)];
@@ -874,6 +876,31 @@ PointInterpolator::PointInterpolator(Interpolate_on_structured_grid* intermediat
     hz = (pENT.z - pWSD.z) / (nz-1);
 
     intermediateGrid->interpolate<double>(*U_from,data,defaultInterpolation_);
+}
+
+PointInterpolator::PointInterpolator(Interpolate_on_structured_grid *intermediateGrid, double defaultInterpolation_)
+{
+    defaultInterpolation = defaultInterpolation_;
+
+    nx = intermediateGrid->nx;
+    ny = intermediateGrid->ny;
+    nz = intermediateGrid->nz;
+
+    data = new double[nx*ny*nz];
+    dataCounter = new int[nx*ny*nz];
+    for (int iter = 0 ; iter < nx*ny*nz;iter++)
+    {
+        data[iter]=0.0;
+        dataCounter[iter]=0;
+    }
+
+    pENT = intermediateGrid->pENT;
+    pWSD = intermediateGrid->pWSD;
+
+    hx = intermediateGrid->getHx();
+    hy = intermediateGrid->getHy();
+    hz = intermediateGrid->getHz();
+
 }
  
  /*
@@ -942,7 +969,6 @@ double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_
   //cout << "id: " << id << " jd: " << jd << " kd: " << kd << endl;
   //cout << "i: " << i << " j: " << j << " k: " << k << endl;
   //cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
-        
 
 
   double uWSD = data[i    +nx*(j    +ny* k)];
@@ -1194,6 +1220,17 @@ void PointInterpolator::smoothGrid()
     }
 }
 
+void PointInterpolator::normToNumberOfWritings()
+{
+    for (int iter = 0 ; iter < nx*ny*nz;iter++)
+    {
+        if (dataCounter[iter] >0)
+        {
+            data[iter] = data[iter] / double(dataCounter[iter]);
+        }
+    }
+}
+
 void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value)
 {
     coord_x-=shiftx;
@@ -1228,23 +1265,39 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y,
 
 
 
-    double uWSD = value * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ);
-    double uESD = value *        locX  * (1.0 - locY) * (1.0 - locZ);
-    double uWND = value * (1.0 - locX) *        locY  * (1.0 - locZ);
-    double uEND = value *        locX  *        locY  * (1.0 - locZ);
-    double uWST = value * (1.0 - locX) * (1.0 - locY) *        locZ ;
-    double uEST = value *        locX  * (1.0 - locY) *        locZ ;
-    double uWNT = value * (1.0 - locX) *        locY  *        locZ ;
-    double uENT = value *        locX  *        locY  *        locZ ;
-
-    data[i    +nx*(j    +ny* k)] = uWSD;
-    data[(i+1)+nx*(j    +ny* k)] = uESD;
-    data[i    +nx*((j+1)+ny* k)] = uWND;
-    data[(i+1)+nx*((j+1)+ny* k)] = uEND;
-    data[i    +nx*(j    +ny*(k+1))] = uWST;
-    data[(i+1)+nx*(j    +ny*(k+1))] = uEST;
-    data[i    +nx*((j+1)+ny*(k+1))] = uWNT;
-    data[(i+1)+nx*((j+1)+ny*(k+1))] = uENT;
+//    double uWSD = value * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ);
+//    double uESD = value *        locX  * (1.0 - locY) * (1.0 - locZ);
+//    double uWND = value * (1.0 - locX) *        locY  * (1.0 - locZ);
+//    double uEND = value *        locX  *        locY  * (1.0 - locZ);
+//    double uWST = value * (1.0 - locX) * (1.0 - locY) *        locZ ;
+//    double uEST = value *        locX  * (1.0 - locY) *        locZ ;
+//    double uWNT = value * (1.0 - locX) *        locY  *        locZ ;
+//    double uENT = value *        locX  *        locY  *        locZ ;
+
+//    data[i    +nx*(j    +ny* k)] = uWSD;
+//    data[(i+1)+nx*(j    +ny* k)] = uESD;
+//    data[i    +nx*((j+1)+ny* k)] = uWND;
+//    data[(i+1)+nx*((j+1)+ny* k)] = uEND;
+//    data[i    +nx*(j    +ny*(k+1))] = uWST;
+//    data[(i+1)+nx*(j    +ny*(k+1))] = uEST;
+//    data[i    +nx*((j+1)+ny*(k+1))] = uWNT;
+//    data[(i+1)+nx*((j+1)+ny*(k+1))] = uENT;
+        data[i    +nx*(j    +ny* k)] += value;
+        data[(i+1)+nx*(j    +ny* k)] += value;
+        data[i    +nx*((j+1)+ny* k)] += value;
+        data[(i+1)+nx*((j+1)+ny* k)] += value;
+        data[i    +nx*(j    +ny*(k+1))] += value;
+        data[(i+1)+nx*(j    +ny*(k+1))] += value;
+        data[i    +nx*((j+1)+ny*(k+1))] += value;
+        data[(i+1)+nx*((j+1)+ny*(k+1))] += value;
+        dataCounter[i    +nx*(j    +ny* k)]++;
+        dataCounter[(i+1)+nx*(j    +ny* k)]++;
+        dataCounter[i    +nx*((j+1)+ny* k)]++;
+        dataCounter[(i+1)+nx*((j+1)+ny* k)]++;
+        dataCounter[i    +nx*(j    +ny*(k+1))]++;
+        dataCounter[(i+1)+nx*(j    +ny*(k+1))]++;
+        dataCounter[i    +nx*((j+1)+ny*(k+1))]++;
+        dataCounter[(i+1)+nx*((j+1)+ny*(k+1))]++;
 
     return;
 
@@ -1339,7 +1392,7 @@ void PointInterpolator::scaleInterpolatedData(double scale, double zeroVal)
 
  
 PointInterpolator::~PointInterpolator() {
-    delete interpolatorStructured;
+    //delete interpolatorStructured;
     delete[] data;
 }
 
@@ -3838,3 +3891,28 @@ int Interpolate_direct::checkForHexaNeighbours(int idHex, int i, int j, int k, D
     counterHexa--;
     return typ;
 }
+
+Interpolate_on_block_grid_from_pointinterpolator::Interpolate_on_block_grid_from_pointinterpolator(PointInterpolator *interp, Blockgrid *blockgrid_to_)
+{
+
+    blockgrid_to = blockgrid_to_;
+    pointInterpolator = interp;
+
+}
+
+Interpolate_on_block_grid_from_pointinterpolator::~Interpolate_on_block_grid_from_pointinterpolator()
+{
+    delete pointInterpolator;
+    delete[] data;
+}
+
+void Interpolate_on_block_grid_from_pointinterpolator::interpolate(Variable<double> *U_to, double defaultInterpolation)
+{
+    Functor3<double,double,PointInterpolator> myFunctor(pointInterpolator);
+
+    X_coordinate Xc(*blockgrid_to);
+    Y_coordinate Yc(*blockgrid_to);
+    Z_coordinate Zc(*blockgrid_to);
+
+    (*U_to) = myFunctor(Xc,Yc,Zc);
+}

program/source/interpol/interpol.h

@@ -76,7 +76,15 @@ class Interpolate_on_structured_grid {
 **/
   template <class DTyp>
     void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
-  public:
+
+    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:
   int nx, ny, nz;
   D3vector pENT,pWSD;
   private:
@@ -200,6 +208,34 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
 /**
  * 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, double defaultInterpolation = 0.0);
+
+ 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:
 /**
@@ -210,7 +246,8 @@ class Interpolate_on_block_grid {
  @param  blockgrid_ of unstructured grid
 **/
   Interpolate_on_block_grid(int nx_, int ny_, int nz_,
-		            Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_);
+                    Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_);
+  Interpolate_on_block_grid(PointInterpolator *interp, Blockgrid* blockgrid_to_ );
   ~Interpolate_on_block_grid();
   
   /**
@@ -220,13 +257,15 @@ class Interpolate_on_block_grid {
    * 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;
   
@@ -280,7 +319,7 @@ class PointInterpolator {
   PointInterpolator(Interpolate_on_structured_grid* intermediateGrid,
 		    Variable<double>* U_from, double defaultInterpolation_ = 0.0);  
   
-
+  PointInterpolator(Interpolate_on_structured_grid* intermediateGrid, double defaultInterpolation_ = -1.0);
 
 
 
@@ -306,20 +345,22 @@ class PointInterpolator {
   std::vector<double> evaluateGradient(double coord_x, double coord_y, double coord_z);
   void smoothGrid();
 		 
-
+  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);
-  
+
+  double shiftx, shifty, shiftz;
+  double rotationx, rotationy, rotationz;
 private:
   int nx, ny, nz;
   double hx, hy, hz;
-  double shiftx, shifty, shiftz;
   Interpolate_on_structured_grid* interpolatorStructured;
   double* data;
+  int* dataCounter;
   
   D3vector pWSD;
   D3vector pENT;