39 template<
class Triangulation>
43 const List<label>& toProc
54 label proci = toProc[i];
65 sendMap[proci].setSize(nSend[proci]);
73 label proci = toProc[i];
75 sendMap[proci][nSend[proci]++] = i;
80 Pstream::exchangeSizes(sendMap, recvSizes);
89 constructMap[Pstream::myProcNo()] =
identity 91 sendMap[Pstream::myProcNo()].size()
94 label constructSize = constructMap[Pstream::myProcNo()].size();
96 forAll(constructMap, proci)
98 if (proci != Pstream::myProcNo())
100 label nRecv = recvSizes[proci];
102 constructMap[proci].setSize(nRecv);
104 for (
label i = 0; i < nRecv; i++)
106 constructMap[proci][i] = constructSize++;
111 return autoPtr<mapDistribute>
125 template<
class Triangulation>
131 DelaunayMesh<Triangulation>(runTime),
132 allBackgroundMeshBounds_()
136 template<
class Triangulation>
143 DelaunayMesh<Triangulation>(runTime, meshName),
144 allBackgroundMeshBounds_()
150 template<
class Triangulation>
157 template<
class Triangulation>
163 allBackgroundMeshBounds_.
reset(
new List<boundBox>(Pstream::nProcs()));
166 allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
168 Pstream::gatherList(allBackgroundMeshBounds_());
169 Pstream::scatterList(allBackgroundMeshBounds_());
175 template<
class Triangulation>
178 const Vertex_handle& v
181 return isLocal(v->procIndex());
185 template<
class Triangulation>
188 const label localProcIndex
191 return localProcIndex == Pstream::myProcNo();
195 template<
class Triangulation>
199 const scalar radiusSqr
202 DynamicList<label> toProc(Pstream::nProcs());
204 forAll(allBackgroundMeshBounds_(), proci)
210 && allBackgroundMeshBounds_()[proci].overlaps(centre, radiusSqr)
213 toProc.append(proci);
221 template<
class Triangulation>
224 const Cell_handle& cit,
225 Map<labelList>& circumsphereOverlaps
230 const scalar crSqr =
magSqr 232 cc -
topoint(cit->vertex(0)->point())
235 labelList circumsphereOverlap = overlapProcessors
241 cit->cellIndex() = this->getNewCellIndex();
243 if (!circumsphereOverlap.empty())
245 circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
254 template<
class Triangulation>
257 Map<labelList>& circumsphereOverlaps
265 Triangulation::number_of_finite_cells()
325 All_cells_iterator cit = Triangulation::all_cells_begin();
326 cit != Triangulation::all_cells_end();
330 if (Triangulation::is_infinite(cit))
333 label i = cit->index(Triangulation::infinite_vertex());
335 Cell_handle c = cit->neighbor(i);
339 c->cellIndex() = this->getNewCellIndex();
341 if (checkProcBoundaryCell(c, circumsphereOverlaps))
343 cellToCheck.insert(c->cellIndex());
347 else if (cit->parallelDualVertex())
349 if (cit->unassigned())
351 if (checkProcBoundaryCell(cit, circumsphereOverlaps))
353 cellToCheck.insert(cit->cellIndex());
361 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
362 cit != Triangulation::finite_cells_end();
366 if (cellToCheck.found(cit->cellIndex()))
369 for (
label adjCelli = 0; adjCelli < 4; ++adjCelli)
371 Cell_handle citNeighbor = cit->neighbor(adjCelli);
376 !citNeighbor->unassigned()
377 || !citNeighbor->internalOrBoundaryDualVertex()
378 || Triangulation::is_infinite(citNeighbor)
386 checkProcBoundaryCell
393 cellToCheck.insert(citNeighbor->cellIndex());
397 cellToCheck.unset(cit->cellIndex());
403 template<
class Triangulation>
406 const Map<labelList>& circumsphereOverlaps,
407 PtrList<labelPairHashSet>& referralVertices,
408 DynamicList<label>& targetProcessor,
409 DynamicList<Vb>& parallelInfluenceVertices
415 Finite_cells_iterator cit = Triangulation::finite_cells_begin();
416 cit != Triangulation::finite_cells_end();
420 if (Triangulation::is_infinite(cit))
425 Map<labelList>::const_iterator iter =
426 circumsphereOverlaps.find(cit->cellIndex());
429 if (iter != circumsphereOverlaps.cend())
435 label proci = citOverlaps[cOI];
437 for (
int i = 0; i < 4; i++)
439 Vertex_handle v = cit->vertex(i);
446 label vProcIndex = v->procIndex();
447 label vIndex = v->index();
449 const labelPair procIndexPair(vProcIndex, vIndex);
454 if (vProcIndex != proci)
456 if (referralVertices[proci].
insert(procIndexPair))
458 targetProcessor.append(proci);
460 parallelInfluenceVertices.append
471 parallelInfluenceVertices.last().targetCellSize() =
473 parallelInfluenceVertices.last().alignment() =
484 template<
class Triangulation>
487 const DynamicList<label>& targetProcessor,
488 DynamicList<Vb>& parallelVertices,
489 PtrList<labelPairHashSet>& referralVertices,
490 labelPairHashSet& receivedVertices
493 DynamicList<Vb> referredVertices(targetProcessor.size());
495 const label preDistributionSize = parallelVertices.size();
497 mapDistribute pointMap = buildMap(targetProcessor);
500 DynamicList<Vb> originalParallelVertices(parallelVertices);
502 pointMap.distribute(parallelVertices);
504 for (
label proci = 0; proci < Pstream::nProcs(); proci++)
506 const labelList& constructMap = pointMap.constructMap()[proci];
508 if (constructMap.size())
512 const Vb& v = parallelVertices[constructMap[i]];
520 referredVertices.append(v);
522 receivedVertices.insert
531 label preInsertionSize = Triangulation::number_of_vertices();
533 labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
535 referredVertices.begin(),
536 referredVertices.end(),
540 if (!pointsNotInserted.empty())
544 typename labelPairHashSet::const_iterator iter
545 = pointsNotInserted.begin();
546 iter != pointsNotInserted.end();
550 if (receivedVertices.found(iter.key()))
552 receivedVertices.erase(iter.key());
557 boolList pointInserted(parallelVertices.size(),
true);
559 forAll(parallelVertices, vI)
563 parallelVertices[vI].procIndex(),
564 parallelVertices[vI].index()
567 if (pointsNotInserted.found(procIndexI))
569 pointInserted[vI] =
false;
573 pointMap.reverseDistribute(preDistributionSize, pointInserted);
575 forAll(originalParallelVertices, vI)
577 const label procIndex = targetProcessor[vI];
579 if (!pointInserted[vI])
581 if (referralVertices[procIndex].size())
585 !referralVertices[procIndex].unset
589 originalParallelVertices[vI].procIndex(),
590 originalParallelVertices[vI].index()
595 Pout<<
"*** not found " 596 << originalParallelVertices[vI].procIndex()
597 <<
" " << originalParallelVertices[vI].index() <<
endl;
603 label postInsertionSize = Triangulation::number_of_vertices();
605 reduce(preInsertionSize, sumOp<label>());
606 reduce(postInsertionSize, sumOp<label>());
608 label nTotalToInsert = referredVertices.size();
610 reduce(nTotalToInsert, sumOp<label>());
612 if (preInsertionSize + nTotalToInsert != postInsertionSize)
617 Info<<
" Inserted = " 618 <<
setw(
name(
label(Triangulation::number_of_finite_cells())).size())
619 << nTotalToInsert - nNotInserted
620 <<
" / " << nTotalToInsert <<
endl;
622 nTotalToInsert -= nNotInserted;
626 Info<<
" Inserted = " << nTotalToInsert <<
endl;
629 return nTotalToInsert;
633 template<
class Triangulation>
637 PtrList<labelPairHashSet>& referralVertices,
638 labelPairHashSet& receivedVertices,
642 if (!Pstream::parRun())
647 if (allBackgroundMeshBounds_.empty())
649 distributeBoundBoxes(bb);
652 label nVerts = Triangulation::number_of_vertices();
653 label nCells = Triangulation::number_of_finite_cells();
655 DynamicList<Vb> parallelInfluenceVertices(0.1*nVerts);
656 DynamicList<label> targetProcessor(0.1*nVerts);
659 DynamicList<Foam::point> circumcentre(0.1*nVerts);
660 DynamicList<scalar> circumradiusSqr(0.1*nVerts);
662 Map<labelList> circumsphereOverlaps(nCells);
664 findProcessorBoundaryCells(circumsphereOverlaps);
666 Info<<
" Influences = " 668 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>()) <<
" / " 673 circumsphereOverlaps,
676 parallelInfluenceVertices
682 parallelInfluenceVertices,
689 label oldNReferred = 0;
690 label nIterations = 1;
693 <<
"Iteratively referring referred vertices..." 697 Info<<
indent <<
"Iteration " << nIterations++ <<
":";
699 circumsphereOverlaps.clear();
700 targetProcessor.clear();
701 parallelInfluenceVertices.clear();
703 findProcessorBoundaryCells(circumsphereOverlaps);
705 nCells = Triangulation::number_of_finite_cells();
707 Info<<
" Influences = " 709 <<
returnReduce(circumsphereOverlaps.size(), sumOp<label>())
715 circumsphereOverlaps,
718 parallelInfluenceVertices
721 label nReferred = referVertices
724 parallelInfluenceVertices,
729 if (nReferred == 0 || nReferred == oldNReferred)
734 oldNReferred = nReferred;
748 template<
class Triangulation>
752 label nRealVertices = 0;
756 Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
757 vit != Triangulation::finite_vertices_end();
762 if (vit->real() && !vit->featurePoint())
776 mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
780 Info<<
" Processor unbalance " << unbalance <<
endl;
786 template<
class Triangulation>
794 if (!Pstream::parRun())
799 distributeBoundBoxes(bb);
805 template<
class Triangulation>
809 const backgroundMeshDecomposition& decomposition,
810 List<Foam::point>& points
813 if (!Pstream::parRun())
815 return autoPtr<mapDistribute>();
818 distributeBoundBoxes(decomposition.procBounds());
820 autoPtr<mapDistribute> mapDist = decomposition.distributePoints(points);
826 template<
class Triangulation>
829 if (!Pstream::parRun())
834 if (allBackgroundMeshBounds_.empty())
836 distributeBoundBoxes(bb);
839 const label nApproxReferred =
840 Triangulation::number_of_vertices()
843 PtrList<labelPairHashSet> referralVertices(Pstream::nProcs());
844 forAll(referralVertices, proci)
848 referralVertices.set(proci,
new labelPairHashSet(nApproxReferred));
852 labelPairHashSet receivedVertices(nApproxReferred);
864 template<
class Triangulation>
865 template<
class Po
intIterator>
874 const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
878 std::pair<scalar, label>
879 > vectorPairPointIndex;
881 vectorPairPointIndex pointsBbDistSqr;
884 for (PointIterator it = begin; it != end; ++it)
888 scalar distFromBbSqr = 0;
890 if (!bb.contains(samplePoint))
892 const Foam::point nearestPoint = bb.nearest(samplePoint);
894 distFromBbSqr =
magSqr(nearestPoint - samplePoint);
897 pointsBbDistSqr.append
899 std::make_pair(distFromBbSqr, count++)
903 std::random_shuffle(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
907 sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
909 typename Triangulation::Vertex_handle hint;
911 typename Triangulation::Locate_type lt;
914 label nNotInserted = 0;
916 labelPairHashSet uninserted
918 Triangulation::number_of_vertices()
924 typename vectorPairPointIndex::const_iterator
p =
925 pointsBbDistSqr.begin();
926 p != pointsBbDistSqr.end();
930 const size_t checkInsertion = Triangulation::number_of_vertices();
932 const Vb& vert = *(begin +
p->second);
933 const Point& pointToInsert = vert.point();
936 Cell_handle c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
938 bool inserted =
false;
940 if (lt == Triangulation::VERTEX)
944 Vertex_handle nearV =
945 Triangulation::nearest_vertex(pointToInsert);
947 Pout<<
"Failed insertion, point already exists" <<
nl 948 <<
"Failed insertion : " << vert.
info()
949 <<
" nearest : " << nearV->info();
952 else if (lt == Triangulation::OUTSIDE_AFFINE_HULL)
955 <<
"Point is outside affine hull! pt = " << pointToInsert
958 else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
971 std::vector<Cell_handle> V;
972 typename Triangulation::Facet
f;
974 Triangulation::find_conflicts
978 CGAL::Oneset_iterator<typename Triangulation::Facet>(f),
979 std::back_inserter(V)
982 for (
size_t i = 0; i < V.size(); ++i)
984 Cell_handle conflictingCell = V[i];
988 Triangulation::dimension() < 3
991 !Triangulation::is_infinite(conflictingCell)
993 conflictingCell->real()
994 || conflictingCell->hasFarPoint()
999 hint = Triangulation::insert_in_hole
1017 if (checkInsertion != Triangulation::number_of_vertices() - 1)
1021 Vertex_handle nearV =
1022 Triangulation::nearest_vertex(pointToInsert);
1024 Pout<<
"Failed insertion : " << vert.
info()
1025 <<
" nearest : " << nearV->info();
1030 hint->index() = vert.
index();
1031 hint->type() = vert.
type();
List< labelList > labelListList
A List of labelList.
scalar calculateLoadUnbalance() const
A HashTable with keys but without contents.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
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)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
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.
bool distribute(const boundBox &bb)
DistributedDelaunayMesh(const Time &runTime)
Construct from components.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Foam::scalar & targetCellSize()
void insert(const scalar, DynamicList< floatScalar > &)
Append scalar to given DynamicList.
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
List< bool > boolList
Bool container classes.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
void reset()
Clear the entire triangulation.
Pair< label > labelPair
Label pair.
List< label > labelList
A List of labels.
Foam::InfoProxy< indexedVertex< Gt, Vb > > info() const
Return info proxy.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Foam::tensor & alignment()
word name(const complex &)
Return a string representation of a complex.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Omanip< int > setw(const int i)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
~DistributedDelaunayMesh()
Destructor.
void sync(const boundBox &bb)
Refer vertices so that the processor interfaces are consistent.