conformalVoronoiMeshCalcDualMesh.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 
26 #include "conformalVoronoiMesh.H"
27 #include "motionSmoother.H"
29 #include "polyMeshGeometry.H"
30 #include "indexedCellChecks.H"
31 #include "OBJstream.H"
32 #include "indexedCellOps.H"
33 #include "ListOps.H"
34 #include "DelaunayMeshTools.H"
35 
36 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
37 
38 void Foam::conformalVoronoiMesh::calcDualMesh
39 (
40  pointField& points,
41  labelList& boundaryPts,
42  faceList& faces,
43  labelList& owner,
44  labelList& neighbour,
45  wordList& patchNames,
46  PtrList<dictionary>& patchDicts,
47  pointField& cellCentres,
48  labelList& cellToDelaunayVertex,
49  labelListList& patchToDelaunayVertex,
50  PackedBoolList& boundaryFacesToRemove
51 )
52 {
53  timeCheck("Start calcDualMesh");
54 
55  setVertexSizeAndAlignment();
56 
57  timeCheck("After setVertexSizeAndAlignment");
58 
59  indexDualVertices(points, boundaryPts);
60 
61  {
62  Info<< nl << "Merging identical points" << endl;
63 
64  // There is no guarantee that a merge of close points is no-risk
65  mergeIdenticalDualVertices(points, boundaryPts);
66  }
67 
68  // Final dual face and owner neighbour construction
69 
70  timeCheck("Before createFacesOwnerNeighbourAndPatches");
71 
72  createFacesOwnerNeighbourAndPatches
73  (
74  points,
75  faces,
76  owner,
77  neighbour,
78  patchNames,
79  patchDicts,
80  patchToDelaunayVertex, // from patch face to Delaunay vertex (slavePp)
81  boundaryFacesToRemove,
82  false
83  );
84 
85  // deferredCollapseFaceSet(owner, neighbour, deferredCollapseFaces);
86 
87  cellCentres = DelaunayMeshTools::allPoints(*this);
88 
89  cellToDelaunayVertex = removeUnusedCells(owner, neighbour);
90 
91  cellCentres = pointField(cellCentres, cellToDelaunayVertex);
92 
93  removeUnusedPoints(faces, points, boundaryPts);
94 
95  timeCheck("End of calcDualMesh");
96 }
97 
98 
99 void Foam::conformalVoronoiMesh::calcTetMesh
100 (
101  pointField& points,
102  labelList& pointToDelaunayVertex,
103  faceList& faces,
104  labelList& owner,
105  labelList& neighbour,
106  wordList& patchNames,
107  PtrList<dictionary>& patchDicts
108 )
109 {
110  labelList vertexMap(number_of_vertices());
111 
112  label vertI = 0;
113 
114  points.setSize(number_of_vertices());
115  pointToDelaunayVertex.setSize(number_of_vertices());
116 
117  for
118  (
119  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
120  vit != finite_vertices_end();
121  ++vit
122  )
123  {
124  if (vit->internalPoint() || vit->boundaryPoint())
125  {
126  vertexMap[vit->index()] = vertI;
127  points[vertI] = topoint(vit->point());
128  pointToDelaunayVertex[vertI] = vit->index();
129  vertI++;
130  }
131  }
132 
133  points.setSize(vertI);
134  pointToDelaunayVertex.setSize(vertI);
135 
136  label celli = 0;
137 
138  for
139  (
140  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
141  cit != finite_cells_end();
142  ++cit
143  )
144  {
145  if (cit->internalOrBoundaryDualVertex())
146  {
147  cit->cellIndex() = celli++;
148  }
149  else
150  {
151  cit->cellIndex() = Cb::ctFar;
152  }
153  }
154 
155  patchNames = geometryToConformTo_.patchNames();
156 
157  patchNames.setSize(patchNames.size() + 1);
158 
159  patchNames[patchNames.size() - 1] = "foamyHexMesh_defaultPatch";
160 
161  label nPatches = patchNames.size();
162 
163  List<DynamicList<face>> patchFaces(nPatches, DynamicList<face>(0));
164 
165  List<DynamicList<label>> patchOwners(nPatches, DynamicList<label>(0));
166 
167  faces.setSize(number_of_finite_facets());
168 
169  owner.setSize(number_of_finite_facets());
170 
171  neighbour.setSize(number_of_finite_facets());
172 
173  label facei = 0;
174 
175  labelList verticesOnTriFace(3, label(-1));
176 
177  face newFace(verticesOnTriFace);
178 
179  for
180  (
181  Delaunay::Finite_facets_iterator fit = finite_facets_begin();
182  fit != finite_facets_end();
183  ++fit
184  )
185  {
186  const Cell_handle c1(fit->first);
187  const label oppositeVertex = fit->second;
188  const Cell_handle c2(c1->neighbor(oppositeVertex));
189 
190  if (c1->hasFarPoint() && c2->hasFarPoint())
191  {
192  // Both tets are outside, skip
193  continue;
194  }
195 
196  label c1I = c1->cellIndex();
197  label c2I = c2->cellIndex();
198 
199  label ownerCell = -1;
200  label neighbourCell = -1;
201 
202  for (label i = 0; i < 3; i++)
203  {
204  verticesOnTriFace[i] = vertexMap
205  [
206  c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
207  ];
208  }
209 
210  newFace = face(verticesOnTriFace);
211 
212  if (c1->hasFarPoint() || c2->hasFarPoint())
213  {
214  // Boundary face...
215  if (c1->hasFarPoint())
216  {
217  //... with c1 outside
218  ownerCell = c2I;
219  }
220  else
221  {
222  // ... with c2 outside
223  ownerCell = c1I;
224 
225  reverse(newFace);
226  }
227 
228  label patchIndex = geometryToConformTo_.findPatch
229  (
230  newFace.centre(points)
231  );
232 
233  if (patchIndex == -1)
234  {
235  patchIndex = patchNames.size() - 1;
236 
238  << newFace.centre(points) << nl
239  << "did not find a surface patch. Adding to "
240  << patchNames[patchIndex]
241  << endl;
242  }
243 
244  patchFaces[patchIndex].append(newFace);
245  patchOwners[patchIndex].append(ownerCell);
246  }
247  else
248  {
249  // Internal face...
250  if (c1I < c2I)
251  {
252  // ...with c1 as the ownerCell
253  ownerCell = c1I;
254  neighbourCell = c2I;
255 
256  reverse(newFace);
257  }
258  else
259  {
260  // ...with c2 as the ownerCell
261  ownerCell = c2I;
262  neighbourCell = c1I;
263  }
264 
265  faces[facei] = newFace;
266  owner[facei] = ownerCell;
267  neighbour[facei] = neighbourCell;
268  facei++;
269  }
270  }
271 
272  label nInternalFaces = facei;
273 
274  faces.setSize(nInternalFaces);
275  owner.setSize(nInternalFaces);
276  neighbour.setSize(nInternalFaces);
277 
278  sortFaces(faces, owner, neighbour);
279 
280 // PackedBoolList boundaryFacesToRemove;
281 // List<DynamicList<bool>> indirectPatchFace;
282 //
283 // addPatches
284 // (
285 // nInternalFaces,
286 // faces,
287 // owner,
288 // patchDicts,
289 // boundaryFacesToRemove,
290 // patchFaces,
291 // patchOwners,
292 // indirectPatchFace
293 // );
294 }
295 
296 
297 void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
298 (
299  const pointField& pts,
300  labelList& boundaryPts
301 )
302 {
303  // Assess close points to be merged
304 
305  label nPtsMerged = 0;
306  label nPtsMergedSum = 0;
307 
308  do
309  {
310  Map<label> dualPtIndexMap;
311 
312  nPtsMerged = mergeIdenticalDualVertices
313  (
314  pts,
315  dualPtIndexMap
316  );
317 
318  reindexDualVertices(dualPtIndexMap, boundaryPts);
319 
320  reduce(nPtsMerged, sumOp<label>());
321 
322  nPtsMergedSum += nPtsMerged;
323 
324  } while (nPtsMerged > 0);
325 
326  if (nPtsMergedSum > 0)
327  {
328  Info<< " Merged " << nPtsMergedSum << " points " << endl;
329  }
330 }
331 
332 
333 Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
334 (
335  const pointField& pts,
336  Map<label>& dualPtIndexMap
337 ) const
338 {
339  label nPtsMerged = 0;
340 
341  for
342  (
343  Delaunay::Finite_facets_iterator fit = finite_facets_begin();
344  fit != finite_facets_end();
345  ++fit
346  )
347  {
348  const Cell_handle c1(fit->first);
349  const label oppositeVertex = fit->second;
350  const Cell_handle c2(c1->neighbor(oppositeVertex));
351 
352  if (is_infinite(c1) || is_infinite(c2))
353  {
354  continue;
355  }
356 
357  label& c1I = c1->cellIndex();
358  label& c2I = c2->cellIndex();
359 
360  if ((c1I != c2I) && !c1->hasFarPoint() && !c2->hasFarPoint())
361  {
362  const Foam::point& p1 = pts[c1I];
363  const Foam::point& p2 = pts[c2I];
364 
365  if (p1 == p2)
366  {
367 // if (c1->parallelDualVertex() || c2->parallelDualVertex())
368 // {
369 // if (c1->vertexLowestProc() < c2->vertexLowestProc())
370 // {
371 // dualPtIndexMap.insert(c1I, c1I);
372 // dualPtIndexMap.insert(c2I, c1I);
373 // }
374 // else
375 // {
376 // dualPtIndexMap.insert(c1I, c2I);
377 // dualPtIndexMap.insert(c2I, c2I);
378 // }
379 // }
380  if (c1I < c2I)
381  {
382  dualPtIndexMap.insert(c1I, c1I);
383  dualPtIndexMap.insert(c2I, c1I);
384  }
385  else
386  {
387  dualPtIndexMap.insert(c1I, c2I);
388  dualPtIndexMap.insert(c2I, c2I);
389  }
390 
391  nPtsMerged++;
392  }
393  }
394  }
395 
396  if (debug)
397  {
398  Info<< "mergeIdenticalDualVertices:" << endl
399  << " zero-length edges : "
400  << returnReduce(nPtsMerged, sumOp<label>()) << endl
401  << endl;
402  }
403 
404  return nPtsMerged;
405 }
406 
407 
408 //void Foam::conformalVoronoiMesh::smoothSurface
409 //(
410 // pointField& pts,
411 // const labelList& boundaryPts
412 //)
413 //{
414 // label nCollapsedFaces = 0;
415 //
416 // label iterI = 0;
417 //
418 // do
419 // {
420 // Map<label> dualPtIndexMap;
421 //
422 // nCollapsedFaces = smoothSurfaceDualFaces
423 // (
424 // pts,
425 // boundaryPts,
426 // dualPtIndexMap
427 // );
428 //
429 // reduce(nCollapsedFaces, sumOp<label>());
430 //
431 // reindexDualVertices(dualPtIndexMap);
432 //
433 // mergeIdenticalDualVertices(pts, boundaryPts);
434 //
435 // if (nCollapsedFaces > 0)
436 // {
437 // Info<< " Collapsed " << nCollapsedFaces << " boundary faces"
438 // << endl;
439 // }
440 //
441 // if (++iterI > foamyHexMeshControls().maxCollapseIterations())
442 // {
443 // Info<< " maxCollapseIterations reached, stopping collapse"
444 // << endl;
445 //
446 // break;
447 // }
448 //
449 // } while (nCollapsedFaces > 0);
450 //
451 // // Force all points of boundary faces to be on the surface
502 //}
503 //
504 //
505 //Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
506 //(
507 // pointField& pts,
508 // const labelList& boundaryPts,
509 // Map<label>& dualPtIndexMap
510 //) const
511 //{
512 // label nCollapsedFaces = 0;
513 //
514 // const scalar cosPerpendicularToleranceAngle = cos
515 // (
516 // degToRad(foamyHexMeshControls().surfaceStepFaceAngle())
517 // );
518 //
519 // for
520 // (
521 // Delaunay::Finite_edges_iterator eit = finite_edges_begin();
522 // eit != finite_edges_end();
523 // ++eit
524 // )
525 // {
526 // Cell_circulator ccStart = incident_cells(*eit);
527 // Cell_circulator cc = ccStart;
528 //
529 // bool skipFace = false;
530 //
531 // do
532 // {
533 // if (dualPtIndexMap.found(cc->cellIndex()))
534 // {
535 // // One of the points of this face has already been
536 // // collapsed this sweep, leave for next sweep
537 //
538 // skipFace = true;
539 //
540 // break;
541 // }
542 //
543 // } while (++cc != ccStart);
544 //
545 // if (skipFace)
546 // {
547 // continue;
548 // }
549 //
550 // if (isBoundaryDualFace(eit))
551 // {
552 // face dualFace = buildDualFace(eit);
553 //
554 // if (dualFace.size() < 3)
555 // {
556 // // This face has been collapsed already
557 // continue;
558 // }
559 //
560 // label maxFC = maxFilterCount(eit);
561 //
562 // if (maxFC > foamyHexMeshControls().filterCountSkipThreshold())
563 // {
564 // // A vertex on this face has been limited too many
565 // // times, skip
566 // continue;
567 // }
568 //
569 // Cell_handle c = eit->first;
570 // Vertex_handle vA = c->vertex(eit->second);
571 // Vertex_handle vB = c->vertex(eit->third);
572 //
573 // if
574 // (
575 // vA->internalBoundaryPoint() && vA->surfacePoint()
576 // && vB->externalBoundaryPoint() && vB->surfacePoint()
577 // )
578 // {
579 // if (vA->index() == vB->index() - 1)
580 // {
581 // continue;
582 // }
583 // }
584 // else if
585 // (
586 // vA->externalBoundaryPoint() && vA->surfacePoint()
587 // && vB->internalBoundaryPoint() && vB->surfacePoint()
588 // )
589 // {
590 // if (vA->index() == vB->index() + 1)
591 // {
592 // continue;
593 // }
594 // }
639 //
640 //
641 // if ((faceNormal & surfaceNormal) < cosPerpendicularToleranceAngle)
642 // {
643 // scalar targetFaceSize = averageAnyCellSize(vA, vB);
644 //
645 // // Selecting faces to collapse based on angle to
646 // // surface, so set collapseSizeLimitCoeff to GREAT to
647 // // allow collapse of all faces
648 //
649 // faceCollapseMode mode = collapseFace
650 // (
651 // dualFace,
652 // pts,
653 // boundaryPts,
654 // dualPtIndexMap,
655 // targetFaceSize,
656 // GREAT,
657 // maxFC
658 // );
659 //
660 // if (mode == fcmPoint || mode == fcmEdge)
661 // {
662 // nCollapsedFaces++;
663 // }
664 // }
665 // }
666 // }
667 //
668 // return nCollapsedFaces;
669 //}
670 
671 
672 void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
673 (
674  labelList& owner,
675  labelList& neighbour,
676  const HashSet<labelPair, labelPair::Hash<>>& deferredCollapseFaces
677 ) const
678 {
679  DynamicList<label> faceLabels;
680 
681  forAll(neighbour, nI)
682  {
683  if (deferredCollapseFaces.found(Pair<label>(owner[nI], neighbour[nI])))
684  {
685  faceLabels.append(nI);
686  }
687  }
688 
689  Pout<< "facesToCollapse" << nl << faceLabels << endl;
690 }
691 
692 
694 Foam::conformalVoronoiMesh::createPolyMeshFromPoints
695 (
696  const pointField& pts
697 ) const
698 {
699  faceList faces;
700  labelList owner;
701  labelList neighbour;
703  PtrList<dictionary> patchDicts;
704  pointField cellCentres;
705  labelListList patchToDelaunayVertex;
706  PackedBoolList boundaryFacesToRemove;
707 
708  timeCheck("Start of checkPolyMeshQuality");
709 
710  Info<< nl << "Creating polyMesh to assess quality" << endl;
711 
712  createFacesOwnerNeighbourAndPatches
713  (
714  pts,
715  faces,
716  owner,
717  neighbour,
718  patchNames,
719  patchDicts,
720  patchToDelaunayVertex,
721  boundaryFacesToRemove,
722  false
723  );
724 
725  cellCentres = DelaunayMeshTools::allPoints(*this);
726 
727  labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
728  cellCentres = pointField(cellCentres, cellToDelaunayVertex);
729 
730  autoPtr<polyMesh> meshPtr
731  (
732  new polyMesh
733  (
734  IOobject
735  (
736  "foamyHexMesh_temporary",
737  runTime_.timeName(),
738  runTime_,
741  ),
742  xferCopy(pts),
743  xferMove(faces),
744  xferMove(owner),
745  xferMove(neighbour)
746  )
747  );
748 
749  polyMesh& pMesh = meshPtr();
750 
751  List<polyPatch*> patches(patchNames.size());
752 
753  label nValidPatches = 0;
754 
755  forAll(patches, p)
756  {
757  label totalPatchSize = readLabel(patchDicts[p].lookup("nFaces"));
758 
759  if
760  (
761  patchDicts.set(p)
762  && word(patchDicts[p].lookup("type")) == processorPolyPatch::typeName
763  )
764  {
765  // Do not create empty processor patches
766  if (totalPatchSize > 0)
767  {
768  patchDicts[p].set("transform", "coincidentFullMatch");
769 
770  patches[nValidPatches] = new processorPolyPatch
771  (
772  patchNames[p],
773  patchDicts[p],
774  nValidPatches,
775  pMesh.boundaryMesh(),
776  processorPolyPatch::typeName
777  );
778 
779  nValidPatches++;
780  }
781  }
782  else
783  {
784  // Check that the patch is not empty on every processor
785  reduce(totalPatchSize, sumOp<label>());
786 
787  if (totalPatchSize > 0)
788  {
789  patches[nValidPatches] = polyPatch::New
790  (
791  patchNames[p],
792  patchDicts[p],
793  nValidPatches,
794  pMesh.boundaryMesh()
795  ).ptr();
796 
797  nValidPatches++;
798  }
799  }
800  }
801 
802  patches.setSize(nValidPatches);
803 
804  pMesh.addPatches(patches);
805 
806  return meshPtr;
807 }
808 
809 
810 void Foam::conformalVoronoiMesh::checkCellSizing()
811 {
812  Info<< "Checking cell sizes..."<< endl;
813 
814  timeCheck("Start of Cell Sizing");
815 
816  labelList boundaryPts(number_of_finite_cells(), internal);
817  pointField ptsField;
818 
819  indexDualVertices(ptsField, boundaryPts);
820 
821  // Merge close dual vertices.
822  mergeIdenticalDualVertices(ptsField, boundaryPts);
823 
824  autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
825  const polyMesh& pMesh = meshPtr();
826 
827  //pMesh.write();
828 
829  // Find cells with poor quality
830  DynamicList<label> checkFaces(identity(pMesh.nFaces()));
831  labelHashSet wrongFaces(pMesh.nFaces()/100);
832 
833  Info<< "Running checkMesh on mesh with " << pMesh.nCells()
834  << " cells "<< endl;
835 
836  const dictionary& dict
838 
839  const dictionary& meshQualityDict
840  = dict.subDict("meshQualityControls");
841 
842  const scalar maxNonOrtho =
843  readScalar(meshQualityDict.lookup("maxNonOrtho", true));
844 
845  label nWrongFaces = 0;
846 
847  if (maxNonOrtho < 180.0 - SMALL)
848  {
850  (
851  false,
852  maxNonOrtho,
853  pMesh,
854  pMesh.cellCentres(),
855  pMesh.faceAreas(),
856  checkFaces,
857  List<labelPair>(),
858  &wrongFaces
859  );
860 
861  label nNonOrthogonal = returnReduce(wrongFaces.size(), sumOp<label>());
862 
863  Info<< " non-orthogonality > " << maxNonOrtho
864  << " degrees : " << nNonOrthogonal << endl;
865 
866  nWrongFaces += nNonOrthogonal;
867  }
868 
869  labelHashSet protrudingCells = findOffsetPatchFaces(pMesh, 0.25);
870 
871  label nProtrudingCells = protrudingCells.size();
872 
873  Info<< " protruding/intruding cells : " << nProtrudingCells << endl;
874 
875  nWrongFaces += nProtrudingCells;
876 
877 // motionSmoother::checkMesh
878 // (
879 // false,
880 // pMesh,
881 // meshQualityDict,
882 // checkFaces,
883 // wrongFaces
884 // );
885 
886  Info<< " Found total of " << nWrongFaces << " bad faces" << endl;
887 
888  {
889  labelHashSet cellsToResizeMap(pMesh.nFaces()/100);
890 
891  // Find cells that are attached to the faces in wrongFaces.
892  forAllConstIter(labelHashSet, wrongFaces, iter)
893  {
894  const label faceOwner = pMesh.faceOwner()[iter.key()];
895  const label faceNeighbour = pMesh.faceNeighbour()[iter.key()];
896 
897  if (!cellsToResizeMap.found(faceOwner))
898  {
899  cellsToResizeMap.insert(faceOwner);
900  }
901 
902  if (!cellsToResizeMap.found(faceNeighbour))
903  {
904  cellsToResizeMap.insert(faceNeighbour);
905  }
906  }
907 
908  cellsToResizeMap += protrudingCells;
909 
910  pointField cellsToResize(cellsToResizeMap.size());
911 
912  label count = 0;
913  for (label celli = 0; celli < pMesh.nCells(); ++celli)
914  {
915  if (cellsToResizeMap.found(celli))
916  {
917  cellsToResize[count++] = pMesh.cellCentres()[celli];
918  }
919  }
920 
921  Info<< " DISABLED: Automatically re-sizing " << cellsToResize.size()
922  << " cells that are attached to the bad faces: " << endl;
923 
924  //cellSizeControl_.setCellSizes(cellsToResize);
925  }
926 
927  timeCheck("End of Cell Sizing");
928 
929  Info<< "Finished checking cell sizes"<< endl;
930 }
931 
932 
933 Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
934 (
935  const polyMesh& mesh,
936  const scalar allowedOffset
937 ) const
938 {
939  timeCheck("Start findRemainingProtrusionSet");
940 
941  const polyBoundaryMesh& patches = mesh.boundaryMesh();
942 
943  cellSet offsetBoundaryCells
944  (
945  mesh,
946  "foamyHexMesh_protrudingCells",
947  mesh.nCells()/1000
948  );
949 
950  forAll(patches, patchi)
951  {
952  const polyPatch& patch = patches[patchi];
953 
954  const faceList& localFaces = patch.localFaces();
955  const pointField& localPoints = patch.localPoints();
956 
957  const labelList& fCell = patch.faceCells();
958 
959  forAll(localFaces, pLFI)
960  {
961  const face& f = localFaces[pLFI];
962 
963  const Foam::point& faceCentre = f.centre(localPoints);
964 
965  const scalar targetSize = targetCellSize(faceCentre);
966 
967  pointIndexHit pHit;
968  label surfHit = -1;
969 
970  geometryToConformTo_.findSurfaceNearest
971  (
972  faceCentre,
973  sqr(targetSize),
974  pHit,
975  surfHit
976  );
977 
978  if
979  (
980  pHit.hit()
981  && (mag(pHit.hitPoint() - faceCentre) > allowedOffset*targetSize)
982  )
983  {
984  offsetBoundaryCells.insert(fCell[pLFI]);
985  }
986  }
987  }
988 
989  if (foamyHexMeshControls().objOutput())
990  {
991  offsetBoundaryCells.write();
992  }
993 
994  return offsetBoundaryCells;
995 }
996 
997 
998 Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
999 (
1000  const pointField& pts
1001 ) const
1002 {
1003  autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(pts);
1004  polyMesh& pMesh = meshPtr();
1005 
1006  timeCheck("polyMesh created, checking quality");
1007 
1008  labelHashSet wrongFaces(pMesh.nFaces()/100);
1009 
1010  DynamicList<label> checkFaces(pMesh.nFaces());
1011 
1012  const vectorField& fAreas = pMesh.faceAreas();
1013 
1014  scalar faceAreaLimit = SMALL;
1015 
1016  forAll(fAreas, fI)
1017  {
1018  if (mag(fAreas[fI]) > faceAreaLimit)
1019  {
1020  checkFaces.append(fI);
1021  }
1022  }
1023 
1024  Info<< nl << "Excluding "
1025  << returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
1026  << " faces from check, < " << faceAreaLimit << " area" << endl;
1027 
1028  const dictionary& dict
1030 
1031  const dictionary& meshQualityDict
1032  = dict.subDict("meshQualityControls");
1033 
1035  (
1036  false,
1037  pMesh,
1038  meshQualityDict,
1039  checkFaces,
1040  wrongFaces
1041  );
1042 
1043  {
1044  // Check for cells with more than 1 but fewer than 4 faces
1045  label nInvalidPolyhedra = 0;
1046 
1047  const cellList& cells = pMesh.cells();
1048 
1049  forAll(cells, cI)
1050  {
1051  if (cells[cI].size() < 4 && cells[cI].size() > 0)
1052  {
1053  // Pout<< "cell " << cI << " " << cells[cI]
1054  // << " has " << cells[cI].size() << " faces."
1055  // << endl;
1056 
1057  nInvalidPolyhedra++;
1058 
1059  forAll(cells[cI], cFI)
1060  {
1061  wrongFaces.insert(cells[cI][cFI]);
1062  }
1063  }
1064  }
1065 
1066  Info<< " cells with more than 1 but fewer than 4 faces : "
1067  << returnReduce(nInvalidPolyhedra, sumOp<label>())
1068  << endl;
1069 
1070  // Check for cells with one internal face only
1071 
1072  labelList nInternalFaces(pMesh.nCells(), label(0));
1073 
1074  for (label fI = 0; fI < pMesh.nInternalFaces(); fI++)
1075  {
1076  nInternalFaces[pMesh.faceOwner()[fI]]++;
1077  nInternalFaces[pMesh.faceNeighbour()[fI]]++;
1078  }
1079 
1080  const polyBoundaryMesh& patches = pMesh.boundaryMesh();
1081 
1082  forAll(patches, patchi)
1083  {
1084  if (patches[patchi].coupled())
1085  {
1086  const labelUList& owners = patches[patchi].faceCells();
1087 
1088  forAll(owners, i)
1089  {
1090  nInternalFaces[owners[i]]++;
1091  }
1092  }
1093  }
1094 
1095  label oneInternalFaceCells = 0;
1096 
1097  forAll(nInternalFaces, cI)
1098  {
1099  if (nInternalFaces[cI] <= 1)
1100  {
1101  oneInternalFaceCells++;
1102 
1103  forAll(cells[cI], cFI)
1104  {
1105  wrongFaces.insert(cells[cI][cFI]);
1106  }
1107  }
1108  }
1109 
1110  Info<< " cells with with zero or one non-boundary face : "
1111  << returnReduce(oneInternalFaceCells, sumOp<label>())
1112  << endl;
1113  }
1114 
1115 
1116  PackedBoolList ptToBeLimited(pts.size(), false);
1117 
1118  forAllConstIter(labelHashSet, wrongFaces, iter)
1119  {
1120  const face f = pMesh.faces()[iter.key()];
1121 
1122  forAll(f, fPtI)
1123  {
1124  ptToBeLimited[f[fPtI]] = true;
1125  }
1126  }
1127 
1128  // // Limit connected cells
1129 
1130  // labelHashSet limitCells(pMesh.nCells()/100);
1131 
1132  // const labelListList& ptCells = pMesh.pointCells();
1133 
1134  // forAllConstIter(labelHashSet, wrongFaces, iter)
1135  // {
1136  // const face f = pMesh.faces()[iter.key()];
1137 
1138  // forAll(f, fPtI)
1139  // {
1140  // label ptI = f[fPtI];
1141 
1142  // const labelList& pC = ptCells[ptI];
1143 
1144  // forAll(pC, pCI)
1145  // {
1146  // limitCells.insert(pC[pCI]);
1147  // }
1148  // }
1149  // }
1150 
1151  // const labelListList& cellPts = pMesh.cellPoints();
1152 
1153  // forAllConstIter(labelHashSet, limitCells, iter)
1154  // {
1155  // label celli = iter.key();
1156 
1157  // const labelList& cP = cellPts[celli];
1158 
1159  // forAll(cP, cPI)
1160  // {
1161  // ptToBeLimited[cP[cPI]] = true;
1162  // }
1163  // }
1164 
1165 
1166  // Apply Delaunay cell filterCounts and determine the maximum
1167  // overall filterCount
1168 
1169  label maxFilterCount = 0;
1170 
1171  for
1172  (
1173  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1174  cit != finite_cells_end();
1175  ++cit
1176  )
1177  {
1178  label cI = cit->cellIndex();
1179 
1180  if (cI >= 0)
1181  {
1182  if (ptToBeLimited[cI] == true)
1183  {
1184  cit->filterCount()++;
1185  }
1186 
1187  if (cit->filterCount() > maxFilterCount)
1188  {
1189  maxFilterCount = cit->filterCount();
1190  }
1191  }
1192  }
1193 
1194  Info<< nl << "Maximum number of filter limits applied: "
1195  << returnReduce(maxFilterCount, maxOp<label>()) << endl;
1196 
1197  return wrongFaces;
1198 }
1199 
1200 
1201 Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
1202 (
1203  Cell_handle cit
1204 ) const
1205 {
1206  if (cit->boundaryDualVertex())
1207  {
1208  if (cit->featurePointDualVertex())
1209  {
1210  return featurePoint;
1211  }
1212  else if (cit->featureEdgeDualVertex())
1213  {
1214  return featureEdge;
1215  }
1216  else
1217  {
1218  return surface;
1219  }
1220  }
1221  else if (cit->baffleSurfaceDualVertex())
1222  {
1223  return surface;
1224  }
1225  else if (cit->baffleEdgeDualVertex())
1226  {
1227  return featureEdge;
1228  }
1229  else
1230  {
1231  return internal;
1232  }
1233 }
1234 
1235 
1236 void Foam::conformalVoronoiMesh::indexDualVertices
1237 (
1238  pointField& pts,
1239  labelList& boundaryPts
1240 )
1241 {
1242  // Indexing Delaunay cells, which are the dual vertices
1243 
1244  this->resetCellCount();
1245 
1246  label nConstrainedVertices = 0;
1247  if (foamyHexMeshControls().guardFeaturePoints())
1248  {
1249  for
1250  (
1251  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1252  vit != finite_vertices_end();
1253  ++vit
1254  )
1255  {
1256  if (vit->constrained())
1257  {
1258  vit->index() = number_of_finite_cells() + nConstrainedVertices;
1259  nConstrainedVertices++;
1260  }
1261  }
1262  }
1263 
1264  pts.setSize(number_of_finite_cells() + nConstrainedVertices);
1265  boundaryPts.setSize
1266  (
1267  number_of_finite_cells() + nConstrainedVertices,
1268  internal
1269  );
1270 
1271  if (foamyHexMeshControls().guardFeaturePoints())
1272  {
1273  nConstrainedVertices = 0;
1274  for
1275  (
1276  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1277  vit != finite_vertices_end();
1278  ++vit
1279  )
1280  {
1281  if (vit->constrained())
1282  {
1283  pts[number_of_finite_cells() + nConstrainedVertices] =
1284  topoint(vit->point());
1285 
1286  boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
1287  constrained;
1288 
1289  nConstrainedVertices++;
1290  }
1291  }
1292  }
1293 
1294  OBJstream snapping1("snapToSurface1.obj");
1295  OBJstream snapping2("snapToSurface2.obj");
1296  OFstream tetToSnapTo("tetsToSnapTo.obj");
1297 
1298  for
1299  (
1300  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1301  cit != finite_cells_end();
1302  ++cit
1303  )
1304  {
1305 // if (tetrahedron(cit).volume() == 0)
1306 // {
1307 // Pout<< "ZERO VOLUME TET" << endl;
1308 // Pout<< cit->info();
1309 // Pout<< "Dual = " << cit->dual();
1310 // }
1311 
1312  if (!cit->hasFarPoint())
1313  {
1314  cit->cellIndex() = getNewCellIndex();
1315 
1316  // For nearly coplanar Delaunay cells that are present on different
1317  // processors the result of the circumcentre calculation depends on
1318  // the ordering of the vertices, so synchronise it across processors
1319 
1320  if (Pstream::parRun() && cit->parallelDualVertex())
1321  {
1322  typedef CGAL::Exact_predicates_exact_constructions_kernel Exact;
1323  typedef CGAL::Point_3<Exact> ExactPoint;
1324 
1325  List<labelPair> cellVerticesPair(4);
1326  List<ExactPoint> cellVertices(4);
1327 
1328  for (label vI = 0; vI < 4; ++vI)
1329  {
1330  cellVerticesPair[vI] = labelPair
1331  (
1332  cit->vertex(vI)->procIndex(),
1333  cit->vertex(vI)->index()
1334  );
1335 
1336  cellVertices[vI] = ExactPoint
1337  (
1338  cit->vertex(vI)->point().x(),
1339  cit->vertex(vI)->point().y(),
1340  cit->vertex(vI)->point().z()
1341  );
1342  }
1343 
1344  // Sort the vertices so that they will be in the same order on
1345  // each processor
1346  labelList oldToNew;
1347  sortedOrder(cellVerticesPair, oldToNew);
1348  oldToNew = invert(oldToNew.size(), oldToNew);
1349  inplaceReorder(oldToNew, cellVertices);
1350 
1351  ExactPoint synchronisedDual = CGAL::circumcenter
1352  (
1353  cellVertices[0],
1354  cellVertices[1],
1355  cellVertices[2],
1356  cellVertices[3]
1357  );
1358 
1359  pts[cit->cellIndex()] = Foam::point
1360  (
1361  CGAL::to_double(synchronisedDual.x()),
1362  CGAL::to_double(synchronisedDual.y()),
1363  CGAL::to_double(synchronisedDual.z())
1364  );
1365  }
1366  else
1367  {
1368  pts[cit->cellIndex()] = cit->dual();
1369  }
1370 
1371  // Feature point snapping
1372  if (foamyHexMeshControls().snapFeaturePoints())
1373  {
1374  if (cit->featurePointDualVertex())
1375  {
1376  pointFromPoint dual = cit->dual();
1377 
1378  pointIndexHit fpHit;
1379  label featureHit;
1380 
1381  // Find nearest feature point and compare
1382  geometryToConformTo_.findFeaturePointNearest
1383  (
1384  dual,
1385  sqr(targetCellSize(dual)),
1386  fpHit,
1387  featureHit
1388  );
1389 
1390  if (fpHit.hit())
1391  {
1392  if (debug)
1393  {
1394  Info<< "Dual = " << dual << nl
1395  << " Nearest = " << fpHit.hitPoint() << endl;
1396  }
1397 
1398  pts[cit->cellIndex()] = fpHit.hitPoint();
1399  }
1400  }
1401  }
1402 
1403 // {
1404 // // Snapping points far outside
1405 // if (cit->boundaryDualVertex() && !cit->parallelDualVertex())
1406 // {
1407 // pointFromPoint dual = cit->dual();
1408 //
1409 // pointIndexHit hitInfo;
1410 // label surfHit;
1411 //
1412 // // Find nearest surface point
1413 // geometryToConformTo_.findSurfaceNearest
1414 // (
1415 // dual,
1416 // sqr(targetCellSize(dual)),
1417 // hitInfo,
1418 // surfHit
1419 // );
1420 //
1421 // if (!hitInfo.hit())
1422 // {
1423 // // Project dual to nearest point on tet
1424 //
1425 // tetPointRef tet
1426 // (
1427 // topoint(cit->vertex(0)->point()),
1428 // topoint(cit->vertex(1)->point()),
1429 // topoint(cit->vertex(2)->point()),
1430 // topoint(cit->vertex(3)->point())
1431 // );
1432 //
1433 // pointFromPoint nearestPointOnTet =
1434 // tet.nearestPoint(dual).rawPoint();
1435 //
1436 // // Get nearest point on surface from tet.
1437 // geometryToConformTo_.findSurfaceNearest
1438 // (
1439 // nearestPointOnTet,
1440 // sqr(targetCellSize(nearestPointOnTet)),
1441 // hitInfo,
1442 // surfHit
1443 // );
1444 //
1445 // vector snapDir = nearestPointOnTet - dual;
1446 // snapDir /= mag(snapDir) + SMALL;
1447 //
1448 // drawDelaunayCell(tetToSnapTo, cit, offset);
1449 // offset += 1;
1450 //
1451 // vectorField norm(1);
1452 // allGeometry_[surfHit].getNormal
1453 // (
1454 // List<pointIndexHit>(1, hitInfo),
1455 // norm
1456 // );
1457 // norm[0] /= mag(norm[0]) + SMALL;
1458 //
1459 // if
1460 // (
1461 // hitInfo.hit()
1462 // && (mag(snapDir & norm[0]) > 0.5)
1463 // )
1464 // {
1465 // snapping1.write
1466 // (
1467 // linePointRef(dual, nearestPointOnTet)
1468 // );
1469 //
1470 // snapping2.write
1471 // (
1472 // linePointRef
1473 // (
1474 // nearestPointOnTet,
1475 // hitInfo.hitPoint()
1476 // )
1477 // );
1478 //
1479 // pts[cit->cellIndex()] = hitInfo.hitPoint();
1480 // }
1481 // }
1482 // }
1483 // }
1484 
1485  boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
1486  }
1487  else
1488  {
1489  cit->cellIndex() = Cb::ctFar;
1490  }
1491  }
1492 
1493  //pts.setSize(this->cellCount());
1494 
1495  //boundaryPts.setSize(this->cellCount());
1496 }
1497 
1498 
1499 void Foam::conformalVoronoiMesh::reindexDualVertices
1500 (
1501  const Map<label>& dualPtIndexMap,
1502  labelList& boundaryPts
1503 )
1504 {
1505  for
1506  (
1507  Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1508  cit != finite_cells_end();
1509  ++cit
1510  )
1511  {
1512  if (dualPtIndexMap.found(cit->cellIndex()))
1513  {
1514  cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
1515  boundaryPts[cit->cellIndex()] =
1516  max
1517  (
1518  boundaryPts[cit->cellIndex()],
1519  boundaryPts[dualPtIndexMap[cit->cellIndex()]]
1520  );
1521  }
1522  }
1523 }
1524 
1525 
1526 Foam::label Foam::conformalVoronoiMesh::createPatchInfo
1527 (
1528  wordList& patchNames,
1529  PtrList<dictionary>& patchDicts
1530 ) const
1531 {
1532  patchNames = geometryToConformTo_.patchNames();
1533 
1534  patchDicts.setSize(patchNames.size() + 1);
1535 
1536  const PtrList<dictionary>& patchInfo = geometryToConformTo_.patchInfo();
1537 
1538  forAll(patchNames, patchi)
1539  {
1540  if (patchInfo.set(patchi))
1541  {
1542  patchDicts.set(patchi, new dictionary(patchInfo[patchi]));
1543  }
1544  else
1545  {
1546  patchDicts.set(patchi, new dictionary());
1547  patchDicts[patchi].set
1548  (
1549  "type",
1550  wallPolyPatch::typeName
1551  );
1552  }
1553  }
1554 
1555  patchNames.setSize(patchNames.size() + 1);
1556  label defaultPatchIndex = patchNames.size() - 1;
1557  patchNames[defaultPatchIndex] = "foamyHexMesh_defaultPatch";
1558  patchDicts.set(defaultPatchIndex, new dictionary());
1559  patchDicts[defaultPatchIndex].set
1560  (
1561  "type",
1562  wallPolyPatch::typeName
1563  );
1564 
1565  label nProcPatches = 0;
1566 
1567  if (Pstream::parRun())
1568  {
1569  List<boolList> procUsedList
1570  (
1571  Pstream::nProcs(),
1572  boolList(Pstream::nProcs(), false)
1573  );
1574 
1575  boolList& procUsed = procUsedList[Pstream::myProcNo()];
1576 
1577  // Determine which processor patches are required
1578  for
1579  (
1580  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1581  vit != finite_vertices_end();
1582  vit++
1583  )
1584  {
1585  // This test is not sufficient if one of the processors does
1586  // not receive a referred vertex from another processor, but does
1587  // send one to the other processor.
1588  if (vit->referred())
1589  {
1590  procUsed[vit->procIndex()] = true;
1591  }
1592  }
1593 
1594  // Because the previous test was insufficient, combine the lists.
1595  Pstream::gatherList(procUsedList);
1596  Pstream::scatterList(procUsedList);
1597 
1598  forAll(procUsedList, proci)
1599  {
1600  if (proci != Pstream::myProcNo())
1601  {
1602  if (procUsedList[proci][Pstream::myProcNo()])
1603  {
1604  procUsed[proci] = true;
1605  }
1606  }
1607  }
1608 
1609  forAll(procUsed, pUI)
1610  {
1611  if (procUsed[pUI])
1612  {
1613  nProcPatches++;
1614  }
1615  }
1616 
1617  label nNonProcPatches = patchNames.size();
1618  label nTotalPatches = nNonProcPatches + nProcPatches;
1619 
1620  patchNames.setSize(nTotalPatches);
1621  patchDicts.setSize(nTotalPatches);
1622  for (label pI = nNonProcPatches; pI < nTotalPatches; ++pI)
1623  {
1624  patchDicts.set(pI, new dictionary());
1625  }
1626 
1627  label procAddI = 0;
1628 
1629  forAll(procUsed, pUI)
1630  {
1631  if (procUsed[pUI])
1632  {
1633  patchNames[nNonProcPatches + procAddI] =
1635 
1636  patchDicts[nNonProcPatches + procAddI].set
1637  (
1638  "type",
1639  processorPolyPatch::typeName
1640  );
1641 
1642  patchDicts[nNonProcPatches + procAddI].set
1643  (
1644  "myProcNo",
1646  );
1647 
1648  patchDicts[nNonProcPatches + procAddI].set("neighbProcNo", pUI);
1649 
1650  procAddI++;
1651  }
1652  }
1653  }
1654 
1655  return defaultPatchIndex;
1656 }
1657 
1658 
1659 Foam::vector Foam::conformalVoronoiMesh::calcSharedPatchNormal
1660 (
1661  Cell_handle c1,
1662  Cell_handle c2
1663 ) const
1664 {
1665  List<Foam::point> patchEdge(2, point::max);
1666 
1667  // Get shared Facet
1668  for (label cI = 0; cI < 4; ++cI)
1669  {
1670  if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
1671  {
1672  if (c1->vertex(cI)->internalBoundaryPoint())
1673  {
1674  patchEdge[0] = topoint(c1->vertex(cI)->point());
1675  }
1676  else
1677  {
1678  patchEdge[1] = topoint(c1->vertex(cI)->point());
1679  }
1680  }
1681  }
1682 
1683  Info<< " " << patchEdge << endl;
1684 
1685  return vector(patchEdge[1] - patchEdge[0]);
1686 }
1687 
1688 
1689 bool Foam::conformalVoronoiMesh::boundaryDualFace
1690 (
1691  Cell_handle c1,
1692  Cell_handle c2
1693 ) const
1694 {
1695  label nInternal = 0;
1696  label nExternal = 0;
1697 
1698  for (label cI = 0; cI < 4; ++cI)
1699  {
1700  if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
1701  {
1702  if (c1->vertex(cI)->internalBoundaryPoint())
1703  {
1704  nInternal++;
1705  }
1706  else if (c1->vertex(cI)->externalBoundaryPoint())
1707  {
1708  nExternal++;
1709  }
1710  }
1711  }
1712 
1713  Info<< "in = " << nInternal << " out = " << nExternal << endl;
1714 
1715  return (nInternal == 1 && nExternal == 1);
1716 }
1717 
1718 
1719 void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
1720 (
1721  const pointField& pts,
1722  faceList& faces,
1723  labelList& owner,
1724  labelList& neighbour,
1725  wordList& patchNames,
1726  PtrList<dictionary>& patchDicts,
1727  labelListList& patchPointPairSlaves,
1728  PackedBoolList& boundaryFacesToRemove,
1729  bool includeEmptyPatches
1730 ) const
1731 {
1732  const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);
1733 
1734  const label nPatches = patchNames.size();
1735 
1736  labelList procNeighbours(nPatches, label(-1));
1737  forAll(procNeighbours, patchi)
1738  {
1739  if (patchDicts[patchi].found("neighbProcNo"))
1740  {
1741  procNeighbours[patchi] =
1742  (
1743  patchDicts[patchi].found("neighbProcNo")
1744  ? readLabel(patchDicts[patchi].lookup("neighbProcNo"))
1745  : -1
1746  );
1747  }
1748  }
1749 
1750  List<DynamicList<face>> patchFaces(nPatches, DynamicList<face>(0));
1751  List<DynamicList<label>> patchOwners(nPatches, DynamicList<label>(0));
1752  // Per patch face the index of the slave node of the point pair
1753  List<DynamicList<label>> patchPPSlaves(nPatches, DynamicList<label>(0));
1754 
1755  List<DynamicList<bool>> indirectPatchFace(nPatches, DynamicList<bool>(0));
1756 
1757 
1758  faces.setSize(number_of_finite_edges());
1759  owner.setSize(number_of_finite_edges());
1760  neighbour.setSize(number_of_finite_edges());
1761  boundaryFacesToRemove.setSize(number_of_finite_edges(), false);
1762 
1763  labelPairPairDynListList procPatchSortingIndex(nPatches);
1764 
1765  label dualFacei = 0;
1766 
1767  if (foamyHexMeshControls().guardFeaturePoints())
1768  {
1769  OBJstream startCellStr("startingCell.obj");
1770  OBJstream featurePointFacesStr("ftPtFaces.obj");
1771  OBJstream featurePointDualsStr("ftPtDuals.obj");
1772  OFstream cellStr("vertexCells.obj");
1773 
1774  label vcount = 1;
1775 
1776  for
1777  (
1778  Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1779  vit != finite_vertices_end();
1780  ++vit
1781  )
1782  {
1783  if (vit->constrained())
1784  {
1785  // Find a starting cell
1786  std::list<Cell_handle> vertexCells;
1787  finite_incident_cells(vit, std::back_inserter(vertexCells));
1788 
1789  Cell_handle startCell;
1790 
1791  for
1792  (
1793  std::list<Cell_handle>::iterator vcit = vertexCells.begin();
1794  vcit != vertexCells.end();
1795  ++vcit
1796  )
1797  {
1798  if ((*vcit)->featurePointExternalCell())
1799  {
1800  startCell = *vcit;
1801  }
1802 
1803  if ((*vcit)->real())
1804  {
1805  featurePointDualsStr.write
1806  (
1807  linePointRef(topoint(vit->point()), (*vcit)->dual())
1808  );
1809  }
1810  }
1811 
1812  // Error if startCell is null
1813  if (startCell == nullptr)
1814  {
1815  Pout<< "Start cell is null!" << endl;
1816  }
1817 
1818  // Need to pick a direction to walk in
1819  Cell_handle vc1 = startCell;
1820  Cell_handle vc2;
1821 
1822  Info<< "c1 index = " << vc1->cellIndex() << " "
1823  << vc1->dual() << endl;
1824 
1825  for (label cI = 0; cI < 4; ++cI)
1826  {
1827  Info<< "c1 = " << cI << " "
1828  << vc1->neighbor(cI)->cellIndex() << " v = "
1829  << vc1->neighbor(cI)->dual() << endl;
1830 
1831  Info<< vc1->vertex(cI)->info();
1832  }
1833 
1834  Cell_handle nextCell;
1835 
1836  for (label cI = 0; cI < 4; ++cI)
1837  {
1838  if (vc1->vertex(cI)->externalBoundaryPoint())
1839  {
1840  vc2 = vc1->neighbor(cI);
1841 
1842  Info<< " c2 is neighbor "
1843  << vc2->cellIndex()
1844  << " of c1" << endl;
1845 
1846  for (label cI = 0; cI < 4; ++cI)
1847  {
1848  Info<< " c2 = " << cI << " "
1849  << vc2->neighbor(cI)->cellIndex() << " v = "
1850  << vc2->vertex(cI)->index() << endl;
1851  }
1852 
1853  face f(3);
1854  f[0] = vit->index();
1855  f[1] = vc1->cellIndex();
1856  f[2] = vc2->cellIndex();
1857 
1858  Info<< "f " << f << endl;
1859  forAll(f, pI)
1860  {
1861  Info<< " " << pts[f[pI]] << endl;
1862  }
1863 
1864  vector correctNormal = calcSharedPatchNormal(vc1, vc2);
1865  correctNormal /= mag(correctNormal);
1866 
1867  Info<< " cN " << correctNormal << endl;
1868 
1869  vector fN = f.normal(pts);
1870 
1871  if (mag(fN) < SMALL)
1872  {
1873  nextCell = vc2;
1874  continue;
1875  }
1876 
1877  fN /= mag(fN);
1878  Info<< " fN " << fN << endl;
1879 
1880  if ((fN & correctNormal) > 0)
1881  {
1882  nextCell = vc2;
1883  break;
1884  }
1885  }
1886  }
1887 
1888  vc2 = nextCell;
1889 
1890  label own = vit->index();
1891  face f(3);
1892  f[0] = own;
1893 
1894  Info<< "Start walk from " << vc1->cellIndex()
1895  << " to " << vc2->cellIndex() << endl;
1896 
1897  // Walk while not at start cell
1898 
1899  label iter = 0;
1900  do
1901  {
1902  Info<< " Walk from " << vc1->cellIndex()
1903  << " " << vc1->dual()
1904  << " to " << vc2->cellIndex()
1905  << " " << vc2->dual()
1906  << endl;
1907 
1908  startCellStr.write(linePointRef(vc1->dual(), vc2->dual()));
1909 
1910  // Get patch by getting face between cells and the two
1911  // points on the face that are not the feature vertex
1912  label patchIndex =
1913  geometryToConformTo_.findPatch
1914  (
1915  topoint(vit->point())
1916  );
1917 
1918  f[1] = vc1->cellIndex();
1919  f[2] = vc2->cellIndex();
1920 
1921  patchFaces[patchIndex].append(f);
1922  patchOwners[patchIndex].append(own);
1923  patchPPSlaves[patchIndex].append(own);
1924 
1925  // Find next cell
1926  Cell_handle nextCell;
1927 
1928  Info<< " c1 vertices " << vc2->dual() << endl;
1929  for (label cI = 0; cI < 4; ++cI)
1930  {
1931  Info<< " " << vc2->vertex(cI)->info();
1932  }
1933  Info<< " c1 neighbour vertices " << endl;
1934  for (label cI = 0; cI < 4; ++cI)
1935  {
1936  if
1937  (
1938  !vc2->vertex(cI)->constrained()
1939  && vc2->neighbor(cI) != vc1
1940  && !is_infinite(vc2->neighbor(cI))
1941  &&
1942  (
1943  vc2->neighbor(cI)->featurePointExternalCell()
1944  || vc2->neighbor(cI)->featurePointInternalCell()
1945  )
1946  && vc2->neighbor(cI)->hasConstrainedPoint()
1947  )
1948  {
1950  (
1951  cellStr,
1952  vc2->neighbor(cI),
1953  vcount++
1954  );
1955 
1956  Info<< " neighbour " << cI << " "
1957  << vc2->neighbor(cI)->dual() << endl;
1958  for (label I = 0; I < 4; ++I)
1959  {
1960  Info<< " "
1961  << vc2->neighbor(cI)->vertex(I)->info();
1962  }
1963  }
1964  }
1965 
1966  for (label cI = 0; cI < 4; ++cI)
1967  {
1968  if
1969  (
1970  !vc2->vertex(cI)->constrained()
1971  && vc2->neighbor(cI) != vc1
1972  && !is_infinite(vc2->neighbor(cI))
1973  &&
1974  (
1975  vc2->neighbor(cI)->featurePointExternalCell()
1976  || vc2->neighbor(cI)->featurePointInternalCell()
1977  )
1978  && vc2->neighbor(cI)->hasConstrainedPoint()
1979  )
1980  {
1981  // check if shared edge is internal/internal
1982  if (boundaryDualFace(vc2, vc2->neighbor(cI)))
1983  {
1984  nextCell = vc2->neighbor(cI);
1985  break;
1986  }
1987  }
1988  }
1989 
1990  vc1 = vc2;
1991  vc2 = nextCell;
1992 
1993  iter++;
1994  } while (vc1 != startCell && iter < 100);
1995  }
1996  }
1997  }
1998 
1999  for
2000  (
2001  Delaunay::Finite_edges_iterator eit = finite_edges_begin();
2002  eit != finite_edges_end();
2003  ++eit
2004  )
2005  {
2006  Cell_handle c = eit->first;
2007  Vertex_handle vA = c->vertex(eit->second);
2008  Vertex_handle vB = c->vertex(eit->third);
2009 
2010  if (vA->constrained() && vB->constrained())
2011  {
2012  continue;
2013  }
2014 
2015  if
2016  (
2017  (vA->constrained() && vB->internalOrBoundaryPoint())
2018  || (vB->constrained() && vA->internalOrBoundaryPoint())
2019  )
2020  {
2021  face newDualFace = buildDualFace(eit);
2022 
2023  label own = -1;
2024  label nei = -1;
2025 
2026  if (ownerAndNeighbour(vA, vB, own, nei))
2027  {
2028  reverse(newDualFace);
2029  }
2030 
2031  // internal face
2032  faces[dualFacei] = newDualFace;
2033  owner[dualFacei] = own;
2034  neighbour[dualFacei] = nei;
2035 
2036  dualFacei++;
2037  }
2038  else if
2039  (
2040  (vA->internalOrBoundaryPoint() && !vA->referred())
2041  || (vB->internalOrBoundaryPoint() && !vB->referred())
2042  )
2043  {
2044  if
2045  (
2046  (vA->internalPoint() && vB->externalBoundaryPoint())
2047  || (vB->internalPoint() && vA->externalBoundaryPoint())
2048  )
2049  {
2050  Cell_circulator ccStart = incident_cells(*eit);
2051  Cell_circulator cc1 = ccStart;
2052  Cell_circulator cc2 = cc1;
2053 
2054  cc2++;
2055 
2056  bool skipEdge = false;
2057 
2058  do
2059  {
2060  if
2061  (
2062  cc1->hasFarPoint() || cc2->hasFarPoint()
2063  || is_infinite(cc1) || is_infinite(cc2)
2064  )
2065  {
2066  Pout<< "Ignoring edge between internal and external: "
2067  << vA->info()
2068  << vB->info();
2069 
2070  skipEdge = true;
2071  break;
2072  }
2073 
2074  cc1++;
2075  cc2++;
2076 
2077  } while (cc1 != ccStart);
2078 
2079 
2080  // Do not create faces if the internal point is outside!
2081  // This occurs because the internal point is not determined to
2082  // be outside in the inside/outside test. This is most likely
2083  // due to the triangle.nearestPointClassify test not returning
2084  // edge/point as the nearest type.
2085 
2086  if (skipEdge)
2087  {
2088  continue;
2089  }
2090  }
2091 
2092  face newDualFace = buildDualFace(eit);
2093 
2094  if (newDualFace.size() >= 3)
2095  {
2096  label own = -1;
2097  label nei = -1;
2098 
2099  if (ownerAndNeighbour(vA, vB, own, nei))
2100  {
2101  reverse(newDualFace);
2102  }
2103 
2104  label patchIndex = -1;
2105 
2106  pointFromPoint ptA = topoint(vA->point());
2107  pointFromPoint ptB = topoint(vB->point());
2108 
2109  if (nei == -1)
2110  {
2111  // boundary face
2112 
2113  if (isProcBoundaryEdge(eit))
2114  {
2115  // One (and only one) of the points is an internal
2116  // point from another processor
2117 
2118  label procIndex = max(vA->procIndex(), vB->procIndex());
2119 
2120  patchIndex = max
2121  (
2122  findIndex(procNeighbours, vA->procIndex()),
2123  findIndex(procNeighbours, vB->procIndex())
2124  );
2125 
2126  // The lower processor index is the owner of the
2127  // two for the purpose of sorting the patch faces.
2128 
2129  if (Pstream::myProcNo() < procIndex)
2130  {
2131  // Use this processor's vertex index as the master
2132  // for sorting
2133 
2134  DynamicList<Pair<labelPair>>& sortingIndex =
2135  procPatchSortingIndex[patchIndex];
2136 
2137  if (vB->internalOrBoundaryPoint() && vB->referred())
2138  {
2139  sortingIndex.append
2140  (
2141  Pair<labelPair>
2142  (
2143  labelPair(vA->index(), vA->procIndex()),
2144  labelPair(vB->index(), vB->procIndex())
2145  )
2146  );
2147  }
2148  else
2149  {
2150  sortingIndex.append
2151  (
2152  Pair<labelPair>
2153  (
2154  labelPair(vB->index(), vB->procIndex()),
2155  labelPair(vA->index(), vA->procIndex())
2156  )
2157  );
2158  }
2159  }
2160  else
2161  {
2162  // Use the other processor's vertex index as the
2163  // master for sorting
2164 
2165  DynamicList<Pair<labelPair>>& sortingIndex =
2166  procPatchSortingIndex[patchIndex];
2167 
2168  if (vA->internalOrBoundaryPoint() && vA->referred())
2169  {
2170  sortingIndex.append
2171  (
2172  Pair<labelPair>
2173  (
2174  labelPair(vA->index(), vA->procIndex()),
2175  labelPair(vB->index(), vB->procIndex())
2176  )
2177  );
2178  }
2179  else
2180  {
2181  sortingIndex.append
2182  (
2183  Pair<labelPair>
2184  (
2185  labelPair(vB->index(), vB->procIndex()),
2186  labelPair(vA->index(), vA->procIndex())
2187  )
2188  );
2189  }
2190  }
2191 
2192 // Pout<< ptA << " " << ptB
2193 // << " proc indices "
2194 // << vA->procIndex() << " " << vB->procIndex()
2195 // << " indices " << vA->index()
2196 // << " " << vB->index()
2197 // << " my proc " << Pstream::myProcNo()
2198 // << " addedIndex "
2199 // << procPatchSortingIndex[patchIndex].last()
2200 // << endl;
2201  }
2202  else
2203  {
2204  patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
2205  }
2206 
2207  if (patchIndex == -1)
2208  {
2209  // Did not find a surface patch between
2210  // between Dv pair, finding nearest patch
2211 
2212 // Pout<< "Did not find a surface patch between "
2213 // << "for face, finding nearest patch to"
2214 // << 0.5*(ptA + ptB) << endl;
2215 
2216  patchIndex = geometryToConformTo_.findPatch
2217  (
2218  0.5*(ptA + ptB)
2219  );
2220  }
2221 
2222  patchFaces[patchIndex].append(newDualFace);
2223  patchOwners[patchIndex].append(own);
2224 
2225  // If the two vertices are a pair, then the patch face is
2226  // a desired one.
2227  if
2228  (
2229  vA->boundaryPoint() && vB->boundaryPoint()
2230  && !ptPairs_.isPointPair(vA, vB)
2231  && !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
2232  )
2233  {
2234  indirectPatchFace[patchIndex].append(true);
2235  }
2236  else
2237  {
2238  indirectPatchFace[patchIndex].append(false);
2239  }
2240 
2241  // Store the non-internal or boundary point
2242  if (vA->internalOrBoundaryPoint())
2243  {
2244  patchPPSlaves[patchIndex].append(vB->index());
2245  }
2246  else
2247  {
2248  patchPPSlaves[patchIndex].append(vA->index());
2249  }
2250  }
2251  else
2252  {
2253 // if
2254 // (
2255 // ptPairs_.isPointPair(vA, vB)
2256 // || ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
2257 // )
2258 // {
2259  patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
2260 // }
2261 
2262  if (patchIndex != -1)
2263  {
2264 // if
2265 // (
2266 // vA->boundaryPoint() && vB->boundaryPoint()
2267 // && !ptPairs_.isPointPair(vA, vB)
2268 // )
2269 // {
2270 // indirectPatchFace[patchIndex].append(true);
2271 // }
2272 // else
2273 // {
2274 // indirectPatchFace[patchIndex].append(false);
2275 // }
2276 // patchFaces[patchIndex].append(newDualFace);
2277 // patchOwners[patchIndex].append(own);
2278 // indirectPatchFace[patchIndex].append(false);
2279 //
2280 // if
2281 // (
2282 // labelPair(vB->index(), vB->procIndex())
2283 // < labelPair(vA->index(), vA->procIndex())
2284 // )
2285 // {
2286 // patchPPSlaves[patchIndex].append(vB->index());
2287 // }
2288 // else
2289 // {
2290 // patchPPSlaves[patchIndex].append(vA->index());
2291 // }
2292  }
2293 // else
2294  {
2295  // internal face
2296  faces[dualFacei] = newDualFace;
2297  owner[dualFacei] = own;
2298  neighbour[dualFacei] = nei;
2299 
2300  dualFacei++;
2301  }
2302  }
2303  }
2304  }
2305  }
2306 
2307  if (!patchFaces[defaultPatchIndex].empty())
2308  {
2309  Pout<< nl << patchFaces[defaultPatchIndex].size()
2310  << " faces were not able to have their patch determined from "
2311  << "the surface. "
2312  << nl << "Adding to patch " << patchNames[defaultPatchIndex]
2313  << endl;
2314  }
2315 
2316  label nInternalFaces = dualFacei;
2317 
2318  faces.setSize(nInternalFaces);
2319  owner.setSize(nInternalFaces);
2320  neighbour.setSize(nInternalFaces);
2321 
2322  timeCheck("polyMesh quality checked");
2323 
2324  sortFaces(faces, owner, neighbour);
2325 
2326  sortProcPatches
2327  (
2328  patchFaces,
2329  patchOwners,
2330  patchPPSlaves,
2331  procPatchSortingIndex
2332  );
2333 
2334  timeCheck("faces, owner, neighbour sorted");
2335 
2336  addPatches
2337  (
2338  nInternalFaces,
2339  faces,
2340  owner,
2341  patchDicts,
2342  boundaryFacesToRemove,
2343  patchFaces,
2344  patchOwners,
2345  indirectPatchFace
2346  );
2347 
2348  // Return patchPointPairSlaves.setSize(nPatches);
2349  patchPointPairSlaves.setSize(nPatches);
2350  forAll(patchPPSlaves, patchi)
2351  {
2352  patchPointPairSlaves[patchi].transfer(patchPPSlaves[patchi]);
2353  }
2354 
2355  if (foamyHexMeshControls().objOutput())
2356  {
2357  Info<< "Writing processor interfaces" << endl;
2358 
2359  forAll(patchDicts, nbI)
2360  {
2361  if (patchFaces[nbI].size() > 0)
2362  {
2363  const label neighbour =
2364  (
2365  patchDicts[nbI].found("neighbProcNo")
2366  ? readLabel(patchDicts[nbI].lookup("neighbProcNo"))
2367  : -1
2368  );
2369 
2370  faceList procPatchFaces = patchFaces[nbI];
2371 
2372  // Reverse faces as it makes it easier to analyse the output
2373  // using a diff
2374  if (neighbour < Pstream::myProcNo())
2375  {
2376  forAll(procPatchFaces, fI)
2377  {
2378  procPatchFaces[fI] = procPatchFaces[fI].reverseFace();
2379  }
2380  }
2381 
2382  if (neighbour != -1)
2383  {
2384  word fName =
2385  "processor_"
2386  + name(Pstream::myProcNo())
2387  + "_to_"
2388  + name(neighbour)
2389  + "_interface.obj";
2390 
2392  (
2393  time().path()/fName,
2394  *this,
2395  procPatchFaces
2396  );
2397  }
2398  }
2399  }
2400  }
2401 }
2402 
2403 
2404 void Foam::conformalVoronoiMesh::sortFaces
2405 (
2406  faceList& faces,
2407  labelList& owner,
2408  labelList& neighbour
2409 ) const
2410 {
2411  // Upper triangular order:
2412  // + owner is sorted in ascending cell order
2413  // + within each block of equal value for owner, neighbour is sorted in
2414  // ascending cell order.
2415  // + faces sorted to correspond
2416  // e.g.
2417  // owner | neighbour
2418  // 0 | 2
2419  // 0 | 23
2420  // 0 | 71
2421  // 1 | 23
2422  // 1 | 24
2423  // 1 | 91
2424 
2425  List<labelPair> ownerNeighbourPair(owner.size());
2426 
2427  forAll(ownerNeighbourPair, oNI)
2428  {
2429  ownerNeighbourPair[oNI] = labelPair(owner[oNI], neighbour[oNI]);
2430  }
2431 
2432  Info<< nl
2433  << "Sorting faces, owner and neighbour into upper triangular order"
2434  << endl;
2435 
2436  labelList oldToNew;
2437 
2438  sortedOrder(ownerNeighbourPair, oldToNew);
2439 
2440  oldToNew = invert(oldToNew.size(), oldToNew);
2441 
2442  inplaceReorder(oldToNew, faces);
2443  inplaceReorder(oldToNew, owner);
2444  inplaceReorder(oldToNew, neighbour);
2445 }
2446 
2447 
2448 void Foam::conformalVoronoiMesh::sortProcPatches
2449 (
2450  List<DynamicList<face>>& patchFaces,
2451  List<DynamicList<label>>& patchOwners,
2452  List<DynamicList<label>>& patchPointPairSlaves,
2453  labelPairPairDynListList& patchSortingIndices
2454 ) const
2455 {
2456  if (!Pstream::parRun())
2457  {
2458  return;
2459  }
2460 
2461  forAll(patchSortingIndices, patchi)
2462  {
2463  faceList& faces = patchFaces[patchi];
2464  labelList& owner = patchOwners[patchi];
2465  DynamicList<label>& slaves = patchPointPairSlaves[patchi];
2466  DynamicList<Pair<labelPair>>& sortingIndices
2467  = patchSortingIndices[patchi];
2468 
2469  if (!sortingIndices.empty())
2470  {
2471  if
2472  (
2473  faces.size() != sortingIndices.size()
2474  || owner.size() != sortingIndices.size()
2475  || slaves.size() != sortingIndices.size()
2476  )
2477  {
2479  << "patch size and size of sorting indices is inconsistent "
2480  << " for patch " << patchi << nl
2481  << " faces.size() " << faces.size() << nl
2482  << " owner.size() " << owner.size() << nl
2483  << " slaves.size() " << slaves.size() << nl
2484  << " sortingIndices.size() "
2485  << sortingIndices.size()
2486  << exit(FatalError) << endl;
2487  }
2488 
2489  labelList oldToNew;
2490 
2491  sortedOrder(sortingIndices, oldToNew);
2492 
2493  oldToNew = invert(oldToNew.size(), oldToNew);
2494 
2495  inplaceReorder(oldToNew, sortingIndices);
2496  inplaceReorder(oldToNew, faces);
2497  inplaceReorder(oldToNew, owner);
2498  inplaceReorder(oldToNew, slaves);
2499  }
2500  }
2501 }
2502 
2503 
2504 void Foam::conformalVoronoiMesh::addPatches
2505 (
2506  const label nInternalFaces,
2507  faceList& faces,
2508  labelList& owner,
2509  PtrList<dictionary>& patchDicts,
2510  PackedBoolList& boundaryFacesToRemove,
2511  const List<DynamicList<face>>& patchFaces,
2512  const List<DynamicList<label>>& patchOwners,
2513  const List<DynamicList<bool>>& indirectPatchFace
2514 ) const
2515 {
2516  label nBoundaryFaces = 0;
2517 
2518  forAll(patchFaces, p)
2519  {
2520  patchDicts[p].set("nFaces", patchFaces[p].size());
2521  patchDicts[p].set("startFace", nInternalFaces + nBoundaryFaces);
2522 
2523  nBoundaryFaces += patchFaces[p].size();
2524  }
2525 
2526  faces.setSize(nInternalFaces + nBoundaryFaces);
2527  owner.setSize(nInternalFaces + nBoundaryFaces);
2528  boundaryFacesToRemove.setSize(nInternalFaces + nBoundaryFaces);
2529 
2530  label facei = nInternalFaces;
2531 
2532  forAll(patchFaces, p)
2533  {
2534  forAll(patchFaces[p], f)
2535  {
2536  faces[facei] = patchFaces[p][f];
2537  owner[facei] = patchOwners[p][f];
2538  boundaryFacesToRemove[facei] = indirectPatchFace[p][f];
2539 
2540  facei++;
2541  }
2542  }
2543 }
2544 
2545 
2546 void Foam::conformalVoronoiMesh::removeUnusedPoints
2547 (
2548  faceList& faces,
2549  pointField& pts,
2550  labelList& boundaryPts
2551 ) const
2552 {
2553  Info<< nl << "Removing unused points" << endl;
2554 
2555  PackedBoolList ptUsed(pts.size(), false);
2556 
2557  // Scan all faces to find all of the points that are used
2558 
2559  forAll(faces, fI)
2560  {
2561  const face& f = faces[fI];
2562 
2563  forAll(f, fPtI)
2564  {
2565  ptUsed[f[fPtI]] = true;
2566  }
2567  }
2568 
2569  label pointi = 0;
2570 
2571  labelList oldToNew(pts.size(), label(-1));
2572 
2573  // Move all of the used points to the start of the pointField and
2574  // truncate it
2575 
2576  forAll(ptUsed, ptUI)
2577  {
2578  if (ptUsed[ptUI] == true)
2579  {
2580  oldToNew[ptUI] = pointi++;
2581  }
2582  }
2583 
2584  inplaceReorder(oldToNew, pts);
2585  inplaceReorder(oldToNew, boundaryPts);
2586 
2587  Info<< " Removing "
2588  << returnReduce(pts.size() - pointi, sumOp<label>())
2589  << " unused points"
2590  << endl;
2591 
2592  pts.setSize(pointi);
2593  boundaryPts.setSize(pointi);
2594 
2595  // Renumber the faces to use the new point numbers
2596 
2597  forAll(faces, fI)
2598  {
2599  inplaceRenumber(oldToNew, faces[fI]);
2600  }
2601 }
2602 
2603 
2604 Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
2605 (
2606  labelList& owner,
2607  labelList& neighbour
2608 ) const
2609 {
2610  Info<< nl << "Removing unused cells" << endl;
2611 
2612  PackedBoolList cellUsed(vertexCount(), false);
2613 
2614  // Scan all faces to find all of the cells that are used
2615 
2616  forAll(owner, oI)
2617  {
2618  cellUsed[owner[oI]] = true;
2619  }
2620 
2621  forAll(neighbour, nI)
2622  {
2623  cellUsed[neighbour[nI]] = true;
2624  }
2625 
2626  label celli = 0;
2627 
2628  labelList oldToNew(cellUsed.size(), label(-1));
2629 
2630  // Move all of the used cellCentres to the start of the pointField and
2631  // truncate it
2632 
2633  forAll(cellUsed, cellUI)
2634  {
2635  if (cellUsed[cellUI] == true)
2636  {
2637  oldToNew[cellUI] = celli++;
2638  }
2639  }
2640 
2641  labelList newToOld(invert(celli, oldToNew));
2642 
2643  // Find all of the unused cells, create a list of them, then
2644  // subtract one from each owner and neighbour entry for each of
2645  // the unused cell indices that it is above.
2646 
2647  DynamicList<label> unusedCells;
2648 
2649  forAll(cellUsed, cUI)
2650  {
2651  if (cellUsed[cUI] == false)
2652  {
2653  unusedCells.append(cUI);
2654  }
2655  }
2656 
2657  if (unusedCells.size() > 0)
2658  {
2659  Info<< " Removing "
2660  << returnReduce(unusedCells.size(), sumOp<label>())
2661  << " unused cell labels" << endl;
2662 
2663  forAll(owner, oI)
2664  {
2665  label& o = owner[oI];
2666 
2667  o -= findLower(unusedCells, o) + 1;
2668  }
2669 
2670  forAll(neighbour, nI)
2671  {
2672  label& n = neighbour[nI];
2673 
2674  n -= findLower(unusedCells, n) + 1;
2675  }
2676  }
2677 
2678  return newToOld;
2679 }
2680 
2681 
2682 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
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.
Delaunay::Cell_handle Cell_handle
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Delaunay::Vertex_handle Vertex_handle
pointFromPoint topoint(const Point &P)
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static fvMesh * meshPtr
Definition: globalFoam.H:52
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:418
List< face > faceList
Definition: faceListFwd.H:43
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
static word newName(const label myProcNo, const label neighbProcNo)
Return the name of a processorPolyPatch.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
void findFeaturePointNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &fpHit, label &featureHit) const
Find the nearest point on any feature edge.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static void timeCheck(const Time &runTime, const string &description=string::null, const bool check=true)
Write the elapsedCpuTime and memory usage, with an optional.
Xfer< T > xferCopy(const T &)
Construct by copying the contents of the arg.
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:644
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
patches[0]
const dictionary & foamyHexMeshDict() const
Return the foamyHexMeshDict.
Definition: cvControlsI.H:28
void drawDelaunayCell(Ostream &os, const CellHandle &c, label offset=0)
Draws a tet cell to an output stream. The offset is supplied as the tet.
UList< label > labelUList
Definition: UList.H:64
Various functions to operate on Lists.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:692
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:116
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
static const Identity< scalar > I
Definition: Identity.H:93
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
wordList patchNames(nPatches)
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
List< label > labelList
A List of labels.
Definition: labelList.H:56
void findSurfaceNearest(const point &sample, scalar nearestDistSqr, pointIndexHit &surfHit, label &hitSurface) const
Find the nearest point to the sample and return it to the.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label readLabel(Istream &is)
Definition: label.H:64
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
static const char nl
Definition: Ostream.H:262
void writeProcessorInterface(const fileName &fName, const Triangulation &t, const faceList &faces)
Write the processor interface to an OBJ file.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
bool isPointPair(const Vertex_handle &vA, const Vertex_handle &vB) const
const List< word > & patchNames() const
Return the patch names.
List< word > wordList
A List of words.
Definition: fileName.H:54
void setSize(const label)
Reset size of List.
Definition: List.C:281
void resetCellCount()
Set the cell count to zero.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:394
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:400
label patchi
label vertexCount() const
Return the vertex count (the next unique vertex index)
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
label getNewCellIndex() const
Create a new unique cell index and return.
Definition: DelaunayMeshI.H:63
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
Field< vector > vectorField
Specialisation of Field<T> for vector.
static autoPtr< polyPatch > New(const word &patchType, const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm)
Return a pointer to a new patch created on freestore from.
Definition: polyPatchNew.C:32
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
virtual Ostream & write(const token &)=0
Write next token to stream.
const cvControls & foamyHexMeshControls() const
Return the foamyHexMeshControls object.
const Time & time() const
Return the Time object.
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:537
volScalarField & p
List< DynamicList< Pair< labelPair > > > labelPairPairDynListList
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
bool found
const pointPairs< Delaunay > & featurePointPairs() const
Return the feature point pair table.
List< cell > cellList
list of cells
Definition: cellList.H:42
InfoProxy< IOstream > info() const
Return info proxy.
Definition: IOstream.H:531
const PtrList< dictionary > & patchInfo() const
Return the patch info.
label findPatch(const point &ptA, const point &ptB) const
Find which patch is intersected by the line from one point to.