diff --git a/program/source/grid/marker.cc b/program/source/grid/marker.cc
index 3ec9ad273ed7cad0f1240c14a1928d42ac25ec2b..1a6f38c75b8623cd1fbd7bdca337cf7dd00e5c33 100644
--- a/program/source/grid/marker.cc
+++ b/program/source/grid/marker.cc
@@ -467,7 +467,7 @@ Boundary_Marker::Boundary_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ )
else {
Set_marker<quadrangleEl> ( id, yes_mark );
for ( int i=0; i<4; ++i )
- Set_marker<edgeEl> ( quad->Give_id_edge ( ( dir2D ) i ), yes_mark );
+ Set_marker<edgeEl> ( quad->Give_id_edge ( ( dir2D ) i ), yes_mark );
for ( int i=0; i<4; ++i )
Set_marker<pointEl> ( quad->Give_id_corner ( ( dir2D_sons ) i ), yes_mark );
}
diff --git a/program/source/interpol/interpol.cc b/program/source/interpol/interpol.cc
index 17d6ecef3629342f71cba0c78f81d83596cc5921..1c4a441db51a44866367654a65b36112b36b4d05 100644
--- a/program/source/interpol/interpol.cc
+++ b/program/source/interpol/interpol.cc
@@ -61,6 +61,14 @@ bool contained_in_tet_strong(D3vector lam) {
if(lam.x + lam.y + lam.z > 1+limit) return false;
return true;
}
+bool contained_in_tet_direct(D3vector lam) {
+ double limit = 0.0; // 0.2
+ if(lam.x < -limit) return false;
+ if(lam.y < -limit) return false;
+ if(lam.z < -limit) return false;
+ if(lam.x + lam.y + lam.z > 1+limit) return false;
+ return true;
+}
bool new_lam_better(D3vector lam_old, D3vector lam_new) {
@@ -1637,13 +1645,47 @@ void PointInterpolator::smoothGrid()
{
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;
+ double div = 0;
+ double sum = 0;
+ if (data[i +nx*(j +ny* k)] != defaultInterpolation)
+ {
+ sum += 2.0*data[i +nx*(j +ny* k)] ;
+ div++;div++;
+ }
+ if (data[i+1 +nx*(j +ny* k)] != defaultInterpolation)
+ {
+ sum += data[i+1 +nx*(j +ny* k)] ;
+ div++;
+ }
+ if (data[i +nx*(j +1 +ny* k)] != defaultInterpolation)
+ {
+ sum += data[i +nx*(j +1 +ny* k)] ;
+ div++;
+ }
+ if (data[i +nx*(j +ny* (k+1 ))] != defaultInterpolation)
+ {
+ sum += data[i +nx*(j +ny* (k+1 ))] ;
+ div++;
+ }
+ if (data[i-1 +nx*(j +ny* k)] != defaultInterpolation)
+ {
+ sum += data[i-1 +nx*(j +ny* k)] ;
+ div++;
+ }
+ if (data[i +nx*(j -1 +ny* k)] != defaultInterpolation)
+ {
+ sum += data[i +nx*(j -1 +ny* k)] ;
+ div++;
+ }
+ if (data[i +nx*(j +ny* (k-1 ))] != defaultInterpolation)
+ {
+ sum += data[i +nx*(j +ny* (k-1 ))] ;
+ div++;
+ }
+ if (div >0)
+ {
+ data[i +nx*(j +ny* k)] = sum / div;
+ }
}
}
@@ -1882,6 +1924,7 @@ void Interpolate_direct::init()
{
arrayBoxWSDENT.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
array_box_boundary.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
+ isOuterBoundaryFlag.resize(blockgrid->Give_unstructured_grid()->Give_number_hexahedra());
int counter = 0;
int counterTwoNeighbours = 0;
@@ -1894,11 +1937,13 @@ void Interpolate_direct::init()
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!)
@@ -2210,7 +2255,56 @@ void Interpolate_direct::init()
}
}
-// cout << "counter " << counter << endl;
+
+ 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);
+ isOuterBoundaryFlag.at(idHex).resize((Nx+1)*(Ny+1)*(Nz+1));
+
+ for (int i = Nx/2 ; i < Nx+1; i = i + Nx) for (int j = Ny/2 ; j < Ny+1; j = j + Ny) for (int k = 0 ; k < Nz+1; k = k + Nz)
+ {
+ int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
+ if (blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(idHex).at(convertedLocalIndex).empty())
+ {
+ for (int iInner = 0 ; iInner < Nx+1; iInner++) for (int jInner = 0 ; jInner < Ny+1; jInner++)
+ {
+ int convertedLocalIndex = iInner +(Nx+1)*(jInner +(Ny+1)* k);
+ isOuterBoundaryFlag.at(idHex).at(convertedLocalIndex) = true;
+ }
+ }
+ }
+ for (int i = Nx/2 ; i < Nx+1; i = i + Nx) for (int j = 0 ; j < Ny+1; j = j + Ny) for (int k = Nz/2 ; k < Nz+1; k = k + Nz)
+ {
+ int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
+ if (blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(idHex).at(convertedLocalIndex).empty())
+ {
+ for (int iInner = 0 ; iInner < Nx+1; iInner++) for (int kInner = 0 ; kInner < Nz+1; kInner++)
+ {
+ int convertedLocalIndex = iInner +(Nx+1)*(j +(Ny+1)* kInner);
+ isOuterBoundaryFlag.at(idHex).at(convertedLocalIndex) = true;
+ }
+ }
+
+ }
+ for (int i = 0 ; i < Nx+1; i = i + Nx) for (int j = Ny/2 ; j < Ny+1; j = j + Ny) for (int k = Nz/2 ; k < Nz+1; k = k + Nz)
+ {
+ int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
+ if (blockgrid->Give_blockgrid_coordinates()->blockgrid_hexa_boundary.at(idHex).at(convertedLocalIndex).empty())
+ {
+ for (int jInner = 0 ; jInner < Ny+1; jInner++) for (int kInner = 0 ; kInner < Nz+1; kInner++)
+ {
+ int convertedLocalIndex = i +(Nx+1)*(jInner +(Ny+1)* kInner);
+ isOuterBoundaryFlag.at(idHex).at(convertedLocalIndex) = true;
+ }
+ }
+
+ }
+ }
+
+
+ cout << "counter " << counter << endl;
// cout << "boundary boxes " << counterBoxesAtBoundary<< endl;
// cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
// cout << "counterTwoNeighbours " << counterTwoNeighbours << endl;
@@ -2915,6 +3009,8 @@ std::vector<int> Interpolate_direct::calculateNeighbourIndex(std::vector<std::ve
iTemp = iTemp + relation.at(0).at(1);
jTemp = jTemp + relation.at(0).at(2);
kTemp = kTemp + relation.at(0).at(3);
+
+
//mark, which one was changed --> invert sign
if (relation.at(0).at(1) != 0)
{ iTemp = -iTemp; }
@@ -3418,7 +3514,186 @@ void Interpolate_direct::find_cell(D3vector v)
setPrevPrevIndex();
setPrevPrevPrevIndex();
counterSlow++;
- // cout << "slow version"<<endl;
+ // cout << "slow version"<<endl;
+}
+
+void Interpolate_direct::find_surface_cell(D3vector v)
+{
+// vPrevPrev = vPrev;
+// vPrev = vNow;
+ if (v == vNow)
+ {
+ counterSamePoint++;
+ return;
+ }
+ badLambdaFound = false;
+ bool badIsBestLambda = false;
+ bool previousBadLambda = false;
+ idHexNow = -1;
+ vNow = v;
+ lambda = D3vector(-1,-1,-1);
+
+
+
+ D3vector boxWSD, boxENT;
+
+ //cout << "looking for new cell!" << endl;
+
+ int Nx, Ny, Nz;
+ // prev
+ if (checkBoxSurface(idHexPrev,iPrev, jPrev, kPrev, v) != -1){
+ if (!badLambdaFound)
+ {
+ counterFastest++;
+ setPrevIndex();
+ return;
+ }
+ else
+ {
+ previousBadLambda = true;
+ //std::cout << "should be in next or is at boundary!" << std::endl;
+ }
+ }
+ //std::cout << idHexPrev<< " " <<iPrev<< " " << jPrev<< " " << kPrev << std::endl;
+ if (checkBoxSurroundingSurface(idHexPrev,iPrev, jPrev, kPrev, v) != -1)
+ {
+ if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
+ {
+ //will stay bad : probably at the boundary!
+ counterFast++;
+ setPrevIndex();
+ return;
+ }
+
+
+ }
+ // prev prev
+ if (checkBoxSurface(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1)
+ {
+ if (!badLambdaFound)
+ {
+ setPrevIndex();
+ setPrevPrevIndex();
+ //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
+ counterSecondTry++;
+ return;
+ }
+ else
+ {
+ previousBadLambda = true;
+ //std::cout << "should be in next or is at boundary!" << std::endl;
+ }
+
+ }
+
+ if (checkBoxSurroundingSurface(idHexPrevPrev,iPrevPrev, jPrevPrev, kPrevPrev, v) != -1)
+ {
+ if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
+ {
+ setPrevIndex();
+ setPrevPrevIndex();
+ //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
+ counterSecondTry++;
+ return;
+ }
+
+ }
+ // prev prev prev
+ if (checkBoxSurface(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1)
+ {
+ if (!badLambdaFound)
+ {
+ setPrevIndex();
+ setPrevPrevIndex();
+ setPrevPrevPrevIndex();
+ //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
+ counterThirdTry++;
+ return;
+ }
+ else
+ {
+ previousBadLambda = true;
+ //std::cout << "should be in next or is at boundary!" << std::endl;
+ }
+
+
+ }
+
+ if (checkBoxSurroundingSurface(idHexPrevPrevPrev,iPrevPrevPrev, jPrevPrevPrev, kPrevPrevPrev, v) != -1)
+ {
+ if ((previousBadLambda && badLambdaFound ) || !badLambdaFound)
+ {
+ setPrevIndex();
+ setPrevPrevIndex();
+ setPrevPrevPrevIndex();
+ //cout << " n.: " << idHexNow << " " << kNow << " " << jNow << " " << kNow << endl;
+ counterThirdTry++;
+ return;
+ }
+
+ }
+
+
+
+ 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+=(Nz-1))
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i) {
+ if (0 < checkBoxSurface(id_hex,i, j, k, v))
+ {
+ if (abs(iNow - iPrev) < 2 && abs(jNow - jPrev) < 2 && abs(kNow - kPrev) < 2)
+ {
+ std::cout<<std::endl;
+ std::cout << "why not found?!\n";
+ }
+ }
+ }
+ for(int k=0;k<Nz;++k)
+ for(int j=0;j<Ny;j+=(Ny-1))
+ for(int i=0;i<Nx;++i) {
+ if (0 < checkBoxSurface(id_hex,i, j, k, v))
+ {
+ if (abs(iNow - iPrev) < 2 && abs(jNow - jPrev) < 2 && abs(kNow - kPrev) < 2)
+ {
+ std::cout<<std::endl;
+ std::cout << "why not found?!\n";
+ }
+ }
+ }
+ for(int k=0;k<Nz;++k)
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;i+=(Nx-1)) {
+ if (0 < checkBoxSurface(id_hex,i, j, k, v))
+ {
+ if (abs(iNow - iPrev) < 2 && abs(jNow - jPrev) < 2 && abs(kNow - kPrev) < 2)
+ {
+ std::cout<<std::endl;
+ std::cout << "why not found?!\n";
+ }
+ }
+ }
+ }
+ if (badLambdaFound)
+ {
+ if (MIN(lambda) < 0 || MAX(lambda)>1)
+ {
+ }
+ else
+ {
+ assert(false && "must not happen!");
+ }
+ //lambda.Print();std::cout<<std::flush;
+ }
+
+ setPrevIndex();
+ setPrevPrevIndex();
+ setPrevPrevPrevIndex();
+ counterSlow++;
+ // cout << "slow version"<<endl;
}
void Interpolate_direct::setPrevIndex()
@@ -3528,7 +3803,6 @@ int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
cENT = blockgrid->Give_coord_hexahedron(id_Hex,i+1,j+1,k+1);
-
D3vector lam, lamTemp ;
@@ -3637,39 +3911,331 @@ int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
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, bool neglectBadLambda)
+int Interpolate_direct::checkBoxSurface(int id_Hex, int i, int j, int k, D3vector v)
{
- bool temp = badLambdaFound;
- badLambdaFound = false;
- int iterationPerSurrounding = 0;
- //std::vector<
- for(int kIter=k+1;kIter>k-2;kIter-=2){
- if (checkBox(idHex,i, j, kIter, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- if (badLambdaFound && !neglectBadLambda)
+ if (id_Hex < 0 )
{
- if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ return -1;
}
-
- for(int iIter=i-1;iIter<i+2;iIter+=2){
+ 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 = checkForHexaNeighboursSurface(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);
+
+ // writeBox(boxWSD,boxENT,std::string("box"));
+
+
+ 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);
+
+ bool cWSDFlag = isOuterBoundaryFlag.at(id_Hex).at((i+0) +(Nx+1)*((j+0) +(Ny+1)* (k+0)));
+ bool cESDFlag = isOuterBoundaryFlag.at(id_Hex).at((i+1) +(Nx+1)*((j+0) +(Ny+1)* (k+0)));
+ bool cWNDFlag = isOuterBoundaryFlag.at(id_Hex).at((i+0) +(Nx+1)*((j+1) +(Ny+1)* (k+0)));
+ bool cENDFlag = isOuterBoundaryFlag.at(id_Hex).at((i+1) +(Nx+1)*((j+1) +(Ny+1)* (k+0)));
+
+ bool cWSTFlag = isOuterBoundaryFlag.at(id_Hex).at((i+0) +(Nx+1)*((j+0) +(Ny+1)* (k+1)));
+ bool cESTFlag = isOuterBoundaryFlag.at(id_Hex).at((i+1) +(Nx+1)*((j+0) +(Ny+1)* (k+1)));
+ bool cWNTFlag = isOuterBoundaryFlag.at(id_Hex).at((i+0) +(Nx+1)*((j+1) +(Ny+1)* (k+1)));
+ bool cENTFlag = isOuterBoundaryFlag.at(id_Hex).at((i+1) +(Nx+1)*((j+1) +(Ny+1)* (k+1)));
+
+
+ bool sideD = false, sideT = false, sideN = false, sideS = false, sideE = false, sideW = false;
+ if ( cWSDFlag && cESDFlag && cWNDFlag && cENDFlag )
+ {
+ sideD = true;
+ }
+ else if ( cWSTFlag && cESTFlag && cWNTFlag && cENTFlag )
+ {
+ sideT = true;
+ }
+ else if ( cWSDFlag && cWNDFlag && cWSTFlag && cWNTFlag )
+ {
+ sideW = true;
+ }
+ else if ( cESDFlag && cENDFlag && cESTFlag && cENTFlag )
+ {
+ sideE = true;
+ }
+ else if ( cWNDFlag && cENDFlag && cWNTFlag && cENTFlag )
+ {
+ sideN = true;
+ }
+ else if ( cWSDFlag && cESDFlag && cWSTFlag && cESTFlag )
+ {
+ sideS = true;
+ }
+ else
+ {
+ std::cout << "error in checkBoxBoundary\n";
+ }
+
+
+
+
+
+ D3vector lam, lamTemp ;
+
+
+ //std::cout << "\n\nstart of tet search\n\n";
+// if (sideW || sideT)
+// {
+// lam = lambda_of_p_in_tet(v,cWND,cWNT,cWST,cEST);
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=0;typCounter0++;lamTemp = lam;}
+// }
+
+// //std::cout << "new best : "; lamTemp.Print();
+// if (sideS)
+// {
+// lam = lambda_of_p_in_tet(v,cEST,cWND,cWST,cESD);
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=1;typCounter1++;lamTemp = lam;}
+
+// }
+// if (sideW || sideD || sideS)
+// {
+// std::cout << std::setprecision(15);
+// lam = lambda_of_p_in_tet(v,cWND,cWSD,cWST,cESD);
+// interpolate_in_tet_D3vector(lam,cWND,cWSD,cWST,cESD).Print();std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,4.0,4.0,5.0,4.0);std::cout<<std::flush;
+// lam.y = 0;
+// interpolate_in_tet_D3vector(lam,cWND,cWSD,cWST,cESD).Print();std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,4.0,4.0,5.0,4.0);std::cout<<std::flush;
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=2;typCounter2++;lamTemp = lam;}
+// }
+
+// if (sideE || sideD)
+// {
+// lam = lambda_of_p_in_tet(v,cEST,cWND,cESD,cEND);
+
+// interpolate_in_tet_D3vector(lam,cEST,cWND,cESD,cEND).Print();std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+// interpolate_in_tet_D3vector(lam,cEST,cWND,cESD,cEND).Print();std::cout<<std::flush;
+// std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=3;typCounter3++;lamTemp = lam;}
+
+// }
+
+// if (sideE || sideT || sideN)
+// {
+// lam = lambda_of_p_in_tet(v,cENT,cWNT,cEST,cEND);
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=4;typCounter4++;lamTemp = lam;}
+
+// }
+
+// if (sideN)
+// {
+// lam = lambda_of_p_in_tet(v,cWNT,cWND,cEST,cEND);
+// if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=5;typCounter5++;lamTemp = lam;}
+// }
+
+ lamTemp = lambda;
+ D3vector zero(0,0,0);
+ D3vector lam2;
+ D3vector lamTrilinear;
+ if (sideN)
+ {
+ lam = lambda_of_p_in_tet(v,cWND,cWNT,cEND,zero);
+ lam2 = lambda_of_p_in_tet(v,cENT,cWNT,cEND,zero);
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=0;typCounter0++;lamTemp = lam;}
+ }
+
+ //std::cout << "new best : "; lamTemp.Print();
+ if (sideS)
+ {
+ lam = lambda_of_p_in_tet(v,cESD,cWSD,cEST,zero);
+ lam2 = lambda_of_p_in_tet(v,cWST,cWSD,cEST,zero);
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=1;typCounter1++;lamTemp = lam;}
+
+ }
+ if (sideE)
+ {
+ lam = lambda_of_p_in_tet(v,cEST,cESD,cENT,zero);
+ lam2 = lambda_of_p_in_tet(v,cEND,cESD,cENT,zero);
+ lam.z = 0;
+ lam2.z = 0;
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=7;typCounter2++;lamTemp = lam;}
+ if (contained_in_tet_strong(lam2) && new_lam_better(lamTemp,lam2)){typ=8;typCounter2++;lamTemp = lam2;}
+ }
+
+ if (sideW)
+ {
+ lam = lambda_of_p_in_tet(v,cEST,cWND,cESD,cEND);
+
+ interpolate_in_tet_D3vector(lam,cEST,cWND,cESD,cEND).Print();std::cout<<std::flush;
+ std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+ interpolate_in_tet_D3vector(lam,cEST,cWND,cESD,cEND).Print();std::cout<<std::flush;
+ std::cout << interpolate_in_tet(lam,5.0,4.0,4.0,4.0);std::cout<<std::flush;
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=3;typCounter3++;lamTemp = lam;}
+
+ }
+
+ if (sideT)
+ {
+
+ lam = lambda_of_p_in_tet(v,cENT,cWNT,cEST,zero);
+ lam2 = lambda_of_p_in_tet(v,cWST,cWNT,cEST,zero);
+ lam.z = 0;
+ lam2.z = 0;
+
+ if (false && (contained_in_tet_direct(lam) || contained_in_tet_direct(lam2)) && false)
+ {
+ lamTrilinear = bilinearInterpolation(v,cENT,cWNT,cEST,cWST);
+ lamTemp = lamTrilinear;
+ typ = 14;
+ }
+ else
+ {
+
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=11;typCounter2++;lamTemp = lam;}
+ if (contained_in_tet_strong(lam2) && new_lam_better(lamTemp,lam2)){typ=12;typCounter2++;lamTemp = lam2;}
+ }
+
+ }
+
+ if (sideD)
+ {
+ lam = lambda_of_p_in_tet(v,cWND,cEND,cWSD,zero);
+ lam2= lambda_of_p_in_tet(v,cESD,cEND,cWSD,zero);
+ lam.z = 0;
+ lam2.z = 0;
+ if (false && (contained_in_tet_direct(lam) || contained_in_tet_direct(lam2)) && false)
+ {
+ lamTrilinear = bilinearInterpolation(v,cWND,cEND,cWSD,cESD);
+ lamTemp = lamTrilinear;
+ typ = 13;
+ }
+ else
+ {
+ if (contained_in_tet_strong(lam) && new_lam_better(lamTemp,lam)){typ=9;typCounter2++;lamTemp = lam;}
+ if (contained_in_tet_strong(lam2) && new_lam_better(lamTemp,lam2)){typ=10;typCounter2++;lamTemp = lam2;}
+ }
+ }
+ //std::cout << "new best : "; lamTemp.Print();
+ if( typ != -1 && new_lam_better(lambda,lamTemp))
+ {
+ if (MIN(lamTemp) < 0.0 || MAX(lamTemp) >1.0 )
+ {
+ badLambdaFound = true;
+ //lamTemp.Print();
+ }
+ else
+ {
+ badLambdaFound = false;
+ }
+
+ 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, bool neglectBadLambda)
+{
+ bool temp = badLambdaFound;
+ badLambdaFound = false;
+ int iterationPerSurrounding = 0;
+ //std::vector<
+ for(int kIter=k+1;kIter>k-2;kIter-=2){
+ if (checkBox(idHex,i, j, kIter, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+ for(int iIter=i-1;iIter<i+2;iIter+=2){
if (checkBox(idHex,iIter, j, k, v) != -1 && !badLambdaFound){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
return 0;}
else { iterationPerSurrounding++;}
@@ -3688,69 +4254,260 @@ int Interpolate_direct::checkBoxSurroundingOptimized(int idHex, int i, int j, in
if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
}
- 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 && !badLambdaFound){
- //counterFast++;
+ 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 && !badLambdaFound){
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ 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 && !badLambdaFound){ //here
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ 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 && !badLambdaFound){
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+
+ 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 && !badLambdaFound){
+ //counterFast++;
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+
+
+ badLambdaFound = temp;
+ return -1;
+}
+
+
+
+
+int Interpolate_direct::checkBoxSurroundingSurface(int idHex, int i, int j, int k, D3vector v, bool neglectBadLambda)
+{
+
+ bool temp = badLambdaFound;
+ badLambdaFound = false;
+ int iterationPerSurrounding = 0;
+ //std::vector<
+
+
+ for(int iIter=i-1;iIter<i+2;iIter+=2){
+ if (checkBoxSurface(idHex,iIter, j, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ for(int jIter=j-1;jIter<j+2;jIter+=2){
+ if (checkBoxSurface(idHex,i, jIter, k, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ for(int jIter=j-1;jIter<j+2;jIter+=2)
+ for(int iIter=i-1;iIter<i+2;iIter+=2) {
+ if (checkBoxSurface(idHex,iIter, jIter, k, v) != -1){
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ for(int kIter=k+1;kIter>k-2;kIter-=2){
+ if (checkBoxSurface(idHex,i, j, kIter, v) != -1){//cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+ for(int kIter=k-1;kIter<k+2;kIter+=2)
+ for(int jIter=j-1;jIter<j+2;jIter+=2) {
+ if (checkBoxSurface(idHex,i, jIter, kIter, v) != -1){
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+ for(int kIter=k-1;kIter<k+2;kIter+=2)
+ for(int iIter=i-1;iIter<i+2;iIter+=2) {
+ if (checkBoxSurface(idHex,iIter, j, kIter, v) != -1){ //here
+ //counterFast++;
+
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+ 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 (checkBoxSurface(idHex,iIter, jIter, kIter, v) != -1){
+ //counterFast++;
+ //setPrevIndex();
+ //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
+ return 0;}
+ else { iterationPerSurrounding++;}
+ }
+ if (badLambdaFound && !neglectBadLambda)
+ {
+ if (checkBoxSurroundingSurface(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ }
+
+
+ return -1;
+
+}
+
+
+int Interpolate_direct::checkCornerSurface(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 +(Ny)*(j +(Nx)* k);
+ index = i +Nx*(j +Ny* 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);
+
+ int NxCorner1 = blockgrid->Give_Nx_hexahedron(idHexCorner1);
+ int NyCorner1 = blockgrid->Give_Ny_hexahedron(idHexCorner1);
+ int NzCorner1 = blockgrid->Give_Nz_hexahedron(idHexCorner1);
+
+ D3vector boxWSDOuter = arrayBoxWSDENT.at(idHexCorner1).at(iCorner1 +NxCorner1*(jCorner1 +NyCorner1* kCorner1)).at(0);
+ D3vector boxENTOuter = arrayBoxWSDENT.at(idHexCorner1).at(iCorner1 +NxCorner1*(jCorner1 +NyCorner1* kCorner1)).at(1);
+ //writeBox(boxWSDInner, boxENTInner, std::string("boxCorner"));
+ //writeBox(boxWSDOuter, boxENTOuter, std::string("boxCorner"));
+
+ if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
+ {
+ writeBox(boxWSDInner, boxENTInner, std::string("box1"));
+ writeBox(boxWSDOuter, boxENTOuter, std::string("box1"));
+ 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 = checkBoxSurface(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 +(Ny)* k)).at(4);
+ int iCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Ny)* k)).at(5);
+ int jCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Ny)* k)).at(6);
+ int kCorner2 = array_box_boundary.at(idHex).at(i +(Nx)*(j +(Ny)* k)).at(7);
+
+ int NxCorner2 = blockgrid->Give_Nx_hexahedron(idHexCorner2);
+ int NyCorner2 = blockgrid->Give_Ny_hexahedron(idHexCorner2);
+ int NzCorner2 = blockgrid->Give_Nz_hexahedron(idHexCorner2);
+
+ boxWSDOuter = arrayBoxWSDENT.at(idHexCorner2).at(iCorner2 +NxCorner2*(jCorner2 +NyCorner2* kCorner2)).at(0);
+ boxENTOuter = arrayBoxWSDENT.at(idHexCorner2).at(iCorner2 +NxCorner2*(jCorner2 +NyCorner2* kCorner2)).at(1);
+
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- if (badLambdaFound && !neglectBadLambda)
- {
- if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
- }
- 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 && !badLambdaFound){ //here
- //counterFast++;
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- if (badLambdaFound && !neglectBadLambda)
+ if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
{
- if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ writeBox(boxWSDInner, boxENTInner, std::string("box1"));
+ writeBox(boxWSDOuter, boxENTOuter, std::string("box1"));
+ cout << "surrounding box is wrong!" << endl;
}
- 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 && !badLambdaFound){
- //counterFast++;
-
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- if (badLambdaFound && !neglectBadLambda)
- {
- if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
- }
- 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 && !badLambdaFound){
- //counterFast++;
- //setPrevIndex();
- //cout << "iterationPerSurrounding " << iterationPerSurrounding+1<< endl;
- return 0;}
- else { iterationPerSurrounding++;}
- }
- if (badLambdaFound && !neglectBadLambda)
+ typ = checkBoxSurface(idHexCorner2, iCorner2, jCorner2, kCorner2, v);
+ if (typ != -1)
{
- if (checkBoxSurroundingOptimized(idHexNow,iNow, jNow, kNow, v,true) != -1) return 0; else return -1;
+ counterCorner++;
+ return typ;
}
+ return typ;
-
- badLambdaFound = temp;
- return -1;
}
int Interpolate_direct::checkCorner(int idHex, int i, int j, int k, D3vector v)
@@ -3840,6 +4597,64 @@ int Interpolate_direct::checkCorner(int idHex, int i, int j, int k, D3vector v)
}
+
+int Interpolate_direct::checkEdgeSurface(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);
+ index = (i)+Nx*(j +Ny*(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);
+
+ int NxInner = blockgrid->Give_Nx_hexahedron(idHexEdge);
+ int NyInner = blockgrid->Give_Ny_hexahedron(idHexEdge);
+ int NzInner = blockgrid->Give_Nz_hexahedron(idHexEdge);
+
+ int indexInner = iEdge +NxInner*(jEdge +NyInner* kEdge);
+ //int indexInnerTest = iEdge +NyInner*(jEdge +NxInner* kEdge);
+
+ D3vector boxWSDOuter = arrayBoxWSDENT.at(idHexEdge).at(indexInner).at(0);
+ D3vector boxENTOuter = arrayBoxWSDENT.at(idHexEdge).at(indexInner).at(1);
+
+ //writeBox(boxWSDInner, boxENTInner, std::string("boxEdge"));
+ //writeBox(boxWSDOuter, boxENTOuter, std::string("boxEdge"));
+
+ if (!checkOverlapOfBoxes(boxWSDInner, boxENTInner, boxWSDOuter, boxENTOuter))
+ {
+ writeBox(boxWSDInner, boxENTInner, std::string("box1"));
+ writeBox(boxWSDOuter, boxENTOuter, std::string("box2"));
+ 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 = checkBoxSurface(idHexEdge, iEdge, jEdge, kEdge, v);
+ if (typ != -1)
+ {
+ counterEdge++;
+ return typ;
+ }
+
+ return typ;
+
+}
+
int Interpolate_direct::checkEdge(int idHex, int i, int j, int k, D3vector v)
{
@@ -4057,6 +4872,127 @@ void Interpolate_direct::writeInterpolationEfficiency()
}
}
+
+D3vector Interpolate_direct::bilinearInterpolation(D3vector v, D3vector x1, D3vector x2, D3vector x3, D3vector x4)
+{
+
+
+ //P = A + x*b + y*C + xy*D
+ D3vector A = x1; // 100 - 000
+ D3vector B = x2-x1; // 100 - 000
+ D3vector C = x3-x1; // 010 - 000
+ D3vector D = x4-x2+x1-x3; // 001 - 000
+ //std::cout << "to be implemented" << std::endl;
+
+ D3vector normalTD = cross_product(x1,B);
+ D3vector normalTDOPP = cross_product(x3 , x4-x3);
+ //D3vector normalTDOPP2 = cross_product(cENT-cEST ,cWNT - cENT);
+ D3vector normalNS = cross_product(x1,C);
+ D3vector normalNSOPP = cross_product(x2,x4-x2);
+ //D3vector normalNSOPP2 = cross_product(cESD-cEST,cESD - cWSD);
+ normalTD = normalTD / D3VectorNorm(normalTD);
+ normalTDOPP = normalTDOPP / D3VectorNorm(normalTDOPP);
+ normalNS = normalNS / D3VectorNorm(normalNS);
+ normalNSOPP = normalNSOPP / D3VectorNorm(normalNSOPP);
+ double eta = 0.4;
+ double xi = 0.5;
+ double phi = 0.5;
+ bool fixEta = false;
+ bool fixXi = false;
+ bool fixPhi = false;
+
+ double distTD0 = product(normalTD,x1);
+ double distTD1 = product(normalTDOPP,x1);
+ double distTDMid = distTD0-product(normalTD,v);
+ double distTDMid2 = distTD1-product(normalTDOPP,v);
+
+ double distNS0 = product(normalNS,x4);
+ double distNS1 = product(normalNSOPP,x4);
+ double distNSMid = distNS0-product(normalNS,v);
+ double distNSMid2 = distNS1 -product(normalNSOPP,v);
+
+
+
+ //coord.Print();std::cout<<std::endl;
+ D3vector coord(0.5);
+ coord.x = 0.45;
+ coord.y = 0.55;
+ if (fabs(angle_between_vectors(normalTD,normalTDOPP)-180) < 0.5)
+ {
+ double delta = (distTD1 - distTD0);
+ phi = (distTDMid + delta - distTDMid2) /delta/ 2.0;
+ coord.z = phi;
+ fixPhi = true;
+ }
+ if (fabs(angle_between_vectors(normalNS,normalNSOPP)-180) < 0.5)
+ {
+ double delta = (distNS1 - distNS0);
+ eta = (distNSMid + delta - distNSMid2) /delta / 2.0;
+ coord.x = eta;
+ fixEta = true;
+ }
+
+
+
+
+
+ D3vector R = A + B * coord.x + C * coord.y + D * coord.x * coord.y;
+
+ {
+ double factorJac = 0.8;
+ for (int iter = 0 ; iter < 20 ; iter++)
+ {
+ D3vector partialX = B + D * coord.y;
+ D3vector partialY = C + D * coord.x;
+
+ D3matrix jac2(partialX,partialY,D3vector(0.0,0.0,1.0));
+ jac2.transpose();
+ jac2.invert_gauss_elimination();
+
+ R = A + B * coord.x + C * coord.y + D * coord.x * coord.y -v;
+
+
+ //std::cout << "residuum f_xn" <<D3VectorNorm(R)<<std::endl;
+ if (D3VectorNormSquared(R) < 1e-10)
+ {
+ iter = 1000;
+ }
+
+ coord = coord - jac2.vectorMultiply(R) * factorJac;
+ if (MAX(coord) > 1.0 || MIN(coord)<0)
+ {
+ factorJac = factorJac * 0.5;
+ if (coord.x < 0.0)
+ coord.x = 0.1;
+ if (coord.y < 0.0)
+ coord.y = 0.15;
+ if (coord.z < 0.0)
+ coord.z = 0.0;
+ if (coord.x > 1.0)
+ coord.x = 0.95;
+ if (coord.y > 1.0)
+ coord.y = 0.9;
+ if (coord.z > 1.0)
+ coord.z = 1.0;
+ }
+ //coord.Print();std::cout<<std::endl;
+
+ }
+ }
+
+
+ D3vector x = A + B * coord.x + C * coord.y + D * coord.x * coord.y;
+// X.Print();
+// std::cout << "tirlinear interp \n";
+// x.Print();
+// std::cout << std::flush;
+
+ if (MAX(coord) > 1.01 || MIN(coord) < -0.01)
+ {
+ return D3vector{-1,-1,-1};
+ }
+ return coord;
+}
bool isInRange(double d0, double d1)
{
if (d0 * d1 >=0 )
@@ -4304,6 +5240,75 @@ D3vector Interpolate_direct::trilinarInterpolation(D3vector X, int id_Hex, int i
}
+
+int Interpolate_direct::checkForHexaNeighboursSurface(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 = checkCornerSurface(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 = checkCornerSurface(idHex, Ti, j, Tk, v);
+ debugTest = false;
+ if (typ != -1 ) return typ;
+ }
+ else
+ {
+ int Ti = (i<0) ? (Ti = i + 1) : (Ti = i - 1);
+ typ = checkEdgeSurface(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 = checkCornerSurface(idHex, i, Tj, Tk, v);
+ debugTest = false;
+ if (typ != -1 ) return typ;
+ }
+ else
+ {
+ int Tj = (j<0) ? (Tj = j + 1) : (Tj = j - 1);
+ typ = checkEdgeSurface(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 = checkEdgeSurface(idHex, i, j, Tk, v);
+ if (typ != -1 ) return typ;
+ }
+ counterHexa--;
+ return typ;
+}
+
int Interpolate_direct::checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v)
{
counterHexa++;
diff --git a/program/source/interpol/interpol.h b/program/source/interpol/interpol.h
index 1f6bd95b774c983432edb0a7a2e96aeecfd29f78..f581b48579140cee3a7569ba56d2a12ae78857e4 100644
--- a/program/source/interpol/interpol.h
+++ b/program/source/interpol/interpol.h
@@ -418,6 +418,7 @@ class Interpolate_direct {
*
**/
void find_cell(D3vector v);
+ void find_surface_cell(D3vector v);
/**
* Print a vtk file with a box, surrounding the cell.
**/
@@ -457,9 +458,12 @@ class Interpolate_direct {
int getTet(){return typ_tet;}
int getHex(){return idHexNow;}
+ D3vector bilinearInterpolation(D3vector v,D3vector a,D3vector b,D3vector c,D3vector d);
D3vector trilinarInterpolation(D3vector v, int id_Hex,int i, int j, int k );
bool vectorInBox(D3vector vWSD, D3vector vENT, D3vector v, double eps = 1e-10);
int checkBox(int idHex, int i, int j, int k, D3vector v);
+ int checkBoxSurface(int idHex, int i, int j, int k, D3vector v);
+
private:
int idHexPrev{-1}, iPrev{-1}, jPrev{-1}, kPrev{-1};
int idHexPrevPrevPrev{-1}, iPrevPrevPrev{-1}, jPrevPrevPrev{-1}, kPrevPrevPrev{-1};
@@ -477,6 +481,7 @@ class Interpolate_direct {
D3vector cWNT, cENT;
Blockgrid *blockgrid;
std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
+ std::vector<std::vector<bool > > isOuterBoundaryFlag;
std::vector<std::vector<std::vector<int> > > array_box_boundary;
std::vector<std::vector<std::vector<int> > > array_point_boundary;
@@ -502,10 +507,14 @@ private:
void setPrevPrevIndex();
void setPrevPrevPrevIndex();
int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
+ int checkForHexaNeighboursSurface(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, bool neglectBadLambda = false);
+ int checkBoxSurroundingSurface(int idHex, int i, int j, int k, D3vector v, bool neglectBadLambda = false);
int checkCorner(int idHex, int i, int j, int k, D3vector v);
+ int checkCornerSurface(int idHex, int i, int j, int k, D3vector v);
int checkEdge(int idHex, int i, int j, int k, D3vector v);
+ int checkEdgeSurface(int idHex, int i, int j, int k, D3vector v);
bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
double calculateOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
@@ -553,6 +562,19 @@ DTyp Interpolate_direct::evaluate(Variable<DTyp> &u)
du.WNT(),du.WND(),du.EST(),du.END());
if(typ==6) returnVal = interpolate_in_tet_trilinear(lambda,du.END(),du.ESD(),du.WSD(),du.WND(),
du.WNT(),du.WST(),du.EST(),du.ENT());
+ if(typ==7) returnVal = interpolate_in_tet(lambda,du.EST(),du.ESD(),du.ENT(),0.0);
+ if(typ==8) returnVal = interpolate_in_tet(lambda,du.END(),du.ESD(),du.ENT(),0.0);
+ if(typ==9) returnVal = interpolate_in_tet(lambda,du.WND(),du.END(),du.WSD(),0.0);
+ if(typ==10) returnVal = interpolate_in_tet(lambda,du.ESD(),du.END(),du.WSD(),0.0);
+ if(typ==11) returnVal = interpolate_in_tet(lambda,du.ENT(),du.WNT(),du.EST(),0.0);
+ if(typ==12) returnVal = interpolate_in_tet(lambda,du.WST(),du.WNT(),du.EST(),0.0);
+ if(typ==13) returnVal = interpolate_in_tet_bilinear(lambda,du.WND(),du.END(),du.WSD(),du.ESD());
+ if(typ==14) returnVal = interpolate_in_tet_bilinear(lambda,du.ENT(),du.WNT(),du.WST(),du.EST());
+
+
+
+
+
}
// D3vector cWSD = {1,1,0}; // 1,1,0 : x2
// D3vector cESD = {2,0,-1}; // 1,0,0 : x1
diff --git a/program/source/math_lib/math_lib.cc b/program/source/math_lib/math_lib.cc
index 350dff42728f40b558290cf65a1b8f67d95045d9..7563dc803aadcd15654ba903371a19e4d0f1f2d1 100644
--- a/program/source/math_lib/math_lib.cc
+++ b/program/source/math_lib/math_lib.cc
@@ -54,6 +54,34 @@ double calc_maximal_face_angle(const D3vector& v0,
angle_between_faces(v0,v1,v3,v2)));
}
+D3matrix rotationMatrix(D3vector vIn, D3vector vOut)
+{
+ //calculates rotation matrix, which rotates vIn, such that is points into direction vOut
+ if (D3VectorNorm(vIn) == 0 || D3VectorNorm(vOut) == 0)
+ {
+ return D3matrix();
+ }
+ D3vector cross = cross_product(vOut,vIn);
+ double normed = D3VectorNorm(cross_product(vOut,vIn));
+ D3vector orthogonal = cross / normed;
+ double alpha = atan2(normed,product(vOut,vIn));
+
+ double s = sin(alpha);
+ double c = cos(alpha);
+ double mc = 1.0 - c;
+ double x = orthogonal.x;
+ double y = orthogonal.y;
+ double z = orthogonal.z;
+
+
+ D3vector cx(c + x * x * mc,x * y * mc + z * s,x * z * mc - y * s); // first column
+ D3vector cy(x * y * mc - z * s,c + y * y * mc,y * z * mc + x * s); // second column
+ D3vector cz(x * z * mc + y * s,y * z * mc - x * s,c + z * z * mc); // third column
+
+ D3matrix ret(cx,cy,cz);
+ ret.transpose();
+ return ret;
+}
////////////////////////////////////
// lambda of p in tet
@@ -72,6 +100,7 @@ D3vector lambda_of_p_in_tet(D3vector p,
}
+
////////////////////////////////////
// D3vector
////////////////////////////////////
diff --git a/program/source/math_lib/math_lib.h b/program/source/math_lib/math_lib.h
index dfc3f6453a9d92c2f9b3582563226f29511058cf..419718c880637f7c911d26eefd369ceb5505c723 100644
--- a/program/source/math_lib/math_lib.h
+++ b/program/source/math_lib/math_lib.h
@@ -153,7 +153,12 @@ class D3matrix {
D3matrix(D3vector cx, D3vector cy, D3vector cz) :
x1(cx.x), y1(cy.x), z1(cz.x),
x2(cx.y), y2(cy.y), z2(cz.y),
- x3(cx.z), y3(cy.z), z3(cz.z) {};
+ x3(cx.z), y3(cy.z), z3(cz.z) {}
+
+ D3matrix() :
+ x1(1), y1(0), z1(0),
+ x2(0), y2(1), z2(0),
+ x3(0), y3(0), z3(1) {}
double Determinante() {
return (x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2));
@@ -225,6 +230,9 @@ class D3matrix {
return D3matrix(cx,cy,cz);
}
+
+
+
D3vector vectorMultiply(D3vector m)
{
// D3vector res = {x1 * m.x + x2 * m.x + x3 * m.x,
@@ -300,7 +308,7 @@ class D3matrix {
for (int i = order - 1; i > 0; i--) {
//Swapping each and every element of the two rows
- if (a[i - 1][0] < a[i][0] && a[i][0]!=0)
+ if (a[i - 1][0] < a[i][0])// && a[i][0]!=0)
for (int j = 0; j < 2 * order; j++) {
// Swapping of the row, if above
@@ -335,6 +343,10 @@ class D3matrix {
for (int k = 0; k < 2 * order; k++) {
a[j][k] -= a[i][k] * temp;
+// if (fabs(a[j][k]) < 1e-20 )
+// {
+// std::cout << "is small???";
+// }
}
}
}
@@ -343,14 +355,18 @@ class D3matrix {
// Multiply each row by a nonzero integer.
// Divide row element by the diagonal element
for (int i = 0; i < order; i++) {
- if (a[i][i] == 0)
- {
- std::cout << "debug here";
- }
+// if (a[i][i] == 0)
+// {
+// std::cout << "debug here";
+// }
temp = a[i][i];
for (int j = 0; j < 2 * order; j++) {
a[i][j] = a[i][j] / temp;
+// if (a[i][j] == 0)
+// {
+// std::cout << "is zero here";
+// }
}
}
// std::cout << "inverted matrix : " << std::endl;
@@ -392,10 +408,17 @@ D3vector lambda_of_p_in_tet(D3vector p,
template <class DTyp>
DTyp interpolate_in_tet(D3vector lambda,
DTyp vA, DTyp vB,
- DTyp vC, DTyp vD) {
+ DTyp vC, DTyp vD) {
+ return vA + (vB-vA) * lambda.x + (vC-vA) * lambda.y + (vD-vA) * lambda.z;
+}
+
+inline D3vector interpolate_in_tet_D3vector(D3vector lambda,
+ D3vector vA, D3vector vB,
+ D3vector vC, D3vector vD) {
return vA + (vB-vA) * lambda.x + (vC-vA) * lambda.y + (vD-vA) * lambda.z;
}
+
template <class DTyp>
DTyp interpolate_in_tet_trilinear(D3vector lambda,
DTyp x0, DTyp x1,
@@ -417,6 +440,20 @@ DTyp interpolate_in_tet_trilinear(D3vector lambda,
// return vA + (vB-vA) * lambda.x + (vC-vA) * lambda.y + (vD-vA) * lambda.z;
}
+template <class DTyp>
+DTyp interpolate_in_tet_bilinear(D3vector lambda,
+ DTyp x1, DTyp x2,
+ DTyp x3, DTyp x4)
+{
+ DTyp A = x1;
+ DTyp B = x2-x1;
+ DTyp C = x3-x1;
+ DTyp D = x4-x2+x1-x3;
+ DTyp R = A + lambda.x * B + lambda.y * C + lambda.x * lambda.y * D;
+ return R;
+}
+
+
//////////////////////////////////////////////////////////////////////
// 3. geometric operators for 3D vectors
@@ -483,6 +520,7 @@ inline double angle_between_vectors_rad(const D3vector& va,
return acos(var);
}
+
inline double max_interior_angel_of_triangle(const D3vector& va,
const D3vector& vb,
const D3vector& vc) {
@@ -513,6 +551,7 @@ double calc_maximal_face_angle(const D3vector& va,
+D3matrix rotationMatrix(D3vector vIn, D3vector vOut);
//////////////////////////////////////////////
// Implementierung einiger Memberfunktionen