Skip to content

GitLab

Commit 8953fc87 authored by Phillip Lino Rall's avatar Phillip Lino Rall
Browse files

increased efficiency of interpolator / finding of adjascent cells

parent 5f99e009
Show whitespace changes
Inline Side-by-side
Showing with 333 additions and 7 deletions
program/source/grid/blockgrid.cc
View file @ 8953fc87
......@@ -1073,6 +1073,7 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
{
blockgrid_hexa_boundary.clear();
int counter(0);
blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
//outer loop
......@@ -1088,7 +1089,7 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
//if ((i > 0 || i < Nx-1) ^ (j > 0 || j < Ny-1) ^ (k > 0 || k < Nz-1))
//remove true here?
if (true || (i == 0 || i == Nx-1) || (j== 0 || j == Ny-1) || (k == 0 || k == Nz-1))
if (false || (i == 0 || i == Nx-1) || (j== 0 || j == Ny-1) || (k == 0 || k == Nz-1))
{
//inner loop
for (int id_hexInner = 0 ; id_hexInner < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hexInner++)
......@@ -1108,7 +1109,7 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
{
// not sure which if loop is correct:
//if ((iInner > 0 || iInner < NxInner-1) ^ (jInner > 0 || jInner < NyInner-1) ^ (kInner > 0 || kInner < NzInner-1))
if (true || (iInner == 0 || iInner == NxInner-1) || (jInner == 0 || jInner == NyInner-1) || (kInner == 0 || kInner == NzInner-1))
if (false || (iInner == 0 || iInner == NxInner-1) || (jInner == 0 || jInner == NyInner-1) || (kInner == 0 || kInner == NzInner-1))
{
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner)) < 1e-10 )
{
......@@ -1139,6 +1140,271 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
blockgrid_hexa_boundaries_calculated = true;
}
void Blockgrid_coordinates::init_blockgrid_coordinates_boundary_fast()
{
blockgrid_hexa_boundary.clear();
counter = 0;
blockgrid_hexa_boundary.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
//outer loop
for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++)
{
int Nx = bg->Give_Nx_hexahedron(id_hex);
int Ny = bg->Give_Ny_hexahedron(id_hex);
int Nz = bg->Give_Nz_hexahedron(id_hex);
blockgrid_hexa_boundary.at(id_hex).resize((Nx+1)*(Ny+1)*(Nz+1));
for(int k=0;k<=Nz;k+=Nz/2) for(int j=0;j<=Ny;j+=Ny/2) for(int i=0;i<=Nx;i+=Nx/2)
{
bool A = (i == 0 || i == Nx);
bool B = (j == 0 || j == Ny);
bool C = (k == 0 || k == Nz);
if ((A&!B&!C)|(!A&B&!C)|(!A&!B&C))
{
//inner loop
for (int id_hexInner = 0 ; id_hexInner < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hexInner++)
{
if (id_hex != id_hexInner)
{
//cout << "id_hexInner " << id_hexInner<< endl;
int NxInner = bg->Give_Nx_hexahedron(id_hexInner);
int NyInner = bg->Give_Ny_hexahedron(id_hexInner);
int NzInner = bg->Give_Nz_hexahedron(id_hexInner);
for(int kInner=0;kInner<=NzInner;kInner+=NzInner/2) for(int jInner=0;jInner<=NyInner;jInner+=NyInner/2) for(int iInner=0;iInner<=NxInner;iInner+=NxInner/2)
{
bool AI = (iInner == 0 || iInner == NxInner);
bool BI = (jInner == 0 || jInner == NyInner);
bool CI = (kInner == 0 || kInner == NzInner);
if ((AI&!BI&!CI)|(!AI&BI&!CI)|(!AI&!BI&CI))
{
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner)) < 1e-10 )
{
//add whole plane :
addBoundaryPlane(id_hex,i,j,k,id_hexInner,iInner,jInner,kInner, Nx, Ny, Nz, NxInner, NyInner, NzInner);
// int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
// counter++;
// blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
// blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner);
// blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner);
// blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner);
}
}
}
}
}
}
}
}
blockgrid_hexa_boundaries_calculated = true;
}
void Blockgrid_coordinates::addBoundaryPlane(int id_hex, int i, int j, int k, int id_hexInner, int iInner, int jInner, int kInner, int Nx, int Ny, int Nz, int NxInner, int NyInner, int NzInner)
{
int PlaneHits{0};
int tempTries{0};
if (i == 0 || i == Nx)
{
//vary j,k
for (j = 0 ; j <= Ny;j++)for (k = 0 ; k <= Nz;k++)
{
if (iInner == 0 || iInner == NxInner)
{
//vary j,k
for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else if (jInner == 0 || jInner == NyInner)
{
//vary i,k
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else
{
//vary i,j
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)
{
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 )
{
tempTries++;
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner);
}
}
}
}
}
else if (j == 0 || j == Ny)
{
//vary i,k
for (i = 0 ; i <= Nx;i++)for (k = 0 ; k <= Nz;k++)
{
if (iInner == 0 || iInner == NxInner)
{
//vary j,k
for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else if (jInner == 0 || jInner == NyInner)
{
//vary i,k
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else
{
//vary i,j
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner);
}
}
}
}
}
else
{
//vary i,j
for (i = 0 ; i <= Nx;i++)for (j = 0 ; j <= Ny;j++)
{
if (iInner == 0 || iInner == NxInner)
{
//vary j,k
for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInnerT,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else if (jInner == 0 || jInner == NyInner)
{
//vary i,k
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int kInnerT = 0 ; kInnerT <= NzInner;kInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInner,kInnerT)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInnerT);
}
}
}
else
{
//vary i,j
for (int iInnerT = 0 ; iInnerT <= NxInner;iInnerT++)for (int jInnerT = 0 ; jInnerT <= NyInner;jInnerT++)
{
tempTries++;
if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInnerT,jInnerT,kInner)) < 1e-10 )
{
int convertedLocalIndex = i +(Nx+1)*(j +(Ny+1)* k);
counter++;
PlaneHits++;
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(id_hexInner);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInnerT);
blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner);
}
}
}
}
}
// std::cout << "tempTries " << tempTries << std::endl;
// std::cout << "PlaneHits " << PlaneHits << std::endl;
// std::cout << "PlaneHits " << PlaneHits << std::endl;
}
void Blockgrid_coordinates::delete_blockgrid_coordinates()
{
......@@ -1149,3 +1415,4 @@ void Blockgrid_coordinates::delete_blockgrid_coordinates()
blockgrid_hexa_boundary.clear();
}
program/source/grid/blockgrid.h
View file @ 8953fc87
......@@ -48,13 +48,22 @@ class Blockgrid_coordinates {
* berechnet indizies der randkoordinaten, sodass direkt die benachabrten konten am Rand von hexahedrons angesprochen werden können.
* */
void init_blockgrid_coordinates_boundary();
void init_blockgrid_coordinates_boundary_fast();
void addBoundaryPlane(int id_hex, int i, int j, int k, int id_hexInner, int iInner, int jInner, int kInner,int Nx, int Ny, int Nz,int NxInner,int NyInner, int NzInner);
void delete_blockgrid_coordinates();
std::vector<std::vector<std::vector<int> > > getHexaCoordiantesBoundary(){return blockgrid_hexa_boundary;}
std::vector<std::vector<std::vector<int> > > blockgrid_hexa_boundary;
std::vector< std::vector< D3vector > > * getBlockgridHexaCoordinates(){return &blockgrid_hexa_coordinates;}
std::vector< std::vector< D3vector > > * getBlockgridEdgeCoordinates(){return &blockgrid_edge_coordinates;}
void setBlockgridHexaCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_hexa_coordinates = coords;}
void setBlockgridEdgeCoordinates(std::vector< std::vector< D3vector > > coords){blockgrid_edge_coordinates = coords;}
std::vector< std::vector< D3vector > > blockgrid_hexa_coordinates;
std::vector< std::vector< D3vector > > blockgrid_edge_coordinates;
Blockgrid *bg;
bool blockgrid_edge_coordinates_calculated;
......@@ -62,6 +71,8 @@ class Blockgrid_coordinates {
bool blockgrid_hexa_boundaries_calculated;
std::vector< std::vector< bool > > blockgrid_edge_coordinates_set;
int counter;
};
......
program/source/interpol/interpol.cc
View file @ 8953fc87
......@@ -1531,11 +1531,6 @@ void Interpolate_direct::init()
correctRelation.at(5) = calculateNeighbourIndexRelation(indexAllFacesInside.at(5), indexAllFacesOutside.at(5));
}
//init surrounding boxes!
for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j) for(int k=0;k<Nz;++k)
{
......@@ -1685,6 +1680,54 @@ void Interpolate_direct::init()
}
void Interpolate_direct::updateBoundaryBoxes()
{
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);
for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++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;
}
}
}
std::vector<std::vector<int> > Interpolate_direct::calculateNeighbourIndexRelation(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer)
{
// 1 : find index which does not change --> not part of the boundary
......
program/source/interpol/interpol.h
View file @ 8953fc87
......@@ -349,6 +349,7 @@ class Interpolate_direct {
* also calculates the neighbour cell, if in other blockgrid.
**/
void init();
void updateBoundaryBoxes();
/**
* Iterates throuh all cells : if D3vector v is inside of a cell, the weighting of the 8 cell points are calculated (D3vector lambda).
......@@ -374,6 +375,10 @@ class Interpolate_direct {
DTyp evaluate(Variable<DTyp>& u);
std::vector<std::vector<std::vector<int> > > getBoundaryBox(){return array_box_boundary;}
std::vector<std::vector<std::vector<D3vector> > > * getArrayBox() {return &arrayBoxWSDENT;}
void setArrayBox(std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT_){arrayBoxWSDENT = arrayBoxWSDENT_;}
int counterFast{}, counterFastest{}, counterSamePoint{}, counterSecondTry{}, counterThirdTry{}, counterSlow{}, counterHexa{}, counterCorner{}, counterEdge{};
int checkCounter{};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment