/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 .
\*---------------------------------------------------------------------------*/
#include "DistributedDelaunayMesh.H"
#include "meshSearch.H"
#include "mapDistribute.H"
#include "zeroGradientFvPatchFields.H"
#include "pointConversion.H"
#include "indexedVertexEnum.H"
#include "IOmanip.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
template
Foam::autoPtr
Foam::DistributedDelaunayMesh::buildMap
(
const List& toProc
)
{
// Determine send map
// ~~~~~~~~~~~~~~~~~~
// 1. Count
labelList nSend(Pstream::nProcs(), 0);
forAll(toProc, i)
{
label procI = toProc[i];
nSend[procI]++;
}
// Send over how many I need to receive
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()] = nSend;
combineReduce(sendSizes, UPstream::listEq());
// 2. Size sendMap
labelListList sendMap(Pstream::nProcs());
forAll(nSend, procI)
{
sendMap[procI].setSize(nSend[procI]);
nSend[procI] = 0;
}
// 3. Fill sendMap
forAll(toProc, i)
{
label procI = toProc[i];
sendMap[procI][nSend[procI]++] = i;
}
// Determine receive map
// ~~~~~~~~~~~~~~~~~~~~~
labelListList constructMap(Pstream::nProcs());
// Local transfers first
constructMap[Pstream::myProcNo()] = identity
(
sendMap[Pstream::myProcNo()].size()
);
label constructSize = constructMap[Pstream::myProcNo()].size();
forAll(constructMap, procI)
{
if (procI != Pstream::myProcNo())
{
label nRecv = sendSizes[procI][Pstream::myProcNo()];
constructMap[procI].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[procI][i] = constructSize++;
}
}
}
return autoPtr
(
new mapDistribute
(
constructSize,
sendMap.xfer(),
constructMap.xfer()
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template
Foam::DistributedDelaunayMesh::DistributedDelaunayMesh
(
const Time& runTime
)
:
DelaunayMesh(runTime),
allBackgroundMeshBounds_()
{}
template
Foam::DistributedDelaunayMesh::DistributedDelaunayMesh
(
const Time& runTime,
const word& meshName
)
:
DelaunayMesh(runTime, meshName),
allBackgroundMeshBounds_()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template
Foam::DistributedDelaunayMesh::~DistributedDelaunayMesh()
{}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template
bool Foam::DistributedDelaunayMesh::distributeBoundBoxes
(
const boundBox& bb
)
{
allBackgroundMeshBounds_.reset(new List(Pstream::nProcs()));
// Give the bounds of every processor to every other processor
allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
Pstream::gatherList(allBackgroundMeshBounds_());
Pstream::scatterList(allBackgroundMeshBounds_());
return true;
}
template
bool Foam::DistributedDelaunayMesh::isLocal
(
const Vertex_handle& v
) const
{
return isLocal(v->procIndex());
}
template
bool Foam::DistributedDelaunayMesh::isLocal
(
const label localProcIndex
) const
{
return localProcIndex == Pstream::myProcNo();
}
template
Foam::labelList Foam::DistributedDelaunayMesh::overlapProcessors
(
const point& centre,
const scalar radiusSqr
) const
{
DynamicList toProc(Pstream::nProcs());
forAll(allBackgroundMeshBounds_(), procI)
{
// Test against the bounding box of the processor
if
(
!isLocal(procI)
&& allBackgroundMeshBounds_()[procI].overlaps(centre, radiusSqr)
)
{
toProc.append(procI);
}
}
return toProc;
}
template
bool Foam::DistributedDelaunayMesh::checkProcBoundaryCell
(
const Cell_handle& cit,
Map& circumsphereOverlaps
) const
{
const Foam::point& cc = cit->dual();
const scalar crSqr = magSqr
(
cc - topoint(cit->vertex(0)->point())
);
labelList circumsphereOverlap = overlapProcessors
(
cc,
sqr(1.01)*crSqr
);
cit->cellIndex() = this->getNewCellIndex();
if (!circumsphereOverlap.empty())
{
circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
return true;
}
return false;
}
template
void Foam::DistributedDelaunayMesh::findProcessorBoundaryCells
(
Map& circumsphereOverlaps
) const
{
// Start by assuming that all the cells have no index
// If they do, they have already been visited so ignore them
labelHashSet cellToCheck
(
Triangulation::number_of_finite_cells()
/Pstream::nProcs()
);
// std::list infinite_cells;
// Triangulation::incident_cells
// (
// Triangulation::infinite_vertex(),
// std::back_inserter(infinite_cells)
// );
//
// for
// (
// typename std::list::iterator vcit
// = infinite_cells.begin();
// vcit != infinite_cells.end();
// ++vcit
// )
// {
// Cell_handle cit = *vcit;
//
// // Index of infinite vertex in this cell.
// label i = cit->index(Triangulation::infinite_vertex());
//
// Cell_handle c = cit->neighbor(i);
//
// if (c->unassigned())
// {
// c->cellIndex() = this->getNewCellIndex();
//
// if (checkProcBoundaryCell(c, circumsphereOverlaps))
// {
// cellToCheck.insert(c->cellIndex());
// }
// }
// }
//
//
// for
// (
// Finite_cells_iterator cit = Triangulation::finite_cells_begin();
// cit != Triangulation::finite_cells_end();
// ++cit
// )
// {
// if (cit->parallelDualVertex())
// {
// if (cit->unassigned())
// {
// if (checkProcBoundaryCell(cit, circumsphereOverlaps))
// {
// cellToCheck.insert(cit->cellIndex());
// }
// }
// }
// }
for
(
All_cells_iterator cit = Triangulation::all_cells_begin();
cit != Triangulation::all_cells_end();
++cit
)
{
if (Triangulation::is_infinite(cit))
{
// Index of infinite vertex in this cell.
label i = cit->index(Triangulation::infinite_vertex());
Cell_handle c = cit->neighbor(i);
if (c->unassigned())
{
c->cellIndex() = this->getNewCellIndex();
if (checkProcBoundaryCell(c, circumsphereOverlaps))
{
cellToCheck.insert(c->cellIndex());
}
}
}
else if (cit->parallelDualVertex())
{
if (cit->unassigned())
{
if (checkProcBoundaryCell(cit, circumsphereOverlaps))
{
cellToCheck.insert(cit->cellIndex());
}
}
}
}
for
(
Finite_cells_iterator cit = Triangulation::finite_cells_begin();
cit != Triangulation::finite_cells_end();
++cit
)
{
if (cellToCheck.found(cit->cellIndex()))
{
// Get the neighbours and check them
for (label adjCellI = 0; adjCellI < 4; ++adjCellI)
{
Cell_handle citNeighbor = cit->neighbor(adjCellI);
// Ignore if has far point or previously visited
if
(
!citNeighbor->unassigned()
|| !citNeighbor->internalOrBoundaryDualVertex()
|| Triangulation::is_infinite(citNeighbor)
)
{
continue;
}
if
(
checkProcBoundaryCell
(
citNeighbor,
circumsphereOverlaps
)
)
{
cellToCheck.insert(citNeighbor->cellIndex());
}
}
cellToCheck.unset(cit->cellIndex());
}
}
}
template
void Foam::DistributedDelaunayMesh::markVerticesToRefer
(
const Map& circumsphereOverlaps,
PtrList& referralVertices,
DynamicList& targetProcessor,
DynamicList& parallelInfluenceVertices
)
{
// Relying on the order of iteration of cells being the same as before
for
(
Finite_cells_iterator cit = Triangulation::finite_cells_begin();
cit != Triangulation::finite_cells_end();
++cit
)
{
if (Triangulation::is_infinite(cit))
{
continue;
}
Map::const_iterator iter =
circumsphereOverlaps.find(cit->cellIndex());
// Pre-tested circumsphere potential influence
if (iter != circumsphereOverlaps.cend())
{
const labelList& citOverlaps = iter();
forAll(citOverlaps, cOI)
{
label procI = citOverlaps[cOI];
for (int i = 0; i < 4; i++)
{
Vertex_handle v = cit->vertex(i);
if (v->farPoint())
{
continue;
}
label vProcIndex = v->procIndex();
label vIndex = v->index();
const labelPair procIndexPair(vProcIndex, vIndex);
// Using the hashSet to ensure that each vertex is only
// referred once to each processor.
// Do not refer a vertex to its own processor.
if (vProcIndex != procI)
{
if (referralVertices[procI].insert(procIndexPair))
{
targetProcessor.append(procI);
parallelInfluenceVertices.append
(
Vb
(
v->point(),
v->index(),
v->type(),
v->procIndex()
)
);
parallelInfluenceVertices.last().targetCellSize() =
v->targetCellSize();
parallelInfluenceVertices.last().alignment() =
v->alignment();
}
}
}
}
}
}
}
template
Foam::label Foam::DistributedDelaunayMesh::referVertices
(
const DynamicList& targetProcessor,
DynamicList& parallelVertices,
PtrList& referralVertices,
labelPairHashSet& receivedVertices
)
{
DynamicList referredVertices(targetProcessor.size());
const label preDistributionSize = parallelVertices.size();
mapDistribute pointMap = buildMap(targetProcessor);
// Make a copy of the original list.
DynamicList originalParallelVertices(parallelVertices);
pointMap.distribute(parallelVertices);
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const labelList& constructMap = pointMap.constructMap()[procI];
if (constructMap.size())
{
forAll(constructMap, i)
{
const Vb& v = parallelVertices[constructMap[i]];
if
(
v.procIndex() != Pstream::myProcNo()
&& !receivedVertices.found(labelPair(v.procIndex(), v.index()))
)
{
referredVertices.append(v);
receivedVertices.insert
(
labelPair(v.procIndex(), v.index())
);
}
}
}
}
label preInsertionSize = Triangulation::number_of_vertices();
labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
(
referredVertices.begin(),
referredVertices.end(),
true
);
if (!pointsNotInserted.empty())
{
for
(
typename labelPairHashSet::const_iterator iter
= pointsNotInserted.begin();
iter != pointsNotInserted.end();
++iter
)
{
if (receivedVertices.found(iter.key()))
{
receivedVertices.erase(iter.key());
}
}
}
boolList pointInserted(parallelVertices.size(), true);
forAll(parallelVertices, vI)
{
const labelPair procIndexI
(
parallelVertices[vI].procIndex(),
parallelVertices[vI].index()
);
if (pointsNotInserted.found(procIndexI))
{
pointInserted[vI] = false;
}
}
pointMap.reverseDistribute(preDistributionSize, pointInserted);
forAll(originalParallelVertices, vI)
{
const label procIndex = targetProcessor[vI];
if (!pointInserted[vI])
{
if (referralVertices[procIndex].size())
{
if
(
!referralVertices[procIndex].unset
(
labelPair
(
originalParallelVertices[vI].procIndex(),
originalParallelVertices[vI].index()
)
)
)
{
Pout<< "*** not found "
<< originalParallelVertices[vI].procIndex()
<< " " << originalParallelVertices[vI].index() << endl;
}
}
}
}
label postInsertionSize = Triangulation::number_of_vertices();
reduce(preInsertionSize, sumOp());
reduce(postInsertionSize, sumOp());
label nTotalToInsert = referredVertices.size();
reduce(nTotalToInsert, sumOp());
if (preInsertionSize + nTotalToInsert != postInsertionSize)
{
label nNotInserted =
returnReduce(pointsNotInserted.size(), sumOp());
Info<< " Inserted = "
<< setw(name(label(Triangulation::number_of_finite_cells())).size())
<< nTotalToInsert - nNotInserted
<< " / " << nTotalToInsert << endl;
nTotalToInsert -= nNotInserted;
}
else
{
Info<< " Inserted = " << nTotalToInsert << endl;
}
return nTotalToInsert;
}
template
void Foam::DistributedDelaunayMesh::sync
(
const boundBox& bb,
PtrList& referralVertices,
labelPairHashSet& receivedVertices,
bool iterateReferral
)
{
if (!Pstream::parRun())
{
return;
}
if (allBackgroundMeshBounds_.empty())
{
distributeBoundBoxes(bb);
}
label nVerts = Triangulation::number_of_vertices();
label nCells = Triangulation::number_of_finite_cells();
DynamicList parallelInfluenceVertices(0.1*nVerts);
DynamicList targetProcessor(0.1*nVerts);
// Some of these values will not be used, i.e. for non-real cells
DynamicList circumcentre(0.1*nVerts);
DynamicList circumradiusSqr(0.1*nVerts);
Map circumsphereOverlaps(nCells);
findProcessorBoundaryCells(circumsphereOverlaps);
Info<< " Influences = "
<< setw(name(nCells).size())
<< returnReduce(circumsphereOverlaps.size(), sumOp()) << " / "
<< returnReduce(nCells, sumOp());
markVerticesToRefer
(
circumsphereOverlaps,
referralVertices,
targetProcessor,
parallelInfluenceVertices
);
referVertices
(
targetProcessor,
parallelInfluenceVertices,
referralVertices,
receivedVertices
);
if (iterateReferral)
{
label oldNReferred = 0;
label nIterations = 1;
Info<< incrIndent << indent
<< "Iteratively referring referred vertices..."
<< endl;
do
{
Info<< indent << "Iteration " << nIterations++ << ":";
circumsphereOverlaps.clear();
targetProcessor.clear();
parallelInfluenceVertices.clear();
findProcessorBoundaryCells(circumsphereOverlaps);
nCells = Triangulation::number_of_finite_cells();
Info<< " Influences = "
<< setw(name(nCells).size())
<< returnReduce(circumsphereOverlaps.size(), sumOp())
<< " / "
<< returnReduce(nCells, sumOp());
markVerticesToRefer
(
circumsphereOverlaps,
referralVertices,
targetProcessor,
parallelInfluenceVertices
);
label nReferred = referVertices
(
targetProcessor,
parallelInfluenceVertices,
referralVertices,
receivedVertices
);
if (nReferred == 0 || nReferred == oldNReferred)
{
break;
}
oldNReferred = nReferred;
} while (true);
Info<< decrIndent;
}
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template
Foam::scalar
Foam::DistributedDelaunayMesh::calculateLoadUnbalance() const
{
label nRealVertices = 0;
for
(
Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
vit != Triangulation::finite_vertices_end();
++vit
)
{
// Only store real vertices that are not feature vertices
if (vit->real() && !vit->featurePoint())
{
nRealVertices++;
}
}
scalar globalNRealVertices = returnReduce
(
nRealVertices,
sumOp()
);
scalar unbalance = returnReduce
(
mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
maxOp()
);
Info<< " Processor unbalance " << unbalance << endl;
return unbalance;
}
template
bool Foam::DistributedDelaunayMesh::distribute
(
const boundBox& bb
)
{
NotImplemented;
if (!Pstream::parRun())
{
return false;
}
distributeBoundBoxes(bb);
return true;
}
template
Foam::autoPtr
Foam::DistributedDelaunayMesh::distribute
(
const backgroundMeshDecomposition& decomposition,
List& points
)
{
if (!Pstream::parRun())
{
return autoPtr();
}
distributeBoundBoxes(decomposition.procBounds());
autoPtr mapDist = decomposition.distributePoints(points);
return mapDist;
}
template
void Foam::DistributedDelaunayMesh::sync(const boundBox& bb)
{
if (!Pstream::parRun())
{
return;
}
if (allBackgroundMeshBounds_.empty())
{
distributeBoundBoxes(bb);
}
const label nApproxReferred =
Triangulation::number_of_vertices()
/Pstream::nProcs();
PtrList referralVertices(Pstream::nProcs());
forAll(referralVertices, procI)
{
if (!isLocal(procI))
{
referralVertices.set(procI, new labelPairHashSet(nApproxReferred));
}
}
labelPairHashSet receivedVertices(nApproxReferred);
sync
(
bb,
referralVertices,
receivedVertices,
true
);
}
template
template
typename Foam::DistributedDelaunayMesh::labelPairHashSet
Foam::DistributedDelaunayMesh::rangeInsertReferredWithInfo
(
PointIterator begin,
PointIterator end,
bool printErrors
)
{
const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
typedef DynamicList
<
std::pair
> vectorPairPointIndex;
vectorPairPointIndex pointsBbDistSqr;
label count = 0;
for (PointIterator it = begin; it != end; ++it)
{
const Foam::point samplePoint(topoint(it->point()));
scalar distFromBbSqr = 0;
if (!bb.contains(samplePoint))
{
const Foam::point nearestPoint = bb.nearest(samplePoint);
distFromBbSqr = magSqr(nearestPoint - samplePoint);
}
pointsBbDistSqr.append
(
std::make_pair(distFromBbSqr, count++)
);
}
std::random_shuffle(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
// Sort in ascending order by the distance of the point from the centre
// of the processor bounding box
sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
typename Triangulation::Vertex_handle hint;
typename Triangulation::Locate_type lt;
int li, lj;
label nNotInserted = 0;
labelPairHashSet uninserted
(
Triangulation::number_of_vertices()
/Pstream::nProcs()
);
for
(
typename vectorPairPointIndex::const_iterator p =
pointsBbDistSqr.begin();
p != pointsBbDistSqr.end();
++p
)
{
const size_t checkInsertion = Triangulation::number_of_vertices();
const Vb& vert = *(begin + p->second);
const Point& pointToInsert = vert.point();
// Locate the point
Cell_handle c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
bool inserted = false;
if (lt == Triangulation::VERTEX)
{
if (printErrors)
{
Vertex_handle nearV =
Triangulation::nearest_vertex(pointToInsert);
Pout<< "Failed insertion, point already exists" << nl
<< "Failed insertion : " << vert.info()
<< " nearest : " << nearV->info();
}
}
else if (lt == Triangulation::OUTSIDE_AFFINE_HULL)
{
WarningIn
(
"Foam::DistributedDelaunayMesh"
"::rangeInsertReferredWithInfo"
) << "Point is outside affine hull! pt = " << pointToInsert
<< endl;
}
else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
{
// @todo Can this be optimised?
//
// Only want to insert if a connection is formed between
// pointToInsert and an internal or internal boundary point.
hint = Triangulation::insert(pointToInsert, c);
inserted = true;
}
else
{
// Get the cells that conflict with p in a vector V,
// and a facet on the boundary of this hole in f.
std::vector V;
typename Triangulation::Facet f;
Triangulation::find_conflicts
(
pointToInsert,
c,
CGAL::Oneset_iterator(f),
std::back_inserter(V)
);
for (size_t i = 0; i < V.size(); ++i)
{
Cell_handle conflictingCell = V[i];
if
(
Triangulation::dimension() < 3 // 2D triangulation
||
(
!Triangulation::is_infinite(conflictingCell)
&& (
conflictingCell->real()
|| conflictingCell->hasFarPoint()
)
)
)
{
hint = Triangulation::insert_in_hole
(
pointToInsert,
V.begin(),
V.end(),
f.first,
f.second
);
inserted = true;
break;
}
}
}
if (inserted)
{
if (checkInsertion != Triangulation::number_of_vertices() - 1)
{
if (printErrors)
{
Vertex_handle nearV =
Triangulation::nearest_vertex(pointToInsert);
Pout<< "Failed insertion : " << vert.info()
<< " nearest : " << nearV->info();
}
}
else
{
hint->index() = vert.index();
hint->type() = vert.type();
hint->procIndex() = vert.procIndex();
hint->targetCellSize() = vert.targetCellSize();
hint->alignment() = vert.alignment();
}
}
else
{
uninserted.insert(labelPair(vert.procIndex(), vert.index()));
nNotInserted++;
}
}
return uninserted;
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
// ************************************************************************* //