snappyHexMesh.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) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-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
27Application
28 snappyHexMesh
29
30Group
31 grpMeshGenerationUtilities
32
33Description
34 Automatic split hex mesher. Refines and snaps to surface.
35
36\*---------------------------------------------------------------------------*/
37
38#include "argList.H"
39#include "Time.H"
40#include "fvMesh.H"
41#include "snappyRefineDriver.H"
42#include "snappySnapDriver.H"
43#include "snappyLayerDriver.H"
44#include "searchableSurfaces.H"
45#include "refinementSurfaces.H"
46#include "refinementFeatures.H"
47#include "shellSurfaces.H"
48#include "decompositionMethod.H"
49#include "fvMeshDistribute.H"
50#include "wallPolyPatch.H"
52#include "snapParameters.H"
53#include "layerParameters.H"
54#include "vtkCoordSetWriter.H"
55#include "faceSet.H"
56#include "motionSmoother.H"
57#include "polyTopoChange.H"
61#include "MeshedSurface.H"
62#include "globalIndex.H"
63#include "IOmanip.H"
64#include "decompositionModel.H"
65#include "fvMeshTools.H"
66#include "profiling.H"
67#include "processorMeshes.H"
69
70using namespace Foam;
71
72// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73
74// Convert size (as fraction of defaultCellSize) to refinement level
75label sizeCoeffToRefinement
76(
77 const scalar level0Coeff, // ratio of hex cell size v.s. defaultCellSize
78 const scalar sizeCoeff
79)
80{
81 return round(::log(level0Coeff/sizeCoeff)/::log(2));
82}
83
84
85autoPtr<refinementSurfaces> createRefinementSurfaces
86(
87 const searchableSurfaces& allGeometry,
88 const dictionary& surfacesDict,
89 const dictionary& shapeControlDict,
90 const label gapLevelIncrement,
91 const scalar level0Coeff
92)
93{
95
96 // Count number of surfaces.
97 label surfi = 0;
98 forAll(allGeometry.names(), geomi)
99 {
100 const word& geomName = allGeometry.names()[geomi];
101
102 if (surfacesDict.found(geomName))
103 {
104 surfi++;
105 }
106 }
107
108 labelList surfaces(surfi);
109 wordList names(surfi);
110 PtrList<surfaceZonesInfo> surfZones(surfi);
111
112 labelList regionOffset(surfi);
113
114 labelList globalMinLevel(surfi, Zero);
115 labelList globalMaxLevel(surfi, Zero);
116 labelList globalLevelIncr(surfi, Zero);
117 PtrList<dictionary> globalPatchInfo(surfi);
118 List<Map<label>> regionMinLevel(surfi);
119 List<Map<label>> regionMaxLevel(surfi);
120 List<Map<label>> regionLevelIncr(surfi);
121 List<Map<scalar>> regionAngle(surfi);
122 List<Map<autoPtr<dictionary>>> regionPatchInfo(surfi);
123
124 wordHashSet unmatchedKeys(surfacesDict.toc());
125
126 surfi = 0;
127 forAll(allGeometry.names(), geomi)
128 {
129 const word& geomName = allGeometry.names()[geomi];
130
131 const entry* ePtr = surfacesDict.findEntry(geomName, keyType::REGEX);
132
133 if (ePtr)
134 {
135 const dictionary& shapeDict = ePtr->dict();
136 unmatchedKeys.erase(ePtr->keyword());
137
138 names[surfi] = geomName;
139 surfaces[surfi] = geomi;
140
141 const searchableSurface& surface = allGeometry[geomi];
142
143 // Find the index in shapeControlDict
144 // Invert surfaceCellSize to get the refinementLevel
145
146 const word scsFuncName =
147 shapeDict.get<word>("surfaceCellSizeFunction");
148
149 const dictionary& scsDict =
150 shapeDict.optionalSubDict(scsFuncName + "Coeffs");
151
152 const scalar surfaceCellSize =
153 scsDict.get<scalar>("surfaceCellSizeCoeff");
154
155 const label refLevel = sizeCoeffToRefinement
156 (
157 level0Coeff,
158 surfaceCellSize
159 );
160
161 globalMinLevel[surfi] = refLevel;
162 globalMaxLevel[surfi] = refLevel;
163 globalLevelIncr[surfi] = gapLevelIncrement;
164
165 // Surface zones
166 surfZones.set
167 (
168 surfi,
170 (
171 surface,
172 shapeDict,
173 allGeometry.regionNames()[surfaces[surfi]]
174 )
175 );
176
177
178 // Global perpendicular angle
179 if (shapeDict.found("patchInfo"))
180 {
181 globalPatchInfo.set
182 (
183 surfi,
184 shapeDict.subDict("patchInfo").clone()
185 );
186 }
187
188
189 // Per region override of patchInfo
190
191 if (shapeDict.found("regions"))
192 {
193 const dictionary& regionsDict = shapeDict.subDict("regions");
194 const wordList& regionNames =
195 allGeometry[surfaces[surfi]].regions();
196
197 forAll(regionNames, regioni)
198 {
199 if (regionsDict.found(regionNames[regioni]))
200 {
201 // Get the dictionary for region
202 const dictionary& regionDict = regionsDict.subDict
203 (
204 regionNames[regioni]
205 );
206
207 if (regionDict.found("patchInfo"))
208 {
209 regionPatchInfo[surfi].insert
210 (
211 regioni,
212 regionDict.subDict("patchInfo").clone()
213 );
214 }
215 }
216 }
217 }
218
219 // Per region override of cellSize
220 if (shapeDict.found("regions"))
221 {
222 const dictionary& shapeControlRegionsDict =
223 shapeDict.subDict("regions");
224 const wordList& regionNames =
225 allGeometry[surfaces[surfi]].regions();
226
227 forAll(regionNames, regioni)
228 {
229 if (shapeControlRegionsDict.found(regionNames[regioni]))
230 {
231 const dictionary& shapeControlRegionDict =
232 shapeControlRegionsDict.subDict
233 (
234 regionNames[regioni]
235 );
236
237 const word scsFuncName =
238 shapeControlRegionDict.get<word>
239 (
240 "surfaceCellSizeFunction"
241 );
242 const dictionary& scsDict =
243 shapeControlRegionDict.subDict
244 (
245 scsFuncName + "Coeffs"
246 );
247
248 const scalar surfaceCellSize =
249 scsDict.get<scalar>("surfaceCellSizeCoeff");
250
251 const label refLevel = sizeCoeffToRefinement
252 (
253 level0Coeff,
254 surfaceCellSize
255 );
256
257 regionMinLevel[surfi].insert(regioni, refLevel);
258 regionMaxLevel[surfi].insert(regioni, refLevel);
259 regionLevelIncr[surfi].insert(regioni, 0);
260 }
261 }
262 }
263
264 surfi++;
265 }
266 }
267
268 // Calculate local to global region offset
269 label nRegions = 0;
270
271 forAll(surfaces, surfi)
272 {
273 regionOffset[surfi] = nRegions;
274 nRegions += allGeometry[surfaces[surfi]].regions().size();
275 }
276
277 // Rework surface specific information into information per global region
278 labelList minLevel(nRegions, Zero);
279 labelList maxLevel(nRegions, Zero);
280 labelList gapLevel(nRegions, -1);
281 PtrList<dictionary> patchInfo(nRegions);
282
283 forAll(globalMinLevel, surfi)
284 {
285 label nRegions = allGeometry[surfaces[surfi]].regions().size();
286
287 // Initialise to global (i.e. per surface)
288 for (label i = 0; i < nRegions; i++)
289 {
290 label globalRegioni = regionOffset[surfi] + i;
291 minLevel[globalRegioni] = globalMinLevel[surfi];
292 maxLevel[globalRegioni] = globalMaxLevel[surfi];
293 gapLevel[globalRegioni] =
294 maxLevel[globalRegioni]
295 + globalLevelIncr[surfi];
296
297 if (globalPatchInfo.set(surfi))
298 {
299 patchInfo.set
300 (
301 globalRegioni,
302 globalPatchInfo[surfi].clone()
303 );
304 }
305 }
306
307 // Overwrite with region specific information
308 forAllConstIters(regionMinLevel[surfi], iter)
309 {
310 label globalRegioni = regionOffset[surfi] + iter.key();
311
312 minLevel[globalRegioni] = iter();
313 maxLevel[globalRegioni] = regionMaxLevel[surfi][iter.key()];
314 gapLevel[globalRegioni] =
315 maxLevel[globalRegioni]
316 + regionLevelIncr[surfi][iter.key()];
317 }
318
319 const Map<autoPtr<dictionary>>& localInfo = regionPatchInfo[surfi];
320 forAllConstIters(localInfo, iter)
321 {
322 label globalRegioni = regionOffset[surfi] + iter.key();
323 patchInfo.set(globalRegioni, iter()().clone());
324 }
325 }
326
327 surfacePtr.reset
328 (
330 (
331 allGeometry,
332 surfaces,
333 names,
334 surfZones,
335 regionOffset,
336 minLevel,
337 maxLevel,
338 gapLevel,
339 scalarField(nRegions, -GREAT), //perpendicularAngle,
340 patchInfo,
341 false //dryRun
342 )
343 );
344
345
346 const refinementSurfaces& rf = surfacePtr();
347
348 // Determine maximum region name length
349 label maxLen = 0;
350 forAll(rf.surfaces(), surfi)
351 {
352 label geomi = rf.surfaces()[surfi];
353 const wordList& regionNames = allGeometry.regionNames()[geomi];
354 forAll(regionNames, regioni)
355 {
356 maxLen = Foam::max(maxLen, label(regionNames[regioni].size()));
357 }
358 }
359
360
361 Info<< setw(maxLen) << "Region"
362 << setw(10) << "Min Level"
363 << setw(10) << "Max Level"
364 << setw(10) << "Gap Level" << nl
365 << setw(maxLen) << "------"
366 << setw(10) << "---------"
367 << setw(10) << "---------"
368 << setw(10) << "---------" << endl;
369
370 forAll(rf.surfaces(), surfi)
371 {
372 label geomi = rf.surfaces()[surfi];
373
374 Info<< rf.names()[surfi] << ':' << nl;
375
376 const wordList& regionNames = allGeometry.regionNames()[geomi];
377
378 forAll(regionNames, regioni)
379 {
380 label globali = rf.globalRegion(surfi, regioni);
381
382 Info<< setw(maxLen) << regionNames[regioni]
383 << setw(10) << rf.minLevel()[globali]
384 << setw(10) << rf.maxLevel()[globali]
385 << setw(10) << rf.gapLevel()[globali] << endl;
386 }
387 }
388
389
390 return surfacePtr;
391}
392
393
394void extractSurface
395(
396 const polyMesh& mesh,
397 const Time& runTime,
398 const labelHashSet& includePatches,
399 const fileName& outFileName
400)
401{
402 const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
403
404 // Collect sizes. Hash on names to handle local-only patches (e.g.
405 // processor patches)
406 HashTable<label> patchSize(1024);
407 label nFaces = 0;
408 for (const label patchi : includePatches)
409 {
410 const polyPatch& pp = bMesh[patchi];
411 patchSize.insert(pp.name(), pp.size());
412 nFaces += pp.size();
413 }
414 Pstream::mapCombineGather(patchSize, plusEqOp<label>());
415
416
417 // Allocate zone/patch for all patches
418 HashTable<label> compactZoneID(1024);
419 forAllConstIters(patchSize, iter)
420 {
421 label sz = compactZoneID.size();
422 compactZoneID.insert(iter.key(), sz);
423 }
424 Pstream::broadcast(compactZoneID);
425
426
427 // Rework HashTable into labelList just for speed of conversion
428 labelList patchToCompactZone(bMesh.size(), -1);
429 forAllConstIters(compactZoneID, iter)
430 {
431 label patchi = bMesh.findPatchID(iter.key());
432 if (patchi != -1)
433 {
434 patchToCompactZone[patchi] = iter();
435 }
436 }
437
438 // Collect faces on zones
439 DynamicList<label> faceLabels(nFaces);
440 DynamicList<label> compactZones(nFaces);
441 for (const label patchi : includePatches)
442 {
443 const polyPatch& pp = bMesh[patchi];
444 forAll(pp, i)
445 {
446 faceLabels.append(pp.start()+i);
447 compactZones.append(patchToCompactZone[pp.index()]);
448 }
449 }
450
451 // Addressing engine for all faces
452 uindirectPrimitivePatch allBoundary
453 (
454 UIndirectList<face>(mesh.faces(), faceLabels),
455 mesh.points()
456 );
457
458
459 // Find correspondence to master points
460 labelList pointToGlobal;
461 labelList uniqueMeshPoints;
462 autoPtr<globalIndex> globalNumbers = mesh.globalData().mergePoints
463 (
464 allBoundary.meshPoints(),
465 allBoundary.meshPointMap(),
466 pointToGlobal,
467 uniqueMeshPoints
468 );
469
470 // Gather all unique points on master
471 List<pointField> gatheredPoints(Pstream::nProcs());
472 gatheredPoints[Pstream::myProcNo()] = pointField
473 (
474 mesh.points(),
475 uniqueMeshPoints
476 );
477 Pstream::gatherList(gatheredPoints);
478
479 // Gather all faces
480 List<faceList> gatheredFaces(Pstream::nProcs());
481 gatheredFaces[Pstream::myProcNo()] = allBoundary.localFaces();
482 forAll(gatheredFaces[Pstream::myProcNo()], i)
483 {
484 inplaceRenumber(pointToGlobal, gatheredFaces[Pstream::myProcNo()][i]);
485 }
486 Pstream::gatherList(gatheredFaces);
487
488 // Gather all ZoneIDs
489 List<labelList> gatheredZones(Pstream::nProcs());
490 gatheredZones[Pstream::myProcNo()].transfer(compactZones);
491 Pstream::gatherList(gatheredZones);
492
493 // On master combine all points, faces, zones
494 if (Pstream::master())
495 {
496 pointField allPoints = ListListOps::combine<pointField>
497 (
498 gatheredPoints,
500 );
501 gatheredPoints.clear();
502
503 faceList allFaces = ListListOps::combine<faceList>
504 (
505 gatheredFaces,
507 );
508 gatheredFaces.clear();
509
510 labelList allZones = ListListOps::combine<labelList>
511 (
512 gatheredZones,
514 );
515 gatheredZones.clear();
516
517
518 // Zones
519 surfZoneIdentifierList surfZones(compactZoneID.size());
520 forAllConstIters(compactZoneID, iter)
521 {
522 surfZones[iter()] = surfZoneIdentifier(iter.key(), iter());
523 Info<< "surfZone " << iter() << " : " << surfZones[iter()].name()
524 << endl;
525 }
526
527 UnsortedMeshedSurface<face> unsortedFace
528 (
529 std::move(allPoints),
530 std::move(allFaces),
531 std::move(allZones),
532 surfZones
533 );
534
535
536 MeshedSurface<face> sortedFace(unsortedFace);
537
538 fileName globalCasePath
539 (
540 runTime.processorCase()
541 ? runTime.globalPath()/outFileName
542 : runTime.path()/outFileName
543 );
544 globalCasePath.clean(); // Remove unneeded ".."
545
546 Info<< "Writing merged surface to " << globalCasePath << endl;
547
548 sortedFace.write(globalCasePath);
549 }
550}
551
552
553label checkAlignment(const polyMesh& mesh, const scalar tol, Ostream& os)
554{
555 // Check all edges aligned with one of the coordinate axes
556 const faceList& faces = mesh.faces();
557 const pointField& points = mesh.points();
558
559 label nUnaligned = 0;
560
561 forAll(faces, facei)
562 {
563 const face& f = faces[facei];
564 forAll(f, fp)
565 {
566 label fp1 = f.fcIndex(fp);
567 const linePointRef e(edge(f[fp], f[fp1]).line(points));
568 const vector v(e.vec());
569 const scalar magV(mag(v));
570 if (magV > ROOTVSMALL)
571 {
572 for
573 (
574 direction dir = 0;
575 dir < pTraits<vector>::nComponents;
576 ++dir
577 )
578 {
579 const scalar s(mag(v[dir]));
580 if (s > magV*tol && s < magV*(1-tol))
581 {
582 ++nUnaligned;
583 break;
584 }
585 }
586 }
587 }
588 }
589
590 reduce(nUnaligned, sumOp<label>());
591
592 if (nUnaligned)
593 {
594 os << "Initial mesh has " << nUnaligned
595 << " edges unaligned with any of the coordinate axes" << nl << endl;
596 }
597 return nUnaligned;
598}
599
600
601// Check writing tolerance before doing any serious work
602scalar getMergeDistance
603(
604 const polyMesh& mesh,
605 const scalar mergeTol,
606 const bool dryRun
607)
608{
609 const boundBox& meshBb = mesh.bounds();
610 scalar mergeDist = mergeTol * meshBb.mag();
611
612 Info<< nl
613 << "Overall mesh bounding box : " << meshBb << nl
614 << "Relative tolerance : " << mergeTol << nl
615 << "Absolute matching distance : " << mergeDist << nl
616 << endl;
617
618 // check writing tolerance
619 if (mesh.time().writeFormat() == IOstream::ASCII && !dryRun)
620 {
621 const scalar writeTol = std::pow
622 (
623 scalar(10),
624 -scalar(IOstream::defaultPrecision())
625 );
626
627 if (mergeTol < writeTol)
628 {
630 << "Your current settings specify ASCII writing with "
631 << IOstream::defaultPrecision() << " digits precision." << nl
632 << "Your merging tolerance (" << mergeTol
633 << ") is finer than this." << nl
634 << "Change to binary writeFormat, "
635 << "or increase the writePrecision" << endl
636 << "or adjust the merge tolerance (mergeTol)."
637 << exit(FatalError);
638 }
639 }
640
641 return mergeDist;
642}
643
644
645void removeZeroSizedPatches(fvMesh& mesh)
646{
647 // Remove any zero-sized ones. Assumes
648 // - processor patches are already only there if needed
649 // - all other patches are available on all processors
650 // - but coupled ones might still be needed, even if zero-size
651 // (e.g. processorCyclic)
652 // See also logic in createPatch.
653 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
654
655 labelList oldToNew(pbm.size(), -1);
656 label newPatchi = 0;
657 forAll(pbm, patchi)
658 {
659 const polyPatch& pp = pbm[patchi];
660
661 if (!isA<processorPolyPatch>(pp))
662 {
663 if
664 (
665 isA<coupledPolyPatch>(pp)
666 || returnReduce(pp.size(), sumOp<label>())
667 )
668 {
669 // Coupled (and unknown size) or uncoupled and used
670 oldToNew[patchi] = newPatchi++;
671 }
672 }
673 }
674
675 forAll(pbm, patchi)
676 {
677 const polyPatch& pp = pbm[patchi];
678
679 if (isA<processorPolyPatch>(pp))
680 {
681 oldToNew[patchi] = newPatchi++;
682 }
683 }
684
685
686 const label nKeepPatches = newPatchi;
687
688 // Shuffle unused ones to end
689 if (nKeepPatches != pbm.size())
690 {
691 Info<< endl
692 << "Removing zero-sized patches:" << endl << incrIndent;
693
694 forAll(oldToNew, patchi)
695 {
696 if (oldToNew[patchi] == -1)
697 {
698 Info<< indent << pbm[patchi].name()
699 << " type " << pbm[patchi].type()
700 << " at position " << patchi << endl;
701 oldToNew[patchi] = newPatchi++;
702 }
703 }
705
706 fvMeshTools::reorderPatches(mesh, oldToNew, nKeepPatches, true);
707 Info<< endl;
708 }
709}
710
711
712// Write mesh and additional information
713void writeMesh
714(
715 const string& msg,
716 const meshRefinement& meshRefiner,
717 const meshRefinement::debugType debugLevel,
718 const meshRefinement::writeType writeLevel
719)
720{
721 const fvMesh& mesh = meshRefiner.mesh();
722
723 meshRefiner.printMeshInfo(debugLevel, msg);
724 Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
725
726 processorMeshes::removeFiles(mesh);
727 if (!debugLevel && !(writeLevel&meshRefinement::WRITELAYERSETS))
728 {
729 topoSet::removeFiles(mesh);
730 }
731 refinementHistory::removeFiles(mesh);
732
733 meshRefiner.write
734 (
735 debugLevel,
736 meshRefinement::writeType(writeLevel | meshRefinement::WRITEMESH),
737 mesh.time().path()/meshRefiner.timeName()
738 );
739 Info<< "Wrote mesh in = "
740 << mesh.time().cpuTimeIncrement() << " s." << endl;
741}
742
743
744int main(int argc, char *argv[])
745{
746 argList::addNote
747 (
748 "Automatic split hex mesher. Refines and snaps to surface"
749 );
750
751 #include "addRegionOption.H"
752 #include "addOverwriteOption.H"
753 #include "addProfilingOption.H"
754 argList::addBoolOption
755 (
756 "checkGeometry",
757 "Check all surface geometry for quality"
758 );
759 argList::addDryRunOption
760 (
761 "Check case set-up only using a single time step"
762 );
763 argList::addOption
764 (
765 "surfaceSimplify",
766 "boundBox",
767 "Simplify the surface using snappyHexMesh starting from a boundBox"
768 );
769 argList::addOption
770 (
771 "patches",
772 "(patch0 .. patchN)",
773 "Only triangulate selected patches (wildcards supported)"
774 );
775 argList::addOption
776 (
777 "outFile",
778 "file",
779 "Name of the file to save the simplified surface to"
780 );
781 argList::addOption("dict", "file", "Alternative snappyHexMeshDict");
782
783 argList::noFunctionObjects(); // Never use function objects
784
785 #include "setRootCase.H"
786 #include "createTime.H"
787
788 const bool overwrite = args.found("overwrite");
789 const bool checkGeometry = args.found("checkGeometry");
790 const bool surfaceSimplify = args.found("surfaceSimplify");
791 const bool dryRun = args.dryRun();
792
793 if (dryRun)
794 {
795 Info<< "Operating in dry-run mode to detect set-up errors"
796 << nl << endl;
797 }
798
799
800 #include "createNamedMesh.H"
801 Info<< "Read mesh in = "
802 << runTime.cpuTimeIncrement() << " s" << endl;
803
804 // Check patches and faceZones are synchronised
805 mesh.boundaryMesh().checkParallelSync(true);
806 meshRefinement::checkCoupledFaceZones(mesh);
807
808 if (dryRun)
809 {
810 // Check if mesh aligned with cartesian axes
811 checkAlignment(mesh, 1e-6, Pout); //FatalIOError);
812 }
813
814
815
816 // Read meshing dictionary
817 const word dictName("snappyHexMeshDict");
820
821
822 // all surface geometry
823 const dictionary& geometryDict =
824 meshRefinement::subDict(meshDict, "geometry", dryRun);
825
826 // refinement parameters
827 const dictionary& refineDict =
828 meshRefinement::subDict(meshDict, "castellatedMeshControls", dryRun);
829
830 // mesh motion and mesh quality parameters
831 const dictionary& motionDict =
832 meshRefinement::subDict(meshDict, "meshQualityControls", dryRun);
833
834 // snap-to-surface parameters
835 const dictionary& snapDict =
836 meshRefinement::subDict(meshDict, "snapControls", dryRun);
837
838 // layer addition parameters
839 const dictionary& layerDict =
840 meshRefinement::subDict(meshDict, "addLayersControls", dryRun);
841
842 // absolute merge distance
843 const scalar mergeDist = getMergeDistance
844 (
845 mesh,
846 meshRefinement::get<scalar>
847 (
848 meshDict,
849 "mergeTolerance",
850 dryRun
851 ),
852 dryRun
853 );
854
855 const bool keepPatches(meshDict.getOrDefault("keepPatches", false));
856
857 // Writer for writing lines
858 autoPtr<coordSetWriter> setFormatter;
859 {
860 const word setFormat
861 (
862 meshDict.getOrDefault<word>
863 (
864 "setFormat",
865 coordSetWriters::vtkWriter::typeName // Default: "vtk"
866 )
867 );
868
869 setFormatter = coordSetWriter::New
870 (
871 setFormat,
872 meshDict.subOrEmptyDict("formatOptions").optionalSubDict(setFormat)
873 );
874 }
875
876
877 const scalar maxSizeRatio
878 (
879 meshDict.getOrDefault<scalar>("maxSizeRatio", 100)
880 );
881
882
883 // Read decomposePar dictionary
884 dictionary decomposeDict;
885 if (Pstream::parRun())
886 {
887 // Ensure demand-driven decompositionMethod finds alternative
888 // decomposeParDict location properly.
889
890 IOdictionary* dictPtr = new IOdictionary
891 (
892 IOobject::selectIO
893 (
895 (
896 decompositionModel::canonicalName,
897 runTime.system(),
898 runTime,
899 IOobject::MUST_READ,
900 IOobject::NO_WRITE
901 ),
902 args.getOrDefault<fileName>("decomposeParDict", "")
903 )
904 );
905
906 // Store it on the object registry, but to be found it must also
907 // have the expected "decomposeParDict" name.
908
909 dictPtr->rename(decompositionModel::canonicalName);
910 runTime.store(dictPtr);
911
912 decomposeDict = *dictPtr;
913 }
914 else
915 {
916 decomposeDict.add("method", "none");
917 decomposeDict.add("numberOfSubdomains", 1);
918 }
919
920
921 // Debug
922 // ~~~~~
923
924 // Set debug level
926 (
927 meshDict.getOrDefault<label>
928 (
929 "debug",
930 0
931 )
932 );
933 {
934 wordList flags;
935 if (meshDict.readIfPresent("debugFlags", flags))
936 {
937 debugLevel = meshRefinement::debugType
938 (
939 meshRefinement::readFlags
940 (
941 meshRefinement::debugTypeNames,
942 flags
943 )
944 );
945 }
946 }
947 if (debugLevel > 0)
948 {
949 meshRefinement::debug = debugLevel;
950 snappyRefineDriver::debug = debugLevel;
951 snappySnapDriver::debug = debugLevel;
952 snappyLayerDriver::debug = debugLevel;
953 }
954
955 // Set file writing level
956 {
957 wordList flags;
958 if (meshDict.readIfPresent("writeFlags", flags))
959 {
960 meshRefinement::writeLevel
961 (
963 (
964 meshRefinement::readFlags
965 (
966 meshRefinement::writeTypeNames,
967 flags
968 )
969 )
970 );
971 }
972 }
973
975 //{
976 // wordList flags;
977 // if (meshDict.readIfPresent("outputFlags", flags))
978 // {
979 // meshRefinement::outputLevel
980 // (
981 // meshRefinement::outputType
982 // (
983 // meshRefinement::readFlags
984 // (
985 // meshRefinement::outputTypeNames,
986 // flags
987 // )
988 // )
989 // );
990 // }
991 //}
992
993 // for the impatient who want to see some output files:
994 profiling::writeNow();
995
996 // Read geometry
997 // ~~~~~~~~~~~~~
998
999 searchableSurfaces allGeometry
1000 (
1001 IOobject
1002 (
1003 "abc", // dummy name
1004 mesh.time().constant(), // instance
1005 //mesh.time().findInstance("triSurface", word::null),// instance
1006 "triSurface", // local
1007 mesh.time(), // registry
1008 IOobject::MUST_READ,
1009 IOobject::NO_WRITE
1010 ),
1011 geometryDict,
1012 meshDict.getOrDefault("singleRegionName", true)
1013 );
1014
1015
1016 // Read refinement surfaces
1017 // ~~~~~~~~~~~~~~~~~~~~~~~~
1018
1019 autoPtr<refinementSurfaces> surfacesPtr;
1020
1021 Info<< "Reading refinement surfaces." << endl;
1022
1023 if (surfaceSimplify)
1024 {
1025 addProfiling(surfaceSimplify, "snappyHexMesh::surfaceSimplify");
1026 IOdictionary foamyHexMeshDict
1027 (
1028 IOobject
1029 (
1030 "foamyHexMeshDict",
1031 runTime.system(),
1032 runTime,
1033 IOobject::MUST_READ_IF_MODIFIED,
1034 IOobject::NO_WRITE
1035 )
1036 );
1037
1038 const dictionary& conformationDict =
1039 foamyHexMeshDict.subDict("surfaceConformation").subDict
1040 (
1041 "geometryToConformTo"
1042 );
1043
1044 const dictionary& motionDict =
1045 foamyHexMeshDict.subDict("motionControl");
1046
1047 const dictionary& shapeControlDict =
1048 motionDict.subDict("shapeControlFunctions");
1049
1050 // Calculate current ratio of hex cells v.s. wanted cell size
1051 const scalar defaultCellSize =
1052 motionDict.get<scalar>("defaultCellSize");
1053
1054 const scalar initialCellSize = ::pow(mesh.V()[0], 1.0/3.0);
1055
1056 //Info<< "Wanted cell size = " << defaultCellSize << endl;
1057 //Info<< "Current cell size = " << initialCellSize << endl;
1058 //Info<< "Fraction = " << initialCellSize/defaultCellSize
1059 // << endl;
1060
1061 surfacesPtr =
1062 createRefinementSurfaces
1063 (
1064 allGeometry,
1065 conformationDict,
1066 shapeControlDict,
1067 refineDict.getOrDefault("gapLevelIncrement", 0),
1068 initialCellSize/defaultCellSize
1069 );
1070
1071 profiling::writeNow();
1072 }
1073 else
1074 {
1075 surfacesPtr.reset
1076 (
1078 (
1079 allGeometry,
1080 meshRefinement::subDict
1081 (
1082 refineDict,
1083 "refinementSurfaces",
1084 dryRun
1085 ),
1086 refineDict.getOrDefault("gapLevelIncrement", 0),
1087 dryRun
1088 )
1089 );
1090
1091 Info<< "Read refinement surfaces in = "
1092 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1093 }
1094
1095 refinementSurfaces& surfaces = surfacesPtr();
1096
1097
1098 // Checking only?
1099
1100 if (checkGeometry)
1101 {
1102 // Check geometry amongst itself (e.g. intersection, size differences)
1103
1104 // Extract patchInfo
1105 List<wordList> patchTypes(allGeometry.size());
1106
1107 const PtrList<dictionary>& patchInfo = surfaces.patchInfo();
1108 const labelList& surfaceGeometry = surfaces.surfaces();
1109 forAll(surfaceGeometry, surfi)
1110 {
1111 label geomi = surfaceGeometry[surfi];
1112 const wordList& regNames = allGeometry.regionNames()[geomi];
1113
1114 patchTypes[geomi].setSize(regNames.size());
1115 forAll(regNames, regioni)
1116 {
1117 label globalRegioni = surfaces.globalRegion(surfi, regioni);
1118
1119 if (patchInfo.set(globalRegioni))
1120 {
1121 patchTypes[geomi][regioni] =
1122 meshRefinement::get<word>
1123 (
1124 patchInfo[globalRegioni],
1125 "type",
1126 dryRun,
1127 keyType::REGEX,
1129 );
1130 }
1131 else
1132 {
1133 patchTypes[geomi][regioni] = wallPolyPatch::typeName;
1134 }
1135 }
1136 }
1137
1138 // Write some stats
1139 allGeometry.writeStats(patchTypes, Info);
1140 // Check topology
1141 allGeometry.checkTopology(true);
1142 // Check geometry
1143 allGeometry.checkGeometry
1144 (
1145 maxSizeRatio, // max size ratio
1146 1e-9, // intersection tolerance
1147 setFormatter,
1148 0.01, // min triangle quality
1149 true
1150 );
1151
1152 if (!dryRun)
1153 {
1154 return 0;
1155 }
1156 }
1157
1158
1159 if (dryRun)
1160 {
1161 // Check geometry to mesh bounding box
1162 Info<< "Checking for geometry size relative to mesh." << endl;
1163 const boundBox& meshBb = mesh.bounds();
1164 forAll(allGeometry, geomi)
1165 {
1166 const searchableSurface& s = allGeometry[geomi];
1167 const boundBox& bb = s.bounds();
1168
1169 scalar ratio = bb.mag() / meshBb.mag();
1170 if (ratio > maxSizeRatio || ratio < 1.0/maxSizeRatio)
1171 {
1172 Warning
1173 << " " << allGeometry.names()[geomi]
1174 << " bounds differ from mesh"
1175 << " by more than a factor " << maxSizeRatio << ":" << nl
1176 << " bounding box : " << bb << nl
1177 << " mesh bounding box : " << meshBb
1178 << endl;
1179 }
1180 if (!meshBb.contains(bb))
1181 {
1182 Warning
1183 << " " << allGeometry.names()[geomi]
1184 << " bounds not fully contained in mesh" << nl
1185 << " bounding box : " << bb << nl
1186 << " mesh bounding box : " << meshBb
1187 << endl;
1188 }
1189 }
1190 Info<< endl;
1191 }
1192
1193
1194
1195
1196 // Read refinement shells
1197 // ~~~~~~~~~~~~~~~~~~~~~~
1198
1199 Info<< "Reading refinement shells." << endl;
1200 shellSurfaces shells
1201 (
1202 allGeometry,
1203 meshRefinement::subDict(refineDict, "refinementRegions", dryRun),
1204 dryRun
1205 );
1206 Info<< "Read refinement shells in = "
1207 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1208
1209
1210 Info<< "Setting refinement level of surface to be consistent"
1211 << " with shells." << endl;
1212 surfaces.setMinLevelFields(shells);
1213 Info<< "Checked shell refinement in = "
1214 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1215
1216
1217 // Optionally read limit shells
1218 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1219
1220 const dictionary limitDict(refineDict.subOrEmptyDict("limitRegions"));
1221
1222 if (!limitDict.empty())
1223 {
1224 Info<< "Reading limit shells." << endl;
1225 }
1226
1227 shellSurfaces limitShells(allGeometry, limitDict, dryRun);
1228
1229 if (!limitDict.empty())
1230 {
1231 Info<< "Read limit shells in = "
1232 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1233 }
1234
1235 if (dryRun)
1236 {
1237 // Check for use of all geometry
1238 const wordList& allGeomNames = allGeometry.names();
1239
1240 labelHashSet unusedGeometries(identity(allGeomNames.size()));
1241 unusedGeometries.erase(surfaces.surfaces());
1242 unusedGeometries.erase(shells.shells());
1243 unusedGeometries.erase(limitShells.shells());
1244
1245 if (unusedGeometries.size())
1246 {
1247 IOWarningInFunction(geometryDict)
1248 << "The following geometry entries are not used:" << nl;
1249 for (const label geomi : unusedGeometries)
1250 {
1251 Info<< " " << allGeomNames[geomi] << nl;
1252 }
1253 Info<< endl;
1254 }
1255 }
1256
1257
1258
1259
1260 // Read feature meshes
1261 // ~~~~~~~~~~~~~~~~~~~
1262
1263 Info<< "Reading features." << endl;
1264 refinementFeatures features
1265 (
1266 mesh,
1268 (
1269 meshRefinement::lookup(refineDict, "features", dryRun)
1270 ),
1271 dryRun
1272 );
1273 Info<< "Read features in = "
1274 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1275
1276
1277 if (dryRun)
1278 {
1279 // Check geometry to mesh bounding box
1280 Info<< "Checking for line geometry size relative to surface geometry."
1281 << endl;
1282
1284 bool hasErrors = features.checkSizes
1285 (
1286 maxSizeRatio, //const scalar maxRatio,
1287 mesh.bounds(),
1288 true, //const bool report,
1289 os //FatalIOError
1290 );
1291 if (hasErrors)
1292 {
1293 Warning<< os.str() << endl;
1294 }
1295 }
1296
1297
1298 // Refinement engine
1299 // ~~~~~~~~~~~~~~~~~
1300
1301 Info<< nl
1302 << "Determining initial surface intersections" << nl
1303 << "-----------------------------------------" << nl
1304 << endl;
1305
1306 // Main refinement engine
1307 meshRefinement meshRefiner
1308 (
1309 mesh,
1310 mergeDist, // tolerance used in sorting coordinates
1311 overwrite, // overwrite mesh files?
1312 surfaces, // for surface intersection refinement
1313 features, // for feature edges/point based refinement
1314 shells, // for volume (inside/outside) refinement
1315 limitShells, // limit of volume refinement
1316 labelList(), // initial faces to test
1317 dryRun
1318 );
1319
1320 if (!dryRun)
1321 {
1322 meshRefiner.updateIntersections(identity(mesh.nFaces()));
1323 Info<< "Calculated surface intersections in = "
1324 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1325 }
1326
1327 // Some stats
1328 meshRefiner.printMeshInfo(debugLevel, "Initial mesh");
1329
1330 meshRefiner.write
1331 (
1332 meshRefinement::debugType(debugLevel&meshRefinement::OBJINTERSECTIONS),
1334 mesh.time().path()/meshRefiner.timeName()
1335 );
1336
1337
1338 // Refinement parameters
1339 const refinementParameters refineParams(refineDict, dryRun);
1340
1341 // Snap parameters
1342 const snapParameters snapParams(snapDict, dryRun);
1343
1344
1345
1346 // Add all the cellZones and faceZones
1347 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1348
1349 // 1. cellZones relating to surface (faceZones added later)
1350
1351 const labelList namedSurfaces
1352 (
1353 surfaceZonesInfo::getNamedSurfaces(surfaces.surfZones())
1354 );
1355
1356 labelList surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
1357 (
1358 surfaces.surfZones(),
1359 namedSurfaces,
1360 mesh
1361 );
1362
1363
1364 // 2. cellZones relating to locations
1365
1366 refineParams.addCellZonesToMesh(mesh);
1367
1368
1369
1370 // Add all the surface regions as patches
1371 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1372
1373 //- Global surface region to patch (non faceZone surface) or patches
1374 // (faceZone surfaces)
1375 labelList globalToMasterPatch;
1376 labelList globalToSlavePatch;
1377
1378
1379 {
1380 Info<< nl
1381 << "Adding patches for surface regions" << nl
1382 << "----------------------------------" << nl
1383 << endl;
1384
1385 // From global region number to mesh patch.
1386 globalToMasterPatch.setSize(surfaces.nRegions(), -1);
1387 globalToSlavePatch.setSize(surfaces.nRegions(), -1);
1388
1389 if (!dryRun)
1390 {
1391 Info<< setf(ios_base::left)
1392 << setw(6) << "Patch"
1393 << setw(20) << "Type"
1394 << setw(30) << "Region" << nl
1395 << setw(6) << "-----"
1396 << setw(20) << "----"
1397 << setw(30) << "------" << endl;
1398 }
1399
1400 const labelList& surfaceGeometry = surfaces.surfaces();
1401 const PtrList<dictionary>& surfacePatchInfo = surfaces.patchInfo();
1402 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
1403
1404 forAll(surfaceGeometry, surfi)
1405 {
1406 label geomi = surfaceGeometry[surfi];
1407
1408 const wordList& regNames = allGeometry.regionNames()[geomi];
1409
1410 if (!dryRun)
1411 {
1412 Info<< surfaces.names()[surfi] << ':' << nl << nl;
1413 }
1414
1415 const wordList& fzNames =
1416 surfaces.surfZones()[surfi].faceZoneNames();
1417
1418 if (fzNames.empty())
1419 {
1420 // 'Normal' surface
1421 forAll(regNames, i)
1422 {
1423 label globalRegioni = surfaces.globalRegion(surfi, i);
1424
1425 label patchi;
1426
1427 if (surfacePatchInfo.set(globalRegioni))
1428 {
1429 patchi = meshRefiner.addMeshedPatch
1430 (
1431 regNames[i],
1432 surfacePatchInfo[globalRegioni]
1433 );
1434 }
1435 else
1436 {
1437 dictionary patchInfo;
1438 patchInfo.set("type", wallPolyPatch::typeName);
1439
1440 patchi = meshRefiner.addMeshedPatch
1441 (
1442 regNames[i],
1443 patchInfo
1444 );
1445 }
1446
1447 if (!dryRun)
1448 {
1449 Info<< setf(ios_base::left)
1450 << setw(6) << patchi
1451 << setw(20) << pbm[patchi].type()
1452 << setw(30) << regNames[i] << nl;
1453 }
1454
1455 globalToMasterPatch[globalRegioni] = patchi;
1456 globalToSlavePatch[globalRegioni] = patchi;
1457 }
1458 }
1459 else
1460 {
1461 // Zoned surface
1462 forAll(regNames, i)
1463 {
1464 label globalRegioni = surfaces.globalRegion(surfi, i);
1465
1466 // Add master side patch
1467 {
1468 label patchi;
1469
1470 if (surfacePatchInfo.set(globalRegioni))
1471 {
1472 patchi = meshRefiner.addMeshedPatch
1473 (
1474 regNames[i],
1475 surfacePatchInfo[globalRegioni]
1476 );
1477 }
1478 else
1479 {
1480 dictionary patchInfo;
1481 patchInfo.set("type", wallPolyPatch::typeName);
1482
1483 patchi = meshRefiner.addMeshedPatch
1484 (
1485 regNames[i],
1486 patchInfo
1487 );
1488 }
1489
1490 if (!dryRun)
1491 {
1492 Info<< setf(ios_base::left)
1493 << setw(6) << patchi
1494 << setw(20) << pbm[patchi].type()
1495 << setw(30) << regNames[i] << nl;
1496 }
1497
1498 globalToMasterPatch[globalRegioni] = patchi;
1499 }
1500 // Add slave side patch
1501 {
1502 const word slaveName = regNames[i] + "_slave";
1503 label patchi;
1504
1505 if (surfacePatchInfo.set(globalRegioni))
1506 {
1507 patchi = meshRefiner.addMeshedPatch
1508 (
1509 slaveName,
1510 surfacePatchInfo[globalRegioni]
1511 );
1512 }
1513 else
1514 {
1515 dictionary patchInfo;
1516 patchInfo.set("type", wallPolyPatch::typeName);
1517
1518 patchi = meshRefiner.addMeshedPatch
1519 (
1520 slaveName,
1521 patchInfo
1522 );
1523 }
1524
1525 if (!dryRun)
1526 {
1527 Info<< setf(ios_base::left)
1528 << setw(6) << patchi
1529 << setw(20) << pbm[patchi].type()
1530 << setw(30) << slaveName << nl;
1531 }
1532
1533 globalToSlavePatch[globalRegioni] = patchi;
1534 }
1535 }
1536
1537 // For now: have single faceZone per surface. Use first
1538 // region in surface for patch for zoning
1539 if (regNames.size())
1540 {
1541 forAll(fzNames, fzi)
1542 {
1543 const word& fzName = fzNames[fzi];
1544 label globalRegioni = surfaces.globalRegion(surfi, fzi);
1545
1546 meshRefiner.addFaceZone
1547 (
1548 fzName,
1549 pbm[globalToMasterPatch[globalRegioni]].name(),
1550 pbm[globalToSlavePatch[globalRegioni]].name(),
1551 surfaces.surfZones()[surfi].faceType()
1552 );
1553 }
1554 }
1555 }
1556
1557 if (!dryRun)
1558 {
1559 Info<< nl;
1560 }
1561 }
1562 Info<< "Added patches in = "
1563 << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
1564 }
1565
1566
1567
1568 // Add all information for all the remaining faceZones
1569 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1570
1571 HashTable<Pair<word>> faceZoneToPatches;
1572 forAll(mesh.faceZones(), zonei)
1573 {
1574 const word& fzName = mesh.faceZones()[zonei].name();
1575
1576 label mpI, spI;
1578 bool hasInfo = meshRefiner.getFaceZoneInfo(fzName, mpI, spI, fzType);
1579
1580 if (!hasInfo)
1581 {
1582 // faceZone does not originate from a surface but presumably
1583 // from a cellZone pair instead
1584 string::size_type i = fzName.find("_to_");
1585 if (i != string::npos)
1586 {
1587 word cz0 = fzName.substr(0, i);
1588 word cz1 = fzName.substr(i+4, fzName.size()-i+4);
1589 word slaveName(cz1 + "_to_" + cz0);
1590 faceZoneToPatches.insert(fzName, Pair<word>(fzName, slaveName));
1591 }
1592 else
1593 {
1594 // Add as fzName + fzName_slave
1595 const word slaveName = fzName + "_slave";
1596 faceZoneToPatches.insert(fzName, Pair<word>(fzName, slaveName));
1597 }
1598 }
1599 }
1600
1601 if (faceZoneToPatches.size())
1602 {
1603 snappyRefineDriver::addFaceZones
1604 (
1605 meshRefiner,
1606 refineParams,
1607 faceZoneToPatches
1608 );
1609 }
1610
1611
1612
1613 // Re-do intersections on meshed boundaries since they use an extrapolated
1614 // other side
1615 {
1616 const labelList adaptPatchIDs(meshRefiner.meshedPatches());
1617
1618 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
1619
1620 label nFaces = 0;
1621 forAll(adaptPatchIDs, i)
1622 {
1623 nFaces += pbm[adaptPatchIDs[i]].size();
1624 }
1625
1626 labelList faceLabels(nFaces);
1627 nFaces = 0;
1628 forAll(adaptPatchIDs, i)
1629 {
1630 const polyPatch& pp = pbm[adaptPatchIDs[i]];
1631 forAll(pp, i)
1632 {
1633 faceLabels[nFaces++] = pp.start()+i;
1634 }
1635 }
1636 meshRefiner.updateIntersections(faceLabels);
1637 }
1638
1639
1640
1641 // Parallel
1642 // ~~~~~~~~
1643
1644 // Construct decomposition engine. Note: cannot use decompositionModel
1645 // MeshObject since we're clearing out the mesh inside the mesh generation.
1646 autoPtr<decompositionMethod> decomposerPtr
1647 (
1648 decompositionMethod::New
1649 (
1650 decomposeDict
1651 )
1652 );
1653 decompositionMethod& decomposer = *decomposerPtr;
1654
1655 if (Pstream::parRun() && !decomposer.parallelAware())
1656 {
1658 << "You have selected decomposition method "
1659 << decomposer.typeName
1660 << " which is not parallel aware." << endl
1661 << "Please select one that is (hierarchical, ptscotch)"
1662 << exit(FatalError);
1663 }
1664
1665 // Mesh distribution engine (uses tolerance to reconstruct meshes)
1666 fvMeshDistribute distributor(mesh);
1667
1668
1669
1670
1671
1672 // Now do the real work -refinement -snapping -layers
1673 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1674
1675 const bool wantRefine
1676 (
1677 meshRefinement::get<bool>(meshDict, "castellatedMesh", dryRun)
1678 );
1679 const bool wantSnap
1680 (
1681 meshRefinement::get<bool>(meshDict, "snap", dryRun)
1682 );
1683 const bool wantLayers
1684 (
1685 meshRefinement::get<bool>(meshDict, "addLayers", dryRun)
1686 );
1687
1688 if (dryRun)
1689 {
1690 string errorMsg(FatalError.message());
1691 string IOerrorMsg(FatalIOError.message());
1692
1693 if (errorMsg.size() || IOerrorMsg.size())
1694 {
1695 //errorMsg = "[dryRun] " + errorMsg;
1696 //errorMsg.replaceAll("\n", "\n[dryRun] ");
1697 //IOerrorMsg = "[dryRun] " + IOerrorMsg;
1698 //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
1699
1700 Warning
1701 << nl
1702 << "Missing/incorrect required dictionary entries:" << nl
1703 << nl
1704 << IOerrorMsg.c_str() << nl
1705 << errorMsg.c_str() << nl << nl
1706 << "Exiting dry-run" << nl << endl;
1707
1708 FatalError.clear();
1710
1711 return 0;
1712 }
1713 }
1714
1715
1716 // How to treat co-planar faces
1718 meshRefinement::FaceMergeType::GEOMETRIC;
1719 {
1720 const bool mergePatchFaces
1721 (
1722 meshDict.getOrDefault("mergePatchFaces", true)
1723 );
1724
1725 if (!mergePatchFaces)
1726 {
1727 Info<< "Not merging patch-faces of cell to preserve"
1728 << " (split)hex cell shape."
1729 << nl << endl;
1730 mergeType = meshRefinement::FaceMergeType::NONE;
1731 }
1732 else
1733 {
1734 const bool mergeAcrossPatches
1735 (
1736 meshDict.getOrDefault("mergeAcrossPatches", false)
1737 );
1738
1739 if (mergeAcrossPatches)
1740 {
1741 Info<< "Merging co-planar patch-faces of cells"
1742 << ", regardless of patch assignment"
1743 << nl << endl;
1744 mergeType = meshRefinement::FaceMergeType::IGNOREPATCH;
1745 }
1746 }
1747 }
1748
1749
1750
1751 if (wantRefine)
1752 {
1753 cpuTime timer;
1754
1755 snappyRefineDriver refineDriver
1756 (
1757 meshRefiner,
1758 decomposer,
1759 distributor,
1760 globalToMasterPatch,
1761 globalToSlavePatch,
1762 setFormatter(),
1763 dryRun
1764 );
1765
1766
1767 if (!overwrite && !debugLevel)
1768 {
1769 const_cast<Time&>(mesh.time())++;
1770 }
1771
1772
1773 refineDriver.doRefine
1774 (
1775 refineDict,
1776 refineParams,
1777 snapParams,
1778 refineParams.handleSnapProblems(),
1779 mergeType,
1780 motionDict
1781 );
1782
1783 // Remove zero sized patches originating from faceZones
1784 if (!keepPatches && !wantSnap && !wantLayers)
1785 {
1786 fvMeshTools::removeEmptyPatches(mesh, true);
1787 }
1788
1789 if (!dryRun)
1790 {
1791 writeMesh
1792 (
1793 "Refined mesh",
1794 meshRefiner,
1795 debugLevel,
1796 meshRefinement::writeLevel()
1797 );
1798 }
1799
1800 Info<< "Mesh refined in = "
1801 << timer.cpuTimeIncrement() << " s." << endl;
1802
1803 profiling::writeNow();
1804 }
1805
1806 if (wantSnap)
1807 {
1808 cpuTime timer;
1809
1810 snappySnapDriver snapDriver
1811 (
1812 meshRefiner,
1813 globalToMasterPatch,
1814 globalToSlavePatch,
1815 dryRun
1816 );
1817
1818 if (!overwrite && !debugLevel)
1819 {
1820 const_cast<Time&>(mesh.time())++;
1821 }
1822
1823 // Use the resolveFeatureAngle from the refinement parameters
1824 scalar curvature = refineParams.curvature();
1825 scalar planarAngle = refineParams.planarAngle();
1826
1827 snapDriver.doSnap
1828 (
1829 snapDict,
1830 motionDict,
1831 mergeType,
1832 curvature,
1833 planarAngle,
1834 snapParams
1835 );
1836
1837 // Remove zero sized patches originating from faceZones
1838 if (!keepPatches && !wantLayers)
1839 {
1840 fvMeshTools::removeEmptyPatches(mesh, true);
1841 }
1842
1843 if (!dryRun)
1844 {
1845 writeMesh
1846 (
1847 "Snapped mesh",
1848 meshRefiner,
1849 debugLevel,
1850 meshRefinement::writeLevel()
1851 );
1852 }
1853
1854 Info<< "Mesh snapped in = "
1855 << timer.cpuTimeIncrement() << " s." << endl;
1856
1857 profiling::writeNow();
1858 }
1859
1860 if (wantLayers)
1861 {
1862 cpuTime timer;
1863
1864 // Layer addition parameters
1865 const layerParameters layerParams
1866 (
1867 layerDict,
1868 mesh.boundaryMesh(),
1869 dryRun
1870 );
1871
1872 snappyLayerDriver layerDriver
1873 (
1874 meshRefiner,
1875 globalToMasterPatch,
1876 globalToSlavePatch,
1877 dryRun
1878 );
1879
1880 // Use the maxLocalCells from the refinement parameters
1881 bool preBalance = returnReduce
1882 (
1883 (mesh.nCells() >= refineParams.maxLocalCells()),
1884 orOp<bool>()
1885 );
1886
1887
1888 if (!overwrite && !debugLevel)
1889 {
1890 const_cast<Time&>(mesh.time())++;
1891 }
1892
1893 layerDriver.doLayers
1894 (
1895 layerDict,
1896 motionDict,
1897 layerParams,
1898 mergeType,
1899 preBalance,
1900 decomposer,
1901 distributor
1902 );
1903
1904 // Remove zero sized patches originating from faceZones
1905 if (!keepPatches)
1906 {
1907 fvMeshTools::removeEmptyPatches(mesh, true);
1908 }
1909
1910 if (!dryRun)
1911 {
1912 writeMesh
1913 (
1914 "Layer mesh",
1915 meshRefiner,
1916 debugLevel,
1917 meshRefinement::writeLevel()
1918 );
1919 }
1920
1921 Info<< "Layers added in = "
1922 << timer.cpuTimeIncrement() << " s." << endl;
1923
1924 profiling::writeNow();
1925 }
1926
1927
1928 {
1929 addProfiling(checkMesh, "snappyHexMesh::checkMesh");
1930
1931 // Check final mesh
1932 Info<< "Checking final mesh ..." << endl;
1933 faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
1934 motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun);
1935 const label nErrors = returnReduce
1936 (
1937 wrongFaces.size(),
1938 sumOp<label>()
1939 );
1940
1941 if (nErrors > 0)
1942 {
1943 Info<< "Finished meshing with " << nErrors << " illegal faces"
1944 << " (concave, zero area or negative cell pyramid volume)"
1945 << endl;
1946 wrongFaces.write();
1947 }
1948 else
1949 {
1950 Info<< "Finished meshing without any errors" << endl;
1951 }
1952
1953 profiling::writeNow();
1954 }
1955
1956
1957 if (surfaceSimplify)
1958 {
1959 addProfiling(surfaceSimplify, "snappyHexMesh::surfaceSimplify");
1960
1961 const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
1962
1963 labelHashSet includePatches(bMesh.size());
1964
1965 if (args.found("patches"))
1966 {
1967 includePatches = bMesh.patchSet
1968 (
1969 args.getList<wordRe>("patches")
1970 );
1971 }
1972 else
1973 {
1974 forAll(bMesh, patchi)
1975 {
1976 const polyPatch& patch = bMesh[patchi];
1977
1978 if (!isA<processorPolyPatch>(patch))
1979 {
1980 includePatches.insert(patchi);
1981 }
1982 }
1983 }
1984
1985 fileName outFileName
1986 (
1988 (
1989 "outFile",
1990 "constant/triSurface/simplifiedSurface.stl"
1991 )
1992 );
1993
1994 extractSurface
1995 (
1996 mesh,
1997 runTime,
1998 includePatches,
1999 outFileName
2000 );
2001
2002 pointIOField cellCentres
2003 (
2004 IOobject
2005 (
2006 "internalCellCentres",
2007 runTime.timeName(),
2008 mesh,
2009 IOobject::NO_READ,
2010 IOobject::AUTO_WRITE
2011 ),
2012 mesh.cellCentres()
2013 );
2014
2015 cellCentres.write();
2016 }
2017
2018 profiling::writeNow();
2019
2020 Info<< "Finished meshing in = "
2021 << runTime.elapsedCpuTime() << " s." << endl;
2022
2023
2024 if (dryRun)
2025 {
2026 string errorMsg(FatalError.message());
2027 string IOerrorMsg(FatalIOError.message());
2028
2029 if (errorMsg.size() || IOerrorMsg.size())
2030 {
2031 //errorMsg = "[dryRun] " + errorMsg;
2032 //errorMsg.replaceAll("\n", "\n[dryRun] ");
2033 //IOerrorMsg = "[dryRun] " + IOerrorMsg;
2034 //IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
2035
2036 Perr<< nl
2037 << "Missing/incorrect required dictionary entries:" << nl
2038 << nl
2039 << IOerrorMsg.c_str() << nl
2040 << errorMsg.c_str() << nl << nl
2041 << "Exiting dry-run" << nl << endl;
2042
2043 FatalError.clear();
2045
2046 return 0;
2047 }
2048 }
2049
2050
2051 Info<< "End\n" << endl;
2052
2053 return 0;
2054}
2055
2056
2057// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
reduce(hasMovingMesh, orOp< bool >())
label size_type
The type to represent the size of a buffer.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
bool erase(T *item)
Remove the specified element from the list and delete it.
Definition: ILList.C:92
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
A surface geometry mesh with zone information, not to be confused with the similarly named surfaceMes...
Definition: MeshedSurface.H:99
Output to string buffer, using a OSstream. Always UNCOMPRESSED.
Definition: StringStream.H:231
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
A list of faces which address into the list of points.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
A surface geometry mesh, in which the surface zone information is conveyed by the 'zoneId' associated...
int dryRun() const noexcept
Return the dry-run flag.
Definition: argListI.H:116
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
List< T > getList(const label index) const
Get a List of values from the argument at index.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:133
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTimeCxx.H:54
Abstract base class for domain decomposition.
virtual bool parallelAware() const =0
Is method parallel aware?
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
Definition: dictionary.C:540
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & optionalSubDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary, otherwise return this dictionary.
Definition: dictionary.C:577
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
autoPtr< dictionary > clone() const
Construct and return clone.
Definition: dictionary.C:172
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
entry * findEntry(const word &keyword, enum keyType::option matchOpt=keyType::REGEX)
Find for an entry (non-const access) with the given keyword.
Definition: dictionaryI.H:97
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
wordList toc() const
Return the table of contents.
Definition: dictionary.C:602
entry * set(entry *entryPtr)
Assign a new entry, overwriting any existing entry.
Definition: dictionary.C:780
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:70
const keyType & keyword() const noexcept
Return keyword.
Definition: entry.H:195
virtual const dictionary & dict() const =0
Return dictionary, if entry is a dictionary.
string message() const
The accumulated error message.
Definition: error.C:319
void clear() const
Clear any messages.
Definition: error.C:325
A list of face labels.
Definition: faceSet.H:54
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Simple container to keep together layer specific information.
A line primitive.
Definition: line.H:68
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
label addFaceZone(const word &fzName, const word &masterPatch, const word &slavePatch, const surfaceZonesInfo::faceZoneType &fzType)
Add/lookup faceZone and update information. Return index of.
bool getFaceZoneInfo(const word &fzName, label &masterPatchID, label &slavePatchID, surfaceZonesInfo::faceZoneType &fzType) const
Lookup faceZone information. Return false if no information.
void updateIntersections(const labelList &changedFaces)
Find any intersection of surface. Store in surfaceIndex_.
writeType
Enumeration for what to write. Used as a bit-pattern.
void printMeshInfo(const bool, const string &) const
Print some mesh stats.
debugType
Enumeration for what to debug. Used as a bit-pattern.
labelList meshedPatches() const
Get patchIDs for patches added in addMeshedPatch.
label addMeshedPatch(const word &name, const dictionary &)
Add patch originating from meshing. Update meshedPatches_.
const fvMesh & mesh() const
Reference to mesh.
FaceMergeType
Enumeration for what to do with co-planar patch faces on a single.
bool write() const
Write mesh and all data.
label index() const noexcept
The index of this patch in the boundaryMesh.
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Lookup type of boundary radiation properties.
Definition: lookup.H:66
Encapsulates queries for features.
Simple container to keep together refinement specific information.
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
const labelList & gapLevel() const
From global region number to small gap refinement level.
const PtrList< dictionary > & patchInfo() const
From global region number to patch type.
const wordList & names() const
Names of surfaces.
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
const labelList & surfaces() const
const labelList & minLevel() const
From global region number to refinement level.
const labelList & maxLevel() const
From global region number to refinement level.
void setMinLevelFields(const shellSurfaces &shells)
Calculate minLevelFields.
const PtrList< surfaceZonesInfo > & surfZones() const
virtual void rename(const word &newName)
Rename.
Definition: regIOobject.C:417
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Container for searchableSurfaces. The collection is specified as a dictionary. For example,...
const wordList & names() const
Surface names, not region names.
label checkTopology(const bool report) const
All topological checks. Return number of failed checks.
label checkGeometry(const scalar maxRatio, const scalar tolerance, autoPtr< coordSetWriter > &setWriter, const scalar minQuality, const bool report) const
All geometric checks. Return number of failed checks.
void writeStats(const List< wordList > &, Ostream &) const
Write some stats.
const List< wordList > & regionNames() const
Region names per surface.
Encapsulates queries for volume refinement ('refine all cells within shell').
Definition: shellSurfaces.H:58
Simple container to keep together snap specific information.
All to do with adding layers.
All to do with snapping to surface.
Identifies a surface patch/zone by name and index, with optional geometric type.
faceZoneType
What to do with faceZone faces.
Implements a timeout mechanism via sigalarm.
Definition: timer.H:84
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A wordRe is a Foam::word, but can contain a regular expression for matching words or strings.
Definition: wordRe.H:83
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const IOdictionary & meshDict
const word dictName("faMeshDefinition")
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
wordList regionNames
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
void inplaceRenumber(const labelUList &oldToNew, IntListType &lists)
Inplace renumber the values (not the indices) of a list of lists.
Definition: ensightOutput.H:90
List< word > names(const UPtrList< T > &list, const UnaryMatchPredicate &matcher)
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
const std::string patch
OpenFOAM patch number as a std::string.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:169
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
label checkGeometry(const polyMesh &mesh, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, autoPtr< coordSetWriter > &setWriter)
prefixOSstream Perr
OSstream wrapped stderr (std::cerr) with parallel prefix.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:349
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:342
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
uint8_t direction
Definition: direction.H:56
IOerror FatalIOError
error FatalError
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.
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
messageStream Warning
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
wordList patchTypes(nPatches)
labelList f(nPoints)
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(rndGen, 1e-3))
IOobject dictIO
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
word setFormat(propsDict.getOrDefault< word >("setFormat", "vtk"))