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