DistributedDelaunayMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "meshSearch.H"
28 #include "mapDistribute.H"
30 #include "pointConversion.H"
31 #include "indexedVertexEnum.H"
32 #include "IOmanip.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 
37 // * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
38 
39 template<class Triangulation>
42 (
43  const List<label>& toProc
44 )
45 {
46  // Determine send map
47  // ~~~~~~~~~~~~~~~~~~
48 
49  // 1. Count
50  labelList nSend(Pstream::nProcs(), 0);
51 
52  forAll(toProc, i)
53  {
54  label proci = toProc[i];
55 
56  nSend[proci]++;
57  }
58 
59 
60  // 2. Size sendMap
61  labelListList sendMap(Pstream::nProcs());
62 
63  forAll(nSend, proci)
64  {
65  sendMap[proci].setSize(nSend[proci]);
66 
67  nSend[proci] = 0;
68  }
69 
70  // 3. Fill sendMap
71  forAll(toProc, i)
72  {
73  label proci = toProc[i];
74 
75  sendMap[proci][nSend[proci]++] = i;
76  }
77 
78  // 4. Send over how many I need to receive
79  labelList recvSizes;
80  Pstream::exchangeSizes(sendMap, recvSizes);
81 
82 
83  // Determine receive map
84  // ~~~~~~~~~~~~~~~~~~~~~
85 
86  labelListList constructMap(Pstream::nProcs());
87 
88  // Local transfers first
89  constructMap[Pstream::myProcNo()] = identity
90  (
91  sendMap[Pstream::myProcNo()].size()
92  );
93 
94  label constructSize = constructMap[Pstream::myProcNo()].size();
95 
96  forAll(constructMap, proci)
97  {
98  if (proci != Pstream::myProcNo())
99  {
100  label nRecv = recvSizes[proci];
101 
102  constructMap[proci].setSize(nRecv);
103 
104  for (label i = 0; i < nRecv; i++)
105  {
106  constructMap[proci][i] = constructSize++;
107  }
108  }
109  }
110 
111  return autoPtr<mapDistribute>
112  (
113  new mapDistribute
114  (
115  constructSize,
116  sendMap.xfer(),
117  constructMap.xfer()
118  )
119  );
120 }
121 
122 
123 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
124 
125 template<class Triangulation>
127 (
128  const Time& runTime
129 )
130 :
131  DelaunayMesh<Triangulation>(runTime),
132  allBackgroundMeshBounds_()
133 {}
134 
135 
136 template<class Triangulation>
138 (
139  const Time& runTime,
140  const word& meshName
141 )
142 :
143  DelaunayMesh<Triangulation>(runTime, meshName),
144  allBackgroundMeshBounds_()
145 {}
146 
147 
148 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
149 
150 template<class Triangulation>
152 {}
153 
154 
155 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
156 
157 template<class Triangulation>
159 (
160  const boundBox& bb
161 )
162 {
163  allBackgroundMeshBounds_.reset(new List<boundBox>(Pstream::nProcs()));
164 
165  // Give the bounds of every processor to every other processor
166  allBackgroundMeshBounds_()[Pstream::myProcNo()] = bb;
167 
168  Pstream::gatherList(allBackgroundMeshBounds_());
169  Pstream::scatterList(allBackgroundMeshBounds_());
170 
171  return true;
172 }
173 
174 
175 template<class Triangulation>
177 (
178  const Vertex_handle& v
179 ) const
180 {
181  return isLocal(v->procIndex());
182 }
183 
184 
185 template<class Triangulation>
187 (
188  const label localProcIndex
189 ) const
190 {
191  return localProcIndex == Pstream::myProcNo();
192 }
193 
194 
195 template<class Triangulation>
197 (
198  const point& centre,
199  const scalar radiusSqr
200 ) const
201 {
202  DynamicList<label> toProc(Pstream::nProcs());
203 
204  forAll(allBackgroundMeshBounds_(), proci)
205  {
206  // Test against the bounding box of the processor
207  if
208  (
209  !isLocal(proci)
210  && allBackgroundMeshBounds_()[proci].overlaps(centre, radiusSqr)
211  )
212  {
213  toProc.append(proci);
214  }
215  }
216 
217  return toProc;
218 }
219 
220 
221 template<class Triangulation>
223 (
224  const Cell_handle& cit,
225  Map<labelList>& circumsphereOverlaps
226 ) const
227 {
228  const Foam::point& cc = cit->dual();
229 
230  const scalar crSqr = magSqr
231  (
232  cc - topoint(cit->vertex(0)->point())
233  );
234 
235  labelList circumsphereOverlap = overlapProcessors
236  (
237  cc,
238  sqr(1.01)*crSqr
239  );
240 
241  cit->cellIndex() = this->getNewCellIndex();
242 
243  if (!circumsphereOverlap.empty())
244  {
245  circumsphereOverlaps.insert(cit->cellIndex(), circumsphereOverlap);
246 
247  return true;
248  }
249 
250  return false;
251 }
252 
253 
254 template<class Triangulation>
256 (
257  Map<labelList>& circumsphereOverlaps
258 ) const
259 {
260  // Start by assuming that all the cells have no index
261  // If they do, they have already been visited so ignore them
262 
263  labelHashSet cellToCheck
264  (
265  Triangulation::number_of_finite_cells()
266  /Pstream::nProcs()
267  );
268 
269 // std::list<Cell_handle> infinite_cells;
270 // Triangulation::incident_cells
271 // (
272 // Triangulation::infinite_vertex(),
273 // std::back_inserter(infinite_cells)
274 // );
275 //
276 // for
277 // (
278 // typename std::list<Cell_handle>::iterator vcit
279 // = infinite_cells.begin();
280 // vcit != infinite_cells.end();
281 // ++vcit
282 // )
283 // {
284 // Cell_handle cit = *vcit;
285 //
286 // // Index of infinite vertex in this cell.
287 // label i = cit->index(Triangulation::infinite_vertex());
288 //
289 // Cell_handle c = cit->neighbor(i);
290 //
291 // if (c->unassigned())
292 // {
293 // c->cellIndex() = this->getNewCellIndex();
294 //
295 // if (checkProcBoundaryCell(c, circumsphereOverlaps))
296 // {
297 // cellToCheck.insert(c->cellIndex());
298 // }
299 // }
300 // }
301 //
302 //
303 // for
304 // (
305 // Finite_cells_iterator cit = Triangulation::finite_cells_begin();
306 // cit != Triangulation::finite_cells_end();
307 // ++cit
308 // )
309 // {
310 // if (cit->parallelDualVertex())
311 // {
312 // if (cit->unassigned())
313 // {
314 // if (checkProcBoundaryCell(cit, circumsphereOverlaps))
315 // {
316 // cellToCheck.insert(cit->cellIndex());
317 // }
318 // }
319 // }
320 // }
321 
322 
323  for
324  (
325  All_cells_iterator cit = Triangulation::all_cells_begin();
326  cit != Triangulation::all_cells_end();
327  ++cit
328  )
329  {
330  if (Triangulation::is_infinite(cit))
331  {
332  // Index of infinite vertex in this cell.
333  label i = cit->index(Triangulation::infinite_vertex());
334 
335  Cell_handle c = cit->neighbor(i);
336 
337  if (c->unassigned())
338  {
339  c->cellIndex() = this->getNewCellIndex();
340 
341  if (checkProcBoundaryCell(c, circumsphereOverlaps))
342  {
343  cellToCheck.insert(c->cellIndex());
344  }
345  }
346  }
347  else if (cit->parallelDualVertex())
348  {
349  if (cit->unassigned())
350  {
351  if (checkProcBoundaryCell(cit, circumsphereOverlaps))
352  {
353  cellToCheck.insert(cit->cellIndex());
354  }
355  }
356  }
357  }
358 
359  for
360  (
361  Finite_cells_iterator cit = Triangulation::finite_cells_begin();
362  cit != Triangulation::finite_cells_end();
363  ++cit
364  )
365  {
366  if (cellToCheck.found(cit->cellIndex()))
367  {
368  // Get the neighbours and check them
369  for (label adjCelli = 0; adjCelli < 4; ++adjCelli)
370  {
371  Cell_handle citNeighbor = cit->neighbor(adjCelli);
372 
373  // Ignore if has far point or previously visited
374  if
375  (
376  !citNeighbor->unassigned()
377  || !citNeighbor->internalOrBoundaryDualVertex()
378  || Triangulation::is_infinite(citNeighbor)
379  )
380  {
381  continue;
382  }
383 
384  if
385  (
386  checkProcBoundaryCell
387  (
388  citNeighbor,
389  circumsphereOverlaps
390  )
391  )
392  {
393  cellToCheck.insert(citNeighbor->cellIndex());
394  }
395  }
396 
397  cellToCheck.unset(cit->cellIndex());
398  }
399  }
400 }
401 
402 
403 template<class Triangulation>
405 (
406  const Map<labelList>& circumsphereOverlaps,
407  PtrList<labelPairHashSet>& referralVertices,
408  DynamicList<label>& targetProcessor,
409  DynamicList<Vb>& parallelInfluenceVertices
410 )
411 {
412  // Relying on the order of iteration of cells being the same as before
413  for
414  (
415  Finite_cells_iterator cit = Triangulation::finite_cells_begin();
416  cit != Triangulation::finite_cells_end();
417  ++cit
418  )
419  {
420  if (Triangulation::is_infinite(cit))
421  {
422  continue;
423  }
424 
426  circumsphereOverlaps.find(cit->cellIndex());
427 
428  // Pre-tested circumsphere potential influence
429  if (iter != circumsphereOverlaps.cend())
430  {
431  const labelList& citOverlaps = iter();
432 
433  forAll(citOverlaps, cOI)
434  {
435  label proci = citOverlaps[cOI];
436 
437  for (int i = 0; i < 4; i++)
438  {
439  Vertex_handle v = cit->vertex(i);
440 
441  if (v->farPoint())
442  {
443  continue;
444  }
445 
446  label vProcIndex = v->procIndex();
447  label vIndex = v->index();
448 
449  const labelPair procIndexPair(vProcIndex, vIndex);
450 
451  // Using the hashSet to ensure that each vertex is only
452  // referred once to each processor.
453  // Do not refer a vertex to its own processor.
454  if (vProcIndex != proci)
455  {
456  if (referralVertices[proci].insert(procIndexPair))
457  {
458  targetProcessor.append(proci);
459 
460  parallelInfluenceVertices.append
461  (
462  Vb
463  (
464  v->point(),
465  v->index(),
466  v->type(),
467  v->procIndex()
468  )
469  );
470 
471  parallelInfluenceVertices.last().targetCellSize() =
472  v->targetCellSize();
473  parallelInfluenceVertices.last().alignment() =
474  v->alignment();
475  }
476  }
477  }
478  }
479  }
480  }
481 }
482 
483 
484 template<class Triangulation>
486 (
487  const DynamicList<label>& targetProcessor,
488  DynamicList<Vb>& parallelVertices,
489  PtrList<labelPairHashSet>& referralVertices,
490  labelPairHashSet& receivedVertices
491 )
492 {
493  DynamicList<Vb> referredVertices(targetProcessor.size());
494 
495  const label preDistributionSize = parallelVertices.size();
496 
497  mapDistribute pointMap = buildMap(targetProcessor);
498 
499  // Make a copy of the original list.
500  DynamicList<Vb> originalParallelVertices(parallelVertices);
501 
502  pointMap.distribute(parallelVertices);
503 
504  for (label proci = 0; proci < Pstream::nProcs(); proci++)
505  {
506  const labelList& constructMap = pointMap.constructMap()[proci];
507 
508  if (constructMap.size())
509  {
510  forAll(constructMap, i)
511  {
512  const Vb& v = parallelVertices[constructMap[i]];
513 
514  if
515  (
517  && !receivedVertices.found(labelPair(v.procIndex(), v.index()))
518  )
519  {
520  referredVertices.append(v);
521 
522  receivedVertices.insert
523  (
524  labelPair(v.procIndex(), v.index())
525  );
526  }
527  }
528  }
529  }
530 
531  label preInsertionSize = Triangulation::number_of_vertices();
532 
533  labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
534  (
535  referredVertices.begin(),
536  referredVertices.end(),
537  true
538  );
539 
540  if (!pointsNotInserted.empty())
541  {
542  for
543  (
545  = pointsNotInserted.begin();
546  iter != pointsNotInserted.end();
547  ++iter
548  )
549  {
550  if (receivedVertices.found(iter.key()))
551  {
552  receivedVertices.erase(iter.key());
553  }
554  }
555  }
556 
557  boolList pointInserted(parallelVertices.size(), true);
558 
559  forAll(parallelVertices, vI)
560  {
561  const labelPair procIndexI
562  (
563  parallelVertices[vI].procIndex(),
564  parallelVertices[vI].index()
565  );
566 
567  if (pointsNotInserted.found(procIndexI))
568  {
569  pointInserted[vI] = false;
570  }
571  }
572 
573  pointMap.reverseDistribute(preDistributionSize, pointInserted);
574 
575  forAll(originalParallelVertices, vI)
576  {
577  const label procIndex = targetProcessor[vI];
578 
579  if (!pointInserted[vI])
580  {
581  if (referralVertices[procIndex].size())
582  {
583  if
584  (
585  !referralVertices[procIndex].unset
586  (
587  labelPair
588  (
589  originalParallelVertices[vI].procIndex(),
590  originalParallelVertices[vI].index()
591  )
592  )
593  )
594  {
595  Pout<< "*** not found "
596  << originalParallelVertices[vI].procIndex()
597  << " " << originalParallelVertices[vI].index() << endl;
598  }
599  }
600  }
601  }
602 
603  label postInsertionSize = Triangulation::number_of_vertices();
604 
605  reduce(preInsertionSize, sumOp<label>());
606  reduce(postInsertionSize, sumOp<label>());
607 
608  label nTotalToInsert = referredVertices.size();
609 
610  reduce(nTotalToInsert, sumOp<label>());
611 
612  if (preInsertionSize + nTotalToInsert != postInsertionSize)
613  {
614  label nNotInserted =
615  returnReduce(pointsNotInserted.size(), sumOp<label>());
616 
617  Info<< " Inserted = "
618  << setw(name(label(Triangulation::number_of_finite_cells())).size())
619  << nTotalToInsert - nNotInserted
620  << " / " << nTotalToInsert << endl;
621 
622  nTotalToInsert -= nNotInserted;
623  }
624  else
625  {
626  Info<< " Inserted = " << nTotalToInsert << endl;
627  }
628 
629  return nTotalToInsert;
630 }
631 
632 
633 template<class Triangulation>
635 (
636  const boundBox& bb,
637  PtrList<labelPairHashSet>& referralVertices,
638  labelPairHashSet& receivedVertices,
639  bool iterateReferral
640 )
641 {
642  if (!Pstream::parRun())
643  {
644  return;
645  }
646 
647  if (allBackgroundMeshBounds_.empty())
648  {
649  distributeBoundBoxes(bb);
650  }
651 
652  label nVerts = Triangulation::number_of_vertices();
653  label nCells = Triangulation::number_of_finite_cells();
654 
655  DynamicList<Vb> parallelInfluenceVertices(0.1*nVerts);
656  DynamicList<label> targetProcessor(0.1*nVerts);
657 
658  // Some of these values will not be used, i.e. for non-real cells
659  DynamicList<Foam::point> circumcentre(0.1*nVerts);
660  DynamicList<scalar> circumradiusSqr(0.1*nVerts);
661 
662  Map<labelList> circumsphereOverlaps(nCells);
663 
664  findProcessorBoundaryCells(circumsphereOverlaps);
665 
666  Info<< " Influences = "
667  << setw(name(nCells).size())
668  << returnReduce(circumsphereOverlaps.size(), sumOp<label>()) << " / "
669  << returnReduce(nCells, sumOp<label>());
670 
671  markVerticesToRefer
672  (
673  circumsphereOverlaps,
674  referralVertices,
675  targetProcessor,
676  parallelInfluenceVertices
677  );
678 
679  referVertices
680  (
681  targetProcessor,
682  parallelInfluenceVertices,
683  referralVertices,
684  receivedVertices
685  );
686 
687  if (iterateReferral)
688  {
689  label oldNReferred = 0;
690  label nIterations = 1;
691 
692  Info<< incrIndent << indent
693  << "Iteratively referring referred vertices..."
694  << endl;
695  do
696  {
697  Info<< indent << "Iteration " << nIterations++ << ":";
698 
699  circumsphereOverlaps.clear();
700  targetProcessor.clear();
701  parallelInfluenceVertices.clear();
702 
703  findProcessorBoundaryCells(circumsphereOverlaps);
704 
705  nCells = Triangulation::number_of_finite_cells();
706 
707  Info<< " Influences = "
708  << setw(name(nCells).size())
709  << returnReduce(circumsphereOverlaps.size(), sumOp<label>())
710  << " / "
711  << returnReduce(nCells, sumOp<label>());
712 
713  markVerticesToRefer
714  (
715  circumsphereOverlaps,
716  referralVertices,
717  targetProcessor,
718  parallelInfluenceVertices
719  );
720 
721  label nReferred = referVertices
722  (
723  targetProcessor,
724  parallelInfluenceVertices,
725  referralVertices,
726  receivedVertices
727  );
728 
729  if (nReferred == 0 || nReferred == oldNReferred)
730  {
731  break;
732  }
733 
734  oldNReferred = nReferred;
735 
736  } while (true);
737 
738  Info<< decrIndent;
739  }
740 }
741 
742 
743 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
744 
745 
746 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
747 
748 template<class Triangulation>
749 Foam::scalar
751 {
752  label nRealVertices = 0;
753 
754  for
755  (
756  Finite_vertices_iterator vit = Triangulation::finite_vertices_begin();
757  vit != Triangulation::finite_vertices_end();
758  ++vit
759  )
760  {
761  // Only store real vertices that are not feature vertices
762  if (vit->real() && !vit->featurePoint())
763  {
764  nRealVertices++;
765  }
766  }
767 
768  scalar globalNRealVertices = returnReduce
769  (
770  nRealVertices,
771  sumOp<label>()
772  );
773 
774  scalar unbalance = returnReduce
775  (
776  mag(1.0 - nRealVertices/(globalNRealVertices/Pstream::nProcs())),
777  maxOp<scalar>()
778  );
779 
780  Info<< " Processor unbalance " << unbalance << endl;
781 
782  return unbalance;
783 }
784 
785 
786 template<class Triangulation>
788 (
789  const boundBox& bb
790 )
791 {
793 
794  if (!Pstream::parRun())
795  {
796  return false;
797  }
798 
799  distributeBoundBoxes(bb);
800 
801  return true;
802 }
803 
804 
805 template<class Triangulation>
808 (
809  const backgroundMeshDecomposition& decomposition,
810  List<Foam::point>& points
811 )
812 {
813  if (!Pstream::parRun())
814  {
815  return autoPtr<mapDistribute>();
816  }
817 
818  distributeBoundBoxes(decomposition.procBounds());
819 
820  autoPtr<mapDistribute> mapDist = decomposition.distributePoints(points);
821 
822  return mapDist;
823 }
824 
825 
826 template<class Triangulation>
828 {
829  if (!Pstream::parRun())
830  {
831  return;
832  }
833 
834  if (allBackgroundMeshBounds_.empty())
835  {
836  distributeBoundBoxes(bb);
837  }
838 
839  const label nApproxReferred =
840  Triangulation::number_of_vertices()
841  /Pstream::nProcs();
842 
843  PtrList<labelPairHashSet> referralVertices(Pstream::nProcs());
844  forAll(referralVertices, proci)
845  {
846  if (!isLocal(proci))
847  {
848  referralVertices.set(proci, new labelPairHashSet(nApproxReferred));
849  }
850  }
851 
852  labelPairHashSet receivedVertices(nApproxReferred);
853 
854  sync
855  (
856  bb,
857  referralVertices,
858  receivedVertices,
859  true
860  );
861 }
862 
863 
864 template<class Triangulation>
865 template<class PointIterator>
868 (
869  PointIterator begin,
870  PointIterator end,
871  bool printErrors
872 )
873 {
874  const boundBox& bb = allBackgroundMeshBounds_()[Pstream::myProcNo()];
875 
876  typedef DynamicList
877  <
878  std::pair<scalar, label>
879  > vectorPairPointIndex;
880 
881  vectorPairPointIndex pointsBbDistSqr;
882 
883  label count = 0;
884  for (PointIterator it = begin; it != end; ++it)
885  {
886  const Foam::point samplePoint(topoint(it->point()));
887 
888  scalar distFromBbSqr = 0;
889 
890  if (!bb.contains(samplePoint))
891  {
892  const Foam::point nearestPoint = bb.nearest(samplePoint);
893 
894  distFromBbSqr = magSqr(nearestPoint - samplePoint);
895  }
896 
897  pointsBbDistSqr.append
898  (
899  std::make_pair(distFromBbSqr, count++)
900  );
901  }
902 
903  std::random_shuffle(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
904 
905  // Sort in ascending order by the distance of the point from the centre
906  // of the processor bounding box
907  sort(pointsBbDistSqr.begin(), pointsBbDistSqr.end());
908 
909  typename Triangulation::Vertex_handle hint;
910 
911  typename Triangulation::Locate_type lt;
912  int li, lj;
913 
914  label nNotInserted = 0;
915 
916  labelPairHashSet uninserted
917  (
918  Triangulation::number_of_vertices()
919  /Pstream::nProcs()
920  );
921 
922  for
923  (
924  typename vectorPairPointIndex::const_iterator p =
925  pointsBbDistSqr.begin();
926  p != pointsBbDistSqr.end();
927  ++p
928  )
929  {
930  const size_t checkInsertion = Triangulation::number_of_vertices();
931 
932  const Vb& vert = *(begin + p->second);
933  const Point& pointToInsert = vert.point();
934 
935  // Locate the point
936  Cell_handle c = Triangulation::locate(pointToInsert, lt, li, lj, hint);
937 
938  bool inserted = false;
939 
940  if (lt == Triangulation::VERTEX)
941  {
942  if (printErrors)
943  {
944  Vertex_handle nearV =
945  Triangulation::nearest_vertex(pointToInsert);
946 
947  Pout<< "Failed insertion, point already exists" << nl
948  << "Failed insertion : " << vert.info()
949  << " nearest : " << nearV->info();
950  }
951  }
952  else if (lt == Triangulation::OUTSIDE_AFFINE_HULL)
953  {
955  << "Point is outside affine hull! pt = " << pointToInsert
956  << endl;
957  }
958  else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
959  {
960  // TODO: Can this be optimised?
961  //
962  // Only want to insert if a connection is formed between
963  // pointToInsert and an internal or internal boundary point.
964  hint = Triangulation::insert(pointToInsert, c);
965  inserted = true;
966  }
967  else
968  {
969  // Get the cells that conflict with p in a vector V,
970  // and a facet on the boundary of this hole in f.
971  std::vector<Cell_handle> V;
972  typename Triangulation::Facet f;
973 
974  Triangulation::find_conflicts
975  (
976  pointToInsert,
977  c,
978  CGAL::Oneset_iterator<typename Triangulation::Facet>(f),
979  std::back_inserter(V)
980  );
981 
982  for (size_t i = 0; i < V.size(); ++i)
983  {
984  Cell_handle conflictingCell = V[i];
985 
986  if
987  (
988  Triangulation::dimension() < 3 // 2D triangulation
989  ||
990  (
991  !Triangulation::is_infinite(conflictingCell)
992  && (
993  conflictingCell->real()
994  || conflictingCell->hasFarPoint()
995  )
996  )
997  )
998  {
999  hint = Triangulation::insert_in_hole
1000  (
1001  pointToInsert,
1002  V.begin(),
1003  V.end(),
1004  f.first,
1005  f.second
1006  );
1007 
1008  inserted = true;
1009 
1010  break;
1011  }
1012  }
1013  }
1014 
1015  if (inserted)
1016  {
1017  if (checkInsertion != Triangulation::number_of_vertices() - 1)
1018  {
1019  if (printErrors)
1020  {
1021  Vertex_handle nearV =
1022  Triangulation::nearest_vertex(pointToInsert);
1023 
1024  Pout<< "Failed insertion : " << vert.info()
1025  << " nearest : " << nearV->info();
1026  }
1027  }
1028  else
1029  {
1030  hint->index() = vert.index();
1031  hint->type() = vert.type();
1032  hint->procIndex() = vert.procIndex();
1033  hint->targetCellSize() = vert.targetCellSize();
1034  hint->alignment() = vert.alignment();
1035  }
1036  }
1037  else
1038  {
1039  uninserted.insert(labelPair(vert.procIndex(), vert.index()));
1040  nNotInserted++;
1041  }
1042  }
1043 
1044  return uninserted;
1045 }
1046 
1047 
1048 // * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
1049 
1050 
1051 // * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * //
1052 
1053 
1054 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
1055 
1056 
1057 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
A HashTable with keys but without contents.
Definition: HashSet.H:59
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
HashTable< nil, Key, Hash >::const_iterator const_iterator
Definition: HashSet.H:67
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
static autoPtr< mapDistribute > buildMap(const List< label > &toProc)
Build a mapDistribute for the supplied destination processor data.
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
pointFromPoint topoint(const Point &P)
label getNewCellIndex() const
Create a new unique cell index and return.
Definition: DelaunayMeshI.H:63
Triangulation::Cell_handle Cell_handle
Definition: DelaunayMesh.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
labelPairHashSet rangeInsertReferredWithInfo(PointIterator begin, PointIterator end, bool printErrors=true)
Inserts points into the triangulation if the point is within.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:417
vertexType & type()
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Foam::label & index()
bool distribute(const boundBox &bb)
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
HashTable< T, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
Foam::scalar & targetCellSize()
Triangulation::Finite_cells_iterator Finite_cells_iterator
Definition: DelaunayMesh.H:75
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:51
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:138
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
static char meshName[]
Definition: globalFoam.H:7
void sort(UList< T > &)
Definition: UList.C:115
timeIndices insert(timeIndex, timeDirs[timeI].value())
void reset()
Clear the entire triangulation.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
List< label > labelList
A List of labels.
Definition: labelList.H:56
scalar calculateLoadUnbalance() const
iterator begin()
Iterator set to the beginning of the HashTable.
Definition: HashTableI.H:419
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
Istream and Ostream manipulators taking arguments.
HashSet< Pair< label >, FixedList< label, 2 >::Hash<> > labelPairHashSet
Definition: DelaunayMesh.H:83
static const char nl
Definition: Ostream.H:262
int procIndex() const
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
labelList f(nPoints)
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.
Definition: complex.C:47
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:399
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
Foam::InfoProxy< indexedVertex< Gt, Vb > > info() const
Return info proxy.
Triangulation::Vertex_handle Vertex_handle
Definition: DelaunayMesh.H:67
messageStream Info
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)
Definition: IOmanip.H:199
volScalarField & p
Triangulation::Point Point
Definition: DelaunayMesh.H:69
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
~DistributedDelaunayMesh()
Destructor.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
void sync(const boundBox &bb)
Refer vertices so that the processor interfaces are consistent.
Triangulation::Finite_vertices_iterator Finite_vertices_iterator
Definition: DelaunayMesh.H:73