PDRMesh.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) 2016-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 PDRMesh
29
30Group
31 grpMeshAdvancedUtilities
32
33Description
34 Mesh and field preparation utility for PDR type simulations.
35
36 Reads
37 - cellSet giving blockedCells
38 - faceSets giving blockedFaces and the patch they should go into
39
40 NOTE: To avoid exposing wrong fields values faceSets should include
41 faces contained in the blockedCells cellset.
42
43 - coupledFaces reads coupledFacesSet to introduces mixed-coupled
44 duplicate baffles
45
46 Subsets out the blocked cells and splits the blockedFaces and updates
47 fields.
48
49 The face splitting is done by duplicating the faces. No duplication of
50 points for now (so checkMesh will show a lot of 'duplicate face' messages)
51
52\*---------------------------------------------------------------------------*/
53
54#include "fvMeshSubset.H"
55#include "argList.H"
56#include "cellSet.H"
57#include "BitOps.H"
58#include "IOobjectList.H"
59#include "volFields.H"
60#include "mapPolyMesh.H"
61#include "faceSet.H"
62#include "cellSet.H"
63#include "pointSet.H"
64#include "syncTools.H"
65#include "ReadFields.H"
66#include "polyTopoChange.H"
67#include "polyModifyFace.H"
68#include "polyAddFace.H"
69#include "regionSplit.H"
70#include "Tuple2.H"
71#include "cyclicFvPatch.H"
72
73using namespace Foam;
74
75// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
76
77void modifyOrAddFace
78(
79 polyTopoChange& meshMod,
80 const face& f,
81 const label facei,
82 const label own,
83 const bool flipFaceFlux,
84 const label newPatchi,
85 const label zoneID,
86 const bool zoneFlip,
87
88 bitSet& modifiedFace
89)
90{
91 if (modifiedFace.set(facei))
92 {
93 // First usage of face. Modify.
94 meshMod.setAction
95 (
97 (
98 f, // modified face
99 facei, // label of face
100 own, // owner
101 -1, // neighbour
102 flipFaceFlux, // face flip
103 newPatchi, // patch for face
104 false, // remove from zone
105 zoneID, // zone for face
106 zoneFlip // face flip in zone
107 )
108 );
109 }
110 else
111 {
112 // Second or more usage of face. Add.
113 meshMod.setAction
114 (
116 (
117 f, // modified face
118 own, // owner
119 -1, // neighbour
120 -1, // master point
121 -1, // master edge
122 facei, // master face
123 flipFaceFlux, // face flip
124 newPatchi, // patch for face
125 zoneID, // zone for face
126 zoneFlip // face flip in zone
127 )
128 );
129 }
130}
131
132
133template<class Type>
134void subsetVolFields
135(
136 const fvMeshSubset& subsetter,
137 const IOobjectList& objects,
138 const label patchi,
139 const Type& exposedValue,
141)
142{
144
145 const fvMesh& baseMesh = subsetter.baseMesh();
146
147 label nFields = 0;
148
149 for (const word& fieldName : objects.sortedNames<GeoField>())
150 {
151 const IOobject* ioptr = objects.findObject(fieldName);
152
153 if (!nFields)
154 {
155 Info<< "Subsetting " << GeoField::typeName << nl;
156 }
157 Info<< " " << fieldName << endl;
158
159 GeoField origField(*ioptr, baseMesh);
160
161 subFields.set(nFields, subsetter.interpolate(origField));
162
163 // Explicitly set exposed faces (in patchi) to exposedValue.
164 if (patchi >= 0)
165 {
167 subFields[nFields].boundaryFieldRef()[patchi];
168
169 const label newStart = fld.patch().patch().start();
170 const label oldPatchi = subsetter.patchMap()[patchi];
171
172 if (oldPatchi == -1)
173 {
174 // New patch. Reset whole value.
175 fld = exposedValue;
176 }
177 else
178 {
179 // Reset faces that originate from different patch
180 // or internal faces.
181
182 const fvPatchField<Type>& origPfld =
183 origField.boundaryField()[oldPatchi];
184
185 const label oldSize = origPfld.size();
186 const label oldStart = origPfld.patch().patch().start();
187
188 forAll(fld, j)
189 {
190 const label oldFacei = subsetter.faceMap()[newStart+j];
191
192 if (oldFacei < oldStart || oldFacei >= oldStart+oldSize)
193 {
194 fld[j] = exposedValue;
195 }
196 }
197 }
198
199 ++nFields;
200 }
201 }
202}
203
204
205template<class Type>
206void subsetSurfaceFields
207(
208 const fvMeshSubset& subsetter,
209 const IOobjectList& objects,
210 const label patchi,
211 const Type& exposedValue,
213)
214{
216
217 const fvMesh& baseMesh = subsetter.baseMesh();
218
219 label nFields = 0;
220
221 for (const word& fieldName : objects.sortedNames<GeoField>())
222 {
223 const IOobject* ioptr = objects.findObject(fieldName);
224
225 if (!nFields)
226 {
227 Info<< "Subsetting " << GeoField::typeName << nl;
228 }
229 Info<< " " << fieldName << endl;
230
231 GeoField origField(*ioptr, baseMesh);
232
233 subFields.set(nFields, subsetter.interpolate(origField));
234
235 // Explicitly set exposed faces (in patchi) to exposedValue.
236 if (patchi >= 0)
237 {
239 subFields[nFields].boundaryFieldRef()[patchi];
240
241 const label newStart = fld.patch().patch().start();
242 const label oldPatchi = subsetter.patchMap()[patchi];
243
244 if (oldPatchi == -1)
245 {
246 // New patch. Reset whole value.
247 fld = exposedValue;
248 }
249 else
250 {
251 // Reset faces that originate from different patch
252 // or internal faces.
253
254 const fvsPatchField<Type>& origPfld =
255 origField.boundaryField()[oldPatchi];
256
257 const label oldSize = origPfld.size();
258 const label oldStart = origPfld.patch().patch().start();
259
260 forAll(fld, j)
261 {
262 const label oldFacei = subsetter.faceMap()[newStart+j];
263
264 if (oldFacei < oldStart || oldFacei >= oldStart+oldSize)
265 {
266 fld[j] = exposedValue;
267 }
268 }
269 }
270 }
271
272 ++nFields;
273 }
274}
275
276
277// Faces introduced into zero-sized patches don't get a value at all.
278// This is hack to set them to an initial value.
279template<class GeoField>
280void initCreatedPatches
281(
282 const fvMesh& mesh,
283 const mapPolyMesh& map,
284 const typename GeoField::value_type initValue
285)
286{
288 (
289 mesh.objectRegistry::lookupClass<GeoField>()
290 );
291
292 forAllIters(fields, fieldIter)
293 {
294 GeoField& field = const_cast<GeoField&>(*fieldIter());
295
296 auto& fieldBf = field.boundaryFieldRef();
297
298 forAll(fieldBf, patchi)
299 {
300 if (map.oldPatchSizes()[patchi] == 0)
301 {
302 // Not mapped.
303 fieldBf[patchi] = initValue;
304
305 if (fieldBf[patchi].fixesValue())
306 {
307 fieldBf[patchi] == initValue;
308 }
309 }
310 }
311 }
312}
313
314
315template<class TopoSet>
316void subsetTopoSets
317(
318 const fvMesh& mesh,
319 const IOobjectList& objects,
320 const labelList& map,
321 const fvMesh& subMesh,
322 PtrList<TopoSet>& subSets
323)
324{
325 // Read original sets
326 PtrList<TopoSet> sets;
327 ReadFields<TopoSet>(objects, sets);
328
329 subSets.resize(sets.size());
330 forAll(sets, seti)
331 {
332 TopoSet& set = sets[seti];
333
334 Info<< "Subsetting " << set.type() << " " << set.name() << endl;
335
336 // Map the data
337 bitSet isSet(set.maxSize(mesh));
338 for (const label id : set)
339 {
340 isSet.set(id);
341 }
342
343 label nSet = 0;
344 for (const label id : map)
345 {
346 if (isSet.test(id))
347 {
348 ++nSet;
349 }
350 }
351
352 subSets.set
353 (
354 seti,
355 new TopoSet(subMesh, set.name(), nSet, IOobject::AUTO_WRITE)
356 );
357 TopoSet& subSet = subSets[seti];
358
359 forAll(map, i)
360 {
361 if (isSet.test(map[i]))
362 {
363 subSet.insert(i);
364 }
365 }
366 }
367}
368
369
370void createCoupledBaffles
371(
372 fvMesh& mesh,
373 const labelList& coupledWantedPatch,
374 polyTopoChange& meshMod,
375 bitSet& modifiedFace
376)
377{
378 const faceZoneMesh& faceZones = mesh.faceZones();
379
380 forAll(coupledWantedPatch, facei)
381 {
382 if (coupledWantedPatch[facei] != -1)
383 {
384 const face& f = mesh.faces()[facei];
385 label zoneID = faceZones.whichZone(facei);
386 bool zoneFlip = false;
387
388 if (zoneID >= 0)
389 {
390 const faceZone& fZone = faceZones[zoneID];
391 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
392 }
393
394 // Use owner side of face
395 modifyOrAddFace
396 (
397 meshMod,
398 f, // modified face
399 facei, // label of face
400 mesh.faceOwner()[facei], // owner
401 false, // face flip
402 coupledWantedPatch[facei], // patch for face
403 zoneID, // zone for face
404 zoneFlip, // face flip in zone
405 modifiedFace // modify or add status
406 );
407
408 if (mesh.isInternalFace(facei))
409 {
410 label zoneID = faceZones.whichZone(facei);
411 bool zoneFlip = false;
412
413 if (zoneID >= 0)
414 {
415 const faceZone& fZone = faceZones[zoneID];
416 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
417 }
418 // Use neighbour side of face
419 modifyOrAddFace
420 (
421 meshMod,
422 f.reverseFace(), // modified face
423 facei, // label of face
424 mesh.faceNeighbour()[facei],// owner
425 false, // face flip
426 coupledWantedPatch[facei], // patch for face
427 zoneID, // zone for face
428 zoneFlip, // face flip in zone
429 modifiedFace // modify or add status
430 );
431 }
432 }
433 }
434}
435
436
437void createCyclicCoupledBaffles
438(
439 fvMesh& mesh,
440 const labelList& cyclicMasterPatch,
441 const labelList& cyclicSlavePatch,
442 polyTopoChange& meshMod,
443 bitSet& modifiedFace
444)
445{
446 const faceZoneMesh& faceZones = mesh.faceZones();
447
448 forAll(cyclicMasterPatch, facei)
449 {
450 if (cyclicMasterPatch[facei] != -1)
451 {
452 const face& f = mesh.faces()[facei];
453
454 label zoneID = faceZones.whichZone(facei);
455 bool zoneFlip = false;
456
457 if (zoneID >= 0)
458 {
459 const faceZone& fZone = faceZones[zoneID];
460 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
461 }
462
463 modifyOrAddFace
464 (
465 meshMod,
466 f.reverseFace(), // modified face
467 facei, // label of face
468 mesh.faceNeighbour()[facei], // owner
469 false, // face flip
470 cyclicMasterPatch[facei], // patch for face
471 zoneID, // zone for face
472 zoneFlip, // face flip in zone
473 modifiedFace // modify or add
474 );
475 }
476 }
477
478 forAll(cyclicSlavePatch, facei)
479 {
480 if (cyclicSlavePatch[facei] != -1)
481 {
482 const face& f = mesh.faces()[facei];
483 if (mesh.isInternalFace(facei))
484 {
485 label zoneID = faceZones.whichZone(facei);
486 bool zoneFlip = false;
487
488 if (zoneID >= 0)
489 {
490 const faceZone& fZone = faceZones[zoneID];
491 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
492 }
493 // Use owner side of face
494 modifyOrAddFace
495 (
496 meshMod,
497 f, // modified face
498 facei, // label of face
499 mesh.faceOwner()[facei], // owner
500 false, // face flip
501 cyclicSlavePatch[facei], // patch for face
502 zoneID, // zone for face
503 zoneFlip, // face flip in zone
504 modifiedFace // modify or add status
505 );
506 }
507 }
508 }
509}
510
511
512void createBaffles
513(
514 fvMesh& mesh,
515 const labelList& wantedPatch,
516 polyTopoChange& meshMod
517)
518{
519 const faceZoneMesh& faceZones = mesh.faceZones();
520 Info << "faceZone:createBaffle " << faceZones << endl;
521 forAll(wantedPatch, facei)
522 {
523 if (wantedPatch[facei] != -1)
524 {
525 const face& f = mesh.faces()[facei];
526
527 label zoneID = faceZones.whichZone(facei);
528 bool zoneFlip = false;
529
530 if (zoneID >= 0)
531 {
532 const faceZone& fZone = faceZones[zoneID];
533 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
534 }
535
536 meshMod.setAction
537 (
539 (
540 f, // modified face
541 facei, // label of face
542 mesh.faceOwner()[facei], // owner
543 -1, // neighbour
544 false, // face flip
545 wantedPatch[facei], // patch for face
546 false, // remove from zone
547 zoneID, // zone for face
548 zoneFlip // face flip in zone
549 )
550 );
551
552 if (mesh.isInternalFace(facei))
553 {
554 label zoneID = faceZones.whichZone(facei);
555 bool zoneFlip = false;
556
557 if (zoneID >= 0)
558 {
559 const faceZone& fZone = faceZones[zoneID];
560 zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
561 }
562
563 meshMod.setAction
564 (
566 (
567 f.reverseFace(), // modified face
568 mesh.faceNeighbour()[facei],// owner
569 -1, // neighbour
570 -1, // masterPointID
571 -1, // masterEdgeID
572 facei, // masterFaceID,
573 false, // face flip
574 wantedPatch[facei], // patch for face
575 zoneID, // zone for face
576 zoneFlip // face flip in zone
577 )
578 );
579 }
580 }
581 }
582}
583
584
585// Wrapper around find patch. Also makes sure same patch in parallel.
586label findPatch(const polyBoundaryMesh& patches, const word& patchName)
587{
588 label patchi = patches.findPatchID(patchName);
589
590 if (patchi == -1)
591 {
593 << "Illegal patch " << patchName
594 << nl << "Valid patches are " << patches.names()
595 << exit(FatalError);
596 }
597
598 // Check same patch for all procs
599 {
600 label newPatch = patchi;
601 reduce(newPatch, minOp<label>());
602
603 if (newPatch != patchi)
604 {
606 << "Patch " << patchName
607 << " should have the same patch index on all processors." << nl
608 << "On my processor it has index " << patchi
609 << " ; on some other processor it has index " << newPatch
610 << exit(FatalError);
611 }
612 }
613 return patchi;
614}
615
616
617
618int main(int argc, char *argv[])
619{
620 argList::addNote
621 (
622 "Mesh and field preparation utility for PDR type simulations."
623 );
624 #include "addOverwriteOption.H"
625
626 argList::noFunctionObjects(); // Never use function objects
627
628 #include "setRootCase.H"
629 #include "createTime.H"
630 #include "createNamedMesh.H"
631
632 // Read control dictionary
633 // ~~~~~~~~~~~~~~~~~~~~~~~
634
636 (
638 (
639 "PDRMeshDict",
640 runTime.system(),
641 mesh,
642 IOobject::MUST_READ_IF_MODIFIED,
643 IOobject::NO_WRITE
644 )
645 );
646
647 // Per faceSet the patch to put the baffles into
648 const List<Pair<word>> setsAndPatches(dict.lookup("blockedFaces"));
649
650 // Per faceSet the patch to put the coupled baffles into
651 DynamicList<FixedList<word, 3>> coupledAndPatches(10);
652
653 const dictionary& functionDicts = dict.subDict("coupledFaces");
654
655 for (const entry& dEntry : functionDicts)
656 {
657 if (!dEntry.isDict()) // Safety
658 {
659 continue;
660 }
661
662 const word& key = dEntry.keyword();
663 const dictionary& dict = dEntry.dict();
664
665 const word cyclicName = dict.get<word>("cyclicMasterPatch");
666 const word wallName = dict.get<word>("wallPatch");
667 FixedList<word, 3> nameAndType;
668 nameAndType[0] = key;
669 nameAndType[1] = wallName;
670 nameAndType[2] = cyclicName;
671 coupledAndPatches.append(nameAndType);
672 }
673
674 forAll(setsAndPatches, setI)
675 {
676 Info<< "Faces in faceSet " << setsAndPatches[setI][0]
677 << " become baffles in patch " << setsAndPatches[setI][1]
678 << endl;
679 }
680
681 forAll(coupledAndPatches, setI)
682 {
683 Info<< "Faces in faceSet " << coupledAndPatches[setI][0]
684 << " become coupled baffles in patch " << coupledAndPatches[setI][1]
685 << endl;
686 }
687
688 // All exposed faces that are not explicitly marked to be put into a patch
689 const word defaultPatch(dict.get<word>("defaultPatch"));
690
691 Info<< "Faces that get exposed become boundary faces in patch "
692 << defaultPatch << endl;
693
694 const word blockedSetName(dict.get<word>("blockedCells"));
695
696 Info<< "Reading blocked cells from cellSet " << blockedSetName
697 << endl;
698
699 const bool overwrite = args.found("overwrite");
700
701
702 // Read faceSets, lookup patches
703 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
704
705 // Check that face sets don't have coincident faces
706 labelList wantedPatch(mesh.nFaces(), -1);
707 forAll(setsAndPatches, setI)
708 {
709 faceSet fSet(mesh, setsAndPatches[setI][0]);
710
711 label patchi = findPatch
712 (
713 mesh.boundaryMesh(),
714 setsAndPatches[setI][1]
715 );
716
717 for (const label facei : fSet)
718 {
719 if (wantedPatch[facei] != -1)
720 {
722 << "Face " << facei
723 << " is in faceSet " << setsAndPatches[setI][0]
724 << " destined for patch " << setsAndPatches[setI][1]
725 << " but also in patch " << wantedPatch[facei]
726 << exit(FatalError);
727 }
728 wantedPatch[facei] = patchi;
729 }
730 }
731
732 // Per face the patch for coupled baffle or -1.
733 labelList coupledWantedPatch(mesh.nFaces(), -1);
734 labelList cyclicWantedPatch_half0(mesh.nFaces(), -1);
735 labelList cyclicWantedPatch_half1(mesh.nFaces(), -1);
736
737 forAll(coupledAndPatches, setI)
738 {
739 const polyBoundaryMesh& patches = mesh.boundaryMesh();
740 const label cyclicId =
741 findPatch(patches, coupledAndPatches[setI][2]);
742
743 const label cyclicSlaveId = findPatch
744 (
745 patches,
746 refCast<const cyclicFvPatch>
747 (
748 mesh.boundary()[cyclicId]
749 ).neighbFvPatch().name()
750 );
751
752 faceSet fSet(mesh, coupledAndPatches[setI][0]);
753 label patchi = findPatch(patches, coupledAndPatches[setI][1]);
754
755 for (const label facei : fSet)
756 {
757 if (coupledWantedPatch[facei] != -1)
758 {
760 << "Face " << facei
761 << " is in faceSet " << coupledAndPatches[setI][0]
762 << " destined for patch " << coupledAndPatches[setI][1]
763 << " but also in patch " << coupledWantedPatch[facei]
764 << exit(FatalError);
765 }
766
767 coupledWantedPatch[facei] = patchi;
768 cyclicWantedPatch_half0[facei] = cyclicId;
769 cyclicWantedPatch_half1[facei] = cyclicSlaveId;
770 }
771 }
772
773 // Exposed faces patch
774 label defaultPatchi = findPatch(mesh.boundaryMesh(), defaultPatch);
775
776
777 //
778 // Removing blockedCells (like subsetMesh)
779 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
780 //
781
782 // Mesh subsetting engine
783 fvMeshSubset subsetter
784 (
785 mesh,
786 BitSetOps::create
787 (
788 mesh.nCells(),
789 cellSet(mesh, blockedSetName), // Blocked cells as labelHashSet
790 false // on=false: invert logic => retain the unblocked cells
791 ),
792 defaultPatchi,
793 true
794 );
795
796
797 // Subset wantedPatch. Note that might also include boundary faces
798 // that have been exposed by subsetting.
799 wantedPatch = IndirectList<label>(wantedPatch, subsetter.faceMap())();
800
801 coupledWantedPatch = IndirectList<label>
802 (
803 coupledWantedPatch,
804 subsetter.faceMap()
805 )();
806
807 cyclicWantedPatch_half0 = IndirectList<label>
808 (
809 cyclicWantedPatch_half0,
810 subsetter.faceMap()
811 )();
812
813 cyclicWantedPatch_half1 = IndirectList<label>
814 (
815 cyclicWantedPatch_half1,
816 subsetter.faceMap()
817 )();
818
819 // Read all fields in time and constant directories
820 IOobjectList objects(mesh, runTime.timeName());
821 {
822 IOobjectList timeObjects(mesh, mesh.facesInstance());
823
824 // Transfer specific types
825 forAllIters(timeObjects, iter)
826 {
827 autoPtr<IOobject> objPtr(timeObjects.remove(iter));
828 const auto& obj = *objPtr;
829
830 if
831 (
832 obj.isHeaderClass<volScalarField>()
842 )
843 {
844 objects.add(objPtr);
845 }
846 }
847 }
848 // Read vol fields and subset.
849
850 wordList scalarNames(objects.sortedNames<volScalarField>());
851 PtrList<volScalarField> scalarFlds(scalarNames.size());
852 subsetVolFields
853 (
854 subsetter,
855 objects,
856 defaultPatchi,
857 scalar(Zero),
858 scalarFlds
859 );
860
861 wordList vectorNames(objects.sortedNames<volVectorField>());
862 PtrList<volVectorField> vectorFlds(vectorNames.size());
863 subsetVolFields
864 (
865 subsetter,
866 objects,
867 defaultPatchi,
868 vector(Zero),
869 vectorFlds
870 );
871
872 wordList sphTensorNames
873 (
875 );
877 (
878 sphTensorNames.size()
879 );
880 subsetVolFields
881 (
882 subsetter,
883 objects,
884 defaultPatchi,
885 sphericalTensor(Zero),
886 sphTensorFlds
887 );
888
889 wordList symmTensorNames(objects.sortedNames<volSymmTensorField>());
890 PtrList<volSymmTensorField> symmTensorFlds(symmTensorNames.size());
891 subsetVolFields
892 (
893 subsetter,
894 objects,
895 defaultPatchi,
896 symmTensor(Zero),
897 symmTensorFlds
898 );
899
900 wordList tensorNames(objects.sortedNames<volTensorField>());
901 PtrList<volTensorField> tensorFlds(tensorNames.size());
902 subsetVolFields
903 (
904 subsetter,
905 objects,
906 defaultPatchi,
907 tensor(Zero),
908 tensorFlds
909 );
910
911 // Read surface fields and subset.
912
913 wordList surfScalarNames(objects.sortedNames<surfaceScalarField>());
914 PtrList<surfaceScalarField> surfScalarFlds(surfScalarNames.size());
915 subsetSurfaceFields
916 (
917 subsetter,
918 objects,
919 defaultPatchi,
920 scalar(Zero),
921 surfScalarFlds
922 );
923
924 wordList surfVectorNames(objects.sortedNames<surfaceVectorField>());
925 PtrList<surfaceVectorField> surfVectorFlds(surfVectorNames.size());
926 subsetSurfaceFields
927 (
928 subsetter,
929 objects,
930 defaultPatchi,
931 vector(Zero),
932 surfVectorFlds
933 );
934
935 wordList surfSphTensorNames
936 (
938 );
939 PtrList<surfaceSphericalTensorField> surfSphericalTensorFlds
940 (
941 surfSphTensorNames.size()
942 );
943 subsetSurfaceFields
944 (
945 subsetter,
946 objects,
947 defaultPatchi,
948 sphericalTensor(Zero),
949 surfSphericalTensorFlds
950 );
951
952 wordList surfSymmTensorNames
953 (
955 );
956
957 PtrList<surfaceSymmTensorField> surfSymmTensorFlds
958 (
959 surfSymmTensorNames.size()
960 );
961
962 subsetSurfaceFields
963 (
964 subsetter,
965 objects,
966 defaultPatchi,
967 symmTensor(Zero),
968 surfSymmTensorFlds
969 );
970
971 wordList surfTensorNames(objects.sortedNames<surfaceTensorField>());
972 PtrList<surfaceTensorField> surfTensorFlds(surfTensorNames.size());
973 subsetSurfaceFields
974 (
975 subsetter,
976 objects,
977 defaultPatchi,
978 tensor(Zero),
979 surfTensorFlds
980 );
981
982
983 // Set handling
984 PtrList<cellSet> cellSets;
985 PtrList<faceSet> faceSets;
986 PtrList<pointSet> pointSets;
987 {
988 IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
989 subsetTopoSets
990 (
991 mesh,
992 objects,
993 subsetter.cellMap(),
994 subsetter.subMesh(),
995 cellSets
996 );
997 subsetTopoSets
998 (
999 mesh,
1000 objects,
1001 subsetter.faceMap(),
1002 subsetter.subMesh(),
1003 faceSets
1004 );
1005 subsetTopoSets
1006 (
1007 mesh,
1008 objects,
1009 subsetter.pointMap(),
1010 subsetter.subMesh(),
1011 pointSets
1012 );
1013 }
1014
1015
1016 if (!overwrite)
1017 {
1018 ++runTime;
1019 }
1020
1021 Info<< "Writing mesh without blockedCells to time " << runTime.value()
1022 << endl;
1023
1024 // Subsetting adds 'subset' prefix. Rename field to be like original.
1025 forAll(scalarFlds, i)
1026 {
1027 scalarFlds[i].rename(scalarNames[i]);
1028 scalarFlds[i].writeOpt(IOobject::AUTO_WRITE);
1029 }
1030 forAll(vectorFlds, i)
1031 {
1032 vectorFlds[i].rename(vectorNames[i]);
1033 vectorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1034 }
1035 forAll(sphTensorFlds, i)
1036 {
1037 sphTensorFlds[i].rename(sphTensorNames[i]);
1038 sphTensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1039 }
1040 forAll(symmTensorFlds, i)
1041 {
1042 symmTensorFlds[i].rename(symmTensorNames[i]);
1043 symmTensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1044 }
1045 forAll(tensorFlds, i)
1046 {
1047 tensorFlds[i].rename(tensorNames[i]);
1048 tensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1049 }
1050
1051 // Surface ones.
1052 forAll(surfScalarFlds, i)
1053 {
1054 surfScalarFlds[i].rename(surfScalarNames[i]);
1055 surfScalarFlds[i].writeOpt(IOobject::AUTO_WRITE);
1056 }
1057 forAll(surfVectorFlds, i)
1058 {
1059 surfVectorFlds[i].rename(surfVectorNames[i]);
1060 surfVectorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1061 }
1062 forAll(surfSphericalTensorFlds, i)
1063 {
1064 surfSphericalTensorFlds[i].rename(surfSphTensorNames[i]);
1065 surfSphericalTensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1066 }
1067 forAll(surfSymmTensorFlds, i)
1068 {
1069 surfSymmTensorFlds[i].rename(surfSymmTensorNames[i]);
1070 surfSymmTensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1071 }
1072 forAll(surfTensorNames, i)
1073 {
1074 surfTensorFlds[i].rename(surfTensorNames[i]);
1075 surfTensorFlds[i].writeOpt(IOobject::AUTO_WRITE);
1076 }
1077
1078 subsetter.subMesh().write();
1079
1080
1081 //
1082 // Splitting blockedFaces
1083 // ~~~~~~~~~~~~~~~~~~~~~~
1084 //
1085
1086 // Synchronize wantedPatch across coupled patches.
1087 syncTools::syncFaceList
1088 (
1089 subsetter.subMesh(),
1090 wantedPatch,
1092 );
1093
1094 // Synchronize coupledWantedPatch across coupled patches.
1095 syncTools::syncFaceList
1096 (
1097 subsetter.subMesh(),
1098 coupledWantedPatch,
1100 );
1101
1102 // Synchronize cyclicWantedPatch across coupled patches.
1103 syncTools::syncFaceList
1104 (
1105 subsetter.subMesh(),
1106 cyclicWantedPatch_half0,
1108 );
1109
1110 // Synchronize cyclicWantedPatch across coupled patches.
1111 syncTools::syncFaceList
1112 (
1113 subsetter.subMesh(),
1114 cyclicWantedPatch_half1,
1116 );
1117
1118 // Topochange container
1119 polyTopoChange meshMod(subsetter.subMesh());
1120
1121
1122 // Whether first use of face (modify) or consecutive (add)
1123 bitSet modifiedFace(mesh.nFaces());
1124
1125 // Create coupled wall-side baffles
1126 createCoupledBaffles
1127 (
1128 subsetter.subMesh(),
1129 coupledWantedPatch,
1130 meshMod,
1131 modifiedFace
1132 );
1133
1134 // Create coupled master/slave cyclic baffles
1135 createCyclicCoupledBaffles
1136 (
1137 subsetter.subMesh(),
1138 cyclicWantedPatch_half0,
1139 cyclicWantedPatch_half1,
1140 meshMod,
1141 modifiedFace
1142 );
1143
1144 // Create wall baffles
1145 createBaffles
1146 (
1147 subsetter.subMesh(),
1148 wantedPatch,
1149 meshMod
1150 );
1151
1152 if (!overwrite)
1153 {
1154 ++runTime;
1155 }
1156
1157 // Change the mesh. Change points directly (no inflation).
1158 autoPtr<mapPolyMesh> mapPtr =
1159 meshMod.changeMesh(subsetter.subMesh(), false);
1160 mapPolyMesh& map = *mapPtr;
1161
1162 // Update fields
1163 subsetter.subMesh().updateMesh(map);
1164
1165 // Fix faces that get mapped to zero-sized patches (these don't get any
1166 // value)
1167 initCreatedPatches<volScalarField>
1168 (
1169 subsetter.subMesh(),
1170 map,
1171 Zero
1172 );
1173 initCreatedPatches<volVectorField>
1174 (
1175 subsetter.subMesh(),
1176 map,
1177 Zero
1178 );
1179 initCreatedPatches<volSphericalTensorField>
1180 (
1181 subsetter.subMesh(),
1182 map,
1183 Zero
1184 );
1185 initCreatedPatches<volSymmTensorField>
1186 (
1187 subsetter.subMesh(),
1188 map,
1189 Zero
1190 );
1191 initCreatedPatches<volTensorField>
1192 (
1193 subsetter.subMesh(),
1194 map,
1195 Zero
1196 );
1197
1198 initCreatedPatches<surfaceScalarField>
1199 (
1200 subsetter.subMesh(),
1201 map,
1202 Zero
1203 );
1204 initCreatedPatches<surfaceVectorField>
1205 (
1206 subsetter.subMesh(),
1207 map,
1208 Zero
1209 );
1210 initCreatedPatches<surfaceSphericalTensorField>
1211 (
1212 subsetter.subMesh(),
1213 map,
1214 Zero
1215 );
1216 initCreatedPatches<surfaceSymmTensorField>
1217 (
1218 subsetter.subMesh(),
1219 map,
1220 Zero
1221 );
1222 initCreatedPatches<surfaceTensorField>
1223 (
1224 subsetter.subMesh(),
1225 map,
1226 Zero
1227 );
1228
1229 // Update numbering of topoSets
1230 topoSet::updateMesh(subsetter.subMesh().facesInstance(), map, cellSets);
1231 topoSet::updateMesh(subsetter.subMesh().facesInstance(), map, faceSets);
1232 topoSet::updateMesh(subsetter.subMesh().facesInstance(), map, pointSets);
1233
1234
1235 // Move mesh (since morphing might not do this)
1236 if (map.hasMotionPoints())
1237 {
1238 subsetter.subMesh().movePoints(map.preMotionPoints());
1239 }
1240
1241 Info<< "Writing mesh with split blockedFaces to time " << runTime.value()
1242 << endl;
1243
1244 subsetter.subMesh().write();
1245
1246
1247 //
1248 // Removing inaccessible regions
1249 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1250 //
1251
1252 // Determine connected regions. regionSplit is the labelList with the
1253 // region per cell.
1254 regionSplit cellRegion(subsetter.subMesh());
1255
1256 if (cellRegion.nRegions() > 1)
1257 {
1259 << "Removing blocked faces and cells created "
1260 << cellRegion.nRegions()
1261 << " regions that are not connected via a face." << nl
1262 << " This is not supported in solvers." << nl
1263 << " Use" << nl << nl
1264 << " splitMeshRegions <root> <case> -largestOnly" << nl << nl
1265 << " to extract a single region of the mesh." << nl
1266 << " This mesh will be written to a new timedirectory"
1267 << " so might have to be moved back to constant/" << nl
1268 << endl;
1269
1270 const word startFrom(runTime.controlDict().get<word>("startFrom"));
1271
1272 if (startFrom != "latestTime")
1273 {
1275 << "To run splitMeshRegions please set your"
1276 << " startFrom entry to latestTime" << endl;
1277 }
1278 }
1279
1280 Info<< "\nEnd\n" << endl;
1281
1282 return 0;
1283}
1284
1285
1286// ************************************************************************* //
Field reading functions for post-processing utilities.
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
reduce(hasMovingMesh, orOp< bool >())
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
Generic GeometricField class.
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
const IOobject * findObject(const word &objName) const
Return const pointer to the object found by name.
Definition: IOobjectList.C:293
wordList sortedNames() const
The sorted names of the IOobjects.
Definition: IOobjectList.C:383
bool add(autoPtr< IOobject > &objectPtr)
Add IOobject to the list.
Definition: IOobjectList.C:187
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
bool isHeaderClass() const
Check if headerClassName() equals Type::typeName.
Definition: IOobjectI.H:156
A List with indirect addressing.
Definition: IndirectList.H:119
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
void resize(const label newLen)
Adjust size of PtrList.
Definition: PtrList.C:103
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
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:288
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
A collection of cell labels.
Definition: cellSet.H:54
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:70
A list of face labels.
Definition: faceSet.H:54
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:372
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Holds a reference to the original mesh (the baseMesh) and optionally to a subset of that mesh (the su...
Definition: fvMeshSubset.H:80
const fvMesh & baseMesh() const noexcept
Original mesh.
Definition: fvMeshSubsetI.H:30
const labelList & faceMap() const
Return face map.
Definition: fvMeshSubsetI.H:72
const labelList & cellMap() const
Return cell map.
Definition: fvMeshSubsetI.H:91
static tmp< DimensionedField< Type, volMesh > > interpolate(const DimensionedField< Type, volMesh > &, const fvMesh &sMesh, const labelUList &cellMap)
Map volume internal (dimensioned) field.
const fvMesh & subMesh() const
Return reference to subset mesh.
Definition: fvMeshSubsetI.H:48
const labelList & patchMap() const
Return patch map.
Definition: fvMeshSubsetI.H:99
const labelList & pointMap() const
Return point map.
Definition: fvMeshSubsetI.H:64
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
virtual bool write(const bool valid=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1079
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
Definition: fvMesh.C:971
virtual void movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:895
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:82
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:362
const polyPatch & patch() const
Return the polyPatch.
Definition: fvPatch.H:167
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:79
const fvPatch & patch() const
Return patch.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:613
const labelList & oldPatchSizes() const
Return list of the old patch sizes.
Definition: mapPolyMesh.H:626
bool hasMotionPoints() const
Has valid preMotionPoints?
Definition: mapPolyMesh.H:619
A face addition data class. A face can be inflated either from a point or from another face and can e...
Definition: polyAddFace.H:56
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
Class describing modification of a face.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Direct mesh changes based on v1.3 polyTopoChange syntax.
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:144
Tensor of scalars, i.e. Tensor<scalar>.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
rDeltaTY field()
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelIOList & zoneID
#define WarningInFunction
Report a warning using Foam::Warning.
void set(List< bool > &bools, const labelUList &locations)
Set the listed locations (assign 'true').
Definition: BitOps.C:38
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
Foam::argList args(argc, argv)
multivariateSurfaceInterpolationScheme< scalar >::fieldTable fields
Definition: createFields.H:97
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:260
static const char *const typeName
The type name used in ensight case files.