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

added interpolator, special for surface only! faster than general...

added interpolator, special for surface only! faster than general interpolator, especially helpfor for interpolation of surface deformation
parent 4f3e9058
......@@ -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 );
}
......
......@@ -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);