41template<
class Triangulation>
45 const List<label>& toProc
56 label proci = toProc[i];
67 sendMap[proci].setSize(nSend[proci]);
75 label proci = toProc[i];
77 sendMap[proci][nSend[proci]++] = i;
82 Pstream::exchangeSizes(sendMap, recvSizes);
91 constructMap[Pstream::myProcNo()] =
identity
93 sendMap[Pstream::myProcNo()].size()
96 label constructSize = constructMap[Pstream::myProcNo()].size();
98 forAll(constructMap, proci)
100 if (proci != Pstream::myProcNo())
102 label nRecv = recvSizes[proci];
104 constructMap[proci].setSize(nRecv);
106 for (label i = 0; i < nRecv; i++)
108 constructMap[proci][i] = constructSize++;
113 return autoPtr<mapDistribute>::New
117 std::move(constructMap)
124template<
class Triangulation>
130 DelaunayMesh<Triangulation>(
runTime),
131 allBackgroundMeshBounds_()
135template<
class Triangulation>
142 DelaunayMesh<Triangulation>(
runTime, meshName),
143 allBackgroundMeshBounds_()
149template<
class Triangulation>
155 allBackgroundMeshBounds_.
reset(
new List<boundBox>(Pstream::nProcs()));
158 allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
160 Pstream::allGatherList(allBackgroundMeshBounds_());
166template<
class Triangulation>
169 const Vertex_handle& v
172 return isLocal(v->procIndex());
176template<
class Triangulation>
179 const label localProcIndex
182 return localProcIndex == Pstream::myProcNo();
186template<
class Triangulation>
190 const scalar radiusSqr
193 DynamicList<label> toProc(Pstream::nProcs());
195 forAll(allBackgroundMeshBounds_(), proci)
201 && allBackgroundMeshBounds_()[proci].overlaps(centre, radiusSqr)
204 toProc.append(proci);
212template<
class Triangulation>
215 const Cell_handle& cit,
216 Map<labelList>& circumsphereOverlaps
221 const scalar crSqr =
magSqr
223 cc -
topoint(cit->vertex(0)->point())
226 labelList circumsphereOverlap = overlapProcessors
232 cit->cellIndex() = this->getNewCellIndex();
234 if (!circumsphereOverlap.empty())
236 circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
245template<
class Triangulation>
248 Map<labelList>& circumsphereOverlaps
256 Triangulation::number_of_finite_cells()
316 All_cells_iterator cit = Triangulation::all_cells_begin();
317 cit != Triangulation::all_cells_end();
321 if (Triangulation::is_infinite(cit))
324 label i = cit->index(Triangulation::infinite_vertex());
326 Cell_handle
c = cit->neighbor(i);
330 c->cellIndex() = this->getNewCellIndex();
332 if (checkProcBoundaryCell(c, circumsphereOverlaps))
334 cellToCheck.insert(
c->cellIndex());
338 else if (cit->parallelDualVertex())
340 if (cit->unassigned())
342 if (checkProcBoundaryCell(cit, circumsphereOverlaps))
344 cellToCheck.insert(cit->cellIndex());
352 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
353 cit != Triangulation::finite_cells_end();
357 if (cellToCheck.found(cit->cellIndex()))
360 for (label adjCelli = 0; adjCelli < 4; ++adjCelli)
362 Cell_handle citNeighbor = cit->neighbor(adjCelli);
367 !citNeighbor->unassigned()
368 || !citNeighbor->internalOrBoundaryDualVertex()
369 || Triangulation::is_infinite(citNeighbor)
377 checkProcBoundaryCell
384 cellToCheck.insert(citNeighbor->cellIndex());
388 cellToCheck.unset(cit->cellIndex());
394template<
class Triangulation>
397 const Map<labelList>& circumsphereOverlaps,
398 PtrList<labelPairHashSet>& referralVertices,
399 DynamicList<label>& targetProcessor,
400 DynamicList<Vb>& parallelInfluenceVertices
406 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
407 cit != Triangulation::finite_cells_end();
411 if (Triangulation::is_infinite(cit))
416 const auto iter = circumsphereOverlaps.cfind(cit->cellIndex());
423 for (
const label proci : citOverlaps)
425 for (
int i = 0; i < 4; i++)
427 Vertex_handle v = cit->vertex(i);
435 label vIndex = v->index();
437 const labelPair procIndexPair(vProcIndex, vIndex);
442 if (vProcIndex != proci)
444 if (referralVertices[proci].
insert(procIndexPair))
446 targetProcessor.append(proci);
448 parallelInfluenceVertices.append
459 parallelInfluenceVertices.last().targetCellSize() =
461 parallelInfluenceVertices.last().alignment() =
472template<
class Triangulation>
475 const DynamicList<label>& targetProcessor,
476 DynamicList<Vb>& parallelVertices,
477 PtrList<labelPairHashSet>& referralVertices,
478 labelPairHashSet& receivedVertices
481 DynamicList<Vb> referredVertices(targetProcessor.size());
483 const label preDistributionSize = parallelVertices.size();
485 autoPtr<mapDistribute> pointMapPtr = buildMap(targetProcessor);
486 mapDistribute& pointMap = *pointMapPtr;
489 DynamicList<Vb> originalParallelVertices(parallelVertices);
491 pointMap.distribute(parallelVertices);
493 for (
const int proci : Pstream::allProcs())
495 const labelList& constructMap = pointMap.constructMap()[proci];
497 if (constructMap.size())
501 const Vb& v = parallelVertices[constructMap[i]];
509 referredVertices.append(v);
511 receivedVertices.insert
520 label preInsertionSize = Triangulation::number_of_vertices();
524 referredVertices.begin(),
525 referredVertices.end(),
529 if (!pointsNotInserted.empty())
533 if (receivedVertices.found(iter.key()))
535 receivedVertices.erase(iter.key());
540 boolList pointInserted(parallelVertices.size(),
true);
542 forAll(parallelVertices, vI)
546 parallelVertices[vI].procIndex(),
547 parallelVertices[vI].index()
550 if (pointsNotInserted.found(procIndexI))
552 pointInserted[vI] =
false;
556 pointMap.reverseDistribute(preDistributionSize, pointInserted);
558 forAll(originalParallelVertices, vI)
560 const label procIndex = targetProcessor[vI];
562 if (!pointInserted[vI])
564 if (referralVertices[procIndex].size())
568 !referralVertices[procIndex].unset
572 originalParallelVertices[vI].procIndex(),
573 originalParallelVertices[vI].index()
578 Pout<<
"*** not found "
579 << originalParallelVertices[vI].procIndex()
580 <<
" " << originalParallelVertices[vI].index() <<
endl;
586 label postInsertionSize = Triangulation::number_of_vertices();
588 reduce(preInsertionSize, sumOp<label>());
589 reduce(postInsertionSize, sumOp<label>());
591 label nTotalToInsert = referredVertices.size();
593 reduce(nTotalToInsert, sumOp<label>());
595 if (preInsertionSize + nTotalToInsert != postInsertionSize)
600 Info<<
" Inserted = "
601 <<
setw(
name(label(Triangulation::number_of_finite_cells())).size())
602 << nTotalToInsert - nNotInserted
603 <<
" / " << nTotalToInsert <<
endl;
605 nTotalToInsert -= nNotInserted;
609 Info<<
" Inserted = " << nTotalToInsert <<
endl;
612 return nTotalToInsert;
616template<
class Triangulation>
620 PtrList<labelPairHashSet>& referralVertices,
621 labelPairHashSet& receivedVertices,
625 if (!Pstream::parRun())
630 if (!allBackgroundMeshBounds_)
632 distributeBoundBoxes(bb);
635 label nVerts = Triangulation::number_of_vertices();
636 label nCells = Triangulation::number_of_finite_cells();
638 DynamicList<Vb> parallelInfluenceVertices(0.1*nVerts);
639 DynamicList<label> targetProcessor(0.1*nVerts);
642 DynamicList<Foam::point> circumcentre(0.1*nVerts);
643 DynamicList<scalar> circumradiusSqr(0.1*nVerts);
645 Map<labelList> circumsphereOverlaps(nCells);
647 findProcessorBoundaryCells(circumsphereOverlaps);
649 Info<<
" Influences = "
651 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>()) <<
" / "
656 circumsphereOverlaps,
659 parallelInfluenceVertices
665 parallelInfluenceVertices,
672 label oldNReferred = 0;
673 label nIterations = 1;
676 <<
"Iteratively referring referred vertices..."
680 Info<<
indent <<
"Iteration " << nIterations++ <<
":";
682 circumsphereOverlaps.clear();
683 targetProcessor.clear();
684 parallelInfluenceVertices.clear();
686 findProcessorBoundaryCells(circumsphereOverlaps);
688 nCells = Triangulation::number_of_finite_cells();
690 Info<<
" Influences = "
692 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>())
698 circumsphereOverlaps,
701 parallelInfluenceVertices
704 label nReferred = referVertices
707 parallelInfluenceVertices,
712 if (nReferred == 0 || nReferred == oldNReferred)
717 oldNReferred = nReferred;
731template<
class Triangulation>
735 label nRealVertices = 0;
739 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
740 vit != Triangulation::finite_vertices_end();
745 if (vit->real() && !vit->featurePoint())
759 mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
763 Info<<
" Processor unbalance " << unbalance <<
endl;
769template<
class Triangulation>
777 if (!Pstream::parRun())
782 distributeBoundBoxes(bb);
788template<
class Triangulation>
792 const backgroundMeshDecomposition& decomposition,
796 if (!Pstream::parRun())
801 distributeBoundBoxes(decomposition.procBounds());
803 return decomposition.distributePoints(
points);
807template<
class Triangulation>
810 if (!Pstream::parRun())
815 if (!allBackgroundMeshBounds_)
817 distributeBoundBoxes(bb);
820 const label nApproxReferred =
821 Triangulation::number_of_vertices()
824 PtrList<labelPairHashSet> referralVertices(Pstream::nProcs());
825 forAll(referralVertices, proci)
845template<
class Triangulation>
846template<
class Po
intIterator>
855 const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
859 std::pair<scalar, label>
860 > vectorPairPointIndex;
862 vectorPairPointIndex pointsBbDistSqr;
865 for (PointIterator it = begin; it !=
end; ++it)
869 scalar distFromBbSqr = 0;
871 if (!bb.contains(samplePoint))
873 const Foam::point nearestPoint = bb.nearest(samplePoint);
875 distFromBbSqr =
magSqr(nearestPoint - samplePoint);
878 pointsBbDistSqr.append
880 std::make_pair(distFromBbSqr, count++)
886 pointsBbDistSqr.begin(),
887 pointsBbDistSqr.end(),
888 std::default_random_engine()
893 sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
897 typename Triangulation::Locate_type lt;
900 label nNotInserted = 0;
904 Triangulation::number_of_vertices()
910 typename vectorPairPointIndex::const_iterator
p =
911 pointsBbDistSqr.begin();
912 p != pointsBbDistSqr.end();
916 const size_t checkInsertion = Triangulation::number_of_vertices();
918 const Vb& vert = *(
begin +
p->second);
919 const Point& pointToInsert = vert.point();
922 Cell_handle
c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
924 bool inserted =
false;
926 if (lt == Triangulation::VERTEX)
930 Vertex_handle nearV =
931 Triangulation::nearest_vertex(pointToInsert);
933 Pout<<
"Failed insertion, point already exists" <<
nl
934 <<
"Failed insertion : " << vert.
info()
935 <<
" nearest : " << nearV->info();
938 else if (lt == Triangulation::OUTSIDE_AFFINE_HULL)
941 <<
"Point is outside affine hull! pt = " << pointToInsert
944 else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
950 hint = Triangulation::insert(pointToInsert, c);
957 std::vector<Cell_handle> V;
960 Triangulation::find_conflicts
964 CGAL::Oneset_iterator<typename Triangulation::Facet>(
f),
965 std::back_inserter(V)
968 for (
size_t i = 0; i < V.size(); ++i)
970 Cell_handle conflictingCell = V[i];
974 Triangulation::dimension() < 3
977 !Triangulation::is_infinite(conflictingCell)
979 conflictingCell->real()
980 || conflictingCell->hasFarPoint()
985 hint = Triangulation::insert_in_hole
1003 if (checkInsertion != Triangulation::number_of_vertices() - 1)
1007 Vertex_handle nearV =
1008 Triangulation::nearest_vertex(pointToInsert);
1010 Pout<<
"Failed insertion : " << vert.
info()
1011 <<
" nearest : " << nearV->info();
1016 hint->index() = vert.
index();
1017 hint->type() = vert.
type();
Istream and Ostream manipulators taking arguments.
reduce(hasMovingMesh, orOp< bool >())
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Foam::scalar & targetCellSize()
Foam::InfoProxy< indexedVertex< Gt, Vb > > info() const
Info proxy, to print information to a stream.
Foam::tensor & alignment()
void reset()
Clear the entire triangulation.
scalar calculateLoadUnbalance() const
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
labelPairHashSet rangeInsertReferredWithInfo(PointIterator begin, PointIterator end, bool printErrors=true)
Inserts points into the triangulation if the point is within.
A HashTable with keys but without contents that is similar to std::unordered_set.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
CellSizeDelaunay::Vertex_handle Vertex_handle
void sync()
Do all: synchronise all IOFields and objectRegistry.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
#define WarningInFunction
Report a warning using Foam::Warning.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
const dimensionedScalar c
Speed of light in a vacuum.
Pair< label > labelPair
A pair of labels.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
pointFromPoint topoint(const Point &P)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setw(const int i)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & indent(Ostream &os)
Indent stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
void sort(UList< T > &list)
Sort the list.
List< bool > boolList
A List of bools.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
HashSet< labelPair, Foam::Hash< labelPair > > labelPairHashSet
A HashSet for a labelPair. The hashing is based on labelPair (FixedList) and is thus non-commutative.
constexpr char nl
The newline '\n' character (0x0a)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.