diff --git a/program/source/interpol/interpol.cc b/program/source/interpol/interpol.cc
index 4d6108ce3ab1aa8179772a538a61104364945bc6..47dc4acd50c68fcbfa6eb8592c1be98b8f7909b2 100644
--- a/program/source/interpol/interpol.cc
+++ b/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()
diff --git a/program/source/interpol/interpol.h b/program/source/interpol/interpol.h
index 0dd1f17009fc3e0b6e64d40a31b5acb36edb6cdb..18cd10ce5065c199270f7e01251ae1a2685e5812 100644
--- a/program/source/interpol/interpol.h
+++ b/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;