conformalVoronoiMeshIO.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-2017 OpenFOAM Foundation
9 Copyright (C) 2015-2018 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 "IOstreams.H"
31#include "OFstream.H"
32#include "pointMesh.H"
33#include "pointFields.H"
34#include "ListOps.H"
35#include "polyMeshFilter.H"
36#include "polyTopoChange.H"
37#include "PrintTable.H"
38#include "indexedVertexOps.H"
39#include "DelaunayMeshTools.H"
40#include "syncTools.H"
41#include "faceSet.H"
42#include "OBJstream.H"
43
44// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
45
47(
48 const string& description
49) const
50{
51 timeCheck(time(), description, foamyHexMeshControls().timeChecks());
52}
53
54
56(
57 const Time& runTime,
58 const string& description,
59 const bool check
60)
61{
62 if (check)
63 {
64 Info<< nl << "--- [ cpuTime "
65 << runTime.elapsedCpuTime() << " s, "
66 << "delta " << runTime.cpuTimeIncrement()<< " s";
67
68 if (description != word::null)
69 {
70 Info<< ", " << description << " ";
71 }
72 else
73 {
74 Info<< " ";
75 }
76
77 Info<< "] --- " << endl;
78
79 memInfo m;
80
81 if (m.valid())
82 {
83 PrintTable<word, label> memoryTable
84 (
85 "Memory Usage (kB): "
86 + description
87 );
88
89 memoryTable.add("mSize", m.size());
90 memoryTable.add("mPeak", m.peak());
91 memoryTable.add("mRss", m.rss());
92
94 memoryTable.print(Info, true, true);
96 }
97 }
98}
99
100
101void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
102{
104
105 // Per cell the Delaunay vertex
106 labelList cellToDelaunayVertex;
107 // Per patch, per face the Delaunay vertex
108 labelListList patchToDelaunayVertex;
109
110 labelList dualPatchStarts;
111
112 {
114 labelList boundaryPts;
115 faceList faces;
116 labelList owner;
117 labelList neighbour;
119 PtrList<dictionary> patchDicts;
120 pointField cellCentres;
121 bitSet boundaryFacesToRemove;
122
123 calcDualMesh
124 (
125 points,
126 boundaryPts,
127 faces,
128 owner,
129 neighbour,
132 cellCentres,
133 cellToDelaunayVertex,
134 patchToDelaunayVertex,
135 boundaryFacesToRemove
136 );
137
138 Info<< nl << "Writing polyMesh to " << instance << endl;
139
140 writeMesh
141 (
143 instance,
144 points,
145 boundaryPts,
146 faces,
147 owner,
148 neighbour,
151 cellCentres,
152 boundaryFacesToRemove
153 );
154
155 dualPatchStarts.setSize(patchDicts.size());
156
157 forAll(dualPatchStarts, patchi)
158 {
159 dualPatchStarts[patchi] =
160 patchDicts[patchi].get<label>("startFace");
161 }
162 }
163
164 if (foamyHexMeshControls().writeCellShapeControlMesh())
165 {
166 cellShapeControls().shapeControlMesh().write();
167 }
168
169 if (foamyHexMeshControls().writeBackgroundMeshDecomposition())
170 {
171 Info<< nl << "Writing " << "backgroundMeshDecomposition" << endl;
172
173 // Have to explicitly update the mesh instance.
174 const_cast<fvMesh&>(decomposition_().mesh()).setInstance
175 (
176 time().timeName()
177 );
178
179 decomposition_().mesh().write();
180 }
181
182 if (foamyHexMeshControls().writeTetDualMesh())
183 {
184 label celli = 0;
185 for
186 (
187 Finite_cells_iterator cit = finite_cells_begin();
188 cit != finite_cells_end();
189 ++cit
190 )
191 {
192 if
193 (
194 !cit->hasFarPoint()
195 && !is_infinite(cit)
196 )
197 {
198 cit->cellIndex() = celli++;
199 }
200 }
201
202 Info<< nl << "Writing " << "tetDualMesh" << endl;
203
204 labelPairLookup vertexMap;
205 labelList cellMap;
206 autoPtr<polyMesh> tetMesh =
207 createMesh("tetDualMesh", vertexMap, cellMap);
208
209 tetMesh().write();
210
211// // Determine map from Delaunay vertex to Dual mesh
212// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213//
214// // From all Delaunay vertices to cell (positive index)
215// // or patch face (negative index)
216// labelList vertexToDualAddressing(number_of_vertices(), Zero);
217//
218// forAll(cellToDelaunayVertex, celli)
219// {
220// label vertI = cellToDelaunayVertex[celli];
221//
222// if (vertexToDualAddressing[vertI] != 0)
223// {
224// FatalErrorInFunction
225// << "Delaunay vertex " << vertI
226// << " from cell " << celli
227// << " is already mapped to "
228// << vertexToDualAddressing[vertI]
229// << exit(FatalError);
230// }
231// vertexToDualAddressing[vertI] = celli+1;
232// }
233//
234// forAll(patchToDelaunayVertex, patchi)
235// {
236// const labelList& patchVertices = patchToDelaunayVertex[patchi];
237//
238// forAll(patchVertices, i)
239// {
240// label vertI = patchVertices[i];
241//
242// if (vertexToDualAddressing[vertI] > 0)
243// {
244// FatalErrorInFunction
245// << "Delaunay vertex " << vertI
246// << " from patch " << patchi
247// << " local index " << i
248// << " is already mapped to cell "
249// << vertexToDualAddressing[vertI]-1
250// << exit(FatalError);
251// }
252//
253// // Vertex might be used by multiple faces. Which one to
254// // use? For now last one wins.
255// label dualFacei = dualPatchStarts[patchi]+i;
256// vertexToDualAddressing[vertI] = -dualFacei-1;
257// }
258// }
259//
260//
261// // Calculate tet mesh addressing
262// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
263//
264// pointField points;
265// labelList boundaryPts(number_of_finite_cells(), -1);
266// // From tet point back to Delaunay vertex index
267// labelList pointToDelaunayVertex;
268// faceList faces;
269// labelList owner;
270// labelList neighbour;
271// wordList patchTypes;
272// wordList patchNames;
273// PtrList<dictionary> patchDicts;
274// pointField cellCentres;
275//
276// calcTetMesh
277// (
278// points,
279// pointToDelaunayVertex,
280// faces,
281// owner,
282// neighbour,
283// patchTypes,
284// patchNames,
285// patchDicts
286// );
287//
288//
289//
290// // Calculate map from tet points to dual mesh cells/patch faces
291// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
292//
293// labelIOList pointDualAddressing
294// (
295// IOobject
296// (
297// "pointDualAddressing",
298// instance,
299// "tetDualMesh"/polyMesh::meshSubDir,
300// runTime_,
301// IOobject::NO_READ,
302// IOobject::AUTO_WRITE,
303// false
304// ),
305// labelUIndList
306// (
307// vertexToDualAddressing,
308// pointToDelaunayVertex
309// )()
310// );
311//
312// label pointi = pointDualAddressing.find(-1);
313// if (pointi != -1)
314// {
315// WarningInFunction
316// << "Delaunay vertex " << pointi
317// << " does not have a corresponding dual cell." << endl;
318// }
319//
320// Info<< "Writing map from tetDualMesh points to Voronoi mesh to "
321// << pointDualAddressing.objectPath() << endl;
322// pointDualAddressing.write();
323//
324//
325//
326// // Write tet points corresponding to the Voronoi cell/face centre
327// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
328// {
329// // Read Voronoi mesh
330// fvMesh mesh
331// (
332// IOobject
333// (
334// polyMesh::defaultRegion,
335// instance,
336// runTime_,
337// IOobject::MUST_READ
338// )
339// );
340// pointIOField dualPoints
341// (
342// IOobject
343// (
344// "dualPoints",
345// instance,
346// "tetDualMesh"/polyMesh::meshSubDir,
347// runTime_,
348// IOobject::NO_READ,
349// IOobject::AUTO_WRITE,
350// false
351// ),
352// points
353// );
354//
355// forAll(pointDualAddressing, pointi)
356// {
357// label index = pointDualAddressing[pointi];
358//
359// if (index > 0)
360// {
361// label celli = index-1;
362// dualPoints[pointi] = mesh.cellCentres()[celli];
363// }
364// else if (index < 0)
365// {
366// label facei = -index-1;
367// if (facei >= mesh.nInternalFaces())
368// {
369// dualPoints[pointi] = mesh.faceCentres()[facei];
370// }
371// }
372// }
373//
374// Info<< "Writing tetDualMesh points mapped onto Voronoi mesh to "
375// << dualPoints.objectPath() << endl
376// << "Replace the polyMesh/points with these." << endl;
377// dualPoints.write();
378// }
379//
380//
381// Info<< nl << "Writing tetDualMesh to " << instance << endl;
382//
383// bitSet boundaryFacesToRemove;
384// writeMesh
385// (
386// "tetDualMesh",
387// instance,
388// points,
389// boundaryPts,
390// faces,
391// owner,
392// neighbour,
393// patchTypes,
394// patchNames,
395// patchDicts,
396// cellCentres,
397// boundaryFacesToRemove
398// );
399 }
400}
401
402
403Foam::autoPtr<Foam::fvMesh> Foam::conformalVoronoiMesh::createDummyMesh
404(
405 const IOobject& io,
406 const wordList& patchNames,
407 const PtrList<dictionary>& patchDicts
408) const
409{
411 fvMesh& mesh = meshPtr();
412
413 List<polyPatch*> patches(patchDicts.size());
414
415 forAll(patches, patchi)
416 {
417 if
418 (
419 patchDicts.set(patchi)
420 && (
421 patchDicts[patchi].get<word>("type")
423 )
424 )
425 {
426 patches[patchi] = new processorPolyPatch
427 (
428 0, //patchSizes[p],
429 0, //patchStarts[p],
430 patchi,
432 patchDicts[patchi].get<label>("myProcNo"),
433 patchDicts[patchi].get<label>("neighbProcNo"),
435 );
436 }
437 else
438 {
439 patches[patchi] = polyPatch::New
440 (
441 patchDicts[patchi].get<word>("type"),
442 patchNames[patchi],
443 0, //patchSizes[p],
444 0, //patchStarts[p],
445 patchi,
447 ).ptr();
448 }
449 }
450
452
453 return meshPtr;
454}
455
456
457void Foam::conformalVoronoiMesh::checkProcessorPatchesMatch
458(
459 const PtrList<dictionary>& patchDicts
460) const
461{
462 // Check patch sizes
463 labelListList procPatchSizes
464 (
467 );
468
469 forAll(patchDicts, patchi)
470 {
471 if
472 (
473 patchDicts.set(patchi)
474 && (
475 patchDicts[patchi].get<word>("type")
477 )
478 )
479 {
480 const label procNeighb =
481 patchDicts[patchi].get<label>("neighbProcNo");
482
483 procPatchSizes[Pstream::myProcNo()][procNeighb]
484 = patchDicts[patchi].get<label>("nFaces");
485 }
486 }
487
488 Pstream::gatherList(procPatchSizes);
489
490 if (Pstream::master())
491 {
492 bool allMatch = true;
493
494 forAll(procPatchSizes, proci)
495 {
496 const labelList& patchSizes = procPatchSizes[proci];
497
498 forAll(patchSizes, patchi)
499 {
500 if (patchSizes[patchi] != procPatchSizes[patchi][proci])
501 {
502 allMatch = false;
503
504 Info<< indent << "Patches " << proci << " and " << patchi
505 << " have different sizes: " << patchSizes[patchi]
506 << " and " << procPatchSizes[patchi][proci] << endl;
507 }
508 }
509 }
510
511 if (allMatch)
512 {
513 Info<< indent << "All processor patches have matching numbers of "
514 << "faces" << endl;
515 }
516 }
517}
518
519
520void Foam::conformalVoronoiMesh::reorderPoints
521(
523 labelList& boundaryPts,
524 faceList& faces,
525 const label nInternalFaces
526) const
527{
528 Info<< incrIndent << indent << "Reordering points into internal/external"
529 << endl;
530
531 labelList oldToNew(points.size(), Zero);
532
533 // Find points that are internal
534 for (label fI = nInternalFaces; fI < faces.size(); ++fI)
535 {
536 const face& f = faces[fI];
537
538 forAll(f, fpI)
539 {
540 oldToNew[f[fpI]] = 1;
541 }
542 }
543
544 const label nInternalPoints = points.size() - sum(oldToNew);
545
546 label countInternal = 0;
547 label countExternal = nInternalPoints;
548
549 forAll(points, pI)
550 {
551 if (oldToNew[pI] == 0)
552 {
553 oldToNew[pI] = countInternal++;
554 }
555 else
556 {
557 oldToNew[pI] = countExternal++;
558 }
559 }
560
561 Info<< indent
562 << "Number of internal points: " << countInternal << nl
563 << indent << "Number of external points: " << countExternal
564 << decrIndent << endl;
565
566 inplaceReorder(oldToNew, points);
567 inplaceReorder(oldToNew, boundaryPts);
568
569 forAll(faces, fI)
570 {
571 face& f = faces[fI];
572
573 forAll(f, fpI)
574 {
575 f[fpI] = oldToNew[f[fpI]];
576 }
577 }
578}
579
580
581void Foam::conformalVoronoiMesh::reorderProcessorPatches
582(
583 const word& meshName,
584 const fileName& instance,
585 const pointField& points,
586 faceList& faces,
587 const wordList& patchNames,
588 const PtrList<dictionary>& patchDicts
589) const
590{
591 Info<< incrIndent << indent << "Reordering processor patches" << endl;
592
594 checkProcessorPatchesMatch(patchDicts);
595
596 // Create dummy mesh with correct proc boundaries to do sorting
597 autoPtr<fvMesh> sortMeshPtr
598 (
599 createDummyMesh
600 (
601 IOobject
602 (
603 meshName,
604 instance,
605 runTime_,
608 false
609 ),
612 )
613 );
614 const fvMesh& sortMesh = sortMeshPtr();
615
616 // Rotation on new faces.
617 labelList rotation(faces.size(), Zero);
618 labelList faceMap(faces.size(), label(-1));
619
620 PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
621
622 // Send ordering
623 forAll(sortMesh.boundaryMesh(), patchi)
624 {
625 const polyPatch& pp = sortMesh.boundaryMesh()[patchi];
626
627 if (isA<processorPolyPatch>(pp))
628 {
629 refCast<const processorPolyPatch>(pp).initOrder
630 (
631 pBufs,
633 (
634 SubList<face>
635 (
636 faces,
637 patchDicts[patchi].get<label>("nFaces"),
638 patchDicts[patchi].get<label>("startFace")
639 ),
640 points
641 )
642 );
643 }
644 }
645
646 pBufs.finishedSends();
647
648 Info<< incrIndent << indent << "Face ordering initialised..." << endl;
649
650 // Receive and calculate ordering
651 bool anyChanged = false;
652
653 forAll(sortMesh.boundaryMesh(), patchi)
654 {
655 const polyPatch& pp = sortMesh.boundaryMesh()[patchi];
656
657 if (isA<processorPolyPatch>(pp))
658 {
659 const label nPatchFaces =
660 patchDicts[patchi].get<label>("nFaces");
661
662 const label patchStartFace =
663 patchDicts[patchi].get<label>("startFace");
664
665 labelList patchFaceMap(nPatchFaces, label(-1));
666 labelList patchFaceRotation(nPatchFaces, Zero);
667
668 bool changed = refCast<const processorPolyPatch>(pp).order
669 (
670 pBufs,
672 (
673 SubList<face>
674 (
675 faces,
676 nPatchFaces,
677 patchStartFace
678 ),
679 points
680 ),
681 patchFaceMap,
682 patchFaceRotation
683 );
684
685 if (changed)
686 {
687 // Merge patch face reordering into mesh face reordering table
688 forAll(patchFaceRotation, patchFacei)
689 {
690 rotation[patchFacei + patchStartFace]
691 = patchFaceRotation[patchFacei];
692 }
693
694 forAll(patchFaceMap, patchFacei)
695 {
696 if (patchFaceMap[patchFacei] != patchFacei)
697 {
698 faceMap[patchFacei + patchStartFace]
699 = patchFaceMap[patchFacei] + patchStartFace;
700 }
701 }
702
703 anyChanged = true;
704 }
705 }
706 }
707
708 Info<< incrIndent << indent << "Faces matched." << endl;
709
710 reduce(anyChanged, orOp<bool>());
711
712 if (anyChanged)
713 {
714 label nReorderedFaces = 0;
715
716 forAll(faceMap, facei)
717 {
718 if (faceMap[facei] != -1)
719 {
720 nReorderedFaces++;
721 }
722 }
723
724 if (nReorderedFaces > 0)
725 {
726 inplaceReorder(faceMap, faces);
727 }
728
729 // Rotate faces (rotation is already in new face indices).
730 label nRotated = 0;
731
732 forAll(rotation, facei)
733 {
734 if (rotation[facei] != 0)
735 {
736 faces[facei] = rotateList(faces[facei], rotation[facei]);
737 nRotated++;
738 }
739 }
740
741 Info<< indent << returnReduce(nReorderedFaces, sumOp<label>())
742 << " faces have been reordered" << nl
743 << indent << returnReduce(nRotated, sumOp<label>())
744 << " faces have been rotated"
746 << decrIndent << decrIndent << endl;
747 }
748}
749
750
752(
753 const word& meshName,
754 const fileName& instance,
756 labelList& boundaryPts,
757 faceList& faces,
758 labelList& owner,
759 labelList& neighbour,
760 const wordList& patchNames,
761 const PtrList<dictionary>& patchDicts,
762 const pointField& cellCentres,
763 bitSet& boundaryFacesToRemove
764) const
765{
766 if (foamyHexMeshControls().objOutput())
767 {
769 (
770 time().path()/word(meshName + ".obj"),
771 points,
772 faces
773 );
774 }
775
776 const label nInternalFaces = patchDicts[0].get<label>("startFace");
777
778 reorderPoints(points, boundaryPts, faces, nInternalFaces);
779
780 if (Pstream::parRun())
781 {
782 reorderProcessorPatches
783 (
784 meshName,
785 instance,
786 points,
787 faces,
790 );
791 }
792
794 Info<< indent << "Constructing mesh" << endl;
795
796 timeCheck("Before fvMesh construction");
797
798 fvMesh mesh
799 (
800 IOobject
801 (
802 meshName,
803 instance,
804 runTime_,
807 ),
808 std::move(points),
809 std::move(faces),
810 std::move(owner),
811 std::move(neighbour)
812 );
813
814 Info<< indent << "Adding patches to mesh" << endl;
815
816 List<polyPatch*> patches(patchNames.size());
817
818 label nValidPatches = 0;
819
821 {
822 label totalPatchSize = patchDicts[p].get<label>("nFaces");
823
824 if
825 (
826 patchDicts.set(p)
827 && (
828 patchDicts[p].get<word>("type")
830 )
831 )
832 {
833 const_cast<dictionary&>(patchDicts[p]).set
834 (
835 "transform",
836 "coincidentFullMatch"
837 );
838
839 // Do not create empty processor patches
840 if (totalPatchSize > 0)
841 {
842 patches[nValidPatches] = new processorPolyPatch
843 (
844 patchNames[p],
845 patchDicts[p],
846 nValidPatches,
849 );
850
851 nValidPatches++;
852 }
853 }
854 else
855 {
856 // Check that the patch is not empty on every processor
857 reduce(totalPatchSize, sumOp<label>());
858
859 if (totalPatchSize > 0)
860 {
861 patches[nValidPatches] = polyPatch::New
862 (
863 patchNames[p],
864 patchDicts[p],
865 nValidPatches,
867 ).ptr();
868
869 nValidPatches++;
870 }
871 }
872 }
873
874 patches.setSize(nValidPatches);
875
877
878
879 // Add zones to the mesh
880 addZones(mesh, cellCentres);
881
882
883 Info<< indent << "Add pointZones" << endl;
884 {
885 label sz = mesh.pointZones().size();
886
887 DynamicList<label> bPts(boundaryPts.size());
888
889 forAll(dualMeshPointTypeNames_, typeI)
890 {
891 const word& znName =
892 dualMeshPointTypeNames_[dualMeshPointType(typeI)];
893
894 forAll(boundaryPts, ptI)
895 {
896 const label& bPtType = boundaryPts[ptI];
897
898 if (bPtType == typeI)
899 {
900 bPts.append(ptI);
901 }
902 }
903
904// syncTools::syncPointList(mesh, bPts, maxEqOp<label>(), -1);
905
907 << "Adding " << bPts.size()
908 << " points of type " << znName
909 << decrIndent << endl;
910
912 (
913 new pointZone
914 (
915 znName,
916 bPts,
917 sz + typeI,
919 )
920 );
921
922 bPts.clear();
923 }
924 }
925
926
927
928 // Add indirectPatchFaces to a face zone
929 Info<< indent << "Adding indirect patch faces set" << endl;
930
932 (
933 mesh,
934 boundaryFacesToRemove,
935 orEqOp<unsigned int>()
936 );
937
938 labelList addr(boundaryFacesToRemove.toc());
939
940 faceSet indirectPatchFaces
941 (
942 mesh,
943 "indirectPatchFaces",
944 addr,
946 );
947
948 indirectPatchFaces.sync(mesh);
949
950
952
953 timeCheck("Before fvMesh filtering");
954
955 autoPtr<polyMeshFilter> meshFilter;
956
957 label nInitialBadFaces = 0;
958
959 if (foamyHexMeshControls().filterEdges())
960 {
961 Info<< nl << "Filtering edges on polyMesh" << nl << endl;
962
963 meshFilter.reset(new polyMeshFilter(mesh, boundaryPts));
964
965 // Filter small edges only. This reduces the number of faces so that
966 // the face filtering is sped up.
967 nInitialBadFaces = meshFilter().filterEdges(0);
968 {
969 const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
970
971 polyTopoChange meshMod(newMesh());
972
973 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
974
975 polyMeshFilter::copySets(newMesh(), mesh);
976 }
977 }
978
979 if (foamyHexMeshControls().filterFaces())
980 {
981 labelIOList boundaryPtsIO
982 (
983 IOobject
984 (
985 "pointPriority",
986 instance,
987 time(),
990 ),
992 );
993
994 forAll(mesh.points(), ptI)
995 {
996 boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
997 }
998
999
1000 Info<< nl << "Filtering faces on polyMesh" << nl << endl;
1001
1002 meshFilter.reset(new polyMeshFilter(mesh, boundaryPtsIO));
1003
1004 meshFilter().filter(nInitialBadFaces);
1005 {
1006 const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
1007
1008 polyTopoChange meshMod(newMesh());
1009
1010 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
1011
1012 polyMeshFilter::copySets(newMesh(), mesh);
1013 }
1014 }
1015
1016 timeCheck("After fvMesh filtering");
1017
1018 mesh.setInstance(instance);
1019
1020 if (!mesh.write())
1021 {
1023 << "Failed writing polyMesh."
1024 << exit(FatalError);
1025 }
1026 else
1027 {
1028 Info<< nl << "Written filtered mesh to "
1029 << mesh.polyMesh::instance() << nl
1030 << endl;
1031 }
1032
1033 {
1034 pointScalarField boundaryPtsScalarField
1035 (
1036 IOobject
1037 (
1038 "boundaryPoints_collapsed",
1039 instance,
1040 time(),
1043 ),
1045 dimensionedScalar("min", dimless, scalar(labelMin))
1046 );
1047
1048 labelIOList boundaryPtsIO
1049 (
1050 IOobject
1051 (
1052 "pointPriority",
1053 instance,
1054 time(),
1057 ),
1059 );
1060
1061 forAll(mesh.points(), ptI)
1062 {
1063 boundaryPtsScalarField[ptI] = mesh.pointZones().whichZone(ptI);
1064 boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
1065 }
1066
1067 boundaryPtsScalarField.write();
1068 boundaryPtsIO.write();
1069 }
1070
1071// writeCellSizes(mesh);
1072
1073// writeCellAlignments(mesh);
1074
1075// writeCellCentres(mesh);
1076
1077 findRemainingProtrusionSet(mesh);
1078}
1079
1080
1082(
1083 const fvMesh& mesh
1084) const
1085{
1086 {
1087 timeCheck("Start writeCellSizes");
1088
1089 Info<< nl << "Create targetCellSize volScalarField" << endl;
1090
1091 volScalarField targetCellSize
1092 (
1093 IOobject
1094 (
1095 "targetCellSize",
1096 mesh.polyMesh::instance(),
1097 mesh,
1100 ),
1101 mesh,
1103 zeroGradientFvPatchScalarField::typeName
1104 );
1105
1106 scalarField& cellSize = targetCellSize.primitiveFieldRef();
1107
1108 const vectorField& C = mesh.cellCentres();
1109
1110 forAll(cellSize, i)
1111 {
1112 cellSize[i] = cellShapeControls().cellSize(C[i]);
1113 }
1114
1115 // Info<< nl << "Create targetCellVolume volScalarField" << endl;
1116
1117 // volScalarField targetCellVolume
1118 // (
1119 // IOobject
1120 // (
1121 // "targetCellVolume",
1122 // mesh.polyMesh::instance(),
1123 // mesh,
1124 // IOobject::NO_READ,
1125 // IOobject::AUTO_WRITE
1126 // ),
1127 // mesh,
1128 // dimensionedScalar(dimLength, Zero),
1129 // zeroGradientFvPatchScalarField::typeName
1130 // );
1131
1132 // targetCellVolume.primitiveFieldRef() = pow3(cellSize);
1133
1134 // Info<< nl << "Create actualCellVolume volScalarField" << endl;
1135
1136 // volScalarField actualCellVolume
1137 // (
1138 // IOobject
1139 // (
1140 // "actualCellVolume",
1141 // mesh.polyMesh::instance(),
1142 // mesh,
1143 // IOobject::NO_READ,
1144 // IOobject::AUTO_WRITE
1145 // ),
1146 // mesh,
1147 // dimensionedScalar(dimVolume, Zero),
1148 // zeroGradientFvPatchScalarField::typeName
1149 // );
1150
1151 // actualCellVolume.primitiveFieldRef() = mesh.cellVolumes();
1152
1153 // Info<< nl << "Create equivalentCellSize volScalarField" << endl;
1154
1155 // volScalarField equivalentCellSize
1156 // (
1157 // IOobject
1158 // (
1159 // "equivalentCellSize",
1160 // mesh.polyMesh::instance(),
1161 // mesh,
1162 // IOobject::NO_READ,
1163 // IOobject::AUTO_WRITE
1164 // ),
1165 // mesh,
1166 // dimensionedScalar(dimLength, Zero),
1167 // zeroGradientFvPatchScalarField::typeName
1168 // );
1169
1170 // equivalentCellSize.primitiveFieldRef() = pow
1171 // (
1172 // actualCellVolume.primitiveField(),
1173 // 1.0/3.0
1174 // );
1175
1176 targetCellSize.correctBoundaryConditions();
1177 // targetCellVolume.correctBoundaryConditions();
1178 // actualCellVolume.correctBoundaryConditions();
1179 // equivalentCellSize.correctBoundaryConditions();
1180
1181 targetCellSize.write();
1182 // targetCellVolume.write();
1183 // actualCellVolume.write();
1184 // equivalentCellSize.write();
1185 }
1186
1187 // {
1188 // polyMesh tetMesh
1189 // (
1190 // IOobject
1191 // (
1192 // "tetDualMesh",
1193 // runTime_.constant(),
1194 // runTime_,
1195 // IOobject::MUST_READ
1196 // )
1197 // );
1198
1199 // pointMesh ptMesh(tetMesh);
1200
1201 // pointScalarField ptTargetCellSize
1202 // (
1203 // IOobject
1204 // (
1205 // "ptTargetCellSize",
1206 // runTime_.timeName(),
1207 // tetMesh,
1208 // IOobject::NO_READ,
1209 // IOobject::AUTO_WRITE
1210 // ),
1211 // ptMesh,
1212 // dimensionedScalar(dimLength, Zero),
1213 // pointPatchVectorField::calculatedType()
1214 // );
1215
1216 // scalarField& cellSize = ptTargetCellSize.primitiveFieldRef();
1217
1218 // const vectorField& P = tetMesh.points();
1219
1220 // forAll(cellSize, i)
1221 // {
1222 // cellSize[i] = cellShapeControls().cellSize(P[i]);
1223 // }
1224
1225 // ptTargetCellSize.write();
1226 // }
1227}
1228
1229
1231(
1232 const fvMesh& mesh
1233) const
1234{
1235// Info<< nl << "Create cellAlignments volTensorField" << endl;
1236//
1237// volTensorField cellAlignments
1238// (
1239// IOobject
1240// (
1241// "cellAlignments",
1242// mesh.polyMesh::instance(),
1243// mesh,
1244// IOobject::NO_READ,
1245// IOobject::AUTO_WRITE
1246// ),
1247// mesh,
1248// tensor::I,
1249// zeroGradientFvPatchTensorField::typeName
1250// );
1251//
1252// tensorField& cellAlignment = cellAlignments.primitiveFieldRef();
1253//
1254// const vectorField& C = mesh.cellCentres();
1255//
1256// vectorField xDir(cellAlignment.size());
1257// vectorField yDir(cellAlignment.size());
1258// vectorField zDir(cellAlignment.size());
1259//
1260// forAll(cellAlignment, i)
1261// {
1262// cellAlignment[i] = cellShapeControls().cellAlignment(C[i]);
1263// xDir[i] = cellAlignment[i] & vector(1, 0, 0);
1264// yDir[i] = cellAlignment[i] & vector(0, 1, 0);
1265// zDir[i] = cellAlignment[i] & vector(0, 0, 1);
1266// }
1267//
1268// OFstream xStr("xDir.obj");
1269// OFstream yStr("yDir.obj");
1270// OFstream zStr("zDir.obj");
1271//
1272// forAll(xDir, i)
1273// {
1274// meshTools::writeOBJ(xStr, C[i], C[i] + xDir[i]);
1275// meshTools::writeOBJ(yStr, C[i], C[i] + yDir[i]);
1276// meshTools::writeOBJ(zStr, C[i], C[i] + zDir[i]);
1277// }
1278//
1279// cellAlignments.correctBoundaryConditions();
1280//
1281// cellAlignments.write();
1282}
1283
1284
1286(
1287 const fvMesh& mesh
1288) const
1289{
1290 Info<< "Writing components of cellCentre positions to volScalarFields"
1291 << " ccx, ccy, ccz in " << runTime_.timeName() << endl;
1292
1293 for (direction i=0; i<vector::nComponents; i++)
1294 {
1295 volScalarField cci
1296 (
1297 IOobject
1298 (
1299 "cc" + word(vector::componentNames[i]),
1300 runTime_.timeName(),
1301 mesh,
1304 ),
1305 mesh.C().component(i)
1306 );
1307
1308 cci.write();
1309 }
1310}
1311
1312
1314(
1315 const polyMesh& mesh
1316) const
1317{
1318 timeCheck("Start findRemainingProtrusionSet");
1319
1320 const polyBoundaryMesh& patches = mesh.boundaryMesh();
1321
1322 labelHashSet protrudingBoundaryPoints;
1323
1324 forAll(patches, patchi)
1325 {
1326 const polyPatch& patch = patches[patchi];
1327
1328 forAll(patch.localPoints(), pLPI)
1329 {
1330 label meshPtI = patch.meshPoints()[pLPI];
1331
1332 const Foam::point& pt = patch.localPoints()[pLPI];
1333
1334 if
1335 (
1336 geometryToConformTo_.wellOutside
1337 (
1338 pt,
1339 sqr(targetCellSize(pt))
1340 )
1341 )
1342 {
1343 protrudingBoundaryPoints.insert(meshPtI);
1344 }
1345 }
1346 }
1347
1348 cellSet protrudingCells
1349 (
1350 mesh,
1351 "foamyHexMesh_remainingProtrusions",
1352 mesh.nCells()/1000
1353 );
1354
1355 for (const label pointi : protrudingBoundaryPoints)
1356 {
1357 const labelList& pCells = mesh.pointCells()[pointi];
1358 protrudingCells.insert(pCells);
1359 }
1360
1361 const label protrudingCellsSize =
1362 returnReduce(protrudingCells.size(), sumOp<label>());
1363
1364 if (foamyHexMeshControls().objOutput() && protrudingCellsSize)
1365 {
1366 Info<< nl << "Found " << protrudingCellsSize
1367 << " cells protruding from the surface, writing cellSet "
1368 << protrudingCells.name()
1369 << endl;
1370
1371 protrudingCells.write();
1372 }
1373
1374 return std::move(protrudingCells);
1375}
1376
1377
1378void Foam::conformalVoronoiMesh::writePointPairs
1379(
1380 const fileName& fName
1381) const
1382{
1383 OBJstream os(fName);
1384
1385 for
1386 (
1387 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1388 eit != finite_edges_end();
1389 ++eit
1390 )
1391 {
1392 Cell_handle c = eit->first;
1393 Vertex_handle vA = c->vertex(eit->second);
1394 Vertex_handle vB = c->vertex(eit->third);
1395
1396 if (ptPairs_.isPointPair(vA, vB))
1397 {
1398 os.write
1399 (
1400 linePointRef(topoint(vA->point()), topoint(vB->point()))
1401 );
1402 }
1403 }
1404}
1405
1406
1407// ************************************************************************* //
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
Various functions to operate on Lists.
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
virtual void print(Ostream &os) const
Print stream description to Ostream.
Definition: IOstream.C:80
virtual Ostream & write(const char c)
Write character.
Definition: OBJstream.C:78
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 gatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
void append(T *ptr)
Append an element to the end of the list.
Definition: PtrListI.H:113
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
@ nonBlocking
"nonBlocking"
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
labelHashSet findRemainingProtrusionSet(const polyMesh &mesh) const
Find the cellSet of the boundary cells which have points that.
const Time & time() const
Return the Time object.
void writeCellSizes(const fvMesh &mesh) const
Calculate and write a field of the target cell size,.
const cvControls & foamyHexMeshControls() const
Return the foamyHexMeshControls object.
void writeCellAlignments(const fvMesh &mesh) const
static void timeCheck(const Time &runTime, const string &description=string::null, const bool check=true)
Write the elapsedCpuTime and memory usage, with an optional.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTimeCxx.C:75
double elapsedCpuTime() const
Return CPU time (in seconds) from the start.
Definition: cpuTimeCxx.C:68
void writeMesh() const
Write equivalent mesh information at the polyMesh faceInstances time.
Writes the cell-centres volVectorField and the three component fields as volScalarFields.
const volVectorField & C() const
Return cell centres as volVectorField.
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1079
void addFvPatches(polyPatchList &plist, const bool validBoundary=true)
Add boundary patches. Constructor helper.
Definition: fvMesh.C:632
static const char *const componentNames[]
Definition: bool.H:104
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
void setInstance(const fileName &instance, const IOobject::writeOption wOpt=IOobject::AUTO_WRITE)
Set the instance for mesh files.
Definition: polyMeshIO.C:36
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:321
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
const pointZoneMesh & pointZones() const noexcept
Return point zone mesh.
Definition: polyMesh.H:492
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
const labelListList & pointCells() const
const vectorField & cellCentres() const
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
int myProcNo() const noexcept
Return processor number.
splitCell * master() const
Definition: splitCell.H:113
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
static const word null
An empty word.
Definition: word.H:80
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
volScalarField & p
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
engineTime & runTime
Foam::autoPtr< Foam::dynamicFvMesh > meshPtr
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const pointField & points
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
void writeInternalDelaunayVertices(const fileName &instance, const Triangulation &t)
Write the internal Delaunay vertices of the tessellation as a.
void writeObjMesh(const fileName &fName, const pointField &points, const faceList &faces)
Write an OBJ mesh consisting of points and faces.
const dimensionedScalar c
Speed of light in a vacuum.
const std::string patch
OpenFOAM patch number as a std::string.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
List< word > wordList
A List of words.
Definition: fileName.H:63
const dimensionSet dimless
Dimensionless.
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
constexpr label labelMin
Definition: label.H:60
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
pointFromPoint topoint(const Point &P)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
ListType rotateList(const ListType &list, const label n)
Rotate a list by n places.
static void check(const int retVal, const char *what)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
HashTable< label, labelPair, Foam::Hash< labelPair > > labelPairLookup
This is a Map of a labelPair to a label. Used for e.g. for face1, face2 to shared edge....
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
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)
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:44
Field< vector > vectorField
Specialisation of Field<T> for vector.
uint8_t direction
Definition: direction.H:56
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
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
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
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:356
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)
volScalarField & C
#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.