interpolator expanded : updating now possible if blockgrid is changed (e.g. rotated)

program/source/interpol/interpol.cc

@@ -111,9 +111,241 @@ Interpolate_on_structured_grid::~Interpolate_on_structured_grid() {
   delete[] ids_k;
 
   delete[] typ_tet;
-  delete[] lambda;
+    delete[] lambda;
 }
 
+void Interpolate_on_structured_grid::update_Interpolate_on_structured_grid(Blockgrid &blockgrid_)
+{
+        int Nx, Ny, Nz;
+        int typ;
+
+        int ilmin, jlmin, klmin;
+        int ilmax, jlmax, klmax;
+
+        double factor = 0.1;
+        //  double factor = 0.00001;
+
+        D3vector lam;
+
+
+
+          //Variable<double> coordXYZ(*blockgrid);
+          X_coordinate Xc(blockgrid_);
+          Y_coordinate Yc(blockgrid_);
+          Z_coordinate Zc(blockgrid_);
+          //D3vector pWSD, pENT;
+          pWSD.x = Minimum(Xc);    pWSD.y = Minimum(Yc);    pWSD.z = Minimum(Zc);
+          pENT.x = Maximum(Xc);    pENT.y = Maximum(Yc);    pENT.z = Maximum(Zc);
+
+        blockgrid = &blockgrid_;
+        ug = blockgrid->Give_unstructured_grid();
+
+
+
+        if(nx>1)
+          hx = (pENT.x - pWSD.x) / (nx-1);
+        else
+          hx = 1.0;
+        if(ny>1)
+          hy = (pENT.y - pWSD.y) / (ny-1);
+        else
+          hy = 1.0;
+        if(nz>1)
+          hz = (pENT.z - pWSD.z) / (nz-1);
+        else
+          hz = 1.0;
+
+        int num_total = nx * ny * nz;
+
+        D3vector cWSD, cESD;
+        D3vector cWND, cEND;
+
+        D3vector cWST, cEST;
+        D3vector cWNT, cENT;
+
+        D3vector boxWSD, boxENT;
+
+        D3vector ploc;
+
+//        ids_hex = new int[num_total];
+
+//        ids_i = new int[num_total];
+//        ids_j = new int[num_total];
+//        ids_k = new int[num_total];
+
+//        typ_tet = new int[num_total];
+
+//        lambda = new D3vector[num_total];
+
+        for(int i=0;i<num_total;++i) ids_hex[i] = -1;
+
+        for(int id_hex=0;id_hex<ug->Give_number_hexahedra();++id_hex) {
+            Nx = blockgrid->Give_Nx_hexahedron(id_hex);
+            Ny = blockgrid->Give_Ny_hexahedron(id_hex);
+            Nz = blockgrid->Give_Nz_hexahedron(id_hex);
+
+            for(int k=0;k<Nz;++k)
+          for(int j=0;j<Ny;++j)
+            for(int i=0;i<Nx;++i) {
+                  // corner points of general hex-cell
+              cWSD = blockgrid->Give_coord_hexahedron(id_hex,i,  j,  k  );
+              cESD = blockgrid->Give_coord_hexahedron(id_hex,i+1,j  ,k  );
+              cWND = blockgrid->Give_coord_hexahedron(id_hex,i,  j+1,k  );
+              cEND = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k  );
+
+              cWST = blockgrid->Give_coord_hexahedron(id_hex,i,  j,  k+1);
+              cEST = blockgrid->Give_coord_hexahedron(id_hex,i+1,j  ,k+1);
+              cWNT = blockgrid->Give_coord_hexahedron(id_hex,i,  j+1,k+1);
+              cENT = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k+1);
+
+                  // bounding box calculation
+              boxWSD.x = MIN(MIN(MIN(cWSD.x,cESD.x),MIN(cWND.x,cEND.x)),
+                     MIN(MIN(cWST.x,cEST.x),MIN(cWNT.x,cENT.x))) - factor *hx;
+              boxWSD.y = MIN(MIN(MIN(cWSD.y,cESD.y),MIN(cWND.y,cEND.y)),
+                     MIN(MIN(cWST.y,cEST.y),MIN(cWNT.y,cENT.y))) - factor *hy;
+              boxWSD.z = MIN(MIN(MIN(cWSD.z,cESD.z),MIN(cWND.z,cEND.z)),
+                     MIN(MIN(cWST.z,cEST.z),MIN(cWNT.z,cENT.z))) - factor *hz;
+
+              boxENT.x = MAX(MAX(MAX(cWSD.x,cESD.x),MAX(cWND.x,cEND.x)),
+                     MAX(MAX(cWST.x,cEST.x),MAX(cWNT.x,cENT.x))) + factor *hx;
+              boxENT.y = MAX(MAX(MAX(cWSD.y,cESD.y),MAX(cWND.y,cEND.y)),
+                     MAX(MAX(cWST.y,cEST.y),MAX(cWNT.y,cENT.y))) + factor *hy;
+              boxENT.z = MAX(MAX(MAX(cWSD.z,cESD.z),MAX(cWND.z,cEND.z)),
+                     MAX(MAX(cWST.z,cEST.z),MAX(cWNT.z,cENT.z))) + factor *hz;
+
+              // calculation of indices of a collection of cells of structured grid which contains bounding box
+              ilmin = Ganzzahliger_Anteil((boxWSD.x - pWSD.x) / hx);
+              jlmin = Ganzzahliger_Anteil((boxWSD.y - pWSD.y) / hy);
+              klmin = Ganzzahliger_Anteil((boxWSD.z - pWSD.z) / hz);
+
+
+              ilmax = Ganzzahliger_Anteil((boxENT.x - pWSD.x) / hx);
+              jlmax = Ganzzahliger_Anteil((boxENT.y - pWSD.y) / hy);
+              klmax = Ganzzahliger_Anteil((boxENT.z - pWSD.z) / hz);
+
+              /*
+              cout << " indices: "
+               << " ilmin: " << ilmin
+               << " jlmin: " << jlmin
+               << " klmin: " << klmin
+               << " ilmax: " << ilmax
+               << " jlmax: " << jlmax
+               << " klmax: " << klmax
+               << " boxWSD.x: " << boxWSD.x
+               << " cWSD.x: " << cWSD.x
+               << " Nx: " <<  Nx
+               << endl;
+              */
+              /*
+      bool now;
+      if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.x < 1.0 && boxENT.x > 1.0 ) {
+        cout << "\n \n WSD   : ";  boxWSD.Print();
+        cout << "\n ENT  : ";  boxENT.Print();
+        now = true;
+       }
+       else now = false;
+
+
+       cout << " tt: " << boxWSD.x << " " << pWSD.x << " " << hx << endl;
+       cout << (boxWSD.x - pWSD.x) << endl;
+
+       cout << " z: " << 0.1 << " g: " << Ganzzahliger_Anteil(0.1) << endl;
+       cout << " z: " << -0.1 << " g: " << Ganzzahliger_Anteil(-0.1) << endl;
+
+       cout << " z: " << 5.1 << " g: " << Ganzzahliger_Anteil(5.1) << endl;
+       cout << " z: " << -5.1 << " g: " << Ganzzahliger_Anteil(-5.1) << endl;
+              */
+
+              if(ilmin<0) ilmin=0;
+              if(jlmin<0) jlmin=0;
+              if(klmin<0) klmin=0;
+
+              for(int il = ilmin; (il <= ilmax) && (il < nx);++il)
+                for(int jl = jlmin; (jl <= jlmax) && (jl < ny);++jl)
+              for(int kl = klmin; (kl <= klmax) && (kl < nz);++kl) {
+
+                ploc = D3vector(il * hx, jl * hy, kl * hz) + pWSD;
+
+                //		  cout << "HI" << endl;
+
+                typ = -1;
+
+                lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST);
+                if(contained_in_tet(lam)) typ=0;
+                else {
+                  lam = lambda_of_p_in_tet(ploc,cEST,cWND,cWST,cESD);
+                  if(contained_in_tet(lam)) typ=1;
+                  else {
+                    lam = lambda_of_p_in_tet(ploc,cWND,cWSD,cWST,cESD);
+                    if(contained_in_tet(lam)) typ=2;
+                    else {
+                  lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);
+                  if(contained_in_tet(lam)) typ=3;
+                  else {
+                    lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND);
+                    if(contained_in_tet(lam)) typ=4;
+                    else {
+                      lam = lambda_of_p_in_tet(ploc,cWNT,cWND,cEST,cEND);
+                      if(contained_in_tet(lam)) typ=5;
+                    }
+                  }
+                    }
+                  }
+                }
+
+                /*
+                cout << " typ " << typ << id_hex
+                     << " il: " << il
+                     << " jl: " << jl
+                     << " kl: " << kl
+                     << endl;
+                */
+
+                if(typ!=-1) {
+                  int ind_global;
+                  ind_global = il+nx*(jl+ny*kl);
+                  bool stop;
+                  stop=false;
+
+                  if(ids_hex[ind_global]!=-1) {
+                    stop=new_lam_worse(lambda[ind_global],lam);
+                  }
+
+                  if(stop==false) {
+                    ids_hex[ind_global] = id_hex;
+                    ids_i[ind_global] = i;
+                    ids_j[ind_global] = j;
+                    ids_k[ind_global] = k;
+
+                    typ_tet[ind_global] = typ;
+
+                    lambda[ind_global] = lam;
+                  }
+                  //go_on = false;
+                }
+
+
+              }
+            }
+        }
+
+
+        for(int i=0;i<num_total;++i) {
+          if(ids_hex[i]==-1) {
+            // wir nehmen default value!!
+            /*
+            cout << i
+             << " Error: Interpolate_on_structured_grid: I cannot interpolate all data!"
+             << endl;
+            ids_hex[i] = 0;
+            */
+          }
+          else {
+            //cout << i << " Interpolate_on_structured_grid: o.k.!" << endl;
+          }
+        }
+      }
+
 Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
                                    D3vector pWSD, D3vector pENT,
                                    Blockgrid& blockgrid_) {
@@ -1220,6 +1452,15 @@ void PointInterpolator::smoothGrid()
     }
 }
 
