diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
index 5aa2a3d5673f87eb43672e3af992a4b823d82aa4..8b0e74a14b1fddfc9785f544b7a83f29c28d87a5 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 209acdec40fa97eb02aff6210de4e059aaf20393..111a0633351e36bf7a18eca1ed11c21f1a9367f5 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 25c8b19b625f95c2da3707679fb45dee7805024c..40bd7ac32f73cee1d97146c2e6578fa0ddad30ad 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 9681ea3b383b59fe8a20558295868d244a6c8129..d4151f806a9f1ef7d2aad9f9eba5e8d7e30b0d53 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 5cae9b1b09139aa3f938748734f1a8a265697010..ada4f48f11d0e10988c8ca9c012f0dd6176ea824 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 8db94718b662a9fa330e8d962142c63d08c64ab0..c5f815be4d5a3cd6ad2a671f43eb0ef891976c0b 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 30b354727c86d4e770f6c4c6d8fc49ae1df5b54f..90f3a7af41467b903fb421d9a9adff59385f1775 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 18733059fad47a5e38088ced398818722b0d0457..d0b5d3f444e56051c363eaa0e43e079a5f6c0fbf 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 ea64d05e2f366d9595cfce4b57fb27f735fbc729..1deb01edbc8b49b3a297035870ec8b002865a86d 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 46b940b89aae9dfcf9ada2ae557033cb421e3604..145df14f47e0083ab5155319e7c67ddcd2c1a3a3 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 096ea09d6effd63d2457060bef366d788f16b2fd..bac39bf08677af486e9536e045fcaf89680d9e7d 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 0e995c62b598549070013ce8504ff1c3fbf9fd67..7ce9b4df3c389d749d2bf654d49348dffcdfd34d 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 e728e71d65577c54145fb3daf8e748c4ec089acf..97b72ec2d7e13cd9dd4d51cb397f2166c74bba39 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 b93f3882674ee80661ddb3f1f58cd4493c09e4a1..7177a1c2a9448f8a667df363dff42fcab179c00c 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 56e3520564ce610bd7a3040df45c1a3b95293aa6..e0677dd0e0ad6dbbac8643a0a9c0ffaa298c5bff 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<const searchableSurface&>(*this), isGlobal)
- ),
+ triSurface(), // construct null
triSurfaceRegionSearch(static_cast<const triSurface&>(*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);
+
+ 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;
+ }
+ }
+ }
- bounds() = boundBox(pts, isGlobal);
+ 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<const searchableSurface&>(*this), dict, isGlobal)
- ),
+ triSurface(), // construct null
triSurfaceRegionSearch(static_cast<const triSurface&>(*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<const searchableSurface&>(*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<const searchableSurface&>(*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<treeDataEdge>::perturbTol();
indexedOctree<treeDataEdge>::perturbTol() = tolerance();
@@ -614,6 +791,14 @@ Foam::triSurfaceMesh::edgeTree() const
);
indexedOctree<treeDataEdge>::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<volumeType> 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<pointIndexHit>& 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<pointIndexHit>& 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<pointIndexHit>& 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<pointIndexHit>& 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<List<pointIndexHit>>& 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<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::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<treeDataTriSurface>::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 abce60f63836895977779a2b764ab0692997174d..b07da194aef75567e8b4237b76e7e3a9b6a55a46 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 655c934f7ade937259895913c5a71db98e14091a..3613206a84bfd32575963d3f538482f64b5647c2 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<pointIndexHit, scalar> 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<treeBoundBox>(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.
+
+ const boundBox& localBb = triSurfaceMesh::bounds();
+
+ 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<volumeType> outsideVolTypes;
+ // triSurfaceMesh::getVolumeType
+ // (
+ // pointField(1, outsidePt),
+ // outsideVolTypes
+ // );
+ // outsideVolType_ = outsideVolTypes[0];
+ //}
+ //dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
+ }
+ else
+ {
+ procBb_[Pstream::myProcNo()] =
+ List<treeBoundBox>(dict_.lookup("bounds"));
+ Pstream::gatherList(procBb_);
+ Pstream::scatterList(procBb_);
- // Distribution type
- distType_ = distributionTypeNames_.get("distributionType", dict_);
+ // Wanted distribution type
+ distType_ = distributionTypeNames_.lookup("distributionType", dict_);
- // Merge distance
- dict_.readEntry("mergeDistance", mergeDist_);
+ // Merge distance
+ mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
+
+ // Distribution type
+ surfaceClosed_ = dict_.lookupOrDefault<bool>("closed", false);
+
+ outsideVolType_ = volumeType::UNKNOWN;
+ word volType;
+ if (dict_.readIfPresent<word>("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<DynamicList<label>> 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<mapDistribute>::New
- (
- segmenti, // size after construction
- std::move(sendMap),
- std::move(constructMap)
- );
+ return autoPtr<mapDistribute>::New(std::move(sendMap));
}
@@ -358,6 +566,14 @@ void Foam::distributedTriSurfaceMesh::findLine
List<pointIndexHit>& info
) const
{
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::findLine :"
+ << " intersecting with "
+ << start.size() << " rays" << endl;
+ }
+ addProfiling(findLine, "distributedTriSurfaceMesh::findLine");
+
const indexedOctree<treeDataTriSurface>& 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]);
}
- }
- }
- else
- {
- // Important:force synchronised construction of indexing
- const globalIndex& triIndexer = globalTris();
-
- // Do any local queries
- // ~~~~~~~~~~~~~~~~~~~~
-
- label nLocal = 0;
-
- forAll(start, i)
- {
- if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
+ if (info[i].hit())
{
- if (nearestIntersection)
- {
- info[i] = octree.findLine(start[i], end[i]);
- }
- else
- {
- info[i] = octree.findLineAny(start[i], end[i]);
- }
-
- if (info[i].hit())
- {
- info[i].setIndex(triIndexer.toGlobal(info[i].index()));
- }
- nLocal++;
+ info[i].setIndex(triIndexer.toGlobal(info[i].index()));
}
+ nLocal++;
}
+ }
- if
- (
- returnReduce(nLocal, sumOp<label>())
- < returnReduce(start.size(), sumOp<label>())
- )
- {
- // Not all can be resolved locally. Build segments and map,
- // send over segments, do intersections, send back and merge.
+ if
+ (
+ returnReduce(nLocal, sumOp<label>())
+ < returnReduce(start.size(), sumOp<label>())
+ )
+ {
+ // Not all can be resolved locally. Build segments and map,
+ // send over segments, do intersections, send back and merge.
- // Construct queries (segments)
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Construct queries (segments)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Segments to test
- List<segment> allSegments(start.size());
- // Original index of segment
- labelList allSegmentMap(start.size());
+ // Segments to test
+ List<segment> allSegments(start.size());
+ // Original index of segment
+ labelList allSegmentMap(start.size());
- const autoPtr<mapDistribute> mapPtr
+ const autoPtr<mapDistribute> mapPtr
+ (
+ distributeSegments
(
- distributeSegments
- (
- start,
- end,
- allSegments,
- allSegmentMap
- )
- );
- const mapDistribute& map = mapPtr();
+ start,
+ end,
+ allSegments,
+ allSegmentMap
+ )
+ );
+ const mapDistribute& map = mapPtr();
- label nOldAllSegments = allSegments.size();
+ label nOldAllSegments = allSegments.size();
- // Exchange the segments
- // ~~~~~~~~~~~~~~~~~~~~~
+ // Exchange the segments
+ // ~~~~~~~~~~~~~~~~~~~~~
- map.distribute(allSegments);
+ map.distribute(allSegments);
- // Do tests i need to do
- // ~~~~~~~~~~~~~~~~~~~~~
+ // Do tests i need to do
+ // ~~~~~~~~~~~~~~~~~~~~~
- // Intersections
- List<pointIndexHit> intersections(allSegments.size());
+ // Intersections
+ List<pointIndexHit> intersections(allSegments.size());
- forAll(allSegments, i)
+ forAll(allSegments, i)
+ {
+ if (nearestIntersection)
{
- if (nearestIntersection)
- {
- intersections[i] = octree.findLine
- (
- allSegments[i].first(),
- allSegments[i].second()
- );
- }
- else
- {
- intersections[i] = octree.findLineAny
- (
- allSegments[i].first(),
- allSegments[i].second()
- );
- }
+ intersections[i] = octree.findLine
+ (
+ allSegments[i].first(),
+ allSegments[i].second()
+ );
+ }
+ else
+ {
+ intersections[i] = octree.findLineAny
+ (
+ allSegments[i].first(),
+ allSegments[i].second()
+ );
+ }
- // Convert triangle index to global numbering
- if (intersections[i].hit())
- {
- intersections[i].setIndex
- (
- triIndexer.toGlobal(intersections[i].index())
- );
- }
+ // Convert triangle index to global numbering
+ if (intersections[i].hit())
+ {
+ intersections[i].setIndex
+ (
+ triIndexer.toGlobal(intersections[i].index())
+ );
}
+ }
- // Exchange the intersections (opposite to segments)
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Exchange the intersections (opposite to segments)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- map.reverseDistribute(nOldAllSegments, intersections);
+ map.reverseDistribute(nOldAllSegments, intersections);
- // Extract the hits
- // ~~~~~~~~~~~~~~~~
+ // Extract the hits
+ // ~~~~~~~~~~~~~~~~
- forAll(intersections, i)
- {
- const pointIndexHit& allInfo = intersections[i];
- label segmenti = allSegmentMap[i];
- pointIndexHit& hitInfo = info[segmenti];
+ forAll(intersections, i)
+ {
+ const pointIndexHit& allInfo = intersections[i];
+ label segmenti = allSegmentMap[i];
+ pointIndexHit& hitInfo = info[segmenti];
- if (allInfo.hit())
+ if (allInfo.hit())
+ {
+ if (!hitInfo.hit())
+ {
+ // No intersection yet so take this one
+ hitInfo = allInfo;
+ }
+ else if (nearestIntersection)
{
- if (!hitInfo.hit())
+ // Nearest intersection
+ if
+ (
+ magSqr(allInfo.hitPoint()-start[segmenti])
+ < magSqr(hitInfo.hitPoint()-start[segmenti])
+ )
{
- // No intersection yet so take this one
hitInfo = allInfo;
}
- else if (nearestIntersection)
- {
- // Nearest intersection
- if
- (
- magSqr(allInfo.hitPoint()-start[segmenti])
- < magSqr(hitInfo.hitPoint()-start[segmenti])
- )
- {
- hitInfo = allInfo;
- }
- }
}
}
}
@@ -530,6 +729,24 @@ void Foam::distributedTriSurfaceMesh::findLine
}
+void Foam::distributedTriSurfaceMesh::convertTriIndices
+(
+ List<pointIndexHit>& info
+) const
+{
+ // Important:force synchronised construction of indexing
+ const globalIndex& triIndexer = globalTris();
+
+ for (pointIndexHit& pi : info)
+ {
+ if (pi.hit())
+ {
+ pi.setIndex(triIndexer.toGlobal(pi.index()));
+ }
+ }
+}
+
+
// Exchanges indices to the processor they come from.
// - calculates exchange map
// - uses map to calculate local triangle index
@@ -540,6 +757,10 @@ Foam::distributedTriSurfaceMesh::calcLocalQueries
labelList& triangleIndex
) const
{
+ // Note: does not filter duplicate queries so a triangle that gets requested
+ // from more than one processor also get local queried more than
+ // once.
+
triangleIndex.setSize(info.size());
const globalIndex& triIndexer = globalTris();
@@ -587,68 +808,147 @@ Foam::distributedTriSurfaceMesh::calcLocalQueries
}
- // Send over how many i need to receive
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Pack into distribution map
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~
- 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);
+ autoPtr<mapDistribute> mapPtr(new mapDistribute(std::move(sendMap)));
- // Determine receive map
- // ~~~~~~~~~~~~~~~~~~~~~
+ // Send over queries
+ // ~~~~~~~~~~~~~~~~~
- labelListList constructMap(Pstream::nProcs());
+ mapPtr().distribute(triangleIndex);
+
+ return mapPtr;
+}
- // My local segments first
- constructMap[Pstream::myProcNo()] = identity
- (
- sendMap[Pstream::myProcNo()].size()
- );
- label segmenti = constructMap[Pstream::myProcNo()].size();
- forAll(constructMap, proci)
+bool Foam::distributedTriSurfaceMesh::contains
+(
+ const List<treeBoundBox>& bbs,
+ const point& sample
+) const
+{
+ forAll(bbs, bbi)
{
- if (proci != Pstream::myProcNo())
+ if (bbs[bbi].contains(sample))
{
- // 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 true;
}
}
+ return false;
+}
- // Pack into distribution map
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+Foam::Tuple2<Foam::label, Foam::scalar>
+Foam::distributedTriSurfaceMesh::findBestProcs
+(
+ const point& centre,
+ const scalar radiusSqr,
+ boolList& procContains,
+ boolList& procOverlaps,
+ label& minProci
+) const
+{
+ // Find processors:
+ // - that contain the centre
+ // - or overlap the search sphere
- autoPtr<mapDistribute> mapPtr
- (
- new mapDistribute
- (
- segmenti, // size after construction
- std::move(sendMap),
- std::move(constructMap)
- )
- );
- const mapDistribute& map = mapPtr();
+ procContains.setSize(Pstream::nProcs());
+ procContains = false;
+ procOverlaps.setSize(Pstream::nProcs());
+ procOverlaps = false;
- // Send over queries
- // ~~~~~~~~~~~~~~~~~
+ minProci = -1;
- map.distribute(triangleIndex);
+ scalar minDistSqr = radiusSqr;
- return mapPtr;
+ label nContain = 0;
+ forAll(procBb_, proci)
+ {
+ const List<treeBoundBox>& bbs = procBb_[proci];
+
+ forAll(bbs, bbi)
+ {
+ if (bbs[bbi].contains(centre))
+ {
+ // We found a bb that contains the centre. There must be
+ // a triangle inside (since otherwise the bb would never
+ // have been created).
+ if (!procContains[proci])
+ {
+ procContains[proci] = true;
+ nContain++;
+ }
+ // Minimum search distance to find the triangle
+ point near, far;
+ bbs[bbi].calcExtremities(centre, near, far);
+ minDistSqr = min(minDistSqr, magSqr(centre-far));
+ }
+ }
+ }
+
+ if (nContain > 0)
+ {
+ return Tuple2<label, scalar>(nContain, minDistSqr);
+ }
+ else
+ {
+ scalar maxDistSqr = radiusSqr;
+
+ // Pass 1: find box with nearest minimum distance. Store its maximum
+ // extent as well. Since box will always contain a triangle
+ // this guarantees at least one hit.
+ forAll(procBb_, proci)
+ {
+ const List<treeBoundBox>& bbs = procBb_[proci];
+
+ forAll(bbs, bbi)
+ {
+ if (bbs[bbi].overlaps(centre, radiusSqr))
+ {
+ point near, far;
+ bbs[bbi].calcExtremities(centre, near, far);
+
+ scalar d2 = magSqr(centre-near);
+ if (d2 < minDistSqr)
+ {
+ minDistSqr = d2;
+ maxDistSqr = min(radiusSqr, magSqr(centre-far));
+ minProci = proci;
+ }
+ }
+ }
+ }
+
+ label nOverlap = 0;
+ if (minProci >= 0)
+ {
+ // Pass 2. Find all bb in range minDistSqr..maxDistSqr
+
+ procOverlaps[minProci] = true;
+ nOverlap++;
+
+ forAll(procBb_, proci)
+ {
+ if (proci != minProci)
+ {
+ const List<treeBoundBox>& bbs = procBb_[proci];
+ forAll(bbs, bbi)
+ {
+ if (bbs[bbi].overlaps(centre, maxDistSqr))
+ {
+ procOverlaps[proci] = true;
+ nOverlap++;
+ break;
+ }
+ }
+ }
+ }
+ }
+ return Tuple2<label, scalar>(nOverlap, maxDistSqr);
+ }
}
@@ -686,6 +986,7 @@ Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
+ const bool includeLocalProcessor,
const pointField& centres,
const scalarField& radiusSqr,
@@ -708,6 +1009,17 @@ Foam::distributedTriSurfaceMesh::calcLocalQueries
// Per processor indices into allSegments to send
List<DynamicList<label>> dynSendMap(Pstream::nProcs());
+ // Pre-size
+ forAll(dynSendMap, proci)
+ {
+ dynSendMap[proci].reserve
+ (
+ (proci == Pstream::myProcNo())
+ ? centres.size()
+ : centres.size()/Pstream::nProcs()
+ );
+ }
+
// Work array - whether processor bb overlaps the bounding sphere.
boolList procBbOverlaps(Pstream::nProcs());
@@ -723,7 +1035,14 @@ Foam::distributedTriSurfaceMesh::calcLocalQueries
forAll(procBbOverlaps, proci)
{
- if (proci != Pstream::myProcNo() && procBbOverlaps[proci])
+ if
+ (
+ procBbOverlaps[proci]
+ && (
+ includeLocalProcessor
+ || proci != Pstream::myProcNo()
+ )
+ )
{
dynSendMap[proci].append(dynAllCentres.size());
dynAllSegmentMap.append(centrei);
@@ -746,229 +1065,1095 @@ Foam::distributedTriSurfaceMesh::calcLocalQueries
allSegmentMap.transfer(dynAllSegmentMap);
}
+ return autoPtr<mapDistribute>::New(std::move(sendMap));
+}
+
+
+Foam::volumeType Foam::distributedTriSurfaceMesh::edgeSide
+(
+ const point& sample,
+ const point& nearestPoint,
+ const label face0,
+ const label face1
+) const
+{
+ const triSurface& surf = static_cast<const triSurface&>(*this);
+ const pointField& points = surf.points();
+
+ // Compare to bisector. This is actually correct since edge is
+ // nearest so there is a knife-edge.
+
+ //const vectorField& faceNormals = surf.faceNormals();
+ //vector n = faceNormals[face0] + faceNormals[face1];
+ vector n = surf[face0].unitNormal(points)+surf[face1].unitNormal(points);
- // Send over how many i need to receive.
- labelListList sendSizes(Pstream::nProcs());
- sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
- forAll(sendMap, proci)
+ if (((sample - nearestPoint) & n) > 0)
+ {
+ return volumeType::OUTSIDE;
+ }
+ else
{
- sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
+ return volumeType::INSIDE;
}
- Pstream::gatherList(sendSizes);
- Pstream::scatterList(sendSizes);
+}
- // Determine order of receiving
- labelListList constructMap(Pstream::nProcs());
+Foam::label Foam::distributedTriSurfaceMesh::findOtherFace
+(
+ const labelListList& pointFaces,
+ const label nearFacei,
+ const label nearLabel
+) const
+{
+ const triSurface& surf = static_cast<const triSurface&>(*this);
+ const triSurface::FaceType& nearF = surf[nearFacei];
- // My local segments first
- constructMap[Pstream::myProcNo()] = identity
- (
- sendMap[Pstream::myProcNo()].size()
- );
+ const edge e(nearF[nearLabel], nearF[nearF.fcIndex(nearLabel)]);
- label segmenti = constructMap[Pstream::myProcNo()].size();
- forAll(constructMap, proci)
+ const labelList& pFaces = pointFaces[e[0]];
+ forAll(pFaces, pFacei)
{
- if (proci != Pstream::myProcNo())
+ const label facei = pFaces[pFacei];
+ if (facei != nearFacei)
{
- // What i need to receive is what other processor is sending to me.
- label nRecv = sendSizes[proci][Pstream::myProcNo()];
- constructMap[proci].setSize(nRecv);
+ const triSurface::FaceType& f = surf[facei];
- for (label i = 0; i < nRecv; i++)
+ int dir = f.edgeDirection(e);
+ if (dir != 0)
{
- constructMap[proci][i] = segmenti++;
+ return facei;
}
}
}
-
- return autoPtr<mapDistribute>::New
- (
- segmenti, // size after construction
- std::move(sendMap),
- std::move(constructMap)
- );
+ return -1;
}
-// Find bounding boxes that guarantee a more or less uniform distribution
-// of triangles. Decomposition in here is only used to get the bounding
-// boxes, actual decomposition is done later on.
-// Returns a per processor a list of bounding boxes that most accurately
-// describe the shape. For now just a single bounding box per processor but
-// optimisation might be to determine a better fitting shape.
-Foam::List<Foam::List<Foam::treeBoundBox>>
-Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
+void Foam::distributedTriSurfaceMesh::calcFaceFaces
(
- const triSurface& s
+ const triSurface& s,
+ const labelListList& pointFaces,
+ labelListList& faceFaces
)
{
- if (!decomposer_.valid())
+ faceFaces.setSize(s.size());
+
+ DynamicList<label> nbrs;
+
+ forAll(faceFaces, facei)
{
- // Use singleton decomposeParDict. Cannot use decompositionModel
- // here since we've only got Time and not a mesh.
+ const labelledTri& f = s[facei];
- const auto* dictPtr =
- searchableSurface::time().findObject<IOdictionary>
- (
- // == decompositionModel::canonicalName
- "decomposeParDict"
- );
+ nbrs.reserve(f.size());
+ nbrs.clear();
- if (dictPtr)
- {
- decomposer_ = decompositionMethod::New(*dictPtr);
- }
- else
+ forAll(f, fp)
{
- if (!decomposeParDict_.valid())
+ const edge e(f[fp], f[f.fcIndex(fp)]);
+ const labelList& pFaces = pointFaces[f[fp]];
+ forAll(pFaces, pFacei)
{
- decomposeParDict_.reset
- (
- new IOdictionary
- (
- IOobject
- (
- // == decompositionModel::canonicalName
- "decomposeParDict",
- searchableSurface::time().system(),
- searchableSurface::time(),
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- )
- );
+ const label otherFacei = pFaces[pFacei];
+ if (otherFacei != facei)
+ {
+ if (s[otherFacei].edgeDirection(e) != 0)
+ {
+ if (!nbrs.find(otherFacei))
+ {
+ nbrs.append(otherFacei);
+ }
+ }
+ }
}
- decomposer_ = decompositionMethod::New(decomposeParDict_());
- }
-
-
- if (!decomposer_().parallelAware())
- {
- FatalErrorInFunction
- << "The decomposition method " << decomposer_().typeName
- << " does not decompose in parallel."
- << " Please choose one that does." << exit(FatalError);
- }
-
- if (!isA<geomDecomp>(decomposer_()))
- {
- FatalErrorInFunction
- << "The decomposition method " << decomposer_().typeName
- << " is not a geometric decomposition method." << endl
- << "Only geometric decomposition methods are currently"
- << " supported."
- << exit(FatalError);
}
+ faceFaces[facei] = std::move(nbrs);
}
+}
- // Do decomposition according to triangle centre
- pointField triCentres(s.size());
- forAll(s, trii)
+
+void Foam::distributedTriSurfaceMesh::surfaceSide
+(
+ const pointField& samples,
+ const List<pointIndexHit>& nearestInfo,
+ List<volumeType>& volType
+) const
+{
+ if (debug)
{
- triCentres[trii] = s[trii].centre(s.points());
+ Pout<< "distributedTriSurfaceMesh::surfaceSide :"
+ << " finding surface side given points on surface for "
+ << samples.size() << " samples" << endl;
}
+ // Use global index to send local tri and nearest back to originating
+ // processor
- geomDecomp& decomposer = refCast<geomDecomp>(decomposer_());
-
- // Do the actual decomposition
- labelList distribution(decomposer.decompose(triCentres));
+ labelList triangleIndex(nearestInfo.size());
+ autoPtr<mapDistribute> mapPtr
+ (
+ calcLocalQueries
+ (
+ nearestInfo,
+ triangleIndex
+ )
+ );
+ const mapDistribute& map = mapPtr();
- // Find bounding box for all triangles on new distribution.
+ // Send over samples
+ pointField localSamples(samples);
+ map.distribute(localSamples);
- // Initialise to inverted box
- List<List<treeBoundBox>> bbs(Pstream::nProcs());
- forAll(bbs, proci)
- {
- bbs[proci].setSize(1, treeBoundBox(boundBox::invertedBox));
- }
- forAll(s, trii)
- {
- const triSurface::FaceType& f = s[trii];
+ // Do my tests
+ // ~~~~~~~~~~~
- treeBoundBox& bb = bbs[distribution[trii]][0];
- bb.add(s.points(), f);
- }
+ volType.setSize(triangleIndex.size());
+ volType = volumeType::UNKNOWN;
- // Now combine for all processors and convert to correct format.
- forAll(bbs, proci)
+ const triSurface& surf = static_cast<const triSurface&>(*this);
+ const pointField& points = surf.points();
{
- Pstream::listCombineGather(bbs[proci], plusEqOp<boundBox>());
- Pstream::listCombineScatter(bbs[proci]);
- }
-
- return bbs;
-}
+ //const labelListList& pointFaces = surf.pointFaces();
+ // Construct pointFaces. Let's hope surface has compact point
+ // numbering ...
+ labelListList pointFaces;
+ invertManyToMany(points.size(), surf, pointFaces);
+ EdgeMap<labelPair> edgeToFaces;
-// Does any part of triangle overlap bb.
-bool Foam::distributedTriSurfaceMesh::overlaps
-(
- const List<treeBoundBox>& bbs,
- const point& p0,
- const point& p1,
- const point& p2
-)
-{
- treeBoundBox triBb(p0);
- triBb.add(p1);
- triBb.add(p2);
+ forAll(triangleIndex, i)
+ {
+ label facei = triangleIndex[i];
+ const triSurface::FaceType& f = surf[facei];
+ const point& sample = localSamples[i];
- forAll(bbs, bbi)
- {
- const treeBoundBox& bb = bbs[bbi];
+ // Find where point is on face
+ label nearType, nearLabel;
+ pointHit pHit =
+ f.nearestPointClassify(sample, points, nearType, nearLabel);
- // Exact test of triangle intersecting bb
+ const point& nearestPoint(pHit.rawPoint());
- // Quick rejection. If whole bounding box of tri is outside cubeBb then
- // there will be no intersection.
- if (bb.overlaps(triBb))
- {
- // Check if one or more triangle point inside
- if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
+ if (nearType == triPointRef::NONE)
{
- // One or more points inside
- return true;
- }
+ const vector sampleNearestVec = (sample - nearestPoint);
- // Now we have the difficult case: all points are outside but
- // connecting edges might go through cube. Use fast intersection
- // of bounding box.
- bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
+ // Nearest to face interior. Use faceNormal to determine side
+ //scalar c = sampleNearestVec & surf.faceNormals()[facei];
+ scalar c = sampleNearestVec & surf[facei].normal(points);
- if (intersect)
+ if (c > 0)
+ {
+ volType[i] = volumeType::OUTSIDE;
+ }
+ else
+ {
+ volType[i] = volumeType::INSIDE;
+ }
+ }
+ else if (nearType == triPointRef::EDGE)
{
- return true;
+ // Nearest to edge nearLabel. Note that this can only be a
+ // knife-edge
+ // situation since otherwise the nearest point could never be
+ // the edge.
+
+ label otherFacei = findOtherFace(pointFaces, facei, nearLabel);
+ if (otherFacei != -1)
+ {
+ volType[i] =
+ edgeSide(sample, nearestPoint, facei, otherFacei);
+ }
+ else
+ {
+ // Open edge. Leave volume type unknown
+ }
}
- }
- }
- return false;
-}
+ else
+ {
+ // Nearest to point. Could use pointNormal here but is not
+ // correct.
+ // Instead determine which edge using point is nearest and
+ // use test above (nearType == triPointRef::EDGE).
+ const label pointi = f[nearLabel];
+ const labelList& pFaces = pointFaces[pointi];
+ const vector sampleNearestVec = (sample - nearestPoint);
-void Foam::distributedTriSurfaceMesh::subsetMeshMap
-(
- const triSurface& s,
- const boolList& include,
- const label nIncluded,
- labelList& newToOldPoints,
- labelList& oldToNewPoints,
- labelList& newToOldFaces
-)
-{
- newToOldFaces.setSize(nIncluded);
- newToOldPoints.setSize(s.points().size());
- oldToNewPoints.setSize(s.points().size());
- oldToNewPoints = -1;
- {
- label facei = 0;
- label pointi = 0;
+ // Loop over all faces. Check if have both edge faces. Do
+ // test.
+ edgeToFaces.clear();
- forAll(include, oldFacei)
- {
+ scalar maxCosAngle = -GREAT;
+ labelPair maxEdgeFaces(-1, -1);
+
+ forAll(pFaces, pFacei)
+ {
+ label facei = pFaces[pFacei];
+ const triSurface::FaceType& f = surf[facei];
+
+ label fp = findIndex(f, pointi);
+ label p1 = f[f.fcIndex(fp)];
+ label pMin1 = f[f.rcIndex(fp)];
+
+ Pair<edge> edges
+ (
+ edge(pointi, p1),
+ edge(pointi, pMin1)
+ );
+
+ // Check edge fp-to-fp+1 and fp-1
+ // determine distance/angle to nearPoint
+ for (const edge& e : edges)
+ {
+ auto iter = edgeToFaces.find(e);
+ if (iter.found())
+ {
+ if (iter().second() == -1)
+ {
+ // Found second face. Now we have edge+faces
+ iter().second() = facei;
+
+ vector eVec(e.vec(points));
+ scalar magEVec = mag(eVec);
+
+ if (magEVec > VSMALL)
+ {
+ eVec /= magEVec;
+
+ // Determine edge most in direction of
+ // sample
+ scalar cosAngle(sampleNearestVec&eVec);
+ if (cosAngle > maxCosAngle)
+ {
+ maxCosAngle = cosAngle;
+ maxEdgeFaces = iter();
+ }
+ }
+ }
+ else
+ {
+ FatalErrorInFunction << "Not closed"
+ << exit(FatalError);
+ }
+ }
+ else
+ {
+ edgeToFaces.insert(e, labelPair(facei, -1));
+ }
+ }
+ }
+
+
+ // Check that surface is closed
+ bool closed = true;
+ for (auto iter : edgeToFaces)
+ {
+ if (iter[0] == -1 || iter[1] == -1)
+ {
+ closed = false;
+ break;
+ }
+ }
+ if (closed)
+ {
+ volType[i] = edgeSide
+ (
+ sample,
+ nearestPoint,
+ maxEdgeFaces[0],
+ maxEdgeFaces[1]
+ );
+ }
+ }
+ }
+ }
+
+
+ // Send back results
+ // ~~~~~~~~~~~~~~~~~
+
+ // Note that we make sure that if multiple processors hold data only
+ // the one with inside/outside wins. (though this should already be
+ // handled by the fact we have a unique nearest triangle so we only
+ // send the volume-query to a single processor)
+
+
+ //forAll(localSamples, i)
+ //{
+ // Pout<< "surfaceSide : for localSample:" << localSamples[i]
+ // << " found volType:" << volumeType::names[volType[i]]
+ // << endl;
+ //}
+
+ const volumeType zero(volumeType::UNKNOWN);
+ mapDistributeBase::distribute
+ (
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(0),
+ nearestInfo.size(),
+ map.constructMap(),
+ map.constructHasFlip(),
+ map.subMap(),
+ map.subHasFlip(),
+ volType,
+ volumeCombineOp(),
+ noOp(), // no flipping
+ zero
+ );
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::surfaceSide :"
+ << " finished finding surface side given points on surface for "
+ << samples.size() << " samples" << endl;
+ }
+}
+
+
+void Foam::distributedTriSurfaceMesh::collectLeafMids
+(
+ const label nodeI,
+ DynamicField<point>& midPoints
+) const
+{
+ const indexedOctree<treeDataTriSurface>::node& nod = tree().nodes()[nodeI];
+
+ for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
+ {
+ const labelBits index = nod.subNodes_[octant];
+
+ if (indexedOctree<treeDataTriSurface>::isNode(index))
+ {
+ // tree node. Recurse.
+ collectLeafMids
+ (
+ indexedOctree<treeDataTriSurface>::getNode(index),
+ midPoints
+ );
+ }
+ else if (indexedOctree<treeDataTriSurface>::isContent(index))
+ {}
+ else
+ {
+ // No data in this octant. Set type for octant acc. to the mid
+ // of its bounding box.
+ const treeBoundBox subBb = nod.bb_.subBbox(octant);
+ midPoints.append(subBb.midpoint());
+ }
+ }
+}
+
+
+Foam::volumeType Foam::distributedTriSurfaceMesh::calcVolumeType
+(
+ const List<volumeType>& midPointTypes,
+ label& midPointi,
+ PackedList<2>& nodeTypes,
+ const label nodeI
+) const
+{
+ // Pre-calculates wherever possible the volume status per node/subnode.
+ // Recurses to determine status of lowest level boxes. Level above is
+ // combination of octants below.
+
+ const indexedOctree<treeDataTriSurface>::node& nod = tree().nodes()[nodeI];
+
+ volumeType myType = volumeType::UNKNOWN;
+
+ for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
+ {
+ volumeType subType;
+
+ const labelBits index = nod.subNodes_[octant];
+
+ if (indexedOctree<treeDataTriSurface>::isNode(index))
+ {
+ // tree node. Recurse.
+ subType = calcVolumeType
+ (
+ midPointTypes,
+ midPointi,
+ nodeTypes,
+ indexedOctree<treeDataTriSurface>::getNode(index)
+ );
+ }
+ else if (indexedOctree<treeDataTriSurface>::isContent(index))
+ {
+ // Contents. Depending on position in box might be on either
+ // side.
+ subType = volumeType::MIXED;
+ }
+ else
+ {
+ // No data in this octant. Set type for octant acc. to the mid
+ // of its bounding box.
+ //Pout<< " for leaf at bb:" << nod.bb_.subBbox(octant)
+ // << " node:" << nodeI
+ // << " octant:" << octant
+ // << " caching volType to:" << midPointTypes[midPointi] << endl;
+ subType = midPointTypes[midPointi++];
+ }
+
+ // Store octant type
+ nodeTypes.set((nodeI<<3)+octant, subType);
+
+ // Combine sub node types into type for treeNode. Result is 'mixed' if
+ // types differ among subnodes.
+ if (myType == volumeType::UNKNOWN)
+ {
+ myType = subType;
+ }
+ else if (subType != myType)
+ {
+ myType = volumeType::MIXED;
+ }
+ }
+ return myType;
+}
+
+
+Foam::volumeType Foam::distributedTriSurfaceMesh::cachedVolumeType
+(
+ const label nodeI,
+ const point& sample
+) const
+{
+ const indexedOctree<treeDataTriSurface>::node& nod = tree().nodes()[nodeI];
+
+ direction octant = nod.bb_.subOctant(sample);
+
+ volumeType octantType = volumeType::type
+ (
+ tree().nodeTypes().get((nodeI<<3)+octant)
+ );
+
+ if (octantType == volumeType::INSIDE)
+ {
+ return octantType;
+ }
+ else if (octantType == volumeType::OUTSIDE)
+ {
+ return octantType;
+ }
+ else if (octantType == volumeType::UNKNOWN)
+ {
+ // Can happen for e.g. non-manifold surfaces.
+ return octantType;
+ }
+ else if (octantType == volumeType::MIXED)
+ {
+ labelBits index = nod.subNodes_[octant];
+
+ if (indexedOctree<treeDataTriSurface>::isNode(index))
+ {
+ // Recurse
+ volumeType subType = cachedVolumeType
+ (
+ indexedOctree<treeDataTriSurface>::getNode(index),
+ sample
+ );
+
+ return subType;
+ }
+ else if (indexedOctree<treeDataTriSurface>::isContent(index))
+ {
+ // Content.
+ return volumeType::UNKNOWN;
+ }
+ else
+ {
+ // Empty node. Cannot have 'mixed' as its type since not divided
+ // up and has no items inside it.
+ FatalErrorInFunction
+ << "Sample:" << sample << " node:" << nodeI
+ << " with bb:" << nod.bb_ << nl
+ << "Empty subnode has invalid volume type MIXED."
+ << abort(FatalError);
+
+ return volumeType::UNKNOWN;
+ }
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Sample:" << sample << " at node:" << nodeI
+ << " octant:" << octant
+ << " with bb:" << nod.bb_.subBbox(octant) << nl
+ << "Node has invalid volume type " << octantType
+ << abort(FatalError);
+
+ return volumeType::UNKNOWN;
+ }
+}
+
+
+// Find bounding boxes that guarantee a more or less uniform distribution
+// of triangles. Decomposition in here is only used to get the bounding
+// boxes, actual decomposition is done later on.
+// Returns a per processor a list of bounding boxes that most accurately
+// describe the shape. For now just a single bounding box per processor but
+// optimisation might be to determine a better fitting shape.
+Foam::List<Foam::List<Foam::treeBoundBox>>
+Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
+(
+ const triSurface& s
+)
+{
+ addProfiling
+ (
+ distribute,
+ "distributedTriSurfaceMesh::independentlyDistributedBbs"
+ );
+
+ if (!decomposer_.valid())
+ {
+ // Use singleton decomposeParDict. Cannot use decompositionModel
+ // here since we've only got Time and not a mesh.
+
+ const auto* dictPtr =
+ searchableSurface::time().findObject<IOdictionary>
+ (
+ // == decompositionModel::canonicalName
+ "decomposeParDict"
+ );
+
+ if (dictPtr)
+ {
+ decomposer_ = decompositionMethod::New(*dictPtr);
+ }
+ else
+ {
+ if (!decomposeParDict_.valid())
+ {
+ decomposeParDict_.reset
+ (
+ new IOdictionary
+ (
+ IOobject
+ (
+ // == decompositionModel::canonicalName
+ "decomposeParDict",
+ searchableSurface::time().system(),
+ searchableSurface::time(),
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ )
+ );
+ }
+ decomposer_ = decompositionMethod::New(decomposeParDict_());
+ }
+ }
+
+
+ // Initialise to inverted box
+ List<List<treeBoundBox>> bbs(Pstream::nProcs());
+ forAll(bbs, proci)
+ {
+ bbs[proci].setSize(1, treeBoundBox(boundBox::invertedBox));
+ }
+
+
+ const globalIndex& triIndexer = globalTris();
+
+ bool masterOnly;
+ {
+ masterOnly = true;
+ for (label proci = 1; proci < Pstream::nProcs(); proci++)
+ {
+ if (triIndexer.localSize(proci) != 0)
+ {
+ masterOnly = false;
+ break;
+ }
+ }
+ }
+
+ if (masterOnly)
+ {
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " determining master-only decomposition for " << s.size()
+ << " centroids for " << searchableSurface::name() << endl;
+ }
+
+ // Triangle centres
+ pointField triCentres(s.size());
+ forAll(s, trii)
+ {
+ triCentres[trii] = s[trii].centre(s.points());
+ }
+
+ labelList distribution;
+ if (!isA<geomDecomp>(decomposer_()))
+ {
+ // Use connectivity
+ labelListList pointFaces;
+ invertManyToMany(s.points().size(), s, pointFaces);
+ labelListList faceFaces(s.size());
+ calcFaceFaces(s, pointFaces, faceFaces);
+
+ // Do the actual decomposition
+ const bool oldParRun = UPstream::parRun();
+ UPstream::parRun() = false;
+ distribution = decomposer_().decompose(faceFaces, triCentres);
+ UPstream::parRun() = oldParRun;
+ }
+ else
+ {
+ // Do the actual decomposition
+ distribution = decomposer_().decompose(triCentres);
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " determining processor bounding boxes" << endl;
+ }
+
+ // Find bounding box for all triangles on new distribution.
+ forAll(s, trii)
+ {
+ const triSurface::FaceType& f = s[trii];
+
+ treeBoundBox& bb = bbs[distribution[trii]][0];
+ bb.add(s.points(), f);
+ }
+
+ // Now combine for all processors and convert to correct format.
+ forAll(bbs, proci)
+ {
+ Pstream::listCombineGather(bbs[proci], plusEqOp<boundBox>());
+ Pstream::listCombineScatter(bbs[proci]);
+ }
+ }
+ else if (distType_ == DISTRIBUTED)
+ {
+ // Fully distributed decomposition. Does not filter duplicate
+ // triangles.
+ if (!decomposer_().parallelAware())
+ {
+ FatalErrorInFunction
+ << "The decomposition method " << decomposer_().typeName
+ << " does not decompose in parallel."
+ << " Please choose one that does." << exit(FatalError);
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " determining decomposition for " << s.size()
+ << " centroids" << endl;
+ }
+
+ // Triangle centres
+ pointField triCentres(s.size());
+ forAll(s, trii)
+ {
+ triCentres[trii] = s[trii].centre(s.points());
+ }
+
+ labelList distribution = decomposer_().decompose(triCentres);
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " determining processor bounding boxes for "
+ << searchableSurface::name() << endl;
+ }
+
+ // Find bounding box for all triangles on new distribution.
+ forAll(s, trii)
+ {
+ const triSurface::FaceType& f = s[trii];
+
+ treeBoundBox& bb = bbs[distribution[trii]][0];
+ bb.add(s.points(), f);
+ }
+
+ // Now combine for all processors and convert to correct format.
+ forAll(bbs, proci)
+ {
+ Pstream::listCombineGather(bbs[proci], plusEqOp<boundBox>());
+ Pstream::listCombineScatter(bbs[proci]);
+ }
+ }
+// //// Tbd. initial duplicate filtering of border points only
+// if (distType_ == DISTRIBUTED)
+// {
+// // Fully distributed decomposition. Does not filter duplicate
+// // triangles.
+// if (!decomposer_().parallelAware())
+// {
+// FatalErrorInFunction
+// << "The decomposition method " << decomposer_().typeName
+// << " does not decompose in parallel."
+// << " Please choose one that does." << exit(FatalError);
+// }
+//
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+// << " determining decomposition for " << s.size()
+// << " centroids" << endl;
+// }
+// const pointField& points = s.points();
+//
+// pointField triCentres(s.size());
+// forAll(s, trii)
+// {
+// triCentres[trii] = s[trii].centre(points);
+// }
+//
+// // Collect all triangles not fully inside the current bb
+// DynamicList<label> borderTris(s.size()/Pstream::nProcs());
+//
+// const List<treeBoundBox>& myBbs = procBb_[Pstream::myProcNo];
+//
+// boolList includedFace;
+// overlappingSurface(s, myBbs, includedFace);
+// boolList internalOrBorderFace(includedFace);
+// forAll(s, trii)
+// {
+// if (includedFace[trii])
+// {
+// // Triangle is either inside or part-inside. Exclude fully
+// // inside triangles.
+// const labelledTri& f = s[trii];
+// const point& p0 = points[f[0]];
+// const point& p1 = points[f[1]];
+// const point& p2 = points[f[2]];
+// if
+// (
+// !contains(myBbs, p0)
+// || !contains(myBbs, p1)
+// || !contains(myBbs, p2)
+// )
+// {
+// borderTris.append(trii);
+// }
+// }
+// }
+//
+// const label nBorderTris = borderTris.size();
+//
+// Pout<< "Have " << borderTris.size() << " border triangles out of "
+// << s.size() << endl;
+//
+// labelListList sendMap(Pstream::nProcs());
+// sendMap[0] = std::move(borderTris);
+//
+// const mapDistribute map(std::move(sendMap));
+//
+// // Gather all borderTris
+// //globalIndex globalBorderTris(borderTris.size());
+// //pointField globalBorderCentres(allCentres, borderTris);
+// //globalBorderTris.gather
+// //(
+// // UPstream::worldComm,
+// // UPstream::procID(Pstream::worldComm),
+// // globalBorderCentres
+// //);
+// pointField globalBorderCentres(allCentres);
+// map.distribute(globalBorderCentres);
+//
+// // Merge on master
+// labelList masterBorder(borderTris.size(), -1);
+// if (Pstream::master())
+// {
+// labelList allToMerged;
+// label nMerged = mergePoints
+// (
+// globalBorderCentres,
+// mergeDist_,
+// false, //const bool verbose,
+// allToMerged
+// // maybe bounds().mid() ?
+// );
+//
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::"
+// << "independentlyDistributedBbs :"
+// << " merged " << globalBorderCentres.size()
+// << " border centroids down to " << nMerged << endl;
+// }
+//
+// labelList mergedMaster(nMerged, -1);
+// isMaster.setSize(globalBorderCentres.size(), false);
+// forAll(allToMerged, i)
+// {
+// label mergedi = allToMerged[i];
+// if (mergedMaster[mergedi] == -1)
+// {
+// mergedMaster[mergedi] = i;
+// isMaster[i] = true;
+// }
+// }
+// forAll(allToMerged, i)
+// {
+// label mergedi = allToMerged[i];
+// masterBorder[i] = mergedMaster[mergedi];
+// }
+// }
+// //globalBorderTris.scatter
+// //(
+// // UPstream::worldComm,
+// // UPstream::procID(Pstream::worldComm),
+// // isMasterPoint
+// //);
+// //boolList isMasterBorder(s.size(), false);
+// //forAll(borderTris, i)
+// //{
+// // isMasterBorder[borderTris[i]] = isMasterPoint[i];
+// //}
+// map.reverseDistribute(s.size(), isMaster);
+//
+// // Collect all non-border or master-border points
+// DynamicList<label> triMap(s.size());
+// forAll(s, trii)
+// {
+// if (includedFace[trii])
+// {
+// // Triangle is either inside or part-inside. Exclude fully
+// // inside triangles.
+// const labelledTri& f = s[trii];
+// const point& p0 = points[f[0]];
+// const point& p1 = points[f[1]];
+// const point& p2 = points[f[2]];
+//
+// if
+// (
+// contains(myBbs, p0)
+// && contains(myBbs, p1)
+// && contains(myBbs, p2)
+// )
+// {
+// // Internal
+// triMap.append(trii);
+// }
+// else if (isMasterBorder[trii])
+// {
+// // Part overlapping and master of overlap
+// triMap.append(trii);
+// }
+// }
+// }
+//
+// pointField masterCentres(allCentres, triMap);
+//
+// // Do the actual decomposition
+// labelList masterDistribution(decomposer_().decompose(masterCentres));
+//
+// // Make map to get the decomposition from the master of each border
+// labelList borderGlobalMaster(borderTris.size());
+// forAll(borderTris, borderi)
+// {
+// borderGlobalMaster[borderi] = ..masterTri
+// }
+// mapDistribute map(globalBorderTris, borderGlobalMaster
+//
+// // Send decomposition
+//
+//
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+// << " determining processor bounding boxes" << endl;
+// }
+//
+// // Find bounding box for all triangles on new distribution.
+// forAll(s, trii)
+// {
+// const triSurface::FaceType& f = s[trii];
+//
+// treeBoundBox& bb = bbs[distribution[trii]][0];
+// bb.add(s.points(), f);
+// }
+//
+// // Now combine for all processors and convert to correct format.
+// forAll(bbs, proci)
+// {
+// Pstream::listCombineGather(bbs[proci], plusEqOp<boundBox>());
+// Pstream::listCombineScatter(bbs[proci]);
+// }
+// }
+ else
+ {
+ // Master-only decomposition. Filters duplicate triangles so repeatable.
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " collecting all centroids" << endl;
+ }
+
+ // Collect all triangle centres
+ pointField allCentres(s.size());
+ {
+ forAll(s, trii)
+ {
+ allCentres[trii] = s[trii].centre(s.points());
+ }
+ globalTris().gather
+ (
+ UPstream::worldComm,
+ UPstream::procID(Pstream::worldComm),
+ allCentres
+ );
+ }
+
+ // Determine local decomposition
+ labelList distribution(s.size());
+ {
+ labelList allDistribution;
+ if (Pstream::master())
+ {
+ labelList allToMerged;
+ label nMerged = mergePoints
+ (
+ allCentres,
+ mergeDist_,
+ false, //const bool verbose,
+ allToMerged
+ // maybe bounds().mid() ?
+ );
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::"
+ << "independentlyDistributedBbs :"
+ << " merged " << allCentres.size()
+ << " centroids down to " << nMerged << endl;
+ }
+
+ // Collect merged points
+ pointField mergedPoints(nMerged);
+ UIndirectList<point>(mergedPoints, allToMerged) = allCentres;
+
+ // Decompose merged centres
+ const bool oldParRun = UPstream::parRun();
+ UPstream::parRun() = false;
+ labelList mergedDist(decomposer_().decompose(mergedPoints));
+ UPstream::parRun() = oldParRun;
+
+ // Convert to all
+ allDistribution = UIndirectList<label>
+ (
+ mergedDist,
+ allToMerged
+ );
+ }
+
+ // Scatter back to processors
+ globalTris().scatter
+ (
+ UPstream::worldComm,
+ UPstream::procID(Pstream::worldComm),
+ allDistribution,
+ distribution
+ );
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::"
+ << "independentlyDistributedBbs :"
+ << " determined decomposition" << endl;
+ }
+ }
+
+ // Find bounding box for all triangles on new distribution.
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
+ << " determining processor bounding boxes for "
+ << searchableSurface::name() << endl;
+ }
+
+ forAll(s, trii)
+ {
+ const triSurface::FaceType& f = s[trii];
+ treeBoundBox& bb = bbs[distribution[trii]][0];
+ bb.add(s.points(), f);
+ }
+
+ // Now combine for all processors and convert to correct format.
+ forAll(bbs, proci)
+ {
+ Pstream::listCombineGather(bbs[proci], plusEqOp<boundBox>());
+ Pstream::listCombineScatter(bbs[proci]);
+ }
+ }
+ return bbs;
+}
+
+
+// Does any part of triangle overlap bb.
+bool Foam::distributedTriSurfaceMesh::overlaps
+(
+ const List<treeBoundBox>& bbs,
+ const point& p0,
+ const point& p1,
+ const point& p2
+)
+{
+ treeBoundBox triBb(p0);
+ triBb.add(p1);
+ triBb.add(p2);
+
+ forAll(bbs, bbi)
+ {
+ const treeBoundBox& bb = bbs[bbi];
+
+ // Exact test of triangle intersecting bb
+
+ // Quick rejection. If whole bounding box of tri is outside cubeBb then
+ // there will be no intersection.
+ if (bb.overlaps(triBb))
+ {
+ // Check if one or more triangle point inside
+ if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
+ {
+ // One or more points inside
+ return true;
+ }
+
+ // Now we have the difficult case: all points are outside but
+ // connecting edges might go through cube. Use fast intersection
+ // of bounding box.
+ bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);
+
+ if (intersect)
+ {
+ return true;
+ }
+ }
+ }
+ return false;
+}
+
+
+void Foam::distributedTriSurfaceMesh::subsetMeshMap
+(
+ const triSurface& s,
+ const boolList& include,
+ const label nIncluded,
+ labelList& newToOldPoints,
+ labelList& oldToNewPoints,
+ labelList& newToOldFaces
+)
+{
+ newToOldFaces.setSize(nIncluded);
+ newToOldPoints.setSize(s.points().size());
+ oldToNewPoints.setSize(s.points().size());
+ oldToNewPoints = -1;
+ {
+ label facei = 0;
+ label pointi = 0;
+
+ forAll(include, oldFacei)
+ {
if (include[oldFacei])
{
// Store new faces compact
@@ -1297,8 +2482,10 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
searchableSurface::registerObject()
),
dict
- )
+ ),
+ currentDistType_(FROZEN) // only used to trigger re-distribution
{
+ // Read from the provided dictionary
read();
bounds().reduce();
@@ -1331,14 +2518,14 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
IOobject
(
io.name(),
- io.time().findInstance(io.local(), word::null),
+ findLocalInstance(io), // findInstance with parent searching
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
io.registerObject()
),
- false
+ triSurfaceMesh::masterOnly // allow parent searching
),
dict_
(
@@ -1348,34 +2535,84 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
searchableSurface::instance(),
searchableSurface::local(),
searchableSurface::db(),
- searchableSurface::readOpt(),
+ (
+ (
+ searchableSurface::readOpt()
+ == IOobject::MUST_READ
+ || searchableSurface::readOpt()
+ == IOobject::MUST_READ_IF_MODIFIED
+ )
+ ? IOobject::READ_IF_PRESENT
+ : searchableSurface::readOpt()
+ ),
searchableSurface::writeOpt(),
searchableSurface::registerObject()
- )
- )
+ ),
+ dictionary()
+ ),
+ currentDistType_(FROZEN) // only used to trigger re-distribution
{
+ // Read from the local, decomposed dictionary
read();
bounds().reduce();
- if (debug)
- {
- InfoInFunction << "Read distributedTriSurface from " << io.objectPath()
- << ':' << endl;
- writeStats(Info);
+ const fileName actualFile(checkFile(io, true));
- labelList nTris(Pstream::nProcs());
- nTris[Pstream::myProcNo()] = triSurface::size();
- Pstream::gatherList(nTris);
- Pstream::scatterList(nTris);
+ if
+ (
+ actualFile != io.localFilePath(triSurfaceMesh::typeName)
+ && (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
+ )
+ {
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name()
+ << " from parent path " << actualFile << endl;
+ }
- Info<< endl<< "\tproc\ttris\tbb" << endl;
- forAll(nTris, proci)
+ if (Pstream::parRun())
{
- Info<< '\t' << proci << '\t' << nTris[proci]
- << '\t' << procBb_[proci] << endl;
+ // Distribute (checks that distType != currentDistType_ so should
+ // always trigger re-distribution)
+ List<treeBoundBox> bbs;
+ autoPtr<mapDistribute> faceMap;
+ autoPtr<mapDistribute> pointMap;
+ distribute
+ (
+ bbs,
+ true, // keep unmapped triangles
+ faceMap,
+ pointMap
+ );
}
- Info<< endl;
+ }
+ else
+ {
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name()
+ << " from actual path " << actualFile << ':' << endl;
+
+ labelList nTris(Pstream::nProcs());
+ nTris[Pstream::myProcNo()] = triSurface::size();
+ Pstream::gatherList(nTris);
+ Pstream::scatterList(nTris);
+
+ Info<< endl<< "\tproc\ttris\tbb" << endl;
+ forAll(nTris, proci)
+ {
+ Info<< '\t' << proci << '\t' << nTris[proci]
+ << '\t' << procBb_[proci] << endl;
+ }
+ Info<< endl;
+ }
+ }
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name() << ':'
+ << endl;
+ writeStats(Info);
}
}
@@ -1391,7 +2628,7 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
IOobject
(
io.name(),
- io.time().findInstance(io.local(), word::null),
+ findLocalInstance(io),
io.local(),
io.db(),
io.readOpt(),
@@ -1399,7 +2636,7 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
io.registerObject()
),
dict,
- false
+ triSurfaceMesh::masterOnly
),
dict_
(
@@ -1409,138 +2646,793 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
searchableSurface::instance(),
searchableSurface::local(),
searchableSurface::db(),
- searchableSurface::readOpt(),
+ (
+ (
+ searchableSurface::readOpt()
+ == IOobject::MUST_READ
+ || searchableSurface::readOpt()
+ == IOobject::MUST_READ_IF_MODIFIED
+ )
+ ? IOobject::READ_IF_PRESENT
+ : searchableSurface::readOpt()
+ ),
searchableSurface::writeOpt(),
searchableSurface::registerObject()
- )
+ ),
+ dictionary()
+ ),
+ currentDistType_(FROZEN) // only used to trigger re-distribution
+{
+ // Read from the local, decomposed dictionary
+ read();
+
+ // Optionally override settings from provided dictionary
+ {
+ // Wanted distribution type
+ distType_ = distributionTypeNames_.lookupOrDefault
+ (
+ "distributionType",
+ dict_,
+ distType_
+ );
+
+ // Merge distance
+ dict_.readIfPresent("mergeDistance", mergeDist_);
+
+ // Distribution type
+ bool closed;
+ if (dict_.readIfPresent<bool>("closed", closed))
+ {
+ surfaceClosed_ = closed;
+ }
+
+ outsideVolType_ = volumeType::names.lookupOrDefault
+ (
+ "outsideVolumeType",
+ dict_,
+ outsideVolType_
+ );
+ }
+
+
+ bounds().reduce();
+
+ const fileName actualFile(checkFile(io, dict, true));
+
+ if
+ (
+ actualFile != io.localFilePath(triSurfaceMesh::typeName)
+ && (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
)
+ {
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name()
+ << " from parent path " << actualFile
+ << " and dictionary" << endl;
+ }
+
+ if (Pstream::parRun())
+ {
+ // Distribute (checks that distType != currentDistType_ so should
+ // always trigger re-distribution)
+ List<treeBoundBox> bbs;
+ autoPtr<mapDistribute> faceMap;
+ autoPtr<mapDistribute> pointMap;
+ distribute
+ (
+ bbs,
+ true, // keep unmapped triangles
+ faceMap,
+ pointMap
+ );
+ }
+ }
+ else
+ {
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name()
+ << " from actual path " << actualFile
+ << " and dictionary:" << endl;
+
+ labelList nTris(Pstream::nProcs());
+ nTris[Pstream::myProcNo()] = triSurface::size();
+ Pstream::gatherList(nTris);
+ Pstream::scatterList(nTris);
+
+ Info<< endl<< "\tproc\ttris\tbb" << endl;
+ forAll(nTris, proci)
+ {
+ Info<< '\t' << proci << '\t' << nTris[proci]
+ << '\t' << procBb_[proci] << endl;
+ }
+ Info<< endl;
+ }
+ }
+ if (debug)
+ {
+ InfoInFunction << "Read distributedTriSurface " << io.name() << ':'
+ << endl;
+ writeStats(Info);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
+{
+ clearOut();
+}
+
+
+void Foam::distributedTriSurfaceMesh::clearOut()
+{
+ globalTris_.clear();
+ triSurfaceMesh::clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
+{
+ if (!globalTris_.valid())
+ {
+ globalTris_.reset(new globalIndex(triSurface::size()));
+ }
+ return *globalTris_;
+}
+
+
+//void Foam::distributedTriSurfaceMesh::findNearest
+//(
+// const pointField& samples,
+// const scalarField& nearestDistSqr,
+// List<pointIndexHit>& info
+//) const
+//{
+// if (!Pstream::parRun())
+// {
+// triSurfaceMesh::findNearest(samples, nearestDistSqr, info);
+// return;
+// }
+//
+// addProfiling
+// (
+// findNearest,
+// "distributedTriSurfaceMesh::findNearest"
+// );
+//
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::findNearest :"
+// << " trying to find nearest for " << samples.size()
+// << " samples with max sphere "
+// << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
+// << endl;
+// }
+//
+//
+// const indexedOctree<treeDataTriSurface>& octree = tree();
+//
+// // Important:force synchronised construction of indexing
+// const globalIndex& triIndexer = globalTris();
+//
+//
+// // Initialise
+// // ~~~~~~~~~~
+//
+// info.setSize(samples.size());
+// forAll(info, i)
+// {
+// info[i].setMiss();
+// }
+//
+//
+//
+// // Do any local queries
+// // ~~~~~~~~~~~~~~~~~~~~
+//
+// label nLocal = 0;
+//
+// {
+// // Work array - whether processor bb overlaps the bounding sphere.
+// boolList procBbOverlaps(Pstream::nProcs());
+//
+// forAll(samples, i)
+// {
+// // Find the processor this sample+radius overlaps.
+// label nProcs = calcOverlappingProcs
+// (
+// samples[i],
+// nearestDistSqr[i],
+// procBbOverlaps
+// );
+//
+// // Overlaps local processor?
+// if (procBbOverlaps[Pstream::myProcNo()])
+// {
+// info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
+// if (info[i].hit())
+// {
+// if
+// (
+// surfaceClosed_
+// && !contains(procBb_[proci], info[i].hitPoint())
+// )
+// {
+// // Nearest point is not on local processor so the
+// // the triangle is only there because some other bit
+// // of it
+// // is on it. Assume there is another processor that
+// // holds
+// // the full surrounding of the triangle so we can
+// // clear this particular nearest.
+// info[i].setMiss();
+// info[i].setIndex(-1);
+// }
+// else
+// {
+// info[i].setIndex
+// (triIndexer.toGlobal(info[i].index()));
+// }
+// }
+// if (nProcs == 1)
+// {
+// // Fully local
+// nLocal++;
+// }
+// }
+// }
+// }
+//
+//
+// if
+// (
+// Pstream::parRun()
+// && (
+// returnReduce(nLocal, sumOp<label>())
+// < returnReduce(samples.size(), sumOp<label>())
+// )
+// )
+// {
+// // Not all can be resolved locally. Build queries and map, send over
+// // queries, do intersections, send back and merge.
+//
+// // Calculate queries and exchange map
+// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+// pointField allCentres;
+// scalarField allRadiusSqr;
+// labelList allSegmentMap;
+// autoPtr<mapDistribute> mapPtr
+// (
+// calcLocalQueries
+// (
+// false, // exclude local processor since already done above
+// samples,
+// nearestDistSqr,
+//
+// allCentres,
+// allRadiusSqr,
+// allSegmentMap
+// )
+// );
+// const mapDistribute& map = mapPtr();
+//
+//
+// // swap samples to local processor
+// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//
+// map.distribute(allCentres);
+// map.distribute(allRadiusSqr);
+//
+//
+// // Do my tests
+// // ~~~~~~~~~~~
+//
+// List<pointIndexHit> allInfo(allCentres.size());
+// forAll(allInfo, i)
+// {
+// allInfo[i] = octree.findNearest
+// (
+// allCentres[i],
+// allRadiusSqr[i]
+// );
+// if (allInfo[i].hit())
+// {
+// // We don't know if the nearest is on an edge/point. If
+// // this is the case we preferentially want to return the
+// // index on the processor that holds all surrounding triangles
+// // so we can do e.g. follow-on inside/outside tests
+// if
+// (
+// surfaceClosed_
+// && !contains
+// (
+// procBb_[Pstream::myProcNo()],
+// allInfo[i].hitPoint()
+// )
+// )
+// {
+// // Nearest point is not on local processor so the
+// // the triangle is only there because some other bit of it
+// // is on it. Assume there is another processor that holds
+// // the full surrounding of the triangle so we can clear
+// // this particular nearest.
+// allInfo[i].setMiss();
+// allInfo[i].setIndex(-1);
+// }
+// else
+// {
+// allInfo[i].setIndex
+// (
+// triIndexer.toGlobal(allInfo[i].index())
+// );
+// }
+// }
+// }
+//
+//
+// // Send back results
+// // ~~~~~~~~~~~~~~~~~
+//
+// map.reverseDistribute(allSegmentMap.size(), allInfo);
+//
+//
+// // Extract information
+// // ~~~~~~~~~~~~~~~~~~~
+//
+// forAll(allInfo, i)
+// {
+// if (allInfo[i].hit())
+// {
+// label pointi = allSegmentMap[i];
+//
+// if (!info[pointi].hit())
+// {
+// // No intersection yet so take this one
+// info[pointi] = allInfo[i];
+// }
+// else
+// {
+// // Nearest intersection
+// if
+// (
+// magSqr(allInfo[i].hitPoint()-samples[pointi])
+// < magSqr(info[pointi].hitPoint()-samples[pointi])
+// )
+// {
+// info[pointi] = allInfo[i];
+// }
+// }
+// }
+// }
+// }
+//}
+
+
+void Foam::distributedTriSurfaceMesh::findNearest
+(
+ const pointField& samples,
+ const scalarField& nearestDistSqr,
+ List<pointIndexHit>& info
+) const
{
- read();
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::findNearest(samples, nearestDistSqr, info);
+ return;
+ }
- bounds().reduce();
+ addProfiling
+ (
+ findNearest,
+ "distributedTriSurfaceMesh::findNearest"
+ );
if (debug)
{
- InfoInFunction << "Read distributedTriSurface from " << io.objectPath()
- << " and dictionary:" << endl;
- writeStats(Info);
+ Pout<< "distributedTriSurfaceMesh::findNearest :"
+ << " trying to find nearest for " << samples.size()
+ << " samples with max sphere "
+ << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
+ << endl;
+ }
- labelList nTris(Pstream::nProcs());
- nTris[Pstream::myProcNo()] = triSurface::size();
- Pstream::gatherList(nTris);
- Pstream::scatterList(nTris);
+ const globalIndex& triIndexer = globalTris();
- Info<< endl<< "\tproc\ttris\tbb" << endl;
- forAll(nTris, proci)
- {
- Info<< '\t' << proci << '\t' << nTris[proci]
- << '\t' << procBb_[proci] << endl;
- }
- Info<< endl;
- }
-}
+ // Two-pass searching:
+ // 1. send the sample to the processor whose bb contains it. This is
+ // most likely also the one that holds the nearest triangle. (In case
+ // there is no containing processor send to nearest processors. Note
+ // that this might cause a lot of traffic if this is likely)
+ // Send the resulting nearest point back.
+ // 2. with the find from 1 look at which other processors might have a
+ // better triangle. Since hopefully step 1) will have produced a tight
+ // bounding box this should limit the amount of points to be retested
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+ // 1. Test samples on processor(s) that contains them
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
-{
- clearOut();
-}
+ autoPtr<mapDistribute> map1Ptr;
+ scalarField distSqr(nearestDistSqr);
+ boolList procContains(Pstream::nProcs(), false);
+ boolList procOverlaps(Pstream::nProcs(), false);
+ label nOutside = 0;
+ {
+ List<DynamicList<label>> dynSendMap(Pstream::nProcs());
+ // Pre-size. Assume samples are uniformly distributed
+ forAll(dynSendMap, proci)
+ {
+ dynSendMap[proci].reserve(samples.size()/Pstream::nProcs());
+ }
-void Foam::distributedTriSurfaceMesh::clearOut()
-{
- globalTris_.clear();
- triSurfaceMesh::clearOut();
-}
+ forAll(samples, samplei)
+ {
+ label minProci = -1;
+ Tuple2<label, scalar> best = findBestProcs
+ (
+ samples[samplei],
+ distSqr[samplei],
+ procContains,
+ procOverlaps,
+ minProci
+ );
+ label nContains = 0;
+ forAll(procBb_, proci)
+ {
+ if (procContains[proci])
+ {
+ nContains++;
+ dynSendMap[proci].append(samplei);
+ distSqr[samplei] = best.second();
+ }
+ }
+ if (nContains == 0)
+ {
+ nOutside++;
+ // Sample is outside all bb. Choices:
+ // - send to all processors
+ // - send to single processor
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+ //forAll(procOverlaps[proci])
+ //{
+ // if (procOverlaps[proci])
+ // {
+ // dynSendMap[proci].append(samplei);
+ // distSqr[samplei] = best.second();
+ // }
+ //}
+ if (minProci != -1)
+ {
+ dynSendMap[minProci].append(samplei);
+ distSqr[samplei] = best.second();
+ }
+ }
+ }
-const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
-{
- if (!globalTris_.valid())
- {
- globalTris_.reset(new globalIndex(triSurface::size()));
+ labelListList sendMap(Pstream::nProcs());
+ forAll(sendMap, proci)
+ {
+ sendMap[proci].transfer(dynSendMap[proci]);
+ }
+ map1Ptr.set(new mapDistribute(std::move(sendMap)));
}
- return *globalTris_;
-}
+ const mapDistribute& map1 = map1Ptr();
-void Foam::distributedTriSurfaceMesh::findNearest
-(
- const pointField& samples,
- const scalarField& nearestDistSqr,
- List<pointIndexHit>& info
-) const
-{
- const indexedOctree<treeDataTriSurface>& octree = tree();
+ if (debug)
+ {
+ Pout<< "Pass1:"
+ << " of " << samples.size() << " samples sending to" << endl;
+ label nSend = 0;
+ forAll(map1.subMap(), proci)
+ {
+ Pout<< " " << proci << "\t" << map1.subMap()[proci].size()
+ << endl;
+ nSend += map1.subMap()[proci].size();
+ }
+ Pout<< " sum\t" << nSend << endl
+ << " outside\t" << nOutside << endl;
+ }
- // Important:force synchronised construction of indexing
- const globalIndex& triIndexer = globalTris();
+ List<nearestAndDist> nearestInfo;
+ {
+ // Get the points I need to test and test locally
+ pointField localPoints(samples);
+ map1.distribute(localPoints);
+ scalarField localDistSqr(distSqr);
+ map1.distribute(localDistSqr);
+ List<pointIndexHit> localInfo;
+ triSurfaceMesh::findNearest(localPoints, localDistSqr, localInfo);
+ convertTriIndices(localInfo);
+
+ // Pack into structure for combining information from multiple
+ // processors
+ nearestInfo.setSize(localInfo.size());
+ nearestInfo = nearestAndDist(pointIndexHit(), Foam::sqr(GREAT));
+
+ label nHit = 0;
+ label nIgnoredHit = 0;
+
+ forAll(nearestInfo, i)
+ {
+ const pointIndexHit& info = localInfo[i];
+ if (info.hit())
+ {
+ nHit++;
- // Initialise
- // ~~~~~~~~~~
+ if
+ (
+ surfaceClosed_
+ && !contains(procBb_[Pstream::myProcNo()], info.hitPoint())
+ )
+ {
+ // Nearest point is not on local processor so the
+ // the triangle is only there because some other bit
+ // of it is on it. Assume there is another processor that
+ // holds the full surrounding of the triangle so we can
+ // ignore this particular nearest.
+ nIgnoredHit++;
+ }
+ else
+ {
+ nearestAndDist& ni = nearestInfo[i];
+ ni.first() = info;
+ ni.second() = magSqr(localPoints[i]-info.hitPoint());
+ }
+ }
+ }
- info.setSize(samples.size());
- forAll(info, i)
- {
- info[i].setMiss();
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::findNearest :"
+ << " searched locally for " << localPoints.size()
+ << " samples with max sphere "
+ << (localDistSqr.size() ? Foam::sqrt(max(localDistSqr)) : Zero)
+ << " found hits:" << nHit
+ << " of which outside local bb:" << nIgnoredHit
+ << endl;
+ }
}
+ // Send back to originating processor. Choose best if sent to multiple
+ // processors. Note that afterwards all unused entries have the unique
+ // value nearestZero (distance < 0). This is used later on to see if
+ // the sample was sent to any processor.
+ mapDistributeBase::distribute
+ (
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(0),
+ samples.size(),
+ map1.constructMap(),
+ map1.constructHasFlip(),
+ map1.subMap(),
+ map1.subHasFlip(),
+ nearestInfo,
+ nearestEqOp(),
+ noOp(), // no flipping
+ nearestZero
+ );
- // Do any local queries
- // ~~~~~~~~~~~~~~~~~~~~
+ // 2. Test samples on other processor(s) that overlap
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- label nLocal = 0;
+ // Now we have (in nearestInfo) for every input sample the current best
+ // hit (on the processor that originates the sample). See if we can
+ // improve it by sending the queries to any other processors
+ autoPtr<mapDistribute> map2Ptr;
{
+ List<DynamicList<label>> dynSendMap(Pstream::nProcs());
+
// Work array - whether processor bb overlaps the bounding sphere.
boolList procBbOverlaps(Pstream::nProcs());
- forAll(samples, i)
+ label nFound = 0;
+
+ forAll(nearestInfo, samplei)
{
- // Find the processor this sample+radius overlaps.
- label nProcs = calcOverlappingProcs
+ const point& sample = samples[samplei];
+ const nearestAndDist& ni = nearestInfo[samplei];
+ const pointIndexHit& info = ni.first();
+
+ if (info.hit())
+ {
+ nFound++;
+ }
+
+ scalar d2 =
(
- samples[i],
- nearestDistSqr[i],
- procBbOverlaps
+ info.hit()
+ ? ni.second()
+ : distSqr[samplei]
+ );
+
+ label hitProci =
+ (
+ info.hit()
+ ? triIndexer.whichProcID(info.index())
+ : -1
);
- // Overlaps local processor?
- if (procBbOverlaps[Pstream::myProcNo()])
+ // Find the processors this sample+radius overlaps.
+ calcOverlappingProcs(sample, d2, procBbOverlaps);
+
+ forAll(procBbOverlaps, proci)
{
- info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
- if (info[i].hit())
+ if (procBbOverlaps[proci])
{
- info[i].setIndex(triIndexer.toGlobal(info[i].index()));
- }
- if (nProcs == 1)
- {
- // Fully local
- nLocal++;
+ // Check this sample wasn't already handled above. This
+ // could be improved since the sample might have been
+ // searched on multiple processors. We now only exclude the
+ // processor where the point was inside.
+ if (proci != hitProci)
+ {
+ dynSendMap[proci].append(samplei);
+ }
}
}
}
+
+ labelListList sendMap(Pstream::nProcs());
+ forAll(sendMap, proci)
+ {
+ sendMap[proci].transfer(dynSendMap[proci]);
+ }
+ map2Ptr.reset(new mapDistribute(std::move(sendMap)));
}
+ const mapDistribute& map2 = map2Ptr();
- if
+ if (debug)
+ {
+ Pout<< "Pass2:"
+ << " of " << samples.size() << " samples sending to" << endl;
+ label nSend = 0;
+ forAll(map2.subMap(), proci)
+ {
+ Pout<< " " << proci << "\t" << map2.subMap()[proci].size()
+ << endl;
+ nSend += map2.subMap()[proci].size();
+ }
+ Pout<< " sum\t" << nSend << endl;
+ }
+
+ // Send samples and current best distance
+ pointField localSamples(samples);
+ map2.distribute(localSamples);
+ scalarField localDistSqr(distSqr);
+ forAll(nearestInfo, samplei)
+ {
+ const nearestAndDist& ni = nearestInfo[samplei];
+ if (ni.first().hit())
+ {
+ localDistSqr[samplei] = ni.second();
+ }
+ }
+ map2.distribute(localDistSqr);
+
+ // Do local test
+ List<pointIndexHit> localInfo;
+ triSurfaceMesh::findNearest(localSamples, localDistSqr, localInfo);
+ convertTriIndices(localInfo);
+
+ // Pack and send back
+ List<nearestAndDist> localBest(localSamples.size());
+ label nHit = 0;
+ label nIgnoredHit = 0;
+ forAll(localInfo, i)
+ {
+ const pointIndexHit& info = localInfo[i];
+ if (info.hit())
+ {
+ nHit++;
+ if
+ (
+ surfaceClosed_
+ && !contains(procBb_[Pstream::myProcNo()], info.hitPoint())
+ )
+ {
+ // See above
+ nIgnoredHit++;
+ }
+ else
+ {
+ nearestAndDist& ni = localBest[i];
+ ni.first() = info;
+ ni.second() = magSqr(info.hitPoint()-localSamples[i]);
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::findNearest :"
+ << " searched locally for " << localSamples.size()
+ << " samples with max sphere "
+ << (localDistSqr.size() ? Foam::sqrt(max(localDistSqr)) : Zero)
+ << " found hits:" << nHit
+ << " of which outside local bb:" << nIgnoredHit
+ << endl;
+ }
+
+ mapDistributeBase::distribute
(
- Pstream::parRun()
- && (
- returnReduce(nLocal, sumOp<label>())
- < returnReduce(samples.size(), sumOp<label>())
- )
- )
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(0),
+ samples.size(),
+ map2.constructMap(),
+ map2.constructHasFlip(),
+ map2.subMap(),
+ map2.subHasFlip(),
+ localBest,
+ nearestEqOp(),
+ noOp(), // no flipping
+ nearestZero
+ );
+
+ // Combine with nearestInfo
+ info.setSize(samples.size());
+ forAll(samples, samplei)
+ {
+ nearestAndDist ni(nearestInfo[samplei]);
+ nearestEqOp()(ni, localBest[samplei]);
+
+ info[samplei] = ni.first();
+ }
+}
+
+
+void Foam::distributedTriSurfaceMesh::findNearest
+(
+ const pointField& samples,
+ const scalarField& nearestDistSqr,
+ const labelList& regionIndices,
+ List<pointIndexHit>& info
+) const
+{
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::findNearest
+ (
+ samples,
+ nearestDistSqr,
+ regionIndices,
+ info
+ );
+ return;
+ }
+
+ addProfiling
+ (
+ findNearestRegion,
+ "distributedTriSurfaceMesh::findNearestRegion"
+ );
+
+ if (debug)
{
- // Not all can be resolved locally. Build queries and map, send over
- // queries, do intersections, send back and merge.
+ Pout<< "distributedTriSurfaceMesh::findNearest :"
+ << " trying to find nearest and region for " << samples.size()
+ << " samples with max sphere "
+ << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
+ << endl;
+ }
+ if (regionIndices.empty())
+ {
+ findNearest(samples, nearestDistSqr, info);
+ }
+ else
+ {
// Calculate queries and exchange map
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1551,6 +3443,7 @@ void Foam::distributedTriSurfaceMesh::findNearest
(
calcLocalQueries
(
+ true, // also send to local processor
samples,
nearestDistSqr,
@@ -1573,16 +3466,37 @@ void Foam::distributedTriSurfaceMesh::findNearest
// ~~~~~~~~~~~
List<pointIndexHit> allInfo(allCentres.size());
+ triSurfaceMesh::findNearest
+ (
+ allCentres,
+ allRadiusSqr,
+ regionIndices,
+ allInfo
+ );
+ convertTriIndices(allInfo);
+
forAll(allInfo, i)
{
- allInfo[i] = octree.findNearest
- (
- allCentres[i],
- allRadiusSqr[i]
- );
if (allInfo[i].hit())
{
- allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
+ if
+ (
+ surfaceClosed_
+ && !contains
+ (
+ procBb_[Pstream::myProcNo()],
+ allInfo[i].hitPoint()
+ )
+ )
+ {
+ // Nearest point is not on local processor so the
+ // the triangle is only there because some other bit of it
+ // is on it. Assume there is another processor that holds
+ // the full surrounding of the triangle so we can clear
+ // this particular nearest.
+ allInfo[i].setMiss();
+ allInfo[i].setIndex(-1);
+ }
}
}
@@ -1632,13 +3546,20 @@ void Foam::distributedTriSurfaceMesh::findLine
List<pointIndexHit>& info
) const
{
- findLine
- (
- true, // nearestIntersection
- start,
- end,
- info
- );
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::findLine(start, end, info);
+ }
+ else
+ {
+ findLine
+ (
+ true, // nearestIntersection
+ start,
+ end,
+ info
+ );
+ }
}
@@ -1649,13 +3570,20 @@ void Foam::distributedTriSurfaceMesh::findLineAny
List<pointIndexHit>& info
) const
{
- findLine
- (
- true, // nearestIntersection
- start,
- end,
- info
- );
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::findLineAny(start, end, info);
+ }
+ else
+ {
+ findLine
+ (
+ true, // nearestIntersection
+ start,
+ end,
+ info
+ );
+ }
}
@@ -1666,6 +3594,25 @@ void Foam::distributedTriSurfaceMesh::findLineAll
List<List<pointIndexHit>>& info
) const
{
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::findLineAll(start, end, info);
+ return;
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::findLineAll :"
+ << " intersecting with "
+ << start.size() << " rays" << endl;
+ }
+
+ addProfiling
+ (
+ findLineAll,
+ "distributedTriSurfaceMesh::findLineAll"
+ );
+
// Reuse fineLine. We could modify all of findLine to do multiple
// intersections but this would mean a lot of data transferred so
// for now we just find nearest intersection and retest from that
@@ -1792,6 +3739,12 @@ void Foam::distributedTriSurfaceMesh::findLineAll
break;
}
}
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::findLineAll :"
+ << " finished intersecting with "
+ << start.size() << " rays" << endl;
+ }
}
@@ -1801,6 +3754,15 @@ void Foam::distributedTriSurfaceMesh::getRegion
labelList& region
) const
{
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getRegion :"
+ << " getting region for "
+ << info.size() << " triangles" << endl;
+ }
+
+ addProfiling(getRegion, "distributedTriSurfaceMesh::getRegion");
+
if (!Pstream::parRun())
{
region.setSize(info.size());
@@ -1815,10 +3777,17 @@ void Foam::distributedTriSurfaceMesh::getRegion
region[i] = -1;
}
}
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getRegion :"
+ << " finished getting region for "
+ << info.size() << " triangles" << endl;
+ }
+
return;
}
-
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
@@ -1852,6 +3821,13 @@ void Foam::distributedTriSurfaceMesh::getRegion
// ~~~~~~~~~~~~~~~~~
map.reverseDistribute(info.size(), region);
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getRegion :"
+ << " finished getting region for "
+ << info.size() << " triangles" << endl;
+ }
}
@@ -1867,6 +3843,15 @@ void Foam::distributedTriSurfaceMesh::getNormal
return;
}
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getNormal :"
+ << " getting normal for "
+ << info.size() << " triangles" << endl;
+ }
+
+ addProfiling(getNormal, "distributedTriSurfaceMesh::getNormal");
+
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
@@ -1901,6 +3886,345 @@ void Foam::distributedTriSurfaceMesh::getNormal
// ~~~~~~~~~~~~~~~~~
map.reverseDistribute(info.size(), normal);
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getNormal :"
+ << " finished getting normal for "
+ << info.size() << " triangles" << endl;
+ }
+}
+
+
+//void Foam::distributedTriSurfaceMesh::getVolumeTypeUncached
+//(
+// const pointField& samples,
+// List<volumeType>& volType
+//) const
+//{
+// if (!Pstream::parRun())
+// {
+// triSurfaceMesh::getVolumeType(samples, volType);
+// return;
+// }
+//
+//
+// if (!hasVolumeType())
+// {
+// FatalErrorInFunction
+// << "Volume type only supported for closed distributed surfaces."
+// << exit(FatalError);
+// }
+//
+// // Trigger (so parallel synchronised) construction of outside type.
+// // Normally this would get triggered from inside individual searches
+// // so would not be parallel synchronised
+// if (outsideVolType_ == volumeType::UNKNOWN)
+// {
+// // Determine nearest (in parallel)
+// const point outsidePt(bounds().max() + 0.5*bounds().span());
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::outsideVolumeType :"
+// << " triggering outsidePoint" << outsidePt
+// << " orientation" << endl;
+// }
+//
+// const pointField outsidePts(1, outsidePt);
+// List<pointIndexHit> nearestInfo;
+// findNearest
+// (
+// outsidePts,
+// scalarField(1, Foam::sqr(GREAT)),
+// nearestInfo
+// );
+//
+// List<volumeType> outsideVolTypes;
+// surfaceSide(outsidePts, nearestInfo, outsideVolTypes);
+// outsideVolType_ = outsideVolTypes[0];
+//
+// if (debug)
+// {
+// Pout<< "distributedTriSurfaceMesh::outsideVolumeType :"
+// << " determined outsidePoint" << outsidePt
+// << " to be " << volumeType::names[outsideVolType_] << endl;
+// }
+// }
+//
+// // Determine nearest (in parallel)
+// List<pointIndexHit> nearestInfo(samples.size());
+// findNearest
+// (
+// samples,
+// scalarField(samples.size(), Foam::sqr(GREAT)),
+// nearestInfo
+// );
+//
+// // Determine orientation (in parallel)
+// surfaceSide(samples, nearestInfo, volType);
+//}
+
+
+void Foam::distributedTriSurfaceMesh::getVolumeType
+(
+ const pointField& samples,
+ List<volumeType>& volType
+) const
+{
+ if (!Pstream::parRun())
+ {
+ triSurfaceMesh::getVolumeType(samples, volType);
+ return;
+ }
+
+
+ if (!hasVolumeType())
+ {
+ FatalErrorInFunction
+ << "Volume type only supported for closed distributed surfaces."
+ << exit(FatalError);
+ }
+
+ // Trigger (so parallel synchronised) construction of outside type.
+ // Normally this would get triggered from inside individual searches
+ // so would not be parallel synchronised
+ if (outsideVolType_ == volumeType::UNKNOWN)
+ {
+ addProfiling
+ (
+ getVolumeType,
+ "distributedTriSurfaceMesh::getCachedVolumeType"
+ );
+
+ // Determine nearest (in parallel)
+ const point outsidePt(bounds().max() + 0.5*bounds().span());
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " triggering outsidePoint" << outsidePt
+ << " orientation" << endl;
+ }
+
+ const pointField outsidePts(1, outsidePt);
+ List<pointIndexHit> nearestInfo;
+ findNearest
+ (
+ outsidePts,
+ scalarField(1, Foam::sqr(GREAT)),
+ nearestInfo
+ );
+
+ List<volumeType> outsideVolTypes;
+ surfaceSide(outsidePts, nearestInfo, outsideVolTypes);
+ outsideVolType_ = outsideVolTypes[0];
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " determined outsidePoint" << outsidePt
+ << " to be " << volumeType::names[outsideVolType_] << endl;
+ }
+
+ if
+ (
+ outsideVolType_ == volumeType::INSIDE
+ || outsideVolType_ == volumeType::OUTSIDE
+ )
+ {
+ // Get local tree
+ const indexedOctree<treeDataTriSurface>& t = tree();
+ PackedList<2>& nt = t.nodeTypes();
+ const List<indexedOctree<treeDataTriSurface>::node>& nodes =
+ t.nodes();
+ nt.setSize(nodes.size());
+ nt = volumeType::UNKNOWN;
+
+ // Collect midpoints
+ DynamicField<point> midPoints(label(0.5*nodes.size()));
+ collectLeafMids(0, midPoints);
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " triggering orientation caching for "
+ << midPoints.size() << " leaf mids" << endl;
+ }
+
+ // Get volume type of mid points
+ List<volumeType> midVolTypes;
+ getVolumeType(midPoints, midVolTypes);
+
+ // Cache on local tree
+ label index = 0;
+ calcVolumeType
+ (
+ midVolTypes,
+ index,
+ nt,
+ 0 // nodeI
+ );
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " done orientation caching for "
+ << midPoints.size() << " leaf mids" << endl;
+ }
+ }
+ }
+
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " finding orientation for " << samples.size()
+ << " samples" << endl;
+ }
+
+ addProfiling
+ (
+ getVolumeType,
+ "distributedTriSurfaceMesh::getVolumeType"
+ );
+
+
+ DynamicList<label> outsideSamples;
+
+ // Distribute samples to relevant processors
+ autoPtr<mapDistribute> mapPtr;
+ {
+ labelListList sendMap(Pstream::nProcs());
+ {
+ // 1. Count
+ labelList nSend(Pstream::nProcs(), 0);
+ forAll(samples, samplei)
+ {
+ // Find the processors this sample overlaps.
+ label nOverlap = 0;
+ forAll(procBb_, proci)
+ {
+ if (contains(procBb_[proci], samples[samplei]))
+ {
+ nSend[proci]++;
+ nOverlap++;
+ }
+ }
+
+ // Special case: point is outside all bbs. These would not
+ // get sent to anyone so handle locally. Note that is the
+ // equivalent of the test in triSurfaceMesh against the local
+ // tree bb
+ if (nOverlap == 0)
+ {
+ outsideSamples.append(samplei);
+ }
+ }
+
+ forAll(nSend, proci)
+ {
+ sendMap[proci].setSize(nSend[proci]);
+ }
+ nSend = 0;
+
+ // 2. Fill
+ forAll(samples, samplei)
+ {
+ // Find the processors this sample overlaps.
+ forAll(procBb_, proci)
+ {
+ if (contains(procBb_[proci], samples[samplei]))
+ {
+ labelList& procSend = sendMap[proci];
+ procSend[nSend[proci]++] = samplei;
+ }
+ }
+ }
+ }
+
+ mapPtr.set(new mapDistribute(std::move(sendMap)));
+ }
+ const mapDistribute& map = mapPtr();
+
+ // Get the points I need to test
+ pointField localPoints(samples);
+ map.distribute(localPoints);
+
+ volType.setSize(localPoints.size());
+ volType = volumeType::UNKNOWN;
+
+ // Split the local queries into those that I can look up on the tree and
+ // those I need to search the nearest for
+ DynamicField<point> fullSearchPoints(localPoints.size());
+ DynamicList<label> fullSearchMap(localPoints.size());
+ forAll(localPoints, i)
+ {
+ volType[i] = cachedVolumeType(0, localPoints[i]);
+ if (volType[i] == volumeType::UNKNOWN)
+ {
+ fullSearchMap.append(i);
+ fullSearchPoints.append(localPoints[i]);
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " original samples:" << samples.size()
+ << " resulting in local queries:"
+ << localPoints.size()
+ << " of which cached:" << localPoints.size()-fullSearchPoints.size()
+ << endl;
+ }
+
+ // Determine nearest (in parallel)
+ List<pointIndexHit> nearestInfo;
+ findNearest
+ (
+ fullSearchPoints,
+ scalarField(fullSearchPoints.size(), Foam::sqr(GREAT)),
+ nearestInfo
+ );
+
+ // Determine orientation (in parallel)
+ List<volumeType> fullSearchType;
+ surfaceSide(fullSearchPoints, nearestInfo, fullSearchType);
+
+ // Combine
+ forAll(fullSearchMap, i)
+ {
+ volType[fullSearchMap[i]] = fullSearchType[i];
+ }
+
+ // Send back to originator. In case of multiple answers choose inside or
+ // outside
+ const volumeType zero(volumeType::UNKNOWN);
+ mapDistributeBase::distribute
+ (
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(0),
+ samples.size(),
+ map.constructMap(),
+ map.constructHasFlip(),
+ map.subMap(),
+ map.subHasFlip(),
+ volType,
+ volumeCombineOp(),
+ noOp(), // no flipping
+ zero
+ );
+
+
+ // Add the points outside the bounding box
+ for (label samplei : outsideSamples)
+ {
+ volType[samplei] = outsideVolType_;
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getVolumeType :"
+ << " finished finding orientation for " << samples.size()
+ << " samples" << endl;
+ }
}
@@ -1916,6 +4240,15 @@ void Foam::distributedTriSurfaceMesh::getField
return;
}
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getField :"
+ << " retrieving field for "
+ << info.size() << " triangles" << endl;
+ }
+
+ addProfiling(getField, "distributedTriSurfaceMesh::getField");
+
const auto* fldPtr = findObject<triSurfaceLabelField>("values");
if (fldPtr)
@@ -1954,33 +4287,26 @@ void Foam::distributedTriSurfaceMesh::getField
map.reverseDistribute(info.size(), values);
}
-}
-
-void Foam::distributedTriSurfaceMesh::getVolumeType
-(
- const pointField& points,
- List<volumeType>& volType
-) const
-{
- FatalErrorInFunction
- << "Volume type not supported for distributed surfaces."
- << exit(FatalError);
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::getField :"
+ << " finished retrieving field for "
+ << info.size() << " triangles" << endl;
+ }
}
-// Subset the part of surface that is overlapping bb.
-Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
+void Foam::distributedTriSurfaceMesh::overlappingSurface
(
const triSurface& s,
const List<treeBoundBox>& bbs,
-
- labelList& subPointMap,
- labelList& subFaceMap
+ boolList& includedFace
)
{
// Determine what triangles to keep.
- boolList includedFace(s.size(), false);
+ includedFace.setSize(s.size());
+ includedFace = false;
// Create a slightly larger bounding box.
List<treeBoundBox> bbsX(bbs.size());
@@ -2006,7 +4332,22 @@ Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
includedFace[trii] = true;
}
}
+}
+
+
+// Subset the part of surface that is overlapping bb.
+Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
+(
+ const triSurface& s,
+ const List<treeBoundBox>& bbs,
+ labelList& subPointMap,
+ labelList& subFaceMap
+)
+{
+ // Determine what triangles to keep.
+ boolList includedFace;
+ overlappingSurface(s, bbs, includedFace);
return subsetMesh(s, includedFace, subPointMap, subFaceMap);
}
@@ -2019,6 +4360,22 @@ void Foam::distributedTriSurfaceMesh::distribute
autoPtr<mapDistribute>& pointMap
)
{
+ if (!Pstream::parRun())
+ {
+ return;
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::distribute :"
+ << " distributing surface according to method:"
+ << distributionTypeNames_[distType_]
+ << " follow bbs:" << flatOutput(bbs) << endl;
+ }
+
+ addProfiling(distribute, "distributedTriSurfaceMesh::distribute");
+
+
// Get bbs of all domains
// ~~~~~~~~~~~~~~~~~~~~~~
@@ -2038,6 +4395,11 @@ void Foam::distributedTriSurfaceMesh::distribute
break;
case INDEPENDENT:
+ case DISTRIBUTED:
+ if (currentDistType_ == distType_)
+ {
+ return;
+ }
newProcBb = independentlyDistributedBbs(*this);
break;
@@ -2097,13 +4459,6 @@ void Foam::distributedTriSurfaceMesh::distribute
pointSendMap[proci],
faceSendMap[proci]
);
-
- if (debug)
- {
- //Pout<< "Overlapping with proc " << proci
- // << " faces:" << faceSendMap[proci].size()
- // << " points:" << pointSendMap[proci].size() << endl << endl;
- }
}
if (keepNonLocal)
@@ -2147,15 +4502,8 @@ void Foam::distributedTriSurfaceMesh::distribute
// Send over how many faces/points i need to receive
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- labelListList faceSendSizes(Pstream::nProcs());
- faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
- forAll(faceSendMap, proci)
- {
- faceSendSizes[Pstream::myProcNo()][proci] = faceSendMap[proci].size();
- }
- Pstream::gatherList(faceSendSizes);
- Pstream::scatterList(faceSendSizes);
-
+ labelList faceRecvSizes;
+ Pstream::exchangeSizes(faceSendMap, faceRecvSizes);
// Exchange surfaces
@@ -2204,11 +4552,11 @@ void Foam::distributedTriSurfaceMesh::distribute
PstreamBuffers pBufs(Pstream::defaultCommsType);
- forAll(faceSendSizes, proci)
+ forAll(faceSendMap, proci)
{
if (proci != Pstream::myProcNo())
{
- if (faceSendSizes[Pstream::myProcNo()][proci] > 0)
+ if (faceSendMap[proci].size() > 0)
{
UOPstream str(proci, pBufs);
@@ -2222,15 +4570,6 @@ void Foam::distributedTriSurfaceMesh::distribute
pointMap
)
);
-
- //if (debug)
- //{
- // Pout<< "Sending to " << proci
- // << " faces:" << faceSendMap[proci].size()
- // << " points:" << subSurface.points().size() << endl
- // << endl;
- //}
-
str << subSurface;
}
}
@@ -2242,25 +4581,17 @@ void Foam::distributedTriSurfaceMesh::distribute
// Receive and merge all
// ~~~~~~~~~~~~~~~~~~~~~
- forAll(faceSendSizes, proci)
+ forAll(faceRecvSizes, proci)
{
if (proci != Pstream::myProcNo())
{
- if (faceSendSizes[proci][Pstream::myProcNo()] > 0)
+ if (faceRecvSizes[proci] > 0)
{
UIPstream str(proci, pBufs);
// Receive
triSurface subSurface(str);
- //if (debug)
- //{
- // Pout<< "Received from " << proci
- // << " faces:" << subSurface.size()
- // << " points:" << subSurface.points().size() << endl
- // << endl;
- //}
-
// Merge into allSurf
merge
(
@@ -2273,12 +4604,6 @@ void Foam::distributedTriSurfaceMesh::distribute
faceConstructMap[proci],
pointConstructMap[proci]
);
-
- //if (debug)
- //{
- // Pout<< "Current merged surface : faces:" << allTris.size()
- // << " points:" << allPoints.size() << endl << endl;
- //}
}
}
}
@@ -2306,11 +4631,13 @@ void Foam::distributedTriSurfaceMesh::distribute
// Construct triSurface. Reuse storage.
triSurface::operator=(triSurface(allTris, patches(), allPoints, true));
+ // Clear trees, preserve topological info (surfaceClosed, outsidePointType)
clearOut();
// Set the bounds() value to the boundBox of the undecomposed surface
bounds() = boundBox(points(), true);
+ currentDistType_ = distType_;
// Regions stays same
// Volume type stays same.
@@ -2338,20 +4665,49 @@ void Foam::distributedTriSurfaceMesh::distribute
}
Info<< endl;
- OBJstream str(searchableSurface::time().path()/"after.obj");
- Info<< "Writing local bounding box to " << str.name() << endl;
- const List<treeBoundBox>& myBbs = procBb_[Pstream::myProcNo()];
- forAll(myBbs, i)
+ if (debug & 2)
+ {
+ OBJstream str(searchableSurface::time().path()/"after.obj");
+ Info<< "Writing local bounding box to " << str.name() << endl;
+ const List<treeBoundBox>& myBbs = procBb_[Pstream::myProcNo()];
+ forAll(myBbs, i)
+ {
+ pointField pts(myBbs[i].points());
+ const edgeList& es = treeBoundBox::edges;
+ forAll(es, ei)
+ {
+ const edge& e = es[ei];
+ str.write(linePointRef(pts[e[0]], pts[e[1]]));
+ }
+ }
+ }
+ if (debug & 2)
{
- pointField pts(myBbs[i].points());
- const edgeList& es = treeBoundBox::edges;
- forAll(es, ei)
+ OBJstream str(searchableSurface::time().path()/"after_all.obj");
+ Info<< "Writing all bounding boxes to " << str.name() << endl;
+ for (auto myBbs : procBb_)
{
- const edge& e = es[ei];
- str.write(linePointRef(pts[e[0]], pts[e[1]]));
+ forAll(myBbs, i)
+ {
+ pointField pts(myBbs[i].points());
+ const edgeList& es = treeBoundBox::edges;
+ forAll(es, ei)
+ {
+ const edge& e = es[ei];
+ str.write(linePointRef(pts[e[0]], pts[e[1]]));
+ }
+ }
}
}
}
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::distribute :"
+ << " done distributing surface according to method:"
+ << distributionTypeNames_[distType_]
+ << " follow bbs:" << flatOutput(bbs) << endl;
+ }
}
@@ -2363,6 +4719,12 @@ bool Foam::distributedTriSurfaceMesh::writeObject
const bool valid
) const
{
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::writeObject :"
+ << " writing surface valid:" << valid << endl;
+ }
+
// Make sure dictionary goes to same directory as surface
const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
@@ -2387,7 +4749,15 @@ bool Foam::distributedTriSurfaceMesh::writeObject
}
// Dictionary needs to be written in ascii - binary output not supported.
- return dict_.writeObject(IOstream::ASCII, ver, cmp, true);
+ bool ok = dict_.writeObject(IOstream::ASCII, ver, cmp, true);
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::writeObject :"
+ << " done writing surface" << endl;
+ }
+
+ return ok;
}
@@ -2401,7 +4771,10 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
os << "Triangles : " << returnReduce(triSurface::size(), sumOp<label>())
<< endl
<< "Vertices : " << returnReduce(nPoints, sumOp<label>()) << endl
- << "Bounding Box : " << bb << endl;
+ << "Bounding Box : " << bb << endl
+ << "Closed : " << surfaceClosed_ << endl
+ << "Outside point: " << volumeType::names[outsideVolType_] << endl
+ << "Distribution : " << distributionTypeNames_[distType_] << endl;
}
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
index 174819bda46f279e929a114ba0ee3aba4887f82d..95ffdce4ea81bd1ca73606d3b60f07bd97b5e376 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 2015-2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -42,6 +42,15 @@ Description
more communication.
- frozen : no change
+ Note: addressing used:
+ distributedTriSurfaceMesh: none
+
+ triSurfaceMesh:
+ - surf.pointFaces() : edge addressing (for volume tests only)
+ - surf.edges() : edgeTree
+ - surf.faceFaces() : only if minQuality > 0
+
+
SourceFiles
distributedTriSurfaceMesh.C
@@ -55,6 +64,7 @@ SourceFiles
#include "IOdictionary.H"
#include "Pair.H"
#include "globalIndex.H"
+#include "DynamicField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -86,7 +96,8 @@ public:
{
FOLLOW = 0,
INDEPENDENT = 1,
- FROZEN = 2
+ DISTRIBUTED = 2,
+ FROZEN = 3
};
static const Enum<distributionType> distributionTypeNames_;
@@ -112,14 +123,22 @@ private:
//- Global triangle numbering
mutable autoPtr<globalIndex> globalTris_;
- //- The distribution type.
+ //- The (wanted) distribution type.
distributionType distType_;
+ //- The (current) distribution type. Used to trigger re-distribution
+ // when starting from undecomposed surface.
+ distributionType currentDistType_;
+
// Private Member Functions
// Read
+ //- Search for io.local directory (=triSurface) either in case
+ // directory or in parent directory
+ static word findLocalInstance(const IOobject& io);
+
//- Read my additional data
bool read();
@@ -183,6 +202,9 @@ private:
// Triangle index
+ //- Helper: convert local triangle indices to global ones
+ void convertTriIndices(List<pointIndexHit>& info) const;
+
//- Obtains global indices from pointIndexHit and swaps them back
// to their original processor. Used to calculate local region
// and normal.
@@ -195,6 +217,24 @@ private:
// Nearest
+ //- Is location inside any of the bounding boxes
+ bool contains
+ (
+ const List<treeBoundBox>& bbs,
+ const point& sample
+ ) const;
+
+ //- Is location inside any of the processors bb or if not
+ // does it overlap
+ Tuple2<label, scalar> findBestProcs
+ (
+ const point& centre,
+ const scalar radiusSqr,
+ boolList& procContains,
+ boolList& procOverlaps,
+ label& minProci
+ ) const;
+
label calcOverlappingProcs
(
const point& centre,
@@ -202,8 +242,10 @@ private:
boolList& overlaps
) const;
+ //- Calculate map to send centres+radius to processors
autoPtr<mapDistribute> calcLocalQueries
(
+ const bool includeLocalProcessor,
const pointField& centres,
const scalarField& radiusSqr,
@@ -213,8 +255,71 @@ private:
) const;
+ // Side
+
+ //- Side of nearest point w.r.t. edge between face0 and face1
+ volumeType edgeSide
+ (
+ const point& sample,
+ const point& nearestPoint,
+ const label face0,
+ const label face1
+ ) const;
+
+ //- Find edge-connected face
+ label findOtherFace
+ (
+ const labelListList& pointFaces,
+ const label nearFacei,
+ const label nearLabel
+ ) const;
+
+ //- Side of nearest point on faces w.r.t. samples. Handles nearest
+ // on edge/point
+ void surfaceSide
+ (
+ const pointField& samples,
+ const List<pointIndexHit>& nearestInfo,
+ List<volumeType>& region
+ ) const;
+
+
+ // Caching of volume type (based on indexedOctree)
+
+ //- Collect mid points of tree boxes
+ void collectLeafMids
+ (
+ const label nodeI,
+ DynamicField<point>& midPoints
+ ) const;
+
+ //- Find volume type of tree boxes
+ volumeType calcVolumeType
+ (
+ const List<volumeType>& midPointTypes,
+ label& index,
+ PackedList<2>& nodeTypes,
+ const label nodeI
+ ) const;
+
+ //- Look up any cached data. Return unknown if cannot be determined.
+ volumeType cachedVolumeType
+ (
+ const label nodeI,
+ const point& sample
+ ) const;
+
+
// Surface redistribution
+ //- Calculate face-faces
+ static void calcFaceFaces
+ (
+ const triSurface& s,
+ const labelListList& pointFaces,
+ labelListList& faceFaces
+ );
+
//- Finds new bounds based on an independent decomposition.
List<List<treeBoundBox>> independentlyDistributedBbs
(
@@ -344,13 +449,6 @@ public:
// searchableSurface implementation
- //- Whether supports volume type below. I.e. whether is closed.
- // Not supported.
- virtual bool hasVolumeType() const
- {
- return false;
- }
-
//- Range of global indices that can be returned.
virtual label globalSize() const
{
@@ -364,6 +462,16 @@ public:
List<pointIndexHit>&
) const;
+ //- Find the nearest locations for the supplied points to a
+ // particular region in the searchable surface.
+ virtual void findNearest
+ (
+ const pointField& samples,
+ const scalarField& nearestDistSqr,
+ const labelList& regionIndices,
+ List<pointIndexHit>& info
+ ) const;
+
virtual void findLine
(
const pointField& start,
@@ -434,6 +542,14 @@ public:
// indices) get the specified field. Misses do not get set.
virtual void getField(const List<pointIndexHit>&, labelList&) const;
+ //- Calculate the triangles that are overlapping bounds.
+ static void overlappingSurface
+ (
+ const triSurface&,
+ const List<treeBoundBox>&,
+ boolList& includedFace
+ );
+
//- Subset the part of surface that is overlapping bounds.
static triSurface overlappingSurface
(
diff --git a/tutorials/incompressible/simpleFoam/motorBike/Allrun b/tutorials/incompressible/simpleFoam/motorBike/Allrun
index e0e6038798e66f16a82874be7f45e11207de7a1c..4654c95afe1538e2aac9e307c42edbf5ecf53ea6 100755
--- a/tutorials/incompressible/simpleFoam/motorBike/Allrun
+++ b/tutorials/incompressible/simpleFoam/motorBike/Allrun
@@ -19,7 +19,8 @@ runApplication $decompDict decomposePar
if foamDictionary -entry geometry -value system/snappyHexMeshDict | \
grep -q distributedTriSurfaceMesh
then
- runParallel $decompDict surfaceRedistributePar motorBike.obj independent
+ echo "surfaceRedistributePar does not need to be run anymore"
+ echo " - distributedTriSurfaceMesh will do on-the-fly redistribution"
fi
runParallel $decompDict snappyHexMesh -overwrite
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allclean b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..a13b9c2ab8de5780209cdab19f8c13f5b07c87c5
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allclean
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions # Tutorial clean functions
+
+cleanCase
+
+(cd constant/triSurface && rm -f box_12_*.obj box.obj box_trans.obj)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allrun b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..f3be65de24cf8279299329cbb0c1cb909311a06f
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/Allrun
@@ -0,0 +1,53 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions
+
+runApplication blockMesh
+
+# Create fine surface
+runApplication -s 1 surfaceRefineRedGreen \
+ constant/triSurface/box_12.obj constant/triSurface/box_12_1.obj
+
+runApplication -s 2 surfaceRefineRedGreen \
+ constant/triSurface/box_12_1.obj constant/triSurface/box_12_2.obj
+
+runApplication -s 3 surfaceRefineRedGreen \
+ constant/triSurface/box_12_2.obj constant/triSurface/box_12_3.obj
+
+runApplication -s 4 surfaceRefineRedGreen \
+ constant/triSurface/box_12_3.obj constant/triSurface/box_12_4.obj
+
+runApplication -s 5 surfaceRefineRedGreen \
+ constant/triSurface/box_12_4.obj constant/triSurface/box_12_5.obj
+
+runApplication -s 6 surfaceRefineRedGreen \
+ constant/triSurface/box_12_5.obj constant/triSurface/box_12_6.obj
+
+runApplication -s 7 surfaceRefineRedGreen \
+ constant/triSurface/box_12_6.obj constant/triSurface/box_12_7.obj
+
+runApplication -s 8 surfaceRefineRedGreen \
+ constant/triSurface/box_12_7.obj constant/triSurface/box_12_8.obj
+
+runApplication -s 9 surfaceRefineRedGreen \
+ constant/triSurface/box_12_8.obj constant/triSurface/box_12_9.obj
+
+runApplication -s 10 surfaceRefineRedGreen \
+ constant/triSurface/box_12_9.obj constant/triSurface/box.obj
+
+
+runApplication surfaceTransformPoints \
+ -rotate '((1 0 0)(1 0.8 0.9))' \
+ -translate '(0.1 0.101 0.103)' \
+ constant/triSurface/box.obj constant/triSurface/box_trans.obj
+
+
+# Run non-parallel
+runApplication snappyHexMesh
+
+# Run parallel
+runApplication decomposePar
+
+runParallel -s parallel snappyHexMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/constant/triSurface/box_12.obj b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/constant/triSurface/box_12.obj
new file mode 100644
index 0000000000000000000000000000000000000000..d3e271f814cec7fa6407c22283f9fdce8c4e65ec
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/constant/triSurface/box_12.obj
@@ -0,0 +1,34 @@
+# Wavefront OBJ file written 2018-09-13T12:59:07
+o box_12
+
+# points : 8
+# faces : 12
+# zones : 1
+# 0 patch0 (nFaces: 12)
+
+# <points count="8">
+v -0.420753 -0.63726 -0.408494
+v 0.228766 -0.76226 0.341506
+v 0.420753 0.63726 0.408494
+v -0.228766 0.76226 -0.341506
+v -0.204247 -0.0122595 0.841506
+v 0.204247 0.0122595 -0.841506
+v -0.853766 0.11274 0.0915064
+v 0.853766 -0.11274 -0.0915064
+# </points>
+
+# <faces count="12">
+g patch0
+f 7 1 5
+f 2 5 1
+f 8 6 3
+f 4 3 6
+f 8 2 6
+f 1 6 2
+f 7 5 4
+f 3 4 5
+f 4 6 7
+f 1 7 6
+f 2 8 5
+f 3 5 8
+# </faces>
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/blockMeshDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..961bf87e809177ffeb2d101ff0158d159db2b44c
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/blockMeshDict
@@ -0,0 +1,98 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale 1;
+
+vertices
+(
+ (-1 -1 -1)
+ ( 1 -1 -1)
+ ( 1 1 -1)
+ (-1 1 -1)
+ (-1 -1 1)
+ ( 1 -1 1)
+ ( 1 1 1)
+ (-1 1 1)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) (4 4 4) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ maxY
+ {
+ type patch;
+ faces
+ (
+ (3 7 6 2)
+ );
+ }
+
+ minX
+ {
+ type patch;
+ faces
+ (
+ (0 4 7 3)
+ );
+ }
+
+ maxX
+ {
+ type patch;
+ faces
+ (
+ (2 6 5 1)
+ );
+ }
+
+ minY
+ {
+ type patch;
+ faces
+ (
+ (1 5 4 0)
+ );
+ }
+
+ minZ
+ {
+ type patch;
+ faces
+ (
+ (0 3 2 1)
+ );
+ }
+
+ maxZ
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ );
+ }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/controlDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a0d442b8909a7a9d79a3acd7d70717aca262f19c
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/controlDict
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 s |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+ //triSurfaceMesh 1;
+ //distributedTriSurfaceMesh 1;
+ //triSurface 1;
+ //PrimitivePatch 1;
+}
+
+application snappyHexMesh;
+
+startFrom startTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 2000;
+
+deltaT 1;
+
+writeControl timeStep;
+
+writeInterval 100;
+
+purgeWrite 0;
+
+writeFormat binary;
+
+writePrecision 6;
+
+writeCompression off;
+
+timeFormat general;
+
+timePrecision 6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/decomposeParDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..daa976064a2064ca6d48ab3fa7b248ebef6dd0bf
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/decomposeParDict
@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method hierarchical;
+
+coeffs
+{
+ n (2 2 2);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSchemes b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..8ef0651ea170f9a6308a9f948564b315776737ee
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSchemes
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default steadyState;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ default none;
+
+ div(phi,U) bounded Gauss upwind;
+ div(phi,T) bounded Gauss upwind;
+ div(phi,k) bounded Gauss upwind;
+ div(phi,epsilon) bounded Gauss upwind;
+ div(phi,R) bounded Gauss upwind;
+ div(R) Gauss linear;
+ div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear limited corrected 0.333;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default limited corrected 0.333;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSolution b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..b8ea24dd90019cdaa866298c9dbbb61e9e57c375
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/fvSolution
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ p_rgh
+ {
+ solver PCG;
+ preconditioner DIC;
+ tolerance 1e-08;
+ relTol 0.01;
+ }
+
+ "(U|T|k|epsilon)"
+ {
+ solver PBiCGStab;
+ preconditioner DILU;
+ tolerance 1e-07;
+ relTol 0.1;
+ }
+}
+
+SIMPLE
+{
+ nNonOrthogonalCorrectors 2;
+ pRefCell 0;
+ pRefValue 0;
+
+ residualControl
+ {
+ p_rgh 1e-2;
+ U 1e-4;
+ T 1e-3;
+
+ // possibly check turbulence fields
+ "(k|epsilon|omega)" 1e-3;
+ }
+}
+
+relaxationFactors
+{
+ fields
+ {
+ p_rgh 0.7;
+ }
+ equations
+ {
+ U 0.2;
+ T 0.5;
+ "(k|epsilon)" 0.7;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/meshQualityDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/meshQualityDict
new file mode 100644
index 0000000000000000000000000000000000000000..f94681c006da542be44fe95a58ebe082ff601116
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/meshQualityDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object meshQualityDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include defaults parameters from master dictionary
+#includeEtc "caseDicts/meshQualityDict"
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..d54c3807f570e1a9063d16edc86e806b0f740f7e
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
@@ -0,0 +1,287 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v1806 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Which of the steps to run
+castellatedMesh true;
+snap true;
+addLayers false;
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
+geometry
+{
+ box
+ {
+ type distributedTriSurfaceMesh;
+ file "box.obj"; //"box_12_2.obj";
+ }
+ box_trans
+ {
+ file "box_trans.obj";
+ type distributedTriSurfaceMesh;
+ }
+};
+
+
+
+// Settings for the castellatedMesh generation.
+castellatedMeshControls
+{
+
+ // Refinement parameters
+ // ~~~~~~~~~~~~~~~~~~~~~
+
+ // If local number of cells is >= maxLocalCells on any processor
+ // switches from from refinement followed by balancing
+ // (current method) to (weighted) balancing before refinement.
+ maxLocalCells 100000;
+
+ // Overall cell limit (approximately). Refinement will stop immediately
+ // upon reaching this number so a refinement level might not complete.
+ // Note that this is the number of cells before removing the part which
+ // is not 'visible' from the keepPoint. The final number of cells might
+ // actually be a lot less.
+ maxGlobalCells 20000000;
+
+ // The surface refinement loop might spend lots of iterations refining just a
+ // few cells. This setting will cause refinement to stop if <= minimumRefine
+ // are selected for refinement. Note: it will at least do one iteration
+ // (unless the number of cells to refine is 0)
+ minRefinementCells 0;
+
+ // Number of buffer layers between different levels.
+ // 1 means normal 2:1 refinement restriction, larger means slower
+ // refinement.
+ nCellsBetweenLevels 1;
+
+
+
+ // Explicit feature edge refinement
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // Specifies a level for any cell intersected by its edges.
+ // This is a featureEdgeMesh, read from constant/triSurface for now.
+ features
+ (
+ );
+
+
+
+ // Surface based refinement
+ // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // Specifies two levels for every surface. The first is the minimum level,
+ // every cell intersecting a surface gets refined up to the minimum level.
+ // The second level is the maximum level. Cells that 'see' multiple
+ // intersections where the intersections make an
+ // angle > resolveFeatureAngle get refined up to the maximum level.
+
+ refinementSurfaces
+ {
+ box
+ {
+ // Surface-wise min and max refinement level
+ level (4 4);
+ }
+ box_trans
+ {
+ level (5 5);
+ faceZone box_trans;
+ cellZone box_trans;
+ cellZoneInside inside;
+ }
+ }
+
+ // Resolve sharp angles
+ resolveFeatureAngle 60;
+
+
+ // Region-wise refinement
+ // ~~~~~~~~~~~~~~~~~~~~~~
+
+ // Specifies refinement level for cells in relation to a surface. One of
+ // three modes
+ // - distance. 'levels' specifies per distance to the surface the
+ // wanted refinement level. The distances need to be specified in
+ // descending order.
+ // - inside. 'levels' is only one entry and only the level is used. All
+ // cells inside the surface get refined up to the level. The surface
+ // needs to be closed for this to be possible.
+ // - outside. Same but cells outside.
+
+ refinementRegions
+ {
+ }
+
+
+ // Mesh selection
+ // ~~~~~~~~~~~~~~
+
+ // After refinement patches get added for all refinementSurfaces and
+ // all cells intersecting the surfaces get put into these patches. The
+ // section reachable from the locationInMesh is kept.
+ // NOTE: This point should never be on a face, always inside a cell, even
+ // after refinement.
+ locationInMesh (1e-5 1e-5 1e-5);
+
+
+ // Whether any faceZones (as specified in the refinementSurfaces)
+ // are only on the boundary of corresponding cellZones or also allow
+ // free-standing zone faces. Not used if there are no faceZones.
+ allowFreeStandingZoneFaces false;
+}
+
+
+
+// Settings for the snapping.
+snapControls
+{
+ //- Number of patch smoothing iterations before finding correspondence
+ // to surface
+ nSmoothPatch 3;
+
+ //- Relative distance for points to be attracted by surface feature point
+ // or edge. True distance is this factor times local
+ // maximum edge length.
+ tolerance 1.0;
+
+ //- Number of mesh displacement relaxation iterations.
+ nSolveIter 30;
+
+ //- Maximum number of snapping relaxation iterations. Should stop
+ // before upon reaching a correct mesh.
+ nRelaxIter 5;
+
+
+ // Feature snapping
+
+ //- Number of feature edge snapping iterations.
+ // Leave out altogether to disable.
+ nFeatureSnapIter 2;
+
+ //- Detect (geometric) features by sampling the surface (default=false)
+ implicitFeatureSnap true;
+
+ //- Use castellatedMeshControls::features (default = true)
+ explicitFeatureSnap false;
+}
+
+
+
+// Settings for the layer addition.
+addLayersControls
+{
+ // Are the thickness parameters below relative to the undistorted
+ // size of the refined cell outside layer (true) or absolute sizes (false).
+ relativeSizes true;
+
+ // Per final patch (so not geometry!) the layer information
+ layers
+ {
+ }
+
+ // Expansion factor for layer mesh
+ expansionRatio 1.0;
+
+ // Wanted thickness of final added cell layer. If multiple layers
+ // is the thickness of the layer furthest away from the wall.
+ // Relative to undistorted size of cell outside layer.
+ // is the thickness of the layer furthest away from the wall.
+ // See relativeSizes parameter.
+ finalLayerThickness 0.5;
+
+ // Minimum thickness of cell layer. If for any reason layer
+ // cannot be above minThickness do not add layer.
+ // Relative to undistorted size of cell outside layer.
+ // See relativeSizes parameter.
+ minThickness 0.25;
+
+ // If points get not extruded do nGrow layers of connected faces that are
+ // also not grown. This helps convergence of the layer addition process
+ // close to features.
+ // Note: changed(corrected) w.r.t 1.7.x! (didn't do anything in 1.7.x)
+ nGrow 0;
+
+
+ // Advanced settings
+
+ // When not to extrude surface. 0 is flat surface, 90 is when two faces
+ // are perpendicular
+ featureAngle 60;
+
+ // Maximum number of snapping relaxation iterations. Should stop
+ // before upon reaching a correct mesh.
+ nRelaxIter 5;
+
+ // Number of smoothing iterations of surface normals
+ nSmoothSurfaceNormals 1;
+
+ // Number of smoothing iterations of interior mesh movement direction
+ nSmoothNormals 3;
+
+ // Smooth layer thickness over surface patches
+ nSmoothThickness 10;
+
+ // Stop layer growth on highly warped cells
+ maxFaceThicknessRatio 0.5;
+
+ // Reduce layer growth where ratio thickness to medial
+ // distance is large
+ maxThicknessToMedialRatio 0.3;
+
+ // Angle used to pick up medial axis points
+ // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
+ minMedianAxisAngle 90;
+
+ // Create buffer region for new layer terminations
+ nBufferCellsNoExtrude 0;
+
+
+ // Overall max number of layer addition iterations. The mesher will exit
+ // if it reaches this number of iterations; possibly with an illegal
+ // mesh.
+ nLayerIter 50;
+}
+
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+ #include "meshQualityDict"
+
+ // Advanced
+
+ //- Number of error distribution iterations
+ nSmoothScale 4;
+ //- amount to scale back displacement at error points
+ errorReduction 0.75;
+}
+
+
+// Advanced
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1e-6;
+
+// ************************************************************************* //