diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index 198be01bcf71f133a77380de31ef7b6ba5ec5ff6..a6547fa1142cc9ce4ce14ffe8a302042d99090c7 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
- Copyright (C) 2015-2022 OpenCFD Ltd.
+ Copyright (C) 2015-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -168,12 +168,7 @@ void Foam::shellSurfaces::orient()
{
const searchableSurface& s = allGeometry_[shells_[shellI]];
- if
- (
- modes_[shellI] != DISTANCE
- && isA<triSurfaceMesh>(s)
- && !isA<distributedTriSurfaceMesh>(s)
- )
+ if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
{
hasSurface = true;
@@ -196,35 +191,39 @@ void Foam::shellSurfaces::orient()
{
const searchableSurface& s = allGeometry_[shells_[shellI]];
- if
- (
- modes_[shellI] != DISTANCE
- && isA<triSurfaceMesh>(s)
- && !isA<distributedTriSurfaceMesh>(s)
- )
+ if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
{
- triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
+ List<pointIndexHit> info;
+ vectorField normal;
+ labelList region;
+ s.findNearest
(
- refCast<const triSurfaceMesh>(s)
+ pointField(1, outsidePt),
+ scalarField(1, GREAT),
+ info,
+ normal,
+ region
);
- // Flip surface so outsidePt is outside.
- bool anyFlipped = orientedSurface::orient
- (
- shell,
- outsidePt,
- true
- );
+ //Pout<< "outsidePt:" << outsidePt << endl;
+ //Pout<< "info :" << info[0] << endl;
+ //Pout<< "normal :" << normal[0] << endl;
+ //Pout<< "region :" << region[0] << endl;
- if (anyFlipped && !dryRun_)
+ bool anyFlipped = false;
+ if ((normal[0] & (info[0].point()-outsidePt)) > 0)
{
- // orientedSurface will have done a clearOut of the surface.
- // we could do a clearout of the triSurfaceMeshes::trees()
- // but these aren't affected by orientation
- // (except for cached
- // sideness which should not be set at this point.
- // !!Should check!)
+ triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
+ (
+ refCast<const triSurfaceMesh>(s)
+ );
+ shell.flip();
+ anyFlipped = true;
+ }
+
+ if (anyFlipped && !dryRun_)
+ {
Info<< "shellSurfaces : Flipped orientation of surface "
<< s.name()
<< " so point " << outsidePt << " is outside." << endl;
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
index 8976887905924d6c8f0fe559955959f44af87b7a..4d13ca4deb75bfa8ff9db1b22e93045bb79a2a45 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2015-2020,2022 OpenCFD Ltd.
+ Copyright (C) 2015-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -101,6 +101,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
if (debug)
{
Pout<< "triSurfaceMesh::isSurfaceClosed:"
+ << " surface:" << searchableSurface::name()
<< " determining closedness for surface with "
<< triSurface::size() << " triangles" << endl;
}
@@ -223,6 +224,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
if (debug)
{
Pout<< "triSurfaceMesh::isSurfaceClosed :"
+ << " surface:" << searchableSurface::name()
<< " surface is non-manifold" << endl;
}
return false;
@@ -240,6 +242,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
if (debug)
{
Pout<< "triSurfaceMesh::isSurfaceClosed :"
+ << " surface:" << searchableSurface::name()
<< " surface is open" << endl;
}
return false;
@@ -274,6 +277,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
if (debug)
{
Pout<< "triSurfaceMesh::isSurfaceClosed :"
+ << " surface:" << searchableSurface::name()
<< " surface is closed" << endl;
}
return true;
@@ -738,6 +742,7 @@ void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
if (debug)
{
Pout<< "triSurfaceMesh::movePoints :"
+ << " surface:" << searchableSurface::name()
<< " moving at time " << objectRegistry::time().timeName()
<< endl;
}
@@ -773,6 +778,7 @@ Foam::triSurfaceMesh::edgeTree() const
if (debug)
{
Pout<< "triSurfaceMesh::edgeTree :"
+ << " surface:" << searchableSurface::name()
<< " constructing tree for " << nEdges() - nInternalEdges()
<< " boundary edges" << endl;
}
@@ -882,6 +888,7 @@ Foam::volumeType Foam::triSurfaceMesh::outsideVolumeType() const
if (debug)
{
Pout<< "triSurfaceMesh::outsideVolumeType :"
+ << " surface:" << searchableSurface::name()
<< " triggering outsidePoint:" << outsidePt
<< " orientation" << endl;
}
@@ -906,6 +913,49 @@ Foam::volumeType Foam::triSurfaceMesh::outsideVolumeType() const
}
+void Foam::triSurfaceMesh::flip()
+{
+ if (debug)
+ {
+ Pout<< "triSurfaceMesh::flip :"
+ << " surface:" << searchableSurface::name()
+ << " with current orientation "
+ << volumeType::names[outsideVolType_]
+ << endl;
+ }
+
+ // Don't bother getting nearest etc. Just flip the triangles.
+
+ // triSurface
+ {
+ triSurface& s = *this;
+ s.clearOut();
+ for (auto& tri : s)
+ {
+ tri.flip();
+ }
+ }
+
+ // triSurfaceRegionSearch (if cached volume type)
+ triSurfaceRegionSearch::flip();
+
+ // edge tree not relevant
+
+ if (hasVolumeType())
+ {
+ // outsideVolType_
+ if (outsideVolType_ == volumeType::INSIDE)
+ {
+ outsideVolType_ = volumeType::OUTSIDE;
+ }
+ else if (outsideVolType_ == volumeType::OUTSIDE)
+ {
+ outsideVolType_ = volumeType::INSIDE;
+ }
+ }
+}
+
+
void Foam::triSurfaceMesh::findNearest
(
const pointField& samples,
@@ -916,6 +966,7 @@ void Foam::triSurfaceMesh::findNearest
if (debug)
{
Pout<< "triSurfaceMesh::findNearest :"
+ << " surface:" << searchableSurface::name()
<< " trying to find nearest for " << samples.size()
<< " samples with max sphere "
<< (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
@@ -942,6 +993,7 @@ void Foam::triSurfaceMesh::findNearest
if (debug)
{
Pout<< "triSurfaceMesh::findNearest :"
+ << " surface:" << searchableSurface::name()
<< " trying to find nearest and region for " << samples.size()
<< " samples with max sphere "
<< (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
@@ -973,6 +1025,7 @@ void Foam::triSurfaceMesh::findLine
if (debug)
{
Pout<< "triSurfaceMesh::findLine :"
+ << " surface:" << searchableSurface::name()
<< " intersecting with "
<< start.size() << " rays" << endl;
}
@@ -996,6 +1049,7 @@ void Foam::triSurfaceMesh::findLineAny
if (debug)
{
Pout<< "triSurfaceMesh::findLineAny :"
+ << " surface:" << searchableSurface::name()
<< " intersecting with "
<< start.size() << " rays" << endl;
}
@@ -1019,6 +1073,7 @@ void Foam::triSurfaceMesh::findLineAll
if (debug)
{
Pout<< "triSurfaceMesh::findLineAll :"
+ << " surface:" << searchableSurface::name()
<< " intersecting with "
<< start.size() << " rays" << endl;
}
@@ -1041,6 +1096,7 @@ void Foam::triSurfaceMesh::getRegion
if (debug)
{
Pout<< "triSurfaceMesh::getRegion :"
+ << " surface:" << searchableSurface::name()
<< " getting region for "
<< info.size() << " triangles" << endl;
}
@@ -1074,6 +1130,7 @@ void Foam::triSurfaceMesh::getNormal
if (debug)
{
Pout<< "triSurfaceMesh::getNormal :"
+ << " surface:" << searchableSurface::name()
<< " getting normal for "
<< info.size() << " triangles" << endl;
}
@@ -1182,6 +1239,7 @@ void Foam::triSurfaceMesh::setField(const labelList& values)
if (debug)
{
Pout<< "triSurfaceMesh::setField :"
+ << " surface:" << searchableSurface::name()
<< " finished setting field for "
<< values.size() << " triangles" << endl;
}
@@ -1213,6 +1271,7 @@ void Foam::triSurfaceMesh::getField
if (debug)
{
Pout<< "triSurfaceMesh::setField :"
+ << " surface:" << searchableSurface::name()
<< " finished getting field for "
<< info.size() << " triangles" << endl;
}
@@ -1231,6 +1290,7 @@ void Foam::triSurfaceMesh::getVolumeType
if (debug)
{
Pout<< "triSurfaceMesh::getVolumeType :"
+ << " surface:" << searchableSurface::name()
<< " finding orientation for " << points.size()
<< " samples" << endl;
}
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
index e5a18bdeb4bbc417cb1c817a79c28dfb58b64cd8..c833e24d2b5b70cc4cd0362140d4ab6a6fe73abf 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
@@ -216,6 +216,9 @@ public:
return triSurface::size();
}
+ //- Flip triangles, outsideVolumeType and all cached inside/outside.
+ virtual void flip();
+
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual tmp<pointField> coordinates() const;
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
index cb8a67db5cc6ae0f4c5970f5a3a56a020b559b3b..eae21fb06e97b34b8b04ec50e7f136d9e5ec8c23 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
@@ -251,4 +251,26 @@ void Foam::triSurfaceRegionSearch::findNearest
}
+void Foam::triSurfaceRegionSearch::flip()
+{
+ triSurfaceSearch::flip();
+
+ for (auto& tree : treeByRegion_)
+ {
+ PackedList<2>& nodeTypes = tree.nodeTypes();
+ forAll(nodeTypes, i)
+ {
+ if (nodeTypes[i] == volumeType::INSIDE)
+ {
+ nodeTypes[i] = volumeType::OUTSIDE;
+ }
+ else if (nodeTypes[i] == volumeType::OUTSIDE)
+ {
+ nodeTypes[i] = volumeType::INSIDE;
+ }
+ }
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
index 724617d1321275ddc266ca0f305fc89170da2e6e..54143930eb55c1f451714210884f6ce5e59d7b99 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
@@ -128,6 +128,11 @@ public:
const labelList& regionIndices,
List<pointIndexHit>& info
) const;
+
+ // Edit
+
+ //- Flip orientation
+ void flip();
};
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
index 33fa142a6149cdb8b53c7d8d097a4e7dbe654ea6..b80944a4910ecf014befa9464c2f0c1bcb2889ca 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
@@ -255,6 +255,26 @@ Foam::triSurfaceSearch::tree() const
}
+void Foam::triSurfaceSearch::flip()
+{
+ if (treePtr_)
+ {
+ PackedList<2>& nodeTypes = treePtr_->nodeTypes();
+ forAll(nodeTypes, i)
+ {
+ if (nodeTypes[i] == volumeType::INSIDE)
+ {
+ nodeTypes[i] = volumeType::OUTSIDE;
+ }
+ else if (nodeTypes[i] == volumeType::OUTSIDE)
+ {
+ nodeTypes[i] = volumeType::INSIDE;
+ }
+ }
+ }
+}
+
+
// Determine inside/outside for samples
Foam::boolList Foam::triSurfaceSearch::calcInside
(
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
index a4ce7c8c986b716d38b1878f7b76dab7eb27470d..2deb9abf9e2ad5e2a0469a8acf18f8279dfc0222 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
@@ -126,6 +126,9 @@ public:
//- Demand driven construction of the octree
const indexedOctree<treeDataTriSurface>& tree() const;
+ //- Flip orientation (if cached on octree)
+ void flip();
+
//- Return reference to the surface.
const triSurface& surface() const
{
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 8c32aface9d13a59d0e1510d601f641950f01d5e..bdb63f152e13184d6d624babf7be5be013612aa8 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -110,6 +110,7 @@ namespace Foam
};
+ ////- Same as above but with some more debugging
//typedef Tuple2<Pair<point>, volumeType> NearType;
//
////- Combine operator for volume types
@@ -337,73 +338,145 @@ Foam::word Foam::distributedTriSurfaceMesh::findLocalInstance
}
-// Read my additional data from the dictionary
-bool Foam::distributedTriSurfaceMesh::read()
+bool Foam::distributedTriSurfaceMesh::readSettings(const bool undecomposed)
{
// Get bb of all domains.
procBb_.resize_nocopy(Pstream::nProcs());
- if (dict_.empty())
+ if (!dict_.readIfPresent("bounds", procBb_[Pstream::myProcNo()]))
{
- // Did not find the distributed version; assume master has loaded the
- // triSurfaceMesh version. Make up some settings.
-
+ // Not found. Use local triangles' bounding box
const boundBox& localBb = triSurfaceMesh::bounds();
procBb_[Pstream::myProcNo()] =
treeBoundBoxList(1, treeBoundBox(localBb));
dict_.add("bounds", procBb_[Pstream::myProcNo()]);
+ }
+
- // Wanted distribution type
+ decomposeUsingBbs_ = dict_.getOrDefault<bool>("decomposeUsingBbs", true);
+
+ // Wanted distribution type
+ if
+ (
+ !distributionTypeNames_.readIfPresent
+ (
+ "distributionType",
+ dict_,
+ distType_
+ )
+ )
+ {
+ // Not found. Make up most sensible data
distType_ = DISTRIBUTED; //INDEPENDENT;
dict_.add("distributionType", distributionTypeNames_[distType_]);
+ }
- // Merge distance
+ // Merge distance
+ if (!dict_.readIfPresent("mergeDistance", mergeDist_))
+ {
mergeDist_ = SMALL;
dict_.add("mergeDistance", mergeDist_);
+ }
- // Force underlying triSurfaceMesh to calculate volume type
- // (is topological walk; does not construct tree)
- surfaceClosed_ = triSurfaceMesh::hasVolumeType();
- Pstream::broadcast(surfaceClosed_);
+ // Is closed?
+ if (!dict_.readIfPresent<label>("closed", surfaceClosed_))
+ {
+ if (undecomposed)
+ {
+ surfaceClosed_ = triSurfaceMesh::hasVolumeType();
+ Pstream::broadcast(surfaceClosed_);
+ }
+ else
+ {
+ surfaceClosed_ = false;
+ }
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.centre());
- // List<volumeType> outsideVolTypes;
- // triSurfaceMesh::getVolumeType
- // (
- // pointField(1, outsidePt),
- // outsideVolTypes
- // );
- // outsideVolType_ = outsideVolTypes[0];
- //}
- //dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
+ // Do we know volume inside/outside?
+ volumeType::type vt;
+ if
+ (
+ volumeType::names.readIfPresent
+ (
+ "outsideVolumeType",
+ dict_,
+ vt
+ )
+ )
+ {
+ outsideVolType_ = vt;
}
else
{
- dict_.readEntry("bounds", procBb_[Pstream::myProcNo()]);
+ if (undecomposed && surfaceClosed_)
+ {
+ if (UPstream::master())
+ {
+ // Determine nearest to do inside/outside testing.
+ // Do linear search to avoid building tree on this large
+ // undecomposed surface. See e.g. triSurfaceTools::surfaceSide
+
+ const triSurface& surf = *this;
+ const pointField& points = surf.points();
+ const boundBox& localBb = triSurfaceMesh::bounds();
+ const point outsidePt(localBb.max()+localBb.centre());
+
+ label nearTrii = -1;
+ pointHit nearInfo; // distance=GREAT
+ forAll(surf, trii)
+ {
+ pointHit info = surf[trii].nearestPoint(outsidePt, points);
- // Wanted distribution type
- distType_ = distributionTypeNames_.get("distributionType", dict_);
+ if (info.distance() < nearInfo.distance())
+ {
+ nearTrii = trii;
+ nearInfo = info;
+ }
+ }
- // Merge distance
- dict_.readEntry("mergeDistance", mergeDist_);
+ if (nearTrii != -1)
+ {
+ const triSurfaceTools::sideType t =
+ triSurfaceTools::surfaceSide
+ (
+ surf,
+ outsidePt,
+ nearTrii
+ );
- // Distribution type
- surfaceClosed_ = dict_.getOrDefault<bool>("closed", false);
+ if (t == triSurfaceTools::UNKNOWN)
+ {
+ outsideVolType_ = volumeType::UNKNOWN;
+ }
+ else if (t == triSurfaceTools::INSIDE)
+ {
+ outsideVolType_ = volumeType::INSIDE;
+ }
+ else if (t == triSurfaceTools::OUTSIDE)
+ {
+ outsideVolType_ = volumeType::OUTSIDE;
+ }
+ else
+ {
+ FatalErrorInFunction << "problem" << abort(FatalError);
+ }
+ }
+ }
+ // Scatter master
+ label volType = outsideVolType_;
+ Pstream::broadcast(volType);
+ outsideVolType_ = volumeType(volType);
+ }
+ else
+ {
+ // Mark as to-be-done
+ outsideVolType_ = volumeType::UNKNOWN;
+ }
- outsideVolType_ = volumeType::names.getOrDefault
- (
- "outsideVolumeType",
- dict_,
- volumeType::UNKNOWN
- );
+ dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
}
Pstream::allGatherList(procBb_);
@@ -412,6 +485,41 @@ bool Foam::distributedTriSurfaceMesh::read()
}
+void Foam::distributedTriSurfaceMesh::calcVertexNormals
+(
+ const triSurface& surf,
+ List<FixedList<vector, 3>>& vn
+) const
+{
+ // For now don't do features so just use the point normals
+ vectorField pointNormals(surf.points().size(), vector::zero);
+ forAll(surf, facei)
+ {
+ const triSurface::face_type& f = surf[facei];
+ const vector n = f.unitNormal(surf.points());
+ forAll(f, fp)
+ {
+ const label pointi = f[fp];
+ pointNormals[pointi] += n;
+ }
+ }
+ pointNormals.normalise();
+
+
+ vn.resize_nocopy(surf.size());
+ vn = Zero;
+
+ forAll(surf, facei)
+ {
+ const triSurface::face_type& f = surf[facei];
+ forAll(f, fp)
+ {
+ vn[facei][fp] = pointNormals[f[fp]];
+ }
+ }
+}
+
+
// Is segment fully local?
bool Foam::distributedTriSurfaceMesh::isLocal
(
@@ -504,28 +612,31 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
List<DynamicList<label>>& sendMap
) const
{
- // 1. Fully local already handled outside. Note: retest is cheap.
- if (isLocal(procBb_[Pstream::myProcNo()], start, end))
+ if (decomposeUsingBbs_)
{
- return;
- }
+ // 1. Fully local already handled outside. Note: retest is cheap.
+ if (isLocal(procBb_[Pstream::myProcNo()], start, end))
+ {
+ return;
+ }
- // 2. If fully inside one other processor, then only need to send
- // to that one processor even if it intersects another. Rare occurrence
- // but cheap to test.
- forAll(procBb_, proci)
- {
- if (proci != Pstream::myProcNo())
+ // 2. If fully inside one other processor, then only need to send
+ // to that one processor even if it intersects another. Rare occurrence
+ // but cheap to test.
+ forAll(procBb_, proci)
{
- const List<treeBoundBox>& bbs = procBb_[proci];
-
- if (isLocal(bbs, start, end))
+ if (proci != Pstream::myProcNo())
{
- sendMap[proci].append(allSegments.size());
- allSegmentMap.append(segmenti);
- allSegments.append(segment(start, end));
- return;
+ const List<treeBoundBox>& bbs = procBb_[proci];
+
+ if (isLocal(bbs, start, end))
+ {
+ sendMap[proci].append(allSegments.size());
+ allSegmentMap.append(segmenti);
+ allSegments.append(segment(start, end));
+ return;
+ }
}
}
}
@@ -644,7 +755,7 @@ void Foam::distributedTriSurfaceMesh::findLine
(
const bool nearestIntersection,
const pointField& start,
- const pointField& end,
+ const pointField& initialEnd,
List<pointIndexHit>& info
) const
{
@@ -659,6 +770,9 @@ void Foam::distributedTriSurfaceMesh::findLine
const indexedOctree<treeDataTriSurface>& octree = tree();
+ // end-point - updated with local hits
+ pointField end(initialEnd);
+
// Initialise
info.setSize(start.size());
forAll(info, i)
@@ -691,6 +805,8 @@ void Foam::distributedTriSurfaceMesh::findLine
if (info[i].hit())
{
info[i].setIndex(triIndexer.toGlobal(info[i].index()));
+ // Update endpoint
+ end[i] = info[i].hitPoint();
}
nLocal++;
}
@@ -699,113 +815,121 @@ void Foam::distributedTriSurfaceMesh::findLine
if
(
- returnReduce(nLocal, sumOp<label>())
- < returnReduce(start.size(), sumOp<label>())
+ decomposeUsingBbs_
+ && (
+ 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.
+ // In decomposeUsingBbs_ mode we're guaranteed to have all the geometry
+ // inside the local bb. Since we've done all tests we're done.
+ return;
+ }
- // Construct queries (segments)
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Not all can be resolved locally. Build segments and map,
+ // send over segments, do intersections, send back and merge.
- // Segments to test
- List<segment> allSegments(start.size());
- // Original index of segment
- labelList allSegmentMap(start.size());
- const autoPtr<mapDistribute> mapPtr
+ // Construct queries (segments)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ // Segments to test
+ List<segment> allSegments(start.size());
+ // Original index of segment
+ labelList allSegmentMap(start.size());
+
+ 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())
{
- if (!hitInfo.hit())
+ // No intersection yet so take this one
+ hitInfo = allInfo;
+ }
+ else if (nearestIntersection)
+ {
+ // Nearest intersection
+ if
+ (
+ start[segmenti].distSqr(allInfo.point())
+ < start[segmenti].distSqr(hitInfo.point())
+ )
{
- // No intersection yet so take this one
hitInfo = allInfo;
}
- else if (nearestIntersection)
- {
- // Nearest intersection
- if
- (
- start[segmenti].distSqr(allInfo.point())
- < start[segmenti].distSqr(hitInfo.point())
- )
- {
- hitInfo = allInfo;
- }
- }
}
}
}
@@ -1154,7 +1278,7 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::edgeSide
const label face1
) const
{
- const triSurface& surf = static_cast<const triSurface&>(*this);
+ const triSurface& surf = *this;
const pointField& points = surf.points();
// Compare to bisector. This is actually correct since edge is
@@ -1182,7 +1306,7 @@ Foam::label Foam::distributedTriSurfaceMesh::findOtherFace
const label nearLabel
) const
{
- const triSurface& surf = static_cast<const triSurface&>(*this);
+ const triSurface& surf = *this;
const triSurface::face_type& nearF = surf[nearFacei];
const edge e(nearF[nearLabel], nearF[nearF.fcIndex(nearLabel)]);
@@ -1255,7 +1379,8 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
Pout<< "distributedTriSurfaceMesh::surfaceSide :"
<< " on surface " << searchableSurface::name()
<< " finding surface side given points on surface for "
- << samples.size() << " samples" << endl;
+ << samples.size() << " samples"
+ << " cached vertexnormals " << vertexNormals_.valid() << endl;
}
// Use global index to send local tri and nearest back to originating
@@ -1276,6 +1401,10 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
pointField localSamples(samples);
map.distribute(localSamples);
+ // Send over nearest point since we have it anyway
+ List<pointIndexHit> localNearestInfo(nearestInfo);
+ map.distribute(localNearestInfo);
+
// Do my tests
// ~~~~~~~~~~~
@@ -1283,8 +1412,43 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
volType.setSize(triangleIndex.size());
volType = volumeType::UNKNOWN;
- const triSurface& surf = static_cast<const triSurface&>(*this);
+ const triSurface& surf = *this;
const pointField& points = surf.points();
+
+
+ if (vertexNormals_)
+ {
+ // Use smooth interpolation for the normal
+
+ forAll(triangleIndex, i)
+ {
+ const label facei = triangleIndex[i];
+ const triSurface::face_type& f = surf[facei];
+ const auto& vn = vertexNormals_()[facei];
+
+ const point& sample = localSamples[i];
+
+ const point& nearest = localNearestInfo[i].point();
+
+ const vector sampleNearestVec(sample - nearest);
+
+ const barycentric2D w(f.tri(points).pointToBarycentric(nearest));
+
+ const vector n = w[0]*vn[0]+w[1]*vn[1]+w[2]*vn[2];
+
+ const scalar c = (sampleNearestVec & n);
+
+ if (c > 0)
+ {
+ volType[i] = volumeType::OUTSIDE;
+ }
+ else
+ {
+ volType[i] = volumeType::INSIDE;
+ }
+ }
+ }
+ else
{
//const labelListList& pointFaces = surf.pointFaces();
// Construct pointFaces. Let's hope surface has compact point
@@ -1313,7 +1477,7 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
// Nearest to face interior. Use faceNormal to determine side
//scalar c = sampleNearestVec & surf.faceNormals()[facei];
- scalar c = sampleNearestVec & surf[facei].areaNormal(points);
+ const scalar c = (sampleNearestVec & f.areaNormal(points));
if (c > 0)
{
@@ -1451,13 +1615,45 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
// send the volume-query to a single processor)
- //forAll(localSamples, i)
//{
- // Pout<< "surfaceSide : for localSample:" << localSamples[i]
- // << " found volType:" << volumeType::names[volType[i]]
- // << endl;
+ //- Alternative which also sends over source sample so we get better
+ //- debug messages.
+ // // Pack sample, nearest and volType so we get nicer error messages
+ // typedef Tuple2<Pair<point>, volumeType> NearType;
+ // List<NearType> nearTypes(volType.size());
+ // forAll(localSamples, i)
+ // {
+ // const point& sample = localSamples[i];
+ // const point& near = localNearestInfo[i].point();
+ // nearTypes[i] = NearType(Pair<point>(sample, near), volType[i]);
+ // }
+ //
+ //
+ // const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
+ // mapDistributeBase::distribute
+ // (
+ // Pstream::commsTypes::nonBlocking,
+ // List<labelPair>::null(),
+ // nearestInfo.size(),
+ // map.constructMap(),
+ // map.constructHasFlip(),
+ // map.subMap(),
+ // map.subHasFlip(),
+ // nearTypes,
+ // zero,
+ // NearTypeCombineOp(),
+ // noOp(), // no flipping
+ // UPstream::msgType(),
+ // map.comm()
+ // );
+ // volType.setSize(nearTypes.size());
+ // forAll(nearTypes, i)
+ // {
+ // volType[i] = nearTypes[i].second();
+ // }
//}
+
const volumeType zero(volumeType::UNKNOWN);
mapDistributeBase::distribute
(
@@ -1476,44 +1672,10 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
map.comm()
);
- //// Pack sample, nearest and volType so we get nicer error messages
- ////typedef Tuple2<Pair<point>, volumeType>> NearType;
- //List<NearType> nearTypes(volType.size());
- //forAll(localSamples, i)
- //{
- // const point& sample = localSamples[i];
- // const point& near = nearest[i];
- // nearTypes[i] = NearType(Pair<point>(sample, near), volType[i]);
- //}
- //
- //
- //const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
- //mapDistributeBase::distribute
- //(
- // Pstream::commsTypes::nonBlocking,
- // List<labelPair>::null(),
- // nearestInfo.size(),
- // map.constructMap(),
- // map.constructHasFlip(),
- // map.subMap(),
- // map.subHasFlip(),
- // nearTypes,
- // zero,
- // NearTypeCombineOp(),
- // noOp(), // no flipping
- // UPstream::msgType(),
- // map.comm()
- //);
- //volType.setSize(nearTypes.size());
- //forAll(nearTypes, i)
- //{
- // volType[i] = nearTypes[i].second();
- //}
-
if (debug)
{
Pout<< "distributedTriSurfaceMesh::surfaceSide :"
- << " finished finding surface" << searchableSurface::name()
+ << " finished finding surface " << searchableSurface::name()
<< " given points on surface "
<< searchableSurface::name() << " for "
<< samples.size() << " samples" << endl;
@@ -1521,6 +1683,196 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
}
+Foam::volumeType Foam::distributedTriSurfaceMesh::markMixed
+(
+ const indexedOctree<treeDataTriSurface>& tree,
+ const label nodei,
+ const triPointRef& tri,
+ PackedList<2>& nodeTypes
+) const
+{
+ const auto& nod = tree.nodes()[nodei];
+ //const auto& contents = tree.contents();
+
+ //Pout<< indent
+ // << "Entering node:" << nodei
+ // << " with bb:" << nod.bb_ << endl;
+ //Pout<< incrIndent;
+
+ volumeType myType = volumeType::UNKNOWN;
+
+ for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
+ {
+ const labelBits index = nod.subNodes_[octant];
+ const treeBoundBox subBb(nod.bb_.subBbox(octant));
+
+ if (subBb.intersects(tri))
+ {
+ volumeType subType;
+
+ if (indexedOctree<treeDataTriSurface>::isNode(index))
+ {
+ //Pout<< indent
+ // << "Testing node octant:" << octant << " with node:"
+ // << indexedOctree<treeDataTriSurface>::getNode(index)
+ // << " subBb:" << subBb << endl;
+
+ // tree node. Recurse.
+ subType = markMixed
+ (
+ tree,
+ indexedOctree<treeDataTriSurface>::getNode(index),
+ tri,
+ nodeTypes
+ );
+ }
+ else if (indexedOctree<treeDataTriSurface>::isContent(index))
+ {
+ // Content inside so we cannot pre-determine whether
+ // whole node box is inside or outside
+
+ //const label contenti =
+ // indexedOctree<treeDataTriSurface>::getContent(index);
+ //Pout<< indent
+ // << "Testing content octant:" << octant
+ // << " with contenti:" << contenti << " subBb:" << subBb
+ // << " holding shapes:" << flatOutput(contents[contenti])
+ // << endl;
+
+ // Store octant type
+ subType = volumeType::MIXED;
+ }
+ else
+ {
+ // Nothing in it but intersects triangle. Can't cache
+
+ //Pout<< indent
+ // << "Testing empty octant:" << octant
+ // << " subBb:" << subBb << endl;
+
+ // Store octant type
+ subType = volumeType::MIXED;
+ }
+
+ // Store octant type
+ nodeTypes.set(labelBits::pack(nodei, 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)
+ {
+ //Pout<< indent
+ // << "for octant:" << octant
+ // << " siwtched from myType:" << myType
+ // << " and subType:" << subType
+ // << " switched to MIXED" << endl;
+ myType = volumeType::MIXED;
+ }
+ }
+ }
+
+ //Pout<< indent
+ // << "Result for node:" << nodei << " myType:" << myType << endl;
+ //Pout<< decrIndent;
+
+ return myType;
+}
+
+
+void Foam::distributedTriSurfaceMesh::markMixedOverlap
+(
+ const indexedOctree<treeDataTriSurface>& tree,
+ PackedList<2>& nodeTypes
+) const
+{
+ // At this point each triangle is on a unique processor so we can just
+ // go and
+ // check what additional triangles would be sent over based on the
+ // bounding boxes. Note that the surface at this point is already
+ // uniquely distributed (every triangle is on a single processor
+ // only) and the tree is only for those.
+
+ PstreamBuffers pBufs;
+
+ // Test & send
+ forAll(procBb_, proci)
+ {
+ if (proci != UPstream::myProcNo())
+ {
+ labelList pointMap;
+ labelList faceMap;
+ overlappingSurface
+ (
+ *this,
+ procBb_[proci],
+ pointMap, // not used
+ faceMap
+ );
+ if (!faceMap.empty())
+ {
+ const triSurface subSurface
+ (
+ subsetMesh
+ (
+ *this,
+ faceMap,
+ pointMap
+ )
+ );
+ UOPstream os(proci, pBufs);
+ os << subSurface;
+ }
+ }
+ }
+
+ pBufs.finishedSends();
+
+ // Receive
+ for (const int proci : pBufs.allProcs())
+ {
+ if (proci != Pstream::myProcNo() && pBufs.recvDataCount(proci))
+ {
+ UIPstream is(proci, pBufs);
+
+ // Receive
+ const triSurface subSurface(is);
+
+ //Pout<< "** from proc:" << proci
+ // << " received:" << subSurface.size()
+ // << " compared to local:" << this->size() << endl;
+ //const fileName fName
+ //(
+ // searchableSurface::time().path()
+ // /searchableSurface::name()
+ // + "_from_proc" + Foam::name(proci) + ".obj"
+ //);
+ //Pout<< "fName:" << fName << endl;
+ //subSurface.write(fName);
+
+
+ // Set all tree nodes overlapping buffer triangles to
+ // mixed. This state will be preserved when caching the
+ // volume type (see below). This uses the actual local
+ // triangles.
+ for (const auto& tri : subSurface)
+ {
+ markMixed
+ (
+ tree,
+ 0, //nodei,
+ tri.tri(subSurface.points()),
+ nodeTypes
+ );
+ }
+ }
+ }
+}
+
+
void Foam::distributedTriSurfaceMesh::collectLeafMids
(
const label nodeI,
@@ -1616,10 +1968,100 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::calcVolumeType
}
else if (subType != myType)
{
- myType = volumeType::MIXED;
+ myType = volumeType::MIXED;
+ }
+ }
+ return myType;
+}
+
+
+void Foam::distributedTriSurfaceMesh::cacheVolumeType(PackedList<2>& nt) const
+{
+ // Get local tree
+ const indexedOctree<treeDataTriSurface>& t = tree();
+ const auto& nodes = t.nodes();
+
+ // Collect midpoints
+ DynamicField<point> midPoints(label(0.5*nodes.size()));
+ if (nodes.size())
+ {
+ collectLeafMids(0, midPoints);
+ }
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::cacheVolumeType :"
+ << " surface " << searchableSurface::name()
+ << " triggering orientation caching for "
+ << midPoints.size() << " leaf mids" << endl;
+ }
+
+ // Get volume type of mid points
+ List<volumeType> midVolTypes;
+
+ // Note: could recurse here (since now outsideVolType_ is set)
+ // but this would use the cached mid point data which we're trying
+ // to calculate. Instead bypass caching and do a full search
+ {
+ List<pointIndexHit> nearestInfo;
+ findNearest
+ (
+ midPoints,
+ scalarField(midPoints.size(), Foam::sqr(GREAT)),
+ nearestInfo
+ );
+ surfaceSide(midPoints, nearestInfo, midVolTypes);
+
+ if (debug & 2)
+ {
+ OBJstream insideStr
+ (
+ searchableSurface::time().path()
+ /searchableSurface::name() + "_cacheVolumeType_inside.obj"
+ );
+ OBJstream outsideStr
+ (
+ searchableSurface::time().path()
+ /searchableSurface::name() + "_cacheVolumeType_outside.obj"
+ );
+ forAll(midVolTypes, i)
+ {
+ const linePointRef ln(midPoints[i], nearestInfo[i].hitPoint());
+ if (midVolTypes[i] == volumeType::INSIDE)
+ {
+ insideStr.write(ln);
+ }
+ else if (midVolTypes[i] == volumeType::OUTSIDE)
+ {
+ outsideStr.write(ln);
+ }
+ }
+ Pout<< "Whilst caching " << searchableSurface::name()
+ << " have inside:" << insideStr.nVertices()
+ << " have outside:" << outsideStr.nVertices()
+ << endl;
}
}
- return myType;
+
+ // Cache on local tree
+ if (nodes.size())
+ {
+ label index = 0;
+ calcVolumeType
+ (
+ midVolTypes,
+ index,
+ nt,
+ 0 // nodeI
+ );
+ }
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::cacheVolumeType :"
+ << " surface " << searchableSurface::name()
+ << " done orientation caching for "
+ << midPoints.size() << " leaf mids" << endl;
+ }
}
@@ -1698,72 +2140,106 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::cachedVolumeType
}
-// 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
-)
+const Foam::decompositionMethod&
+Foam::distributedTriSurfaceMesh::decomposer() const
{
- addProfiling
- (
- distribute,
- "distributedTriSurfaceMesh::independentlyDistributedBbs"
- );
-
if (!decomposer_)
{
- // 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)
+ // Check if local dictionary contains decomposition entries
+ if (dict_.found("method"))
{
- decomposer_ = decompositionMethod::New(*dictPtr);
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::decomposer() :"
+ << " constructing decomposer from provided dictionary "
+ << dict_ << endl;
+ }
+
+ decomposer_ = decompositionMethod::New(dict_);
}
else
{
- if (!decomposeParDict_)
- {
- decomposeParDict_.reset
+ // Use singleton decomposeParDict. Cannot use decompositionModel
+ // here since we've only got Time and not a mesh.
+
+ const auto* dictPtr =
+ searchableSurface::time().findObject<IOdictionary>
(
- new IOdictionary
+ // == decompositionModel::canonicalName
+ "decomposeParDict"
+ );
+
+ if (dictPtr)
+ {
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::decomposer() :"
+ << " constructing decomposer from registered"
+ << " dictionary " << *dictPtr << endl;
+ }
+
+ decomposer_ = decompositionMethod::New(*dictPtr);
+ }
+ else
+ {
+ if (!decomposeParDict_)
+ {
+ decomposeParDict_.reset
(
- IOobject
+ new IOdictionary
(
- // == decompositionModel::canonicalName
- "decomposeParDict",
- searchableSurface::time().system(),
- searchableSurface::time(),
- IOobject::MUST_READ,
- IOobject::NO_WRITE
+ IOobject
+ (
+ // == decompositionModel::canonicalName
+ "decomposeParDict",
+ searchableSurface::time().system(),
+ searchableSurface::time(),
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
)
- )
- );
+ );
+
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::decomposer() :"
+ << " reading "
+ << decomposeParDict_().objectRelPath()
+ << endl;
+ }
+ }
+ decomposer_ = decompositionMethod::New(*decomposeParDict_);
}
- decomposer_ = decompositionMethod::New(*decomposeParDict_);
}
}
+ return decomposer_();
+}
+
- // Initialise to inverted box
- List<List<treeBoundBox>> bbs
+// 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.
+void Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
+(
+ const triSurface& s,
+ labelList& distribution, // per local triangle the destination processor
+ List<List<treeBoundBox>>& bbs
+) const
+{
+ addProfiling
(
- Pstream::nProcs(),
- List<treeBoundBox>(1, treeBoundBox::null())
+ distribute,
+ "distributedTriSurfaceMesh::independentlyDistributedBbs"
);
+ // Initialise to inverted box
+ bbs.resize_nocopy(Pstream::nProcs());
+ bbs = List<treeBoundBox>(1, treeBoundBox::null());
+
const globalIndex& triIndexer = globalTris();
bool masterOnly;
@@ -1795,8 +2271,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
triCentres[trii] = s[trii].centre(s.points());
}
- labelList distribution;
- if (!isA<geomDecomp>(decomposer_()))
+ if (!isA<geomDecomp>(decomposer()))
{
// Use connectivity
labelListList pointFaces;
@@ -1806,19 +2281,19 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
// Do the actual decomposition
const bool oldParRun = UPstream::parRun(false);
- distribution = decomposer_().decompose(faceFaces, triCentres);
+ distribution = decomposer().decompose(faceFaces, triCentres);
UPstream::parRun(oldParRun); // Restore parallel state
}
else
{
// Do the actual decomposition
- distribution = decomposer_().decompose(triCentres);
+ distribution = decomposer().decompose(triCentres);
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
- << " determining processor bounding boxes for surface"
+ << " determining processor bounding boxes for surface "
<< searchableSurface::name() << endl;
}
@@ -1845,10 +2320,10 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
{
// Fully distributed decomposition. Does not filter duplicate
// triangles.
- if (!decomposer_().parallelAware())
+ if (!decomposer().parallelAware())
{
FatalErrorInFunction
- << "The decomposition method " << decomposer_().typeName
+ << "The decomposition method " << decomposer().typeName
<< " does not decompose in parallel."
<< " Please choose one that does." << exit(FatalError);
}
@@ -1868,7 +2343,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
triCentres[trii] = s[trii].centre(s.points());
}
- labelList distribution = decomposer_().decompose(triCentres);
+ distribution = decomposer().decompose(triCentres);
if (debug)
{
@@ -1901,10 +2376,10 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
// {
// // Fully distributed decomposition. Does not filter duplicate
// // triangles.
-// if (!decomposer_().parallelAware())
+// if (!decomposer().parallelAware())
// {
// FatalErrorInFunction
-// << "The decomposition method " << decomposer_().typeName
+// << "The decomposition method " << decomposer().typeName
// << " does not decompose in parallel."
// << " Please choose one that does." << exit(FatalError);
// }
@@ -2060,7 +2535,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
// pointField masterCentres(allCentres, triMap);
//
// // Do the actual decomposition
-// labelList masterDistribution(decomposer_().decompose(masterCentres));
+// labelList masterDistribution(decomposer().decompose(masterCentres));
//
// // Make map to get the decomposition from the master of each border
// labelList borderGlobalMaster(borderTris.size());
@@ -2120,7 +2595,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
}
// Determine local decomposition
- labelList distribution(s.size());
{
labelList allDistribution;
if (Pstream::master())
@@ -2149,7 +2623,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
// Decompose merged centres
const bool oldParRun = UPstream::parRun(false);
- labelList mergedDist(decomposer_().decompose(mergedPoints));
+ labelList mergedDist(decomposer().decompose(mergedPoints));
UPstream::parRun(oldParRun); // Restore parallel state
// Convert to all
@@ -2202,7 +2676,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
}
Pstream::broadcast(bbs);
}
- return bbs;
}
@@ -2576,8 +3049,10 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
),
currentDistType_(FROZEN) // only used to trigger re-distribution
{
- // Read from the provided dictionary
- read();
+ // Construct from components. If empty dictionary get all the info from the
+ // supplied (undecomposed) triSurface.
+
+ readSettings(dict_.empty());
bounds().reduce();
@@ -2639,21 +3114,34 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
),
currentDistType_(FROZEN) // only used to trigger re-distribution
{
- // Read from the local, decomposed dictionary
- read();
-
- bounds().reduce();
-
// Try and find out where the triSurface was loaded from. On slave this
// might give empty filename.
const fileName actualFile(triSurfaceMesh::findFile(io, true));
+ // Was it read undecomposed?
+ const bool readFromMaster =
+ (
+ actualFile.empty()
+ || actualFile != io.localFilePath(triSurfaceMesh::typeName)
+ );
+
+ readSettings(readFromMaster);
+
+ bounds().reduce();
+
+ if (readFromMaster && !decomposeUsingBbs_)
+ {
+ // No fill-in so store normals
+ DebugInFunction
+ << "Determining vertex based normals since undecomposed" << endl;
+ const triSurface& surf = *this;
+ vertexNormals_ = autoPtr<List<FixedList<vector, 3>>>::New();
+ calcVertexNormals(surf, vertexNormals_());
+ }
+
if
(
- (
- actualFile.empty()
- || actualFile != io.localFilePath(triSurfaceMesh::typeName)
- )
+ readFromMaster
&& (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
)
{
@@ -2748,50 +3236,37 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
),
currentDistType_(FROZEN) // only used to trigger re-distribution
{
- // Read from the local, decomposed dictionary
- read();
-
- // Optionally override settings from provided dictionary
- {
- // Wanted distribution type
- distributionTypeNames_.readIfPresent
- (
- "distributionType",
- dict,
- distType_
- );
-
- // Merge distance
- dict.readIfPresent("mergeDistance", mergeDist_);
+ // Try and find out where the triSurface was loaded from. On slave this
+ // might give empty filename.
+ const fileName actualFile(triSurfaceMesh::findFile(io, dict, true));
- // Distribution type
- bool closed;
- if (dict.readIfPresent<bool>("closed", closed))
- {
- surfaceClosed_ = closed;
- }
+ // Was it read undecomposed?
+ const bool readFromMaster =
+ (
+ actualFile.empty()
+ || actualFile != io.localFilePath(triSurfaceMesh::typeName)
+ );
- outsideVolType_ = volumeType::names.getOrDefault
- (
- "outsideVolumeType",
- dict,
- outsideVolType_
- );
- }
+ // Add supplied dictionary entries (override existing)
+ dict_ <<= dict;
+ readSettings(readFromMaster);
bounds().reduce();
- // Try and find out where the triSurface was loaded from. On slave this
- // might give empty filename.
- const fileName actualFile(triSurfaceMesh::findFile(io, dict, true));
+ if (readFromMaster && !decomposeUsingBbs_)
+ {
+ // No fill-in so store enough normals to calculate
+ DebugInFunction
+ << "Determining vertex based normals since undecomposed" << endl;
+ const triSurface& surf = *this;
+ vertexNormals_ = autoPtr<List<FixedList<vector, 3>>>::New();
+ calcVertexNormals(surf, vertexNormals_());
+ }
if
(
- (
- actualFile.empty()
- || actualFile != io.localFilePath(triSurfaceMesh::typeName)
- )
+ readFromMaster
&& (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
)
{
@@ -2878,6 +3353,35 @@ const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
}
+void Foam::distributedTriSurfaceMesh::flip()
+{
+ if (debug)
+ {
+ Pout<< "distributedTriSurfaceMesh::flip :"
+ << " surface " << searchableSurface::name()
+ << " with current orientation "
+ << volumeType::names[outsideVolType_]
+ << " vertexNormals_ " << vertexNormals_.good()
+ << endl;
+ }
+
+ triSurfaceMesh::flip();
+
+ if (vertexNormals_)
+ {
+ List<FixedList<vector, 3>>& vn = vertexNormals_();
+
+ for (auto& normals : vn)
+ {
+ for (auto& n : normals)
+ {
+ n = -n;
+ }
+ }
+ }
+}
+
+
//void Foam::distributedTriSurfaceMesh::findNearest
//(
// const pointField& samples,
@@ -3902,7 +4406,7 @@ void Foam::distributedTriSurfaceMesh::getRegion
// Do my tests
// ~~~~~~~~~~~
- const triSurface& s = static_cast<const triSurface&>(*this);
+ const triSurface& s = *this;
region.setSize(triangleIndex.size());
@@ -3969,14 +4473,36 @@ void Foam::distributedTriSurfaceMesh::getNormal
// Do my tests
// ~~~~~~~~~~~
- const triSurface& s = static_cast<const triSurface&>(*this);
+ const triSurface& s = *this;
normal.setSize(triangleIndex.size());
- forAll(triangleIndex, i)
+ if (vertexNormals_)
{
- label trii = triangleIndex[i];
- normal[i] = s[trii].unitNormal(s.points());
+ // Use smooth interpolation for the normal
+
+ // Send over nearest point since we have it anyway
+ List<pointIndexHit> localNearestInfo(info);
+ map.distribute(localNearestInfo);
+
+ forAll(triangleIndex, i)
+ {
+ const point& nearest = localNearestInfo[i].point();
+
+ const label trii = triangleIndex[i];
+ const triPointRef tri(s[trii].tri(s.points()));
+ const auto& vn = vertexNormals_()[trii];
+ const barycentric2D w(tri.pointToBarycentric(nearest));
+ normal[i] = w[0]*vn[0]+w[1]*vn[1]+w[2]*vn[2];
+ }
+ }
+ else
+ {
+ forAll(triangleIndex, i)
+ {
+ const label trii = triangleIndex[i];
+ normal[i] = s[trii].unitNormal(s.points());
+ }
}
@@ -4138,64 +4664,18 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
|| outsideVolType_ == volumeType::OUTSIDE
)
{
- // Get local tree
const indexedOctree<treeDataTriSurface>& t = tree();
const auto& nodes = t.nodes();
PackedList<2>& nt = t.nodeTypes();
nt.resize(nodes.size());
nt = volumeType::UNKNOWN;
- // Collect midpoints
- DynamicField<point> midPoints(label(0.5*nodes.size()));
- if (nodes.size())
- {
- collectLeafMids(0, midPoints);
- }
-
- if (debug)
- {
- Pout<< "distributedTriSurfaceMesh::getVolumeType :"
- << " surface " << searchableSurface::name()
- << " triggering orientation caching for "
- << midPoints.size() << " leaf mids" << endl;
- }
-
- // Get volume type of mid points
- List<volumeType> midVolTypes;
-
- // Note: could recurse here (since now outsideVolType_ is set)
- // but this would use the cached mid point data which we're trying
- // to calculate. Instead bypass caching and do a full search
+ if (!decomposeUsingBbs_)
{
- List<pointIndexHit> nearestInfo;
- findNearest
- (
- midPoints,
- scalarField(midPoints.size(), Foam::sqr(GREAT)),
- nearestInfo
- );
- surfaceSide(midPoints, nearestInfo, midVolTypes);
+ markMixedOverlap(t, nt);
}
- // Cache on local tree
- if (nodes.size())
- {
- label index = 0;
- calcVolumeType
- (
- midVolTypes,
- index,
- nt,
- 0 // nodeI
- );
- }
- if (debug)
- {
- Pout<< "distributedTriSurfaceMesh::getVolumeType :"
- << " surface " << searchableSurface::name()
- << " done orientation caching for "
- << midPoints.size() << " leaf mids" << endl;
- }
+ cacheVolumeType(nt);
}
}
@@ -4282,9 +4762,18 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
// those I need to search the nearest for
DynamicField<point> fullSearchPoints(localPoints.size());
DynamicList<label> fullSearchMap(localPoints.size());
+
forAll(localPoints, i)
{
- if (tree().nodes().size())
+ // If no-fill in:
+ // - tree nodes on the outside will not include fill-in triangles
+ // - so might decide the wrong state (might be 'mixed' with fill-in
+ // triangles)
+ // - for now do not cache
+ // - TBD. Note that the cached volume types already include remote
+ // fill-in triangles. There is still a problem there though. TBD.
+
+ if (decomposeUsingBbs_ && tree().nodes().size())
{
volType[i] = cachedVolumeType(0, localPoints[i]);
}
@@ -4322,9 +4811,68 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
// Combine
forAll(fullSearchMap, i)
{
- volType[fullSearchMap[i]] = fullSearchType[i];
+ volumeType& vt = volType[fullSearchMap[i]];
+ if (vt == volumeType::UNKNOWN)
+ {
+ vt = fullSearchType[i];
+ }
+ else if (vt != fullSearchType[i])
+ {
+ Pout<< "At point:" << fullSearchPoints[i]
+ << " or point:" << localPoints[fullSearchMap[i]]
+ << " have cached status:" << vt
+ << " have full-search status:" << fullSearchType[i]
+ << endl;
+ }
+ //volType[fullSearchMap[i]] = fullSearchType[i];
}
+ //{
+ // //- Variant of below that produces better error messages
+ // typedef Tuple2<Pair<point>, volumeType> NearType;
+ // const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
+ //
+ // // Combine nearest and volType so we can give nicer error messages
+ // List<NearType> nearTypes(volType.size(), zero);
+ //
+ // NearTypeCombineOp cop;
+ //
+ // // For all volType we do have the originating point but not the
+ // // nearest
+ // forAll(localPoints, i)
+ // {
+ // nearTypes[i] =
+ // NearType(Pair<point>(localPoints[i], Zero), volType[i]);
+ // }
+ // // But for all full-search we also have the nearest
+ // forAll(fullSearchMap, i)
+ // {
+ // const NearType nt
+ // (
+ // Pair<point>(fullSearchPoints[i], nearestInfo[i].hitPoint()),
+ // fullSearchType[i]
+ // );
+ //
+ // cop(nearTypes[fullSearchMap[i]], nt);
+ // }
+ // mapDistributeBase::distribute
+ // (
+ // Pstream::commsTypes::nonBlocking,
+ // List<labelPair>::null(),
+ // samples.size(),
+ // map.constructMap(),
+ // map.constructHasFlip(),
+ // map.subMap(),
+ // map.subHasFlip(),
+ // nearTypes,
+ // zero,
+ // cop, //NearTypeCombineOp(),
+ // noOp(), // no flipping
+ // UPstream::msgType(),
+ // map.comm()
+ // );
+ //}
+
// Send back to originator. In case of multiple answers choose inside or
// outside
const volumeType zero(volumeType::UNKNOWN);
@@ -4345,47 +4893,6 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
map.comm()
);
-
- //// Combine nearest and volType so we can give nicer error messages
- //List<NearType> nearTypes(volType.size());
- //// For all volType we do have the originating point but not the nearest
- //forAll(localPoints, i)
- //{
- // nearTypes[i] =
- // NearType(Pair<point>(localPoints[i], Zero), volType[i]);
- //}
- //// But for all full-search we also have the nearest
- //forAll(fullSearchMap, i)
- //{
- // nearTypes[fullSearchMap[i]] = NearType
- // (
- // Pair<point>(fullSearchPoints[i], nearestInfo[i].hitPoint()),
- // fullSearchType[i]
- // );
- //}
- //const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
- //mapDistributeBase::distribute
- //(
- // Pstream::commsTypes::nonBlocking,
- // List<labelPair>::null(),
- // samples.size(),
- // map.constructMap(),
- // map.constructHasFlip(),
- // map.subMap(),
- // map.subHasFlip(),
- // nearTypes,
- // zero,
- // NearTypeCombineOp(),
- // noOp(), // no flipping
- // UPstream::msgType(),
- // map.comm()
- //);
- //volType.setSize(nearTypes.size());
- //forAll(nearTypes, i)
- //{
- // volType[i] = nearTypes[i].second();
- //}
-
// Add the points outside the bounding box
for (label samplei : outsideSamples)
{
@@ -4613,6 +5120,7 @@ void Foam::distributedTriSurfaceMesh::distribute
<< " surface " << searchableSurface::name()
<< " distributing surface according to method:"
<< distributionTypeNames_[distType_]
+ << " fill-in:" << decomposeUsingBbs_
<< " follow bbs:" << flatOutput(bbs) << endl;
}
@@ -4622,7 +5130,10 @@ void Foam::distributedTriSurfaceMesh::distribute
// Get bbs of all domains
// ~~~~~~~~~~~~~~~~~~~~~~
+ // Per triangle the destination processor
+ labelList distribution;
{
+ // Per processor the bounding box of all (destination) triangles
List<List<treeBoundBox>> newProcBb(Pstream::nProcs());
switch (distType_)
@@ -4638,7 +5149,7 @@ void Foam::distributedTriSurfaceMesh::distribute
{
return;
}
- newProcBb = independentlyDistributedBbs(*this);
+ independentlyDistributedBbs(*this, distribution, newProcBb);
break;
case FROZEN:
@@ -4691,15 +5202,37 @@ void Foam::distributedTriSurfaceMesh::distribute
labelListList faceSendMap(Pstream::nProcs());
labelListList pointSendMap(Pstream::nProcs());
- forAll(procBb_, proci)
+ if (decomposeUsingBbs_)
{
- overlappingSurface
- (
- *this,
- procBb_[proci],
- pointSendMap[proci],
- faceSendMap[proci]
- );
+ forAll(procBb_, proci)
+ {
+ overlappingSurface
+ (
+ *this,
+ procBb_[proci],
+ pointSendMap[proci],
+ faceSendMap[proci]
+ );
+ }
+ }
+ else
+ {
+ const triSurface& s = *this;
+ forAll(procBb_, proci)
+ {
+ boolList includedFace(s.size(), false);
+ forAll(s, trii)
+ {
+ includedFace[trii] = (distribution[trii] == proci);
+ }
+ subsetMesh
+ (
+ s,
+ includedFace,
+ pointSendMap[proci],
+ faceSendMap[proci]
+ );
+ }
}
if (keepNonLocal)
@@ -4707,7 +5240,7 @@ void Foam::distributedTriSurfaceMesh::distribute
// Include in faceSendMap/pointSendMap the triangles that are
// not mapped to any processor so they stay local.
- const triSurface& s = static_cast<const triSurface&>(*this);
+ const triSurface& s = *this;
boolList includedFace(s.size(), true);
@@ -4819,7 +5352,7 @@ void Foam::distributedTriSurfaceMesh::distribute
UIPstream is(proci, pBufs);
// Receive
- triSurface subSurface(is);
+ const triSurface subSurface(is);
// Merge into allSurf
merge
@@ -4876,6 +5409,63 @@ void Foam::distributedTriSurfaceMesh::distribute
distributeFields<sphericalTensor>(faceMap());
distributeFields<symmTensor>(faceMap());
distributeFields<tensor>(faceMap());
+ if (vertexNormals_)
+ {
+ List<FixedList<vector, 3>>& vn = vertexNormals_();
+ faceMap().distribute(vn);
+
+ if (debug & 2)
+ {
+ OBJstream str
+ (
+ searchableSurface::time().path()
+ /searchableSurface::name() + "_vertex_normals.obj"
+ );
+ const triSurface& s = *this;
+ forAll(s, trii)
+ {
+ const auto& f = s[trii];
+ const scalar l = Foam::sqrt(f.mag(s.points()));
+ forAll(f, fp)
+ {
+ const point& pt = s.points()[f[fp]];
+ const vector& n = vn[trii][fp];
+ str.write(linePointRef(pt, pt+l*n));
+ }
+ }
+ Pout<< "Dumped local vertex normals to " << str.name() << endl;
+ }
+ }
+
+
+ if
+ (
+ outsideVolType_ == volumeType::INSIDE
+ || outsideVolType_ == volumeType::OUTSIDE
+ )
+ {
+ // Re-build tree & inside/outside for closed surfaces
+
+ const indexedOctree<treeDataTriSurface>& t = tree();
+ const auto& nodes = t.nodes();
+ PackedList<2>& nt = t.nodeTypes();
+ nt.resize(nodes.size());
+ nt = volumeType::UNKNOWN;
+
+ // Send over any overlapping but not included triangles
+ if (!decomposeUsingBbs_)
+ {
+ // We do not backfill the bounding box so the local tree might
+ // not have all the triangles to cache inside/outside correctly.
+ // Do a pre-pass to mark 'mixed' the tree leaves that intersect
+ // the triangles of the (would be) fill-in triangles.
+ markMixedOverlap(t, nt);
+ }
+
+ cacheVolumeType(nt);
+ }
+
+
if (debug)
{
@@ -4995,16 +5585,17 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
{
boundBox bb;
label nPoints;
- PatchTools::calcBounds(static_cast<const triSurface&>(*this), bb, nPoints);
+ PatchTools::calcBounds(*this, bb, nPoints);
bb.reduce();
- os << "Triangles : " << returnReduce(triSurface::size(), sumOp<label>())
- << endl
- << "Vertices : " << returnReduce(nPoints, sumOp<label>()) << endl
- << "Bounding Box : " << bb << endl
- << "Closed : " << surfaceClosed_ << endl
- << "Outside point: " << volumeType::names[outsideVolType_] << endl
- << "Distribution : " << distributionTypeNames_[distType_] << endl;
+ os << "Triangles : "
+ << returnReduce(triSurface::size(), sumOp<label>()) << endl
+ << "Vertices : " << returnReduce(nPoints, sumOp<label>()) << endl
+ << "Vertex normals : " << vertexNormals_.valid() << endl
+ << "Bounding Box : " << bb << endl
+ << "Closed : " << surfaceClosed_ << endl
+ << "Outside type : " << 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 05adb82d5c8d08746b7a500b339c410b6cfa4af4..5c4498994b535a82275a2d6dfb16d35242db3f41 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2015-2022 OpenCFD Ltd.
+ Copyright (C) 2015-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -109,20 +109,28 @@ private:
//- Merging distance
scalar mergeDist_;
- autoPtr<IOdictionary> decomposeParDict_;
+ mutable autoPtr<IOdictionary> decomposeParDict_;
//- Decomposition used when independently decomposing surface.
- autoPtr<decompositionMethod> decomposer_;
+ mutable autoPtr<decompositionMethod> decomposer_;
//- Bounding box settings
localIOdictionary dict_;
+ //- Use bounding boxes (default) or unique decomposition of triangles
+ //- (i.e. do not duplicate triangles)
+ bool decomposeUsingBbs_;
+
//- Bounding boxes of all processors
List<List<treeBoundBox>> procBb_;
//- Global triangle numbering
mutable autoPtr<globalIndex> globalTris_;
+ //- Optional per-vertex normals. TBD: move to triSurface? or have
+ //- per-triangle 3 normals so we can interpolate and have features
+ mutable autoPtr<List<FixedList<vector, 3>>> vertexNormals_;
+
//- The (wanted) distribution type.
distributionType distType_;
@@ -139,8 +147,17 @@ private:
// directory or in parent directory
static word findLocalInstance(const IOobject& io);
- //- Read my additional data
- bool read();
+ //- Read my additional data from dictionary. Additional flag to
+ //- say whether we can use master-only geometric tests.
+ bool readSettings(const bool isUndecomposed);
+
+ //- Construction helper: generate vertex normals upon reading
+ //- undecomposed surface
+ void calcVertexNormals
+ (
+ const triSurface& surf,
+ List<FixedList<vector, 3>>& vn
+ ) const;
// Line intersection
@@ -287,6 +304,23 @@ private:
// Caching of volume type (based on indexedOctree)
+ //- Set node type on any node containing the triangle
+ volumeType markMixed
+ (
+ const indexedOctree<treeDataTriSurface>& tree,
+ const label nodei,
+ const triPointRef& tri,
+ PackedList<2>& nodeTypes
+ ) const;
+
+ //- Set node type on any node overlapping any remote triangles.
+ //- Only valid if using unique decomposition.
+ void markMixedOverlap
+ (
+ const indexedOctree<treeDataTriSurface>& tree,
+ PackedList<2>& nodeTypes
+ ) const;
+
//- Collect mid points of tree boxes
void collectLeafMids
(
@@ -303,6 +337,9 @@ private:
const label nodeI
) const;
+ //- Calculate inside/outside of midpoint of tree nodes
+ void cacheVolumeType(PackedList<2>& nt) const;
+
//- Look up any cached data. Return unknown if cannot be determined.
volumeType cachedVolumeType
(
@@ -321,11 +358,16 @@ private:
labelListList& faceFaces
);
+ //- Helper: get decompositionMethod
+ const decompositionMethod& decomposer() const;
+
//- Finds new bounds based on an independent decomposition.
- List<List<treeBoundBox>> independentlyDistributedBbs
+ void independentlyDistributedBbs
(
- const triSurface&
- );
+ const triSurface& s,
+ labelList& distribution,
+ List<List<treeBoundBox>>& bbs
+ ) const;
//- Does any part of triangle overlap bb.
static bool overlaps
@@ -413,7 +455,8 @@ public:
// Constructors
- //- Construct from triSurface
+ //- Construct from components. Assumes triSurface is already decomposed
+ //- and dictionary contains corresponding information
distributedTriSurfaceMesh
(
const IOobject&,
@@ -421,11 +464,17 @@ public:
const dictionary& dict
);
- //- Construct read. Does findInstance to find io.local().
+ //- Construct read. Does findInstance to find io.local()
+ //- - if found local : assume distributed
+ //- - if found in parent : assume undistributed. Can e.g. check for
+ //- closedness.
distributedTriSurfaceMesh(const IOobject& io);
//- Construct from dictionary (used by searchableSurface).
// Does read. Does findInstance to find io.local().
+ //- - if found local : assume distributed
+ //- - if found in parent : assume undistributed. Can e.g. check for
+ //- closedness.
distributedTriSurfaceMesh
(
const IOobject& io,
@@ -454,6 +503,9 @@ public:
return globalTris().totalSize();
}
+ //- Flip triangles, outsideVolumeType and all cached inside/outside.
+ virtual void flip();
+
virtual void findNearest
(
const pointField& sample,
diff --git a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
index 73fe2e15fe4d2750c286255c0922fe4d31ff62eb..2fce36be064a19bbf7782dcab39539732edacee0 100644
--- a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
+++ b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2018-2022 OpenCFD Ltd.
+ Copyright (C) 2018-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -61,10 +61,9 @@ Foam::patchDistMethods::exact::patchSurface() const
localBb.extend(rndGen, 1E-3)
);
- // Dummy bounds dictionary
- dictionary dict;
- dict.add("bounds", meshBb);
- dict.add
+ // Add any missing properties (but not override existing ones)
+ dict_.add("bounds", meshBb);
+ dict_.add
(
"distributionType",
distributedTriSurfaceMesh::distributionTypeNames_
@@ -74,8 +73,7 @@ Foam::patchDistMethods::exact::patchSurface() const
distributedTriSurfaceMesh::DISTRIBUTED // parallel decomp
]
);
- dict.add("mergeDistance", 1e-6*localBb.mag());
-
+ dict_.add("mergeDistance", 1e-6*localBb.mag());
Info<< "Triangulating local patch faces" << nl << endl;
@@ -100,7 +98,7 @@ Foam::patchDistMethods::exact::patchSurface() const
patchIDs_,
mapTriToGlobal
),
- dict
+ dict_
)
);
@@ -131,7 +129,8 @@ Foam::patchDistMethods::exact::exact
const labelHashSet& patchIDs
)
:
- patchDistMethod(mesh, patchIDs)
+ patchDistMethod(mesh, patchIDs),
+ dict_(dict.subOrEmptyDict(typeName + "Coeffs"))
{}
diff --git a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
index 4fd1f82fb1ac3a1fc3f18a3f56a4565b411a3cfb..3ac84561ce1cbeb27d14b4f14f234bbe5db098a2 100644
--- a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
+++ b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2018-2019 OpenCFD Ltd.
+ Copyright (C) 2018-2019,2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -30,9 +30,32 @@ Description
Calculation of exact distance to nearest patch for all cells and
boundary by constructing a search tree for all patch faces.
+Usage
+
+ \begin{verbatim}
+ wallDist
+ {
+ method exactDistance;
+
+ // Optional entries (currently for distributedTriSurfaceMesh)
+ exactDistanceCoeffs
+ {
+ // Optional decomposition method. If not supplied will use the
+ // system/decomposeParDict.
+ method hierarchical;
+ numberOfSubdomains 8;
+ n (2 2 2);
+
+ // (not)add fill-in triangles
+ decomposeUsingBbs false;
+ }
+ }
+ \end{verbatim}
+
See also
Foam::patchDistMethod::meshWave
Foam::wallDist
+ Foam::distributedTriSurfaceMesh
SourceFiles
exactPatchDistMethod.C
@@ -55,7 +78,7 @@ namespace patchDistMethods
{
/*---------------------------------------------------------------------------*\
- Class exact Declaration
+ Class exact Declaration
\*---------------------------------------------------------------------------*/
class exact
@@ -64,6 +87,9 @@ class exact
{
// Private Member Data
+ //- Dictionary contents from the dictionary constructor
+ mutable dictionary dict_;
+
//- Cache surface+searching of patch
mutable autoPtr<distributedTriSurfaceMesh> patchSurfPtr_;
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
index 9fb6072860498dd34dfa7daea8d5eeedab832bfe..bb4f915e9e8c3ffe5e49756ed7736c20786dcfb4 100644
--- a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
@@ -20,6 +20,18 @@ snap true;
addLayers false;
+// Define common distributedTriSurface properties
+_decomp
+{
+ numberOfSubdomains 8;
+ method hierarchical;
+ n (2 2 2);
+
+ // unique triangle decomposition; no fill-in
+ decomposeUsingBbs false;
+}
+
+
// Geometry. Definition of all surfaces. All surfaces are of class
// searchableSurface.
// Surfaces are used
@@ -32,16 +44,19 @@ geometry
{
type distributedTriSurfaceMesh;
file "box.obj"; //"box_12_2.obj";
+ ${_decomp}
}
box_trans
{
- file "box_trans.obj";
type distributedTriSurfaceMesh;
+ file "box_trans.obj";
+ ${_decomp}
}
shell
{
- file "box_scaled.obj";
type distributedTriSurfaceMesh;
+ file "box_scaled.obj";
+ ${_decomp}
}
};