From a65c6f234d670937f15fd32123ec32901e18ae89 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 20 Aug 2008 12:10:43 +0100 Subject: [PATCH 01/13] make SortableList::sort always use stable sort --- .../Lists/SortableList/SortableList.C | 23 +------------------ .../Lists/SortableList/SortableList.H | 5 +--- 2 files changed, 2 insertions(+), 26 deletions(-) diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C index 7cde863a63..a4264578d1 100644 --- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C +++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C @@ -86,28 +86,7 @@ void Foam::SortableList::sort() indices_[i] = i; } - Foam::sort(indices_, less(*this)); - - List tmpValues(this->size()); - - forAll(indices_, i) - { - tmpValues[i] = this->operator[](indices_[i]); - } - - List::transfer(tmpValues); -} - - - -template -void Foam::SortableList::stableSort() -{ - forAll(indices_, i) - { - indices_[i] = i; - } - + //Foam::sort(indices_, less(*this)); Foam::stableSort(indices_, less(*this)); List tmpValues(this->size()); diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H index a9ab2ae829..dd0f0af0d4 100644 --- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H +++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H @@ -109,12 +109,9 @@ public: //- Size the list. If grow can cause undefined indices (until next sort) void setSize(const label); - //- Sort the list (if changed after construction time) + //- (stable) sort the list (if changed after construction time) void sort(); - //- Sort the list (if changed after construction time) - void stableSort(); - // Member Operators From ce16ef250eeddebc70e2f5fd7a311ec818c8573a Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 20 Aug 2008 12:11:13 +0100 Subject: [PATCH 02/13] unnecessary include --- .../polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C | 1 - src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C | 2 -- src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C | 1 - src/meshTools/searchableSurface/searchableSphere.C | 1 - .../booleanOps/intersectedSurface/intersectedSurface.C | 2 -- 5 files changed, 7 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C index af86f37fed..f4e073d4e1 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "coupledPolyPatch.H" -#include "SortableList.H" #include "ListOps.H" #include "transform.H" #include "OFstream.H" diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C index fedd050a78..da7e6368c1 100644 --- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C +++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C @@ -24,7 +24,6 @@ License \*---------------------------------------------------------------------------*/ -#include "SortableList.H" #include "dynamicRefineFvMesh.H" #include "addToRunTimeSelectionTable.H" #include "volFields.H" @@ -32,7 +31,6 @@ License #include "surfaceFields.H" #include "fvCFD.H" #include "syncTools.H" -#include "ListListOps.H" #include "pointFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C index 9c6a50177a..6d97285d73 100644 --- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C +++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C @@ -32,7 +32,6 @@ License #include "octreeDataFace.H" #include "octree.H" #include "OFstream.H" -#include "SortableList.H" #include "IndirectList.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/src/meshTools/searchableSurface/searchableSphere.C b/src/meshTools/searchableSurface/searchableSphere.C index 47abfbfff5..04f3f2f2b5 100644 --- a/src/meshTools/searchableSurface/searchableSphere.C +++ b/src/meshTools/searchableSurface/searchableSphere.C @@ -26,7 +26,6 @@ License #include "searchableSphere.H" #include "addToRunTimeSelectionTable.H" -#include "SortableList.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C index f31b30602b..7e716c1679 100644 --- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C +++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C @@ -27,8 +27,6 @@ License #include "intersectedSurface.H" #include "surfaceIntersection.H" #include "faceList.H" -#include "SortableList.H" -#include "triSurfaceTools.H" #include "faceTriangulation.H" #include "treeBoundBox.H" #include "OFstream.H" From b0e8eaf614ec612c5545507c8ebc6333c2cce24e Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 20 Aug 2008 12:11:30 +0100 Subject: [PATCH 03/13] make SortableList::sort always use stable sort --- src/triSurface/triSurface/triSurface.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C index adc28176c5..336987ade8 100644 --- a/src/triSurface/triSurface/triSurface.C +++ b/src/triSurface/triSurface/triSurface.C @@ -626,7 +626,7 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const { sortedRegion[faceI] = operator[](faceI).region(); } - sortedRegion.stableSort(); + sortedRegion.sort(); faceMap = sortedRegion.indices(); From d88bc072a1490739ee58df236a35ef88f1435079 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 17:58:07 +0100 Subject: [PATCH 04/13] schedule calculation; skipping unused procs --- .../mapPolyMesh/mapDistribute/mapDistribute.C | 45 ++++++------ .../mapPolyMesh/mapDistribute/mapDistribute.H | 20 +++--- .../mapDistribute/mapDistributeTemplates.C | 70 +++++++++++++++---- 3 files changed, 93 insertions(+), 42 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C index 796cb0983a..5421b328d8 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C @@ -31,7 +31,11 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::mapDistribute::calcSchedule() const +Foam::List Foam::mapDistribute::schedule +( + const labelListList& subMap, + const labelListList& constructMap +) { // Communications: send and receive processor List allComms; @@ -40,16 +44,16 @@ void Foam::mapDistribute::calcSchedule() const HashSet > commsSet(Pstream::nProcs()); // Find what communication is required - forAll(subMap_, procI) + forAll(subMap, procI) { if (procI != Pstream::myProcNo()) { - if (subMap_[procI].size() > 0) + if (subMap[procI].size() > 0) { // I need to send to procI commsSet.insert(labelPair(Pstream::myProcNo(), procI)); } - if (constructMap_[procI].size() > 0) + if (constructMap[procI].size() > 0) { // I need to receive from procI commsSet.insert(labelPair(procI, Pstream::myProcNo())); @@ -120,13 +124,7 @@ void Foam::mapDistribute::calcSchedule() const ); // Processors involved in my schedule - schedulePtr_.reset - ( - new List - ( - IndirectList(allComms, mySchedule) - ) - ); + return IndirectList(allComms, mySchedule); //if (debug) @@ -152,6 +150,22 @@ void Foam::mapDistribute::calcSchedule() const } +const Foam::List& Foam::mapDistribute::schedule() const +{ + if (!schedulePtr_.valid()) + { + schedulePtr_.reset + ( + new List + ( + schedule(subMap_, constructMap_) + ) + ); + } + return schedulePtr_(); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // //- Construct from components @@ -257,13 +271,4 @@ Foam::mapDistribute::mapDistribute } -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // - - // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H index bdc71c286d..65888ea79e 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H @@ -36,6 +36,8 @@ Note: Schedule is a list of processor pairs (one send, one receive. One of them will be myself) which forms a scheduled (i.e. non-buffered) exchange. See distribute on how to use it. + Note2: number of items send on one processor have to equal the number + of items received on the other processor. SourceFiles @@ -80,8 +82,6 @@ class mapDistribute // Private Member Functions - void calcSchedule() const; - //- Disallow default bitwise copy construct mapDistribute(const mapDistribute&); @@ -142,15 +142,15 @@ public: return constructMap_; } + //- Calculate a schedule. See above. + static List schedule + ( + const labelListList& subMap, + const labelListList& constructMap + ); + //- Return a schedule. Demand driven. See above. - const List& schedule() const - { - if (!schedulePtr_.valid()) - { - calcSchedule(); - } - return schedulePtr_(); - } + const List& schedule() const; // Other diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C index 86ee3d96cf..6ad85c7971 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C @@ -48,15 +48,28 @@ void Foam::mapDistribute::distribute // Send sub field to neighbour for (label domain = 0; domain < Pstream::nProcs(); domain++) { - if (domain != Pstream::myProcNo()) + const labelList& map = subMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) { + List subField(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } OPstream toNbr(Pstream::blocking, domain); - toNbr << IndirectList(field, subMap[domain])(); + toNbr << subField; } } // Subset myself - List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + const labelList& mySubMap = subMap[Pstream::myProcNo()]; + + List subField(mySubMap.size()); + forAll(mySubMap, i) + { + subField[i] = field[mySubMap[i]]; + } // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -71,7 +84,11 @@ void Foam::mapDistribute::distribute // Receive sub field from neighbour for (label domain = 0; domain < Pstream::nProcs(); domain++) { - if (domain != Pstream::myProcNo()) + if + ( + domain != Pstream::myProcNo() + && constructMap[domain].size() > 0 + ) { IPstream fromNbr(Pstream::blocking, domain); List subField(fromNbr); @@ -93,7 +110,13 @@ void Foam::mapDistribute::distribute List newField(constructSize); // Subset myself - List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + const labelList& mySubMap = subMap[Pstream::myProcNo()]; + + List subField(mySubMap.size()); + forAll(mySubMap, i) + { + subField[i] = field[mySubMap[i]]; + } // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -112,8 +135,16 @@ void Foam::mapDistribute::distribute if (Pstream::myProcNo() == sendProc) { // I am sender. Send to recvProc. + const labelList& map = subMap[recvProc]; + + List subField(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + OPstream toNbr(Pstream::scheduled, recvProc); - toNbr << IndirectList(field, subMap[recvProc])(); + toNbr << subField; } else { @@ -136,7 +167,13 @@ void Foam::mapDistribute::distribute List newField(constructSize); // Subset myself - List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + const labelList& mySubMap = subMap[Pstream::myProcNo()]; + + List subField(mySubMap.size()); + forAll(mySubMap, i) + { + subField[i] = field[mySubMap[i]]; + } // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -149,10 +186,19 @@ void Foam::mapDistribute::distribute // Send sub field to neighbour for (label domain = 0; domain < Pstream::nProcs(); domain++) { - if (domain != Pstream::myProcNo()) + const labelList& map = subMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) { + + List subField(map.size()); + forAll(map, i) + { + subField[i] = field[map[i]]; + } + OPstream toNbr(Pstream::nonBlocking, domain); - toNbr << IndirectList(field, subMap[domain])(); + toNbr << subField; } } @@ -160,13 +206,13 @@ void Foam::mapDistribute::distribute // Receive sub field from neighbour for (label domain = 0; domain < Pstream::nProcs(); domain++) { - if (domain != Pstream::myProcNo()) + const labelList& map = constructMap[domain]; + + if (domain != Pstream::myProcNo() && map.size() > 0) { IPstream fromNbr(Pstream::nonBlocking, domain); List subField(fromNbr); - const labelList& map = constructMap[domain]; - forAll(map, i) { newField[map[i]] = subField[i]; From a688017797d8928541dbab4bba9c32650c845936 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 17:59:28 +0100 Subject: [PATCH 05/13] constructor from data with reuse --- .../PrimitivePatch/PrimitivePatch.C | 36 +++++++++++++++++++ .../PrimitivePatch/PrimitivePatch.H | 8 +++++ 2 files changed, 44 insertions(+) diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C index a47c9c7d8c..835f884399 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C +++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C @@ -68,6 +68,42 @@ PrimitivePatch::PrimitivePatch {} +// Construct from components +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +PrimitivePatch::PrimitivePatch +( + FaceList& faces, + Field& points, + const bool reUse +) +: + FaceList(faces, reUse), + points_(points, reUse), + edgesPtr_(NULL), + nInternalEdges_(-1), + boundaryPointsPtr_(NULL), + faceFacesPtr_(NULL), + edgeFacesPtr_(NULL), + faceEdgesPtr_(NULL), + pointEdgesPtr_(NULL), + pointFacesPtr_(NULL), + localFacesPtr_(NULL), + meshPointsPtr_(NULL), + meshPointMapPtr_(NULL), + edgeLoopsPtr_(NULL), + localPointsPtr_(NULL), + localPointOrderPtr_(NULL), + faceNormalsPtr_(NULL), + pointNormalsPtr_(NULL) +{} + + // Construct as copy template < diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H index c22514d25f..49b3be8260 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H +++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H @@ -235,6 +235,14 @@ public: const Field& points ); + //- Construct from components, reuse storage + PrimitivePatch + ( + FaceList& faces, + Field& points, + const bool reUse + ); + //- Construct as copy PrimitivePatch ( From 084473cc0c6953442051a245d9ef73f682c5dd43 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 18:00:02 +0100 Subject: [PATCH 06/13] delayed faceEdges construction to save memory --- .../meshes/primitiveMesh/primitiveMeshEdges.C | 52 ++++++++++++++++--- 1 file changed, 46 insertions(+), 6 deletions(-) diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C index 8b5fa7bb03..6b3039ff30 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C @@ -467,7 +467,8 @@ const edgeList& primitiveMesh::edges() const { if (!edgesPtr_) { - calcEdges(true); + //calcEdges(true); + calcEdges(false); } return *edgesPtr_; @@ -477,10 +478,8 @@ const labelListList& primitiveMesh::pointEdges() const { if (!pePtr_) { - //// Invert edges - //pePtr_ = new labelListList(nPoints()); - //invertManyToMany(nPoints(), edges(), *pePtr_); - calcEdges(true); + //calcEdges(true); + calcEdges(false); } return *pePtr_; @@ -491,12 +490,53 @@ const labelListList& primitiveMesh::faceEdges() const { if (!fePtr_) { - calcEdges(true); + if (debug) + { + Pout<< "primitiveMesh::faceEdges() : " + << "calculating faceEdges" << endl; + } + + //calcEdges(true); + const faceList& fcs = faces(); + const labelListList& pe = pointEdges(); + const edgeList& es = edges(); + + fePtr_ = new labelListList(fcs.size()); + labelListList& faceEdges = *fePtr_; + + forAll(fcs, faceI) + { + const face& f = fcs[faceI]; + + labelList& fEdges = faceEdges[faceI]; + fEdges.setSize(f.size()); + + forAll(f, fp) + { + label pointI = f[fp]; + label nextPointI = f[f.fcIndex(fp)]; + + // Find edge between pointI, nextPontI + const labelList& pEdges = pe[pointI]; + + forAll(pEdges, i) + { + label edgeI = pEdges[i]; + + if (es[edgeI].otherVertex(pointI) == nextPointI) + { + fEdges[fp] = edgeI; + break; + } + } + } + } } return *fePtr_; } + void primitiveMesh::clearOutEdges() { deleteDemandDrivenData(edgesPtr_); From 05e733ce3a0149953aed1a10fa3d428eb58ae8e0 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 18:01:17 +0100 Subject: [PATCH 07/13] extra checks --- .../parMetisDecomp/parMetisDecomp.C | 47 +++++++++++++++++-- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C index a5778bcdf0..8b3092f08f 100644 --- a/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C +++ b/src/decompositionAgglomeration/parMetisDecomp/parMetisDecomp.C @@ -41,8 +41,7 @@ extern "C" # include "parmetis.h" } - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { @@ -57,6 +56,8 @@ namespace Foam } +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + //- Does prevention of 0 cell domains and calls parmetis. Foam::label Foam::parMetisDecomp::decompose ( @@ -76,6 +77,16 @@ Foam::label Foam::parMetisDecomp::decompose // Number of dimensions int nDims = 3; + + if (cellCentres.size() != xadj.size()-1) + { + FatalErrorIn("parMetisDecomp::decompose(..)") + << "cellCentres:" << cellCentres.size() + << " xadj:" << xadj.size() + << abort(FatalError); + } + + // Get number of cells on all processors List nLocalCells(Pstream::nProcs()); nLocalCells[Pstream::myProcNo()] = xadj.size()-1; @@ -106,12 +117,12 @@ Foam::label Foam::parMetisDecomp::decompose // Make sure every domain has at least one cell // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // (Metis falls over with zero sized domains) - // Trickle cells from processors that have them down to those that + // Trickle cells from processors that have them up to those that // don't. - // Number of cells to send down (is same as number of cells next processor - // has to receive) + // Number of cells to send to the next processor + // (is same as number of cells next processor has to receive) List nSendCells(Pstream::nProcs(), 0); for (label procI = nLocalCells.size()-1; procI >=1; procI--) @@ -135,6 +146,15 @@ Foam::label Foam::parMetisDecomp::decompose Field prevCellWeights(fromPrevProc); Field prevFaceWeights(fromPrevProc); + if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) + { + FatalErrorIn("parMetisDecomp::decompose(..)") + << "Expected from processor " << Pstream::myProcNo()-1 + << " connectivity for " << nSendCells[Pstream::myProcNo()-1] + << " nCells but only received " << prevXadj.size() + << abort(FatalError); + } + // Insert adjncy prepend(prevAdjncy, adjncy); // Adapt offsets and prepend xadj @@ -222,6 +242,14 @@ Foam::label Foam::parMetisDecomp::decompose } + if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1)) + { + FatalErrorIn("parMetisDecomp::decompose(..)") + << "Have connectivity for " << xadj.size()-1 + << " cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()] + << abort(FatalError); + } + // Weight info int wgtFlag = 0; int* vwgtPtr = NULL; @@ -292,6 +320,15 @@ Foam::label Foam::parMetisDecomp::decompose List nextFinalDecomp(fromNextProc); + if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) + { + FatalErrorIn("parMetisDecomp::decompose(..)") + << "Expected from processor " << Pstream::myProcNo()+1 + << " decomposition for " << nSendCells[Pstream::myProcNo()] + << " nCells but only received " << nextFinalDecomp.size() + << abort(FatalError); + } + append(nextFinalDecomp, finalDecomp); } From bd27a662989b6bc5dedc56fcc2b764f73e2c47bc Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 18:01:31 +0100 Subject: [PATCH 08/13] binary output --- src/meshTools/octree/PointIndexHit.H | 44 ++++++++++++++++++++++++++-- src/meshTools/octree/pointIndexHit.H | 9 +++++- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/meshTools/octree/PointIndexHit.H b/src/meshTools/octree/PointIndexHit.H index 9b72d99605..f29965748c 100644 --- a/src/meshTools/octree/PointIndexHit.H +++ b/src/meshTools/octree/PointIndexHit.H @@ -94,6 +94,12 @@ public: index_(-1) {} + //- Construct from Istream + PointIndexHit(Istream& is) + { + is >> *this; + } + // Member Functions @@ -193,13 +199,45 @@ public: friend Ostream& operator<< (Ostream& os, const PointIndexHit& pHit) { - return os << pHit.hit_ << token::SPACE << pHit.hitPoint_ - << token::SPACE << pHit.index_; + if (os.format() == IOstream::ASCII) + { + os << pHit.hit_ << token::SPACE << pHit.hitPoint_ + << token::SPACE << pHit.index_; + } + else + { + os.write + ( + reinterpret_cast(&pHit), + sizeof(PointIndexHit) + ); + } + + // Check state of Ostream + os.check("Ostream& operator<<(Ostream&, const PointIndexHit&)"); + + return os; } friend Istream& operator>>(Istream& is, PointIndexHit& pHit) { - return is >> pHit.hit_ >> pHit.hitPoint_ >> pHit.index_; + if (is.format() == IOstream::ASCII) + { + return is >> pHit.hit_ >> pHit.hitPoint_ >> pHit.index_; + } + else + { + is.read + ( + reinterpret_cast(&pHit), + sizeof(PointIndexHit) + ); + } + + // Check state of Istream + is.check("Istream& operator>>(Istream&, PointIndexHit&)"); + + return is; } }; diff --git a/src/meshTools/octree/pointIndexHit.H b/src/meshTools/octree/pointIndexHit.H index e41fd28499..b43f12d551 100644 --- a/src/meshTools/octree/pointIndexHit.H +++ b/src/meshTools/octree/pointIndexHit.H @@ -39,7 +39,14 @@ Description namespace Foam { - typedef PointIndexHit pointIndexHit; + +typedef PointIndexHit pointIndexHit; + + +//- Specify data associated with pointIndexHit type is contiguous +template<> +inline bool contiguous() {return true;} + } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From 6e02125e05709d0444edbd8dd318554dacd4b308 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 18:09:32 +0100 Subject: [PATCH 09/13] stats without mesh addressing; reuse constructor --- src/triSurface/triSurface/triSurface.C | 49 +++++++++++++++++++------- src/triSurface/triSurface/triSurface.H | 11 +++++- 2 files changed, 47 insertions(+), 13 deletions(-) diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C index 336987ade8..89a48a30d6 100644 --- a/src/triSurface/triSurface/triSurface.C +++ b/src/triSurface/triSurface/triSurface.C @@ -31,6 +31,7 @@ License #include "Time.H" #include "boundBox.H" #include "SortableList.H" +#include "PackedList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -743,6 +744,26 @@ triSurface::triSurface {} +triSurface::triSurface +( + List& triangles, + const geometricSurfacePatchList& patches, + pointField& points, + const bool reUse +) +: + PrimitivePatch + ( + triangles, + points, + reUse + ), + patches_(patches), + sortedEdgeFacesPtr_(NULL), + edgeOwnerPtr_(NULL) +{} + + triSurface::triSurface ( const List& triangles, @@ -1148,9 +1169,7 @@ triSurface triSurface::subsetMesh } // Construct subsurface - triSurface subSurface(newTriangles, patches(), newPoints); - - return subSurface; + return triSurface(newTriangles, patches(), newPoints, true); } @@ -1187,30 +1206,36 @@ void triSurface::write(const Time& d) const void triSurface::writeStats(Ostream& os) const { - // Calculate bounding box without any additional addressing - // Copy of treeBoundBox code. Cannot use meshTools from triSurface... + // Unfortunately nPoints constructs meshPoints() so do compact version + // ourselves. + PackedList<1> pointIsUsed(points().size()); + pointIsUsed = 0U; + + label nPoints = 0; boundBox bb ( point(VGREAT, VGREAT, VGREAT), point(-VGREAT, -VGREAT, -VGREAT) ); + forAll(*this, triI) { const labelledTri& f = operator[](triI); forAll(f, fp) { - const point& pt = points()[f[fp]]; - bb.min() = ::Foam::min(bb.min(), pt); - bb.max() = ::Foam::max(bb.max(), pt); + label pointI = f[fp]; + if (pointIsUsed.set(pointI, 1)) + { + bb.min() = ::Foam::min(bb.min(), points()[pointI]); + bb.max() = ::Foam::max(bb.max(), points()[pointI]); + nPoints++; + } } } - // Unfortunately nPoints constructs meshPoints() ... - os << "Triangles : " << size() << endl - //<< "Edges : " << nEdges() << endl - << "Vertices : " << nPoints() << endl + << "Vertices : " << nPoints << endl << "Bounding Box : " << bb << endl; } diff --git a/src/triSurface/triSurface/triSurface.H b/src/triSurface/triSurface/triSurface.H index 56b1abc849..624b793117 100644 --- a/src/triSurface/triSurface/triSurface.H +++ b/src/triSurface/triSurface/triSurface.H @@ -36,9 +36,9 @@ SourceFiles #ifndef triSurface_H #define triSurface_H +#include "PrimitivePatch.H" #include "pointField.H" #include "labelledTri.H" -#include "PrimitivePatch.H" #include "boolList.H" #include "geometricSurfacePatchList.H" #include "surfacePatchList.H" @@ -215,6 +215,15 @@ public: const pointField& ); + //- Construct from triangles, patches, points. Reuse storage. + triSurface + ( + List&, + const geometricSurfacePatchList&, + pointField&, + const bool reUse + ); + //- Construct from triangles, points. Set patchnames to default. triSurface(const List&, const pointField&); From 4dc1c78a281f105a0f44bdd01d38cbe91a6216d4 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 27 Aug 2008 18:09:53 +0100 Subject: [PATCH 10/13] handling tiny triangles --- src/meshTools/triSurface/triangleFuncs/triangleFuncs.C | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C index 56f8799aec..e8c9a29552 100644 --- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C +++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C @@ -101,6 +101,7 @@ bool Foam::triangleFuncs::intersectAxesBundle // Since direction is coordinate axis there is no need to do projection, // we can directly check u,v components for inclusion in triangle. + scalar localScale = max(max(magSqr(V10), magSqr(V20)), 1.0); // Get other components label i1 = (i0 + 1) % 3; @@ -114,8 +115,11 @@ bool Foam::triangleFuncs::intersectAxesBundle scalar det = v2*u1 - u2*v1; - // Fix for V0:(-31.71428 0 -15.10714) V10:(-1.285715 8.99165e-16 -1.142858) V20:(0 0 -1.678573) i0:0 - if (Foam::mag(det)/max(max(mag(V10),mag(V20)),1) < SMALL) + // Fix for V0:(-31.71428 0 -15.10714) + // V10:(-1.285715 8.99165e-16 -1.142858) + // V20:(0 0 -1.678573) + // i0:0 + if (Foam::mag(det)/localScale < SMALL) { // Triangle parallel to dir return false; @@ -132,7 +136,7 @@ bool Foam::triangleFuncs::intersectAxesBundle scalar beta = 0; bool inter = false; - if (Foam::mag(u1) < SMALL) + if (Foam::mag(u1)/localScale < SMALL) { beta = u0/u2; if ((beta >= 0) && (beta <= 1)) From 212a70a17b17f9657b3e1bdfbb728d4ce3554869 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 Aug 2008 15:50:11 +0100 Subject: [PATCH 11/13] boundbox intersect renaming --- .../preProcessing/mapFields/mapFields.C | 4 +- src/OpenFOAM/meshes/boundBox/boundBox.H | 2 +- .../boundaryMesh/octreeDataFaceList.C | 2 +- src/meshTools/indexedOctree/indexedOctree.C | 20 +++++----- src/meshTools/indexedOctree/indexedOctree.H | 6 +-- src/meshTools/indexedOctree/treeDataCell.C | 4 +- src/meshTools/indexedOctree/treeDataEdge.C | 4 +- src/meshTools/indexedOctree/treeDataFace.C | 4 +- .../indexedOctree/treeDataTriSurface.C | 4 +- src/meshTools/octree/octreeDataCell.C | 2 +- src/meshTools/octree/octreeDataEdges.C | 2 +- src/meshTools/octree/octreeDataFace.C | 4 +- src/meshTools/octree/treeBoundBox.C | 39 +++++++++++++++++++ src/meshTools/octree/treeBoundBox.H | 7 +++- src/meshTools/octree/treeBoundBoxI.H | 6 +-- src/meshTools/octree/treeNode.C | 16 ++++---- .../octreeData/octreeDataTriSurface.C | 4 +- .../octreeData/octreeDataTriSurfaceTreeLeaf.C | 2 +- 18 files changed, 86 insertions(+), 46 deletions(-) diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C index 0fb8e40fff..5d4a6471e4 100644 --- a/applications/utilities/preProcessing/mapFields/mapFields.C +++ b/applications/utilities/preProcessing/mapFields/mapFields.C @@ -504,7 +504,7 @@ int main(int argc, char *argv[]) !bbsTargetSet[procITarget] || ( bbsTargetSet[procITarget] - && bbsTarget[procITarget].intersects(bbSource) + && bbsTarget[procITarget].overlaps(bbSource) ) ) { @@ -533,7 +533,7 @@ int main(int argc, char *argv[]) bbsTarget[procITarget] = meshTarget.bounds(); bbsTargetSet[procITarget] = true; - if (bbsTarget[procITarget].intersects(bbSource)) + if (bbsTarget[procITarget].overlaps(bbSource)) { if (consistent) { diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H index a11aa68266..40287e53eb 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.H +++ b/src/OpenFOAM/meshes/boundBox/boundBox.H @@ -118,7 +118,7 @@ public: // Query //- Intersects other boundingbox? - bool intersects(const boundBox& bb) const + bool overlaps(const boundBox& bb) const { if ( diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C index 2249bc006e..2ce870af9d 100644 --- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C +++ b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C @@ -399,7 +399,7 @@ bool Foam::octreeDataFaceList::overlaps const treeBoundBox& sampleBb ) const { - return sampleBb.intersects(allBb_[index]); + return sampleBb.overlaps(allBb_[index]); } diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C index f0f8bd94d0..1a55a88547 100644 --- a/src/meshTools/indexedOctree/indexedOctree.C +++ b/src/meshTools/indexedOctree/indexedOctree.C @@ -40,7 +40,7 @@ namespace Foam // Does bb intersect a sphere around sample? Or is any corner point of bb // closer than nearestDistSqr to sample. template -bool indexedOctree::intersects +bool indexedOctree::overlaps ( const point& p0, const point& p1, @@ -84,7 +84,7 @@ bool indexedOctree::intersects // Does bb intersect a sphere around sample? Or is any corner point of bb // closer than nearestDistSqr to sample. template -bool indexedOctree::intersects +bool indexedOctree::overlaps ( const treeBoundBox& parentBb, const direction octant, @@ -94,7 +94,7 @@ bool indexedOctree::intersects { //- Speeded up version of // treeBoundBox subBb(parentBb.subBbox(mid, octant)) - // intersects + // overlaps // ( // subBb.min(), // subBb.max(), @@ -136,7 +136,7 @@ bool indexedOctree::intersects const point mid(0.5*(min+max)); - return intersects(mid, other, nearestDistSqr, sample); + return overlaps(mid, other, nearestDistSqr, sample); } @@ -567,7 +567,7 @@ void indexedOctree::findNearest const treeBoundBox& subBb = nodes_[subNodeI].bb_; - if (intersects(subBb.min(), subBb.max(), nearestDistSqr, sample)) + if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample)) { findNearest ( @@ -584,7 +584,7 @@ void indexedOctree::findNearest { if ( - intersects + overlaps ( nod.bb_, octant, @@ -639,7 +639,7 @@ void indexedOctree::findNearest { const treeBoundBox& subBb = nodes_[getNode(index)].bb_; - if (subBb.intersects(tightest)) + if (subBb.overlaps(tightest)) { findNearest ( @@ -657,7 +657,7 @@ void indexedOctree::findNearest { const treeBoundBox subBb(nodeBb.subBbox(octant)); - if (subBb.intersects(tightest)) + if (subBb.overlaps(tightest)) { shapes_.findNearest ( @@ -1121,7 +1121,7 @@ void indexedOctree::findBox { const treeBoundBox& subBb = nodes_[getNode(index)].bb_; - if (subBb.intersects(searchBox)) + if (subBb.overlaps(searchBox)) { findBox(getNode(index), searchBox, elements); } @@ -1130,7 +1130,7 @@ void indexedOctree::findBox { const treeBoundBox subBb(nodeBb.subBbox(octant)); - if (subBb.intersects(searchBox)) + if (subBb.overlaps(searchBox)) { const labelList& indices = contents_[getContent(index)]; diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H index 249eb1b5f1..97b6ab607c 100644 --- a/src/meshTools/indexedOctree/indexedOctree.H +++ b/src/meshTools/indexedOctree/indexedOctree.H @@ -36,9 +36,9 @@ SourceFiles #ifndef indexedOctree_H #define indexedOctree_H +#include "treeBoundBox.H" #include "pointIndexHit.H" #include "FixedList.H" -#include "treeBoundBox.H" #include "Ostream.H" #include "labelHashSet.H" #include "labelBits.H" @@ -146,7 +146,7 @@ private: //- Like above but now bb is implicitly provided as parent bb + mid // + octant - static bool intersects + static bool overlaps ( const treeBoundBox& parentBb, const direction octant, @@ -497,7 +497,7 @@ public: //- Helper: does bb intersect a sphere around sample? Or is any // corner point of bb closer than nearestDistSqr to sample. - static bool intersects + static bool overlaps ( const point& bbMin, const point& bbMax, diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/meshTools/indexedOctree/treeDataCell.C index 2ade66780c..cb41899306 100644 --- a/src/meshTools/indexedOctree/treeDataCell.C +++ b/src/meshTools/indexedOctree/treeDataCell.C @@ -137,11 +137,11 @@ bool Foam::treeDataCell::overlaps { if (cacheBb_) { - return cubeBb.intersects(bbs_[index]); + return cubeBb.overlaps(bbs_[index]); } else { - return cubeBb.intersects(calcCellBb(cellLabels_[index])); + return cubeBb.overlaps(calcCellBb(cellLabels_[index])); } } diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C index d918fab944..6b8d54399e 100644 --- a/src/meshTools/indexedOctree/treeDataEdge.C +++ b/src/meshTools/indexedOctree/treeDataEdge.C @@ -110,11 +110,11 @@ bool Foam::treeDataEdge::overlaps { if (cacheBb_) { - return cubeBb.intersects(bbs_[index]); + return cubeBb.overlaps(bbs_[index]); } else { - return cubeBb.intersects(calcBb(edgeLabels_[index])); + return cubeBb.overlaps(calcBb(edgeLabels_[index])); } } diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C index e8d148ec23..5764841bbd 100644 --- a/src/meshTools/indexedOctree/treeDataFace.C +++ b/src/meshTools/indexedOctree/treeDataFace.C @@ -412,14 +412,14 @@ bool Foam::treeDataFace::overlaps // 1. Quick rejection: bb does not intersect face bb at all if (cacheBb_) { - if (!cubeBb.intersects(bbs_[index])) + if (!cubeBb.overlaps(bbs_[index])) { return false; } } else { - if (!cubeBb.intersects(calcBb(faceLabels_[index]))) + if (!cubeBb.overlaps(calcBb(faceLabels_[index]))) { return false; } diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C index 23ba7f5116..c226de9499 100644 --- a/src/meshTools/indexedOctree/treeDataTriSurface.C +++ b/src/meshTools/indexedOctree/treeDataTriSurface.C @@ -284,13 +284,13 @@ bool Foam::treeDataTriSurface::overlaps triBb.max() = max(triBb.max(), p2); //- For testing: robust one - //return cubeBb.intersects(triBb); + //return cubeBb.overlaps(triBb); //- Exact test of triangle intersecting bb // Quick rejection. If whole bounding box of tri is outside cubeBb then // there will be no intersection. - if (!cubeBb.intersects(triBb)) + if (!cubeBb.overlaps(triBb)) { return false; } diff --git a/src/meshTools/octree/octreeDataCell.C b/src/meshTools/octree/octreeDataCell.C index 81a3600be9..17bba8b2c9 100644 --- a/src/meshTools/octree/octreeDataCell.C +++ b/src/meshTools/octree/octreeDataCell.C @@ -113,7 +113,7 @@ bool Foam::octreeDataCell::overlaps const treeBoundBox& cubeBb ) const { - return cubeBb.intersects(bbs_[index]); + return cubeBb.overlaps(bbs_[index]); } diff --git a/src/meshTools/octree/octreeDataEdges.C b/src/meshTools/octree/octreeDataEdges.C index 3974a0b91d..9828599e61 100644 --- a/src/meshTools/octree/octreeDataEdges.C +++ b/src/meshTools/octree/octreeDataEdges.C @@ -107,7 +107,7 @@ bool Foam::octreeDataEdges::overlaps const treeBoundBox& sampleBb ) const { - return sampleBb.intersects(allBb_[index]); + return sampleBb.overlaps(allBb_[index]); } diff --git a/src/meshTools/octree/octreeDataFace.C b/src/meshTools/octree/octreeDataFace.C index 46f03993cc..6a8c3f875a 100644 --- a/src/meshTools/octree/octreeDataFace.C +++ b/src/meshTools/octree/octreeDataFace.C @@ -507,12 +507,12 @@ bool Foam::octreeDataFace::overlaps const treeBoundBox& sampleBb ) const { - //return sampleBb.intersects(allBb_[index]); + //return sampleBb.overlaps(allBb_[index]); //- Exact test of face intersecting bb // 1. Quick rejection: bb does not intersect face bb at all - if (!sampleBb.intersects(allBb_[index])) + if (!sampleBb.overlaps(allBb_[index])) { return false; } diff --git a/src/meshTools/octree/treeBoundBox.C b/src/meshTools/octree/treeBoundBox.C index 042209ec7c..26530b5a3a 100644 --- a/src/meshTools/octree/treeBoundBox.C +++ b/src/meshTools/octree/treeBoundBox.C @@ -297,6 +297,45 @@ Foam::treeBoundBox Foam::treeBoundBox::subBbox } +bool Foam::treeBoundBox::overlaps +( + const point& centre, + const scalar radiusSqr +) const +{ + // Find out where centre is in relation to bb. + // Find nearest point on bb. + scalar distSqr = 0; + + for (direction dir = 0; dir < vector::nComponents; dir++) + { + scalar d0 = min()[dir] - centre[dir]; + scalar d1 = max()[dir] - centre[dir]; + + if ((d0 > 0) != (d1 > 0)) + { + // centre inside both extrema. This component does not add any + // distance. + } + else if (Foam::mag(d0) < Foam::mag(d1)) + { + distSqr += d0*d0; + } + else + { + distSqr += d1*d1; + } + + if (distSqr > radiusSqr) + { + return false; + } + } + + return true; +} + + // line intersection. Returns true if line (start to end) inside // bb or intersects bb. Sets pt to intersection. // diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H index 01d7a9a6bd..02dbac38af 100644 --- a/src/meshTools/octree/treeBoundBox.H +++ b/src/meshTools/octree/treeBoundBox.H @@ -263,8 +263,11 @@ public: FixedList& octantOrder ) const; - //- Intersects other boundingbox? - inline bool intersects(const treeBoundBox&) const; + //- Overlaps other boundingbox? + inline bool overlaps(const treeBoundBox&) const; + + //- Overlaps boundingSphere (centre + sqr(radius))? + bool overlaps(const point&, const scalar radiusSqr) const; //- Intersects segment; set point to intersection position, // return true if intersection found. diff --git a/src/meshTools/octree/treeBoundBoxI.H b/src/meshTools/octree/treeBoundBoxI.H index 431eb8a04d..c444b2808a 100644 --- a/src/meshTools/octree/treeBoundBoxI.H +++ b/src/meshTools/octree/treeBoundBoxI.H @@ -24,9 +24,7 @@ License \*---------------------------------------------------------------------------*/ -#include "error.H" #include "treeBoundBox.H" -#include "point.H" #include "Random.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -414,9 +412,9 @@ inline void treeBoundBox::searchOrder // true if bb's intersect or overlap. // Note: <= to make sure we catch all. -inline bool treeBoundBox::intersects(const treeBoundBox& bb) const +inline bool treeBoundBox::overlaps(const treeBoundBox& bb) const { - return boundBox::intersects(bb); + return boundBox::overlaps(bb); } diff --git a/src/meshTools/octree/treeNode.C b/src/meshTools/octree/treeNode.C index 678bcfb694..771ef2cdd9 100644 --- a/src/meshTools/octree/treeNode.C +++ b/src/meshTools/octree/treeNode.C @@ -798,7 +798,7 @@ bool treeNode::findTightest // Node: recurse into subnodes const treeNode* subNodePtr = getNodePtr(octant); - if (subNodePtr->bb().intersects(tightest)) + if (subNodePtr->bb().overlaps(tightest)) { // there might be a better fit inside this subNode changed |= @@ -815,7 +815,7 @@ bool treeNode::findTightest // Leaf: let leaf::find handle this const treeLeaf* subLeafPtr = getLeafPtr(octant); - if (subLeafPtr->bb().intersects(tightest)) + if (subLeafPtr->bb().overlaps(tightest)) { // there might be a better fit inside this subLeaf changed |= @@ -884,7 +884,7 @@ bool treeNode::findNearest // Node const treeNode* subNodePtr = getNodePtr(octant); - if (subNodePtr->bb().intersects(tightest)) + if (subNodePtr->bb().overlaps(tightest)) { // there might be a better fit inside this subNode changed |= @@ -903,7 +903,7 @@ bool treeNode::findNearest // Leaf: let leaf::find handle this const treeLeaf* subLeafPtr = getLeafPtr(octant); - if (subLeafPtr->bb().intersects(tightest)) + if (subLeafPtr->bb().overlaps(tightest)) { // there might be a better fit inside this subNode changed |= @@ -975,7 +975,7 @@ bool treeNode::findNearest // Node const treeNode* subNodePtr = getNodePtr(octant); - if (subNodePtr->bb().intersects(tightest)) + if (subNodePtr->bb().overlaps(tightest)) { // there might be a better fit inside this subNode changed |= @@ -995,7 +995,7 @@ bool treeNode::findNearest // Leaf: let leaf::find handle this const treeLeaf* subLeafPtr = getLeafPtr(octant); - if (subLeafPtr->bb().intersects(tightest)) + if (subLeafPtr->bb().overlaps(tightest)) { // there might be a better fit inside this subNode changed |= @@ -1060,7 +1060,7 @@ bool treeNode::findBox // Node const treeNode* subNodePtr = getNodePtr(octant); - if (subNodePtr->bb().intersects(box)) + if (subNodePtr->bb().overlaps(box)) { // Visit sub node. changed |= subNodePtr->findBox(shapes, box, elements); @@ -1071,7 +1071,7 @@ bool treeNode::findBox // Leaf: let leaf::find handle this const treeLeaf* subLeafPtr = getLeafPtr(octant); - if (subLeafPtr->bb().intersects(box)) + if (subLeafPtr->bb().overlaps(box)) { // Visit sub leaf. changed |= subLeafPtr->findBox(shapes, box, elements); diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C index 6a21879f35..f3dd4f3e03 100644 --- a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C +++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C @@ -385,12 +385,12 @@ bool Foam::octreeDataTriSurface::overlaps const treeBoundBox& cubeBb ) const { - //return cubeBb.intersects(allBb_[index]); + //return cubeBb.overlaps(allBb_[index]); //- Exact test of triangle intersecting bb // Quick rejection. - if (!cubeBb.intersects(allBb_[index])) + if (!cubeBb.overlaps(allBb_[index])) { return false; } diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C index 6582cfd239..ca4a34975c 100644 --- a/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C +++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C @@ -54,7 +54,7 @@ bool Foam::treeLeaf::findNearest label faceI = indices_[i]; // Quick rejection test. - if (tightest.intersects(allBb[faceI])) + if (tightest.overlaps(allBb[faceI])) { // Full calculation scalar dist = shapes.calcNearest(faceI, sample, nearest); From 48a74535dbcd329b2bb699c9885074dfaeeb82e8 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 28 Aug 2008 15:50:46 +0100 Subject: [PATCH 12/13] added check for non-aligned cell centres --- .../directMappedPolyPatch.C | 90 +++++++++++++++++-- .../directMappedPolyPatch.H | 6 +- 2 files changed, 87 insertions(+), 9 deletions(-) diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C index 83f621363e..e519da59d9 100644 --- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C +++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C @@ -29,6 +29,8 @@ License #include "ListListOps.H" #include "meshSearch.H" #include "mapDistribute.H" +#include "meshTools.H" +#include "OFstream.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -47,16 +49,18 @@ void Foam::directMappedPolyPatch::collectSamples ( pointField& samples, labelList& patchFaceProcs, - labelList& patchFaces + labelList& patchFaces, + pointField& patchFc ) const { - const vectorField::subField fc = this->faceCentres(); // Collect all sample points and the faces they come from. + List globalFc(Pstream::nProcs()); List globalSamples(Pstream::nProcs()); labelListList globalFaces(Pstream::nProcs()); - globalSamples[Pstream::myProcNo()] = fc+offset_; + globalFc[Pstream::myProcNo()] = this->faceCentres(); + globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_; globalFaces[Pstream::myProcNo()] = identity(size()); // Distribute to all processors @@ -64,6 +68,8 @@ void Foam::directMappedPolyPatch::collectSamples Pstream::scatterList(globalSamples); Pstream::gatherList(globalFaces); Pstream::scatterList(globalFaces); + Pstream::gatherList(globalFc); + Pstream::scatterList(globalFc); // Rework into straight list samples = ListListOps::combine @@ -76,6 +82,11 @@ void Foam::directMappedPolyPatch::collectSamples globalFaces, accessOp() ); + patchFc = ListListOps::combine + ( + globalFc, + accessOp() + ); patchFaceProcs.setSize(patchFaces.size()); labelList nPerProc @@ -103,11 +114,14 @@ void Foam::directMappedPolyPatch::findSamples ( const pointField& samples, labelList& sampleCellProcs, - labelList& sampleCells + labelList& sampleCells, + pointField& sampleCc ) const { sampleCellProcs.setSize(samples.size()); sampleCells.setSize(samples.size()); + sampleCc.setSize(samples.size()); + sampleCc = point(-GREAT, -GREAT, -GREAT); { // Octree based search engine @@ -124,6 +138,8 @@ void Foam::directMappedPolyPatch::findSamples else { sampleCellProcs[sampleI] = Pstream::myProcNo(); + sampleCc[sampleI] = + boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]]; } } } @@ -136,6 +152,9 @@ void Foam::directMappedPolyPatch::findSamples Pstream::listCombineGather(sampleCellProcs, maxEqOp