41 template<
class Triangulation>
45 const List<label>& toProc
56 label proci = toProc[i];
66 sendMap[proci].resize_nocopy(nSend[proci]);
73 label proci = toProc[i];
74 sendMap[proci][nSend[proci]++] = i;
83 template<
class Triangulation>
89 DelaunayMesh<Triangulation>(
runTime),
90 allBackgroundMeshBounds_()
94 template<
class Triangulation>
101 DelaunayMesh<Triangulation>(
runTime, meshName),
102 allBackgroundMeshBounds_()
108 template<
class Triangulation>
114 allBackgroundMeshBounds_.
reset(
new List<boundBox>(Pstream::nProcs()));
117 allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
119 Pstream::allGatherList(allBackgroundMeshBounds_());
125 template<
class Triangulation>
128 const Vertex_handle& v
131 return isLocal(v->procIndex());
135 template<
class Triangulation>
138 const label localProcIndex
141 return localProcIndex == Pstream::myProcNo();
145 template<
class Triangulation>
149 const scalar radiusSqr
152 DynamicList<label> toProc(Pstream::nProcs());
154 forAll(allBackgroundMeshBounds_(), proci)
160 && allBackgroundMeshBounds_()[proci].overlaps(centre, radiusSqr)
163 toProc.append(proci);
171 template<
class Triangulation>
174 const Cell_handle& cit,
175 Map<labelList>& circumsphereOverlaps
180 const scalar crSqr =
magSqr 182 cc -
topoint(cit->vertex(0)->point())
185 labelList circumsphereOverlap = overlapProcessors
191 cit->cellIndex() = this->getNewCellIndex();
193 if (!circumsphereOverlap.empty())
195 circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
204 template<
class Triangulation>
207 Map<labelList>& circumsphereOverlaps
215 Triangulation::number_of_finite_cells()
275 All_cells_iterator cit = Triangulation::all_cells_begin();
276 cit != Triangulation::all_cells_end();
280 if (Triangulation::is_infinite(cit))
283 label i = cit->index(Triangulation::infinite_vertex());
285 Cell_handle
c = cit->neighbor(i);
289 c->cellIndex() = this->getNewCellIndex();
291 if (checkProcBoundaryCell(
c, circumsphereOverlaps))
293 cellToCheck.insert(
c->cellIndex());
297 else if (cit->parallelDualVertex())
299 if (cit->unassigned())
301 if (checkProcBoundaryCell(cit, circumsphereOverlaps))
303 cellToCheck.insert(cit->cellIndex());
311 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
312 cit != Triangulation::finite_cells_end();
316 if (cellToCheck.found(cit->cellIndex()))
319 for (label adjCelli = 0; adjCelli < 4; ++adjCelli)
321 Cell_handle citNeighbor = cit->neighbor(adjCelli);
326 !citNeighbor->unassigned()
327 || !citNeighbor->internalOrBoundaryDualVertex()
328 || Triangulation::is_infinite(citNeighbor)
336 checkProcBoundaryCell
343 cellToCheck.insert(citNeighbor->cellIndex());
347 cellToCheck.unset(cit->cellIndex());
353 template<
class Triangulation>
356 const Map<labelList>& circumsphereOverlaps,
357 PtrList<labelPairHashSet>& referralVertices,
358 DynamicList<label>& targetProcessor,
359 DynamicList<Vb>& parallelInfluenceVertices
365 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
366 cit != Triangulation::finite_cells_end();
370 if (Triangulation::is_infinite(cit))
375 const auto iter = circumsphereOverlaps.cfind(cit->cellIndex());
382 for (
const label proci : citOverlaps)
384 for (
int i = 0; i < 4; i++)
386 Vertex_handle v = cit->vertex(i);
393 label vProcIndex = v->procIndex();
394 label vIndex = v->index();
396 const labelPair procIndexPair(vProcIndex, vIndex);
401 if (vProcIndex != proci)
403 if (referralVertices[proci].
insert(procIndexPair))
405 targetProcessor.append(proci);
407 parallelInfluenceVertices.append
418 parallelInfluenceVertices.last().targetCellSize() =
420 parallelInfluenceVertices.last().alignment() =
431 template<
class Triangulation>
434 const DynamicList<label>& targetProcessor,
435 DynamicList<Vb>& parallelVertices,
436 PtrList<labelPairHashSet>& referralVertices,
440 DynamicList<Vb> referredVertices(targetProcessor.size());
442 const label preDistributionSize = parallelVertices.size();
444 autoPtr<mapDistribute> pointMapPtr = buildMap(targetProcessor);
445 mapDistribute& pointMap = *pointMapPtr;
448 DynamicList<Vb> originalParallelVertices(parallelVertices);
450 pointMap.distribute(parallelVertices);
452 for (
const int proci : Pstream::allProcs())
454 const labelList& constructMap = pointMap.constructMap()[proci];
456 if (constructMap.size())
460 const Vb& v = parallelVertices[constructMap[i]];
468 referredVertices.append(v);
470 receivedVertices.insert
479 label preInsertionSize = Triangulation::number_of_vertices();
483 referredVertices.begin(),
484 referredVertices.end(),
488 if (!pointsNotInserted.empty())
492 if (receivedVertices.found(iter.key()))
494 receivedVertices.erase(iter.key());
499 boolList pointInserted(parallelVertices.size(),
true);
501 forAll(parallelVertices, vI)
505 parallelVertices[vI].procIndex(),
506 parallelVertices[vI].index()
509 if (pointsNotInserted.found(procIndexI))
511 pointInserted[vI] =
false;
515 pointMap.reverseDistribute(preDistributionSize, pointInserted);
517 forAll(originalParallelVertices, vI)
519 const label procIndex = targetProcessor[vI];
521 if (!pointInserted[vI])
523 if (referralVertices[procIndex].size())
527 !referralVertices[procIndex].
unset 531 originalParallelVertices[vI].procIndex(),
532 originalParallelVertices[vI].index()
537 Pout<<
"*** not found " 538 << originalParallelVertices[vI].procIndex()
539 <<
" " << originalParallelVertices[vI].index() <<
endl;
545 label postInsertionSize = Triangulation::number_of_vertices();
547 reduce(preInsertionSize, sumOp<label>());
548 reduce(postInsertionSize, sumOp<label>());
550 label nTotalToInsert =
553 if (preInsertionSize + nTotalToInsert != postInsertionSize)
558 Info<<
" Inserted = " 559 <<
setw(
name(label(Triangulation::number_of_finite_cells())).size())
560 << nTotalToInsert - nNotInserted
561 <<
" / " << nTotalToInsert <<
endl;
563 nTotalToInsert -= nNotInserted;
567 Info<<
" Inserted = " << nTotalToInsert <<
endl;
570 return nTotalToInsert;
574 template<
class Triangulation>
578 PtrList<labelPairHashSet>& referralVertices,
583 if (!Pstream::parRun())
588 if (!allBackgroundMeshBounds_)
590 distributeBoundBoxes(bb);
593 label nVerts = Triangulation::number_of_vertices();
594 label nCells = Triangulation::number_of_finite_cells();
596 DynamicList<Vb> parallelInfluenceVertices(0.1*nVerts);
597 DynamicList<label> targetProcessor(0.1*nVerts);
600 DynamicList<Foam::point> circumcentre(0.1*nVerts);
601 DynamicList<scalar> circumradiusSqr(0.1*nVerts);
603 Map<labelList> circumsphereOverlaps(nCells);
605 findProcessorBoundaryCells(circumsphereOverlaps);
607 Info<<
" Influences = " 609 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>()) <<
" / " 614 circumsphereOverlaps,
617 parallelInfluenceVertices
623 parallelInfluenceVertices,
630 label oldNReferred = 0;
631 label nIterations = 1;
634 <<
"Iteratively referring referred vertices..." 638 Info<<
indent <<
"Iteration " << nIterations++ <<
":";
640 circumsphereOverlaps.clear();
641 targetProcessor.clear();
642 parallelInfluenceVertices.clear();
644 findProcessorBoundaryCells(circumsphereOverlaps);
646 nCells = Triangulation::number_of_finite_cells();
648 Info<<
" Influences = " 650 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>())
656 circumsphereOverlaps,
659 parallelInfluenceVertices
662 label nReferred = referVertices
665 parallelInfluenceVertices,
670 if (nReferred == 0 || nReferred == oldNReferred)
675 oldNReferred = nReferred;
689 template<
class Triangulation>
693 label nRealVertices = 0;
697 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
698 vit != Triangulation::finite_vertices_end();
703 if (vit->real() && !vit->featurePoint())
717 mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
721 Info<<
" Processor unbalance " << unbalance <<
endl;
727 template<
class Triangulation>
735 if (!Pstream::parRun())
740 distributeBoundBoxes(bb);
746 template<
class Triangulation>
750 const backgroundMeshDecomposition& decomposition,
754 if (!Pstream::parRun())
759 distributeBoundBoxes(decomposition.procBounds());
761 return decomposition.distributePoints(
points);
765 template<
class Triangulation>
768 if (!Pstream::parRun())
773 if (!allBackgroundMeshBounds_)
775 distributeBoundBoxes(bb);
778 const label nApproxReferred =
779 Triangulation::number_of_vertices()
782 PtrList<labelPairHashSet> referralVertices(Pstream::nProcs());
783 forAll(referralVertices, proci)
803 template<
class Triangulation>
804 template<
class Po
intIterator>
813 const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
817 std::pair<scalar, label>
818 > vectorPairPointIndex;
820 vectorPairPointIndex pointsBbDistSqr;
823 for (PointIterator it =
begin; it !=
end; ++it)
827 scalar distFromBbSqr = 0;
829 if (!bb.contains(samplePoint))
831 distFromBbSqr = bb.nearest(samplePoint).distSqr(samplePoint);
834 pointsBbDistSqr.append
836 std::make_pair(distFromBbSqr,
count++)
842 pointsBbDistSqr.begin(),
843 pointsBbDistSqr.end(),
844 std::default_random_engine()
849 sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
851 typename Triangulation::Vertex_handle hint;
853 typename Triangulation::Locate_type lt;
860 Triangulation::number_of_vertices()
866 typename vectorPairPointIndex::const_iterator
p =
867 pointsBbDistSqr.begin();
868 p != pointsBbDistSqr.end();
872 const size_t checkInsertion = Triangulation::number_of_vertices();
874 const Vb& vert = *(
begin +
p->second);
875 const Point& pointToInsert = vert.point();
878 Cell_handle
c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
880 bool inserted =
false;
882 if (lt == Triangulation::VERTEX)
886 Vertex_handle nearV =
887 Triangulation::nearest_vertex(pointToInsert);
889 Pout<<
"Failed insertion, point already exists" <<
nl 890 <<
"Failed insertion : " << vert.
info()
891 <<
" nearest : " << nearV->info();
894 else if (lt == Triangulation::OUTSIDE_AFFINE_HULL)
897 <<
"Point is outside affine hull! pt = " << pointToInsert
900 else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
913 std::vector<Cell_handle> V;
914 typename Triangulation::Facet
f;
916 Triangulation::find_conflicts
920 CGAL::Oneset_iterator<typename Triangulation::Facet>(
f),
921 std::back_inserter(V)
924 for (
size_t i = 0; i < V.size(); ++i)
926 Cell_handle conflictingCell = V[i];
930 Triangulation::dimension() < 3
933 !Triangulation::is_infinite(conflictingCell)
935 conflictingCell->real()
936 || conflictingCell->hasFarPoint()
941 hint = Triangulation::insert_in_hole
959 if (checkInsertion != Triangulation::number_of_vertices() - 1)
963 Vertex_handle nearV =
964 Triangulation::nearest_vertex(pointToInsert);
966 Pout<<
"Failed insertion : " << vert.
info()
967 <<
" nearest : " << nearV->info();
972 hint->index() = vert.
index();
973 hint->type() = vert.
type();
scalar calculateLoadUnbalance() const
A HashTable with keys but without contents that is similar to std::unordered_set. ...
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
Ostream & indent(Ostream &os)
Indent stream.
pointFromPoint topoint(const Point &P)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
constexpr char nl
The newline '\n' character (0x0a)
labelPairHashSet rangeInsertReferredWithInfo(PointIterator begin, PointIterator end, bool printErrors=true)
Inserts points into the triangulation if the point is within.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
bool distribute(const boundBox &bb)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
List< labelList > labelListList
List of labelList.
Foam::scalar & targetCellSize()
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
#define forAll(list, i)
Loop across all elements in list.
void unset(List< bool > &bools, const labelUList &locations)
Unset the listed locations (assign 'false').
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Foam::InfoProxy< indexedVertex< Gt, Vb > > info() const noexcept
Return info proxy, used to print information to a stream.
void sort(UList< T > &list)
Sort the list.
void reset()
Clear the entire triangulation.
Istream and Ostream manipulators taking arguments.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
Pair< label > labelPair
A pair of labels.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Foam::tensor & alignment()
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setw(const int i)
List< label > labelList
A List of labels.
HashSet< labelPair, Foam::Hash< labelPair > > labelPairHashSet
A HashSet for a labelPair. The hashing is based on labelPair (FixedList) and is thus non-commutative...
void shuffle(UList< T > &list)
Randomise the list order.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
List< bool > boolList
A List of bools.
void sync(const boundBox &bb)
Refer vertices so that the processor interfaces are consistent.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
forAllConstIters(mixture.phases(), phase)
static constexpr const zero Zero
Global zero (0)