+void PointInterpolator::resetInterpolator()
+{
+    for (int iter = 0 ; iter < nx*ny*nz;iter++)
+    {
+        dataCounter[iter] = 0;
+        data[iter] = defaultInterpolation;
+    }
+}
+
 void PointInterpolator::normToNumberOfWritings()
 {
     for (int iter = 0 ; iter < nx*ny*nz;iter++)
@@ -1396,6 +1637,27 @@ PointInterpolator::~PointInterpolator() {
     delete[] data;
 }
 
+void PointInterpolator::updatePointInterpolator(Interpolate_on_structured_grid *intermediateGrid)
+{
+
+    nx = intermediateGrid->nx;
+    ny = intermediateGrid->ny;
+    nz = intermediateGrid->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();
+}
+
 
 
 void Interpolate_direct::init()

program/source/interpol/interpol.h

@@ -77,6 +77,8 @@ class Interpolate_on_structured_grid {
   template <class DTyp>
     void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
 
+    void update_Interpolate_on_structured_grid(Blockgrid& blockgrid_);
+
     int getNx(){return nx;}
     int getNy(){return ny;}
     int getNz(){return nz;}
@@ -341,9 +343,12 @@ class PointInterpolator {
    * @param coord_x, coord_y, coord_z, Koordinaten des Punktes
    * @return interpolierter Wert
    **/
+
+  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);
@@ -353,9 +358,9 @@ class PointInterpolator {
 
   void QPrint_VTK(QString DateiName);
 
-  double shiftx, shifty, shiftz;
-  double rotationx, rotationy, rotationz;
 private:
+  double shiftx{0}, shifty{0}, shiftz{0};
+  double rotationx{0}, rotationy{0}, rotationz{0};
   int nx, ny, nz;
   double hx, hy, hz;
   Interpolate_on_structured_grid* interpolatorStructured;