/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2015-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see 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;
}
}
};
}
const Foam::Enum
<
Foam::distributedTriSurfaceMesh::distributionType
>
Foam::distributedTriSurfaceMesh::distributionTypeNames_
({
{ distributionType::FOLLOW, "follow" },
{ distributionType::INDEPENDENT, "independent" },
{ 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());
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 outsideVolTypes;
// triSurfaceMesh::getVolumeType
// (
// pointField(1, outsidePt),
// outsideVolTypes
// );
// outsideVolType_ = outsideVolTypes[0];
//}
//dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
}
else
{
dict_.readEntry("bounds", procBb_[Pstream::myProcNo()]);
Pstream::gatherList(procBb_);
Pstream::scatterList(procBb_);
// Wanted distribution type
distType_ = distributionTypeNames_.get("distributionType", dict_);
// Merge distance
dict_.readEntry("mergeDistance", mergeDist_);
// Distribution type
surfaceClosed_ = dict_.getOrDefault("closed", false);
outsideVolType_ = volumeType::names.getOrDefault
(
"outsideVolumeType",
dict_,
volumeType::UNKNOWN
);
}
return true;
}
// Is segment fully local?
bool Foam::distributedTriSurfaceMesh::isLocal
(
const List& myBbs,
const point& start,
const point& end
)
{
forAll(myBbs, bbi)
{
if (myBbs[bbi].contains(start) && myBbs[bbi].contains(end))
{
return true;
}
}
return false;
}
//void Foam::distributedTriSurfaceMesh::splitSegment
//(
// const label segmenti,
// const point& start,
// const point& end,
// const treeBoundBox& bb,
//
// DynamicList& allSegments,
// DynamicList& allSegmentMap,
// DynamicList sendMap
//) const
//{
// // Work points
// point clipPt0, clipPt1;
//
// if (bb.contains(start))
// {
// // start within, trim end to bb
// bool clipped = bb.intersects(end, start, clipPt0);
//
// if (clipped)
// {
// // segment from start to clippedStart passes
// // through proc.
// sendMap[proci].append(allSegments.size());
// allSegmentMap.append(segmenti);
// allSegments.append(segment(start, clipPt0));
// }
// }
// else if (bb.contains(end))
// {
// // end within, trim start to bb
// bool clipped = bb.intersects(start, end, clipPt0);
//
// if (clipped)
// {
// sendMap[proci].append(allSegments.size());
// allSegmentMap.append(segmenti);
// allSegments.append(segment(clipPt0, end));
// }
// }
// else
// {
// // trim both
// bool clippedStart = bb.intersects(start, end, clipPt0);
//
// if (clippedStart)
// {
// bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
//
// if (clippedEnd)
// {
// // middle part of segment passes through proc.
// sendMap[proci].append(allSegments.size());
// allSegmentMap.append(segmenti);
// allSegments.append(segment(clipPt0, clipPt1));
// }
// }
// }
//}
void Foam::distributedTriSurfaceMesh::distributeSegment
(
const label segmenti,
const point& start,
const point& end,
DynamicList& allSegments,
DynamicList& allSegmentMap,
List>& sendMap
) const
{
// 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())
{
const List& bbs = procBb_[proci];
if (isLocal(bbs, start, end))
{
sendMap[proci].append(allSegments.size());
allSegmentMap.append(segmenti);
allSegments.append(segment(start, end));
return;
}
}
}
// 3. If not contained in single processor send to all intersecting
// processors.
forAll(procBb_, proci)
{
const List& bbs = procBb_[proci];
forAll(bbs, bbi)
{
const treeBoundBox& bb = bbs[bbi];
// Scheme a: any processor that intersects the segment gets
// the whole segment.
// Intersection point
point clipPt;
if (bb.intersects(start, end, clipPt))
{
sendMap[proci].append(allSegments.size());
allSegmentMap.append(segmenti);
allSegments.append(segment(start, end));
}
// Alternative: any processor only gets clipped bit of
// segment. This gives small problems with additional
// truncation errors.
//splitSegment
//(
// segmenti,
// start,
// end,
// bb,
//
// allSegments,
// allSegmentMap,
// sendMap[proci]
//);
}
}
}
Foam::autoPtr
Foam::distributedTriSurfaceMesh::distributeSegments
(
const pointField& start,
const pointField& end,
List& allSegments,
labelList& allSegmentMap
) const
{
// Determine send map
// ~~~~~~~~~~~~~~~~~~
labelListList sendMap(Pstream::nProcs());
{
// Since intersection test is quite expensive compared to memory
// allocation we use DynamicList to immediately store the segment
// in the correct bin.
// Segments to test
DynamicList dynAllSegments(start.size());
// Original index of segment
DynamicList dynAllSegmentMap(start.size());
// Per processor indices into allSegments to send
List> dynSendMap(Pstream::nProcs());
// Pre-size
forAll(dynSendMap, proci)
{
dynSendMap[proci].reserve
(
(proci == Pstream::myProcNo())
? start.size()
: start.size()/Pstream::nProcs()
);
}
forAll(start, segmenti)
{
distributeSegment
(
segmenti,
start[segmenti],
end[segmenti],
dynAllSegments,
dynAllSegmentMap,
dynSendMap
);
}
// Convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
dynSendMap[proci].shrink();
sendMap[proci].transfer(dynSendMap[proci]);
}
allSegments.transfer(dynAllSegments);
allSegmentMap.transfer(dynAllSegmentMap);
}
return autoPtr::New(std::move(sendMap));
}
void Foam::distributedTriSurfaceMesh::findLine
(
const bool nearestIntersection,
const pointField& start,
const pointField& end,
List& info
) const
{
if (debug)
{
Pout<< "distributedTriSurfaceMesh::findLine :"
<< " intersecting with "
<< start.size() << " rays" << endl;
}
addProfiling(findLine, "distributedTriSurfaceMesh::findLine");
const indexedOctree& octree = tree();
// Initialise
info.setSize(start.size());
forAll(info, i)
{
info[i].setMiss();
}
// 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 (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++;
}
}
if
(
returnReduce(nLocal, sumOp())
< returnReduce(start.size(), sumOp())
)
{
// Not all can be resolved locally. Build segments and map,
// send over segments, do intersections, send back and merge.
// Construct queries (segments)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Segments to test
List allSegments(start.size());
// Original index of segment
labelList allSegmentMap(start.size());
const autoPtr mapPtr
(
distributeSegments
(
start,
end,
allSegments,
allSegmentMap
)
);
const mapDistribute& map = mapPtr();
label nOldAllSegments = allSegments.size();
// Exchange the segments
// ~~~~~~~~~~~~~~~~~~~~~
map.distribute(allSegments);
// Do tests i need to do
// ~~~~~~~~~~~~~~~~~~~~~
// Intersections
List intersections(allSegments.size());
forAll(allSegments, i)
{
if (nearestIntersection)
{
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())
);
}
}
// Exchange the intersections (opposite to segments)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
map.reverseDistribute(nOldAllSegments, intersections);
// Extract the hits
// ~~~~~~~~~~~~~~~~
forAll(intersections, i)
{
const pointIndexHit& allInfo = intersections[i];
label segmenti = allSegmentMap[i];
pointIndexHit& hitInfo = info[segmenti];
if (allInfo.hit())
{
if (!hitInfo.hit())
{
// 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;
}
}
}
}
}
}
void Foam::distributedTriSurfaceMesh::convertTriIndices
(
List& 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
Foam::autoPtr
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
const List& info,
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();
// Determine send map
// ~~~~~~~~~~~~~~~~~~
// Since determining which processor the query should go to is
// cheap we do a multi-pass algorithm to save some memory temporarily.
// 1. Count
labelList nSend(Pstream::nProcs(), Zero);
forAll(info, i)
{
if (info[i].hit())
{
label proci = triIndexer.whichProcID(info[i].index());
nSend[proci]++;
}
}
// 2. Size sendMap
labelListList sendMap(Pstream::nProcs());
forAll(nSend, proci)
{
sendMap[proci].setSize(nSend[proci]);
nSend[proci] = 0;
}
// 3. Fill sendMap
forAll(info, i)
{
if (info[i].hit())
{
label proci = triIndexer.whichProcID(info[i].index());
triangleIndex[i] = triIndexer.toLocal(proci, info[i].index());
sendMap[proci][nSend[proci]++] = i;
}
else
{
triangleIndex[i] = -1;
}
}
// Pack into distribution map
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr mapPtr(new mapDistribute(std::move(sendMap)));
// Send over queries
// ~~~~~~~~~~~~~~~~~
mapPtr().distribute(triangleIndex);
return mapPtr;
}
bool Foam::distributedTriSurfaceMesh::contains
(
const List& bbs,
const point& sample
) const
{
forAll(bbs, bbi)
{
if (bbs[bbi].contains(sample))
{
return true;
}
}
return false;
}
Foam::Tuple2
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
procContains.setSize(Pstream::nProcs());
procContains = false;
procOverlaps.setSize(Pstream::nProcs());
procOverlaps = false;
minProci = -1;
scalar minDistSqr = radiusSqr;
label nContain = 0;
forAll(procBb_, proci)
{
const List& 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(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& 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& bbs = procBb_[proci];
forAll(bbs, bbi)
{
if (bbs[bbi].overlaps(centre, maxDistSqr))
{
procOverlaps[proci] = true;
nOverlap++;
break;
}
}
}
}
}
return Tuple2(nOverlap, maxDistSqr);
}
}
Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
(
const point& centre,
const scalar radiusSqr,
boolList& overlaps
) const
{
overlaps = false;
label nOverlaps = 0;
forAll(procBb_, proci)
{
const List& bbs = procBb_[proci];
forAll(bbs, bbi)
{
if (bbs[bbi].overlaps(centre, radiusSqr))
{
overlaps[proci] = true;
nOverlaps++;
break;
}
}
}
return nOverlaps;
}
// Generate queries for parallel distance calculation
// - calculates exchange map
// - uses map to exchange points and radius
Foam::autoPtr
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
const bool includeLocalProcessor,
const pointField& centres,
const scalarField& radiusSqr,
pointField& allCentres,
scalarField& allRadiusSqr,
labelList& allSegmentMap
) const
{
// Determine queries
// ~~~~~~~~~~~~~~~~~
labelListList sendMap(Pstream::nProcs());
{
// Queries
DynamicList dynAllCentres(centres.size());
DynamicList dynAllRadiusSqr(centres.size());
// Original index of segment
DynamicList dynAllSegmentMap(centres.size());
// Per processor indices into allSegments to send
List> 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());
forAll(centres, centrei)
{
// Find the processor this sample+radius overlaps.
calcOverlappingProcs
(
centres[centrei],
radiusSqr[centrei],
procBbOverlaps
);
forAll(procBbOverlaps, proci)
{
if
(
procBbOverlaps[proci]
&& (
includeLocalProcessor
|| proci != Pstream::myProcNo()
)
)
{
dynSendMap[proci].append(dynAllCentres.size());
dynAllSegmentMap.append(centrei);
dynAllCentres.append(centres[centrei]);
dynAllRadiusSqr.append(radiusSqr[centrei]);
}
}
}
// Convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, proci)
{
dynSendMap[proci].shrink();
sendMap[proci].transfer(dynSendMap[proci]);
}
allCentres.transfer(dynAllCentres);
allRadiusSqr.transfer(dynAllRadiusSqr);
allSegmentMap.transfer(dynAllSegmentMap);
}
return autoPtr::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(*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);
if (((sample - nearestPoint) & n) > 0)
{
return volumeType::OUTSIDE;
}
else
{
return volumeType::INSIDE;
}
}
Foam::label Foam::distributedTriSurfaceMesh::findOtherFace
(
const labelListList& pointFaces,
const label nearFacei,
const label nearLabel
) const
{
const triSurface& surf = static_cast(*this);
const triSurface::face_type& nearF = surf[nearFacei];
const edge e(nearF[nearLabel], nearF[nearF.fcIndex(nearLabel)]);
const labelList& pFaces = pointFaces[e[0]];
for (const label facei : pFaces)
{
if (facei != nearFacei)
{
int dir = surf[facei].edgeDirection(e);
if (dir != 0)
{
return facei;
}
}
}
return -1;
}
void Foam::distributedTriSurfaceMesh::calcFaceFaces
(
const triSurface& s,
const labelListList& pointFaces,
labelListList& faceFaces
)
{
faceFaces.setSize(s.size());
DynamicList nbrs;
forAll(faceFaces, facei)
{
const labelledTri& f = s[facei];
nbrs.reserve(f.size());
nbrs.clear();
forAll(f, fp)
{
const edge e(f[fp], f[f.fcIndex(fp)]);
const labelList& pFaces = pointFaces[f[fp]];
for (const label otherFacei : pFaces)
{
if (otherFacei != facei)
{
if (s[otherFacei].edgeDirection(e) != 0)
{
if (!nbrs.find(otherFacei))
{
nbrs.append(otherFacei);
}
}
}
}
}
faceFaces[facei] = std::move(nbrs);
}
}
void Foam::distributedTriSurfaceMesh::surfaceSide
(
const pointField& samples,
const List& nearestInfo,
List& volType
) const
{
if (debug)
{
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
labelList triangleIndex(nearestInfo.size());
autoPtr mapPtr
(
calcLocalQueries
(
nearestInfo,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Send over samples
pointField localSamples(samples);
map.distribute(localSamples);
// Do my tests
// ~~~~~~~~~~~
volType.setSize(triangleIndex.size());
volType = volumeType::UNKNOWN;
const triSurface& surf = static_cast(*this);
const pointField& points = surf.points();
{
//const labelListList& pointFaces = surf.pointFaces();
// Construct pointFaces. Let's hope surface has compact point
// numbering ...
labelListList pointFaces;
invertManyToMany(points.size(), surf, pointFaces);
EdgeMap edgeToFaces;
forAll(triangleIndex, i)
{
const label facei = triangleIndex[i];
const triSurface::face_type& f = surf[facei];
const point& sample = localSamples[i];
// Find where point is on face
label nearType, nearLabel;
pointHit pHit =
f.nearestPointClassify(sample, points, nearType, nearLabel);
const point& nearestPoint(pHit.rawPoint());
if (nearType == triPointRef::NONE)
{
const vector sampleNearestVec = (sample - nearestPoint);
// Nearest to face interior. Use faceNormal to determine side
//scalar c = sampleNearestVec & surf.faceNormals()[facei];
scalar c = sampleNearestVec & surf[facei].areaNormal(points);
if (c > 0)
{
volType[i] = volumeType::OUTSIDE;
}
else
{
volType[i] = volumeType::INSIDE;
}
}
else if (nearType == triPointRef::EDGE)
{
// 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
}
}
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);
// Loop over all faces. Check if have both edge faces. Do
// test.
edgeToFaces.clear();
scalar maxCosAngle = -GREAT;
labelPair maxEdgeFaces(-1, -1);
for (const label facei : pFaces)
{
const triSurface::face_type& f = surf[facei];
label fp = f.find(pointi);
label p1 = f[f.fcIndex(fp)];
label pMin1 = f[f.rcIndex(fp)];
Pair 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(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& midPoints
) const
{
const indexedOctree::node& nod = tree().nodes()[nodeI];
for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
{
const labelBits index = nod.subNodes_[octant];
if (indexedOctree::isNode(index))
{
// tree node. Recurse.
collectLeafMids
(
indexedOctree::getNode(index),
midPoints
);
}
else if (indexedOctree::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& 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::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::isNode(index))
{
// tree node. Recurse.
subType = calcVolumeType
(
midPointTypes,
midPointi,
nodeTypes,
indexedOctree::getNode(index)
);
}
else if (indexedOctree::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::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::isNode(index))
{
// Recurse
volumeType subType = cachedVolumeType
(
indexedOctree::getNode(index),
sample
);
return subType;
}
else if (indexedOctree::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::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
(
// == 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> 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(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)
{
treeBoundBox& bb = bbs[distribution[trii]][0];
bb.add(s.points(), s[trii]);
}
// Now combine for all processors and convert to correct format.
forAll(bbs, proci)
{
Pstream::listCombineGather(bbs[proci], plusEqOp());
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)
{
treeBoundBox& bb = bbs[distribution[trii]][0];
bb.add(s.points(), s[trii]);
}
// Now combine for all processors and convert to correct format.
forAll(bbs, proci)
{
Pstream::listCombineGather(bbs[proci], plusEqOp());
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 borderTris(s.size()/Pstream::nProcs());
//
// const List& 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 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)
// {
// treeBoundBox& bb = bbs[distribution[trii]][0];
// bb.add(s.points(), s[trii]);
// }
//
// // Now combine for all processors and convert to correct format.
// forAll(bbs, proci)
// {
// Pstream::listCombineGather(bbs[proci], plusEqOp());
// 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(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
(
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)
{
treeBoundBox& bb = bbs[distribution[trii]][0];
bb.add(s.points(), s[trii]);
}
// Now combine for all processors and convert to correct format.
forAll(bbs, proci)
{
Pstream::listCombineGather(bbs[proci], plusEqOp());
Pstream::listCombineScatter(bbs[proci]);
}
}
return bbs;
}
// Does any part of triangle overlap bb.
bool Foam::distributedTriSurfaceMesh::overlaps
(
const List& 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
newToOldFaces[facei++] = oldFacei;
// Renumber labels for face
for (const label oldPointi : s[oldFacei])
{
if (oldToNewPoints[oldPointi] == -1)
{
oldToNewPoints[oldPointi] = pointi;
newToOldPoints[pointi++] = oldPointi;
}
}
}
}
newToOldPoints.setSize(pointi);
}
}
Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
const triSurface& s,
const labelList& newToOldPoints,
const labelList& oldToNewPoints,
const labelList& newToOldFaces
)
{
// Extract points
pointField newPoints(newToOldPoints.size());
forAll(newToOldPoints, i)
{
newPoints[i] = s.points()[newToOldPoints[i]];
}
// Extract faces
List newTriangles(newToOldFaces.size());
forAll(newToOldFaces, i)
{
// Get old vertex labels
const labelledTri& tri = s[newToOldFaces[i]];
newTriangles[i][0] = oldToNewPoints[tri[0]];
newTriangles[i][1] = oldToNewPoints[tri[1]];
newTriangles[i][2] = oldToNewPoints[tri[2]];
newTriangles[i].region() = tri.region();
}
// Reuse storage.
return triSurface(newTriangles, s.patches(), newPoints, true);
}
Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
const triSurface& s,
const boolList& include,
labelList& newToOldPoints,
labelList& newToOldFaces
)
{
label n = 0;
forAll(include, i)
{
if (include[i])
{
n++;
}
}
labelList oldToNewPoints;
subsetMeshMap
(
s,
include,
n,
newToOldPoints,
oldToNewPoints,
newToOldFaces
);
return subsetMesh
(
s,
newToOldPoints,
oldToNewPoints,
newToOldFaces
);
}
Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
const triSurface& s,
const labelList& newToOldFaces,
labelList& newToOldPoints
)
{
const boolList include
(
ListOps::createWithValue(s.size(), newToOldFaces, true, false)
);
newToOldPoints.setSize(s.points().size());
labelList oldToNewPoints(s.points().size(), -1);
{
label pointi = 0;
forAll(include, oldFacei)
{
if (include[oldFacei])
{
// Renumber labels for face
for (const label oldPointi : s[oldFacei])
{
if (oldToNewPoints[oldPointi] == -1)
{
oldToNewPoints[oldPointi] = pointi;
newToOldPoints[pointi++] = oldPointi;
}
}
}
}
newToOldPoints.setSize(pointi);
}
return subsetMesh
(
s,
newToOldPoints,
oldToNewPoints,
newToOldFaces
);
}
Foam::label Foam::distributedTriSurfaceMesh::findTriangle
(
const List& allFaces,
const labelListList& allPointFaces,
const labelledTri& otherF
)
{
// allFaces connected to otherF[0]
const labelList& pFaces = allPointFaces[otherF[0]];
forAll(pFaces, i)
{
const labelledTri& f = allFaces[pFaces[i]];
if (f.region() == otherF.region())
{
// Find index of otherF[0]
label fp0 = f.find(otherF[0]);
// Check rest of triangle in same order
label fp1 = f.fcIndex(fp0);
label fp2 = f.fcIndex(fp1);
if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
{
return pFaces[i];
}
}
}
return -1;
}
// Merge into allSurf.
void Foam::distributedTriSurfaceMesh::merge
(
const scalar mergeDist,
const List& subTris,
const pointField& subPoints,
List& allTris,
pointField& allPoints,
labelList& faceConstructMap,
labelList& pointConstructMap
)
{
labelList subToAll;
matchPoints
(
subPoints,
allPoints,
scalarField(subPoints.size(), mergeDist), // match distance
false, // verbose
pointConstructMap
);
label nOldAllPoints = allPoints.size();
// Add all unmatched points
// ~~~~~~~~~~~~~~~~~~~~~~~~
label allPointi = nOldAllPoints;
forAll(pointConstructMap, pointi)
{
if (pointConstructMap[pointi] == -1)
{
pointConstructMap[pointi] = allPointi++;
}
}
if (allPointi > nOldAllPoints)
{
allPoints.setSize(allPointi);
forAll(pointConstructMap, pointi)
{
if (pointConstructMap[pointi] >= nOldAllPoints)
{
allPoints[pointConstructMap[pointi]] = subPoints[pointi];
}
}
}
// To check whether triangles are same we use pointFaces.
labelListList allPointFaces;
invertManyToMany(nOldAllPoints, allTris, allPointFaces);
// Add all unmatched triangles
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
label allTrii = allTris.size();
allTris.setSize(allTrii + subTris.size());
faceConstructMap.setSize(subTris.size());
forAll(subTris, trii)
{
const labelledTri& subTri = subTris[trii];
// Get triangle in new numbering
labelledTri mappedTri
(
pointConstructMap[subTri[0]],
pointConstructMap[subTri[1]],
pointConstructMap[subTri[2]],
subTri.region()
);
// Check if all points were already in surface
bool fullMatch = true;
forAll(mappedTri, fp)
{
if (mappedTri[fp] >= nOldAllPoints)
{
fullMatch = false;
break;
}
}
if (fullMatch)
{
// All three points are mapped to old points. See if same
// triangle.
label i = findTriangle
(
allTris,
allPointFaces,
mappedTri
);
if (i == -1)
{
// Add
faceConstructMap[trii] = allTrii;
allTris[allTrii] = mappedTri;
allTrii++;
}
else
{
faceConstructMap[trii] = i;
}
}
else
{
// Add
faceConstructMap[trii] = allTrii;
allTris[allTrii] = mappedTri;
allTrii++;
}
}
allTris.setSize(allTrii);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
(
const IOobject& io,
const triSurface& s,
const dictionary& dict
)
:
triSurfaceMesh(io, s),
dict_
(
IOobject
(
searchableSurface::name() + "Dict",
searchableSurface::instance(),
searchableSurface::local(),
searchableSurface::db(),
searchableSurface::NO_READ,
searchableSurface::writeOpt(),
searchableSurface::registerObject()
),
dict
),
currentDistType_(FROZEN) // only used to trigger re-distribution
{
// Read from the provided dictionary
read();
bounds().reduce();
if (debug)
{
InfoInFunction << "Constructed from triSurface:" << endl;
writeStats(Info);
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;
}
}
Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
:
triSurfaceMesh
(
IOobject
(
io.name(),
findLocalInstance(io), // findInstance with parent searching
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
io.registerObject()
),
triSurfaceMesh::masterOnly // allow parent searching
),
dict_
(
IOobject
(
searchableSurface::name() + "Dict",
searchableSurface::instance(),
searchableSurface::local(),
searchableSurface::db(),
(
(
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();
const fileName actualFile(triSurfaceMesh::checkFile(io, true));
if
(
actualFile != io.localFilePath(triSurfaceMesh::typeName)
&& (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
)
{
if (debug)
{
InfoInFunction << "Read distributedTriSurface " << io.name()
<< " from parent path " << actualFile << endl;
}
if (Pstream::parRun())
{
// Distribute (checks that distType != currentDistType_ so should
// always trigger re-distribution)
List bbs;
autoPtr faceMap;
autoPtr pointMap;
distribute
(
bbs,
true, // keep unmapped triangles
faceMap,
pointMap
);
}
}
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);
}
}
Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
(
const IOobject& io,
const dictionary& dict
)
:
triSurfaceMesh
(
IOobject
(
io.name(),
findLocalInstance(io),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
io.registerObject()
),
dict,
triSurfaceMesh::masterOnly
),
dict_
(
IOobject
(
searchableSurface::name() + "Dict",
searchableSurface::instance(),
searchableSurface::local(),
searchableSurface::db(),
(
(
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
distributionTypeNames_.readIfPresent
(
"distributionType",
dict_,
distType_
);
// Merge distance
dict_.readIfPresent("mergeDistance", mergeDist_);
// Distribution type
bool closed;
if (dict_.readIfPresent("closed", closed))
{
surfaceClosed_ = closed;
}
outsideVolType_ = volumeType::names.getOrDefault
(
"outsideVolumeType",
dict_,
outsideVolType_
);
}
bounds().reduce();
const fileName actualFile(triSurfaceMesh::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 bbs;
autoPtr faceMap;
autoPtr 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& 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& 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())
// < returnReduce(samples.size(), sumOp())
// )
// )
// {
// // 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 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 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& 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 globalIndex& triIndexer = globalTris();
// 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
// 1. Test samples on processor(s) that contains them
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr map1Ptr;
scalarField distSqr(nearestDistSqr);
boolList procContains(Pstream::nProcs(), false);
boolList procOverlaps(Pstream::nProcs(), false);
label nOutside = 0;
{
List> dynSendMap(Pstream::nProcs());
// Pre-size. Assume samples are uniformly distributed
forAll(dynSendMap, proci)
{
dynSendMap[proci].reserve(samples.size()/Pstream::nProcs());
}
forAll(samples, samplei)
{
label minProci = -1;
Tuple2 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
//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();
}
}
}
labelListList sendMap(Pstream::nProcs());
forAll(sendMap, proci)
{
sendMap[proci].transfer(dynSendMap[proci]);
}
map1Ptr.set(new mapDistribute(std::move(sendMap)));
}
const mapDistribute& map1 = map1Ptr();
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;
}
List nearestInfo;
{
// Get the points I need to test and test locally
pointField localPoints(samples);
map1.distribute(localPoints);
scalarField localDistSqr(distSqr);
map1.distribute(localDistSqr);
List 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++;
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());
}
}
}
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(0),
samples.size(),
map1.constructMap(),
map1.constructHasFlip(),
map1.subMap(),
map1.subHasFlip(),
nearestInfo,
nearestEqOp(),
noOp(), // no flipping
nearestZero
);
// 2. Test samples on other processor(s) that overlap
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 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 map2Ptr;
{
List> dynSendMap(Pstream::nProcs());
// Work array - whether processor bb overlaps the bounding sphere.
boolList procBbOverlaps(Pstream::nProcs());
label nFound = 0;
forAll(nearestInfo, samplei)
{
const point& sample = samples[samplei];
const nearestAndDist& ni = nearestInfo[samplei];
const pointIndexHit& info = ni.first();
if (info.hit())
{
nFound++;
}
scalar d2 =
(
info.hit()
? ni.second()
: distSqr[samplei]
);
label hitProci =
(
info.hit()
? triIndexer.whichProcID(info.index())
: -1
);
// Find the processors this sample+radius overlaps.
calcOverlappingProcs(sample, d2, procBbOverlaps);
forAll(procBbOverlaps, proci)
{
if (procBbOverlaps[proci])
{
// 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 (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 localInfo;
triSurfaceMesh::findNearest(localSamples, localDistSqr, localInfo);
convertTriIndices(localInfo);
// Pack and send back
List 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::commsTypes::nonBlocking,
List(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& info
) const
{
if (!Pstream::parRun())
{
triSurfaceMesh::findNearest
(
samples,
nearestDistSqr,
regionIndices,
info
);
return;
}
addProfiling
(
findNearestRegion,
"distributedTriSurfaceMesh::findNearestRegion"
);
if (debug)
{
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
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pointField allCentres;
scalarField allRadiusSqr;
labelList allSegmentMap;
autoPtr mapPtr
(
calcLocalQueries
(
true, // also send to local processor
samples,
nearestDistSqr,
allCentres,
allRadiusSqr,
allSegmentMap
)
);
const mapDistribute& map = mapPtr();
// swap samples to local processor
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
map.distribute(allCentres);
map.distribute(allRadiusSqr);
// Do my tests
// ~~~~~~~~~~~
List allInfo(allCentres.size());
triSurfaceMesh::findNearest
(
allCentres,
allRadiusSqr,
regionIndices,
allInfo
);
convertTriIndices(allInfo);
forAll(allInfo, i)
{
if (allInfo[i].hit())
{
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);
}
}
}
// 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::findLine
(
const pointField& start,
const pointField& end,
List& info
) const
{
if (!Pstream::parRun())
{
triSurfaceMesh::findLine(start, end, info);
}
else
{
findLine
(
true, // nearestIntersection
start,
end,
info
);
}
}
void Foam::distributedTriSurfaceMesh::findLineAny
(
const pointField& start,
const pointField& end,
List& info
) const
{
if (!Pstream::parRun())
{
triSurfaceMesh::findLineAny(start, end, info);
}
else
{
findLine
(
true, // nearestIntersection
start,
end,
info
);
}
}
void Foam::distributedTriSurfaceMesh::findLineAll
(
const pointField& start,
const pointField& end,
List>& 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
// intersection onwards.
// Work array.
List hitInfo(start.size());
findLine
(
true, // nearestIntersection
start,
end,
hitInfo
);
// Tolerances:
// To find all intersections we add a small vector to the last intersection
// This is chosen such that
// - it is significant (SMALL is smallest representative relative tolerance;
// we need something bigger since we're doing calculations)
// - if the start-end vector is zero we still progress
const vectorField dirVec(end-start);
const scalarField magSqrDirVec(magSqr(dirVec));
const vectorField smallVec
(
ROOTSMALL*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
// Copy to input and compact any hits
labelList pointMap(start.size());
pointField e0(start.size());
pointField e1(start.size());
label compacti = 0;
info.setSize(hitInfo.size());
forAll(hitInfo, pointi)
{
if (hitInfo[pointi].hit())
{
info[pointi].setSize(1);
info[pointi][0] = hitInfo[pointi];
point pt = hitInfo[pointi].hitPoint() + smallVec[pointi];
if (((pt-start[pointi])&dirVec[pointi]) <= magSqrDirVec[pointi])
{
e0[compacti] = pt;
e1[compacti] = end[pointi];
pointMap[compacti] = pointi;
compacti++;
}
}
else
{
info[pointi].clear();
}
}
e0.setSize(compacti);
e1.setSize(compacti);
pointMap.setSize(compacti);
label iter = 0;
while (returnReduce(e0.size(), sumOp()) > 0)
{
findLine
(
true, // nearestIntersection
e0,
e1,
hitInfo
);
// Extract
label compacti = 0;
forAll(hitInfo, i)
{
if (hitInfo[i].hit())
{
label pointi = pointMap[i];
label sz = info[pointi].size();
info[pointi].setSize(sz+1);
info[pointi][sz] = hitInfo[i];
point pt = hitInfo[i].hitPoint() + smallVec[pointi];
// Check current coordinate along ray
scalar d = ((pt-start[pointi])&dirVec[pointi]);
// Note check for d>0. Very occasionally the octree will find
// an intersection to the left of the ray due to tolerances.
if (d > 0 && d <= magSqrDirVec[pointi])
{
e0[compacti] = pt;
e1[compacti] = end[pointi];
pointMap[compacti] = pointi;
compacti++;
}
}
}
// Trim
e0.setSize(compacti);
e1.setSize(compacti);
pointMap.setSize(compacti);
iter++;
if (iter == 1000)
{
Pout<< "distributedTriSurfaceMesh::findLineAll :"
<< " Exiting loop due to excessive number of"
<< " intersections along ray"
<< " start:" << UIndirectList(start, pointMap)
<< " end:" << UIndirectList(end, pointMap)
<< " e0:" << UIndirectList(e0, pointMap)
<< " e1:" << UIndirectList(e1, pointMap)
<< endl;
break;
}
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::findLineAll :"
<< " finished intersecting with "
<< start.size() << " rays" << endl;
}
}
void Foam::distributedTriSurfaceMesh::getRegion
(
const List& info,
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());
forAll(info, i)
{
if (info[i].hit())
{
region[i] = triSurface::operator[](info[i].index()).region();
}
else
{
region[i] = -1;
}
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::getRegion :"
<< " finished getting region for "
<< info.size() << " triangles" << endl;
}
return;
}
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr mapPtr
(
localQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
const triSurface& s = static_cast(*this);
region.setSize(triangleIndex.size());
forAll(triangleIndex, i)
{
label trii = triangleIndex[i];
region[i] = s[trii].region();
}
// Send back results
// ~~~~~~~~~~~~~~~~~
map.reverseDistribute(info.size(), region);
if (debug)
{
Pout<< "distributedTriSurfaceMesh::getRegion :"
<< " finished getting region for "
<< info.size() << " triangles" << endl;
}
}
void Foam::distributedTriSurfaceMesh::getNormal
(
const List& info,
vectorField& normal
) const
{
if (!Pstream::parRun())
{
triSurfaceMesh::getNormal(info, normal);
return;
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::getNormal :"
<< " getting normal for "
<< info.size() << " triangles" << endl;
}
addProfiling(getNormal, "distributedTriSurfaceMesh::getNormal");
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr mapPtr
(
localQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
const triSurface& s = static_cast(*this);
normal.setSize(triangleIndex.size());
forAll(triangleIndex, i)
{
label trii = triangleIndex[i];
normal[i] = s[trii].unitNormal(s.points());
}
// Send back results
// ~~~~~~~~~~~~~~~~~
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& 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 nearestInfo;
// findNearest
// (
// outsidePts,
// scalarField(1, Foam::sqr(GREAT)),
// nearestInfo
// );
//
// List 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 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& 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 nearestInfo;
findNearest
(
outsidePts,
scalarField(1, Foam::sqr(GREAT)),
nearestInfo
);
List 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& t = tree();
PackedList<2>& nt = t.nodeTypes();
const List::node>& nodes =
t.nodes();
nt.setSize(nodes.size());
nt = volumeType::UNKNOWN;
// Collect midpoints
DynamicField 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 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 outsideSamples;
// Distribute samples to relevant processors
autoPtr 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 fullSearchPoints(localPoints.size());
DynamicList 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 nearestInfo;
findNearest
(
fullSearchPoints,
scalarField(fullSearchPoints.size(), Foam::sqr(GREAT)),
nearestInfo
);
// Determine orientation (in parallel)
List 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(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;
}
}
void Foam::distributedTriSurfaceMesh::getField
(
const List& info,
labelList& values
) const
{
if (!Pstream::parRun())
{
triSurfaceMesh::getField(info, values);
return;
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::getField :"
<< " retrieving field for "
<< info.size() << " triangles" << endl;
}
addProfiling(getField, "distributedTriSurfaceMesh::getField");
const auto* fldPtr = findObject("values");
if (fldPtr)
{
const triSurfaceLabelField& fld = *fldPtr;
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr mapPtr
(
localQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
values.setSize(triangleIndex.size());
forAll(triangleIndex, i)
{
label trii = triangleIndex[i];
values[i] = fld[trii];
}
// Send back results
// ~~~~~~~~~~~~~~~~~
map.reverseDistribute(info.size(), values);
}
if (debug)
{
Pout<< "distributedTriSurfaceMesh::getField :"
<< " finished retrieving field for "
<< info.size() << " triangles" << endl;
}
}
void Foam::distributedTriSurfaceMesh::overlappingSurface
(
const triSurface& s,
const List& bbs,
boolList& includedFace
)
{
// Determine what triangles to keep.
includedFace.setSize(s.size());
includedFace = false;
// Create a slightly larger bounding box.
List bbsX(bbs.size());
const scalar eps = 1.0e-4;
forAll(bbs, i)
{
const point mid = bbs[i].centre();
const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);
bbsX[i].min() = mid - halfSpan;
bbsX[i].max() = mid + halfSpan;
}
forAll(s, trii)
{
const labelledTri& f = s[trii];
const point& p0 = s.points()[f[0]];
const point& p1 = s.points()[f[1]];
const point& p2 = s.points()[f[2]];
if (overlaps(bbsX, p0, p1, p2))
{
includedFace[trii] = true;
}
}
}
// Subset the part of surface that is overlapping bb.
Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
(
const triSurface& s,
const List& bbs,
labelList& subPointMap,
labelList& subFaceMap
)
{
// Determine what triangles to keep.
boolList includedFace;
overlappingSurface(s, bbs, includedFace);
return subsetMesh(s, includedFace, subPointMap, subFaceMap);
}
// Exchanges indices to the processor they come from.
// - calculates exchange map
// - uses map to calculate local triangle index
Foam::autoPtr
Foam::distributedTriSurfaceMesh::localQueries
(
const List& info,
labelList& triangleIndex
) const
{
triangleIndex.setSize(info.size());
const globalIndex& triIndexer = globalTris();
// Determine send map
// ~~~~~~~~~~~~~~~~~~
// Since determining which processor the query should go to is
// cheap we do a multi-pass algorithm to save some memory temporarily.
// 1. Count
labelList nSend(Pstream::nProcs(), 0);
forAll(info, i)
{
if (info[i].hit())
{
label proci = triIndexer.whichProcID(info[i].index());
nSend[proci]++;
}
}
// 2. Size sendMap
labelListList sendMap(Pstream::nProcs());
forAll(nSend, proci)
{
sendMap[proci].setSize(nSend[proci]);
nSend[proci] = 0;
}
// 3. Fill sendMap
forAll(info, i)
{
if (info[i].hit())
{
label proci = triIndexer.whichProcID(info[i].index());
triangleIndex[i] = triIndexer.toLocal(proci, info[i].index());
sendMap[proci][nSend[proci]++] = i;
}
else
{
triangleIndex[i] = -1;
}
}
// 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 receive map
// ~~~~~~~~~~~~~~~~~~~~~
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++;
}
}
}
// Pack into distribution map
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr mapPtr
(
new mapDistribute
(
segmenti, // size after construction
std::move(sendMap),
std::move(constructMap)
)
);
const mapDistribute& map = mapPtr();
// Send over queries
// ~~~~~~~~~~~~~~~~~
map.distribute(triangleIndex);
return mapPtr;
}
void Foam::distributedTriSurfaceMesh::distribute
(
const List& bbs,
const bool keepNonLocal,
autoPtr& faceMap,
autoPtr& 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
// ~~~~~~~~~~~~~~~~~~~~~~
{
List> newProcBb(Pstream::nProcs());
switch(distType_)
{
case FOLLOW:
newProcBb[Pstream::myProcNo()].setSize(bbs.size());
forAll(bbs, i)
{
newProcBb[Pstream::myProcNo()][i] = bbs[i];
}
Pstream::gatherList(newProcBb);
Pstream::scatterList(newProcBb);
break;
case INDEPENDENT:
case DISTRIBUTED:
if (currentDistType_ == distType_)
{
return;
}
newProcBb = independentlyDistributedBbs(*this);
break;
case FROZEN:
return;
break;
default:
FatalErrorInFunction
<< "Unsupported distribution type." << exit(FatalError);
break;
}
if (newProcBb == procBb_)
{
return;
}
else
{
procBb_.transfer(newProcBb);
dict_.set("bounds", procBb_[Pstream::myProcNo()]);
}
}
// Debug information
if (debug)
{
labelList nTris(Pstream::nProcs());
nTris[Pstream::myProcNo()] = triSurface::size();
Pstream::gatherList(nTris);
Pstream::scatterList(nTris);
InfoInFunction
<< "before distribution:" << endl << "\tproc\ttris" << endl;
forAll(nTris, proci)
{
Info<< '\t' << proci << '\t' << nTris[proci] << endl;
}
Info<< endl;
}
// Use procBbs to determine which faces go where
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelListList faceSendMap(Pstream::nProcs());
labelListList pointSendMap(Pstream::nProcs());
forAll(procBb_, proci)
{
overlappingSurface
(
*this,
procBb_[proci],
pointSendMap[proci],
faceSendMap[proci]
);
}
if (keepNonLocal)
{
// Include in faceSendMap/pointSendMap the triangles that are
// not mapped to any processor so they stay local.
const triSurface& s = static_cast(*this);
boolList includedFace(s.size(), true);
forAll(faceSendMap, proci)
{
if (proci != Pstream::myProcNo())
{
forAll(faceSendMap[proci], i)
{
includedFace[faceSendMap[proci][i]] = false;
}
}
}
// Combine includedFace (all triangles that are not on any neighbour)
// with overlap.
forAll(faceSendMap[Pstream::myProcNo()], i)
{
includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
}
subsetMesh
(
s,
includedFace,
pointSendMap[Pstream::myProcNo()],
faceSendMap[Pstream::myProcNo()]
);
}
// Send over how many faces/points i need to receive
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList faceRecvSizes;
Pstream::exchangeSizes(faceSendMap, faceRecvSizes);
// Exchange surfaces
// ~~~~~~~~~~~~~~~~~
// Storage for resulting surface
List allTris;
pointField allPoints;
labelListList faceConstructMap(Pstream::nProcs());
labelListList pointConstructMap(Pstream::nProcs());
// My own surface first
// ~~~~~~~~~~~~~~~~~~~~
{
labelList pointMap;
triSurface subSurface
(
subsetMesh
(
*this,
faceSendMap[Pstream::myProcNo()],
pointMap
)
);
allTris = subSurface;
allPoints = subSurface.points();
faceConstructMap[Pstream::myProcNo()] = identity
(
faceSendMap[Pstream::myProcNo()].size()
);
pointConstructMap[Pstream::myProcNo()] = identity
(
pointSendMap[Pstream::myProcNo()].size()
);
}
// Send all
// ~~~~~~~~
PstreamBuffers pBufs(Pstream::defaultCommsType);
forAll(faceSendMap, proci)
{
if (proci != Pstream::myProcNo())
{
if (faceSendMap[proci].size() > 0)
{
UOPstream str(proci, pBufs);
labelList pointMap;
triSurface subSurface
(
subsetMesh
(
*this,
faceSendMap[proci],
pointMap
)
);
str << subSurface;
}
}
}
pBufs.finishedSends(); // no-op for blocking
// Receive and merge all
// ~~~~~~~~~~~~~~~~~~~~~
forAll(faceRecvSizes, proci)
{
if (proci != Pstream::myProcNo())
{
if (faceRecvSizes[proci] > 0)
{
UIPstream str(proci, pBufs);
// Receive
triSurface subSurface(str);
// Merge into allSurf
merge
(
mergeDist_,
subSurface,
subSurface.points(),
allTris,
allPoints,
faceConstructMap[proci],
pointConstructMap[proci]
);
}
}
}
faceMap.reset
(
new mapDistribute
(
allTris.size(),
std::move(faceSendMap),
std::move(faceConstructMap)
)
);
pointMap.reset
(
new mapDistribute
(
allPoints.size(),
std::move(pointSendMap),
std::move(pointConstructMap)
)
);
// 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.
distributeFields(faceMap());
distributeFields(faceMap());
distributeFields(faceMap());
distributeFields(faceMap());
distributeFields(faceMap());
distributeFields(faceMap());
if (debug)
{
labelList nTris(Pstream::nProcs());
nTris[Pstream::myProcNo()] = triSurface::size();
Pstream::gatherList(nTris);
Pstream::scatterList(nTris);
InfoInFunction
<< "after distribution:" << endl << "\tproc\ttris" << endl;
forAll(nTris, proci)
{
Info<< '\t' << proci << '\t' << nTris[proci] << endl;
}
Info<< endl;
if (debug & 2)
{
OBJstream str(searchableSurface::time().path()/"after.obj");
Info<< "Writing local bounding box to " << str.name() << endl;
const List& 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)
{
OBJstream str(searchableSurface::time().path()/"after_all.obj");
Info<< "Writing all bounding boxes to " << str.name() << endl;
for (auto myBbs : procBb_)
{
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;
}
}
bool Foam::distributedTriSurfaceMesh::writeObject
(
IOstreamOption streamOpt,
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(dict_.instance()) = searchableSurface::instance();
// Copy of triSurfaceMesh::writeObject except for the sorting of
// triangles by region. This is done so we preserve region names,
// even if locally we have zero triangles.
{
fileName fullPath(searchableSurface::objectPath());
if (!mkDir(fullPath.path()))
{
return false;
}
// Important: preserve any zero-sized patches
triSurface::write(fullPath, true);
if (!isFile(fullPath))
{
return false;
}
}
// Dictionary needs to be written in ascii - binary output not supported.
streamOpt.format(IOstream::ASCII);
bool ok = dict_.writeObject(streamOpt, true);
if (debug)
{
Pout<< "distributedTriSurfaceMesh::writeObject :"
<< " done writing surface" << endl;
}
return ok;
}
void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
{
boundBox bb;
label nPoints;
PatchTools::calcBounds(static_cast(*this), bb, nPoints);
bb.reduce();
os << "Triangles : " << returnReduce(triSurface::size(), sumOp())
<< endl
<< "Vertices : " << returnReduce(nPoints, sumOp()) << endl
<< "Bounding Box : " << bb << endl
<< "Closed : " << surfaceClosed_ << endl
<< "Outside point: " << volumeType::names[outsideVolType_] << endl
<< "Distribution : " << distributionTypeNames_[distType_] << endl;
}
// ************************************************************************* //