Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PABLO: rework the function that evaluates the octant owner of a point #452

Merged
merged 4 commits into from
May 16, 2024
Merged

PABLO: rework the function that evaluates the octant owner of a point #452

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e625106
patchkernel: add function to get a constant VTK object
andrea-iob Feb 22, 2024
a097a30
PABLO: add a function to evaluate the Morton number of the anchor of …
andrea-iob Feb 22, 2024
ce02fa9
PABLO: simpilfy the function that evaluates the rank of the octant th…
andrea-iob Feb 22, 2024
a4ab652
PABLO: rework the function that evaluates the octant owner of a point
andrea-iob Feb 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
312 changes: 86 additions & 226 deletions src/PABLO/ParaTree.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3041,76 +3041,27 @@ namespace bitpit {
*/
uint32_t
ParaTree::getPointOwnerIdx(const double * point) const {
uint32_t noctants = m_octree.m_octants.size();
if(noctants==0)
return numeric_limits<uint32_t>::max();
uint32_t idxtry = noctants/2;
uint32_t x, y, z;
uint64_t morton;
int powner = 0;
//ParaTree works in [0,1] domain
if (point[0] > 1+m_tol || point[1] > 1+m_tol || point[2] > 1+m_tol
|| point[0] < -m_tol || point[1] < -m_tol || point[2] < -m_tol){
// Evaluate the Morton associated with the point
uint64_t morton = evalPointAnchorMorton(point);
if (morton == PABLO::INVALID_MORTON) {
return numeric_limits<uint32_t>::max();
}

x = m_trans.mapX(std::min(std::max(point[0], 0.0), 1.0));
y = m_trans.mapY(std::min(std::max(point[1], 0.0), 1.0));
z = m_trans.mapZ(std::min(std::max(point[2], 0.0), 1.0));

uint32_t maxLength = getMaxLength();
if (x == maxLength) x = x - 1;
if (y == maxLength) y = y - 1;
if (z == maxLength) z = z - 1;
morton = PABLO::computeMorton(m_dim,x,y,z);


powner = 0;
if(!m_serial) powner = findOwner(morton);

if ((powner!=m_rank) && (!m_serial))
// Early return if the point belongs to a ghost
if (findOwner(morton) != m_rank && (!m_serial)) {
return numeric_limits<uint32_t>::max();

int32_t jump = idxtry;
while(abs(jump) > 0){

uint64_t mortontry = m_octree.m_octants[idxtry].getMorton();
jump = ((mortontry<morton)-(mortontry>morton))*abs(jump)/2;
idxtry += jump;
if (idxtry > noctants-1){
if (jump > 0){
idxtry = noctants - 1;
jump = 0;
}
else if (jump < 0){
idxtry = 0;
jump = 0;
}
}
}
if(m_octree.m_octants[idxtry].getMorton() == morton){
return idxtry;
}
else{
// Step until the mortontry lower than morton (one idx of distance)
{
while(m_octree.m_octants[idxtry].getMorton() < morton){
idxtry++;
if(idxtry > noctants-1){
idxtry = noctants-1;
break;
}
}
while(m_octree.m_octants[idxtry].getMorton() > morton){
idxtry--;
if(idxtry > noctants-1){
idxtry = 0;
break;
}
}
}
return idxtry;
}

// Identify the octant that contains the point
//
// If the owner of the point is an internal octant, decrementing the Morton upper bound
// by one will give us the requested octant.
uint32_t pointOwnerIdx;
uint64_t pointOwnerMorton;
m_octree.findMortonUpperBound(morton, m_octree.m_octants, &pointOwnerIdx, &pointOwnerMorton);
--pointOwnerIdx;

return pointOwnerIdx;
};

/** Get the octant owner of an input point.
Expand All @@ -3120,193 +3071,102 @@ namespace bitpit {
*/
uint32_t
ParaTree::getPointOwnerIdx(const double * point, bool & isghost) const {
uint32_t noctants = m_octree.m_octants.size();
if(noctants==0)
return numeric_limits<uint32_t>::max();
uint32_t idxtry = noctants/2;
uint32_t x, y, z;
uint64_t morton, mortontry;
int powner = 0;
isghost = false;
//ParaTree works in [0,1] domain
if (point[0] > 1+m_tol || point[1] > 1+m_tol || point[2] > 1+m_tol
|| point[0] < -m_tol || point[1] < -m_tol || point[2] < -m_tol){
uint64_t morton = evalPointAnchorMorton(point);
if (morton == PABLO::INVALID_MORTON) {
isghost = false;

return numeric_limits<uint32_t>::max();
}

x = m_trans.mapX(std::min(std::max(point[0], 0.0), 1.0));
y = m_trans.mapY(std::min(std::max(point[1], 0.0), 1.0));
z = m_trans.mapZ(std::min(std::max(point[2], 0.0), 1.0));

uint32_t maxLength = getMaxLength();
if (x == maxLength) x = x - 1;
if (y == maxLength) y = y - 1;
if (z == maxLength) z = z - 1;
morton = PABLO::computeMorton(m_dim,x,y,z);

// Check if the point belongs to a ghost
isghost = (findOwner(morton) != m_rank);

powner = 0;
if(!m_serial) powner = findOwner(morton);

if (powner==m_rank){

int32_t jump = idxtry;
while(abs(jump) > 0){

mortontry = m_octree.m_octants[idxtry].getMorton();
jump = ((mortontry<morton)-(mortontry>morton))*abs(jump)/2;
idxtry += jump;
if (idxtry > noctants-1){
if (jump > 0){
idxtry = noctants - 1;
jump = 0;
}
else if (jump < 0){
idxtry = 0;
jump = 0;
// Identify the octant that contains the point
//
// If the owner of the point is an internal octant, decrementing the Morton upper bound
// by one will give us the requested octant. If the owner of the point is a ghost octant,
// we need to test all the octants before the Morton upper bound.
uint32_t pointOwnerIdx;
uint64_t pointOwnerMorton;
if (!isghost) {
m_octree.findMortonUpperBound(morton, m_octree.m_octants, &pointOwnerIdx, &pointOwnerMorton);
--pointOwnerIdx;
} else {
m_octree.findMortonUpperBound(morton, m_octree.m_ghosts, &pointOwnerIdx, &pointOwnerMorton);
while (pointOwnerIdx > 0) {
--pointOwnerIdx;

const Octant *octant = getGhostOctant(pointOwnerIdx);
darray3 octantAnchor = getCoordinates(octant);
double octantSize = getSize(octant);

bool found = true;
for (int i = 0; i < m_dim; ++i) {
if (point[i] < octantAnchor[i] - m_tol || point[i] > (octantAnchor[i] + octantSize + m_tol)) {
found = false;
break;
}
}
}
if(m_octree.m_octants[idxtry].getMorton() == morton){
return idxtry;
}
else{
// Step until the mortontry lower than morton (one idx of distance)
{
while(m_octree.m_octants[idxtry].getMorton() < morton){
idxtry++;
if(idxtry > noctants-1){
idxtry = noctants-1;
break;
}
}
while(m_octree.m_octants[idxtry].getMorton() > morton){
idxtry--;
if(idxtry > noctants-1){
idxtry = 0;
break;
}
}

if (found) {
return pointOwnerIdx;
}
return idxtry;
}
}
else if((powner != m_rank) && m_serial){

return numeric_limits<uint32_t>::max();
}
else{
//GHOST SEARCH
uint32_t nghosts = m_octree.m_ghosts.size();
idxtry = nghosts/2;
int32_t jump = idxtry;
while(abs(jump) > 0){

mortontry = m_octree.m_ghosts[idxtry].getMorton();
jump = ((mortontry<morton)-(mortontry>morton))*abs(jump)/2;
idxtry += jump;
if (idxtry > nghosts-1){
if (jump > 0){
idxtry = nghosts - 1;
jump = 0;
}
else if (jump < 0){
idxtry = 0;
jump = 0;
}
}
}
if(m_octree.m_ghosts[idxtry].getMorton() == morton){
isghost = true;
return idxtry;
}
else{
// Step until the mortontry lower than morton (one idx of distance)
{
while(m_octree.m_ghosts[idxtry].getMorton() < morton){
idxtry++;
if(idxtry > nghosts-1){
idxtry = nghosts-1;
break;
}
}
while(m_octree.m_ghosts[idxtry].getMorton() > morton){
idxtry--;
if(idxtry > nghosts-1){
idxtry = 0;
break;
}
}
}

const Octant* octtry = getGhostOctant(idxtry);
dvector anchor_idxtry = {{getX(octtry),getY(octtry),getZ(octtry)}};
double size_try = getSize(octtry);
bool isInIdxtry = true;

for(int i = 0; i < m_dim; ++i){
isInIdxtry = isInIdxtry && (point[i] >= anchor_idxtry[i] && point[i] <= (anchor_idxtry[i] + size_try));
}

if( isInIdxtry){
isghost = true;
return idxtry;
}
else{
return numeric_limits<uint32_t>::max();
}
}
}///end ghosts search
return pointOwnerIdx;
};

/** Get the octant owner rank of an input point.
/** Get the rank of the octant that contains the specified point.
* \param[in] point Coordinates of target point.
* \return Owner rank of target point (negative if out of global domain).
*/
int
ParaTree::getPointOwnerRank(darray3 point){

uint32_t x,y,z;
uint64_t morton;

if (point[0] > 1+m_tol || point[1] > 1+m_tol || point[2] > 1+m_tol
|| point[0] < -m_tol || point[1] < -m_tol || point[2] < -m_tol){
ParaTree::getPointOwnerRank(const darray3 &point){
// Early return if the point is not associatd with a valid Morton
uint64_t morton = evalPointAnchorMorton(point.data());
if (morton == PABLO::INVALID_MORTON) {
return -1;
}
point[0] = min(max(point[0],0.0),1.0);
point[1] = min(max(point[1],0.0),1.0);
point[2] = min(max(point[2],0.0),1.0);

x = m_trans.mapX(point[0]);
y = m_trans.mapY(point[1]);
z = m_trans.mapZ(point[2]);
// Identify that partition that contains the point
return findOwner(morton);
};

uint32_t maxLength = getMaxLength();
if ((x > maxLength) || (y > maxLength) || (z > maxLength)
|| (point[0] < m_trans.m_origin[0]) || (point[1] < m_trans.m_origin[1]) || (point[2] < m_trans.m_origin[2])){
return -1;
/** Evaluate the Morton number of the anchor associated with the specified point.
* The anchor of a point is the lower-left-back vertex of the smallest octant that
* contains the point.
* \param[in] point Coordinates of target point.
* \return The Morton number of the anchor associated with the specified point.
*/
uint64_t
ParaTree::evalPointAnchorMorton(const double * point) const {
// Early return if the tree is empty
if (m_octree.m_octants.empty()) {
return PABLO::INVALID_MORTON;
}

if (m_serial)
return m_rank;

if (x == maxLength)
x = x - 1;

if (y == maxLength)
y = y - 1;

if (z == maxLength)
z = z - 1;
// Early return if the point is outside the domain
if (point[0] > 1+m_tol || point[1] > 1+m_tol || point[2] > 1+m_tol) {
return PABLO::INVALID_MORTON;
} else if (point[0] < -m_tol || point[1] < -m_tol || point[2] < -m_tol) {
return PABLO::INVALID_MORTON;
}

morton = PABLO::computeMorton(m_dim, x, y, z);
// Evaluate the Morton associated to the point
uint32_t x = m_trans.mapX(std::min(std::max(point[0], 0.0), 1.0));
uint32_t y = m_trans.mapY(std::min(std::max(point[1], 0.0), 1.0));
uint32_t z = m_trans.mapZ(std::min(std::max(point[2], 0.0), 1.0));

for (int p = 0; p < m_nproc; ++p){
if (morton <= m_partitionLastDesc[p] && morton >= m_partitionFirstDesc[p])
return p;
}
uint32_t maxLength = getMaxLength();
if (x == maxLength) x = x - 1;
if (y == maxLength) y = y - 1;
if (z == maxLength) z = z - 1;

return -1;
};
return PABLO::computeMorton(m_dim, x, y, z);
}

/** Get the local index of the node of a target octant, corresponding to the splitting node of its family; i.e. the index of the local node
* coincident with the center point of its father.
Expand Down
4 changes: 3 additions & 1 deletion src/PABLO/ParaTree.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -469,6 +469,8 @@ namespace bitpit {

void reset(bool createRoot);

uint64_t evalPointAnchorMorton(const double * point) const;

// =================================================================================== //
// OTHER OCTANT BASED METHODS //
// =================================================================================== //
Expand Down Expand Up @@ -503,7 +505,7 @@ namespace bitpit {
bool isNodeOnOctant(const Octant* nodeOctant, uint8_t nodeIndex, const Octant* octant) const;
bool isEdgeOnOctant(const Octant* edgeOctant, uint8_t edgeIndex, const Octant* octant) const;
bool isFaceOnOctant(const Octant* faceOctant, uint8_t faceIndex, const Octant* octant) const;
int getPointOwnerRank(darray3 point);
int getPointOwnerRank(const darray3 &point);
uint8_t getFamilySplittingNode(const Octant*) const;
void expectedOctantAdapt(const Octant* octant, int8_t marker, octvector* result) const;

Expand Down
10 changes: 10 additions & 0 deletions src/patchkernel/patch_kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7811,6 +7811,16 @@ VTKUnstructuredGrid & PatchKernel::getVTK()
return m_vtk;
}

/*!
Get the VTK object.

\result The VTK object.
*/
const VTKUnstructuredGrid & PatchKernel::getVTK() const
{
return m_vtk;
}

/*!
Get the VTK write target.

Expand Down
1 change: 1 addition & 0 deletions src/patchkernel/patch_kernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -678,6 +678,7 @@ friend class PatchManager;
void displayInterfaces(std::ostream &out, unsigned int padding = 0) const;

VTKUnstructuredGrid & getVTK();
const VTKUnstructuredGrid & getVTK() const;
WriteTarget getVTKWriteTarget() const;
void setVTKWriteTarget(WriteTarget targetCells);
const CellConstRange getVTKCellWriteRange() const;
Expand Down
Loading