From 77e51249c75e6cbfbeb314af8a6b84cfed86421c Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 18 Feb 2019 13:29:01 +0000 Subject: [PATCH] ENH: distributedTriSurfaceMesh: auto-decomposition; inside/outside support --- .../algorithms/indexedOctree/indexedOctree.H | 6 + .../mapPolyMesh/mapDistribute/mapDistribute.H | 4 +- .../searchableBox/searchableBox.H | 5 + .../searchableCone/searchableCone.H | 6 + .../searchableCylinder/searchableCylinder.H | 6 + .../searchableDisk/searchableDisk.H | 6 + .../searchableExtrudedCircle.H | 6 + .../searchablePlane/searchablePlane.H | 6 + .../searchablePlate/searchablePlate.H | 6 + .../searchableRotatedBox.H | 6 + .../searchableSphere/searchableSphere.H | 6 + .../searchableSurface/searchableSurface.H | 8 +- .../searchableSurfaceCollection.H | 6 + .../searchableSurfaceWithGaps.H | 6 + .../triSurfaceMesh/triSurfaceMesh.C | 381 +- .../triSurfaceMesh/triSurfaceMesh.H | 24 +- .../distributedTriSurfaceMesh.C | 3323 ++++++++++++++--- .../distributedTriSurfaceMesh.H | 136 +- .../simpleFoam/motorBike/Allrun | 3 +- .../distributedTriSurfaceMesh/Allclean | 9 + .../distributedTriSurfaceMesh/Allrun | 53 + .../constant/triSurface/box_12.obj | 34 + .../system/blockMeshDict | 98 + .../system/controlDict | 57 + .../system/decomposeParDict | 27 + .../system/fvSchemes | 57 + .../system/fvSolution | 69 + .../system/meshQualityDict | 21 + .../system/snappyHexMeshDict | 287 ++ 29 files changed, 4142 insertions(+), 520 deletions(-) create mode 100755 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allclean create mode 100755 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allrun create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/constant/triSurface/box_12.obj create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/blockMeshDict create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/controlDict create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/decomposeParDict create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSchemes create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSolution create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/meshQualityDict create mode 100644 tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H index 2557e42b50..e0ccc214c8 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H @@ -471,6 +471,12 @@ public: return nodes_[0].bb_; } + //- Per node, per octant whether is fully inside/outside/mixed. + PackedList<2>& nodeTypes() const + { + return nodeTypes_; + } + // Conversions for entries in subNodes_. diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H index 6dff36ae39..1cca640641 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H @@ -330,10 +330,10 @@ public: mapDistribute(); //- Copy construct - mapDistribute(const mapDistribute& map); + explicit mapDistribute(const mapDistribute& map); //- Move construct - mapDistribute(mapDistribute&& map); + explicit mapDistribute(mapDistribute&& map); //- Move construct from components mapDistribute diff --git a/src/meshTools/searchableSurfaces/searchableBox/searchableBox.H b/src/meshTools/searchableSurfaces/searchableBox/searchableBox.H index b0f0aff690..4d9a7c1759 100644 --- a/src/meshTools/searchableSurfaces/searchableBox/searchableBox.H +++ b/src/meshTools/searchableSurfaces/searchableBox/searchableBox.H @@ -128,6 +128,11 @@ public: { return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } //- Range of local indices that can be returned. virtual label size() const diff --git a/src/meshTools/searchableSurfaces/searchableCone/searchableCone.H b/src/meshTools/searchableSurfaces/searchableCone/searchableCone.H index 9681ea3b38..d4151f806a 100644 --- a/src/meshTools/searchableSurfaces/searchableCone/searchableCone.H +++ b/src/meshTools/searchableSurfaces/searchableCone/searchableCone.H @@ -186,6 +186,12 @@ public: return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableCylinder/searchableCylinder.H b/src/meshTools/searchableSurfaces/searchableCylinder/searchableCylinder.H index 4fd34f197e..8ebdf3a1ef 100644 --- a/src/meshTools/searchableSurfaces/searchableCylinder/searchableCylinder.H +++ b/src/meshTools/searchableSurfaces/searchableCylinder/searchableCylinder.H @@ -158,6 +158,12 @@ public: return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H index 8db94718b6..c5f815be4d 100644 --- a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H +++ b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H @@ -141,6 +141,12 @@ public: return false; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::UNKNOWN; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.H b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.H index 30b354727c..90f3a7af41 100644 --- a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.H +++ b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.H @@ -124,6 +124,12 @@ public: return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } + //- Range of local indices that can be returned. virtual label size() const; diff --git a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H index 85823245bb..099adf0653 100644 --- a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H +++ b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H @@ -126,6 +126,12 @@ public: return false; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::UNKNOWN; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.H b/src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.H index 546ff2a3e9..e4c7bdfdbe 100644 --- a/src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.H +++ b/src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.H @@ -152,6 +152,12 @@ public: return false; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::UNKNOWN; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.H b/src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.H index 46b940b89a..145df14f47 100644 --- a/src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.H +++ b/src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.H @@ -131,6 +131,12 @@ public: return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.H b/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.H index 4d4a8c7c71..505c726e8a 100644 --- a/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.H +++ b/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.H @@ -144,6 +144,12 @@ public: return true; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::OUTSIDE; + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurfaces/searchableSurface/searchableSurface.H index 774cb16932..e1c6e87783 100644 --- a/src/meshTools/searchableSurfaces/searchableSurface/searchableSurface.H +++ b/src/meshTools/searchableSurfaces/searchableSurface/searchableSurface.H @@ -176,13 +176,13 @@ public: } //- Return const reference to boundBox - const boundBox& bounds() const + virtual const boundBox& bounds() const { return bounds_; } //- Return non-const access to the boundBox to allow it to be set. - boundBox& bounds() + virtual boundBox& bounds() { return bounds_; } @@ -194,6 +194,10 @@ public: // This is false for the base class. virtual bool hasVolumeType() const; + //- If surface supports volume queries, what is type of points outside + // bounds + virtual volumeType outsideVolumeType() const = 0; + //- Range of local indices that can be returned virtual label size() const = 0; diff --git a/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.H b/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.H index 4b69f19ae0..0bf2a369c0 100644 --- a/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.H +++ b/src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.H @@ -181,6 +181,12 @@ public: return false; } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return volumeType::UNKNOWN; + } + //- Range of local indices that can be returned. virtual label size() const; diff --git a/src/meshTools/searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.H b/src/meshTools/searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.H index 8b47b79f25..55065dc80e 100644 --- a/src/meshTools/searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.H +++ b/src/meshTools/searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.H @@ -166,6 +166,12 @@ public: return surface().hasVolumeType(); } + //- What is type of points outside bounds + virtual volumeType outsideVolumeType() const + { + return surface().outsideVolumeType(); + } + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C index edcf88d9ec..6da6dc6a7a 100644 --- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C @@ -71,7 +71,7 @@ Foam::fileName Foam::triSurfaceMesh::checkFile Foam::fileName Foam::triSurfaceMesh::relativeFilePath ( - const regIOobject& io, + const IOobject& io, const fileName& f, const bool isGlobal ) @@ -96,7 +96,7 @@ Foam::fileName Foam::triSurfaceMesh::relativeFilePath Foam::fileName Foam::triSurfaceMesh::checkFile ( - const regIOobject& io, + const IOobject& io, const dictionary& dict, const bool isGlobal ) @@ -156,6 +156,13 @@ bool Foam::triSurfaceMesh::addFaceToEdge bool Foam::triSurfaceMesh::isSurfaceClosed() const { + if (debug) + { + Pout<< "triSurfaceMesh::isSurfaceClosed:" + << " determining closedness for surface with " + << triSurface::size() << " triangles" << endl; + } + const pointField& pts = triSurface::points(); // Construct pointFaces. Let's hope surface has compact point @@ -196,6 +203,11 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const if (!okFace) { + if (debug) + { + Pout<< "triSurfaceMesh::isSurfaceClosed :" + << " surface is open" << endl; + } return false; } } @@ -213,6 +225,11 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const if (!okFace) { + if (debug) + { + Pout<< "triSurfaceMesh::isSurfaceClosed :" + << " surface is open" << endl; + } return false; } } @@ -223,11 +240,21 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const { if (iter.val() != 2) { + if (debug) + { + Pout<< "triSurfaceMesh::isSurfaceClosed :" + << " surface is open" << endl; + } return false; } } } + if (debug) + { + Pout<< "triSurfaceMesh::isSurfaceClosed :" + << " surface is closed" << endl; + } return true; } @@ -358,7 +385,7 @@ Foam::triSurfaceMesh::triSurfaceMesh } -Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const bool isGlobal) +Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const readAction r) : // Find instance for triSurfaceMesh searchableSurface(io), @@ -376,18 +403,79 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const bool isGlobal) false // searchableSurface already registered under name ) ), - triSurface - ( - checkFile(static_cast(*this), isGlobal) - ), + triSurface(), // construct null triSurfaceRegionSearch(static_cast(*this)), minQuality_(-1), surfaceClosed_(-1), outsideVolType_(volumeType::UNKNOWN) { - const pointField& pts = triSurface::points(); + // Check IO flags + if (io.readOpt() != IOobject::NO_READ) + { + const bool searchGlobal(r == localOrGlobal || r == masterOnly); - bounds() = boundBox(pts, isGlobal); + const fileName actualFile + ( + searchGlobal + ? io.globalFilePath(typeName) + : io.localFilePath(typeName) + ); + + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loading surface " << io.objectPath() + << " local filePath:" << io.localFilePath(typeName) + << " from:" << actualFile << endl; + } + + if (searchGlobal && Pstream::parRun()) + { + // Check where surface was found + const fileName localFile(io.localFilePath(typeName)); + + if (r == masterOnly && (actualFile != localFile)) + { + // Found undecomposed surface. Load on master only + if (Pstream::master()) + { + triSurface s2(actualFile); + triSurface::transfer(s2); + } + Pstream::scatter(triSurface::patches()); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } + else + { + // Read on all processors + triSurface s2(actualFile); + triSurface::transfer(s2); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } + } + else + { + // Read on all processors + triSurface s2(actualFile); + triSurface::transfer(s2); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } + } + + const pointField& pts = triSurface::points(); + bounds() = boundBox(pts, false); } @@ -395,7 +483,7 @@ Foam::triSurfaceMesh::triSurfaceMesh ( const IOobject& io, const dictionary& dict, - const bool isGlobal + const readAction r ) : searchableSurface(io), @@ -413,26 +501,91 @@ Foam::triSurfaceMesh::triSurfaceMesh false // searchableSurface already registered under name ) ), - triSurface - ( - checkFile(static_cast(*this), dict, isGlobal) - ), + triSurface(), // construct null triSurfaceRegionSearch(static_cast(*this), dict), minQuality_(-1), surfaceClosed_(-1), outsideVolType_(volumeType::UNKNOWN) { - // Reading from supplied file name instead of objectPath/filePath - if (dict.readIfPresent("file", fName_, keyType::LITERAL)) + // Check IO flags + if (io.readOpt() != IOobject::NO_READ) { - fName_ = relativeFilePath + const bool searchGlobal(r == localOrGlobal || r == masterOnly); + + fileName actualFile ( - static_cast(*this), - fName_, - isGlobal + searchGlobal + ? io.globalFilePath(typeName) + : io.localFilePath(typeName) ); + + // Reading from supplied file name instead of objectPath/filePath + if (dict.readIfPresent("file", fName_, keyType::LITERAL)) + { + fName_ = relativeFilePath + ( + static_cast(*this), + fName_, + searchGlobal + ); + actualFile = fName_; + } + + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io, const dictionary&) :" + << " loading surface " << io.objectPath() + << " local filePath:" << io.localFilePath(typeName) + << " from:" << actualFile << endl; + } + + if (searchGlobal && Pstream::parRun()) + { + // Check where surface was found + const fileName localFile(io.localFilePath(typeName)); + + if (r == masterOnly && (actualFile != localFile)) + { + // Surface not loaded from processor directories -> undecomposed + // surface. Load on master only + if (Pstream::master()) + { + triSurface s2(actualFile); + triSurface::transfer(s2); + } + Pstream::scatter(triSurface::patches()); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } + else + { + // Read on all processors + triSurface s2(actualFile); + triSurface::transfer(s2); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } + } + else + { + // Read on all processors + triSurface s2(actualFile); + triSurface::transfer(s2); + if (debug) + { + Pout<< "triSurfaceMesh(const IOobject& io) :" + << " loaded triangles:" << triSurface::size() << endl; + } + } } + scalar scaleFactor = 0; // Allow rescaling of the surface points @@ -445,8 +598,7 @@ Foam::triSurfaceMesh::triSurfaceMesh } const pointField& pts = triSurface::points(); - - bounds() = boundBox(pts, isGlobal); + bounds() = boundBox(pts, false); // Have optional minimum quality for normal calculation if (dict.readIfPresent("minQuality", minQuality_) && minQuality_ > 0) @@ -468,7 +620,7 @@ Foam::triSurfaceMesh::~triSurfaceMesh() void Foam::triSurfaceMesh::clearOut() { - outsideVolType_ = volumeType::UNKNOWN; + // Do not clear closedness status triSurfaceRegionSearch::clearOut(); edgeTree_.clear(); triSurface::clearOut(); @@ -539,7 +691,14 @@ bool Foam::triSurfaceMesh::overlaps(const boundBox& bb) const void Foam::triSurfaceMesh::movePoints(const pointField& newPoints) { - outsideVolType_ = volumeType::UNKNOWN; + if (debug) + { + Pout<< "triSurfaceMesh::movePoints :" + << " moving at time " << objectRegistry::time().timeName() + << endl; + } + + // Preserve topological point status (surfaceClosed_, outsideVolType_) // Update local information (instance, event number) searchableSurface::instance() = objectRegistry::time().timeName(); @@ -555,6 +714,10 @@ void Foam::triSurfaceMesh::movePoints(const pointField& newPoints) triSurface::movePoints(newPoints); bounds() = boundBox(triSurface::points(), false); + if (debug) + { + Pout<< "triSurfaceMesh::movePoints: finished moving points" << endl; + } } @@ -563,6 +726,13 @@ Foam::triSurfaceMesh::edgeTree() const { if (edgeTree_.empty()) { + if (debug) + { + Pout<< "triSurfaceMesh::edgeTree :" + << " constructing tree for " << nEdges() - nInternalEdges() + << " boundary edges" << endl; + } + // Boundary edges labelList bEdges ( @@ -592,6 +762,13 @@ Foam::triSurfaceMesh::edgeTree() const bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); } + + if (debug) + { + Pout<< "triSurfaceMesh::edgeTree : " + << "calculating edge tree for bb:" << bb << endl; + } + scalar oldTol = indexedOctree::perturbTol(); indexedOctree::perturbTol() = tolerance(); @@ -614,6 +791,14 @@ Foam::triSurfaceMesh::edgeTree() const ); indexedOctree::perturbTol() = oldTol; + + if (debug) + { + Pout<< "triSurfaceMesh::edgeTree :" + << " finished constructing tree for " + << nEdges() - nInternalEdges() + << " boundary edges" << endl; + } } return *edgeTree_; @@ -652,6 +837,40 @@ bool Foam::triSurfaceMesh::hasVolumeType() const } +Foam::volumeType Foam::triSurfaceMesh::outsideVolumeType() const +{ + if (outsideVolType_ == volumeType::UNKNOWN) + { + // Get point outside bounds() + const point outsidePt(bounds().max() + 0.5*bounds().span()); + + if (debug) + { + Pout<< "triSurfaceMesh::outsideVolumeType :" + << " triggering outsidePoint:" << outsidePt + << " orientation" << endl; + } + + //outsideVolType_ = tree().shapes().getVolumeType(tree(), outsidePt); + // Note: do not use tree directly so e.g. distributedTriSurfaceMesh + // has opportunity to intercept + List outsideVolTypes; + getVolumeType(pointField(1, outsidePt), outsideVolTypes); + outsideVolType_ = outsideVolTypes[0]; + + if (debug) + { + Pout<< "triSurfaceMesh::outsideVolumeType :" + << " finished outsidePoint:" << outsidePt + << " orientation:" << volumeType::names[outsideVolType_] + << endl; + } + } + + return outsideVolType_; +} + + void Foam::triSurfaceMesh::findNearest ( const pointField& samples, @@ -659,7 +878,21 @@ void Foam::triSurfaceMesh::findNearest List& info ) const { + if (debug) + { + Pout<< "triSurfaceMesh::findNearest :" + << " trying to find nearest for " << samples.size() + << " samples with max sphere " + << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero) + << endl; + } triSurfaceSearch::findNearest(samples, nearestDistSqr, info); + if (debug) + { + Pout<< "triSurfaceMesh::findNearest :" + << " finished trying to find nearest for " << samples.size() + << " samples" << endl; + } } @@ -671,6 +904,14 @@ void Foam::triSurfaceMesh::findNearest List& info ) const { + if (debug) + { + Pout<< "triSurfaceMesh::findNearest :" + << " trying to find nearest and region for " << samples.size() + << " samples with max sphere " + << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero) + << endl; + } triSurfaceRegionSearch::findNearest ( samples, @@ -678,6 +919,12 @@ void Foam::triSurfaceMesh::findNearest regionIndices, info ); + if (debug) + { + Pout<< "triSurfaceMesh::findNearest :" + << " finished trying to find nearest and region for " + << samples.size() << " samples" << endl; + } } @@ -688,7 +935,19 @@ void Foam::triSurfaceMesh::findLine List& info ) const { + if (debug) + { + Pout<< "triSurfaceMesh::findLine :" + << " intersecting with " + << start.size() << " rays" << endl; + } triSurfaceSearch::findLine(start, end, info); + if (debug) + { + Pout<< "triSurfaceMesh::findLine :" + << " finished intersecting with " + << start.size() << " rays" << endl; + } } @@ -699,7 +958,19 @@ void Foam::triSurfaceMesh::findLineAny List& info ) const { + if (debug) + { + Pout<< "triSurfaceMesh::findLineAny :" + << " intersecting with " + << start.size() << " rays" << endl; + } triSurfaceSearch::findLineAny(start, end, info); + if (debug) + { + Pout<< "triSurfaceMesh::findLineAny :" + << " finished intersecting with " + << start.size() << " rays" << endl; + } } @@ -710,7 +981,19 @@ void Foam::triSurfaceMesh::findLineAll List>& info ) const { + if (debug) + { + Pout<< "triSurfaceMesh::findLineAll :" + << " intersecting with " + << start.size() << " rays" << endl; + } triSurfaceSearch::findLineAll(start, end, info); + if (debug) + { + Pout<< "triSurfaceMesh::findLineAll :" + << " finished intersecting with " + << start.size() << " rays" << endl; + } } @@ -720,6 +1003,12 @@ void Foam::triSurfaceMesh::getRegion labelList& region ) const { + if (debug) + { + Pout<< "triSurfaceMesh::getRegion :" + << " getting region for " + << info.size() << " triangles" << endl; + } region.setSize(info.size()); forAll(info, i) { @@ -732,6 +1021,12 @@ void Foam::triSurfaceMesh::getRegion region[i] = -1; } } + if (debug) + { + Pout<< "triSurfaceMesh::getRegion :" + << " finished getting region for " + << info.size() << " triangles" << endl; + } } @@ -741,6 +1036,13 @@ void Foam::triSurfaceMesh::getNormal vectorField& normal ) const { + if (debug) + { + Pout<< "triSurfaceMesh::getNormal :" + << " getting normal for " + << info.size() << " triangles" << endl; + } + const triSurface& s = *this; const pointField& pts = s.points(); @@ -803,6 +1105,12 @@ void Foam::triSurfaceMesh::getNormal } } } + if (debug) + { + Pout<< "triSurfaceMesh::getNormal :" + << " finished getting normal for " + << info.size() << " triangles" << endl; + } } @@ -835,6 +1143,12 @@ void Foam::triSurfaceMesh::setField(const labelList& values) // Store field on triMesh fldPtr->store(); } + if (debug) + { + Pout<< "triSurfaceMesh::setField :" + << " finished setting field for " + << values.size() << " triangles" << endl; + } } @@ -860,6 +1174,12 @@ void Foam::triSurfaceMesh::getField } } } + if (debug) + { + Pout<< "triSurfaceMesh::setField :" + << " finished getting field for " + << info.size() << " triangles" << endl; + } } @@ -872,6 +1192,13 @@ void Foam::triSurfaceMesh::getVolumeType const scalar oldTol = indexedOctree::perturbTol(); indexedOctree::perturbTol() = tolerance(); + if (debug) + { + Pout<< "triSurfaceMesh::getVolumeType :" + << " finding orientation for " << points.size() + << " samples" << endl; + } + volType.setSize(points.size()); forAll(points, pointi) @@ -900,6 +1227,12 @@ void Foam::triSurfaceMesh::getVolumeType } indexedOctree::perturbTol() = oldTol; + if (debug) + { + Pout<< "triSurfaceMesh::getVolumeType :" + << " finished finding orientation for " << points.size() + << " samples" << endl; + } } diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H index 006c4b9189..39ea790307 100644 --- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H @@ -52,8 +52,8 @@ SourceFiles #ifndef triSurfaceMesh_H #define triSurfaceMesh_H -#include "treeBoundBox.H" #include "searchableSurface.H" +#include "treeBoundBox.H" #include "objectRegistry.H" #include "indexedOctree.H" #include "treeDataTriSurface.H" @@ -79,6 +79,8 @@ class triSurfaceMesh public triSurface, public triSurfaceRegionSearch { +protected: + // Private member data //- Supplied fileName override @@ -110,7 +112,7 @@ class triSurfaceMesh // IOobject static fileName relativeFilePath ( - const regIOobject&, + const IOobject&, const fileName&, const bool isGlobal ); @@ -118,7 +120,7 @@ class triSurfaceMesh //- Return fileName to load IOobject from. Optional override of fileName static fileName checkFile ( - const regIOobject&, + const IOobject&, const dictionary&, const bool isGlobal ); @@ -179,15 +181,22 @@ public: // Special constructors for use by distributedTriSurface. File search - // status (local/global) supplied. + // method supplied: - triSurfaceMesh(const IOobject& io, const bool isGlobal); + enum readAction + { + localOnly, // load from (processor-)local file only + localOrGlobal, // load from (processor-)local or global file + masterOnly // as localOrGlobal but only load on master + }; + + triSurfaceMesh(const IOobject& io, const readAction r); triSurfaceMesh ( const IOobject& io, const dictionary& dict, - const bool isGlobal + const readAction r ); @@ -216,6 +225,9 @@ public: //- Whether supports volume type (below) - i.e. whether is closed. virtual bool hasVolumeType() const; + //- If surface is closed, what is type of points outside bounds + virtual volumeType outsideVolumeType() const; + //- Range of local indices that can be returned. virtual label size() const { diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C index 59e21e39b0..9b3cf90221 100644 --- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C +++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C @@ -41,6 +41,8 @@ License #include "bitSet.H" #include "PatchTools.H" #include "OBJstream.H" +#include "labelBits.H" +#include "profiling.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -53,6 +55,89 @@ namespace Foam distributedTriSurfaceMesh, dict ); + + + //- Combine operator for volume types + class volumeCombineOp + { + public: + void operator()(volumeType& x, const volumeType& y) const + { + if (x == volumeType::MIXED || y == volumeType::MIXED) + { + FatalErrorInFunction << "Illegal volume type " + << volumeType::names[x] + << "," << volumeType::names[y] << exit(FatalError); + } + else + { + switch (x) + { + case volumeType::UNKNOWN: + { + if (y == volumeType::INSIDE || y == volumeType::OUTSIDE) + { + x = y; + } + } + break; + case volumeType::INSIDE: + { + if (y == volumeType::OUTSIDE) + { + FatalErrorInFunction << "Conflicting volume types " + << volumeType::names[x] << "," + << volumeType::names[y] << exit(FatalError); + } + } + break; + case volumeType::OUTSIDE: + { + if (y == volumeType::INSIDE) + { + FatalErrorInFunction << "Conflicting volume types " + << volumeType::names[x] << "," + << volumeType::names[y] << exit(FatalError); + } + } + break; + case volumeType::MIXED: + break; + } + } + } + }; + + + //- Combine operator for nearest + typedef Tuple2 nearestAndDist; + + const nearestAndDist nearestZero + ( + nearestAndDist + ( + pointIndexHit(), + -GREAT + ) + ); + class nearestEqOp + { + public: + void operator()(nearestAndDist& x, const nearestAndDist& y) const + { + if (x.first().hit()) + { + if (y.first().hit() && y.second() < x.second()) + { + x = y; + } + } + else if (y.first().hit()) + { + x = y; + } + } + }; } @@ -64,27 +149,181 @@ Foam::distributedTriSurfaceMesh::distributionTypeNames_ ({ { distributionType::FOLLOW, "follow" }, { distributionType::INDEPENDENT, "independent" }, - { distributionType::FROZEN, "frozen" }, + { distributionType::DISTRIBUTED, "distributed" }, + { distributionType::FROZEN, "frozen" } }); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +Foam::word Foam::distributedTriSurfaceMesh::findLocalInstance +( + const IOobject& io +) +{ + // Modified findInstance which also looks in parent directory + word instance + ( + io.time().findInstance + ( + io.local(), + word::null, + IOobject::READ_IF_PRESENT + ) + ); + + if (instance.size()) + { + return instance; + } + + + // Replicate findInstance operation but now on parent directory + + // Search in parent directory + fileName parentDir = + io.rootPath()/io.time().globalCaseName() + /io.instance()/io.db().dbDir()/io.local()/io.name(); + + if (fileHandler().isDir(parentDir)) + { + return io.instance(); + } + + instantList ts = io.time().times(); + label instanceI; + + const scalar startValue = io.time().timeOutputValue(); + + for (instanceI = ts.size()-1; instanceI >= 0; --instanceI) + { + if (ts[instanceI].value() <= startValue) + { + break; + } + } + + // continue searching from here + for (; instanceI >= 0; --instanceI) + { + // Shortcut: if actual directory is the timeName we've already tested it + if (ts[instanceI].name() == io.instance()) + { + continue; + } + + fileName parentDir = + io.rootPath()/io.time().globalCaseName() + /ts[instanceI].name()/io.db().dbDir()/io.local()/io.name(); + + if (fileHandler().isDir(parentDir)) + { + return ts[instanceI].name(); + } + } + + // times() usually already includes the constant() so would have been + // checked above. Re-test if + // - times() is empty. Sometimes this can happen (e.g. decomposePar with + // collated) + // - times()[0] is not constant + if (!ts.size() || ts[0].name() != io.time().constant()) + { + // Note. This needs to be a hard-coded constant, rather than the + // constant function of the time, because the latter points to + // the case constant directory in parallel cases + + fileName parentDir = + io.rootPath()/io.time().globalCaseName() + /io.time().constant()/io.db().dbDir()/io.local()/io.name(); + + if (fileHandler().isDir(parentDir)) + { + return io.time().constant(); + } + } + + FatalErrorInFunction + << "Cannot find directory " << io.local() << " in times " << ts + << exit(FatalError); + + return word::null; +} + + // Read my additional data from the dictionary bool Foam::distributedTriSurfaceMesh::read() { // Get bb of all domains. procBb_.setSize(Pstream::nProcs()); - procBb_[Pstream::myProcNo()] = List(dict_.lookup("bounds")); - Pstream::gatherList(procBb_); - Pstream::scatterList(procBb_); + if (dict_.empty()) + { + // Did not find the distributed version; assume master has loaded the + // triSurfaceMesh version. Make up some settings. - // Distribution type - distType_ = distributionTypeNames_.get("distributionType", dict_); + const boundBox& localBb = triSurfaceMesh::bounds(); - // Merge distance - dict_.readEntry("mergeDistance", mergeDist_); + procBb_[Pstream::myProcNo()] = + treeBoundBoxList(1, treeBoundBox(localBb)); + Pstream::gatherList(procBb_); + Pstream::scatterList(procBb_); + + dict_.add("bounds", procBb_[Pstream::myProcNo()]); + + // Wanted distribution type + distType_ = DISTRIBUTED; //INDEPENDENT; + dict_.add("distributionType", distributionTypeNames_[distType_]); + + // Merge distance + mergeDist_ = SMALL; + dict_.add("mergeDistance", mergeDist_); + + // Force underlying triSurfaceMesh to calculate volume type + // (is topological walk; does not construct tree) + surfaceClosed_ = triSurfaceMesh::hasVolumeType(); + Pstream::scatter(surfaceClosed_); + dict_.add("closed", surfaceClosed_); + + // Delay calculating outside vol type since constructs tree. Is ok + // after distributing since then local surfaces much smaller + //outsideVolType_ = volumeType::UNKNOWN; + //if (surfaceClosed_) + //{ + // point outsidePt(localBb.max()+localBb.midpoint()); + // List outsideVolTypes; + // triSurfaceMesh::getVolumeType + // ( + // pointField(1, outsidePt), + // outsideVolTypes + // ); + // outsideVolType_ = outsideVolTypes[0]; + //} + //dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]); + } + else + { + procBb_[Pstream::myProcNo()] = + List(dict_.lookup("bounds")); + Pstream::gatherList(procBb_); + Pstream::scatterList(procBb_); + + // Wanted distribution type + distType_ = distributionTypeNames_.lookup("distributionType", dict_); + + // Merge distance + mergeDist_ = readScalar(dict_.lookup("mergeDistance")); + + // Distribution type + surfaceClosed_ = dict_.lookupOrDefault("closed", false); + + outsideVolType_ = volumeType::UNKNOWN; + word volType; + if (dict_.readIfPresent("outsideVolumeType", volType)) + { + outsideVolType_ = volumeType::names[volType]; + } + } return true; } @@ -220,7 +459,7 @@ void Foam::distributedTriSurfaceMesh::distributeSegment const treeBoundBox& bb = bbs[bbi]; // Scheme a: any processor that intersects the segment gets - // the segment. + // the whole segment. // Intersection point point clipPt; @@ -278,6 +517,17 @@ Foam::distributedTriSurfaceMesh::distributeSegments // Per processor indices into allSegments to send List> dynSendMap(Pstream::nProcs()); + // Pre-size + forAll(dynSendMap, proci) + { + dynSendMap[proci].reserve + ( + (proci == Pstream::myProcNo()) + ? start.size() + : start.size()/Pstream::nProcs() + ); + } + forAll(start, segmenti) { distributeSegment @@ -304,49 +554,7 @@ Foam::distributedTriSurfaceMesh::distributeSegments allSegmentMap.transfer(dynAllSegmentMap); } - - // Send over how many i need to receive. - labelListList sendSizes(Pstream::nProcs()); - sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs()); - forAll(sendMap, proci) - { - sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size(); - } - Pstream::gatherList(sendSizes); - Pstream::scatterList(sendSizes); - - - // Determine order of receiving - labelListList constructMap(Pstream::nProcs()); - - // My local segments first - constructMap[Pstream::myProcNo()] = identity - ( - sendMap[Pstream::myProcNo()].size() - ); - - label segmenti = constructMap[Pstream::myProcNo()].size(); - forAll(constructMap, proci) - { - if (proci != Pstream::myProcNo()) - { - // What i need to receive is what other processor is sending to me. - label nRecv = sendSizes[proci][Pstream::myProcNo()]; - constructMap[proci].setSize(nRecv); - - for (label i = 0; i < nRecv; i++) - { - constructMap[proci][i] = segmenti++; - } - } - } - - return autoPtr::New - ( - segmenti, // size after construction - std::move(sendMap), - std::move(constructMap) - ); + return autoPtr::New(std::move(sendMap)); } @@ -358,6 +566,14 @@ void Foam::distributedTriSurfaceMesh::findLine List& info ) const { + if (debug) + { + Pout<< "distributedTriSurfaceMesh::findLine :" + << " intersecting with " + << start.size() << " rays" << endl; + } + addProfiling(findLine, "distributedTriSurfaceMesh::findLine"); + const indexedOctree& octree = tree(); // Initialise @@ -367,9 +583,18 @@ void Foam::distributedTriSurfaceMesh::findLine info[i].setMiss(); } - if (!Pstream::parRun()) + // Important:force synchronised construction of indexing + const globalIndex& triIndexer = globalTris(); + + + // Do any local queries + // ~~~~~~~~~~~~~~~~~~~~ + + label nLocal = 0; + + forAll(start, i) { - forAll(start, i) + if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i])) { if (nearestIntersection) { @@ -379,150 +604,124 @@ void Foam::distributedTriSurfaceMesh::findLine { info[i] = octree.findLineAny(start[i], end[i]); } + + if (info[i].hit()) + { + info[i].setIndex(triIndexer.toGlobal(info[i].index())); + } + nLocal++; } } - else + + + if + ( + returnReduce(nLocal, sumOp