Commit eadd7402 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

interpolator now trilinear

parent 5ebf0313
...@@ -19,5 +19,8 @@ source/program.pro ...@@ -19,5 +19,8 @@ source/program.pro
*.stash *.stash
*pro.user* *pro.user*
*pro.user *pro.user
Makefile
*.opt
*.user
/program2D/Makefile /program2D/Makefile
/program2D/debug/ /program2D/debug/
\ No newline at end of file
...@@ -1288,6 +1288,22 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary() ...@@ -1288,6 +1288,22 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
void Blockgrid_coordinates::init_blockgrid_coordinates_boundary_fast() void Blockgrid_coordinates::init_blockgrid_coordinates_boundary_fast()
{ {
bool cantUseFast = false;
int numfree = bg->Give_unstructured_grid()->degree_of_freedom();
for (int iter = 0 ; iter < numfree; iter++)
{
if(bg->getNumberPointsDegree(iter)%2 !=0)
{
cantUseFast = true;
}
}
if (cantUseFast)
{
std::cout << "odd gridpoint number, cant use fast method to initialize boundary points. Slow version used instead.\n";
init_blockgrid_coordinates_boundary();
return;
}
blockgrid_hexa_boundary.clear(); blockgrid_hexa_boundary.clear();
counter = 0; counter = 0;
blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra()); blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
......
...@@ -382,13 +382,13 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i ...@@ -382,13 +382,13 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i
// X.x = X.x +15; // X.x = X.x +15;
// X.y = X.y +15; // X.y = X.y +15;
D3vector normalTD = cross_product(-1*B,C); D3vector normalTD = cross_product(-1*B,C);
D3vector normalTDOPP = cross_product(cWST-cEST , cWNT - cWST); D3vector normalTDOPP = -1*cross_product(cWST-cEST , cWNT - cWST);
//D3vector normalTDOPP2 = cross_product(cENT-cEST ,cWNT - cENT); //D3vector normalTDOPP2 = cross_product(cENT-cEST ,cWNT - cENT);
D3vector normalNS = cross_product(-1*C,D); D3vector normalNS = cross_product(-1*C,D);
D3vector normalNSOPP = cross_product(cWST-cEST,cWST - cWSD); D3vector normalNSOPP = -1*cross_product(cWST-cEST,cWST - cWSD);
//D3vector normalNSOPP2 = cross_product(cESD-cEST,cESD - cWSD); //D3vector normalNSOPP2 = cross_product(cESD-cEST,cESD - cWSD);
D3vector normalEW = cross_product(-1*D,B); D3vector normalEW = cross_product(-1*D,B);
D3vector normalEWOPP = cross_product(cWST-cWSD,cWST - cWNT); D3vector normalEWOPP = -1*cross_product(cWST-cWSD,cWST - cWNT);
//D3vector normalEWOPP2 = cross_product(cWND-cWSD,cWND - cWNT); //D3vector normalEWOPP2 = cross_product(cWND-cWSD,cWND - cWNT);
normalTD = normalTD / D3VectorNorm(normalTD); normalTD = normalTD / D3VectorNorm(normalTD);
normalTDOPP = normalTDOPP / D3VectorNorm(normalTDOPP); normalTDOPP = normalTDOPP / D3VectorNorm(normalTDOPP);
...@@ -405,6 +405,9 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i ...@@ -405,6 +405,9 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i
D3vector coord(eta,xi,phi); D3vector coord(eta,xi,phi);
double distTD0 = product(normalTD,cEND); double distTD0 = product(normalTD,cEND);
double TEST1 = product(normalTD,cWND);
double TEST2 = product(normalTD,cESD);
double TEST3 = product(normalTD,cWSD);
double distTD1 = product(normalTDOPP,cWST); double distTD1 = product(normalTDOPP,cWST);
double distTDMid = distTD0-product(normalTD,X); double distTDMid = distTD0-product(normalTD,X);
double distTDMid2 = distTD1-product(normalTDOPP,X); double distTDMid2 = distTD1-product(normalTDOPP,X);
...@@ -437,37 +440,79 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i ...@@ -437,37 +440,79 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i
D3vector F = cWNT - cWND - cENT + cEND; // 011 - 010 - 001 + 000 D3vector F = cWNT - cWND - cENT + cEND; // 011 - 010 - 001 + 000
D3vector G = cEST - cESD - cENT + cEND; // 101 - 100 - 001 + 000 D3vector G = cEST - cESD - cENT + cEND; // 101 - 100 - 001 + 000
D3vector H = cWST + cESD + cWND + cENT - cEND - cWSD - cWNT - cEST; // 111 + 100 + 010 + 001 - 000 - 110 - 011 - 101 D3vector H = cWST + cESD + cWND + cENT - cEND - cWSD - cWNT - cEST; // 111 + 100 + 010 + 001 - 000 - 110 - 011 - 101
if (fabs(angle_between_vectors(normalTD,normalTDOPP)-180) < 0.5) if (fabs(angle_between_vectors(normalTD,normalTDOPP)-180*0) < 0.5)
{ {
double delta = (distTD1 - distTD0); double delta = (distTD1 - distTD0);
phi = (distTDMid + delta - distTDMid2) /delta/ 2.0; //phi = (distTDMid + delta - distTDMid2) /delta/ 2.0;
//std::cout << "phi 1 "<< phi << std::endl;
phi = -distTDMid / delta;
//std::cout << "phi 2 "<< phi << std::endl;
//phi = distTDMid2 / delta;
//std::cout << "phi 3 "<< phi << std::endl;
coord.z = phi; coord.z = phi;
fixPhi = true; fixPhi = true;
} }
if (fabs(angle_between_vectors(normalNS,normalNSOPP)-180) < 0.5) if (fabs(angle_between_vectors(normalNS,normalNSOPP)-180*0) < 0.5)
{ {
double delta = (distNS1 - distNS0); double delta = (distNS1 - distNS0);
eta = (distNSMid + delta - distNSMid2) /delta / 2.0; //eta = (distNSMid + delta - distNSMid2) /delta / 2.0;
//std::cout << "eta 1 "<< eta << std::endl;
eta = -distNSMid / delta;
//std::cout << "eta 2 "<< eta << std::endl;
//eta = distNSMid2 / delta;
//std::cout << "eta 3 "<< eta << std::endl;
coord.x = eta; coord.x = eta;
fixEta = true; fixEta = true;
} }
if (fabs(angle_between_vectors(normalEW,normalEWOPP)-180) < 0.5) if (fabs(angle_between_vectors(normalEW,normalEWOPP)-180*0) < 0.5)
{ {
double delta = (distEW1 - distEW0); double delta = (distEW1 - distEW0);
xi = (distEWMid + delta - distEWMid2) /delta/ 2.0; //xi = (distEWMid + delta - distEWMid2) /delta/ 2.0;
//std::cout << "xi 1 "<< xi << std::endl;
xi = -distEWMid / delta;
//std::cout << "xi 2 "<< xi << std::endl;
//xi = distEWMid2 / delta;
//std::cout << "xi 3 "<< xi << std::endl;
//xi = (distEWMid + distEWMid2)/ 2.0 / (distEW1 - distEW0); //xi = (distEWMid + distEWMid2)/ 2.0 / (distEW1 - distEW0);
coord.y = xi; coord.y = xi;
fixXi = true; fixXi = true;
} }
//coord.Print();std::cout<<std::endl; // std::cout << "initial guess\n";
// coord.Print();std::cout<<std::endl;
D3vector R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z; D3vector R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
//
// coord.x = 1; coord.y = 0; coord.z = 0;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 0; coord.y = 1; coord.z = 0;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 0; coord.y = 0; coord.z = 1;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 1; coord.y = 1; coord.z = 0;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 0; coord.y = 1; coord.z = 1;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 1; coord.y = 0; coord.z = 1;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 0; coord.y = 0; coord.z = 0;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();
// coord.x = 1; coord.y = 1; coord.z = 1;
// R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
// R.Print();std::cout<<std::endl; // R.Print();std::cout<<std::endl;
// X.Print(); // X.Print();
D3vector x;
for (int iter = 0 ; iter < 20 ; iter++) D3vector coordPrev;
for (int iter = 0 ; iter < 50 ; iter++)
{ {
coordPrev = coord;
D3vector partialEta = B + E * coord.y + G * coord.z + H * coord.y*coord.z; D3vector partialEta = B + E * coord.y + G * coord.z + H * coord.y*coord.z;
D3vector partialXi = C + E * coord.x + F * coord.z + H * coord.x*coord.z; D3vector partialXi = C + E * coord.x + F * coord.z + H * coord.x*coord.z;
D3vector partialPhi = D + F * coord.y + G * coord.x + H * coord.x*coord.y; D3vector partialPhi = D + F * coord.y + G * coord.x + H * coord.x*coord.y;
...@@ -475,25 +520,30 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i ...@@ -475,25 +520,30 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i
D3matrix jac2(partialEta,partialXi,partialPhi); D3matrix jac2(partialEta,partialXi,partialPhi);
jac2.transpose(); //jac2.transpose();
jac2.invert_gauss_elimination(); //jac2.invert_gauss_elimination();
R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z -X; R = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z -X;
//std::cout << "residuum f_xn" <<D3VectorNorm(R)<<std::endl; //std::cout << "residuum f_xn" <<D3VectorNorm(R)<<std::endl;
if (D3VectorNormSquared(R) < 1e-10) // if (D3VectorNormSquared(R) < 1e-10)
// {
// iter = 1000;
// }
// coord.Print(); std::cout << std::flush;
coord = coord - jac2.invert_apply(R) * 1.0;
if (D3VectorNorm(coordPrev-coord) < 1e-25)
{ {
iter = 1000; iter = 1000;
} }
// coord.Print();
coord = coord - jac2.vectorMultiply(R) * 0.8;
if (std::isnan(coord.x) || std::isnan(coord.y) || std::isnan(coord.z)) if (std::isnan(coord.x) || std::isnan(coord.y) || std::isnan(coord.z))
{ {
std::cout << "trilinear interpolation failed!!!\n"<< std::endl; std::cout << "trilinear interpolation failed!!!\n"<< std::endl;
} }
D3vector x = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z; x = A + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
if (coord.x > 2.0) if (coord.x > 2.0)
coord.x = 1.1; coord.x = 1.1;
if (coord.y > 2.0) if (coord.y > 2.0)
...@@ -517,7 +567,7 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i ...@@ -517,7 +567,7 @@ D3vector Interpolate_on_structured_grid::trilinarInterpolation(D3vector X, int i
// std::cout << std::flush; // std::cout << std::flush;
if (MAX(coord) > 1.01 || MIN(coord) < -0.01) if (MAX(coord) > 1.2 || MIN(coord) < -0.2)
{ {
return D3vector{-1,-1,-1}; return D3vector{-1,-1,-1};
} }
...@@ -530,7 +580,7 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, ...@@ -530,7 +580,7 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
// id = trilinearInterpolationFlag_; // id = trilinearInterpolationFlag_;
trilinearInterpolationFlag = trilinearInterpolationFlag_; trilinearInterpolationFlag = trilinearInterpolationFlag_;
int Nx, Ny, Nz; int Nx, Ny, Nz;
// int typ; int typ;
assert(nx_ > 1); assert(nx_ > 1);
assert(ny_ > 1); assert(ny_ > 1);
...@@ -542,7 +592,7 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_, ...@@ -542,7 +592,7 @@ Interpolate_on_structured_grid::Interpolate_on_structured_grid(int nx_, int ny_,
double factor = 0.1; double factor = 0.1;
// double factor = 0.00001; // double factor = 0.00001;
// D3vector lam; D3vector lam;
blockgrid = &blockgrid_; blockgrid = &blockgrid_;
ug = blockgrid->Give_unstructured_grid(); ug = blockgrid->Give_unstructured_grid();
...@@ -681,9 +731,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -681,9 +731,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
// cout << "HI" << endl; // cout << "HI" << endl;
int typ = -1;
D3vector lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST); typ = -1;
lam = lambda_of_p_in_tet(ploc,cWND,cWNT,cWST,cEST);
if(contained_in_tet(lam)) typ=0; if(contained_in_tet(lam)) typ=0;
else { else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cWST,cESD); lam = lambda_of_p_in_tet(ploc,cEST,cWND,cWST,cESD);
...@@ -693,7 +744,8 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -693,7 +744,8 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
if(contained_in_tet(lam)) typ=2; if(contained_in_tet(lam)) typ=2;
else { else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND); lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);
if(contained_in_tet(lam)) typ=3; if(contained_in_tet(lam))
{typ=3;lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);}
else { else {
lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND); lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND);
if(contained_in_tet(lam)) typ=4; if(contained_in_tet(lam)) typ=4;
...@@ -705,21 +757,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -705,21 +757,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
} }
} }
} }
bool useTrilinear = false; if (trilinearInterpolationFlag && MIN(lam)>= -0.1 && MAX(lam) <= 1.1 && SUM(lam) <= 1.2 && SUM(lam)>=-0.2)
if (trilinearInterpolationFlag && MIN(lam)>= -0.2 && MAX(lam) <= 1.2 )
{ {
// std::cout << "lam before; "; lam = trilinarInterpolation(ploc, id_hex,i, j, k);
// lam.Print(); typ = 6;
// std::cout<<std::endl;
D3vector lamTri = trilinarInterpolation(ploc, id_hex,i, j, k);
if (! (std::isnan(lamTri.x) || std::isnan(lamTri.y) || std::isnan(lamTri.z)) && MIN(lamTri)>= 0.0 && MAX(lamTri))
{
useTrilinear = true;
lam = lamTri;
typ = 6;
}
// lam.Print();
// std::cout<<std::endl;
} }
/* /*
cout << " typ " << typ << id_hex cout << " typ " << typ << id_hex
...@@ -736,11 +777,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -736,11 +777,10 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
stop=false; stop=false;
if(ids_hex[ind_global]!=-1) { if(ids_hex[ind_global]!=-1) {
stop=(new_lam_worse(lambda[ind_global],lam) || !useTrilinear); stop=new_lam_worse(lambda[ind_global],lam);
} }
#pragma omp critical if(stop==false) {
if(stop==false && typ_tet[ind_global] != 6) {
ids_hex[ind_global] = id_hex; ids_hex[ind_global] = id_hex;
ids_i[ind_global] = i; ids_i[ind_global] = i;
ids_j[ind_global] = j; ids_j[ind_global] = j;
...@@ -982,7 +1022,8 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -982,7 +1022,8 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
if(contained_in_tet(lam)) typ=2; if(contained_in_tet(lam)) typ=2;
else { else {
lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND); lam = lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);
if(contained_in_tet(lam)) typ=3; if(contained_in_tet(lam))
{typ=3;lambda_of_p_in_tet(ploc,cEST,cWND,cESD,cEND);}
else { else {
lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND); lam = lambda_of_p_in_tet(ploc,cENT,cWNT,cEST,cEND);
if(contained_in_tet(lam)) typ=4; if(contained_in_tet(lam)) typ=4;
...@@ -994,7 +1035,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -994,7 +1035,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
} }
} }
} }
if (trilinearInterpolationFlag && MIN(lam)>= 0.0 && MAX(lam) <= 1.0 ) if (trilinearInterpolationFlag && MIN(lam)>= -0.2 && MAX(lam) <= 1.2 && SUM(lam) <= 1.2 && SUM(lam)>=-0.2)
{ {
lam = trilinarInterpolation(ploc, id_hex,i, j, k); lam = trilinarInterpolation(ploc, id_hex,i, j, k);
typ = 6; typ = 6;
...@@ -1022,7 +1063,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD. ...@@ -1022,7 +1063,7 @@ if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.
ids_i[ind_global] = i; ids_i[ind_global] = i;
ids_j[ind_global] = j; ids_j[ind_global] = j;
ids_k[ind_global] = k; ids_k[ind_global] = k;
typ_tet[ind_global] = typ; typ_tet[ind_global] = typ;
lambda[ind_global] = lam; lambda[ind_global] = lam;
...@@ -1707,13 +1748,20 @@ void PointInterpolator::resetInterpolator() ...@@ -1707,13 +1748,20 @@ void PointInterpolator::resetInterpolator()
} }
} }
void PointInterpolator::normToNumberOfWritings() void PointInterpolator::normToNumberOfWritings(bool p2norm)
{ {
for (int iter = 0 ; iter < nx*ny*nz;iter++) for (int iter = 0 ; iter < nx*ny*nz;iter++)
{ {
if (dataCounter[iter] >0) if (dataCounter[iter] >0)
{ {
data[iter] = data[iter] / double(dataCounter[iter]); if(p2norm)
{
data[iter] = sqrt(data[iter] / double(dataCounter[iter]));
}
else
{
data[iter] = data[iter] / double(dataCounter[iter]);
}
} }
} }
} }
...@@ -1753,24 +1801,6 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y, ...@@ -1753,24 +1801,6 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y,
double locZ = posZ / hz - k; 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;
data[i +nx*(j +ny* k)] += value; data[i +nx*(j +ny* k)] += value;
data[(i+1)+nx*(j +ny* k)] += value; data[(i+1)+nx*(j +ny* k)] += value;
data[i +nx*((j+1)+ny* k)] += value; data[i +nx*((j+1)+ny* k)] += value;
...@@ -1789,8 +1819,20 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y, ...@@ -1789,8 +1819,20 @@ void PointInterpolator::writeOnInterpolatedGrid(double coord_x, double coord_y,
dataCounter[(i+1)+nx*((j+1)+ny*(k+1))]++; dataCounter[(i+1)+nx*((j+1)+ny*(k+1))]++;
return; return;
}
void PointInterpolator::writeOnInterpolatedGridPoint(int i, int j, int k, double value)
{
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;
data[i +nx*(j +ny* k)] += value;
dataCounter[i +nx*(j +ny* k)]++;
return;
} }
void PointInterpolator::subtractOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value) void PointInterpolator::subtractOnInterpolatedGrid(double coord_x, double coord_y, double coord_z, double value)
...@@ -2310,7 +2352,7 @@ void Interpolate_direct::init() ...@@ -2310,7 +2352,7 @@ void Interpolate_direct::init()
} }
cout << "counter " << counter << endl; // cout << "counter " << counter << endl;
// cout << "boundary boxes " << counterBoxesAtBoundary<< endl; // cout << "boundary boxes " << counterBoxesAtBoundary<< endl;
// cout << "counterTwoNeighbours " << counterTwoNeighbours << endl; // cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
// cout << "counterTwoNeighbours " << counterTwoNeighbours << endl; // cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
...@@ -2994,8 +3036,11 @@ std::vector<std::vector<int> > Interpolate_direct::calculateNeighbourIndexRelati ...@@ -2994,8 +3036,11 @@ std::vector<std::vector<int> > Interpolate_direct::calculateNeighbourIndexRelati
} }
cout << "checking ... Finished! " << endl; cout << "checking ... Finished! " << endl;
cout << endl; cout << endl;
// std::vector<int> emptyVectorInt(4);
// std::vector<std::vector<int> > emptyVector(3);
// emptyVector.at(0) = emptyVectorInt; emptyVector.at(1) = emptyVectorInt; emptyVector.at(2) = emptyVectorInt;
// return emptyVector;
return std::vector< std::vector< int > >(0); return std::vector< std::vector< int > >(0);
} }
...@@ -3335,6 +3380,7 @@ void Interpolate_direct::find_cell(D3vector v) ...@@ -3335,6 +3380,7 @@ void Interpolate_direct::find_cell(D3vector v)
badLambdaFound = false; badLambdaFound = false;
bool badIsBestLambda = false; bool badIsBestLambda = false;
bool previousBadLambda = false; bool previousBadLambda = false;
idHexNowSurface = -1;
idHexNow = -1; idHexNow = -1;
vNow = v; vNow = v;
lambda = D3vector(-1,-1,-1); lambda = D3vector(-1,-1,-1);
...@@ -3535,6 +3581,7 @@ void Interpolate_direct::find_surface_cell(D3vector v) ...@@ -3535,6 +3581,7 @@ void Interpolate_direct::find_surface_cell(D3vector v)
badLambdaFound = false; badLambdaFound = false;
bool badIsBestLambda = false; bool badIsBestLambda = false;
bool previousBadLambda = false; bool previousBadLambda = false;
idHexNowSurface = -1;
idHexNow = -1; idHexNow = -1;
vNow = v; vNow = v;
lambda = D3vector(-1,-1,-1); lambda = D3vector(-1,-1,-1);
...@@ -3547,11 +3594,14 @@ void Interpolate_direct::find_surface_cell(D3vector v) ...@@ -3547,11 +3594,14 @@ void Interpolate_direct::find_surface_cell(D3vector v)
int Nx, Ny, Nz; int Nx, Ny, Nz;
// prev // prev
if (checkBoxSurface(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
//set kPrev to boundary value! Only works here, if k=0 or k=Nz is the refracting surface
// (double(kPrev)/double(blockgrid->Give_Nz_hexahedron(0)))
if (checkBoxSurface(idHexPrevSurface,iPrevSurface, jPrevSurface, kPrevSurface, v) != -1){
if (!badLambdaFound) if (!badLambdaFound)
{ {
counterFastest++; counterFastest++;
setPrevIndex(); setPrevIndexSurface();
return; return;
} }
else else
...@@ -3561,13 +3611,13 @@ void Interpolate_direct::find_surface_cell(D3vector v) ...@@ -3561,13 +3611,13 @@ void Interpolate_direct::find_surface_cell(D3vector v)
} }
} }
//std::cout << idHexPrev<< " " <<iPrev<< " " << jPrev<< " " << kPrev << std::endl; //std::cout << idHexPrev<< " " <<iPrev<< " " << jPrev<< " " << kPrev << std::endl;
if (checkBoxSurroundingSurface(idHexPrev,iPrev, jPrev, kPrev, v) != -1) if (checkBoxSurroundingSurface(idHexPrevSurface,iPrevSurface, jPrevSurface, kPrevSurface, v) != -1)
{ {
if ((previousBadLambda && badLambdaFound ) || !badLambdaFound) if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
{ {
//will stay bad : probably at the boundary! //will stay bad : probably at the boundary!
counterFast++; counterFast++;
setPrevIndex(); setPrevIndexSurface();
return; return;
} }