diff --git a/program/source/interpol/interpol.cc b/program/source/interpol/interpol.cc
index 89c65fb8d76549647e84f3050775dc8365f058d7..556a186faa7b0c0e1be924343f201f2f6cc1193a 100644
--- a/program/source/interpol/interpol.cc
+++ b/program/source/interpol/interpol.cc
@@ -30,32 +30,10 @@
#include "../extemp/extemp.h"
#include "../extemp/parallel.h"
#include "../extemp/variable.h"
-
-
-
-#include "../mympi.h"
-
-#include "../abbrevi.h"
-#include "../parameter.h"
-#include "../math_lib/math_lib.h"
-#include "../math_lib/mg_array.h"
-#include "../basics/basic.h"
-#include "../grid/elements.h"
-#include "../grid/parti.h"
-#include "../grid/ug.h"
-#include "../grid/examples_ug.h"
-#include "../grid/blockgrid.h"
-#include "../grid/marker.h"
-#include "../extemp/extemp.h"
-#include "../extemp/parallel.h"
-#include "../extemp/variable.h"
+#include "../extemp/cellvar.h"
#include "../extemp/co_fu.h"
#include "../extemp/functor.h"
-
-
#include "interpol.h"
-#include <sstream>
-#include <iomanip>
#include "assert.h"
@@ -72,73 +50,40 @@ bool contained_in_tet(D3vector lam) {
return true;
}
-bool contained_in_tet_strong(D3vector lam) {
- if(lam.x < -0.1) return false;
- if(lam.y < -0.1) return false;
- if(lam.z < -0.1) return false;
- if(lam.x + lam.y + lam.z > 1.1) return false;
- return true;
-}
-
-bool new_lam_better(D3vector lam_old, D3vector lam_new) {
-
- if (MIN(lam_new ) < -0.1 || MAX(lam_new) > 1.1) return false;
- if (MIN(lam_old ) < -0.1 || MAX(lam_old) > 1.1) return true;
-
- if (MIN(lam_new) < MIN(lam_old)) return true;
- if (MAX(lam_new) < MAX(lam_old)) return true;
-}
-
bool new_lam_worse(D3vector lam_old, D3vector lam_new) {
if(MIN(lam_new) < MIN(lam_old) && MIN(lam_old) < -0.1) return true;
if(MAX(lam_new) > MAX(lam_old) && MAX(lam_old) > 1.1) return true;
return false;
}
-/*
+Interpolate_on_structured_grid::~Interpolate_on_structured_grid() {
+ delete[] ids_hex;
+ delete[] ids_i;
+ delete[] ids_j;
+ delete[] ids_k;
-Intermadiate_grid_for_PointInterpolator::Intermadiate_grid_for_PointInterpolator(int nx_, int ny_, int nz_, Variable<double>* U_from)
-{
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- if(nx<=2) nx = 3;
- if(ny<=2) ny = 3;
- if(nz<=2) nz = 3;
-
- Blockgrid* blockgrid_from = U_from->Give_blockgrid();
-
- //Variable<double> coordXYZ(*blockgrid);
- X_coordinate Xc(*blockgrid_from);
- Y_coordinate Yc(*blockgrid_from);
- Z_coordinate Zc(*blockgrid_from);
- pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
- pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
+ delete[] typ_tet;
-
- interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
-
+ delete[] lambda;
}
-*/
Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
D3vector pWSD, D3vector pENT,
Blockgrid& blockgrid_) {
int Nx, Ny, Nz;
- int typ;
+ // int typ;
assert(nx_ > 1);
assert(ny_ > 1);
assert(nz_ > 1);
- int ilmin, jlmin, klmin;
- int ilmax, jlmax, klmax;
+ //int ilmin, jlmin, klmin;
+ //int ilmax, jlmax, klmax;
double factor = 0.1;
// double factor = 0.00001;
- D3vector lam;
+ // D3vector lam;
blockgrid = &blockgrid_;
ug = blockgrid->Give_unstructured_grid();
@@ -162,15 +107,17 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
int num_total = nx * ny * nz;
+ /*
D3vector cWSD, cESD;
D3vector cWND, cEND;
D3vector cWST, cEST;
D3vector cWNT, cENT;
+*/
+
+ // D3vector boxWSD, boxENT;
- D3vector boxWSD, boxENT;
-
- D3vector ploc;
+ // D3vector ploc;
ids_hex = new int[num_total];
@@ -188,22 +135,24 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
Nx = blockgrid->Give_Nx_hexahedron(id_hex);
Ny = blockgrid->Give_Ny_hexahedron(id_hex);
Nz = blockgrid->Give_Nz_hexahedron(id_hex);
-
+
+#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
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 );
+ D3vector cWSD = blockgrid->Give_coord_hexahedron(id_hex,i, j, k );
+ D3vector cESD = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k );
+ D3vector cWND = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k );
+ D3vector 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);
+ D3vector cWST = blockgrid->Give_coord_hexahedron(id_hex,i, j, k+1);
+ D3vector cEST = blockgrid->Give_coord_hexahedron(id_hex,i+1,j ,k+1);
+ D3vector cWNT = blockgrid->Give_coord_hexahedron(id_hex,i, j+1,k+1);
+ D3vector cENT = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k+1);
// bounding box calculation
+ D3vector boxWSD, boxENT;
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)),
@@ -219,14 +168,14 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
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);
+ int ilmin = Ganzzahliger_Anteil((boxWSD.x - pWSD.x) / hx);
+ int jlmin = Ganzzahliger_Anteil((boxWSD.y - pWSD.y) / hy);
+ int 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);
+ int ilmax = Ganzzahliger_Anteil((boxENT.x - pWSD.x) / hx);
+ int jlmax = Ganzzahliger_Anteil((boxENT.y - pWSD.y) / hy);
+ int klmax = Ganzzahliger_Anteil((boxENT.z - pWSD.z) / hz);
/*
cout << " indices: "
@@ -269,13 +218,13 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
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;
+ D3vector ploc = D3vector(il * hx, jl * hy, kl * hz) + pWSD;
// cout << "HI" << endl;
- typ = -1;
+ int typ = -1;
- lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST);
+ D3vector 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);
@@ -316,6 +265,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
stop=new_lam_worse(lambda[ind_global],lam);
}
+ #pragma omp critical
if(stop==false) {
ids_hex[ind_global] = id_hex;
ids_i[ind_global] = i;
@@ -373,2015 +323,126 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
}
}
-Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
- Blockgrid& blockgrid_) {
- int Nx, Ny, Nz;
- int typ;
-
-
-
-
- assert(nx_ > 1);
- assert(ny_ > 1);
- assert(nz_ > 1);
-
- int ilmin, jlmin, klmin;
- int ilmax, jlmax, klmax;
-
- double factor = 0.1;
- // double factor = 0.00001;
+/////////////////////////////////////////////////////////////
+// 2. Interpolate from blockgrid to blockgrid
+/////////////////////////////////////////////////////////////
- D3vector lam;
+
+Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
+ Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_) {
+ nx = nx_;
+ ny = ny_;
+ nz = nz_;
+
+ if(nx<=2) nx = 3;
+ if(ny<=2) ny = 3;
+ if(nz<=2) nz = 3;
+ blockgrid_to = blockgrid_to_;
//Variable<double> coordXYZ(*blockgrid);
- X_coordinate Xc(blockgrid_);
- Y_coordinate Yc(blockgrid_);
- Z_coordinate Zc(blockgrid_);
- //D3vector pWSD, pENT;
+ X_coordinate Xc(*blockgrid_to);
+ Y_coordinate Yc(*blockgrid_to);
+ Z_coordinate Zc(*blockgrid_to);
+
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();
-
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- 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;
+ interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
+ data = new double[nx*ny*nz];
+
+
+ hx = (pENT.x - pWSD.x) / (nx-1); // geandert nach Rall-Vorschlag von nx+1 zu nx-1
+ hy = (pENT.y - pWSD.y) / (ny-1);
+ hz = (pENT.z - pWSD.z) / (nz-1);
+
+
+ /*
+
+ // test GGGG
+ cout << "\n WSD: " ; pWSD.Print();
+ cout << "\n ENT: " ; pENT.Print();
+ cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
+ */
+}
+
- 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;
+
+void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<double>* U_to,
+ double defaultInterpolation) {
+/*
+ //test GGGG
+ X_coordinate Xfrom(*U_from->Give_blockgrid());
+ (*U_from) = Xfrom;
+*/
+ interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation);
- // 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);
-
+ /*
+ //test GGGG
+ for(int i=0;i<Nx;++i) for(int j=0;j<nz;++j) for(int k=0;k<nz;++k)
+ data[i +nx*(j +ny* k)] = hx * i;
+ */
+
+ Functor3<double,double,Interpolate_on_block_grid> myFunctor(this);
+
+ X_coordinate Xc(*blockgrid_to);
+ Y_coordinate Yc(*blockgrid_to);
+ Z_coordinate Zc(*blockgrid_to);
+
+ (*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;
+ if(coord_x < pWSD.x) return 0.0;
+ if(coord_y > pENT.y) return 0.0;
+ if(coord_y < pWSD.y) return 0.0;
+ if(coord_z > pENT.z) return 0.0;
+ if(coord_z < pWSD.z) return 0.0;
+
+ int i = (coord_x - pWSD.x) / hx;
+ int j = (coord_y - pWSD.y) / hy;
+ int k = (coord_z - pWSD.z) / hz;
+
+ if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
+ 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)];
+ double uEND = data[(i+1)+nx*((j+1)+ny* k)];
+ double uWST = data[i +nx*(j +ny*(k+1))];
+ double uEST = data[(i+1)+nx*(j +ny*(k+1))];
+ double uWNT = data[i +nx*((j+1)+ny*(k+1))];
+ double uENT = data[(i+1)+nx*((j+1)+ny*(k+1))];
+
+
+ // assert( (i+1)+nx*((j+1)+ny*(k+1)) < nx*ny*nz);
+
+ double locX = (coord_x - pWSD.x) / hx - i;
+ double locY = (coord_y - pWSD.y) / hy - j;
+ double locZ = (coord_z - pWSD.z) / hz - k;
+
+ //return uWSD;
+
+ return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
+ uESD * locX * (1.0 - locY) * (1.0 - locZ) +
+ uWND * (1.0 - locX) * locY * (1.0 - locZ) +
+ uEND * locX * locY * (1.0 - locZ) +
+ uWST * (1.0 - locX) * (1.0 - locY) * locZ +
+ uEST * locX * (1.0 - locY) * locZ +
+ uWNT * (1.0 - locX) * locY * locZ +
+ uENT * locX * locY * locZ;
+}
+
+Interpolate_on_block_grid::~Interpolate_on_block_grid() {
+ delete interpolatorStructured;
+ delete[] data;
+}
- 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;
- }
-
- /*
- cout << " out "
- << " ilmin: " << ilmin
- << " ilmax: " << ilmax
- << " jlmin: " << jlmin
- << " jlmax: " << jlmax
- << " klmin: " << klmin
- << " klmax: " << klmax;
- cout << "\n "; cWSD.Print();
- cout << "\n "; cESD.Print();
- cout << "\n "; cWND.Print();
- cout << "\n "; cEND.Print();
- cout << "\n "; cWST.Print();
- cout << "\n "; cEST.Print();
- cout << "\n "; cWNT.Print();
- cout << "\n "; cENT.Print();
- cout << "\n p: "; ploc.Print();
-
- cout << "\n : "; boxWSD.Print();
- cout << "\n : "; boxENT.Print();
-
- cout << endl;
- */
- }
- }
- }
-
-
- 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;
- }
- }
-}
-
-/////////////////////////////////////////////////////////////
-// 2. Interpolate from blockgrid to blockgrid
-/////////////////////////////////////////////////////////////
-
-
-Interpolate_on_block_grid::Interpolate_on_block_grid(int nx_, int ny_, int nz_,
- Blockgrid* blockgrid_from, Blockgrid* blockgrid_to_) {
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- if(nx<=2) nx = 3;
- if(ny<=2) ny = 3;
- if(nz<=2) nz = 3;
-
- blockgrid_to = blockgrid_to_;
-
- //Variable<double> coordXYZ(*blockgrid);
- X_coordinate Xc(*blockgrid_to);
- Y_coordinate Yc(*blockgrid_to);
- Z_coordinate Zc(*blockgrid_to);
-
- pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
- pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
-
- interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
- data = new double[nx*ny*nz];
-
-
- hx = (pENT.x - pWSD.x) / (nx+1);
- hy = (pENT.y - pWSD.y) / (ny+1); // CHANGES OF RALL MISSING; HAS TO BE FIXED. -> still old version in interpol.cc file used!
- hz = (pENT.z - pWSD.z) / (nz+1);
-
-
- /*
-
- // test GGGG
- cout << "\n WSD: " ; pWSD.Print();
- cout << "\n ENT: " ; pENT.Print();
- cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
- */
-}
-
-
-
-void Interpolate_on_block_grid::interpolate(Variable<double>* U_from, Variable<double>* U_to,
- double defaultInterpolation) {
-/*
- //test GGGG
- X_coordinate Xfrom(*U_from->Give_blockgrid());
- (*U_from) = Xfrom;
-*/
- interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation);
-
- /*
- //test GGGG
- for(int i=0;i<Nx;++i) for(int j=0;j<nz;++j) for(int k=0;k<nz;++k)
- data[i +nx*(j +ny* k)] = hx * i;
- */
-
- Functor3<double,double,Interpolate_on_block_grid> myFunctor(this);
-
- X_coordinate Xc(*blockgrid_to);
- Y_coordinate Yc(*blockgrid_to);
- Z_coordinate Zc(*blockgrid_to);
-
- (*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;
- if(coord_x < pWSD.x) return 0.0;
- if(coord_y > pENT.y) return 0.0;
- if(coord_y < pWSD.y) return 0.0;
- if(coord_z > pENT.z) return 0.0;
- if(coord_z < pWSD.z) return 0.0;
-
- int i = (coord_x - pWSD.x) / hx;
- int j = (coord_y - pWSD.y) / hy;
- int k = (coord_z - pWSD.z) / hz;
-
- if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
- 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)];
- double uEND = data[(i+1)+nx*((j+1)+ny* k)];
- double uWST = data[i +nx*(j +ny*(k+1))];
- double uEST = data[(i+1)+nx*(j +ny*(k+1))];
- double uWNT = data[i +nx*((j+1)+ny*(k+1))];
- double uENT = data[(i+1)+nx*((j+1)+ny*(k+1))];
-
-
- // assert( (i+1)+nx*((j+1)+ny*(k+1)) < nx*ny*nz);
-
- double locX = (coord_x - pWSD.x) / hx - i;
- double locY = (coord_y - pWSD.y) / hy - j;
- double locZ = (coord_z - pWSD.z) / hz - k;
-
-
- return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
- uESD * locX * (1.0 - locY) * (1.0 - locZ) +
- uWND * (1.0 - locX) * locY * (1.0 - locZ) +
- uEND * locX * locY * (1.0 - locZ) +
- uWST * (1.0 - locX) * (1.0 - locY) * locZ +
- uEST * locX * (1.0 - locY) * locZ +
- uWNT * (1.0 - locX) * locY * locZ +
- uENT * locX * locY * locZ;
-}
-
-Interpolate_on_block_grid::~Interpolate_on_block_grid() {
- delete interpolatorStructured;
- delete[] data;
-}
-
-
-/////////////////////////////////////////////////////////////
-// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid
-/////////////////////////////////////////////////////////////
-
-PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
- Variable<double>* U_from, double defaultInterpolation_) {
- defaultInterpolation = defaultInterpolation_;
- shiftx = 0.0;
- shifty = 0.0;
- shiftz = 0.0;
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- if(nx<=2) nx = 3;
- if(ny<=2) ny = 3;
- if(nz<=2) nz = 3;
-
- Blockgrid* blockgrid_from = U_from->Give_blockgrid();
-
- //Variable<double> coordXYZ(*blockgrid);
- X_coordinate Xc(*blockgrid_from);
- Y_coordinate Yc(*blockgrid_from);
- Z_coordinate Zc(*blockgrid_from);
- pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
- pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
-
-
- interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
-
- data = new double[nx*ny*nz];
-
-
- hx = (pENT.x - pWSD.x) / (nx-1);
- hy = (pENT.y - pWSD.y) / (ny-1);
- hz = (pENT.z - pWSD.z) / (nz-1);
-
- interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation_);
-
-
- /*
-
- // test GGGG
- cout << "\n WSD: " ; pWSD.Print();
- cout << "\n ENT: " ; pENT.Print();
- cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
- */
-}
-
-PointInterpolator::PointInterpolator(int nx_, int ny_, int nz_,
- D3vector pWSD_, D3vector pENT_,
- Variable<double>* U_from, double defaultInterpolation_) {
- defaultInterpolation = defaultInterpolation_;
-
- shiftx = 0.0;
- shifty = 0.0;
- shiftz = 0.0;
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- if(nx<=2) nx = 3;
- if(ny<=2) ny = 3;
- if(nz<=2) nz = 3;
-
- Blockgrid* blockgrid_from = U_from->Give_blockgrid();
-
- pWSD = pWSD_;
- pENT = pENT_;
-
- interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
- data = new double[nx*ny*nz];
-
-
- hx = (pENT.x - pWSD.x) / (nx-1);
- hy = (pENT.y - pWSD.y) / (ny-1);
- hz = (pENT.z - pWSD.z) / (nz-1);
-
-
- interpolatorStructured->interpolate<double>(*U_from,data,defaultInterpolation);
-
-
-
-
- /*
-
- // test GGGG
- cout << "\n WSD: " ; pWSD.Print();
- cout << "\n ENT: " ; pENT.Print();
- cout << "nx: " << nx << " ny: " << ny << " nz: " << nz << endl;
- */
-}
-
-PointInterpolator::PointInterpolator(Interpolate_on_structured_grid* intermediateGrid, Variable<double>* U_from, double defaultInterpolation_)
-{
- defaultInterpolation = defaultInterpolation_;
-
- nx = intermediateGrid->nx;
- ny = intermediateGrid->ny;
- nz = intermediateGrid->nz;
-
- data = new double[nx*ny*nz];
-
- pENT = intermediateGrid->pENT;
- pWSD = intermediateGrid->pWSD;
-
- shiftx = 0.0;
- shifty = 0.0;
- shiftz = 0.0;
- hx = (pENT.x - pWSD.x) / (nx-1);
- hy = (pENT.y - pWSD.y) / (ny-1);
- hz = (pENT.z - pWSD.z) / (nz-1);
-
- intermediateGrid->interpolate<double>(*U_from,data,defaultInterpolation_);
-}
-
- /*
-Interpolate_on_structured_grid* PointInterpolator::intermediateGrid(int nx_, int ny_, int nz_, Variable<double>* U_from)
-{
- nx = nx_;
- ny = ny_;
- nz = nz_;
-
- if(nx<=2) nx = 3;
- if(ny<=2) ny = 3;
- if(nz<=2) nz = 3;
-
- Blockgrid* blockgrid_from = U_from->Give_blockgrid();
-
- //Variable<double> coordXYZ(*blockgrid);
- X_coordinate Xc(*blockgrid_from);
- Y_coordinate Yc(*blockgrid_from);
- Z_coordinate Zc(*blockgrid_from);
- pWSD.x = Minimum(Xc); pWSD.y = Minimum(Yc); pWSD.z = Minimum(Zc);
- pENT.x = Maximum(Xc); pENT.y = Maximum(Yc); pENT.z = Maximum(Zc);
-
-
- interpolatorStructured = new Interpolate_on_structured_grid(nx,ny,nz, pWSD, pENT, *blockgrid_from);
-
- return interpolatorStructured;
-}
-*/
-
-
-
-double PointInterpolator::evaluate(double coord_x, double coord_y, double coord_z) {
-
- coord_x-=shiftx;
- coord_y-=shifty;
- coord_z-=shiftz;
-
- if(coord_x > pENT.x) return defaultInterpolation;
- if(coord_x < pWSD.x) return defaultInterpolation;
- if(coord_y > pENT.y) return defaultInterpolation;
- if(coord_y < pWSD.y) return defaultInterpolation;
- if(coord_z > pENT.z) return defaultInterpolation;
- if(coord_z < pWSD.z) return defaultInterpolation;
-
- //cout << "coord_z " << coord_z << " pWSD.z " << pWSD.z << endl;
- /*
- int i = (coord_x - pWSD.x) / hx;
- int j = (coord_y - pWSD.y) / hy;
- int k = (coord_z - pWSD.z) / hz;
- */
- double id = (coord_x - pWSD.x) / hx;
- double jd = (coord_y - pWSD.y) / hy;
- double kd = (coord_z - pWSD.z) / hz;
-
-
- int i = int(id);
- int j = int(jd);
- int k = int(kd);
-
- if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
- if(i>=nx-1) i=nx-2; if(j>=ny-1) j=ny-2; if(k>=nz-1) k=nz-2;
-
-
-
- //cout << "hx " << hx << " hy "<< hy << " hz " << hz << endl;
- //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)];
- double uESD = data[(i+1)+nx*(j +ny* k)];
- double uWND = data[i +nx*((j+1)+ny* k)];
- double uEND = data[(i+1)+nx*((j+1)+ny* k)];
- double uWST = data[i +nx*(j +ny*(k+1))];
- double uEST = data[(i+1)+nx*(j +ny*(k+1))];
- double uWNT = data[i +nx*((j+1)+ny*(k+1))];
- double uENT = data[(i+1)+nx*((j+1)+ny*(k+1))];
- //i++;
- //j++;
- //k++;
- //k++;
- //cout << "uWSD " << uWSD << endl;
- //cout << "uESD " << uESD << endl;
- //cout << "uWND " << uWND << endl;
- //cout << "uEND " << uEND << endl;
- //cout << "uWST " << uWST << endl;
- //cout << "uEST " << uEST << endl;
- //cout << "uWNT " << uWNT << endl;
- //cout << "uENT " << uENT << endl;
-
-
- //cout << "x+1 "<< data[(i+2)+nx*(j +ny* k)] << endl;
- //cout << "x-1 " <<data[i-1 +nx*(j +ny* k)] << endl;
- //cout << "x-1, y-1 " <<data[i-1 +nx*(j-1 +ny* k)] << endl;
-
- // assert( (i+1)+nx*((j+1)+ny*(k+1)) < nx*ny*nz);
-
- double posX = (coord_x - pWSD.x) ;
- double locX = posX / hx - i;
- double posY = (coord_y - pWSD.y);
- double locY = posY / hy - j;
- double posZ = (coord_z - pWSD.z);
- double locZ = posZ / hz - k;
-
-
-
- //cout << "locX, Y, Z: " << locX << " " << locY << " " << locZ << endl;
- //return uWSD;
-
-
- //cout << "uPOS : " << uWSD << " , " << uESD << " , " << uWND << " , " << uEND << " , " << uWST << " , " << uEST << " , " << uWNT << " , " << uENT << endl;
- double uTOT(0);
- double uET, uWT, uWD, uED;
- double uT, uD;
-
- if ( (uEST != defaultInterpolation) == (uENT != defaultInterpolation) ) { uET = uEST * (1.0 - locY) + uENT * locY ;}
- else if ( (uEST != defaultInterpolation) && (uENT == defaultInterpolation) ) { uET = uEST;}
- else { uET = uENT;}
-
- if ( (uWST != defaultInterpolation) == (uWNT != defaultInterpolation) ) {uWT = uWST * (1.0 - locY) + uWNT * locY ;}
- else if ( (uWST != defaultInterpolation) && (uWNT == defaultInterpolation) ) {uWT = uWST;}
- else {uWT = uWNT;}
-
- if ( (uESD != defaultInterpolation) == (uEND != defaultInterpolation) ) {uED = uESD * (1.0 - locY) + uEND * locY ;}
- else if ( (uESD != defaultInterpolation) && (uEND == defaultInterpolation) ) {uED = uESD;}
- else {uED = uEND;}
-
- if ( (uWSD != defaultInterpolation) == (uWND != defaultInterpolation) ) {uWD = uWSD * (1.0 - locY) + uWND * locY ;}
- else if ( (uWSD != defaultInterpolation) && (uWND == defaultInterpolation) ) {uWD = uWSD;}
- else {uWD = uWND;}
-
- if ( (uET != defaultInterpolation) == (uWT != defaultInterpolation) ) {uT = uWT * (1.0 - locX) + uET * locX ;}
- else if ( (uET != defaultInterpolation) && (uWT == defaultInterpolation) ) {uT = uET;}
- else {uT = uWT;}
-
- if ( (uED != defaultInterpolation) == (uWD != defaultInterpolation) ) {uD = uWD * (1.0 - locX) + uED * locX ;}
- else if ( (uED != defaultInterpolation) && (uWD == defaultInterpolation) ) {uD = uED;}
- else {uD = uWD;}
-
- if ( (uT != defaultInterpolation) == (uD != defaultInterpolation) ) {uTOT = uD * (1.0 - locZ) + uT * locZ ;}
- else if ( (uT != defaultInterpolation) && (uD == defaultInterpolation) ) {uTOT = uT;}
- else {uTOT = uD;}
-
-
-
-// if (uWSD != defaultInterpolation)
-// { uTOT += uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ); }
-// if (uESD != defaultInterpolation)
-// { uTOT += uESD * locX * (1.0 - locY) * (1.0 - locZ); }
-// if (uWND != defaultInterpolation)
-// { uTOT += uWND * (1.0 - locX) * locY * (1.0 - locZ); }
-// if (uEND != defaultInterpolation)
-// { uTOT += uEND * locX * locY * (1.0 - locZ); }
-// if (uWST != defaultInterpolation)
-// { uTOT += uWST * (1.0 - locX) * (1.0 - locY) * locZ; }
-// if (uEST != defaultInterpolation)
-// { uTOT += uEST * locX * (1.0 - locY) * locZ; }
-// if (uWNT != defaultInterpolation)
-// { uTOT += uWNT * (1.0 - locX) * locY * locZ; }
-// if (uENT != defaultInterpolation)
-// { uTOT += uENT * locX * locY * locZ; }
-
- //cout << "my method, other method " << uTOT << " , " << uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
- // uESD * locX * (1.0 - locY) * (1.0 - locZ) +
- // uWND * (1.0 - locX) * locY * (1.0 - locZ) +
- // uEND * locX * locY * (1.0 - locZ) +
- // uWST * (1.0 - locX) * (1.0 - locY) * locZ +
- // uEST * locX * (1.0 - locY) * locZ +
- // uWNT * (1.0 - locX) * locY * locZ +
- // uENT * locX * locY * locZ << endl;
-
-
- //cout <<endl<< "RESULT: " << uTOT<<endl<<endl;
- return uTOT;
-
-
-
-
- /*
- return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
- uESD * locX * (1.0 - locY) * (1.0 - locZ) +
- uWND * (1.0 - locX) * locY * (1.0 - locZ) +
- uEND * locX * locY * (1.0 - locZ) +
- uWST * (1.0 - locX) * (1.0 - locY) * locZ +
- uEST * locX * (1.0 - locY) * locZ +
- uWNT * (1.0 - locX) * locY * locZ +
- uENT * locX * locY * locZ;
- */
-}
-
-std::vector<double> PointInterpolator::evaluateGradient(double coord_x, double coord_y, double coord_z){
- coord_x-=shiftx;
- coord_y-=shifty;
- coord_z-=shiftz;
- //cout << "coord_x y z " << coord_x << " " << coord_y << " " << coord_z << endl;
- if(coord_x > pENT.x) return std::vector<double>({0.0 , 0.0 , 0.0});
- if(coord_x < pWSD.x) return std::vector<double>({0.0 , 0.0 , 0.0});
- if(coord_y > pENT.y) return std::vector<double>({0.0 , 0.0 , 0.0});
- if(coord_y < pWSD.y) return std::vector<double>({0.0 , 0.0 , 0.0});
- if(coord_z > pENT.z) return std::vector<double>({0.0 , 0.0 , 0.0});
- if(coord_z < pWSD.z) return std::vector<double>({0.0 , 0.0 , 0.0});
-
-
- double id = (coord_x - pWSD.x) / hx;
- double jd = (coord_y - pWSD.y) / hy;
- double kd = (coord_z - pWSD.z) / hz;
-
-
- int i = int(id);
- int j = int(jd);
- int k = int(kd);
-
- if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
- if(i>=nx-1) i=nx-2; if(j>=ny-1) j=ny-2; if(k>=nz-1) k=nz-2;
-
-
- //cout << "coord_i j k " << i << " " << j << " " << 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)];
- double uEND = data[(i+1)+nx*((j+1)+ny* k)];
- double uWST = data[i +nx*(j +ny*(k+1))];
- double uEST = data[(i+1)+nx*(j +ny*(k+1))];
- double uWNT = data[i +nx*((j+1)+ny*(k+1))];
- double uENT = data[(i+1)+nx*((j+1)+ny*(k+1))];
-
- double posX = (coord_x - pWSD.x) ;
- double locX = posX / hx - i;
- double posY = (coord_y - pWSD.y);
- double locY = posY / hy - j;
- double posZ = (coord_z - pWSD.z);
- double locZ = posZ / hz - k;
-
-
-
- //cout << "coord_locX locY locZ " << locX << " " << locY << " " << locZ << endl;
-
- if (uWSD == 0.0 || uESD == 0.0 || uWND == 0.0 || uEND == 0.0 || uWST == 0.0 || uEST == 0.0 || uWNT == 0.0 || uENT == 0.0)
- {
- std::vector<double> gradient = { 0.0 , 0.0 , 0.0};
- return gradient;
- }
- double uTOT(0);
- double uET, uWT, uWD, uED;
- double uT, uD, uN, uS, uW, uE;
-
- double uGradientZ, uGradientX, uGradientY;
-
- //assume: all values are != defaultInterpolation
- //does not hold for curved interfaces
- /*
- return uWSD * (1.0 - locX) * (1.0 - locY) * (1.0 - locZ) +
- uESD * locX * (1.0 - locY) * (1.0 - locZ) +
- uWND * (1.0 - locX) * locY * (1.0 - locZ) +
- uEND * locX * locY * (1.0 - locZ) +
- uWST * (1.0 - locX) * (1.0 - locY) * locZ +
- uEST * locX * (1.0 - locY) * locZ +
- uWNT * (1.0 - locX) * locY * locZ +
- uENT * locX * locY * locZ;
- */
- uGradientX = uWSD * (0.0 - 1.0) * (1.0 - locY) * (1.0 - locZ) +
- uESD * 1.0 * (1.0 - locY) * (1.0 - locZ) +
- uWND * (0.0 - 1.0) * locY * (1.0 - locZ) +
- uEND * 1.0 * locY * (1.0 - locZ) +
- uWST * (0.0 - 1.0) * (1.0 - locY) * locZ +
- uEST * 1.0 * (1.0 - locY) * locZ +
- uWNT * (0.0 - 1.0) * locY * locZ +
- uENT * 1.0 * locY * locZ;
- uGradientY = uWSD * (1.0 - locX) * (0.0 - 1.0) * (1.0 - locZ) +
- uESD * locX * (0.0 - 1.0) * (1.0 - locZ) +
- uWND * (1.0 - locX) * 1.0 * (1.0 - locZ) +
- uEND * locX * 1.0 * (1.0 - locZ) +
- uWST * (1.0 - locX) * (0.0 - 1.0) * locZ +
- uEST * locX * (0.0 - 1.0) * locZ +
- uWNT * (1.0 - locX) * 1.0 * locZ +
- uENT * locX * 1.0 * locZ;
-
- uGradientZ = uWSD * (1.0 - locX) * (1.0 - locY) * (0.0 - 1.0) +
- uESD * locX * (1.0 - locY) * (0.0 - 1.0) +
- uWND * (1.0 - locX) * locY * (0.0 - 1.0) +
- uEND * locX * locY * (0.0 - 1.0) +
- uWST * (1.0 - locX) * (1.0 - locY) * 1.0 +
- uEST * locX * (1.0 - locY) * 1.0 +
- uWNT * (1.0 - locX) * locY * 1.0 +
- uENT * locX * locY * 1.0;
-
- std::vector<double> gradient = {uGradientX / hx, uGradientY / hy, uGradientZ / hz};
-
- return gradient;
-
-}
-
-void PointInterpolator::smoothGrid()
-{
- for (int i = 1 ; i < nx -1 ; i++)
- {
- for (int j = 1 ; j < ny -1 ; j++)
- {
- for (int k = 1 ; k < nz -1 ; k++)
- {
- data[i +nx*(j +ny* k)] = (2.0 * data[i +nx*(j +ny* k)]
- + data[i+1 +nx*(j +ny* k)]
- + data[i-1 +nx*(j +ny* k)]
- + data[i+1 +nx*(j+1 +ny* k)]
- + data[i+1 +nx*(j-1 +ny* k)]
- + data[i+1 +nx*(j +ny* k+1)]
- + data[i+1 +nx*(j +ny* k-1)]) / 8.0;
-
- }
- }
- }
-}
-
-void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value)
-{
- coord_x-=shiftx;
- coord_y-=shifty;
- coord_z-=shiftz;
-
- if(coord_x > pENT.x) return;
- if(coord_x < pWSD.x) return;
- if(coord_y > pENT.y) return;
- if(coord_y < pWSD.y) return;
- if(coord_z > pENT.z) return;
- if(coord_z < pWSD.z) return;
-
- double id = (coord_x - pWSD.x) / hx;
- double jd = (coord_y - pWSD.y) / hy;
- double kd = (coord_z - pWSD.z) / hz;
-
-
- int i = int(id);
- int j = int(jd);
- int k = int(kd);
-
- if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
- if(i>=nx-1) i=nx-2; if(j>=ny-1) j=ny-2; if(k>=nz-1) k=nz-2;
-
- double posX = (coord_x - pWSD.x) ;
- double locX = posX / hx - i;
- double posY = (coord_y - pWSD.y);
- double locY = posY / hy - j;
- double posZ = (coord_z - pWSD.z);
- double locZ = posZ / hz - k;
-
-
-
- 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;
-
- return;
-
-
-}
-
-void PointInterpolator::subtractOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value)
-{
- coord_x+=shiftx;
- coord_y+=shifty;
- coord_z+=shiftz;
- if(coord_x > pENT.x) return;
- if(coord_x < pWSD.x) return;
- if(coord_y > pENT.y) return;
- if(coord_y < pWSD.y) return;
- if(coord_z > pENT.z) return;
- if(coord_z < pWSD.z) return;
-
- double id = (coord_x - pWSD.x) / hx;
- double jd = (coord_y - pWSD.y) / hy;
- double kd = (coord_z - pWSD.z) / hz;
-
-
- int i = int(id);
- int j = int(jd);
- int k = int(kd);
-
- if(i < 0) i=0; if(j <0 ) j=0; if(k<0) k=0;
- if(i>=nx-1) i=nx-2; if(j>=ny-1) j=ny-2; if(k>=nz-1) k=nz-2;
-
- double posX = (coord_x - pWSD.x) ;
- double locX = posX / hx - i;
- double posY = (coord_y - pWSD.y);
- double locY = posY / hy - j;
- double posZ = (coord_z - pWSD.z);
- double locZ = posZ / hz - k;
-
-
-
- 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;
-
- return;
-}
-
-void PointInterpolator::shiftInterpolatedGrid(double shift_x, double shift_y, double shift_z)
-{
- shiftx = shift_x;
- shifty = shift_y;
- shiftz = shift_z;
-}
-
-void PointInterpolator::scaleInterpolatedData(double scale, double zeroVal)
-{
- if (scale == 1.0)
- {
- return;
- }
- for (int k = 0; k < nz; k++)
- {
- for (int j = 0; j < ny; j++)
- {
- for (int i = 0; i < nx; i++)
- {
- if (data[i +nx*(j +ny* k)] != defaultInterpolation )
- {
- data[i +nx*(j +ny* k)] = (data[i +nx*(j +ny* k)] - zeroVal) * scale + zeroVal;
- if (data[i +nx*(j +ny* k)] <= 1.0 & zeroVal != 0.0)
- {
- cout << "Warning: " << data[i +nx*(j +ny* k)] << endl;
- }
- }
- }
- }
- }
-}
-
-
-PointInterpolator::~PointInterpolator() {
- delete interpolatorStructured;
- delete[] data;
-}
-
-
-
-void Interpolate_direct::init()
-{
- arrayBoxWSDENT.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
- array_box_boundary.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
-
- int counter = 0;
- int counterTwoNeighbours = 0;
- int counterBoxesAtBoundary = 0;
-
- for (int id_hex = 0 ; id_hex < blockgrid->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++)
- {
- int Nx = blockgrid->Give_Nx_hexahedron(id_hex);
- int Ny = blockgrid->Give_Ny_hexahedron(id_hex);
- int Nz = blockgrid->Give_Nz_hexahedron(id_hex);
- arrayBoxWSDENT.at(id_hex).resize(Nx*Ny*Nz);
- array_box_boundary.at(id_hex).resize(Nx*Ny*Nz);
- }
- std::vector<std::vector<std::vector<int> > >addToIndex;
- addToIndex.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
-
-
- for (int id_hex = 0 ; id_hex < blockgrid->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++)
- {
- int Nx = blockgrid->Give_Nx_hexahedron(id_hex); // now without +1, since for ( n gridpoints, only n-1 surrounding blocks exist!)
- int Ny = blockgrid->Give_Ny_hexahedron(id_hex);
- int Nz = blockgrid->Give_Nz_hexahedron(id_hex);
- addToIndex.at(id_hex).resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
-
- //if (!blockgrid->blockgrid_hexa_boundary.at(id_hex).at(0).empty() || true)
- //{
-
- for (int i = 0 ; i < Nx+1; i = i + Nx) for (int j = 0 ; j < Ny+1; j = j + Ny) for (int k = 0 ; k < Nz+1; k = k + Nz)
- {
- for (int iterIndexToAdd = 0 ; iterIndexToAdd <(int)blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).size() ; iterIndexToAdd+=4 )
- {
- //cout << "i "<< i << " j "<< j << " k "<< k << endl;
-
- int idHexOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+0);
- int iOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+1);
- int jOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+2);
- int kOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+3);
-
- int iAdd = (iOpp==i) ? (0) : (-1);
- int jAdd = (jOpp==j) ? (0) : (-1);
- int kAdd = (kOpp==k) ? (0) : (-1);
-
- //cout << "convert: = " <<idHexOpp<< " , iAdd = " << iAdd << " , jAdd = " << jAdd << " , kAdd = " << kAdd << endl;
- //blockgrid->blockgrid_hexa_coordinates.at(idHexOpp).at(iOpp +(Nx+1)*(jOpp +(Ny+1)* kOpp)).Print();cout<<endl;
- //blockgrid->blockgrid_hexa_coordinates.at(id_hex ).at(i +(Nx+1)*(j +(Ny+1)* k)).Print();cout<<endl;
-
-
- if (addToIndex.at(id_hex).at(idHexOpp).empty())
- {
- //cout << "added!" << endl<< endl;
- addToIndex.at(id_hex).at(idHexOpp).resize(4);
- addToIndex.at(id_hex).at(idHexOpp).at(0) = iAdd;
- addToIndex.at(id_hex).at(idHexOpp).at(1) = jAdd;
- addToIndex.at(id_hex).at(idHexOpp).at(2) = kAdd;
- addToIndex.at(id_hex).at(idHexOpp).at(3) = 0;
-
- }
- else
- {
- //check is same:
- if (addToIndex.at(id_hex).at(idHexOpp).at(0) == iAdd && addToIndex.at(id_hex).at(idHexOpp).at(1) == jAdd && addToIndex.at(id_hex).at(idHexOpp).at(2) == kAdd)
- {
- //cout << "ok!" << endl<< endl;
- }
- else
- {
-
- //add : invert ij, jk, ik
-
- addToIndex.at(id_hex).at(idHexOpp).at(3) = 1;
-
- }
- }
-
- }
-
- }
- //check again, maybe invert changes the need for invert everything
- for (int i = 0 ; i < Nx+1; i = i + Nx) for (int j = 0 ; j < Ny+1; j = j + Ny) for (int k = 0 ; k < Nz+1; k = k + Nz)
- {
- for (int iterIndexToAdd = 0 ; iterIndexToAdd <(int)blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).size() ; iterIndexToAdd+=4 )
- {
- int idHexOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+0);
- int iOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+1);
- int jOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+2);
- int kOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+3);
-
- int iAdd = (iOpp==j) ? (0) : (-1); // invert here already!
- int jAdd = (jOpp==i) ? (0) : (-1);
- int kAdd = (kOpp==k) ? (0) : (-1);
-
- // compare to previous :
- if (addToIndex.at(id_hex).at(idHexOpp).at(3) == 1)
- {
- if (addToIndex.at(id_hex).at(idHexOpp).at(0) == jAdd && addToIndex.at(id_hex).at(idHexOpp).at(1) == iAdd )
- {
-
- }
- else
- {
-
- addToIndex.at(id_hex).at(idHexOpp).at(0) = (addToIndex.at(id_hex).at(idHexOpp).at(0) == 0) ? -1 : 0 ;
- addToIndex.at(id_hex).at(idHexOpp).at(1) = (addToIndex.at(id_hex).at(idHexOpp).at(1) == 0) ? -1 : 0 ;
-
- }
- }
-
-
-
- }
-
- }
-
- for (int i = 0 ; i < Nx+1; i = i + Nx) for (int j = 0 ; j < Ny+1; j = j + Ny) for (int k = 0 ; k < Nz+1; k = k + Nz)
- {
- for (int iterIndexToAdd = 0 ; iterIndexToAdd <(int)blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).size() ; iterIndexToAdd+=4 )
- {
- int idHexOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+0);
- int iOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+1);
- int jOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+2);
- int kOpp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(i +(Nx+1)*(j +(Ny+1)* k)).at(iterIndexToAdd+3);
-
- int iAdd = (iOpp==j) ? (0) : (-1); // invert here already!
- int jAdd = (jOpp==i) ? (0) : (-1);
- int kAdd = (kOpp==k) ? (0) : (-1);
-
- if (addToIndex.at(id_hex).at(idHexOpp).at(3) == 1)
- {
- if (addToIndex.at(id_hex).at(idHexOpp).at(0) == jAdd && addToIndex.at(id_hex).at(idHexOpp).at(1) == iAdd )
- {
-
- }
- else
- {
- addToIndex.at(id_hex).at(idHexOpp).at(0) = (addToIndex.at(id_hex).at(idHexOpp).at(0) == 0) ? -1 : 0 ;
- addToIndex.at(id_hex).at(idHexOpp).at(1) = (addToIndex.at(id_hex).at(idHexOpp).at(1) == 0) ? -1 : 0 ;
-
- }
- }
-
-
-
- }
-
- }
-
-
-
-
-
-//cout << "new try "<<endl;
- for(int i=0;i<Nx;++i) for(int j=0;j<Nz;++j) for(int k=0;k<Nz;++k)
- {
- arrayBoxWSDENT.at(id_hex).at(i +Nx*(j +Ny* k)).resize(2);
- D3vector cWSD, cESD;
- D3vector cWND, cEND;
-
- D3vector cWST, cEST;
- D3vector cWNT, cENT;
-
- D3vector boxWSD, boxENT;
- 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;
- arrayBoxWSDENT.at(id_hex).at(i +Nx*(j +Ny* k)).at(0) = boxWSD;
- arrayBoxWSDENT.at(id_hex).at(i +Nx*(j +Ny* k)).at(1) = boxENT;
-
-
-
- }
- for(int i=0;i<Nx+1;++i) for(int j=0;j<Nz+1;++j) for(int k=0;k<Nz+1;++k)
- {
- if ( ( (i == 0 || i == Nx) || (j == 0 || j == Ny) || (k == 0 || k == Nz) ) )
- {
- int convertedI, convertedJ, convertedK;
- convertedI = i;
- convertedJ = j;
- convertedK = k;
- //cout << "before: i = " << i << " , j = " << j << " , k = " << k << endl;
-
-// if ((i == 0 || i == Nx-1))
-// {convertedI = (i == 0) ? 0 : Nx;}
-// if ((j == 0 || j == Ny-1))
-// {convertedJ = (j == 0) ? 0 : Ny;}
-// if ((k == 0 || k == Nz-1))
-// {convertedK = (k == 0) ? 0 : Nx;}
-
- //cout << "after : i = " << convertedI << " , j = " << convertedJ << " , k = " << convertedK << endl;
- int convertedLocalIndex = convertedI +(Nx+1)*(convertedJ +(Ny+1)* convertedK);
- if (!blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).empty())
- {
- for (int iterHexaTemp = 0 ; iterHexaTemp < (int)blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).size() ; iterHexaTemp +=4)
- {
- convertedI = i;
- convertedJ = j;
- convertedK = k;
-
-
- int idHexTemp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).at(iterHexaTemp+0);
- int iTemp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).at(iterHexaTemp+1);
- int jTemp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).at(iterHexaTemp+2);
- int kTemp = blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).at(iterHexaTemp+3);
-
- int iAdd = addToIndex.at(id_hex).at(idHexTemp).at(0);
- int jAdd = addToIndex.at(id_hex).at(idHexTemp).at(1);
- int kAdd = addToIndex.at(id_hex).at(idHexTemp).at(2);
-
- // cout << "idHex " << idHexTemp << " i " << iTemp << " j " << jTemp << " k " << kTemp <<endl;
- if (addToIndex.at(id_hex).at(idHexTemp).at(2) == 1)
- {
- int tempVar = iTemp;
- iTemp = jTemp;
- jTemp = tempVar;
- iTemp = iTemp+ iAdd;
- jTemp = jTemp+ jAdd;
- kTemp = kTemp+ kAdd;
- }
- else
- {
- iTemp = iTemp+ iAdd;
- jTemp = jTemp+ jAdd;
- kTemp = kTemp+ kAdd;
- }
- //
- //if (iTemp >=0 && jTemp >=0 && kTemp >=0 && iTemp < Nx && jTemp < Ny && kTemp < Nz)
- if (true || iTemp >=0 && jTemp >=0 && kTemp >=0 && i < Nx+1 && j < Ny+1 && k < Nz+1)
- {
-
- if (kTemp != k)
- {
- cout << "why? should not happen " << endl;
- }
-// if (i == 4 && j == 0 && k == 4)
-// {
-// cout << "stop " << endl;
-// }
-
- //cout << endl;
- //cout << "On blockgrid : idHex i j k " << id_hex << " "<< i << " "<< j << " "<< k << endl;
- //cout << "Neighbour of blockgrid: idHex i j k " << idHexTemp << " " << iTemp << " " << jTemp << " " << kTemp <<endl;
- if (addToIndex.at(id_hex).at(idHexTemp).at(3) == 0)
- {
- if (i == Nx)
- {
- convertedI = Nx-1;
- if (iAdd == 0)
- iTemp--;
- else
- iTemp++;
- }
-
- if (j == Ny)
- {
- convertedJ = Ny-1;
- if (jAdd == 0)
- jTemp--;
- else
- jTemp++;
-
- }
- }
- else
- {
- if (j == Nx)
- {
- convertedJ = Nx-1;
- if (jAdd == 0)
- iTemp--;
- else
- iTemp++;
- }
-
- if (i == Ny)
- {
- convertedI = Ny-1;
- if (iAdd == 0)
- jTemp--;
- else
- jTemp++;
-
- }
- }
-
- if (k == Nz)
- {convertedK = Nz-1;
- kTemp--;}
-
-
-
- int indexLokal = convertedI +(Nx-0)*(convertedJ +(Nx-0)* convertedK);
-
- bool addToBoundary = true;
- for (int iterArrayTest = 0 ; iterArrayTest < (int)this->array_box_boundary.at(id_hex).at(indexLokal).size() ; iterArrayTest +=4)
- {
- if (this->array_box_boundary.at(id_hex).at(indexLokal).at(iterArrayTest+0) == idHexTemp)
- {
- addToBoundary = false;
- }
- }
- if (addToBoundary)
- {
- counterBoxesAtBoundary++;
- this->array_box_boundary.at(id_hex).at(indexLokal).push_back(idHexTemp);
- this->array_box_boundary.at(id_hex).at(indexLokal).push_back(iTemp);
- this->array_box_boundary.at(id_hex).at(indexLokal).push_back(jTemp);
- this->array_box_boundary.at(id_hex).at(indexLokal).push_back(kTemp);
- }
-
-
- }
- else
- {
-
- }
-
- }
-
-
- counter += blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedI +(Nx+1)*(convertedJ +(Ny+1)* convertedK)).size() / 4;
- if (blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(id_hex).at(convertedI +(Nx+1)*(convertedJ +(Ny+1)* convertedK)).size()>4)
- counterTwoNeighbours++;
- // blockgrid->blockgrid_hexa_boundary.at(id_hex).at(convertedI +(Nx+1)*(convertedJ +(Ny+1)* convertedK)).clear();
- }
- for (int iterBoundaryBoxes = 0 ; iterBoundaryBoxes < blockgrid->Give_unstructured_grid()->Give_number_hexahedra(); iterBoundaryBoxes++)
- {
- addToIndex.at(id_hex).at(iterBoundaryBoxes);
- addToIndex.at(id_hex).at(iterBoundaryBoxes);
-
-
- }
- }
- }
- }
- for (int idHex = 0 ; idHex < blockgrid->Give_unstructured_grid()->Give_number_hexahedra() ; idHex++)
- {
- int Nx = blockgrid->Give_Nx_hexahedron(idHex);
- int Ny = blockgrid->Give_Ny_hexahedron(idHex);
- int Nz = blockgrid->Give_Nz_hexahedron(idHex);
- for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j) for(int k=0;k<Nz;++k)
- {
-
- for (int iterHexaTemp = 0 ; iterHexaTemp < array_box_boundary.at(idHex).at(i +(Nx)*(j +(Ny)* k)).size() ; iterHexaTemp +=4)
- {
- int index = i +(Nx)*(j +(Ny)* k) ;
- int idHexTemp = array_box_boundary.at(idHex).at(index).at(iterHexaTemp+0);
- int iTemp = array_box_boundary.at(idHex).at(index).at(iterHexaTemp+1);
- int jTemp = array_box_boundary.at(idHex).at(index).at(iterHexaTemp+2);
- int kTemp = array_box_boundary.at(idHex).at(index).at(iterHexaTemp+3);
- if (iterHexaTemp == 4)
- {
- //cout << "4!" << endl;
- }
-
-// D3vector t1 = arrayBoxWSDENT.at(idHexTemp).at(iTemp +(Nx-0)*(jTemp +(Nx-0)* kTemp)).at(0);
-// D3vector t2 = arrayBoxWSDENT.at(idHexTemp).at(iTemp +(Nx-0)*(jTemp +(Nx-0)* kTemp)).at(1);
-// writeBox(t1,t2,std::string("testBox"));
- }
-
- }
- }
-
- //cout << "counter " << counter << endl;
- //cout << "boundary boxes " << counterBoxesAtBoundary<< endl;
- //cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
- //cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
-
-}
-
-void Interpolate_direct::initBoundary()
-{
-
-}
-
-void Interpolate_direct::find_cell(D3vector v)
-{
-// vPrevPrev = vPrev;
-// vPrev = vNow;
- vNow = v;
- lambda = D3vector(-1,-1,-1);
-
-
- D3vector cWSD, cESD;
- D3vector cWND, cEND;
-
- D3vector cWST, cEST;
- D3vector cWNT, cENT;
-
- D3vector boxWSD, boxENT;
-
- //cout << "looking for new cell!" << endl;
-
- int Nx, Ny, Nz;
- // prev
- if (checkBox(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
- counterFastest++;
- setPrevIndex();
-
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- return;
- }
-
- if (checkBoxSurroundingOptimized(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
- counterFast++;
- setPrevIndex();
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- return;}
- // prev prev
- if (checkBox(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1){
- //counterFastest++;
- setPrevIndex();
- setPrevPrevIndex();
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- counterSecondTry++;
- return;}
-
- if (checkBoxSurroundingOptimized(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1){
- //counterFast++;
- setPrevIndex();
- setPrevPrevIndex();
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- counterSecondTry++;
- return;}
- // prev prev prev
- if (checkBox(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1){
- //counterFastest++;
- setPrevIndex();
- setPrevPrevPrevIndex();
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- counterSecondTry++;
- return;}
-
- if (checkBoxSurroundingOptimized(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1){
- //counterFast++;
- setPrevIndex();
- setPrevPrevIndex();
- setPrevPrevPrevIndex();
- //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
- counterSecondTry++;
- return;}
-
-
-
- idHexNow = -1;
-
-
- // cout << "Point looking for :";v.Print();cout<<endl;
-
- for(int id_hex=0;id_hex<blockgrid->Give_unstructured_grid()->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) {
- if (0 < checkBox(id_hex,i, j, k, v))
- {
-// k =Nz;
-// j =Ny;
-// i =Nx;
-// id_hex = blockgrid->Give_unstructured_grid()->Give_number_hexahedra();
- }
-
- }
- }
-// cout << " pp: " << idHexPrevPrevPrev << " " << iPrevPrevPrev << " " << jPrevPrevPrev << " " << kPrevPrevPrev << endl;
-// cout << " pp: " << idHexPrevPrev << " " << iPrevPrev << " " << jPrevPrev << " " << kPrevPrev << endl;
-// cout << " p.: " << idHexPrev << " " << iPrev << " " << jPrev << " " << kPrev << endl;
-// cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
-
- if (idHexPrevPrev != -1)
- {
- //vPrevPrevPrev.Print() ; cout << endl;
- //vPrevPrev.Print() ; cout << endl;
- //vPrev.Print() ; cout << endl;
- //vNow.Print() ; cout << endl;
-
-// writePoint(vPrevPrevPrev, std::string("testPoint"));
-// D3vector t1 = arrayBoxWSDENT.at(idHexPrevPrevPrev).at(iPrevPrevPrev +(Nx-0)*(jPrevPrevPrev +(Nx-0)* kPrevPrevPrev)).at(0);
-// D3vector t2 = arrayBoxWSDENT.at(idHexPrevPrevPrev).at(iPrevPrevPrev +(Nx-0)*(jPrevPrevPrev +(Nx-0)* kPrevPrevPrev)).at(1);
-// writeBox(t1,t2,std::string("testBox"));
-
-// writePoint(vPrevPrev, std::string("testPoint"));
-// t1 = arrayBoxWSDENT.at(idHexPrevPrev).at(iPrevPrev +(Nx-0)*(jPrevPrev +(Nx-0)* kPrevPrev)).at(0);
-// t2 = arrayBoxWSDENT.at(idHexPrevPrev).at(iPrevPrev +(Nx-0)*(jPrevPrev +(Nx-0)* kPrevPrev)).at(1);
-// writeBox(t1,t2,std::string("testBox"));
-
-
-// writePoint(vPrev, std::string("testPoint"));
-// t1 = arrayBoxWSDENT.at(idHexPrev).at(iPrev +(Nx-0)*(jPrev +(Nx-0)* kPrev)).at(0);
-// t2 = arrayBoxWSDENT.at(idHexPrev).at(iPrev +(Nx-0)*(jPrev +(Nx-0)* kPrev)).at(1);
-// writeBox(t1,t2,std::string("testBox"));
-
-
-// writePoint(vNow, std::string("testPoint"));
-// if (idHexNow != -1)
-// {
-// t1 = arrayBoxWSDENT.at(idHexNow).at(kNow +(Nx-0)*(jNow +(Nx-0)* kNow)).at(0);
-// t2 = arrayBoxWSDENT.at(idHexNow).at(kNow +(Nx-0)*(jNow +(Nx-0)* kNow)).at(1);
-// }
-// writeBox(t1,t2,std::string("testBox"));
-// boxCounter = 0;
- }
-
-
- setPrevIndex();
- setPrevPrevIndex();
- setPrevPrevPrevIndex();
- counterSlow++;
- // cout << "slow version"<<endl;
-}
-
-void Interpolate_direct::setPrevIndex()
-{
- vPrev = vNow;
- idHexPrev = idHexNow;
- iPrev = iNow ;
- jPrev = jNow ;
- kPrev = kNow ;
-}
-
-void Interpolate_direct::setPrevPrevIndex()
-{
- vPrevPrev = vNow;
- idHexPrevPrev = idHexNow;
- iPrevPrev = iNow ;
- jPrevPrev = jNow ;
- kPrevPrev = kNow ;
-}
-
-void Interpolate_direct::setPrevPrevPrevIndex()
-{
- vPrevPrevPrev = vNow;
- idHexPrevPrevPrev = idHexNow;
- iPrevPrevPrev = iNow ;
- jPrevPrevPrev = jNow ;
- kPrevPrevPrev = kNow ;
-}
-
-int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
-{
- if (id_Hex < 0 )
- {
- return -1;
- }
- int Nx = blockgrid->Give_Nx_hexahedron(id_Hex);
- int Ny = blockgrid->Give_Ny_hexahedron(id_Hex);
- int Nz = blockgrid->Give_Nz_hexahedron(id_Hex);
- if ((i < -1 && j < -1 && k < -1 ) || (i > Nx && j > Ny && k > Nz ) || ( (i == -1 || i == Nx ) && (j == -1 || j == Ny) && (k == -1 || k == Nz) ) )
- {
- return -1;
- }
- int typ = -1;
- if ( (i == -1 || j == -1 || k == -1 ) || (i == Nx || j == Ny || k == Nz ) )
- {
- typ = checkForHexaNeighbours(id_Hex, i, j, k, v);
- if ( -1 !=typ)
- {
- return typ;
- }
- }
-
- // figure out, whether intersection with other hex is true;
-
-
-
- if ( id_Hex< 0 || i < 0 || j < 0 || k < 0)
- {
- return -1;
- }
- if (i >= Nx || j >= Ny || k >= Nz )
- {
- return -1;
- }
- checkCounter++;
- D3vector cWSD, cESD;
- D3vector cWND, cEND;
-
- D3vector cWST, cEST;
- D3vector cWNT, cENT;
-
- D3vector boxWSD, boxENT;
-
- D3vector ploc;
-
- if ( i +Nx*(j +Ny* k) >= Nx*Ny*Nz)
- {
- cout << "Nx " << Nx << endl;
- cout << "Ny " << Ny << endl;
- cout << "i " << i << endl;
- cout << "j " << j << endl;
- cout << "k " << k << endl;
- cout << "id of array: " << i +Nx*(j +Ny* k) << endl;
- }
- boxWSD = arrayBoxWSDENT.at(id_Hex).at(i +Nx*(j +Ny* k)).at(0);
- boxENT = arrayBoxWSDENT.at(id_Hex).at(i +Nx*(j +Ny* k)).at(1);
-
-
- if (boxENT>v && boxWSD<v)
- {
-
- if (debugTest)
- {
- //cout << "debug " << endl;
- }
- // cout << "inside box!" << endl;
- 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);
-
- D3vector lam, lamTemp ;
-
-
-
- lam = lambda_of_p_in_tet(v,cWND,cWNT,cWST,cEST);
- lamTemp = lam;
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=0;typCounter0++;}
-
- lam = lambda_of_p_in_tet(v,cEST,cWND,cWST,cESD);
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=1;typCounter1++;lamTemp = lam;}
-
- lam = lambda_of_p_in_tet(v,cWND,cWSD,cWST,cESD);
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=2;typCounter2++;lamTemp = lam;}
-
- lam = lambda_of_p_in_tet(v,cEST,cWND,cESD,cEND);
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=3;typCounter3++;lamTemp = lam;}
-
- lam = lambda_of_p_in_tet(v,cENT,cWNT,cEST,cEND);
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=4;typCounter4++;lamTemp = lam;}
-
- lam = lambda_of_p_in_tet(v,cWNT,cWND,cEST,cEND);
- if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=5;typCounter5++;lamTemp = lam;}
-
- if( typ != -1 && new_lam_better(lambda,lamTemp))
- {
- if (MIN(lamTemp) < -0.1)
- {
- cout << "error !" << endl;
- }
- //lambda = lam;
- lambda = lamTemp;
- typ_tet = typ;
- idHexNow = id_Hex;
- iNow = i;
- jNow = j;
- kNow = k;
- }
- }
- return typ;
-}
-
-int Interpolate_direct::checkBoxSurrounding(int idHex, int i, int j, int k, D3vector v)
-{
- int iterationPerSurrounding = 0;
- for(int kIter=k-1;kIter<k+2;++kIter)
- for(int jIter=j-1;jIter<j+2;++jIter)
- for(int iIter=i-1;iIter<i+2;++iIter) {
- if (checkBox(idHex,iIter, jIter, kIter, v) != -1){
- //counterFast++;
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- return -1;
-}
-
-int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v)
-{
- int iterationPerSurrounding = 0;
- //std::vector<
- for(int iIter=i-1;iIter<i+2;iIter+=2){
- if (checkBox(idHex,iIter, j, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- for(int jIter=j-1;jIter<j+2;jIter+=2){
- if (checkBox(idHex,i, jIter, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- for(int kIter=k-1;kIter<k+2;kIter+=2){
- if (checkBox(idHex,i, j, kIter, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- for(int kIter=k-1;kIter<k+2;kIter+=2)
- for(int jIter=j-1;jIter<j+2;jIter+=2) {
- if (checkBox(idHex,i, jIter, kIter, v) != -1){
- //counterFast++;
-
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- for(int kIter=k-1;kIter<k+2;kIter+=2)
- for(int iIter=i-1;iIter<i+2;iIter+=2) {
- if (checkBox(idHex,iIter, j, kIter, v) != -1){
- //counterFast++;
-
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- for(int jIter=j-1;jIter<j+2;jIter+=2)
- for(int iIter=i-1;iIter<i+2;iIter+=2) {
- if (checkBox(idHex,iIter, jIter, k, v) != -1){
- //counterFast++;
-
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
-
-
-
- for(int kIter=k-1;kIter<k+2;kIter+=2)
- for(int jIter=j-1;jIter<j+2;jIter+=2)
- for(int iIter=i-1;iIter<i+2;iIter+=2) {
- if (checkBox(idHex,iIter, jIter, kIter, v) != -1){
- //counterFast++;
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
-
-
-
-
- return -1;
-}
-
-int Interpolate_direct::checkCorner(int idHex, int i, int j, int k, D3vector v)
-{
- int typ;
-
- int Nx = blockgrid->Give_Nx_hexahedron(idHex);
- int Ny = blockgrid->Give_Ny_hexahedron(idHex);
- int Nz = blockgrid->Give_Nz_hexahedron(idHex);
-
- int index = i +(Nx)*(j +(Nx)* k);
- if (array_box_boundary.at(idHex).at(index).empty())
- return -1;
- int idHexCorner1 = array_box_boundary.at(idHex).at(index).at(0);
- int iCorner1 = array_box_boundary.at(idHex).at(index).at(1);
- int jCorner1 = array_box_boundary.at(idHex).at(index).at(2);
- int kCorner1 = array_box_boundary.at(idHex).at(index).at(3);
-
-
- D3vector boxWSDInner = arrayBoxWSDENT.at(idHex).at(i +Nx*(j +Ny* k)).at(0);
- D3vector boxENTInner = arrayBoxWSDENT.at(idHex).at(i +Nx*(j +Ny* k)).at(1);
-
- D3vector boxWSDOuter = arrayBoxWSDENT.at(idHexCorner1).at(iCorner1 +Nx*(jCorner1 +Ny* kCorner1)).at(0);
- D3vector boxENTOuter = arrayBoxWSDENT.at(idHexCorner1).at(iCorner1 +Nx*(jCorner1 +Ny* kCorner1)).at(1);
- //writeBox(boxWSDInner, boxENTInner, std::string("boxCorner"));
- //writeBox(boxWSDOuter, boxENTOuter, std::string("boxCorner"));
-
- if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
- {
- cout << "surrounding box is wrong!" << endl;
- }
-
- //cout << "corner box 1 : ";
-
- //cout << "id: " << idHexCorner1 << " " << iCorner1 << " " << jCorner1 << " " << kCorner1 << endl;
- //boxWSDOuter.Print();cout << endl;
- //boxENTOuter.Print();cout << endl;
-
- typ = checkBox(idHexCorner1, iCorner1, jCorner1, kCorner1, v);
- if (typ != -1)
- {
- counterCorner++;
- return typ;
- }
-
- if (array_box_boundary.at(idHex).at(index).size() < 5)
- return -1;
-
- int idHexCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Nx)* k)).at(4);
- int iCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Nx)* k)).at(5);
- int jCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Nx)* k)).at(6);
- int kCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Nx)* k)).at(7);
-
- boxWSDOuter = arrayBoxWSDENT.at(idHexCorner2).at(iCorner2 +Nx*(jCorner2 +Ny* kCorner2)).at(0);
- boxENTOuter = arrayBoxWSDENT.at(idHexCorner2).at(iCorner2 +Nx*(jCorner2 +Ny* kCorner2)).at(1);
-
-
-
- if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
- {
- cout << "surrounding box is wrong!" << endl;
- }
-
-
- typ = checkBox(idHexCorner2, iCorner2, jCorner2, kCorner2, v);
- if (typ != -1)
- {
- counterCorner++;
- return typ;
- }
-
- return typ;
-
-}
-
-int Interpolate_direct::checkEdge(int idHex, int i, int j, int k, D3vector v)
-{
-
- int typ;
-
- int Nx = blockgrid->Give_Nx_hexahedron(idHex);
- int Ny = blockgrid->Give_Ny_hexahedron(idHex);
- int Nz = blockgrid->Give_Nz_hexahedron(idHex);
-
- int index = i +(Nx)*(j +(Nx)* k);
- if (array_box_boundary.at(idHex).at(index).empty())
- return -1;
- int idHexEdge = array_box_boundary.at(idHex).at(index).at(0);
- int iEdge = array_box_boundary.at(idHex).at(index).at(1);
- int jEdge = array_box_boundary.at(idHex).at(index).at(2);
- int kEdge = array_box_boundary.at(idHex).at(index).at(3);
-
- D3vector boxWSDInner = arrayBoxWSDENT.at(idHex).at(i +Nx*(j +Ny* k)).at(0);
- D3vector boxENTInner = arrayBoxWSDENT.at(idHex).at(i +Nx*(j +Ny* k)).at(1);
-
- D3vector boxWSDOuter = arrayBoxWSDENT.at(idHexEdge).at(iEdge +Nx*(jEdge +Ny* kEdge)).at(0);
- D3vector boxENTOuter = arrayBoxWSDENT.at(idHexEdge).at(iEdge +Nx*(jEdge +Ny* kEdge)).at(1);
-
- //writeBox(boxWSDInner, boxENTInner, std::string("boxEdge"));
- //writeBox(boxWSDOuter, boxENTOuter, std::string("boxEdge"));
-
- if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
- {
- cout << "surrounding box is wrong!" << endl;
- }
-
- //cout << "edge box : ";
- //cout << "id: " << idHexEdge << " " << iEdge << " " << jEdge << " " << kEdge << endl;
- //boxWSDOuter.Print();cout << endl;
- //boxENTOuter.Print();cout << endl;
-
- typ = checkBox(idHexEdge, iEdge, jEdge, kEdge, v);
- if (typ != -1)
- {
- counterEdge++;
- return typ;
- }
-
- return typ;
-
-}
-
-bool Interpolate_direct::checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT)
-{
- if ( (( (vWSD.x <= wENT.x) && (vWSD.x >= wWSD.x) ) || ( (vENT.x <= wENT.x) && (vENT.x >= wWSD.x) ) )
- && (( (vWSD.y <= wENT.y) && (vWSD.y >= wWSD.y) ) || ( (vENT.y <= wENT.y) && (vENT.y >= wWSD.y) ) )
- && (( (vWSD.z <= wENT.z) && (vWSD.z >= wWSD.z) ) || ( (vENT.z <= wENT.z) && (vENT.z >= wWSD.z) ) ) )
- {
- return true;
- }
-
- if ( (( (wWSD.x <= vENT.x) && (wWSD.x >= vWSD.x) ) || ( (wENT.x <= vENT.x) && (wENT.x >= vWSD.x) ) )
- && (( (wWSD.y <= vENT.y) && (wWSD.y >= vWSD.y) ) || ( (wENT.y <= vENT.y) && (wENT.y >= vWSD.y) ) )
- && (( (wWSD.z <= vENT.z) && (wWSD.z >= vWSD.z) ) || ( (wENT.z <= vENT.z) && (wENT.z >= vWSD.z) ) ) )
- {
- return true;
- }
-
- return false;
-}
-
-void Interpolate_direct::writeBox(D3vector vWSD, D3vector vENT, std::string str)
-{
-
- std::stringstream ss;
- ss << std::setw(5) << std::setfill('0') << boxCounter;
- std::string s = ss.str();
- ofstream myfile;
- myfile.open ("C:/Users/rall/Documents/UG_Blocks_thermal/testergebnisse/" + s + str + ".vtk");
-myfile << "# vtk DataFile Version 4.2\n";
-myfile << "vtk output\n";
-myfile << "ASCII\n";
-myfile << "DATASET POLYDATA\n";
-myfile << "POINTS 24 float\n";
-myfile << vWSD.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << vWSD.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << vWSD.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vWSD.z << "\n" ;
-myfile << vWSD.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vWSD.y << " " << vENT.z << "\n" ;
-myfile << vWSD.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << vENT.x << " " << vENT.y << " " << vENT.z << "\n" ;
-myfile << "\n" ;
-myfile << "POLYGONS 6 30\n";
-myfile << "4 0 1 3 2 \n";
-myfile << "4 4 6 7 5 \n";
-myfile << "4 8 10 11 9 \n";
-myfile << "4 12 13 15 14 \n";
-myfile << "4 16 18 19 17 \n";
-myfile << "4 20 21 23 22 \n";
-
- myfile.close();
-
- boxCounter++;
-}
-
-void Interpolate_direct::writePoint(D3vector v, std::string str)
-{
- std::stringstream ss;
- ss << std::setw(5) << std::setfill('0') << boxCounter;
- std::string s = ss.str();
- ofstream myfile;
- myfile.open ("C:/Users/rall/Documents/UG_Blocks_thermal/testergebnisse/" + s + str + ".vtk");
- myfile << "# vtk DataFile Version 4.2\n";
- myfile << "vtk output\n";
- myfile << "ASCII\n";
- myfile << "DATASET POLYDATA\n";
- myfile << "POINTS 1 float\n";
- myfile << v.x << " " << v.y << " " << v.z << "\n" ;
- myfile << "VERTICES 1 2\n";
- myfile << "1 0 \n";
- myfile.close();
-}
-
-int Interpolate_direct::checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v)
-{
- counterHexa++;
- int Nx = blockgrid->Give_Nx_hexahedron(idHex);
- int Ny = blockgrid->Give_Ny_hexahedron(idHex);
- int Nz = blockgrid->Give_Nz_hexahedron(idHex);
-
- int typ = -1;
- if ( i == -1 || i == Nx )
- {
- if ( j == -1 || j == Ny )
- {
- int Tj = (j<0) ? (Tj = j + 1 ) : (Tj = j - 1 );
- int Ti = (i<0) ? (Ti = i + 1) : (Ti = i - 1);
- debugTest = true;
- typ = checkCorner(idHex, Ti, Tj, k, v);
- debugTest = false;
- if (typ != -1 ) return typ;
-
- }
- else if (k == -1 || k == Nz )
- {
- int Tk = (k<0) ? (Tk = k + 1) : (Tk = k - 1);
- int Ti = (i<0) ? (Ti = i + 1) : (Ti = i - 1);
- debugTest = true;
- typ = checkCorner(idHex, Ti, j, Tk, v);
- debugTest = false;
- if (typ != -1 ) return typ;
- }
- else
- {
- int Ti = (i<0) ? (Ti = i + 1) : (Ti = i - 1);
- typ = checkEdge(idHex, Ti, j, k, v);
- if (typ != -1 ) return typ;
- }
- }
- else if ( j == -1 || j == Ny )
- {
- /*if ( i == -1)
- {
- checkCorner(id_Hex, i, j, k, v);
- }
- else */if (k == -1 || k == Nz )
- {
- int Tk = (k<0) ? (Tk = k + 1) : (Tk = k - 1);
- int Tj = (j<0) ? (Tj = j + 1) : (Tj = j - 1);
- debugTest = true;
- typ = checkCorner(idHex, i, Tj, Tk, v);
- debugTest = false;
- if (typ != -1 ) return typ;
- }
- else
- {
- int Tj = (j<0) ? (Tj = j + 1) : (Tj = j - 1);
- typ = checkEdge(idHex, i, Tj, k, v);
- if (typ != -1 ) return typ;
- }
- }
- else if ( k == -1 || k == Nz )
- {
- int Tk = (k<0) ? (Tk = k + 1) : (Tk = k - 1);
- typ = checkEdge(idHex, i, j, Tk, v);
- if (typ != -1 ) return typ;
- }
- counterHexa--;
- return typ;
-}
diff --git a/program/source/interpol/interpol.h b/program/source/interpol/interpol.h
index 68256edd46deb95526db4d0090ecf30a13c53acb..1fac9a84f70291b499a999d8605cf05381c3bd1b 100644
--- a/program/source/interpol/interpol.h
+++ b/program/source/interpol/interpol.h
@@ -18,7 +18,7 @@
// ------------------------------------------------------------
//
-// variable.h
+// interpolate.h
//
// ------------------------------------------------------------
@@ -38,23 +38,23 @@
/*
data structure on structured grid
-ny * * * * * *
-ny-1 * + + + + *
+ny-1 * * * * * *
+ * + + + + *
3 * + + + + *
2 * + + + + *
1 * + + + + *
0 * * * * * *
- 0 1 2 3 nx-1 nx
+ 0 1 2 3 nx-1
*/
-/***
+/**
* interpolates data from blockgrid_ to rectangular block [pWSD, pENT]
-***/
+**/
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
@@ -62,38 +62,42 @@ class Interpolate_on_structured_grid {
@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_);
- Interpolate_on_structured_grid(int nx_, int ny_, int nz_,
- Blockgrid& blockgrid_);
-/***
+/**
* 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);
- public:
- int nx, ny, nz;
- D3vector pENT,pWSD;
- private:
+
+/**
+ * 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(Cell_variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
+
+
+ ~Interpolate_on_structured_grid();
+private:
int* ids_hex;
int *ids_i, *ids_j, *ids_k;
int *typ_tet;
D3vector *lambda;
-
+ int nx, ny, nz;
double hx, hy, hz;
Blockgrid* blockgrid;
Unstructured_grid* ug;
};
-
-
template <class DTyp>
class Give_corner_data_of_cube {
public:
@@ -140,22 +144,22 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
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];
- }
-*/
+ 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);
@@ -174,38 +178,70 @@ void Interpolate_on_structured_grid::interpolate(Variable<DTyp>& u, DTyp* data,
<< endl;
*/
- if(typ==0) data[ind_global] = interpolate_in_tet(lambda[ind_global],
+ 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],
+ 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],
+ 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],
+ 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],
+ if(typ==5) data[ind_global] = interpolate_in_tet(lambda[ind_global],
du.WNT(),du.WND(),du.EST(),du.END());
}
}
}
+
+
+template <class DTyp>
+void Interpolate_on_structured_grid::interpolate(Cell_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];
+
+ if(id_hex < 0) data[ind_global] = defaultInterpolation;
+ else {
+ Blockgrid* blockgrid = u.getBlockgrid();
+ int Nx = blockgrid->Give_Nx_hexahedron(id_hex);
+ int Ny = blockgrid->Give_Ny_hexahedron(id_hex);
+ int Nz = blockgrid->Give_Nz_hexahedron(id_hex);
+
+ data[ind_global] = u.Give_cell_hexahedra(id_hex,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
+ }
+ }
+}
+
/////////////////////////////////////////////////////////////
// 2. Interpolate from blockgrid to blockgrid
/////////////////////////////////////////////////////////////
-/***
- * interpolates data from blockgrid_from to blockgrid_to_ using an internal rectangular grid if size nx,ny,nz
-***/
+/**
+ * interpolates data from blockgrid_ to rectangular block [pWSD, 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();
@@ -233,211 +269,4 @@ private:
};
-
-/////////////////////////////////////////////////////////////
-// 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;
-};
-
- * 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);
-
- 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);
-
-
-
-
-
- /**
- * 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
- ***/
- 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 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);
-
-private:
- int nx, ny, nz;
- double hx, hy, hz;
- double shiftx, shifty, shiftz;
- Interpolate_on_structured_grid* interpolatorStructured;
- double* data;
-
- D3vector pWSD;
- D3vector pENT;
-
- double defaultInterpolation;
-
-};
-
-
-//template <class DTyp>
-class Interpolate_direct {
- public:
-
- Interpolate_direct (Blockgrid* bg):blockgrid(bg)
- ,idHexPrev(-1)
- ,iPrev(-1)
- ,jPrev(-1)
- ,kPrev(-1)
- ,idHexPrevPrev(-1)
- ,iPrevPrev(-1)
- ,jPrevPrev(-1)
- ,kPrevPrev(-1)
- ,idHexPrevPrevPrev(-1)
- ,iPrevPrevPrev(-1)
- ,jPrevPrevPrev(-1)
- ,kPrevPrevPrev(-1)
- ,counterFast(0)
- ,counterFastest(0)
- ,counterSecondTry(0)
- ,counterSlow(0)
- ,counterHexa(0)
- ,checkCounter(0)
- ,counterEdge(0)
- ,counterCorner(0)
- ,typCounter0(0)
- ,typCounter1(0)
- ,typCounter2(0)
- ,typCounter3(0)
- ,typCounter4(0)
- ,typCounter5(0){}
-
- void init();
- void initBoundary();
-
- void find_cell(D3vector v);
- void setPrevIndex();
- void setPrevPrevIndex();
- void setPrevPrevPrevIndex();
- int checkBox(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);
- int checkCorner(int idHex, int i, int j, int k, D3vector v);
- int checkEdge(int idHex, int i, int j, int k, D3vector v);
- bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
- void writeBox(D3vector vWSD, D3vector vENT, std::string str);
- void writePoint(D3vector v, std::string str);
- int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
- // void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
- template <class DTyp>
- DTyp evaluate(Variable<DTyp>& u);
-
- int counterFast, counterFastest, counterSecondTry, counterSlow, counterHexa, counterCorner, counterEdge;
- int checkCounter;
- int boxCounter = 0;
- int typCounter0,typCounter1,typCounter2,typCounter3,typCounter4,typCounter5;
- bool debugTest;
- D3vector lambda;
- D3vector vNow, vPrev, vPrevPrev, vPrevPrevPrev;
- private:
- int idHexPrevPrevPrev, iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev;
- int idHexPrevPrev, iPrevPrev, jPrevPrev, kPrevPrev;
- int idHexPrev, iPrev, jPrev, kPrev;
- int idHexNow, iNow, jNow, kNow;
-
- int typ_tet;
-
-
- Blockgrid *blockgrid;
- std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
-
- std::vector<std::vector<std::vector<int> > > array_box_boundary;
- // template <class DTyp>
- // Variable<DTyp>* u;
-};
-
-
-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<ug->Give_number_hexahedra();++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());
- }
- return returnVal;
-};
-
-
#endif // INTERPOL_H
diff --git a/program/source/interpol/interpolCellVar.cc b/program/source/interpol/interpolCellVar.cc
index 702cd42616214c8ea7770efefb3042a495a0e372..5f999213350244d1c7310ad0c4fce03b0caf2286 100644
--- a/program/source/interpol/interpolCellVar.cc
+++ b/program/source/interpol/interpolCellVar.cc
@@ -15,18 +15,20 @@
* See the License for the specific language governing permissions and
* limitations under the License.
**********************************************************************************/
-#include "../../../Source_Colsamm/Colsamm.h"
-using namespace ::_COLSAMM_;
+
+
#include "interpolCellVar.h"
#include "../extemp/variable_cc.h"
#include "../pde_op/diffop_cc.h"
-
+#include "../../../Source_Colsamm/Colsamm.h"
#include "assert.h"
+using namespace ::_COLSAMM_;
+
template<>
void Interpolate_from_cell_to_point<Stencil_variable>::init(Blockgrid* blockgrid_) {