diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index 07dd39f93f..25a4f8fbce 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/src/OpenFOAM/meshes/meshShapes/face/face.H @@ -222,16 +222,16 @@ public: ) const; //- Fast intersection with a ray. - // For a hit, the pointHit.distance() is the line parameter t : - // intersection=p+t*q. Only defined for FULL_RAY or - // HALF_RAY. + // Does face-center decomposition and returns triangle intersection + // point closest to p. See triangle::intersection for details. pointHit intersection ( const point& p, const vector& q, const point& ctr, const pointField& meshPoints, - const intersection::algorithm alg + const intersection::algorithm alg, + const scalar tol = 0.0 ) const; //- Return nearest point to face diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C index b514af7749..c5f9a80250 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C @@ -136,7 +136,8 @@ Foam::pointHit Foam::face::intersection const vector& q, const point& ctr, const pointField& meshPoints, - const intersection::algorithm alg + const intersection::algorithm alg, + const scalar tol ) const { scalar nearestHitDist = VGREAT; @@ -154,7 +155,7 @@ Foam::pointHit Foam::face::intersection meshPoints[f[pI]], meshPoints[f[fcIndex(pI)]], ctr - ).intersection(p, q, alg); + ).intersection(p, q, alg, tol); if (curHit.hit()) { diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index ecc04cbb8b..4d7e8d5549 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -153,7 +153,7 @@ public: inline scalar sweptVol(const triangle& t) const; //- Return point intersection with a ray. - // For a hit, the distance is signed. Positive number + // For a hit, the distance is signed. Positive number // represents the point in front of triangle. // In case of miss pointHit is set to nearest point // on triangle and its distance to the distance between @@ -166,12 +166,14 @@ public: const intersection::direction dir = intersection::VECTOR ) const; - //- Fast intersection with a ray. - // For a hit, the pointHit.distance() is the line parameter t : - // intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or + //- Fast intersection with a ray. Distance is normal distance + // to the intersection. + // For a hit, the distance is signed. Positive number + // represents the point in front of triangle. In case of miss + // pointHit position is undefined. + // Only defined for VISIBLE, FULL_RAY or // HALF_RAY. tol increases the virtual size of the triangle - // by a relative factor - can be used to compensate for - // limited precision. + // by a relative factor. inline pointHit intersection ( const point& p, diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 876389057e..33b3cab15e 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -585,7 +585,7 @@ inline bool triangle::classify bool hit = false; - if (Foam::mag(u1) < SMALL) + if (Foam::mag(u1) < ROOTVSMALL) { beta = u0/u2; diff --git a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C index bcf92e9747..2f58e5cc03 100644 --- a/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C +++ b/src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C @@ -73,7 +73,7 @@ Foam::scalar Foam::ExactParticle::trackToFace const labelList& cFaces = mesh.cells()[this->celli_]; point intersection(vector::zero); - scalar trackFraction = VGREAT; + scalar distanceSqr = VGREAT; label hitFacei = -1; const vector vec = endPosition-this->position_; @@ -93,32 +93,21 @@ Foam::scalar Foam::ExactParticle::trackToFace intersection::HALF_RAY ); - if (inter.hit() && inter.distance() < trackFraction) + if (inter.hit()) { - trackFraction = inter.distance(); - hitFacei = facei; - intersection = inter.hitPoint(); + scalar s = magSqr(inter.hitPoint()-this->position_); + + if (s < distanceSqr) + { + distanceSqr = s; + hitFacei = facei; + intersection = inter.hitPoint(); + } } } } - - if (hitFacei != -1) - { - if (trackFraction > 1.0) - { - // Nearest intersection beyond endPosition so we hit endPosition. - this->position_ = endPosition; - this->facei_ = -1; - return 1.0; - } - else - { - this->position_ = intersection; - this->facei_ = hitFacei; - } - } - else + if (hitFacei == -1) { // Did not find any intersection. Fall back to original approximate // algorithm @@ -130,7 +119,24 @@ Foam::scalar Foam::ExactParticle::trackToFace } - // Normal situation (trackFraction 0..1) + scalar trackFraction = Foam::sqrt(distanceSqr/magSqr(vec)); + + if (trackFraction >= (1.0-SMALL)) + { + // Nearest intersection beyond endPosition so we hit endPosition. + this->position_ = endPosition; + this->facei_ = -1; + return 1.0; + } + else + { + this->position_ = intersection; + this->facei_ = hitFacei; + } + + + // Normal situation (trackFraction 0..1). Straight copy + // of Particle::trackToFace. bool internalFace = this->cloud().internalFace(this->facei_); diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C index 2131c3aaf6..4f55564999 100644 --- a/src/meshTools/indexedOctree/indexedOctree.C +++ b/src/meshTools/indexedOctree/indexedOctree.C @@ -29,6 +29,12 @@ License #include "meshTools.H" #include "OFstream.H" +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +Foam::scalar Foam::indexedOctree::perturbTol_ = 10*SMALL; + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Does bb intersect a sphere around sample? Or is any corner point of bb @@ -203,7 +209,7 @@ Foam::indexedOctree::divide || bb.min()[2] >= bb.max()[2] ) { - FatalErrorIn("indexedOctree::divide") + FatalErrorIn("indexedOctree::divide(..)") << "Badly formed bounding box:" << bb << abort(FatalError); } @@ -674,61 +680,768 @@ void Foam::indexedOctree::findNearest } +template +Foam::treeBoundBox Foam::indexedOctree::subBbox +( + const label parentNodeI, + const direction octant +) const +{ + // Get type of node at octant + const node& nod = nodes_[parentNodeI]; + labelBits index = nod.subNodes_[octant]; + + if (isNode(index)) + { + // Use stored bb + return nodes_[getNode(index)].bb_; + } + else + { + // Calculate subBb + return nod.bb_.subBbox(octant); + } +} + + +// Takes a bb and a point on/close to the edge of the bb and pushes the point +// inside by a small fraction. +template +Foam::point Foam::indexedOctree::pushPoint +( + const treeBoundBox& bb, + const point& pt, + const bool pushInside +) +{ + // Get local length scale. + const vector perturbVec = perturbTol_*(bb.span()); + + point perturbedPt(pt); + + // Modify all components which are close to any face of the bb to be + // well inside/outside them. + + if (pushInside) + { + for (direction dir = 0; dir < vector::nComponents; dir++) + { + if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir])) + { + // Close to 'left' side. Push well beyond left side. + scalar perturbDist = perturbVec[dir] + ROOTVSMALL; + perturbedPt[dir] = bb.min()[dir] + perturbDist; + } + else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir])) + { + // Close to 'right' side. Push well beyond right side. + scalar perturbDist = perturbVec[dir] + ROOTVSMALL; + perturbedPt[dir] = bb.max()[dir] - perturbDist; + } + } + } + else + { + for (direction dir = 0; dir < vector::nComponents; dir++) + { + if (mag(pt[dir]-bb.min()[dir]) < mag(perturbVec[dir])) + { + scalar perturbDist = perturbVec[dir] + ROOTVSMALL; + perturbedPt[dir] = bb.min()[dir] - perturbDist; + } + else if (mag(pt[dir]-bb.max()[dir]) < mag(perturbVec[dir])) + { + scalar perturbDist = perturbVec[dir] + ROOTVSMALL; + perturbedPt[dir] = bb.max()[dir] + perturbDist; + } + } + } + + if (debug) + { + if (pushInside != bb.contains(perturbedPt)) + { + FatalErrorIn("indexedOctree::pushPoint(..)") + << "pushed point:" << pt + << " to:" << perturbedPt + << " wanted side:" << pushInside + << " obtained side:" << bb.contains(perturbedPt) + << " of bb:" << bb + << abort(FatalError); + } + } + + return perturbedPt; +} + + +// Takes a bb and a point on the edge of the bb and pushes the point +// outside by a small fraction. +template +Foam::point Foam::indexedOctree::pushPoint +( + const treeBoundBox& bb, + const direction faceID, + const point& pt, + const bool pushInside +) +{ + // Get local length scale. + const vector perturbVec = perturbTol_*bb.span(); + + point perturbedPt(pt); + + // Modify all components which are close to any face of the bb to be + // well outside them. + + if (faceID == 0) + { + FatalErrorIn("indexedOctree::pushPoint(..)") + << abort(FatalError); + } + + if (faceID & treeBoundBox::LEFTBIT) + { + if (pushInside) + { + perturbedPt[0] = bb.min()[0] + (perturbVec[0] + ROOTVSMALL); + } + else + { + perturbedPt[0] = bb.min()[0] - (perturbVec[0] + ROOTVSMALL); + } + } + else if (faceID & treeBoundBox::RIGHTBIT) + { + if (pushInside) + { + perturbedPt[0] = bb.max()[0] - (perturbVec[0] + ROOTVSMALL); + } + else + { + perturbedPt[0] = bb.max()[0] + (perturbVec[0] + ROOTVSMALL); + } + } + + if (faceID & treeBoundBox::BOTTOMBIT) + { + if (pushInside) + { + perturbedPt[1] = bb.min()[1] + (perturbVec[1] + ROOTVSMALL); + } + else + { + perturbedPt[1] = bb.min()[1] - (perturbVec[1] + ROOTVSMALL); + } + } + else if (faceID & treeBoundBox::TOPBIT) + { + if (pushInside) + { + perturbedPt[1] = bb.max()[1] - (perturbVec[1] + ROOTVSMALL); + } + else + { + perturbedPt[1] = bb.max()[1] + (perturbVec[1] + ROOTVSMALL); + } + } + + if (faceID & treeBoundBox::BACKBIT) + { + if (pushInside) + { + perturbedPt[2] = bb.min()[2] + (perturbVec[2] + ROOTVSMALL); + } + else + { + perturbedPt[2] = bb.min()[2] - (perturbVec[2] + ROOTVSMALL); + } + } + else if (faceID & treeBoundBox::FRONTBIT) + { + if (pushInside) + { + perturbedPt[2] = bb.max()[2] - (perturbVec[2] + ROOTVSMALL); + } + else + { + perturbedPt[2] = bb.max()[2] + (perturbVec[2] + ROOTVSMALL); + } + } + + if (debug) + { + if (pushInside != bb.contains(perturbedPt)) + { + FatalErrorIn("indexedOctree::pushPoint(..)") + << "pushed point:" << pt << " on face:" << faceString(faceID) + << " to:" << perturbedPt + << " wanted side:" << pushInside + << " obtained side:" << bb.contains(perturbedPt) + << " of bb:" << bb + << abort(FatalError); + } + } + + return perturbedPt; +} + + +// Guarantees that if pt is on a face it gets perturbed so it is away +// from the face edges. +// If pt is not on a face does nothing. +template +Foam::point Foam::indexedOctree::pushPointIntoFace +( + const treeBoundBox& bb, + const vector& dir, // end-start + const point& pt +) +{ + if (debug) + { + if (bb.posBits(pt) != 0) + { + FatalErrorIn("indexedOctree::pushPointIntoFace(..)") + << " bb:" << bb << endl + << "does not contain point " << pt << abort(FatalError); + } + } + + + // Handle two cases: + // - point exactly on multiple faces. Push away from all but one. + // - point on a single face. Push away from edges of face. + + direction ptFaceID = bb.faceBits(pt); + + direction nFaces = 0; + FixedList faceIndices; + + if (ptFaceID & treeBoundBox::LEFTBIT) + { + faceIndices[nFaces++] = treeBoundBox::LEFT; + } + else if (ptFaceID & treeBoundBox::RIGHTBIT) + { + faceIndices[nFaces++] = treeBoundBox::RIGHT; + } + + if (ptFaceID & treeBoundBox::BOTTOMBIT) + { + faceIndices[nFaces++] = treeBoundBox::BOTTOM; + } + else if (ptFaceID & treeBoundBox::TOPBIT) + { + faceIndices[nFaces++] = treeBoundBox::TOP; + } + + if (ptFaceID & treeBoundBox::BACKBIT) + { + faceIndices[nFaces++] = treeBoundBox::BACK; + } + else if (ptFaceID & treeBoundBox::FRONTBIT) + { + faceIndices[nFaces++] = treeBoundBox::FRONT; + } + + + // Determine the face we want to keep the point on + + direction keepFaceID; + + if (nFaces == 0) + { + // Return original point + return pt; + } + else if (nFaces == 1) + { + keepFaceID = faceIndices[0]; + } + else + { + // Determine best face out of faceIndices[0 .. nFaces-1]. + // The best face is the one most perpendicular to the ray direction. + + keepFaceID = faceIndices[0]; + scalar maxInproduct = mag(treeBoundBox::faceNormals[keepFaceID] & dir); + + for (direction i = 1; i < nFaces; i++) + { + direction face = faceIndices[i]; + scalar s = mag(treeBoundBox::faceNormals[face] & dir); + if (s > maxInproduct) + { + maxInproduct = s; + keepFaceID = face; + } + } + } + + + // 1. Push point into bb, away from all corners + + point facePoint(pushPoint(bb, pt, true)); + direction faceID = 0; + + // 2. Snap it back onto the preferred face + + if (keepFaceID == treeBoundBox::LEFT) + { + facePoint.x() = bb.min().x(); + faceID = treeBoundBox::LEFTBIT; + } + else if (keepFaceID == treeBoundBox::RIGHT) + { + facePoint.x() = bb.max().x(); + faceID = treeBoundBox::RIGHTBIT; + } + else if (keepFaceID == treeBoundBox::BOTTOM) + { + facePoint.y() = bb.min().y(); + faceID = treeBoundBox::BOTTOMBIT; + } + else if (keepFaceID == treeBoundBox::TOP) + { + facePoint.y() = bb.max().y(); + faceID = treeBoundBox::TOPBIT; + } + else if (keepFaceID == treeBoundBox::BACK) + { + facePoint.z() = bb.min().z(); + faceID = treeBoundBox::BACKBIT; + } + else if (keepFaceID == treeBoundBox::FRONT) + { + facePoint.z() = bb.max().z(); + faceID = treeBoundBox::FRONTBIT; + } + + + if (debug) + { + if (faceID != bb.faceBits(facePoint)) + { + FatalErrorIn("indexedOctree::pushPointIntoFace(..)") + << "Pushed point from " << pt + << " on face:" << ptFaceID << " of bb:" << bb << endl + << "onto " << facePoint + << " on face:" << faceID + << " which is not consistent with geometric face " + << bb.faceBits(facePoint) + << abort(FatalError); + } + if (bb.posBits(facePoint) != 0) + { + FatalErrorIn("indexedOctree::pushPointIntoFace(..)") + << " bb:" << bb << endl + << "does not contain perturbed point " + << facePoint << abort(FatalError); + } + } + + + return facePoint; +} + + +//// Takes a bb and a point on the outside of the bb. Checks if on multiple faces +//// and if so perturbs point so it is only on one face. +//template +//void Foam::indexedOctree::checkMultipleFaces +//( +// const treeBoundBox& bb, +// const vector& dir, // end-start +// pointIndexHit& faceHitInfo, +// direction& faceID +//) +//{ +// // Do the quick elimination of no or one face. +// if +// ( +// (faceID == 0) +// || (faceID == treeBoundBox::LEFTBIT) +// || (faceID == treeBoundBox::RIGHTBIT) +// || (faceID == treeBoundBox::BOTTOMBIT) +// || (faceID == treeBoundBox::TOPBIT) +// || (faceID == treeBoundBox::BACKBIT) +// || (faceID == treeBoundBox::FRONTBIT) +// ) +// { +// return; +// } +// +// +// // Check the direction of vector w.r.t. faces being intersected. +// FixedList inproducts(-GREAT); +// +// direction nFaces = 0; +// +// if (faceID & treeBoundBox::LEFTBIT) +// { +// inproducts[treeBoundBox::LEFT] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::LEFT] +// & dir +// ); +// nFaces++; +// } +// if (faceID & treeBoundBox::RIGHTBIT) +// { +// inproducts[treeBoundBox::RIGHT] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::RIGHT] +// & dir +// ); +// nFaces++; +// } +// +// if (faceID & treeBoundBox::BOTTOMBIT) +// { +// inproducts[treeBoundBox::BOTTOM] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::BOTTOM] +// & dir +// ); +// nFaces++; +// } +// if (faceID & treeBoundBox::TOPBIT) +// { +// inproducts[treeBoundBox::TOP] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::TOP] +// & dir +// ); +// nFaces++; +// } +// +// if (faceID & treeBoundBox::BACKBIT) +// { +// inproducts[treeBoundBox::BACK] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::BACK] +// & dir +// ); +// nFaces++; +// } +// if (faceID & treeBoundBox::FRONTBIT) +// { +// inproducts[treeBoundBox::FRONT] = mag +// ( +// treeBoundBox::faceNormals[treeBoundBox::FRONT] +// & dir +// ); +// nFaces++; +// } +// +// if (nFaces == 0 || nFaces == 1 || nFaces > 3) +// { +// FatalErrorIn("indexedOctree::checkMultipleFaces(..)") +// << "Problem : nFaces:" << nFaces << abort(FatalError); +// } +// +// // Keep point on most perpendicular face; shift it away from the aligned +// // ones. +// // E.g. line hits top and left face: +// // a +// // ----+----+ +// // | | +// // | | +// // +----+ +// // Shift point down (away from top): +// // +// // a+----+ +// // ----| | +// // | | +// // +----+ +// +// label maxIndex = -1; +// scalar maxInproduct = -GREAT; +// +// for (direction i = 0; i < 6; i++) +// { +// if (inproducts[i] > maxInproduct) +// { +// maxInproduct = inproducts[i]; +// maxIndex = i; +// } +// } +// +// if (maxIndex == -1) +// { +// FatalErrorIn("indexedOctree::checkMultipleFaces(..)") +// << "Problem maxIndex:" << maxIndex << " inproducts:" << inproducts +// << abort(FatalError); +// } +// +// const point oldPoint(faceHitInfo.rawPoint()); +// const direction oldFaceID = faceID; +// +// // 1. Push point into bb, away from all corners +// +// faceHitInfo.rawPoint() = pushPoint(bb, oldFaceID, oldPoint, true); +// +// // 2. Snap it back onto the preferred face +// +// if (maxIndex == treeBoundBox::LEFT) +// { +// faceHitInfo.rawPoint().x() = bb.min().x(); +// faceID = treeBoundBox::LEFTBIT; +// } +// else if (maxIndex == treeBoundBox::RIGHT) +// { +// faceHitInfo.rawPoint().x() = bb.max().x(); +// faceID = treeBoundBox::RIGHTBIT; +// } +// else if (maxIndex == treeBoundBox::BOTTOM) +// { +// faceHitInfo.rawPoint().y() = bb.min().y(); +// faceID = treeBoundBox::BOTTOMBIT; +// } +// else if (maxIndex == treeBoundBox::TOP) +// { +// faceHitInfo.rawPoint().y() = bb.max().y(); +// faceID = treeBoundBox::TOPBIT; +// } +// else if (maxIndex == treeBoundBox::BACK) +// { +// faceHitInfo.rawPoint().z() = bb.min().z(); +// faceID = treeBoundBox::BACKBIT; +// } +// else if (maxIndex == treeBoundBox::FRONT) +// { +// faceHitInfo.rawPoint().z() = bb.max().z(); +// faceID = treeBoundBox::FRONTBIT; +// } +// +// Pout<< "From ray:" << dir +// << " from point:" << oldPoint +// << " on faces:" << faceString(oldFaceID) +// << " of bb:" << bb +// << " with inprods:" << inproducts +// << " maxIndex:" << maxIndex << endl +// << "perturbed to point:" << faceHitInfo.rawPoint() +// << " on face:" << faceString(faceID) +// << endl; +// +// +// if (debug) +// { +// if (faceID != bb.faceBits(faceHitInfo.rawPoint())) +// { +// FatalErrorIn("indexedOctree::checkMultipleFaces(..)") +// << "Pushed point from " << oldPoint +// << " on face:" << oldFaceID << " of bb:" << bb << endl +// << "onto " << faceHitInfo.rawPoint() +// << " on face:" << faceID +// << " which is not consistent with geometric face " +// << bb.faceBits(faceHitInfo.rawPoint()) +// << abort(FatalError); +// } +// } +//} + + +// Get parent node and octant. Return false if top of tree reached. +template +bool Foam::indexedOctree::walkToParent +( + const label nodeI, + const direction octant, + + label& parentNodeI, + label& parentOctant +) const +{ + parentNodeI = nodes_[nodeI].parent_; + + if (parentNodeI == -1) + { + // Reached edge of tree + return false; + } + + const node& parentNode = nodes_[parentNodeI]; + + // Find octant nodeI is in. + parentOctant = 255; + + for (direction i = 0; i < parentNode.subNodes_.size(); i++) + { + labelBits index = parentNode.subNodes_[i]; + + if (isNode(index) && getNode(index) == nodeI) + { + parentOctant = i; + break; + } + } + + if (parentOctant == 255) + { + FatalErrorIn("walkToParent(..)") + << "Problem: no parent found for octant:" << octant + << " node:" << nodeI + << abort(FatalError); + } + + return true; +} + + // Walk tree to neighbouring node. Gets current position as // node and octant in this node and walks in the direction given by -// the faceID (one of treeBoundBox::LEFTBIT, RIGHTBIT etc.) +// the facePointBits (combination of treeBoundBox::LEFTBIT, TOPBIT etc.) // Returns false if edge of tree hit. template bool Foam::indexedOctree::walkToNeighbour ( const point& facePoint, - const direction faceID, // direction to walk in + const direction faceID, // face(s) that facePoint is on label& nodeI, direction& octant ) const { + label oldNodeI = nodeI; + direction oldOctant = octant; + // Find out how to test for octant. Say if we want to go left we need // to traverse up the tree until we hit a node where our octant is // on the right. + // Coordinate direction to test + const direction X = treeBoundBox::RIGHTHALF; + const direction Y = treeBoundBox::TOPHALF; + const direction Z = treeBoundBox::FRONTHALF; + direction octantMask = 0; - direction valueMask = 0; + direction wantedValue = 0; if ((faceID & treeBoundBox::LEFTBIT) != 0) { - // We want to go left so check if in right octant. - octantMask |= treeBoundBox::RIGHTHALF; - valueMask |= treeBoundBox::RIGHTHALF; + // We want to go left so check if in right octant (i.e. x-bit is set) + octantMask |= X; + wantedValue |= X; } else if ((faceID & treeBoundBox::RIGHTBIT) != 0) { - octantMask |= treeBoundBox::RIGHTHALF; // valueMask already 0 + octantMask |= X; // wantedValue already 0 } if ((faceID & treeBoundBox::BOTTOMBIT) != 0) { - octantMask |= treeBoundBox::TOPHALF; - valueMask |= treeBoundBox::TOPHALF; + // Want to go down so check for y-bit set. + octantMask |= Y; + wantedValue |= Y; } else if ((faceID & treeBoundBox::TOPBIT) != 0) { - octantMask |= treeBoundBox::TOPHALF; + // Want to go up so check for y-bit not set. + octantMask |= Y; } if ((faceID & treeBoundBox::BACKBIT) != 0) { - octantMask |= treeBoundBox::FRONTHALF; - valueMask |= treeBoundBox::FRONTHALF; + octantMask |= Z; + wantedValue |= Z; } else if ((faceID & treeBoundBox::FRONTBIT) != 0) { - octantMask |= treeBoundBox::FRONTHALF; + octantMask |= Z; } + // So now we have the coordinate directions in the octant we need to check + // and the resulting value. + + /* + // +---+---+ + // | | | + // | | | + // | | | + // +---+-+-+ + // | | | | + // | a+-+-+ + // | |\| | + // +---+-+-+ + // \ + // + // e.g. ray is at (a) in octant 0(or 4) with faceIDs : LEFTBIT+TOPBIT. + // If we would be in octant 1(or 5) we could go to the correct octant + // in the same node by just flipping the x and y bits (exoring). + // But if we are not in octant 1/5 we have to go up until we are. + // In general for leftbit+topbit: + // - we have to check for x and y : octantMask = 011 + // - the checked bits have to be : wantedValue = ?01 + */ + + //Pout<< "For point " << facePoint << endl; + // Go up until we have chance to cross to the wanted direction - while (valueMask != (octant & octantMask)) + while (wantedValue != (octant & octantMask)) { - label parentNodeI = nodes_[nodeI].parent_; + // Go up to the parent. + + // Remove the directions that are not on the boundary of the parent. + // See diagram above + if (wantedValue & X) // && octantMask&X + { + // Looking for right octant. + if (octant & X) + { + // My octant is one of the left ones so punch out the x check + octantMask &= ~X; + wantedValue &= ~X; + } + } + else + { + if (!(octant & X)) + { + // My octant is right but I want to go left. + octantMask &= ~X; + wantedValue &= ~X; + } + } + + if (wantedValue & Y) + { + if (octant & Y) + { + octantMask &= ~Y; + wantedValue &= ~Y; + } + } + else + { + if (!(octant & Y)) + { + octantMask &= ~Y; + wantedValue &= ~Y; + } + } + + if (wantedValue & Z) + { + if (octant & Z) + { + octantMask &= ~Z; + wantedValue &= ~Z; + } + } + else + { + if (!(octant & Z)) + { + octantMask &= ~Z; + wantedValue &= ~Z; + } + } + + + label parentNodeI; + label parentOctant; + walkToParent(nodeI, octant, parentNodeI, parentOctant); if (parentNodeI == -1) { @@ -736,38 +1449,46 @@ bool Foam::indexedOctree::walkToNeighbour return false; } - const node& parentNode = nodes_[parentNodeI]; + //Pout<< " walked from node:" << nodeI << " octant:" << octant + // << " bb:" << nodes_[nodeI].bb_.subBbox(octant) << endl + // << " to:" << parentNodeI << " octant:" << parentOctant + // << " bb:" << nodes_[parentNodeI].bb_.subBbox(parentOctant) + // << endl; + // + //Pout<< " octantMask:" << octantMask + // << " wantedValue:" << wantedValue << endl; - // Find octant nodeI is in. - direction parentOctant = 255; - - for (direction i = 0; i < parentNode.subNodes_.size(); i++) - { - labelBits index = parentNode.subNodes_[i]; - - if (isNode(index) && getNode(index) == nodeI) - { - parentOctant = i; - break; - } - } - - if (parentOctant == 255) - { - FatalErrorIn("walkToNeighbour") - << abort(FatalError); - } nodeI = parentNodeI; octant = parentOctant; } - // So now we hit the correct parent node. Switch to the correct octant + // So now we hit the correct parent node. Switch to the correct octant. + // Is done by jumping to the 'other half' so if e.g. in x direction in + // right half we now jump to the left half. octant ^= octantMask; + //Pout<< " to node:" << nodeI << " octant:" << octant + // << " subBb:" <::walkToNeighbour(..)") + << "When searching for " << facePoint + << " ended up in node:" << nodeI + << " octant:" << octant + << " with bb:" << subBb + << abort(FatalError); + } + } + + // See if we need to travel down. Note that we already go into the - // the first level ourselves (instead of having findNode decide) since - // the facePoint is now exactly on the mid of the node so there could - // be truncation problems. + // the first level ourselves (instead of having findNode decide) labelBits index = nodes_[nodeI].subNodes_[octant]; if (isNode(index)) @@ -778,49 +1499,82 @@ bool Foam::indexedOctree::walkToNeighbour octant = getOctant(node); } + + if (debug) + { + const treeBoundBox subBb(subBbox(nodeI, octant)); + + if (nodeI == oldNodeI && octant == oldOctant) + { + FatalErrorIn("indexedOctree::walkToNeighbour(..)") + << "Did not go to neighbour when searching for " << facePoint + << endl + << " starting from face:" << faceString(faceID) + << " node:" << nodeI + << " octant:" << octant + << " bb:" << subBb + << abort(FatalError); + } + + if (!subBb.contains(facePoint)) + { + FatalErrorIn("indexedOctree::walkToNeighbour(..)") + << "When searching for " << facePoint + << " ended up in node:" << nodeI + << " octant:" << octant + << " bb:" << subBb + << abort(FatalError); + } + } + + return true; } -// Return unique number for the face of a bounding box a point is on. -// (number is single bit but not really nessecary) -// Return 0 if point not on any face of bb. template -Foam::direction Foam::indexedOctree::getFace +Foam::word Foam::indexedOctree::faceString ( - const treeBoundBox& bb, - const point& pt + const direction faceID ) { - direction faceID = 0; + word desc; - if (pt.x() <= bb.min().x()) + if (faceID == 0) { - faceID |= treeBoundBox::LEFTBIT; + desc = "noFace"; } - if (pt.x() >= bb.max().x()) + if (faceID & treeBoundBox::LEFTBIT) { - faceID |= treeBoundBox::RIGHTBIT; + if (!desc.empty()) desc += "+"; + desc += "left"; } - - if (pt.y() <= bb.min().y()) + if (faceID & treeBoundBox::RIGHTBIT) { - faceID |= treeBoundBox::BOTTOMBIT; + if (!desc.empty()) desc += "+"; + desc += "right"; } - if (pt.y() >= bb.max().y()) + if (faceID & treeBoundBox::BOTTOMBIT) { - faceID |= treeBoundBox::TOPBIT; + if (!desc.empty()) desc += "+"; + desc += "bottom"; } - - if (pt.z() <= bb.min().z()) + if (faceID & treeBoundBox::TOPBIT) { - faceID |= treeBoundBox::BACKBIT; + if (!desc.empty()) desc += "+"; + desc += "top"; } - if (pt.z() >= bb.max().z()) + if (faceID & treeBoundBox::BACKBIT) { - faceID |= treeBoundBox::FRONTBIT; + if (!desc.empty()) desc += "+"; + desc += "back"; } - return faceID; + if (faceID & treeBoundBox::FRONTBIT) + { + if (!desc.empty()) desc += "+"; + desc += "front"; + } + return desc; } @@ -829,21 +1583,37 @@ Foam::direction Foam::indexedOctree::getFace // hitInfo.point = point on shape // Else return a miss and the bounding box face hit: // hitInfo.point = coordinate of intersection of ray with bounding box -// faceID = index of bounding box face +// hitBits = posbits of point on bounding box template void Foam::indexedOctree::traverseNode ( const bool findAny, + const point& treeStart, + const vector& treeVec, + const point& start, const point& end, - const vector& dir, const label nodeI, const direction octant, pointIndexHit& hitInfo, - direction& faceID + direction& hitBits ) const { + if (debug) + { + const treeBoundBox octantBb(subBbox(nodeI, octant)); + + if (octantBb.posBits(start) != 0) + { + FatalErrorIn("indexedOctree::traverseNode(..)") + << "Node:" << nodeI << " octant:" << octant + << " bb:" << octantBb << endl + << "does not contain point " << start << abort(FatalError); + } + } + + const node& nod = nodes_[nodeI]; labelBits index = nod.subNodes_[octant]; @@ -909,49 +1679,72 @@ void Foam::indexedOctree::traverseNode // Nothing intersected in this node // Traverse to end of node. Do by ray tracing back from end and finding // intersection with bounding box of node. + // start point is now considered to be inside bounding box. + + const treeBoundBox octantBb(subBbox(nodeI, octant)); point pt; + bool intersected = octantBb.intersects + ( + end, //treeStart, + (start-end), //treeVec, - if (isNode(index)) + end, + start, + + pt, + hitBits + ); + + + if (intersected) { - const treeBoundBox& subBb = nodes_[getNode(index)].bb_; - - if (!subBb.intersects(end, start, pt)) - { - faceID = 0; - - FatalErrorIn("indexedOctree::traverseNode(..)") - << "Did not hit side of box " << subBb - << " with ray from " << start << " to " << end - << abort(FatalError); - } - else - { - faceID = getFace(subBb, pt); - } + // Return miss. Set misspoint to face. + hitInfo.setPoint(pt); } else { - const treeBoundBox subBb(nod.bb_.subBbox(octant)); + // Rare case: did not find intersection of ray with octantBb. Can + // happen if end is on face/edge of octantBb. Do like + // lagrangian and re-track with perturbed vector (slightly pushed out + // of bounding box) - if (!subBb.intersects(end, start, pt)) - { - faceID = 0; + point perturbedEnd(pushPoint(octantBb, end, false)); - FatalErrorIn("indexedOctree::traverseNode(..)") - << "Did not hit side of box " << subBb - << " with ray from " << start << " to " << end - << abort(FatalError); - } - else + + //if (debug) { - faceID = getFace(subBb, pt); + // Dump octantBb to obj + writeOBJ(nodeI, octant); + // Dump ray to obj as well + { + OFstream str("ray.obj"); + meshTools::writeOBJ(str, start); + meshTools::writeOBJ(str, end); + str << "l 1 2" << nl; + } + WarningIn("indexedOctree::traverseNode(..)") + << "Did not intersect ray from endpoint:" << end + << " to startpoint:" << start + << " with bounding box:" << octantBb << nl + << "Re-intersecting with perturbed endpoint:" << perturbedEnd + << endl; } + + traverseNode + ( + findAny, + treeStart, + treeVec, + start, + perturbedEnd, + nodeI, + octant, + + hitInfo, + hitBits + ); } - - - // Return miss. Set misspoint to face. - hitInfo.setPoint(pt); } @@ -963,40 +1756,75 @@ Foam::pointIndexHit Foam::indexedOctree::findLine const point& treeStart, const point& treeEnd, const label startNodeI, - const direction startOctant + const direction startOctant, + const bool verbose ) const { - const vector dir(treeEnd - treeStart); + const vector treeVec(treeEnd - treeStart); // Current node as parent+octant label nodeI = startNodeI; direction octant = startOctant; + if (verbose) + { + Pout<< "findLine : treeStart:" << treeStart + << " treeEnd:" << treeEnd << endl + << "node:" << nodeI + << " octant:" << octant + << " bb:" << subBbox(nodeI, octant) << endl; + } + // Current position. Initialize to miss pointIndexHit hitInfo(false, treeStart, -1); - // Side of current bounding box current point is on - direction faceID = 0; - - while (true) + //while (true) + label i = 0; + for (; i < 100000; i++) { + if (verbose) + { + Pout<< "iter:" << i + << " at startPoint:" << hitInfo.rawPoint() << endl + << " node:" << nodeI + << " octant:" << octant + << " bb:" << subBbox(nodeI, octant) << endl; + } + + // Ray-trace to end of current node. Updates point (either on triangle // in case of hit or on node bounding box in case of miss) - point startPoint(hitInfo.rawPoint()); + const treeBoundBox octantBb(subBbox(nodeI, octant)); + + // Make sure point is away from any edges/corners + point startPoint + ( + pushPointIntoFace + ( + octantBb, + treeVec, + hitInfo.rawPoint() + ) + ); + + // Faces of current bounding box current point is on + direction hitFaceID = 0; traverseNode ( findAny, + treeStart, + treeVec, + startPoint, // Note: pass in copy since hitInfo // also used in return value. - treeEnd, - dir, + treeEnd, // pass in overall end so is nicely outside bb nodeI, octant, hitInfo, - faceID + hitFaceID ); // Did we hit a triangle? @@ -1005,29 +1833,48 @@ Foam::pointIndexHit Foam::indexedOctree::findLine break; } - if (faceID == 0) + if (hitFaceID == 0) { - // Was never inside the tree. Return miss. + // endpoint inside the tree. Return miss. break; } - //Pout<< "findLine : treeStart:" << treeStart - // << " treeEnd:" << treeEnd - // << " tracked from " << startPoint << " to edge of nodeBb:" - // << hitInfo.missPoint() - // << " node:" << nodeI << " octant:" << octant - // << " subBb:" - // << nodes_[nodeI].bb_.subBbox(octant) - // << endl; + + if (verbose) + { + Pout<< " iter:" << i + << " hit face:" << faceString(hitFaceID) + << " at:" << hitInfo.rawPoint() << nl + << " node:" << nodeI + << " octant:" << octant + << " bb:" << subBbox(nodeI, octant) << nl + << " walking to neighbour containing:" + << pushPoint + ( + octantBb, + hitFaceID, + hitInfo.rawPoint(), + false // push outside of octantBb + ) + << endl; + } // Nothing hit so we are on face of bounding box (given as node+octant+ - // faceID). Traverse to neighbouring node. + // position bits). Traverse to neighbouring node. Use slightly perturbed + // point. bool ok = walkToNeighbour ( - hitInfo.rawPoint(), // point on face - faceID, // direction to walk in + pushPoint + ( + octantBb, + hitFaceID, + hitInfo.rawPoint(), + false // push outside of octantBb + ), + hitFaceID, // face(s) that hitInfo is on + nodeI, octant ); @@ -1037,7 +1884,53 @@ Foam::pointIndexHit Foam::indexedOctree::findLine // Hit the edge of the tree. Return miss. break; } + + if (verbose) + { + const treeBoundBox octantBb(subBbox(nodeI, octant)); + Pout<< " walked for point:" << hitInfo.rawPoint() << endl + << " to neighbour node:" << nodeI + << " octant:" << octant + << " face:" << faceString(octantBb.faceBits(hitInfo.rawPoint())) + << " of octantBb:" << octantBb << endl + << endl; + } } + + if (i == 100000) + { + // Probably in loop. + if (!verbose) + { + // Redo intersection but now with verbose flag switched on. + return findLine + ( + findAny, + treeStart, + treeEnd, + startNodeI, + startOctant, + true //verbose + ); + } + if (debug) + { + FatalErrorIn("indexedOctree::findLine(..)") + << "Got stuck in loop raytracing from:" << treeStart + << " to:" << treeEnd << endl + << "inside top box:" << subBbox(startNodeI, startOctant) + << abort(FatalError); + } + else + { + WarningIn("indexedOctree::findLine(..)") + << "Got stuck in loop raytracing from:" << treeStart + << " to:" << treeEnd << endl + << "inside top box:" << subBbox(startNodeI, startOctant) + << endl; + } + } + return hitInfo; } @@ -1057,6 +1950,9 @@ Foam::pointIndexHit Foam::indexedOctree::findLine { const treeBoundBox& treeBb = nodes_[0].bb_; + // No effort is made to deal with points which are on edge of tree + // bounding box for now. + direction startBit = treeBb.posBits(start); direction endBit = treeBb.posBits(end); @@ -1066,6 +1962,9 @@ Foam::pointIndexHit Foam::indexedOctree::findLine return pointIndexHit(false, vector::zero, -1); } + + // Trim segment to treeBb + point trackStart(start); point trackEnd(end); @@ -1087,6 +1986,7 @@ Foam::pointIndexHit Foam::indexedOctree::findLine } } + // Find lowest level tree node that start is in. labelBits index = findNode(0, trackStart); @@ -1234,32 +2134,6 @@ void Foam::indexedOctree::writeOBJ str<< "l " << e[0]+pointVertI+1 << ' ' << e[1]+pointVertI+1 << nl; } - - - //// Dump triangles - //if (isContent(index)) - //{ - // const labelList& indices = contents_[getContent(index)]; - // const triSurface& surf = shapes_.surface(); - // const pointField& points = surf.points(); - // - // forAll(indices, i) - // { - // label shapeI = indices[i]; - // - // const labelledTri& f = surf[shapeI]; - // - // meshTools::writeOBJ(str, points[f[0]]); - // vertI++; - // meshTools::writeOBJ(str, points[f[1]]); - // vertI++; - // meshTools::writeOBJ(str, points[f[2]]); - // vertI++; - // - // str<< "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << ' ' - // << vertI-2 << nl; - // } - //} } @@ -1305,6 +2179,14 @@ Foam::indexedOctree::indexedOctree contents_(0), nodeTypes_(0) { + if (debug) + { + Pout<< "indexedOctree::indexedOctree:" << nl + << " shapes:" << shapes.size() << nl + << " bb:" << bb << nl + << endl; + } + if (shapes.size() == 0) { return; @@ -1409,7 +2291,6 @@ Foam::indexedOctree::indexedOctree nodes_.transfer(nodes); nodes.clear(); - if (debug) { label nEntries = 0; @@ -1418,14 +2299,18 @@ Foam::indexedOctree::indexedOctree nEntries += contents_[i].size(); } - Pout<< "indexedOctree::indexedOctree : finished construction:" + Pout<< "indexedOctree::indexedOctree" + << " : finished construction of tree of:" << shapes.typeName << nl + << " bb:" << this->bb() << nl << " shapes:" << shapes.size() << nl << " nLevels:" << nLevels << nl << " treeNodes:" << nodes_.size() << nl << " nEntries:" << nEntries << nl - << " per treeLeaf:" << nEntries/contents.size() << nl - << " per shape (duplicity):" << nEntries/shapes.size() << nl + << " per treeLeaf:" + << scalar(nEntries)/contents.size() << nl + << " per shape (duplicity):" + << scalar(nEntries)/shapes.size() << nl << endl; } } @@ -1447,6 +2332,13 @@ Foam::indexedOctree::indexedOctree // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +Foam::scalar& Foam::indexedOctree::perturbTol() +{ + return perturbTol_; +} + + template Foam::pointIndexHit Foam::indexedOctree::findNearest ( @@ -1570,6 +2462,16 @@ Foam::labelBits Foam::indexedOctree::findNode const node& nod = nodes_[nodeI]; + if (debug) + { + if (!nod.bb_.contains(sample)) + { + FatalErrorIn("findNode(..)") + << "Cannot find " << sample << " in node " << nodeI + << abort(FatalError); + } + } + direction octant = nod.bb_.subOctant(sample); labelBits index = nod.subNodes_[octant]; diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H index b4aad9a305..6c2bfa3c6e 100644 --- a/src/meshTools/indexedOctree/indexedOctree.H +++ b/src/meshTools/indexedOctree/indexedOctree.H @@ -128,6 +128,15 @@ public: private: + // Static data + + //- Relative peturbation tolerance. Determines when point is + // considered to be close to face/edge of bb of node. + // The tolerance is relative to the bounding box of the smallest + // node. + static scalar perturbTol_; + + // Private data //- Underlying shapes for geometric queries. @@ -144,8 +153,9 @@ private: // Private Member Functions - //- Like above but now bb is implicitly provided as parent bb + mid - // + octant + //- Helper: does bb intersect a sphere around sample? Or is any + // corner point of bb closer than nearestDistSqr to sample. + // (bb is implicitly provided as parent bb + octant) static bool overlaps ( const treeBoundBox& parentBb, @@ -215,6 +225,60 @@ private: point& nearestPoint ) const; + //- Return bbox of octant + treeBoundBox subBbox + ( + const label parentNodeI, + const direction octant + ) const; + + //- Helper: take a point on/close to face of bb and push it + // inside or outside of bb. + static point pushPoint + ( + const treeBoundBox&, + const point&, + const bool pushInside + ); + + //- Helper: take a point on face of bb and push it + // inside or outside of bb. + static point pushPoint + ( + const treeBoundBox&, + const direction, + const point&, + const bool pushInside + ); + + //- Helper: take point on face(s) of bb and push it away from + // edges of face. + static point pushPointIntoFace + ( + const treeBoundBox& bb, + const vector& dir, // end-start + const point& pt + ); + + ////- Push point on multiple faces away from any corner/edge. + //static void checkMultipleFaces + //( + // const treeBoundBox& bb, + // const vector& dir, // end-start + // pointIndexHit& faceHitInfo, + // direction& faceID + //); + + //- Walk to parent of node+octant. + bool walkToParent + ( + const label nodeI, + const direction octant, + + label& parentNodeI, + label& parentOctant + ) const; + //- Walk tree to neighbouring node. Return false if edge of tree // hit. bool walkToNeighbour @@ -225,8 +289,8 @@ private: direction& octant ) const; - //- Categorize point on bounding box. - static direction getFace(const treeBoundBox&, const point&); + //- Debug: return verbose the bounding box faces + static word faceString(const direction faceID); //- Traverse a node. If intersects a triangle return first // intersection point. @@ -235,9 +299,11 @@ private: void traverseNode ( const bool findAny, + const point& treeStart, + const vector& treeVec, + const point& start, const point& end, - const vector& dir, const label nodeI, const direction octantI, @@ -252,7 +318,8 @@ private: const point& treeStart, const point& treeEnd, const label startNodeI, - const direction startOctantI + const direction startOctantI, + const bool verbose = false ) const; //- Find any or nearest intersection of line between start and end. @@ -310,6 +377,10 @@ private: public: + //- Get the perturbation tolerance + static scalar& perturbTol(); + + // Constructors //- Construct null @@ -505,6 +576,7 @@ public: const point& sample ); + // Write //- Print tree. Either print all indices (printContent = true) or diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C index f86c487f4e..ffcbd0c06c 100644 --- a/src/meshTools/indexedOctree/treeDataFace.C +++ b/src/meshTools/indexedOctree/treeDataFace.C @@ -528,7 +528,7 @@ bool Foam::treeDataFace::intersects intersection::HALF_RAY ); - if (inter.hit() && inter.distance() <= 1) + if (inter.hit() && magSqr(inter.hitPoint()-start) <= magSqr(dir)) { // Note: no extra test on whether intersection is in front of us // since using half_ray diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C index 06a496ee88..c1a16f7128 100644 --- a/src/meshTools/indexedOctree/treeDataTriSurface.C +++ b/src/meshTools/indexedOctree/treeDataTriSurface.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "treeDataTriSurface.H" @@ -226,7 +224,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType if (!pHit.hit()) { - FatalErrorIn("treeDataTriSurface::getVolumeType") + FatalErrorIn("treeDataTriSurface::getVolumeType(..)") << "treeBb:" << treeBb << " sample:" << sample << " pHit:" << pHit @@ -238,7 +236,8 @@ Foam::label Foam::treeDataTriSurface::getVolumeType surface_, sample, pHit.index(), - pHit.hitPoint() + pHit.hitPoint(), + indexedOctree::perturbTol() ); if (t == triSurfaceTools::UNKNOWN) @@ -255,12 +254,13 @@ Foam::label Foam::treeDataTriSurface::getVolumeType } else { - FatalErrorIn("treeDataTriSurface::getVolumeType") + FatalErrorIn("treeDataTriSurface::getVolumeType(..)") << "problem" << abort(FatalError); return indexedOctree::UNKNOWN; } } + // Check if any point on triangle is inside cubeBb. bool Foam::treeDataTriSurface::overlaps ( @@ -305,6 +305,7 @@ bool Foam::treeDataTriSurface::overlaps // Now we have the difficult case: all points are outside but connecting // edges might go through cube. Use fast intersection of bounding box. + //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb); return triangleFuncs::intersectBb(p0, p1, p2, cubeBb); } @@ -445,9 +446,17 @@ bool Foam::treeDataTriSurface::intersects const vector dir(end - start); - pointHit inter = tri.intersection(start, dir, intersection::HALF_RAY); + // Use relative tolerance (from octree) to determine intersection. - if (inter.hit() && inter.distance() <= 1) + pointHit inter = tri.intersection + ( + start, + dir, + intersection::HALF_RAY, + indexedOctree::perturbTol() + ); + + if (inter.hit() && magSqr(inter.hitPoint()-start) <= magSqr(dir)) { // Note: no extra test on whether intersection is in front of us // since using half_ray. @@ -459,6 +468,16 @@ bool Foam::treeDataTriSurface::intersects { return false; } + + + //- Using exact intersections + //pointHit info = f.tri(points).intersectionExact(start, end); + // + //if (info.hit()) + //{ + // intersectionPoint = info.hitPoint(); + //} + //return info.hit(); } diff --git a/src/meshTools/octree/treeBoundBox.C b/src/meshTools/octree/treeBoundBox.C index ce62c47192..6133bc3106 100644 --- a/src/meshTools/octree/treeBoundBox.C +++ b/src/meshTools/octree/treeBoundBox.C @@ -309,12 +309,14 @@ bool Foam::treeBoundBox::overlaps // This makes sure that posBits tests 'inside' bool Foam::treeBoundBox::intersects ( + const point& overallStart, + const vector& overallVec, const point& start, const point& end, - point& pt + point& pt, + direction& ptOnFaces ) const { - const vector vec(end - start); const direction endBits = posBits(end); pt = start; @@ -325,80 +327,106 @@ bool Foam::treeBoundBox::intersects if (ptBits == 0) { // pt inside bb + ptOnFaces = faceBits(pt); return true; } if ((ptBits & endBits) != 0) { // pt and end in same block outside of bb + ptOnFaces = faceBits(pt); return false; } if (ptBits & LEFTBIT) { // Intersect with plane V=min, n=-1,0,0 - if (Foam::mag(vec.x()) > VSMALL) + if (Foam::mag(overallVec.x()) > VSMALL) { - scalar s = (min().x() - pt.x())/vec.x(); + scalar s = (min().x() - overallStart.x())/overallVec.x(); + pt.x() = min().x(); + pt.y() = overallStart.y() + overallVec.y()*s; + pt.z() = overallStart.z() + overallVec.z()*s; + } + else + { + // Vector not in x-direction. But still intersecting bb planes. + // So must be close - just snap to bb. pt.x() = min().x(); - pt.y() = pt.y() + vec.y()*s; - pt.z() = pt.z() + vec.z()*s; } } - if (ptBits & RIGHTBIT) + else if (ptBits & RIGHTBIT) { // Intersect with plane V=max, n=1,0,0 - if (Foam::mag(vec.x()) > VSMALL) + if (Foam::mag(overallVec.x()) > VSMALL) + { + scalar s = (max().x() - overallStart.x())/overallVec.x(); + pt.x() = max().x(); + pt.y() = overallStart.y() + overallVec.y()*s; + pt.z() = overallStart.z() + overallVec.z()*s; + } + else { - scalar s = (max().x() - pt.x())/vec.x(); pt.x() = max().x(); - pt.y() = pt.y() + vec.y()*s; - pt.z() = pt.z() + vec.z()*s; } } - - if (ptBits & BOTTOMBIT) + else if (ptBits & BOTTOMBIT) { // Intersect with plane V=min, n=0,-1,0 - if (Foam::mag(vec.y()) > VSMALL) + if (Foam::mag(overallVec.y()) > VSMALL) { - scalar s = (min().y() - pt.y())/vec.y(); - pt.x() = pt.x() + vec.x()*s; + scalar s = (min().y() - overallStart.y())/overallVec.y(); + pt.x() = overallStart.x() + overallVec.x()*s; pt.y() = min().y(); - pt.z() = pt.z() + vec.z()*s; + pt.z() = overallStart.z() + overallVec.z()*s; + } + else + { + pt.x() = min().y(); } } - if (ptBits & TOPBIT) + else if (ptBits & TOPBIT) { // Intersect with plane V=max, n=0,1,0 - if (Foam::mag(vec.y()) > VSMALL) + if (Foam::mag(overallVec.y()) > VSMALL) + { + scalar s = (max().y() - overallStart.y())/overallVec.y(); + pt.x() = overallStart.x() + overallVec.x()*s; + pt.y() = max().y(); + pt.z() = overallStart.z() + overallVec.z()*s; + } + else { - scalar s = (max().y() - pt.y())/vec.y(); - pt.x() = pt.x() + vec.x()*s; pt.y() = max().y(); - pt.z() = pt.z() + vec.z()*s; } } - - if (ptBits & BACKBIT) + else if (ptBits & BACKBIT) { // Intersect with plane V=min, n=0,0,-1 - if (Foam::mag(vec.z()) > VSMALL) + if (Foam::mag(overallVec.z()) > VSMALL) + { + scalar s = (min().z() - overallStart.z())/overallVec.z(); + pt.x() = overallStart.x() + overallVec.x()*s; + pt.y() = overallStart.y() + overallVec.y()*s; + pt.z() = min().z(); + } + else { - scalar s = (min().z() - pt.z())/vec.z(); - pt.x() = pt.x() + vec.x()*s; - pt.y() = pt.y() + vec.y()*s; pt.z() = min().z(); } } - if (ptBits & FRONTBIT) + else if (ptBits & FRONTBIT) { // Intersect with plane V=max, n=0,0,1 - if (Foam::mag(vec.z()) > VSMALL) + if (Foam::mag(overallVec.z()) > VSMALL) + { + scalar s = (max().z() - overallStart.z())/overallVec.z(); + pt.x() = overallStart.x() + overallVec.x()*s; + pt.y() = overallStart.y() + overallVec.y()*s; + pt.z() = max().z(); + } + else { - scalar s = (max().z() - pt.z())/vec.z(); - pt.x() = pt.x() + vec.x()*s; - pt.y() = pt.y() + vec.y()*s; pt.z() = max().z(); } } @@ -406,6 +434,18 @@ bool Foam::treeBoundBox::intersects } +bool Foam::treeBoundBox::intersects +( + const point& start, + const point& end, + point& pt +) const +{ + direction ptBits; + return intersects(start, end-start, start, end, pt, ptBits); +} + + // this.bb fully contains bb bool Foam::treeBoundBox::contains(const treeBoundBox& bb) const { @@ -452,6 +492,40 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const } +// Code position of pt on bounding box faces +Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const +{ + direction faceBits = 0; + if (pt.x() == min().x()) + { + faceBits |= LEFTBIT; + } + else if (pt.x() == max().x()) + { + faceBits |= RIGHTBIT; + } + + if (pt.y() == min().y()) + { + faceBits |= BOTTOMBIT; + } + else if (pt.y() == max().y()) + { + faceBits |= TOPBIT; + } + + if (pt.z() == min().z()) + { + faceBits |= BACKBIT; + } + else if (pt.z() == max().z()) + { + faceBits |= FRONTBIT; + } + return faceBits; +} + + // Code position of point relative to box Foam::direction Foam::treeBoundBox::posBits(const point& pt) const { @@ -461,7 +535,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const { posBits |= LEFTBIT; } - if (pt.x() > max().x()) + else if (pt.x() > max().x()) { posBits |= RIGHTBIT; } @@ -470,7 +544,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const { posBits |= BOTTOMBIT; } - if (pt.y() > max().y()) + else if (pt.y() > max().y()) { posBits |= TOPBIT; } @@ -479,7 +553,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const { posBits |= BACKBIT; } - if (pt.z() > max().z()) + else if (pt.z() > max().z()) { posBits |= FRONTBIT; } diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H index 77498b2935..69791e04a8 100644 --- a/src/meshTools/octree/treeBoundBox.H +++ b/src/meshTools/octree/treeBoundBox.H @@ -120,12 +120,12 @@ public: enum faceBit { NOFACE = 0, - LEFTBIT = 0x1 << LEFT, - RIGHTBIT = 0x1 << RIGHT, - BOTTOMBIT = 0x1 << BOTTOM, - TOPBIT = 0x1 << TOP, - BACKBIT = 0x1 << BACK, - FRONTBIT = 0x1 << FRONT, + LEFTBIT = 0x1 << LEFT, //1 + RIGHTBIT = 0x1 << RIGHT, //2 + BOTTOMBIT = 0x1 << BOTTOM, //4 + TOPBIT = 0x1 << TOP, //8 + BACKBIT = 0x1 << BACK, //16 + FRONTBIT = 0x1 << FRONT, //32 }; //- Edges codes. @@ -265,9 +265,25 @@ public: //- Overlaps boundingSphere (centre + sqr(radius))? bool overlaps(const point&, const scalar radiusSqr) const; - //- Intersects segment; set point to intersection position, + //- Intersects segment; set point to intersection position and face, // return true if intersection found. - // (intPt argument used during calculation even if not intersecting) + // (pt argument used during calculation even if not intersecting). + // Calculates intersections from outside supplied vector + // (overallStart, overallVec). This is so when + // e.g. tracking through lots of consecutive boxes + // (typical octree) we're not accumulating truncation errors. Set + // to start, (end-start) if not used. + bool intersects + ( + const point& overallStart, + const vector& overallVec, + const point& start, + const point& end, + point& pt, + direction& ptBits + ) const; + + //- Like above but does not return faces point is on bool intersects ( const point& start, @@ -285,6 +301,9 @@ public: // dir would cause it to go inside. bool contains(const vector& dir, const point&) const; + //- Code position of pt on bounding box faces + direction faceBits(const point& pt) const; + //- Position of point relative to bounding box direction posBits(const point&) const; diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C index 10e1fa8f2f..113bdf58f4 100644 --- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C +++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C @@ -238,7 +238,8 @@ void Foam::orientedSurface::propagateOrientation s, samplePoint, nearestFaceI, - nearestPt + nearestPt, + 10*SMALL ); if (side == triSurfaceTools::UNKNOWN) diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 08138f450b..b89afdab63 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -2203,7 +2203,8 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide const triSurface& surf, const point& sample, const label nearestFaceI, // nearest face - const point& nearestPoint // nearest point on nearest face + const point& nearestPoint, // nearest point on nearest face + const scalar tol ) { const labelledTri& f = surf[nearestFaceI]; @@ -2217,7 +2218,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide points[f[0]], points[f[1]], points[f[2]] - ).classify(nearestPoint, 1E-6, nearType, nearLabel); + ).classify(nearestPoint, tol, nearType, nearLabel); if (nearType == triPointRef::NONE) { diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H index 8ddac22a29..b88d460b77 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H @@ -454,13 +454,14 @@ public: //- Given nearest point (to sample) on surface determines which side // sample is. Uses either face normal, edge normal or point normal - // (non-trivial). Feed in output of e.g. triangle::classify. + // (non-trivial). Uses triangle::classify. static sideType surfaceSide ( const triSurface& surf, const point& sample, const label nearestFaceI, // nearest face - const point& nearestPt // nearest point on nearest face + const point& nearestPt, // nearest point on nearest face + const scalar tol // tolerance for nearness test. ); // Triangulation of faces diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C index be83b6326c..1bef529b03 100644 --- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C +++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "triangleFuncs.H" @@ -31,8 +29,6 @@ Description #include "treeBoundBox.H" #include "SortableList.H" #include "boolList.H" -#include "scalar.H" - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -170,7 +166,6 @@ bool Foam::triangleFuncs::intersectAxesBundle } - // Intersect triangle with bounding box. Return true if // any of the faces of bb intersect triangle. // Note: so returns false if triangle inside bb. @@ -262,6 +257,272 @@ bool Foam::triangleFuncs::intersectBb } +//// Intersect triangle with bounding box. Return true if +//// any of the faces of bb intersect triangle. +//// Note: so returns false if triangle inside bb. +//bool Foam::triangleFuncs::intersectBbExact +//( +// const point& p0, +// const point& p1, +// const point& p2, +// const treeBoundBox& cubeBb +//) +//{ +// const point& min = cubeBb.min(); +// const point& max = cubeBb.max(); +// +// const point& cube0 = min; +// const point cube1(min.x(), min.y(), max.z()); +// const point cube2(max.x(), min.y(), max.z()); +// const point cube3(max.x(), min.y(), min.z()); +// +// const point cube4(min.x(), max.y(), min.z()); +// const point cube5(min.x(), max.y(), max.z()); +// const point& cube6 = max; +// const point cube7(max.x(), max.y(), min.z()); +// +// // Test intersection of triangle with twelve edges of box. +// { +// triPointRef tri(p0, p1, p2); +// if (tri.intersectionExact(cube0, cube1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube1, cube2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube2, cube3).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube3, cube0).hit()) +// { +// return true; +// } +// +// if (tri.intersectionExact(cube4, cube5).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube5, cube6).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube6, cube7).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube7, cube4).hit()) +// { +// return true; +// } +// +// if (tri.intersectionExact(cube0, cube4).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube1, cube5).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube2, cube6).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(cube3, cube7).hit()) +// { +// return true; +// } +// } +// // Test intersection of triangle edges with bounding box +// { +// triPointRef tri(cube0, cube1, cube2); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube2, cube3, cube0); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube4, cube5, cube6); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube6, cube7, cube4); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// +// +// { +// triPointRef tri(cube4, cube5, cube1); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube1, cube0, cube4); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube7, cube6, cube2); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube2, cube3, cube7); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// +// { +// triPointRef tri(cube0, cube4, cube7); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube7, cube3, cube0); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube1, cube5, cube6); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// { +// triPointRef tri(cube6, cube2, cube1); +// if (tri.intersectionExact(p0, p1).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p1, p2).hit()) +// { +// return true; +// } +// if (tri.intersectionExact(p2, p0).hit()) +// { +// return true; +// } +// } +// return false; +//} + + bool Foam::triangleFuncs::intersect ( const point& va0, @@ -470,145 +731,4 @@ bool Foam::triangleFuncs::intersect } -bool Foam::triangleFuncs::classify -( - const point& baseVertex, - const vector& E0, - const vector& E1, - const vector& n, - const point& pInter, - const scalar tol, - label& nearType, - label& nearLabel -) -{ - // Get largest component of normal - scalar magX = Foam::mag(n.x()); - scalar magY = Foam::mag(n.y()); - scalar magZ = Foam::mag(n.z()); - - label i0 = -1; - if ((magX >= magY) && (magX >= magZ)) - { - i0 = 0; - } - else if ((magY >= magX) && (magY >= magZ)) - { - i0 = 1; - } - else - { - i0 = 2; - } - - // Get other components - label i1 = (i0 + 1) % 3; - label i2 = (i1 + 1) % 3; - - scalar u1 = E0[i1]; - scalar v1 = E0[i2]; - - scalar u2 = E1[i1]; - scalar v2 = E1[i2]; - - scalar det = v2*u1 - u2*v1; - - scalar u0 = pInter[i1] - baseVertex[i1]; - scalar v0 = pInter[i2] - baseVertex[i2]; - - scalar alpha = 0; - scalar beta = 0; - - bool hit = false; - - if (Foam::mag(u1) < ROOTVSMALL) - { - beta = u0/u2; - - alpha = (v0 - beta*v2)/v1; - - hit = ((alpha >= 0) && ((alpha + beta) <= 1)); - } - else - { - beta = (v0*u1 - u0*v1)/det; - - alpha = (u0 - beta*u2)/u1; - - hit = ((alpha >= 0) && ((alpha + beta) <= 1)); - } - - // - // Now alpha, beta are the coordinates in the triangle local coordinate - // system E0, E1 - // - - nearType = NONE; - nearLabel = -1; - - if (mag(alpha+beta-1) < tol) - { - // On edge between vert 1-2 (=edge 1) - - if (mag(alpha) < tol) - { - nearType = POINT; - nearLabel = 2; - } - else if (mag(beta) < tol) - { - nearType = POINT; - nearLabel = 1; - } - else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1)) - { - nearType = EDGE; - nearLabel = 1; - } - } - else if (mag(alpha) < tol) - { - // On edge between vert 2-0 (=edge 2) - - if (mag(beta) < tol) - { - nearType = POINT; - nearLabel = 0; - } - else if (mag(beta-1) < tol) - { - nearType = POINT; - nearLabel = 2; - } - else if ((beta >= 0) && (beta <= 1)) - { - nearType = EDGE; - nearLabel = 2; - } - } - else if (mag(beta) < tol) - { - // On edge between vert 0-1 (= edge 0) - - if (mag(alpha) < tol) - { - nearType = POINT; - nearLabel = 0; - } - else if (mag(alpha-1) < tol) - { - nearType = POINT; - nearLabel = 1; - } - else if ((alpha >= 0) && (alpha <= 1)) - { - nearType = EDGE; - nearLabel = 0; - } - } - - return hit; -} - - // ************************************************************************* // diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H index 1180e10195..0aa2a4c1c4 100644 --- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H +++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.H @@ -147,25 +147,6 @@ public: point& pInter0, point& pInter1 ); - - - //- Classify point on triangle plane w.r.t. triangle edges. - // - inside/outside (returned) - // - near point (0, 1, 2) - // - near edge (0, 1, 2). Note: edges are counted from starting vertex - // so e.g. edge 2 is from f[2] to f[0] - static bool classify - ( - const point& baseVertex, - const vector& E0, - const vector& E1, - const vector& n, - const point& pInter, - const scalar tol, - label& nearType, - label& nearLabel - ); - };