splitMeshRegions.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-2017 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 splitMeshRegions
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 Splits mesh into multiple regions.
35
36 Each region is defined as a domain whose cells can all be reached by
37 cell-face-cell walking without crossing
38 - boundary faces
39 - additional faces from faceset (-blockedFaces faceSet).
40 - any face between differing cellZones (-cellZones)
41
42 Output is:
43 - volScalarField with regions as different scalars (-detectOnly)
44 or
45 - mesh with multiple regions and mapped patches. These patches
46 either cover the whole interface between two region (default) or
47 only part according to faceZones (-useFaceZones)
48 or
49 - mesh with cells put into cellZones (-makeCellZones)
50
51 Note:
52
53 - multiple cellZones can be combined into a single region (cluster)
54 for further analysis using the 'addZones' or 'combineZones' option:
55 -addZones '((allSolids zoneA "zoneB.*")(allFluids none otherZone))'
56 or
57 -combineZones '((zoneA "zoneB.*")(none otherZone))
58 This can be combined with e.g. 'cellZones' or 'cellZonesOnly'. The
59 addZones option supplies the destination region name as first element in
60 the list. The combineZones option synthesises the region name e.g.
61 zoneA_zoneB0_zoneB1
62
63 - cellZonesOnly does not do a walk and uses the cellZones only. Use
64 this if you don't mind having disconnected domains in a single region.
65 This option requires all cells to be in one (and one only) cellZone.
66
67 - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
68 from the specified file. This allows one to explicitly specify the region
69 distribution and still have multiple cellZones per region.
70
71 - prefixRegion prefixes all normal patches with region name (interface
72 (patches already have region name prefix)
73
74 - Should work in parallel.
75 cellZones can differ on either side of processor boundaries in which case
76 the faces get moved from processor patch to mapped patch. Not very well
77 tested.
78
79 - If a cell zone gets split into more than one region it can detect
80 the largest matching region (-sloppyCellZones). This will accept any
81 region that covers more than 50% of the zone. It has to be a subset
82 so cannot have any cells in any other zone.
83
84 - If explicitly a single region has been selected (-largestOnly or
85 -insidePoint) its region name will be either
86 - name of a cellZone it matches to or
87 - "largestOnly" respectively "insidePoint" or
88 - polyMesh::defaultRegion if additionally -overwrite
89 (so it will overwrite the input mesh!)
90
91 - writes maps like decomposePar back to original mesh:
92 - pointRegionAddressing : for every point in this region the point in
93 the original mesh
94 - cellRegionAddressing : ,, cell ,, cell ,,
95 - faceRegionAddressing : ,, face ,, face in
96 the original mesh + 'turning index'. For a face in the same orientation
97 this is the original facelabel+1, for a turned face this is -facelabel-1
98 - boundaryRegionAddressing : for every patch in this region the
99 patch in the original mesh (or -1 if added patch)
100
101\*---------------------------------------------------------------------------*/
102
103#include "argList.H"
104#include "regionSplit.H"
105#include "fvMeshSubset.H"
106#include "IOobjectList.H"
107#include "volFields.H"
108#include "faceSet.H"
109#include "cellSet.H"
110#include "polyTopoChange.H"
111#include "removeCells.H"
112#include "edgeHashes.H"
113#include "syncTools.H"
114#include "ReadFields.H"
115#include "mappedWallPolyPatch.H"
116#include "fvMeshTools.H"
118#include "processorMeshes.H"
119
120using namespace Foam;
121
122// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
123
124// Prepend prefix to selected patches.
125void renamePatches
126(
127 fvMesh& mesh,
128 const word& prefix,
129 const labelList& patchesToRename
130)
131{
132 polyBoundaryMesh& polyPatches =
133 const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
134 forAll(patchesToRename, i)
135 {
136 label patchi = patchesToRename[i];
137 polyPatch& pp = polyPatches[patchi];
138
139 if (isA<coupledPolyPatch>(pp))
140 {
142 << "Encountered coupled patch " << pp.name()
143 << ". Will only rename the patch itself,"
144 << " not any referred patches."
145 << " This might have to be done by hand."
146 << endl;
147 }
148
149 pp.name() = prefix + '_' + pp.name();
150 }
151 // Recalculate any demand driven data (e.g. group to name lookup)
152 polyPatches.updateMesh();
153}
154
155
156template<class GeoField>
157void subsetVolFields
158(
159 const fvMesh& mesh,
160 const fvMesh& subMesh,
161 const labelList& cellMap,
162 const labelList& faceMap,
163 const labelHashSet& addedPatches
164)
165{
166 const labelList patchMap(identity(mesh.boundaryMesh().size()));
167
169 (
170 mesh.objectRegistry::lookupClass<GeoField>()
171 );
173 {
174 const GeoField& fld = *iter.val();
175
176 Info<< "Mapping field " << fld.name() << endl;
177
178 tmp<GeoField> tSubFld
179 (
180 fvMeshSubset::interpolate
181 (
182 fld,
183 subMesh,
184 patchMap,
185 cellMap,
186 faceMap
187 )
188 );
189
190 // Hack: set value to 0 for introduced patches (since don't
191 // get initialised.
192 forAll(tSubFld().boundaryField(), patchi)
193 {
194 if (addedPatches.found(patchi))
195 {
196 tSubFld.ref().boundaryFieldRef()[patchi] ==
197 typename GeoField::value_type(Zero);
198 }
199 }
200
201 // Store on subMesh
202 GeoField* subFld = tSubFld.ptr();
203 subFld->rename(fld.name());
204 subFld->writeOpt(IOobject::AUTO_WRITE);
205 subFld->store();
206 }
207}
208
209
210template<class GeoField>
211void subsetSurfaceFields
212(
213 const fvMesh& mesh,
214 const fvMesh& subMesh,
215 const labelList& cellMap,
216 const labelList& faceMap,
217 const labelHashSet& addedPatches
218)
219{
220 const labelList patchMap(identity(mesh.boundaryMesh().size()));
221
223 (
224 mesh.objectRegistry::lookupClass<GeoField>()
225 );
227 {
228 const GeoField& fld = *iter.val();
229
230 Info<< "Mapping field " << fld.name() << endl;
231
232 tmp<GeoField> tSubFld
233 (
234 fvMeshSubset::interpolate
235 (
236 fld,
237 subMesh,
238 patchMap,
239 cellMap,
240 faceMap
241 )
242 );
243
244 // Hack: set value to 0 for introduced patches (since don't
245 // get initialised.
246 forAll(tSubFld().boundaryField(), patchi)
247 {
248 if (addedPatches.found(patchi))
249 {
250 tSubFld.ref().boundaryFieldRef()[patchi] ==
251 typename GeoField::value_type(Zero);
252 }
253 }
254
255 // Store on subMesh
256 GeoField* subFld = tSubFld.ptr();
257 subFld->rename(fld.name());
258 subFld->writeOpt(IOobject::AUTO_WRITE);
259 subFld->store();
260 }
261}
262
263// Select all cells not in the region
264labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
265{
266 DynamicList<label> nonRegionCells(cellRegion.size());
267 forAll(cellRegion, celli)
268 {
269 if (cellRegion[celli] != regionI)
270 {
271 nonRegionCells.append(celli);
272 }
273 }
274 return nonRegionCells.shrink();
275}
276
277
278void addToInterface
279(
280 const polyMesh& mesh,
281 const label zoneID,
282 const label ownRegion,
283 const label neiRegion,
284 EdgeMap<Map<label>>& regionsToSize
285)
286{
288 (
289 min(ownRegion, neiRegion),
290 max(ownRegion, neiRegion)
291 );
292
293 auto iter = regionsToSize.find(interface);
294
295 if (iter.found())
296 {
297 // Check if zone present
298 auto zoneIter = iter().find(zoneID);
299 if (zoneIter.found())
300 {
301 // Found zone. Increment count.
302 ++(*zoneIter);
303 }
304 else
305 {
306 // New or no zone. Insert with count 1.
307 iter().insert(zoneID, 1);
308 }
309 }
310 else
311 {
312 // Create new interface of size 1.
313 Map<label> zoneToSize;
314 zoneToSize.insert(zoneID, 1);
315 regionsToSize.insert(interface, zoneToSize);
316 }
317}
318
319
320// Get region-region interface name and sizes.
321// Returns interfaces as straight list for looping in identical order.
322void getInterfaceSizes
323(
324 const polyMesh& mesh,
325 const bool useFaceZones,
326 const labelList& cellRegion,
327 const wordList& regionNames,
328
329 edgeList& interfaces,
330 List<Pair<word>>& interfaceNames,
331 labelList& interfaceSizes,
332 labelList& faceToInterface
333)
334{
335 // From region-region to faceZone (or -1) to number of faces.
336
337 EdgeMap<Map<label>> regionsToSize;
338
339
340 // Internal faces
341 // ~~~~~~~~~~~~~~
342
343 forAll(mesh.faceNeighbour(), facei)
344 {
345 label ownRegion = cellRegion[mesh.faceOwner()[facei]];
346 label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
347
348 if (ownRegion != neiRegion)
349 {
350 addToInterface
351 (
352 mesh,
353 (useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
354 ownRegion,
355 neiRegion,
356 regionsToSize
357 );
358 }
359 }
360
361 // Boundary faces
362 // ~~~~~~~~~~~~~~
363
364 // Neighbour cellRegion.
365 labelList coupledRegion(mesh.nBoundaryFaces());
366
367 forAll(coupledRegion, i)
368 {
369 label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
370 coupledRegion[i] = cellRegion[celli];
371 }
372 syncTools::swapBoundaryFaceList(mesh, coupledRegion);
373
374 forAll(coupledRegion, i)
375 {
376 label facei = i+mesh.nInternalFaces();
377 label ownRegion = cellRegion[mesh.faceOwner()[facei]];
378 label neiRegion = coupledRegion[i];
379
380 if (ownRegion != neiRegion)
381 {
382 addToInterface
383 (
384 mesh,
385 (useFaceZones ? mesh.faceZones().whichZone(facei) : -1),
386 ownRegion,
387 neiRegion,
388 regionsToSize
389 );
390 }
391 }
392
393
394 if (Pstream::parRun())
395 {
396 if (Pstream::master())
397 {
398 // Receive and add to my sizes
399 for (const int slave : Pstream::subProcs())
400 {
401 IPstream fromSlave(Pstream::commsTypes::blocking, slave);
402
403 EdgeMap<Map<label>> slaveSizes(fromSlave);
404
405 forAllConstIters(slaveSizes, slaveIter)
406 {
407 const Map<label>& slaveInfo = *slaveIter;
408
409 auto masterIter = regionsToSize.find(slaveIter.key());
410
411 if (masterIter.found())
412 {
413 // Same inter-region
414 Map<label>& masterInfo = *masterIter;
415
416 forAllConstIters(slaveInfo, iter)
417 {
418 const label zoneID = iter.key();
419 const label slaveSize = iter.val();
420
421 auto zoneIter = masterInfo.find(zoneID);
422 if (zoneIter.found())
423 {
424 *zoneIter += slaveSize;
425 }
426 else
427 {
428 masterInfo.insert(zoneID, slaveSize);
429 }
430 }
431 }
432 else
433 {
434 regionsToSize.insert(slaveIter.key(), slaveInfo);
435 }
436 }
437 }
438 }
439 else
440 {
441 // Send to master
442 {
443 OPstream toMaster
444 (
445 Pstream::commsTypes::blocking,
446 Pstream::masterNo()
447 );
448 toMaster << regionsToSize;
449 }
450 }
451 }
452
453 // Rework
454
455 Pstream::broadcast(regionsToSize);
456
457
458 // Now we have the global sizes of all inter-regions.
459 // Invert this on master and distribute.
460 label nInterfaces = 0;
461 forAllConstIters(regionsToSize, iter)
462 {
463 const Map<label>& info = iter.val();
464 nInterfaces += info.size();
465 }
466
467 interfaces.setSize(nInterfaces);
468 interfaceNames.setSize(nInterfaces);
469 interfaceSizes.setSize(nInterfaces);
470 EdgeMap<Map<label>> regionsToInterface(nInterfaces);
471
472 nInterfaces = 0;
473 forAllConstIters(regionsToSize, iter)
474 {
475 const edge& e = iter.key();
476 const Map<label>& info = iter.val();
477
478 const word& name0 = regionNames[e[0]];
479 const word& name1 = regionNames[e[1]];
480
481 forAllConstIters(info, infoIter)
482 {
483 interfaces[nInterfaces] = iter.key();
484 label zoneID = infoIter.key();
485 if (zoneID == -1)
486 {
487 interfaceNames[nInterfaces] = Pair<word>
488 (
489 name0 + "_to_" + name1,
490 name1 + "_to_" + name0
491 );
492 }
493 else
494 {
495 const word& zoneName = mesh.faceZones()[zoneID].name();
496 interfaceNames[nInterfaces] = Pair<word>
497 (
498 zoneName + "_" + name0 + "_to_" + name1,
499 zoneName + "_" + name1 + "_to_" + name0
500 );
501 }
502 interfaceSizes[nInterfaces] = infoIter();
503
504 if (regionsToInterface.found(e))
505 {
506 regionsToInterface[e].insert(zoneID, nInterfaces);
507 }
508 else
509 {
510 Map<label> zoneAndInterface;
511 zoneAndInterface.insert(zoneID, nInterfaces);
512 regionsToInterface.insert(e, zoneAndInterface);
513 }
514 nInterfaces++;
515 }
516 }
517
518
519 // Consistent interface information for all processors
520 Pstream::broadcasts
521 (
522 UPstream::worldComm,
523 interfaces,
524 interfaceNames,
525 interfaceSizes,
526 regionsToInterface
527 );
528
529 // Mark all inter-region faces.
530 faceToInterface.setSize(mesh.nFaces(), -1);
531
532 forAll(mesh.faceNeighbour(), facei)
533 {
534 label ownRegion = cellRegion[mesh.faceOwner()[facei]];
535 label neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
536
537 if (ownRegion != neiRegion)
538 {
539 label zoneID = -1;
540 if (useFaceZones)
541 {
542 zoneID = mesh.faceZones().whichZone(facei);
543 }
544
546 (
547 min(ownRegion, neiRegion),
548 max(ownRegion, neiRegion)
549 );
550
551 faceToInterface[facei] = regionsToInterface[interface][zoneID];
552 }
553 }
554 forAll(coupledRegion, i)
555 {
556 label facei = i+mesh.nInternalFaces();
557 label ownRegion = cellRegion[mesh.faceOwner()[facei]];
558 label neiRegion = coupledRegion[i];
559
560 if (ownRegion != neiRegion)
561 {
562 label zoneID = -1;
563 if (useFaceZones)
564 {
565 zoneID = mesh.faceZones().whichZone(facei);
566 }
567
569 (
570 min(ownRegion, neiRegion),
571 max(ownRegion, neiRegion)
572 );
573
574 faceToInterface[facei] = regionsToInterface[interface][zoneID];
575 }
576 }
577}
578
579
580// Create mesh for region.
581autoPtr<mapPolyMesh> createRegionMesh
582(
583 const fvMesh& mesh,
584 // Region info
585 const labelList& cellRegion,
586 const label regionI,
587 const word& regionName,
588 // Interface info
589 const labelList& interfacePatches,
590 const labelList& faceToInterface,
591
592 autoPtr<fvMesh>& newMesh
593)
594{
595 // Create dummy system/fv*
596 fvMeshTools::createDummyFvMeshFiles(mesh, regionName, true);
597
598 // Neighbour cellRegion.
599 labelList coupledRegion(mesh.nBoundaryFaces());
600
601 forAll(coupledRegion, i)
602 {
603 label celli = mesh.faceOwner()[i+mesh.nInternalFaces()];
604 coupledRegion[i] = cellRegion[celli];
605 }
606 syncTools::swapBoundaryFaceList(mesh, coupledRegion);
607
608
609 // Topology change container. Start off from existing mesh.
610 polyTopoChange meshMod(mesh);
611
612 // Cell remover engine
613 removeCells cellRemover(mesh);
614
615 // Select all but region cells
616 labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
617
618 // Find out which faces will get exposed. Note that this
619 // gets faces in mesh face order. So both regions will get same
620 // face in same order (important!)
621 labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
622
623 labelList exposedPatchIDs(exposedFaces.size());
624 forAll(exposedFaces, i)
625 {
626 label facei = exposedFaces[i];
627 label interfacei = faceToInterface[facei];
628
629 label ownRegion = cellRegion[mesh.faceOwner()[facei]];
630 label neiRegion = -1;
631
632 if (mesh.isInternalFace(facei))
633 {
634 neiRegion = cellRegion[mesh.faceNeighbour()[facei]];
635 }
636 else
637 {
638 neiRegion = coupledRegion[facei-mesh.nInternalFaces()];
639 }
640
641
642 // Check which side is being kept - determines which of the two
643 // patches will be used.
644
645 label otherRegion = -1;
646
647 if (ownRegion == regionI && neiRegion != regionI)
648 {
649 otherRegion = neiRegion;
650 }
651 else if (ownRegion != regionI && neiRegion == regionI)
652 {
653 otherRegion = ownRegion;
654 }
655 else
656 {
658 << "Exposed face:" << facei
659 << " fc:" << mesh.faceCentres()[facei]
660 << " has owner region " << ownRegion
661 << " and neighbour region " << neiRegion
662 << " when handling region:" << regionI
663 << exit(FatalError);
664 }
665
666 // Find the patch.
667 if (regionI < otherRegion)
668 {
669 exposedPatchIDs[i] = interfacePatches[interfacei];
670 }
671 else
672 {
673 exposedPatchIDs[i] = interfacePatches[interfacei]+1;
674 }
675 }
676
677 // Remove faces
678 cellRemover.setRefinement
679 (
680 cellsToRemove,
681 exposedFaces,
682 exposedPatchIDs,
683 meshMod
684 );
685
686 autoPtr<mapPolyMesh> map = meshMod.makeMesh
687 (
688 newMesh,
690 (
692 mesh.time().timeName(),
693 mesh.time(),
694 IOobject::READ_IF_PRESENT, // read fv* if present
695 IOobject::AUTO_WRITE
696 ),
697 mesh
698 );
699
700 return map;
701}
702
703
704void createAndWriteRegion
705(
706 const fvMesh& mesh,
707 const labelList& cellRegion,
708 const wordList& regionNames,
709 const bool prefixRegion,
710 const labelList& faceToInterface,
711 const labelList& interfacePatches,
712 const label regionI,
713 const word& newMeshInstance
714)
715{
716 Info<< "Creating mesh for region " << regionI
717 << ' ' << regionNames[regionI] << endl;
718
719 autoPtr<fvMesh> newMesh;
720 autoPtr<mapPolyMesh> map = createRegionMesh
721 (
722 mesh,
723 cellRegion,
724 regionI,
725 regionNames[regionI],
726 interfacePatches,
727 faceToInterface,
728 newMesh
729 );
730
731
732 // Make map of all added patches
733 labelHashSet addedPatches(2*interfacePatches.size());
734 forAll(interfacePatches, interfacei)
735 {
736 addedPatches.insert(interfacePatches[interfacei]);
737 addedPatches.insert(interfacePatches[interfacei]+1);
738 }
739
740
741 Info<< "Mapping fields" << endl;
742
743 // Map existing fields
744 newMesh().updateMesh(map());
745
746 // Add subsetted fields
747 subsetVolFields<volScalarField>
748 (
749 mesh,
750 newMesh(),
751 map().cellMap(),
752 map().faceMap(),
753 addedPatches
754 );
755 subsetVolFields<volVectorField>
756 (
757 mesh,
758 newMesh(),
759 map().cellMap(),
760 map().faceMap(),
761 addedPatches
762 );
763 subsetVolFields<volSphericalTensorField>
764 (
765 mesh,
766 newMesh(),
767 map().cellMap(),
768 map().faceMap(),
769 addedPatches
770 );
771 subsetVolFields<volSymmTensorField>
772 (
773 mesh,
774 newMesh(),
775 map().cellMap(),
776 map().faceMap(),
777 addedPatches
778 );
779 subsetVolFields<volTensorField>
780 (
781 mesh,
782 newMesh(),
783 map().cellMap(),
784 map().faceMap(),
785 addedPatches
786 );
787
788 subsetSurfaceFields<surfaceScalarField>
789 (
790 mesh,
791 newMesh(),
792 map().cellMap(),
793 map().faceMap(),
794 addedPatches
795 );
796 subsetSurfaceFields<surfaceVectorField>
797 (
798 mesh,
799 newMesh(),
800 map().cellMap(),
801 map().faceMap(),
802 addedPatches
803 );
804 subsetSurfaceFields<surfaceSphericalTensorField>
805 (
806 mesh,
807 newMesh(),
808 map().cellMap(),
809 map().faceMap(),
810 addedPatches
811 );
812 subsetSurfaceFields<surfaceSymmTensorField>
813 (
814 mesh,
815 newMesh(),
816 map().cellMap(),
817 map().faceMap(),
818 addedPatches
819 );
820 subsetSurfaceFields<surfaceTensorField>
821 (
822 mesh,
823 newMesh(),
824 map().cellMap(),
825 map().faceMap(),
826 addedPatches
827 );
828
829
830 const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
831 newPatches.checkParallelSync(true);
832
833 // Delete empty patches
834 // ~~~~~~~~~~~~~~~~~~~~
835
836 // Create reordering list to move patches-to-be-deleted to end
837 labelList oldToNew(newPatches.size(), -1);
838 DynamicList<label> sharedPatches(newPatches.size());
839 label newI = 0;
840
841 Info<< "Deleting empty patches" << endl;
842
843 // Assumes all non-proc boundaries are on all processors!
844 forAll(newPatches, patchi)
845 {
846 const polyPatch& pp = newPatches[patchi];
847
848 if (!isA<processorPolyPatch>(pp))
849 {
850 if (returnReduce(pp.size(), sumOp<label>()) > 0)
851 {
852 oldToNew[patchi] = newI;
853 if (!addedPatches.found(patchi))
854 {
855 sharedPatches.append(newI);
856 }
857 newI++;
858 }
859 }
860 }
861
862 // Same for processor patches (but need no reduction)
863 forAll(newPatches, patchi)
864 {
865 const polyPatch& pp = newPatches[patchi];
866
867 if (isA<processorPolyPatch>(pp) && pp.size())
868 {
869 oldToNew[patchi] = newI++;
870 }
871 }
872
873 const label nNewPatches = newI;
874
875 // Move all deleteable patches to the end
876 forAll(oldToNew, patchi)
877 {
878 if (oldToNew[patchi] == -1)
879 {
880 oldToNew[patchi] = newI++;
881 }
882 }
883
884 //reorderPatches(newMesh(), oldToNew, nNewPatches);
885 fvMeshTools::reorderPatches(newMesh(), oldToNew, nNewPatches, true);
886
887 // Rename shared patches with region name
888 if (prefixRegion)
889 {
890 Info<< "Prefixing patches with region name" << endl;
891
892 renamePatches(newMesh(), regionNames[regionI], sharedPatches);
893 }
894
895
896 Info<< "Writing new mesh" << endl;
897
898 newMesh().setInstance(newMeshInstance);
899 newMesh().write();
900 topoSet::removeFiles(newMesh());
901 processorMeshes::removeFiles(newMesh());
902
903 // Write addressing files like decomposePar
904 Info<< "Writing addressing to base mesh" << endl;
905
906 labelIOList pointProcAddressing
907 (
909 (
910 "pointRegionAddressing",
911 newMesh().facesInstance(),
912 newMesh().meshSubDir,
913 newMesh(),
914 IOobject::NO_READ,
915 IOobject::NO_WRITE,
916 false
917 ),
918 map().pointMap()
919 );
920 Info<< "Writing map " << pointProcAddressing.name()
921 << " from region" << regionI
922 << " points back to base mesh." << endl;
923 pointProcAddressing.write();
924
925 labelIOList faceProcAddressing
926 (
928 (
929 "faceRegionAddressing",
930 newMesh().facesInstance(),
931 newMesh().meshSubDir,
932 newMesh(),
933 IOobject::NO_READ,
934 IOobject::NO_WRITE,
935 false
936 ),
937 newMesh().nFaces()
938 );
939 forAll(faceProcAddressing, facei)
940 {
941 // face + turning index. (see decomposePar)
942 // Is the face pointing in the same direction?
943 label oldFacei = map().faceMap()[facei];
944
945 if
946 (
947 map().cellMap()[newMesh().faceOwner()[facei]]
948 == mesh.faceOwner()[oldFacei]
949 )
950 {
951 faceProcAddressing[facei] = oldFacei+1;
952 }
953 else
954 {
955 faceProcAddressing[facei] = -(oldFacei+1);
956 }
957 }
958 Info<< "Writing map " << faceProcAddressing.name()
959 << " from region" << regionI
960 << " faces back to base mesh." << endl;
961 faceProcAddressing.write();
962
963 labelIOList cellProcAddressing
964 (
966 (
967 "cellRegionAddressing",
968 newMesh().facesInstance(),
969 newMesh().meshSubDir,
970 newMesh(),
971 IOobject::NO_READ,
972 IOobject::NO_WRITE,
973 false
974 ),
975 map().cellMap()
976 );
977 Info<< "Writing map " <<cellProcAddressing.name()
978 << " from region" << regionI
979 << " cells back to base mesh." << endl;
980 cellProcAddressing.write();
981
982 labelIOList boundaryProcAddressing
983 (
985 (
986 "boundaryRegionAddressing",
987 newMesh().facesInstance(),
988 newMesh().meshSubDir,
989 newMesh(),
990 IOobject::NO_READ,
991 IOobject::NO_WRITE,
992 false
993 ),
994 labelList(nNewPatches, -1)
995 );
996 forAll(oldToNew, i)
997 {
998 if (!addedPatches.found(i))
999 {
1000 label newI = oldToNew[i];
1001 if (newI >= 0 && newI < nNewPatches)
1002 {
1003 boundaryProcAddressing[oldToNew[i]] = i;
1004 }
1005 }
1006 }
1007 Info<< "Writing map " << boundaryProcAddressing.name()
1008 << " from region" << regionI
1009 << " boundary back to base mesh." << endl;
1010 boundaryProcAddressing.write();
1011}
1012
1013
1014// Create for every region-region interface with non-zero size two patches.
1015// First one is for minimumregion to maximumregion.
1016// Note that patches get created in same order on all processors (if parallel)
1017// since looping over synchronised 'interfaces'.
1018labelList addRegionPatches
1019(
1020 fvMesh& mesh,
1021 const wordList& regionNames,
1022 const edgeList& interfaces,
1023 const List<Pair<word>>& interfaceNames
1024)
1025{
1026 Info<< nl << "Adding patches" << nl << endl;
1027
1028 labelList interfacePatches(interfaces.size());
1029
1030 forAll(interfaces, interI)
1031 {
1032 const edge& e = interfaces[interI];
1033 const Pair<word>& names = interfaceNames[interI];
1034
1035 //Info<< "For interface " << interI
1036 // << " between regions " << e
1037 // << " trying to add patches " << names << endl;
1038
1039
1040 mappedWallPolyPatch patch1
1041 (
1042 names[0],
1043 0, // overridden
1044 0, // overridden
1045 0, // overridden
1046 regionNames[e[1]], // sampleRegion
1047 mappedPatchBase::NEARESTPATCHFACE,
1048 names[1], // samplePatch
1049 point::zero, // offset
1050 mesh.boundaryMesh()
1051 );
1052
1053 interfacePatches[interI] = fvMeshTools::addPatch
1054 (
1055 mesh,
1056 patch1,
1057 dictionary(), //optional per field value
1059 true //validBoundary
1060 );
1061
1062 mappedWallPolyPatch patch2
1063 (
1064 names[1],
1065 0,
1066 0,
1067 0,
1068 regionNames[e[0]], // sampleRegion
1069 mappedPatchBase::NEARESTPATCHFACE,
1070 names[0],
1071 point::zero, // offset
1072 mesh.boundaryMesh()
1073 );
1074 fvMeshTools::addPatch
1075 (
1076 mesh,
1077 patch2,
1078 dictionary(), //optional per field value
1080 true //validBoundary
1081 );
1082
1083 Info<< "For interface between region " << regionNames[e[0]]
1084 << " and " << regionNames[e[1]] << " added patches" << endl
1085 << " " << interfacePatches[interI]
1086 << "\t" << mesh.boundaryMesh()[interfacePatches[interI]].name()
1087 << endl
1088 << " " << interfacePatches[interI]+1
1089 << "\t" << mesh.boundaryMesh()[interfacePatches[interI]+1].name()
1090 << endl;
1091 }
1092 return interfacePatches;
1093}
1094
1095
1096// Find region that covers most of cell zone
1097label findCorrespondingRegion
1098(
1099 const labelList& existingZoneID, // per cell the (unique) zoneID
1100 const labelList& cellRegion,
1101 const label nCellRegions,
1102 const label zoneI,
1103 const label minOverlapSize
1104)
1105{
1106 // Per region the number of cells in zoneI
1107 labelList cellsInZone(nCellRegions, Zero);
1108
1109 forAll(cellRegion, celli)
1110 {
1111 if (existingZoneID[celli] == zoneI)
1112 {
1113 cellsInZone[cellRegion[celli]]++;
1114 }
1115 }
1116
1117 Pstream::listCombineAllGather(cellsInZone, plusEqOp<label>());
1118
1119 // Pick region with largest overlap of zoneI
1120 label regionI = findMax(cellsInZone);
1121
1122
1123 if (cellsInZone[regionI] < minOverlapSize)
1124 {
1125 // Region covers too little of zone. Not good enough.
1126 regionI = -1;
1127 }
1128 else
1129 {
1130 // Check that region contains no cells that aren't in cellZone.
1131 forAll(cellRegion, celli)
1132 {
1133 if (cellRegion[celli] == regionI && existingZoneID[celli] != zoneI)
1134 {
1135 // celli in regionI but not in zoneI
1136 regionI = -1;
1137 break;
1138 }
1139 }
1140 // If one in error, all should be in error. Note that branch gets taken
1141 // on all procs.
1142 reduce(regionI, minOp<label>());
1143 }
1144
1145 return regionI;
1146}
1147
1148
1149void getClusterID
1150(
1151 const polyMesh& mesh,
1152 const cellZoneMesh& cellZones,
1153 const wordList& clusterNames,
1154 const labelListList& clusterToZones,
1155 labelList& clusterID,
1156 labelList& neiClusterID
1157)
1158{
1159 // Existing zoneID
1160 clusterID.setSize(mesh.nCells());
1161 clusterID = -1;
1162
1163 forAll(clusterToZones, clusterI)
1164 {
1165 for (const label zoneI : clusterToZones[clusterI])
1166 {
1167 const cellZone& cz = cellZones[zoneI];
1168
1169 forAll(cz, i)
1170 {
1171 label celli = cz[i];
1172 if (clusterID[celli] == -1)
1173 {
1174 clusterID[celli] = clusterI;
1175 }
1176 else
1177 {
1179 << "Cell " << celli << " with cell centre "
1180 << mesh.cellCentres()[celli]
1181 << " is multiple zones. This is not allowed." << endl
1182 << "It is in zone " << clusterNames[clusterID[celli]]
1183 << " and in zone " << clusterNames[clusterI]
1184 << exit(FatalError);
1185 }
1186 }
1187 }
1188 }
1189
1190 // Neighbour zoneID.
1191 syncTools::swapBoundaryCellList(mesh, clusterID, neiClusterID);
1192}
1193
1194
1195word makeRegionName
1196(
1197 const cellZoneMesh& czs,
1198 const label regioni,
1199 const labelList& zoneIDs
1200)
1201{
1202 // Synthesise region name. Equals the zone name if cluster consist of only
1203 // one zone
1204
1205 if (zoneIDs.empty())
1206 {
1207 return word("domain") + Foam::name(regioni);
1208 }
1209 else
1210 {
1211 // Zone indices are in cellZone order ...
1212 word regionName(czs[zoneIDs[0]].name());
1213
1214 // Synthesize name from appended cellZone names
1215 for (label i = 1; i < zoneIDs.size(); i++)
1216 {
1217 regionName += "_" + czs[zoneIDs[i]].name();
1218 }
1219 return regionName;
1220 }
1221}
1222
1223
1224void makeClusters
1225(
1226 const List<wordRes>& zoneClusters,
1227 const wordList& zoneClusterNames,
1228 const cellZoneMesh& cellZones,
1229 wordList& clusterNames,
1230 labelListList& clusterToZones,
1231 labelList& zoneToCluster
1232)
1233{
1234 // Check if there are clustering for zones. If none every zone goes into
1235 // its own cluster.
1236
1237 clusterNames.clear();
1238 clusterToZones.clear();
1239 zoneToCluster.setSize(cellZones.size());
1240 zoneToCluster = -1;
1241
1242 if (zoneClusters.size())
1243 {
1244 forAll(zoneClusters, clusteri)
1245 {
1246 const labelList zoneIDs(cellZones.indices(zoneClusters[clusteri]));
1247 UIndirectList<label>(zoneToCluster, zoneIDs) = clusteri;
1248 clusterNames.append
1249 (
1250 zoneClusterNames[clusteri].size()
1251 ? zoneClusterNames[clusteri]
1252 : makeRegionName
1253 (
1254 cellZones,
1255 clusteri,
1256 zoneIDs
1257 )
1258 );
1259 clusterToZones.append(std::move(zoneIDs));
1260 }
1261
1262 // Unclustered zone
1263 forAll(zoneToCluster, zonei)
1264 {
1265 if (zoneToCluster[zonei] == -1)
1266 {
1267 clusterNames.append(cellZones[zonei].name());
1268 clusterToZones.append(labelList(1, zonei));
1269 zoneToCluster[zonei] = clusterToZones.size();
1270 }
1271 }
1272 }
1273 else
1274 {
1275 for (const auto& cellZone : cellZones)
1276 {
1277 const label nClusters = clusterToZones.size();
1278 clusterNames.append(cellZone.name());
1279 clusterToZones.append(labelList(1, cellZone.index()));
1280 zoneToCluster[cellZone.index()] = nClusters;
1281 }
1282 }
1283}
1284
1285
1286void matchRegions
1287(
1288 const bool sloppyCellZones,
1289 const polyMesh& mesh,
1290
1291 const wordList& clusterNames,
1292 const labelListList& clusterToZones,
1293 const labelList& clusterID,
1294
1295 const label nCellRegions,
1296 const labelList& cellRegion,
1297
1298 labelListList& regionToZones,
1300 labelList& zoneToRegion
1301)
1302{
1303 const cellZoneMesh& cellZones = mesh.cellZones();
1304
1305 regionToZones.setSize(nCellRegions);
1306 regionToZones = labelList();
1307 regionNames.setSize(nCellRegions);
1309 zoneToRegion.setSize(cellZones.size(), -1);
1310
1311
1312 // Sizes per cluster
1313 labelList clusterSizes(clusterToZones.size(), Zero);
1314 forAll(clusterToZones, clusterI)
1315 {
1316 for (const label zoneI : clusterToZones[clusterI])
1317 {
1318 clusterSizes[clusterI] += cellZones[zoneI].size();
1319 }
1320 reduce(clusterSizes[clusterI], sumOp<label>());
1321 }
1322
1323 if (sloppyCellZones)
1324 {
1325 Info<< "Trying to match regions to existing cell zones;"
1326 << " region can be subset of cell zone." << nl << endl;
1327
1328 forAll(clusterToZones, clusterI)
1329 {
1330 label regionI = findCorrespondingRegion
1331 (
1332 clusterID,
1333 cellRegion,
1334 nCellRegions,
1335 clusterI,
1336 label(0.5*clusterSizes[clusterI]) // minimum overlap
1337 );
1338
1339 if (regionI != -1)
1340 {
1341 Info<< "Sloppily matched region " << regionI
1342 //<< " size " << regionSizes[regionI]
1343 << " to cluster " << clusterI
1344 << " size " << clusterSizes[clusterI]
1345 << endl;
1347 (
1348 zoneToRegion,
1349 clusterToZones[clusterI]
1350 ) = regionI;
1351 regionToZones[regionI] = clusterToZones[clusterI];
1352 regionNames[regionI] = clusterNames[clusterI];
1353 }
1354 }
1355 }
1356 else
1357 {
1358 Info<< "Trying to match regions to existing cell zones." << nl << endl;
1359
1360 forAll(clusterToZones, clusterI)
1361 {
1362 label regionI = findCorrespondingRegion
1363 (
1364 clusterID,
1365 cellRegion,
1366 nCellRegions,
1367 clusterI,
1368 clusterSizes[clusterI] // require exact match
1369 );
1370
1371 if (regionI != -1)
1372 {
1374 (
1375 zoneToRegion,
1376 clusterToZones[clusterI]
1377 ) = regionI;
1378 regionToZones[regionI] = clusterToZones[clusterI];
1379 regionNames[regionI] = clusterNames[clusterI];
1380 }
1381 }
1382 }
1383 // Allocate region names for unmatched regions.
1384 forAll(regionNames, regionI)
1385 {
1386 if (regionNames[regionI].empty())
1387 {
1388 regionNames[regionI] = makeRegionName
1389 (
1390 cellZones,
1391 regionI,
1392 regionToZones[regionI]
1393 );
1394 }
1395 }
1396}
1397
1398
1399void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
1400{
1401 // Write to manual decomposition option
1402 {
1403 labelIOList cellToRegion
1404 (
1405 IOobject
1406 (
1407 "cellToRegion",
1408 mesh.facesInstance(),
1409 mesh,
1410 IOobject::NO_READ,
1411 IOobject::NO_WRITE,
1412 false
1413 ),
1414 cellRegion
1415 );
1416 cellToRegion.write();
1417
1418 Info<< "Writing region per cell file (for manual decomposition) to "
1419 << cellToRegion.objectPath() << nl << endl;
1420 }
1421 // Write for postprocessing
1422 {
1423 volScalarField cellToRegion
1424 (
1425 IOobject
1426 (
1427 "cellToRegion",
1428 mesh.time().timeName(),
1429 mesh,
1430 IOobject::NO_READ,
1431 IOobject::NO_WRITE,
1432 false
1433 ),
1434 mesh,
1435 dimensionedScalar(dimless, Zero),
1437 );
1438 forAll(cellRegion, celli)
1439 {
1440 cellToRegion[celli] = cellRegion[celli];
1441 }
1442 cellToRegion.write();
1443
1444 Info<< "Writing region per cell as volScalarField to "
1445 << cellToRegion.objectPath() << nl << endl;
1446 }
1447}
1448
1449
1450
1451
1452int main(int argc, char *argv[])
1453{
1454 argList::addNote
1455 (
1456 "Split mesh into multiple regions (detected by walking across faces)"
1457 );
1458 #include "addRegionOption.H"
1459 #include "addOverwriteOption.H"
1460 argList::addBoolOption
1461 (
1462 "cellZones",
1463 "Additionally split cellZones off into separate regions"
1464 );
1465 argList::addBoolOption
1466 (
1467 "cellZonesOnly",
1468 "Use cellZones only to split mesh into regions; do not use walking"
1469 );
1470 argList::addOption
1471 (
1472 "cellZonesFileOnly",
1473 "file",
1474 "Like -cellZonesOnly, but use specified file"
1475 );
1476 argList::addOption
1477 (
1478 "combineZones",
1479 "lists of zones",
1480 "Combine zones in follow-on analysis"
1481 );
1482 argList::addOption
1483 (
1484 "addZones",
1485 "lists of zones",
1486 "Combine zones in follow-on analysis"
1487 );
1488 argList::addOption
1489 (
1490 "blockedFaces",
1491 "faceSet",
1492 "Specify additional region boundaries that walking does not cross"
1493 );
1494 argList::addBoolOption
1495 (
1496 "makeCellZones",
1497 "Place cells into cellZones instead of splitting mesh"
1498 );
1499 argList::addBoolOption
1500 (
1501 "largestOnly",
1502 "Only write largest region"
1503 );
1504 argList::addOption
1505 (
1506 "insidePoint",
1507 "point",
1508 "Only write region containing point"
1509 );
1510 argList::addBoolOption
1511 (
1512 "detectOnly",
1513 "Do not write mesh"
1514 );
1515 argList::addBoolOption
1516 (
1517 "sloppyCellZones",
1518 "Try to match heuristically regions to existing cell zones"
1519 );
1520 argList::addBoolOption
1521 (
1522 "useFaceZones",
1523 "Use faceZones to patch inter-region faces instead of single patch"
1524 );
1525 argList::addBoolOption
1526 (
1527 "prefixRegion",
1528 "Prefix region name to all patches, not just coupling patches"
1529 );
1530
1531 argList::noFunctionObjects(); // Never use function objects
1532
1533 #include "setRootCase.H"
1534 #include "createTime.H"
1535 #include "createNamedMesh.H"
1536
1537 const word oldInstance = mesh.pointsInstance();
1538
1539 word blockedFacesName;
1540 if (args.readIfPresent("blockedFaces", blockedFacesName))
1541 {
1542 Info<< "Reading blocked internal faces from faceSet "
1543 << blockedFacesName << nl << endl;
1544 }
1545
1546 const bool makeCellZones = args.found("makeCellZones");
1547 const bool largestOnly = args.found("largestOnly");
1548 const bool insidePoint = args.found("insidePoint");
1549 const bool useCellZones = args.found("cellZones");
1550 const bool useCellZonesOnly = args.found("cellZonesOnly");
1551 const bool useCellZonesFile = args.found("cellZonesFileOnly");
1552 const bool combineZones = args.found("combineZones");
1553 const bool addZones = args.found("addZones");
1554 const bool overwrite = args.found("overwrite");
1555 const bool detectOnly = args.found("detectOnly");
1556 const bool sloppyCellZones = args.found("sloppyCellZones");
1557 const bool useFaceZones = args.found("useFaceZones");
1558 const bool prefixRegion = args.found("prefixRegion");
1559
1560
1561 if
1562 (
1563 (useCellZonesOnly || useCellZonesFile)
1564 && (useCellZones || blockedFacesName.size())
1565 )
1566 {
1568 << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
1569 << " (which specify complete split)"
1570 << " in combination with -blockedFaces or -cellZones"
1571 << " (which imply a split based on topology)"
1572 << exit(FatalError);
1573 }
1574
1575
1576 if (useFaceZones)
1577 {
1578 Info<< "Using current faceZones to divide inter-region interfaces"
1579 << " into multiple patches."
1580 << nl << endl;
1581 }
1582 else
1583 {
1584 Info<< "Creating single patch per inter-region interface."
1585 << nl << endl;
1586 }
1587
1588
1589
1590 if (insidePoint && largestOnly)
1591 {
1593 << "You cannot specify both -largestOnly"
1594 << " (keep region with most cells)"
1595 << " and -insidePoint (keep region containing point)"
1596 << exit(FatalError);
1597 }
1598
1599
1600 // Make sure cellZone names consistent across processors
1601 mesh.cellZones().checkParallelSync(true);
1602
1603 List<wordRes> zoneClusters;
1604 wordList zoneClusterNames;
1605 if (combineZones)
1606 {
1607 if (addZones)
1608 {
1610 << "Cannot specify both combineZones and addZones"
1611 << exit(FatalError);
1612 }
1613 zoneClusters = args.get<List<wordRes>>("combineZones");
1614 zoneClusterNames.setSize(zoneClusters.size());
1615 }
1616 else if (addZones)
1617 {
1618 zoneClusters = args.get<List<wordRes>>("addZones");
1619 zoneClusterNames.setSize(zoneClusters.size());
1620 forAll(zoneClusters, clusteri)
1621 {
1622 // Pop of front - is name
1623
1624 wordRes& wrs = zoneClusters[clusteri];
1625
1626 zoneClusterNames[clusteri] = wrs[0];
1627
1628 for (label i = 1; i < wrs.size(); i++)
1629 {
1630 wrs[i-1] = wrs[i];
1631 }
1632 wrs.setSize(wrs.size()-1);
1633 }
1634 }
1635
1636
1637 // Determine per cell the region it belongs to
1638 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1639
1640 // cellRegion is the labelList with the region per cell.
1641 labelList cellRegion;
1642 // Region to zone(s)
1643 labelListList regionToZones;
1644 // Name of region
1646 // Zone to region
1647 labelList zoneToRegion;
1648
1649 label nCellRegions = 0;
1650 if (useCellZonesOnly)
1651 {
1652 Info<< "Using current cellZones to split mesh into regions."
1653 << " This requires all"
1654 << " cells to be in one and only one cellZone." << nl << endl;
1655
1656 // Collect sets of zones into clusters. If no cluster is just an identity
1657 // list (cluster 0 is cellZone 0 etc.)
1658 wordList clusterNames;
1659 labelListList clusterToZones;
1660 labelList zoneToCluster;
1661 makeClusters
1662 (
1663 zoneClusters,
1664 zoneClusterNames,
1665 mesh.cellZones(),
1666 clusterNames,
1667 clusterToZones,
1668 zoneToCluster
1669 );
1670
1671 // Existing clusterID
1672 labelList clusterID(mesh.nCells(), -1);
1673 // Neighbour clusterID.
1674 labelList neiClusterID(mesh.nBoundaryFaces());
1675 getClusterID
1676 (
1677 mesh,
1678 mesh.cellZones(),
1679 clusterNames,
1680 clusterToZones,
1681 clusterID,
1682 neiClusterID
1683 );
1684
1685 label unzonedCelli = clusterID.find(-1);
1686 if (unzonedCelli != -1)
1687 {
1689 << "For the cellZonesOnly option all cells "
1690 << "have to be in a cellZone." << endl
1691 << "Cell " << unzonedCelli
1692 << " at" << mesh.cellCentres()[unzonedCelli]
1693 << " is not in a cellZone. There might be more unzoned cells."
1694 << exit(FatalError);
1695 }
1696 cellRegion = clusterID;
1697 nCellRegions = gMax(cellRegion)+1;
1698 zoneToRegion = zoneToCluster;
1699 regionToZones = clusterToZones;
1700 regionNames = clusterNames;
1701 }
1702 else if (useCellZonesFile)
1703 {
1704 const word zoneFile(args["cellZonesFileOnly"]);
1705 Info<< "Reading split from cellZones file " << zoneFile << endl
1706 << "This requires all"
1707 << " cells to be in one and only one cellZone." << nl << endl;
1708
1709 cellZoneMesh newCellZones
1710 (
1711 IOobject
1712 (
1713 zoneFile,
1714 mesh.facesInstance(),
1715 polyMesh::meshSubDir,
1716 mesh,
1717 IOobject::MUST_READ,
1718 IOobject::NO_WRITE,
1719 false
1720 ),
1721 mesh
1722 );
1723
1724 wordList clusterNames;
1725 labelListList clusterToZones;
1726 labelList zoneToCluster;
1727 makeClusters
1728 (
1729 zoneClusters,
1730 zoneClusterNames,
1731 newCellZones,
1732 clusterNames,
1733 clusterToZones,
1734 zoneToCluster
1735 );
1736
1737
1738 // Existing clusterID
1739 labelList clusterID(mesh.nCells(), -1);
1740 // Neighbour clusterID.
1741 labelList neiClusterID(mesh.nBoundaryFaces());
1742 getClusterID
1743 (
1744 mesh,
1745 newCellZones,
1746 clusterNames,
1747 clusterToZones,
1748 clusterID,
1749 neiClusterID
1750 );
1751
1752
1753 label unzonedCelli = clusterID.find(-1);
1754 if (unzonedCelli != -1)
1755 {
1757 << "For the cellZonesFileOnly option all cells "
1758 << "have to be in a cellZone." << endl
1759 << "Cell " << unzonedCelli
1760 << " at" << mesh.cellCentres()[unzonedCelli]
1761 << " is not in a cellZone. There might be more unzoned cells."
1762 << exit(FatalError);
1763 }
1764 cellRegion = clusterID;
1765 nCellRegions = gMax(cellRegion)+1;
1766 zoneToRegion = zoneToCluster;
1767 regionToZones = clusterToZones;
1768 regionNames = clusterNames;
1769 }
1770 else
1771 {
1772 // Determine connected regions
1773 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1774
1775 // Mark additional faces that are blocked
1776 boolList blockedFace;
1777
1778 // Read from faceSet
1779 if (blockedFacesName.size())
1780 {
1781 faceSet blockedFaceSet(mesh, blockedFacesName);
1782 Info<< "Read "
1783 << returnReduce(blockedFaceSet.size(), sumOp<label>())
1784 << " blocked faces from set " << blockedFacesName << nl << endl;
1785
1786 blockedFace.setSize(mesh.nFaces(), false);
1787
1788 for (const label facei : blockedFaceSet)
1789 {
1790 blockedFace[facei] = true;
1791 }
1792 }
1793
1794 // Collect sets of zones into clusters. If no cluster is just an
1795 // identity list (cluster 0 is cellZone 0 etc.)
1796 wordList clusterNames;
1797 labelListList clusterToZones;
1798 labelList zoneToCluster;
1799 makeClusters
1800 (
1801 zoneClusters,
1802 zoneClusterNames,
1803 mesh.cellZones(),
1804 clusterNames,
1805 clusterToZones,
1806 zoneToCluster
1807 );
1808
1809 // Existing clusterID
1810 labelList clusterID(mesh.nCells(), -1);
1811 // Neighbour clusterID.
1812 labelList neiClusterID(mesh.nBoundaryFaces());
1813 getClusterID
1814 (
1815 mesh,
1816 mesh.cellZones(),
1817 clusterNames,
1818 clusterToZones,
1819 clusterID,
1820 neiClusterID
1821 );
1822
1823
1824 // Imply from differing cellZones
1825 if (useCellZones)
1826 {
1827 blockedFace.setSize(mesh.nFaces(), false);
1828
1829 for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
1830 {
1831 label ownCluster = clusterID[mesh.faceOwner()[facei]];
1832 label neiCluster = clusterID[mesh.faceNeighbour()[facei]];
1833
1834 if (ownCluster != neiCluster)
1835 {
1836 blockedFace[facei] = true;
1837 }
1838 }
1839
1840 // Different cellZones on either side of processor patch.
1841 forAll(neiClusterID, i)
1842 {
1843 label facei = i+mesh.nInternalFaces();
1844 label ownCluster = clusterID[mesh.faceOwner()[facei]];
1845 label neiCluster = neiClusterID[i];
1846
1847 if (ownCluster != neiCluster)
1848 {
1849 blockedFace[facei] = true;
1850 }
1851 }
1852 }
1853
1854 // Do a topological walk to determine regions
1855 regionSplit regions(mesh, blockedFace);
1856 nCellRegions = regions.nRegions();
1857 cellRegion.transfer(regions);
1858
1859 // Match regions to zones
1860 matchRegions
1861 (
1862 sloppyCellZones,
1863 mesh,
1864
1865 // cluster-to-zone and cluster-to-cell addressing
1866 clusterNames,
1867 clusterToZones,
1868 clusterID,
1869
1870 // cell-to-region addressing
1871 nCellRegions,
1872 cellRegion,
1873
1874 // matched zones
1875 regionToZones,
1877 zoneToRegion
1878 );
1879
1880 // Override any default region names if single region selected
1881 if (largestOnly || insidePoint)
1882 {
1883 forAll(regionToZones, regionI)
1884 {
1885 if (regionToZones[regionI].empty())
1886 {
1887 if (overwrite)
1888 {
1889 regionNames[regionI] = polyMesh::defaultRegion;
1890 }
1891 else if (insidePoint)
1892 {
1893 regionNames[regionI] = "insidePoint";
1894 }
1895 else if (largestOnly)
1896 {
1897 regionNames[regionI] = "largestOnly";
1898 }
1899 }
1900 }
1901 }
1902 }
1903
1904 Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
1905
1906
1907 // Write decomposition to file
1908 writeCellToRegion(mesh, cellRegion);
1909
1910
1911
1912 // Sizes per region
1913 // ~~~~~~~~~~~~~~~~
1914
1915 labelList regionSizes(nCellRegions, Zero);
1916
1917 forAll(cellRegion, celli)
1918 {
1919 regionSizes[cellRegion[celli]]++;
1920 }
1921 forAll(regionSizes, regionI)
1922 {
1923 reduce(regionSizes[regionI], sumOp<label>());
1924 }
1925
1926 Info<< "Region\tCells" << nl
1927 << "------\t-----" << endl;
1928
1929 forAll(regionSizes, regionI)
1930 {
1931 Info<< regionI << "\t\t" << regionSizes[regionI] << nl;
1932 }
1933 Info<< endl;
1934
1935
1936
1937 // Print region to zone
1938 Info<< "Region\tZone\tName" << nl
1939 << "------\t----\t----" << endl;
1940 forAll(regionToZones, regionI)
1941 {
1942 Info<< regionI << "\t\t" << flatOutput(regionToZones[regionI])
1943 << '\t'
1944 << regionNames[regionI] << nl;
1945 }
1946 Info<< endl;
1947
1948
1949
1950 // Since we're going to mess with patches and zones make sure all
1951 // is synchronised
1952 mesh.boundaryMesh().checkParallelSync(true);
1953 mesh.faceZones().checkParallelSync(true);
1954
1955
1956 // Interfaces
1957 // ----------
1958 // per interface:
1959 // - the two regions on either side
1960 // - the name
1961 // - the (global) size
1962 edgeList interfaces;
1963 List<Pair<word>> interfaceNames;
1964 labelList interfaceSizes;
1965 // per face the interface
1966 labelList faceToInterface;
1967
1968 getInterfaceSizes
1969 (
1970 mesh,
1971 useFaceZones,
1972 cellRegion,
1974
1975 interfaces,
1976 interfaceNames,
1977 interfaceSizes,
1978 faceToInterface
1979 );
1980
1981 Info<< "Sizes of interfaces between regions:" << nl << nl
1982 << "Interface\tRegion\tRegion\tFaces" << nl
1983 << "---------\t------\t------\t-----" << endl;
1984
1985 forAll(interfaces, interI)
1986 {
1987 const edge& e = interfaces[interI];
1988
1989 Info<< interI
1990 << "\t\t\t" << e[0] << "\t\t" << e[1]
1991 << "\t\t" << interfaceSizes[interI] << nl;
1992 }
1993 Info<< endl;
1994
1995
1996 if (detectOnly)
1997 {
1998 Info<< "End\n" << endl;
1999
2000 return 0;
2001 }
2002
2003
2004 // Read objects in time directory
2005 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2006
2007 IOobjectList objects(mesh, runTime.timeName());
2008
2009 // Read vol fields.
2010
2012 ReadFields(mesh, objects, vsFlds);
2013
2015 ReadFields(mesh, objects, vvFlds);
2016
2018 ReadFields(mesh, objects, vstFlds);
2019
2021 ReadFields(mesh, objects, vsymtFlds);
2022
2024 ReadFields(mesh, objects, vtFlds);
2025
2026 // Read surface fields.
2027
2029 ReadFields(mesh, objects, ssFlds);
2030
2032 ReadFields(mesh, objects, svFlds);
2033
2035 ReadFields(mesh, objects, sstFlds);
2036
2038 ReadFields(mesh, objects, ssymtFlds);
2039
2041 ReadFields(mesh, objects, stFlds);
2042
2043 Info<< endl;
2044
2045
2046 // Remove any demand-driven fields ('S', 'V' etc)
2047 mesh.clearOut();
2048
2049
2050 if (nCellRegions == 1)
2051 {
2052 Info<< "Only one region. Doing nothing." << endl;
2053 }
2054 else if (makeCellZones)
2055 {
2056 Info<< "Putting cells into cellZones instead of splitting mesh."
2057 << endl;
2058
2059 // Check if region overlaps with existing zone. If so keep.
2060
2061 for (label regionI = 0; regionI < nCellRegions; regionI++)
2062 {
2063 const labelList& zones = regionToZones[regionI];
2064
2065 if (zones.size() == 1 && zones[0] != -1)
2066 {
2067 // Existing zone
2068 const label zoneI = zones[0];
2069 Info<< " Region " << regionI << " : corresponds to existing"
2070 << " cellZone "
2071 << zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
2072 }
2073 else
2074 {
2075 // Create new cellZone.
2076 const labelList regionCells(findIndices(cellRegion, regionI));
2077
2078 const word zoneName = "region" + Foam::name(regionI);
2079
2080 label zoneI = mesh.cellZones().findZoneID(zoneName);
2081
2082 if (zoneI == -1)
2083 {
2084 zoneI = mesh.cellZones().size();
2085 mesh.cellZones().setSize(zoneI+1);
2086 mesh.cellZones().set
2087 (
2088 zoneI,
2089 new cellZone
2090 (
2091 zoneName, //name
2092 std::move(regionCells), //addressing
2093 zoneI, //index
2094 mesh.cellZones() //cellZoneMesh
2095 )
2096 );
2097 }
2098 else
2099 {
2100 mesh.cellZones()[zoneI].clearAddressing();
2101 mesh.cellZones()[zoneI] = regionCells;
2102 }
2103 Info<< " Region " << regionI << " : created new cellZone "
2104 << zoneI << ' ' << mesh.cellZones()[zoneI].name() << endl;
2105 }
2106 }
2107 mesh.cellZones().writeOpt(IOobject::AUTO_WRITE);
2108
2109 if (!overwrite)
2110 {
2111 ++runTime;
2112 mesh.setInstance(runTime.timeName());
2113 }
2114 else
2115 {
2116 mesh.setInstance(oldInstance);
2117 }
2118
2119 Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
2120 << nl << endl;
2121
2122 mesh.write();
2123
2124
2125 // Write cellSets for convenience
2126 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2127
2128 Info<< "Writing cellSets corresponding to cellZones." << nl << endl;
2129
2130 for (const auto& cz : mesh.cellZones())
2131 {
2132 cellSet(mesh, cz.name(), cz).write();
2133 }
2134 }
2135 else
2136 {
2137 // Add patches for interfaces
2138 // ~~~~~~~~~~~~~~~~~~~~~~~~~~
2139
2140 // Add all possible patches. Empty ones get filtered later on.
2141 Info<< nl << "Adding patches" << nl << endl;
2142
2143 labelList interfacePatches
2144 (
2145 addRegionPatches
2146 (
2147 mesh,
2149 interfaces,
2150 interfaceNames
2151 )
2152 );
2153
2154
2155 if (!overwrite)
2156 {
2157 ++runTime;
2158 }
2159
2160
2161 // Create regions
2162 // ~~~~~~~~~~~~~~
2163
2164 if (insidePoint)
2165 {
2166 const point insidePoint = args.get<point>("insidePoint");
2167
2168 label regionI = -1;
2169
2170 (void)mesh.tetBasePtIs();
2171
2172 label celli = mesh.findCell(insidePoint);
2173
2174 Info<< nl << "Found point " << insidePoint << " in cell " << celli
2175 << endl;
2176
2177 if (celli != -1)
2178 {
2179 regionI = cellRegion[celli];
2180 }
2181
2182 reduce(regionI, maxOp<label>());
2183
2184 Info<< nl
2185 << "Subsetting region " << regionI
2186 << " containing point " << insidePoint << endl;
2187
2188 if (regionI == -1)
2189 {
2191 << "Point " << insidePoint
2192 << " is not inside the mesh." << nl
2193 << "Bounding box of the mesh:" << mesh.bounds()
2194 << exit(FatalError);
2195 }
2196
2197 createAndWriteRegion
2198 (
2199 mesh,
2200 cellRegion,
2202 prefixRegion,
2203 faceToInterface,
2204 interfacePatches,
2205 regionI,
2206 (overwrite ? oldInstance : runTime.timeName())
2207 );
2208 }
2209 else if (largestOnly)
2210 {
2211 label regionI = findMax(regionSizes);
2212
2213 Info<< nl
2214 << "Subsetting region " << regionI
2215 << " of size " << regionSizes[regionI]
2216 << " as named region " << regionNames[regionI] << endl;
2217
2218 createAndWriteRegion
2219 (
2220 mesh,
2221 cellRegion,
2223 prefixRegion,
2224 faceToInterface,
2225 interfacePatches,
2226 regionI,
2227 (overwrite ? oldInstance : runTime.timeName())
2228 );
2229 }
2230 else
2231 {
2232 // Split all
2233 for (label regionI = 0; regionI < nCellRegions; regionI++)
2234 {
2235 Info<< nl
2236 << "Region " << regionI << nl
2237 << "-------- " << endl;
2238
2239 createAndWriteRegion
2240 (
2241 mesh,
2242 cellRegion,
2244 prefixRegion,
2245 faceToInterface,
2246 interfacePatches,
2247 regionI,
2248 (overwrite ? oldInstance : runTime.timeName())
2249 );
2250 }
2251 }
2252 }
2253
2254 Info<< "End\n" << endl;
2255
2256 return 0;
2257}
2258
2259
2260// ************************************************************************* //
Field reading functions for post-processing utilities.
Y[inertIndex] max(0.0)
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
Map from edge (expressed as its endpoints) to value. For easier forward declaration it is currently i...
Definition: EdgeMap.H:54
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
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
iterator find(const Key &key)
Find and return an iterator set at the hashed entry.
Definition: HashTableI.H:114
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
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
Input inter-processor communications stream.
Definition: IPstream.H:57
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
Output inter-processor communications stream.
Definition: OPstream.H:57
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 pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
void append(T *ptr)
Append an element to the end of the list.
Definition: PtrListI.H:113
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
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
labelList indices(const wordRe &matcher, const bool useGroups=true) const
Return (sorted) zone indices for all matches.
Definition: ZoneMesh.C:377
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
A collection of cell labels.
Definition: cellSet.H:54
A subset of mesh cells.
Definition: cellZone.H:65
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A list of face labels.
Definition: faceSet.H:54
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Determines a mapping between patch face centres and mesh cell or face centres and processors they're ...
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
bool checkParallelSync(const bool report=false) const
Check whether all procs have all patches and in same order.
void updateMesh()
Correct polyBoundaryMesh after topology update.
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
Direct mesh changes based on v1.3 polyTopoChange syntax.
virtual bool write(const bool valid=true) const
Write using setting from DB.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:144
Given list of cells to remove, insert all the topology changes.
Definition: removeCells.H:64
A class for managing temporary objects.
Definition: tmp.H:65
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
label index() const noexcept
The index of this zone in the zone list.
const word & name() const noexcept
The zone name.
dynamicFvMesh & mesh
engineTime & runTime
Foam::word regionName(Foam::polyMesh::defaultRegion)
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
wordList regionNames
const labelIOList & zoneIDs
Definition: correctPhi.H:59
const labelIOList & zoneID
#define WarningInFunction
Report a warning using Foam::Warning.
List< word > names(const UPtrList< T > &list, const UnaryMatchPredicate &matcher)
Namespace for OpenFOAM.
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
Definition: FlatOutput.H:215
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.
label findMax(const ListType &input, label start=0)
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
multivariateSurfaceInterpolationScheme< scalar >::fieldTable fields
Definition: createFields.H:97
volScalarField & e
Definition: createFields.H:11
interfaceProperties interface(alpha1, U, thermo->transportPropertiesDict())
#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
static const char *const typeName
The type name used in ensight case files.
A non-counting (dummy) refCount.
Definition: refCount.H:59