From d7e880289647c3bad384aaa0907239002aa7ea79 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 11 Nov 2009 17:19:52 +0000 Subject: [PATCH] removed octreeDataFaceList; created generic PrimitivePatch search one --- src/dynamicMesh/Make/files | 1 - src/dynamicMesh/boundaryMesh/boundaryMesh.C | 138 +++-- .../boundaryMesh/octreeDataFaceList.C | 573 ------------------ .../boundaryMesh/octreeDataFaceList.H | 220 ------- src/meshTools/Make/files | 1 + .../indexedOctree/treeDataPrimitivePatch.C | 567 +++++++++++++++++ .../indexedOctree/treeDataPrimitivePatch.H | 210 +++++++ .../treeDataPrimitivePatchName.C | 33 + 8 files changed, 888 insertions(+), 855 deletions(-) delete mode 100644 src/dynamicMesh/boundaryMesh/octreeDataFaceList.C delete mode 100644 src/dynamicMesh/boundaryMesh/octreeDataFaceList.H create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatch.C create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatch.H create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatchName.C diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files index e59bbe3800..8201b7236d 100644 --- a/src/dynamicMesh/Make/files +++ b/src/dynamicMesh/Make/files @@ -49,7 +49,6 @@ perfectInterface/perfectInterface.C boundaryMesh/boundaryMesh.C boundaryPatch/boundaryPatch.C -boundaryMesh/octreeDataFaceList.C setUpdater/setUpdater.C meshModifiers = meshCut/meshModifiers diff --git a/src/dynamicMesh/boundaryMesh/boundaryMesh.C b/src/dynamicMesh/boundaryMesh/boundaryMesh.C index 494cda9295..b71e0a569e 100644 --- a/src/dynamicMesh/boundaryMesh/boundaryMesh.C +++ b/src/dynamicMesh/boundaryMesh/boundaryMesh.C @@ -29,11 +29,12 @@ License #include "polyMesh.H" #include "repatchPolyTopoChanger.H" #include "faceList.H" -#include "octree.H" -#include "octreeDataFaceList.H" +#include "indexedOctree.H" +#include "treeDataPrimitivePatch.H" #include "triSurface.H" #include "SortableList.H" #include "OFstream.H" +#include "uindirectPrimitivePatch.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -892,6 +893,17 @@ Foam::labelList Foam::boundaryMesh::getNearest << endl; } + uindirectPrimitivePatch leftPatch + ( + UIndirectList(mesh(), leftFaces), + mesh().points() + ); + uindirectPrimitivePatch rightPatch + ( + UIndirectList(mesh(), rightFaces), + mesh().points() + ); + // Overall bb treeBoundBox overallBb(mesh().localPoints()); @@ -911,29 +923,35 @@ Foam::labelList Foam::boundaryMesh::getNearest bbMax.z() += 2*tol; // Create the octrees - octree leftTree + indexedOctree + < + treeDataPrimitivePatch + > leftTree ( - overallBb, - octreeDataFaceList + treeDataPrimitivePatch ( - mesh(), - leftFaces + false, // cacheBb + leftPatch ), - 1, - 20, - 10 + overallBb, + 10, // maxLevel + 10, // leafSize + 3.0 // duplicity ); - octree rightTree + indexedOctree + < + treeDataPrimitivePatch + > rightTree ( - overallBb, - octreeDataFaceList + treeDataPrimitivePatch ( - mesh(), - rightFaces + false, // cacheBb + rightPatch ), - 1, - 20, - 10 + overallBb, + 10, // maxLevel + 10, // leafSize + 3.0 // duplicity ); if (debug) @@ -953,7 +971,7 @@ Foam::labelList Foam::boundaryMesh::getNearest treeBoundBox tightest; - const scalar searchDim = mag(searchSpan); + const scalar searchDimSqr = magSqr(searchSpan); forAll(nearestBFaceI, patchFaceI) { @@ -982,50 +1000,25 @@ Foam::labelList Foam::boundaryMesh::getNearest } // Search right tree - tightest.min() = ctr - searchSpan; - tightest.max() = ctr + searchSpan; - scalar rightDist = searchDim; - label rightI = rightTree.findNearest(ctr, tightest, rightDist); - + pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr); // Search left tree. Note: could start from rightDist bounding box // instead of starting from top. - tightest.min() = ctr - searchSpan; - tightest.max() = ctr + searchSpan; - scalar leftDist = searchDim; - label leftI = leftTree.findNearest(ctr, tightest, leftDist); + pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr); - - if (rightI == -1) + if (rightInfo.hit()) { - // No face found in right tree. - - if (leftI == -1) - { - // No face found in left tree. - nearestBFaceI[patchFaceI] = -1; - } - else - { - // Found in left but not in right. Choose left regardless if - // correct sign. Note: do we want this? - nearestBFaceI[patchFaceI] = leftFaces[leftI]; - } - } - else - { - if (leftI == -1) - { - // Found in right but not in left. Choose right regardless if - // correct sign. Note: do we want this? - nearestBFaceI[patchFaceI] = rightFaces[rightI]; - } - else + if (leftInfo.hit()) { // Found in both trees. Compare normals. + label rightFaceI = rightFaces[rightInfo.index()]; + label leftFaceI = leftFaces[leftInfo.index()]; - scalar rightSign = n & ns[rightFaces[rightI]]; - scalar leftSign = n & ns[leftFaces[leftI]]; + label rightDist = mag(rightInfo.hitPoint()-ctr); + label leftDist = mag(leftInfo.hitPoint()-ctr); + + scalar rightSign = n & ns[rightFaceI]; + scalar leftSign = n & ns[leftFaceI]; if ( @@ -1036,11 +1029,11 @@ Foam::labelList Foam::boundaryMesh::getNearest // Both same sign. Choose nearest. if (rightDist < leftDist) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } else @@ -1059,11 +1052,11 @@ Foam::labelList Foam::boundaryMesh::getNearest // Different sign and nearby. Choosing matching normal if (rightSign > 0) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } else @@ -1071,15 +1064,38 @@ Foam::labelList Foam::boundaryMesh::getNearest // Different sign but faraway. Choosing nearest. if (rightDist < leftDist) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } } } + else + { + // Found in right but not in left. Choose right regardless if + // correct sign. Note: do we want this? + label rightFaceI = rightFaces[rightInfo.index()]; + nearestBFaceI[patchFaceI] = rightFaceI; + } + } + else + { + // No face found in right tree. + + if (leftInfo.hit()) + { + // Found in left but not in right. Choose left regardless if + // correct sign. Note: do we want this? + nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()]; + } + else + { + // No face found in left tree. + nearestBFaceI[patchFaceI] = -1; + } } } diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C deleted file mode 100644 index c44c455c74..0000000000 --- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C +++ /dev/null @@ -1,573 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public 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 "octreeDataFaceList.H" -#include "octree.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -defineTypeNameAndDebug(Foam::octreeDataFaceList, 0); - -Foam::scalar Foam::octreeDataFaceList::tol = 1E-6; - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::octreeDataFaceList::calcBb() -{ - allBb_.setSize(faceLabels_.size()); - allBb_ = treeBoundBox - ( - vector(GREAT, GREAT, GREAT), - vector(-GREAT, -GREAT, -GREAT) - ); - - forAll (faceLabels_, faceLabelI) - { - label faceI = faceLabels_[faceLabelI]; - - // Update bb of face - treeBoundBox& myBb = allBb_[faceLabelI]; - - const face& f = mesh_.localFaces()[faceI]; - - forAll(f, fp) - { - const point& coord = mesh_.localPoints()[f[fp]]; - - myBb.min() = min(myBb.min(), coord); - myBb.max() = max(myBb.max(), coord); - } - } -} - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from all faces in bMesh -Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh) -: - mesh_(mesh), - faceLabels_(mesh_.size()), - allBb_(mesh_.size()) -{ - forAll(faceLabels_, faceI) - { - faceLabels_[faceI] = faceI; - } - - // Generate tight fitting bounding box - calcBb(); -} - - -// Construct from selected faces in bMesh -Foam::octreeDataFaceList::octreeDataFaceList -( - const bMesh& mesh, - const labelList& faceLabels -) -: - mesh_(mesh), - faceLabels_(faceLabels), - allBb_(faceLabels.size()) -{ - // Generate tight fitting bounding box - calcBb(); -} - - - -// Construct as copy -Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes) -: - mesh_(shapes.mesh()), - faceLabels_(shapes.faceLabels()), - allBb_(shapes.allBb()) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::octreeDataFaceList::~octreeDataFaceList() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -Foam::label Foam::octreeDataFaceList::getSampleType -( - const octree& oc, - const point& sample -) const -{ - // Need to determine whether sample is 'inside' or 'outside' - // Done by finding nearest face. This gives back a face which is - // guaranteed to contain nearest point. This point can be - // - in interior of face: compare to face normal - // - on edge of face: compare to edge normal - // - on point of face: compare to point normal - // Unfortunately the octree does not give us back the intersection point - // or where on the face it has hit so we have to recreate all that - // information. - - - // Find nearest face to sample - treeBoundBox tightest(treeBoundBox::greatBox); - - scalar tightestDist = GREAT; - - label index = oc.findNearest(sample, tightest, tightestDist); - - if (index == -1) - { - FatalErrorIn - ( - "octreeDataFaceList::getSampleType" - "(octree&, const point&)" - ) << "Could not find " << sample << " in octree." - << abort(FatalError); - } - - label faceI = faceLabels_[index]; - - // Get actual intersection point on face - - if (debug & 2) - { - Pout<< "getSampleType : sample:" << sample - << " nearest face:" << faceI; - } - - const face& f = mesh_.localFaces()[faceI]; - - const pointField& points = mesh_.localPoints(); - - pointHit curHit = f.nearestPoint(sample, points); - - // - // 1] Check whether sample is above face - // - - if (curHit.hit()) - { - // Simple case. Compare to face normal. - - if (debug & 2) - { - Pout<< " -> face hit:" << curHit.hitPoint() - << " comparing to face normal " << mesh_.faceNormals()[faceI] - << endl; - } - return octree::getVolType - ( - mesh_.faceNormals()[faceI], - sample - curHit.hitPoint() - ); - } - - if (debug & 2) - { - Pout<< " -> face miss:" << curHit.missPoint(); - } - - // - // 2] Check whether intersection is on one of the face vertices or - // face centre - // - - // typical dimension as sqrt of face area. - scalar typDim = sqrt(mag(f.normal(points))) + VSMALL; - - forAll(f, fp) - { - if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol) - { - // Face intersection point equals face vertex fp - - if (debug & 2) - { - Pout<< " -> face point hit :" << points[f[fp]] - << " point normal:" << mesh_.pointNormals()[f[fp]] - << " distance:" - << mag(points[f[fp]] - curHit.missPoint())/typDim - << endl; - } - return octree::getVolType - ( - mesh_.pointNormals()[f[fp]], - sample - curHit.missPoint() - ); - } - } - - // Get face centre - point ctr(f.centre(points)); - - if ((mag(ctr - curHit.missPoint())/typDim) < tol) - { - // Face intersection point equals face centre. Normal at face centre - // is already average of face normals - - if (debug & 2) - { - Pout<< " -> centre hit:" << ctr - << " distance:" - << mag(ctr - curHit.missPoint())/typDim - << endl; - } - - return octree::getVolType - ( - mesh_.faceNormals()[faceI], - sample - curHit.missPoint() - ); - } - - - // - // 3] Get the 'real' edge the face intersection is on - // - - const labelList& myEdges = mesh_.faceEdges()[faceI]; - - forAll(myEdges, myEdgeI) - { - const edge& e = mesh_.edges()[myEdges[myEdgeI]]; - - pointHit edgeHit = - line - ( - points[e.start()], - points[e.end()] - ).nearestDist(sample); - - point edgePoint; - if (edgeHit.hit()) - { - edgePoint = edgeHit.hitPoint(); - } - else - { - edgePoint = edgeHit.missPoint(); - } - - - if ((mag(edgePoint - curHit.missPoint())/typDim) < tol) - { - // Face intersection point lies on edge e - - // Calculate edge normal (wrong: uses face normals instead of - // triangle normals) - const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]]; - - vector edgeNormal(vector::zero); - - forAll(myFaces, myFaceI) - { - edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]]; - } - - if (debug & 2) - { - Pout<< " -> real edge hit point:" << edgePoint - << " comparing to edge normal:" << edgeNormal - << endl; - } - - // Found face intersection point on this edge. Compare to edge - // normal - return octree::getVolType - ( - edgeNormal, - sample - curHit.missPoint() - ); - } - } - - - // - // 4] Get the internal edge (vertex - faceCentre) the face intersection - // is on - // - - forAll(f, fp) - { - pointHit edgeHit = - line - ( - points[f[fp]], - ctr - ).nearestDist(sample); - - point edgePoint; - if (edgeHit.hit()) - { - edgePoint = edgeHit.hitPoint(); - } - else - { - edgePoint = edgeHit.missPoint(); - } - - if ((mag(edgePoint - curHit.missPoint())/typDim) < tol) - { - // Face intersection point lies on edge between two face triangles - - // Calculate edge normal as average of the two triangle normals - label fpPrev = f.rcIndex(fp); - label fpNext = f.fcIndex(fp); - - vector e = points[f[fp]] - ctr; - vector ePrev = points[f[fpPrev]] - ctr; - vector eNext = points[f[fpNext]] - ctr; - - vector nLeft = ePrev ^ e; - nLeft /= mag(nLeft) + VSMALL; - - vector nRight = e ^ eNext; - nRight /= mag(nRight) + VSMALL; - - if (debug & 2) - { - Pout<< " -> internal edge hit point:" << edgePoint - << " comparing to edge normal " - << 0.5*(nLeft + nRight) - << endl; - } - - // Found face intersection point on this edge. Compare to edge - // normal - return octree::getVolType - ( - 0.5*(nLeft + nRight), - sample - curHit.missPoint() - ); - } - } - - if (debug & 2) - { - Pout<< "Did not find sample " << sample - << " anywhere related to nearest face " << faceI << endl - << "Face:"; - - forAll(f, fp) - { - Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]] - << endl; - } - } - - // Can't determine status of sample with respect to nearest face. - // Either - // - tolerances are wrong. (if e.g. face has zero area) - // - or (more likely) surface is not closed. - - return octree::UNKNOWN; -} - - -bool Foam::octreeDataFaceList::overlaps -( - const label index, - const treeBoundBox& sampleBb -) const -{ - return sampleBb.overlaps(allBb_[index]); -} - - -bool Foam::octreeDataFaceList::contains -( - const label, - const point& -) const -{ - notImplemented - ( - "octreeDataFaceList::contains(const label, const point&)" - ); - return false; -} - - -bool Foam::octreeDataFaceList::intersects -( - const label index, - const point& start, - const point& end, - point& intersectionPoint -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - const vector dir(end - start); - - // Disable picking up intersections behind us. - scalar oldTol = intersection::setPlanarTol(0.0); - - pointHit inter = - f.ray - ( - start, - dir, - mesh_.localPoints(), - intersection::HALF_RAY, - intersection::VECTOR - ); - - intersection::setPlanarTol(oldTol); - - if (inter.hit() && inter.distance() <= mag(dir)) - { - intersectionPoint = inter.hitPoint(); - - return true; - } - else - { - return false; - } -} - - -bool Foam::octreeDataFaceList::findTightest -( - const label index, - const point& sample, - treeBoundBox& tightest -) const -{ - // Get nearest and furthest away vertex - point myNear, myFar; - allBb_[index].calcExtremities(sample, myNear, myFar); - - const point dist = myFar - sample; - scalar myFarDist = mag(dist); - - point tightestNear, tightestFar; - tightest.calcExtremities(sample, tightestNear, tightestFar); - - scalar tightestFarDist = mag(tightestFar - sample); - - if (tightestFarDist < myFarDist) - { - // Keep current tightest. - return false; - } - else - { - // Construct bb around sample and myFar - const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z())); - - tightest.min() = sample - dist2; - tightest.max() = sample + dist2; - - return true; - } -} - - -// Determine numerical value of sign of sample compared to shape at index -Foam::scalar Foam::octreeDataFaceList::calcSign -( - const label index, - const point& sample, - vector& -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - point ctr = f.centre(mesh_.localPoints()); - - vector vec = sample - ctr; - - vec /= mag(vec) + VSMALL; - - return mesh_.faceNormals()[faceI] & vec; -} - - -// Calculate nearest point on/in shapei -Foam::scalar Foam::octreeDataFaceList::calcNearest -( - const label index, - const point& sample, - point& nearest -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints()); - - if (nearHit.hit()) - { - nearest = nearHit.hitPoint(); - } - else - { - nearest = nearHit.missPoint(); - } - - if (debug & 1) - { - point ctr = f.centre(mesh_.localPoints()); - - scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest); - - Pout<< "octreeDataFaceList::calcNearest : " - << "sample:" << sample - << " faceI:" << faceI - << " ctr:" << ctr - << " sign:" << sign - << " nearest point:" << nearest - << " distance to face:" << nearHit.distance() - << endl; - } - return nearHit.distance(); -} - - -void Foam::octreeDataFaceList::write -( - Ostream& os, - const label index -) const -{ - os << faceLabels_[index] << " " << allBb_[index]; -} - - -// ************************************************************************* // diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H deleted file mode 100644 index addd71d79b..0000000000 --- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H +++ /dev/null @@ -1,220 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::octreeDataFaceList - -Description - Holds data for octree to work on list of faces on a bMesh - (= PrimitivePatch which holds faces, not references them) - Same as octreeDataFace except for that. - -SourceFiles - octreeDataFaceList.C - -\*---------------------------------------------------------------------------*/ - -#ifndef octreeDataFaceList_H -#define octreeDataFaceList_H - -#include "treeBoundBoxList.H" -#include "faceList.H" -#include "point.H" -#include "className.H" -#include "bMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -template class octree; - -/*---------------------------------------------------------------------------*\ - Class octreeDataFaceList Declaration -\*---------------------------------------------------------------------------*/ - -class octreeDataFaceList -{ - // Static data - - //- tolerance on linear dimensions - static scalar tol; - - - // Static function - - static inline label nexti(label max, label i) - { - return (i + 1) % max; - } - - - // Private data - - //- the mesh - const bMesh& mesh_; - - //- labels (in mesh indexing) of faces - labelList faceLabels_; - - //- bbs for all above faces - treeBoundBoxList allBb_; - - - // Private Member Functions - - //- Set allBb to tight fitting bounding box - void calcBb(); - -public: - - // Declare name of the class and its debug switch - ClassName("octreeDataFaceList"); - - // Constructors - - //- Construct from all faces in bMesh. - octreeDataFaceList(const bMesh& mesh); - - //- Construct from selected faces in bMesh. - octreeDataFaceList(const bMesh& mesh, const labelList& faceLabels); - - //- Construct as copy - octreeDataFaceList(const octreeDataFaceList&); - - - // Destructor - - ~octreeDataFaceList(); - - - // Member Functions - - // Access - - const bMesh& mesh() const - { - return mesh_; - } - - const labelList& faceLabels() const - { - return faceLabels_; - } - - const treeBoundBoxList& allBb() const - { - return allBb_; - } - - label size() const - { - return allBb_.size(); - } - - // Search - - //- Get type of sample - label getSampleType - ( - const octree&, - const point& - ) const; - - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; - - //- Does shape at index contain sample - bool contains - ( - const label index, - const point& sample - ) const; - - //- Segment (from start to end) intersection with shape - // at index. If intersects returns true and sets intersectionPoint - bool intersects - ( - const label index, - const point& start, - const point& end, - point& intersectionPoint - ) const; - - //- Sets newTightest to bounding box (and returns true) if - // nearer to sample than tightest bounding box. Otherwise - // returns false. - bool findTightest - ( - const label index, - const point& sample, - treeBoundBox& tightest - ) const; - - //- Given index get unit normal and calculate (numerical) sign - // of sample. - // Used to determine accuracy of calcNearest or inside/outside. - scalar calcSign - ( - const label index, - const point& sample, - vector& n - ) const; - - //- Calculates nearest (to sample) point in shape. - // Returns point and mag(nearest - sample). Returns GREAT if - // sample does not project onto (triangle decomposition) of face. - scalar calcNearest - ( - const label index, - const point& sample, - point& nearest - ) const; - - - // Edit - - // Write - - //- Write shape at index - void write(Ostream& os, const label index) const; -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -#endif - -// ************************************************************************* // diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 272f3e067c..ada844b5a5 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -50,6 +50,7 @@ indexedOctree/treeDataCell.C indexedOctree/treeDataEdge.C indexedOctree/treeDataFace.C indexedOctree/treeDataPoint.C +indexedOctree/treeDataPrimitivePatchName.C indexedOctree/treeDataTriSurface.C searchableSurface = searchableSurface diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C new file mode 100644 index 0000000000..97b7821b1e --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C @@ -0,0 +1,567 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "treeDataPrimitivePatch.H" +#include "indexedOctree.H" +#include "triangleFuncs.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::scalar +Foam::treeDataPrimitivePatch:: +tolSqr = sqr(1E-6); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::treeBoundBox +Foam::treeDataPrimitivePatch:: +calcBb +( + const pointField& points, + const face& f +) +{ + treeBoundBox bb(points[f[0]], points[f[0]]); + + for (label fp = 1; fp < f.size(); fp++) + { + const point& p = points[f[fp]]; + + bb.min() = min(bb.min(), p); + bb.max() = max(bb.max(), p); + } + return bb; +} + + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +void Foam::treeDataPrimitivePatch:: +update() +{ + if (cacheBb_) + { + bbs_.setSize(patch_.size()); + + forAll(patch_, i) + { + bbs_[i] = calcBb(patch_.points(), patch_[i]); + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::treeDataPrimitivePatch:: +treeDataPrimitivePatch +( + const bool cacheBb, + const PrimitivePatch& patch +) +: + patch_(patch), + cacheBb_(cacheBb) +{ + update(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::pointField +Foam::treeDataPrimitivePatch:: +points() const +{ + pointField cc(patch_.size()); + + forAll(patch_, i) + { + cc[i] = patch_[i].centre(patch_.points()); + } + + return cc; +} + + +//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. +// Only makes sense for closed surfaces. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::label +Foam::treeDataPrimitivePatch:: +getVolumeType +( + const indexedOctree + < + treeDataPrimitivePatch + < + Face, + FaceList, + PointField, + PointType + > + >& oc, + const point& sample +) const +{ + // Need to determine whether sample is 'inside' or 'outside' + // Done by finding nearest face. This gives back a face which is + // guaranteed to contain nearest point. This point can be + // - in interior of face: compare to face normal + // - on edge of face: compare to edge normal + // - on point of face: compare to point normal + // Unfortunately the octree does not give us back the intersection point + // or where on the face it has hit so we have to recreate all that + // information. + + + // Find nearest face to sample + pointIndexHit info = oc.findNearest(sample, sqr(GREAT)); + + if (info.index() == -1) + { + FatalErrorIn + ( + "treeDataPrimitivePatch::getSampleType" + "(indexedOctree&, const point&)" + ) << "Could not find " << sample << " in octree." + << abort(FatalError); + } + + + // Get actual intersection point on face + label faceI = info.index(); + + if (debug & 2) + { + Pout<< "getSampleType : sample:" << sample + << " nearest face:" << faceI; + } + + const pointField& points = patch_.localPoints(); + const face& f = patch_.localFaces()[faceI]; + + // Retest to classify where on face info is. Note: could be improved. We + // already have point. + + pointHit curHit = f.nearestPoint(sample, points); + const vector area = f.normal(points); + const point& curPt = curHit.rawPoint(); + + // + // 1] Check whether sample is above face + // + + if (curHit.hit()) + { + // Nearest point inside face. Compare to face normal. + + if (debug & 2) + { + Pout<< " -> face hit:" << curPt + << " comparing to face normal " << area << endl; + } + return indexedOctree::getSide + ( + area, + sample - curPt + ); + } + + if (debug & 2) + { + Pout<< " -> face miss:" << curPt; + } + + // + // 2] Check whether intersection is on one of the face vertices or + // face centre + // + + const scalar typDimSqr = mag(area) + VSMALL; + + forAll(f, fp) + { + if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point equals face vertex fp + + // Calculate point normal (wrong: uses face normals instead of + // triangle normals) + + return indexedOctree::getSide + ( + patch_.pointNormals()[f[fp]], + sample - curPt + ); + } + } + + const point fc(f.centre(points)); + + if ((magSqr(fc - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point equals face centre. Normal at face centre + // is already average of face normals + + if (debug & 2) + { + Pout<< " -> centre hit:" << fc + << " distance:" << magSqr(fc - curPt)/typDimSqr << endl; + } + + return indexedOctree::getSide + ( + area, + sample - curPt + ); + } + + + + // + // 3] Get the 'real' edge the face intersection is on + // + + const labelList& fEdges = patch_.faceEdges()[faceI]; + + forAll(fEdges, fEdgeI) + { + label edgeI = fEdges[fEdgeI]; + const edge& e = patch_.edges()[edgeI]; + + pointHit edgeHit = e.line(points).nearestDist(sample); + + if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point lies on edge e + + // Calculate edge normal (wrong: uses face normals instead of + // triangle normals) + const labelList& eFaces = patch_.edgeFaces()[edgeI]; + + vector edgeNormal(vector::zero); + + forAll(eFaces, i) + { + edgeNormal += patch_.faceNormal()[eFaces[i]]; + } + + if (debug & 2) + { + Pout<< " -> real edge hit point:" << edgeHit.rawPoint() + << " comparing to edge normal:" << edgeNormal + << endl; + } + + // Found face intersection point on this edge. Compare to edge + // normal + return indexedOctree::getSide + ( + edgeNormal, + sample - curPt + ); + } + } + + + // + // 4] Get the internal edge the face intersection is on + // + + forAll(f, fp) + { + pointHit edgeHit = linePointRef + ( + points[f[fp]], + fc + ).nearestDist(sample); + + if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point lies on edge between two face triangles + + // Calculate edge normal as average of the two triangle normals + vector e = points[f[fp]] - fc; + vector ePrev = points[f[f.rcIndex(fp)]] - fc; + vector eNext = points[f[f.fcIndex(fp)]] - fc; + + vector nLeft = ePrev ^ e; + nLeft /= mag(nLeft) + VSMALL; + + vector nRight = e ^ eNext; + nRight /= mag(nRight) + VSMALL; + + if (debug & 2) + { + Pout<< " -> internal edge hit point:" << edgeHit.rawPoint() + << " comparing to edge normal " + << 0.5*(nLeft + nRight) + << endl; + } + + // Found face intersection point on this edge. Compare to edge + // normal + return indexedOctree::getSide + ( + 0.5*(nLeft + nRight), + sample - curPt + ); + } + } + + if (debug & 2) + { + Pout<< "Did not find sample " << sample + << " anywhere related to nearest face " << faceI << endl + << "Face:"; + + forAll(f, fp) + { + Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]] + << endl; + } + } + + // Can't determine status of sample with respect to nearest face. + // Either + // - tolerances are wrong. (if e.g. face has zero area) + // - or (more likely) surface is not closed. + + return indexedOctree::UNKNOWN; +} + + +// Check if any point on shape is inside cubeBb. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +bool +Foam::treeDataPrimitivePatch:: +overlaps +( + const label index, + const treeBoundBox& cubeBb +) const +{ + // 1. Quick rejection: bb does not intersect face bb at all + if (cacheBb_) + { + if (!cubeBb.overlaps(bbs_[index])) + { + return false; + } + } + else + { + if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index]))) + { + return false; + } + } + + + // 2. Check if one or more face points inside + + const pointField& points = patch_.points(); + const face& f = patch_[index]; + + forAll(f, fp) + { + if (cubeBb.contains(points[f[fp]])) + { + return true; + } + } + + // 3. Difficult case: all points are outside but connecting edges might + // go through cube. Use triangle-bounding box intersection. + const point fc = f.centre(points); + + forAll(f, fp) + { + bool triIntersects = triangleFuncs::intersectBb + ( + points[f[fp]], + points[f[f.fcIndex(fp)]], + fc, + cubeBb + ); + + if (triIntersects) + { + return true; + } + } + return false; +} + + +// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex, +// nearestPoint. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +void +Foam::treeDataPrimitivePatch:: +findNearest +( + const labelList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + const pointField& points = patch_.points(); + + forAll(indices, i) + { + label index = indices[i]; + + const face& f = patch_[index]; + + pointHit nearHit = f.nearestPoint(sample, points); + scalar distSqr = sqr(nearHit.distance()); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + minIndex = index; + nearestPoint = nearHit.rawPoint(); + } + } +} + + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +bool +Foam::treeDataPrimitivePatch:: +intersects +( + const label index, + const point& start, + const point& end, + point& intersectionPoint +) const +{ + // Do quick rejection test + if (cacheBb_) + { + const treeBoundBox& faceBb = bbs_[index]; + + if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0) + { + // start and end in same block outside of faceBb. + return false; + } + } + + const pointField& points = patch_.points(); + const face& f = patch_[index]; + const point fc = f.centre(points); + const vector dir(end - start); + + pointHit inter = patch_[index].intersection + ( + start, + dir, + fc, + points, + intersection::HALF_RAY + ); + + if (inter.hit() && inter.distance() <= 1) + { + // Note: no extra test on whether intersection is in front of us + // since using half_ray + intersectionPoint = inter.hitPoint(); + return true; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H new file mode 100644 index 0000000000..9fd4da2fd5 --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H @@ -0,0 +1,210 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::treeDataPrimitivePatch + +Description + Encapsulation of data needed to search on PrimitivePatches + +SourceFiles + treeDataPrimitivePatch.C + +\*---------------------------------------------------------------------------*/ + +#ifndef treeDataPrimitivePatch_H +#define treeDataPrimitivePatch_H + +#include "PrimitivePatch.H" +//#include "indexedOctree.H" +#include "treeBoundBoxList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +template class indexedOctree; + + +/*---------------------------------------------------------------------------*\ + Class treeDataPrimitivePatchName Declaration +\*---------------------------------------------------------------------------*/ + +TemplateName(treeDataPrimitivePatch); + + +/*---------------------------------------------------------------------------*\ + Class treeDataPrimitivePatch Declaration +\*---------------------------------------------------------------------------*/ + +template +< + class Face, + template class FaceList, + class PointField, + class PointType=point +> +class treeDataPrimitivePatch +: + public treeDataPrimitivePatchName +{ + // Static data + + //- tolerance on linear dimensions + static scalar tolSqr; + + // Private data + + //- Underlying geometry + const PrimitivePatch& patch_; + + //- Whether to precalculate and store face bounding box + const bool cacheBb_; + + //- face bounding boxes (valid only if cacheBb_) + treeBoundBoxList bbs_; + + + // Private Member Functions + + //- Calculate face bounding box + static treeBoundBox calcBb(const pointField&, const face&); + + //- Initialise all member data + void update(); + +public: + + // Constructors + + //- Construct from patch. + treeDataPrimitivePatch + ( + const bool cacheBb, + const PrimitivePatch& + ); + + + // Member Functions + + // Access + + label size() const + { + return patch_.size(); + } + + //- Get representative point cloud for all shapes inside + // (one point per shape) + pointField points() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + label getVolumeType + ( + const indexedOctree + < + treeDataPrimitivePatch + < + Face, + FaceList, + PointField, + PointType + > + >&, + const point& + ) const; + + //- Does (bb of) shape at index overlap bb + bool overlaps + ( + const label index, + const treeBoundBox& sampleBb + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; + + //- Calculates nearest (to line) point in shape. + // Returns point and distance (squared) + void findNearest + ( + const labelList& indices, + const linePointRef& ln, + + treeBoundBox& tightest, + label& minIndex, + point& linePoint, + point& nearestPoint + ) const + { + notImplemented + ( + "treeDataPrimitivePatch::findNearest" + "(const labelList&, const linePointRef&, ..)" + ); + } + + //- Calculate intersection of shape with ray. Sets result + // accordingly + bool intersects + ( + const label index, + const point& start, + const point& end, + point& result + ) const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "treeDataPrimitivePatch.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C b/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C new file mode 100644 index 0000000000..59cf44b8a6 --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C @@ -0,0 +1,33 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "treeDataPrimitivePatch.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::treeDataPrimitivePatchName, 0); + +// ************************************************************************* //