1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29 Splits mesh into multiple regions.
31 Each region is defined as a domain whose cells can all be reached by
32 cell-face-cell walking without crossing
34 - additional faces from faceset (-blockedFaces faceSet).
35 - any face inbetween differing cellZones (-cellZones)
38 - volScalarField with regions as different scalars (-detectOnly) or
39 - mesh with multiple regions or
40 - mesh with cells put into cellZones (-makeCellZones)
43 - cellZonesOnly does not do a walk and uses the cellZones only. Use
44 this if you don't mind having disconnected domains in a single region.
45 This option requires all cells to be in one (and one only) cellZone.
47 - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
48 from the specified file. This allows one to explicitly specify the region
49 distribution and still have multiple cellZones per region.
51 - useCellZonesOnly does not do a walk and uses the cellZones only. Use
52 this if you don't mind having disconnected domains in a single region.
53 This option requires all cells to be in one (and one only) cellZone.
56 - Should work in parallel.
57 cellZones can differ on either side of processor boundaries in which case
58 the faces get moved from processor patch to directMapped patch. Not
61 - If a cell zone gets split into more than one region it can detect
62 the largest matching region (-sloppyCellZones). This will accept any
63 region that covers more than 50% of the zone. It has to be a subset
64 so cannot have any cells in any other zone.
66 - writes maps like decomposePar back to original mesh:
67 - pointRegionAddressing : for every point in this region the point in
69 - cellRegionAddressing : ,, cell ,, cell ,,
70 - faceRegionAddressing : ,, face ,, face in
71 the original mesh + 'turning index'. For a face in the same orientation
72 this is the original facelabel+1, for a turned face this is -facelabel-1
73 \*---------------------------------------------------------------------------*/
75 #include "SortableList.H"
77 #include "regionSplit.H"
78 #include "fvMeshSubset.H"
79 #include "IOobjectList.H"
80 #include "volFields.H"
83 #include "polyTopoChange.H"
84 #include "removeCells.H"
86 #include "syncTools.H"
87 #include "ReadFields.H"
88 #include "directMappedWallPolyPatch.H"
89 #include "zeroGradientFvPatchFields.H"
93 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
95 template<class GeoField>
96 void addPatchFields(fvMesh& mesh, const word& patchFieldType)
98 HashTable<const GeoField*> flds
100 mesh.objectRegistry::lookupClass<GeoField>()
105 typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
110 const GeoField& fld = *iter();
112 typename GeoField::GeometricBoundaryField& bfld =
113 const_cast<typename GeoField::GeometricBoundaryField&>
118 label sz = bfld.size();
123 GeoField::PatchFieldType::New
127 fld.dimensionedInternalField()
134 // Remove last patch field
135 template<class GeoField>
136 void trimPatchFields(fvMesh& mesh, const label nPatches)
138 HashTable<const GeoField*> flds
140 mesh.objectRegistry::lookupClass<GeoField>()
145 typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
150 const GeoField& fld = *iter();
152 const_cast<typename GeoField::GeometricBoundaryField&>
160 // Reorder patch field
161 template<class GeoField>
162 void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
164 HashTable<const GeoField*> flds
166 mesh.objectRegistry::lookupClass<GeoField>()
171 typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
176 const GeoField& fld = *iter();
178 typename GeoField::GeometricBoundaryField& bfld =
179 const_cast<typename GeoField::GeometricBoundaryField&>
184 bfld.reorder(oldToNew);
189 // Adds patch if not yet there. Returns patchID.
190 label addPatch(fvMesh& mesh, const polyPatch& patch)
192 polyBoundaryMesh& polyPatches =
193 const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
195 label patchI = polyPatches.findPatchID(patch.name());
198 if (polyPatches[patchI].type() == patch.type())
205 FatalErrorIn("addPatch(fvMesh&, const polyPatch*)")
206 << "Already have patch " << patch.name()
207 << " but of type " << patch.type()
213 label insertPatchI = polyPatches.size();
214 label startFaceI = mesh.nFaces();
216 forAll(polyPatches, patchI)
218 const polyPatch& pp = polyPatches[patchI];
220 if (isA<processorPolyPatch>(pp))
222 insertPatchI = patchI;
223 startFaceI = pp.start();
229 // Below is all quite a hack. Feel free to change once there is a better
230 // mechanism to insert and reorder patches.
232 // Clear local fields and e.g. polyMesh parallelInfo.
235 label sz = polyPatches.size();
237 fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
239 // Add polyPatch at the end
240 polyPatches.setSize(sz+1);
247 insertPatchI, //index
252 fvPatches.setSize(sz+1);
258 polyPatches[sz], // point to newly added polyPatch
263 addPatchFields<volScalarField>
266 calculatedFvPatchField<scalar>::typeName
268 addPatchFields<volVectorField>
271 calculatedFvPatchField<vector>::typeName
273 addPatchFields<volSphericalTensorField>
276 calculatedFvPatchField<sphericalTensor>::typeName
278 addPatchFields<volSymmTensorField>
281 calculatedFvPatchField<symmTensor>::typeName
283 addPatchFields<volTensorField>
286 calculatedFvPatchField<tensor>::typeName
291 addPatchFields<surfaceScalarField>
294 calculatedFvPatchField<scalar>::typeName
296 addPatchFields<surfaceVectorField>
299 calculatedFvPatchField<vector>::typeName
301 addPatchFields<surfaceSphericalTensorField>
304 calculatedFvPatchField<sphericalTensor>::typeName
306 addPatchFields<surfaceSymmTensorField>
309 calculatedFvPatchField<symmTensor>::typeName
311 addPatchFields<surfaceTensorField>
314 calculatedFvPatchField<tensor>::typeName
317 // Create reordering list
318 // patches before insert position stay as is
319 labelList oldToNew(sz+1);
320 for (label i = 0; i < insertPatchI; i++)
324 // patches after insert position move one up
325 for (label i = insertPatchI; i < sz; i++)
329 // appended patch gets moved to insert position
330 oldToNew[sz] = insertPatchI;
332 // Shuffle into place
333 polyPatches.reorder(oldToNew);
334 fvPatches.reorder(oldToNew);
336 reorderPatchFields<volScalarField>(mesh, oldToNew);
337 reorderPatchFields<volVectorField>(mesh, oldToNew);
338 reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
339 reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
340 reorderPatchFields<volTensorField>(mesh, oldToNew);
341 reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
342 reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
343 reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
344 reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
345 reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
351 // Reorder and delete patches.
355 const labelList& oldToNew,
356 const label nNewPatches
359 polyBoundaryMesh& polyPatches =
360 const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
361 fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
363 // Shuffle into place
364 polyPatches.reorder(oldToNew);
365 fvPatches.reorder(oldToNew);
367 reorderPatchFields<volScalarField>(mesh, oldToNew);
368 reorderPatchFields<volVectorField>(mesh, oldToNew);
369 reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
370 reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
371 reorderPatchFields<volTensorField>(mesh, oldToNew);
372 reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
373 reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
374 reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
375 reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
376 reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
379 polyPatches.setSize(nNewPatches);
380 fvPatches.setSize(nNewPatches);
381 trimPatchFields<volScalarField>(mesh, nNewPatches);
382 trimPatchFields<volVectorField>(mesh, nNewPatches);
383 trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
384 trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
385 trimPatchFields<volTensorField>(mesh, nNewPatches);
386 trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
387 trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
388 trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
389 trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
390 trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
394 template<class GeoField>
398 const fvMesh& subMesh,
399 const labelList& cellMap,
400 const labelList& faceMap,
401 const labelHashSet& addedPatches
404 const labelList patchMap(identity(mesh.boundaryMesh().size()));
406 HashTable<const GeoField*> fields
408 mesh.objectRegistry::lookupClass<GeoField>()
410 forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
412 const GeoField& fld = *iter();
414 Info<< "Mapping field " << fld.name() << endl;
416 tmp<GeoField> tSubFld
418 fvMeshSubset::interpolate
428 // Hack: set value to 0 for introduced patches (since don't
430 forAll(tSubFld().boundaryField(), patchI)
432 if (addedPatches.found(patchI))
434 tSubFld().boundaryField()[patchI] ==
435 pTraits<typename GeoField::value_type>::zero;
440 GeoField* subFld = tSubFld.ptr();
441 subFld->rename(fld.name());
442 subFld->writeOpt() = IOobject::AUTO_WRITE;
448 template<class GeoField>
449 void subsetSurfaceFields
452 const fvMesh& subMesh,
453 const labelList& faceMap,
454 const labelHashSet& addedPatches
457 const labelList patchMap(identity(mesh.boundaryMesh().size()));
459 HashTable<const GeoField*> fields
461 mesh.objectRegistry::lookupClass<GeoField>()
463 forAllConstIter(typename HashTable<const GeoField*>, fields, iter)
465 const GeoField& fld = *iter();
467 Info<< "Mapping field " << fld.name() << endl;
469 tmp<GeoField> tSubFld
471 fvMeshSubset::interpolate
480 // Hack: set value to 0 for introduced patches (since don't
482 forAll(tSubFld().boundaryField(), patchI)
484 if (addedPatches.found(patchI))
486 tSubFld().boundaryField()[patchI] ==
487 pTraits<typename GeoField::value_type>::zero;
492 GeoField* subFld = tSubFld.ptr();
493 subFld->rename(fld.name());
494 subFld->writeOpt() = IOobject::AUTO_WRITE;
499 // Select all cells not in the region
500 labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
502 DynamicList<label> nonRegionCells(cellRegion.size());
503 forAll(cellRegion, cellI)
505 if (cellRegion[cellI] != regionI)
507 nonRegionCells.append(cellI);
510 return nonRegionCells.shrink();
514 // Get per region-region interface the sizes. If sumParallel sums sizes.
515 // Returns interfaces as straight list for looping in identical order.
516 void getInterfaceSizes
518 const polyMesh& mesh,
519 const labelList& cellRegion,
520 const bool sumParallel,
522 edgeList& interfaces,
523 EdgeMap<label>& interfaceSizes
530 forAll(mesh.faceNeighbour(), faceI)
532 label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
533 label neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
535 if (ownRegion != neiRegion)
539 min(ownRegion, neiRegion),
540 max(ownRegion, neiRegion)
543 EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
545 if (iter != interfaceSizes.end())
551 interfaceSizes.insert(interface, 1);
559 // Neighbour cellRegion.
560 labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
562 forAll(coupledRegion, i)
564 label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
565 coupledRegion[i] = cellRegion[cellI];
567 syncTools::swapBoundaryFaceList(mesh, coupledRegion, false);
569 forAll(coupledRegion, i)
571 label faceI = i+mesh.nInternalFaces();
572 label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
573 label neiRegion = coupledRegion[i];
575 if (ownRegion != neiRegion)
579 min(ownRegion, neiRegion),
580 max(ownRegion, neiRegion)
583 EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
585 if (iter != interfaceSizes.end())
591 interfaceSizes.insert(interface, 1);
597 if (sumParallel && Pstream::parRun())
599 if (Pstream::master())
601 // Receive and add to my sizes
604 int slave=Pstream::firstSlave();
605 slave<=Pstream::lastSlave();
609 IPstream fromSlave(Pstream::blocking, slave);
611 EdgeMap<label> slaveSizes(fromSlave);
613 forAllConstIter(EdgeMap<label>, slaveSizes, slaveIter)
615 EdgeMap<label>::iterator masterIter =
616 interfaceSizes.find(slaveIter.key());
618 if (masterIter != interfaceSizes.end())
620 masterIter() += slaveIter();
624 interfaceSizes.insert(slaveIter.key(), slaveIter());
632 int slave=Pstream::firstSlave();
633 slave<=Pstream::lastSlave();
637 // Receive the edges using shared points from the slave.
638 OPstream toSlave(Pstream::blocking, slave);
639 toSlave << interfaceSizes;
646 OPstream toMaster(Pstream::blocking, Pstream::masterNo());
647 toMaster << interfaceSizes;
649 // Receive from master
651 IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
652 fromMaster >> interfaceSizes;
657 // Make sure all processors have interfaces in same order
658 interfaces = interfaceSizes.toc();
661 Pstream::scatter(interfaces);
666 // Create mesh for region.
667 autoPtr<mapPolyMesh> createRegionMesh
669 const labelList& cellRegion,
670 const EdgeMap<label>& interfaceToPatch,
673 const word& regionName,
674 autoPtr<fvMesh>& newMesh
677 // Create dummy system/fv*
682 mesh.time().system(),
690 Info<< "Testing:" << io.objectPath() << endl;
693 // if (!exists(io.objectPath()))
695 Info<< "Writing dummy " << regionName/io.name() << endl;
696 dictionary dummyDict;
698 dummyDict.add("divSchemes", divDict);
700 dummyDict.add("gradSchemes", gradDict);
702 dummyDict.add("laplacianSchemes", laplDict);
704 IOdictionary(io, dummyDict).regIOobject::write();
711 mesh.time().system(),
720 //if (!exists(io.objectPath()))
722 Info<< "Writing dummy " << regionName/io.name() << endl;
723 dictionary dummyDict;
724 IOdictionary(io, dummyDict).regIOobject::write();
729 // Neighbour cellRegion.
730 labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
732 forAll(coupledRegion, i)
734 label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
735 coupledRegion[i] = cellRegion[cellI];
737 syncTools::swapBoundaryFaceList(mesh, coupledRegion, false);
740 // Topology change container. Start off from existing mesh.
741 polyTopoChange meshMod(mesh);
743 // Cell remover engine
744 removeCells cellRemover(mesh);
746 // Select all but region cells
747 labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
749 // Find out which faces will get exposed. Note that this
750 // gets faces in mesh face order. So both regions will get same
751 // face in same order (important!)
752 labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
754 labelList exposedPatchIDs(exposedFaces.size());
755 forAll(exposedFaces, i)
757 label faceI = exposedFaces[i];
759 label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
760 label neiRegion = -1;
762 if (mesh.isInternalFace(faceI))
764 neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
768 neiRegion = coupledRegion[faceI-mesh.nInternalFaces()];
771 label otherRegion = -1;
773 if (ownRegion == regionI && neiRegion != regionI)
775 otherRegion = neiRegion;
777 else if (ownRegion != regionI && neiRegion == regionI)
779 otherRegion = ownRegion;
783 FatalErrorIn("createRegionMesh(..)")
784 << "Exposed face:" << faceI
785 << " fc:" << mesh.faceCentres()[faceI]
786 << " has owner region " << ownRegion
787 << " and neighbour region " << neiRegion
788 << " when handling region:" << regionI
792 if (otherRegion != -1)
794 edge interface(regionI, otherRegion);
797 if (regionI < otherRegion)
799 exposedPatchIDs[i] = interfaceToPatch[interface];
803 exposedPatchIDs[i] = interfaceToPatch[interface]+1;
809 cellRemover.setRefinement
817 autoPtr<mapPolyMesh> map = meshMod.makeMesh
823 mesh.time().timeName(),
835 void createAndWriteRegion
838 const labelList& cellRegion,
839 const wordList& regionNames,
840 const EdgeMap<label>& interfaceToPatch,
842 const word& newMeshInstance
845 Info<< "Creating mesh for region " << regionI
846 << ' ' << regionNames[regionI] << endl;
848 autoPtr<fvMesh> newMesh;
849 autoPtr<mapPolyMesh> map = createRegionMesh
855 regionNames[regionI],
860 // Make map of all added patches
861 labelHashSet addedPatches(2*interfaceToPatch.size());
862 forAllConstIter(EdgeMap<label>, interfaceToPatch, iter)
864 addedPatches.insert(iter());
865 addedPatches.insert(iter()+1);
868 Info<< "Mapping fields" << endl;
870 // Map existing fields
871 newMesh().updateMesh(map());
873 // Add subsetted fields
874 subsetVolFields<volScalarField>
882 subsetVolFields<volVectorField>
890 subsetVolFields<volSphericalTensorField>
898 subsetVolFields<volSymmTensorField>
906 subsetVolFields<volTensorField>
915 subsetSurfaceFields<surfaceScalarField>
922 subsetSurfaceFields<surfaceVectorField>
929 subsetSurfaceFields<surfaceSphericalTensorField>
936 subsetSurfaceFields<surfaceSymmTensorField>
943 subsetSurfaceFields<surfaceTensorField>
952 const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
953 newPatches.checkParallelSync(true);
955 // Delete empty patches
956 // ~~~~~~~~~~~~~~~~~~~~
958 // Create reordering list to move patches-to-be-deleted to end
959 labelList oldToNew(newPatches.size(), -1);
962 Info<< "Deleting empty patches" << endl;
964 // Assumes all non-proc boundaries are on all processors!
965 forAll(newPatches, patchI)
967 const polyPatch& pp = newPatches[patchI];
969 if (!isA<processorPolyPatch>(pp))
971 if (returnReduce(pp.size(), sumOp<label>()) > 0)
973 oldToNew[patchI] = newI++;
978 // Same for processor patches (but need no reduction)
979 forAll(newPatches, patchI)
981 const polyPatch& pp = newPatches[patchI];
983 if (isA<processorPolyPatch>(pp) && pp.size())
985 oldToNew[patchI] = newI++;
989 const label nNewPatches = newI;
991 // Move all deleteable patches to the end
992 forAll(oldToNew, patchI)
994 if (oldToNew[patchI] == -1)
996 oldToNew[patchI] = newI++;
1000 reorderPatches(newMesh(), oldToNew, nNewPatches);
1003 Info<< "Writing new mesh" << endl;
1005 newMesh().setInstance(newMeshInstance);
1008 // Write addressing files like decomposePar
1009 Info<< "Writing addressing to base mesh" << endl;
1011 labelIOList pointProcAddressing
1015 "pointRegionAddressing",
1016 newMesh().facesInstance(),
1017 newMesh().meshSubDir,
1025 Info<< "Writing map " << pointProcAddressing.name()
1026 << " from region" << regionI
1027 << " points back to base mesh." << endl;
1028 pointProcAddressing.write();
1030 labelIOList faceProcAddressing
1034 "faceRegionAddressing",
1035 newMesh().facesInstance(),
1036 newMesh().meshSubDir,
1044 forAll(faceProcAddressing, faceI)
1046 // face + turning index. (see decomposePar)
1047 // Is the face pointing in the same direction?
1048 label oldFaceI = map().faceMap()[faceI];
1052 map().cellMap()[newMesh().faceOwner()[faceI]]
1053 == mesh.faceOwner()[oldFaceI]
1056 faceProcAddressing[faceI] = oldFaceI+1;
1060 faceProcAddressing[faceI] = -(oldFaceI+1);
1063 Info<< "Writing map " << faceProcAddressing.name()
1064 << " from region" << regionI
1065 << " faces back to base mesh." << endl;
1066 faceProcAddressing.write();
1068 labelIOList cellProcAddressing
1072 "cellRegionAddressing",
1073 newMesh().facesInstance(),
1074 newMesh().meshSubDir,
1082 Info<< "Writing map " <<cellProcAddressing.name()
1083 << " from region" << regionI
1084 << " cells back to base mesh." << endl;
1085 cellProcAddressing.write();
1089 // Create for every region-region interface with non-zero size two patches.
1090 // First one is for minimumregion to maximumregion.
1091 // Note that patches get created in same order on all processors (if parallel)
1092 // since looping over synchronised 'interfaces'.
1093 EdgeMap<label> addRegionPatches
1096 const labelList& cellRegion,
1097 const label nCellRegions,
1098 const edgeList& interfaces,
1099 const EdgeMap<label>& interfaceSizes,
1100 const wordList& regionNames
1103 // Check that all patches are present in same order.
1104 mesh.boundaryMesh().checkParallelSync(true);
1106 Info<< nl << "Adding patches" << nl << endl;
1108 EdgeMap<label> interfaceToPatch(nCellRegions);
1110 forAll(interfaces, interI)
1112 const edge& e = interfaces[interI];
1114 if (interfaceSizes[e] > 0)
1116 const word inter1 = regionNames[e[0]] + "_to_" + regionNames[e[1]];
1117 const word inter2 = regionNames[e[1]] + "_to_" + regionNames[e[0]];
1119 directMappedWallPolyPatch patch1
1125 regionNames[e[1]], // sampleRegion
1126 directMappedPatchBase::NEARESTPATCHFACE,
1127 inter2, // samplePatch
1128 point::zero, // offset
1132 label patchI = addPatch(mesh, patch1);
1134 directMappedWallPolyPatch patch2
1140 regionNames[e[0]], // sampleRegion
1141 directMappedPatchBase::NEARESTPATCHFACE,
1143 point::zero, // offset
1146 addPatch(mesh, patch2);
1148 Info<< "For interface between region " << e[0]
1149 << " and " << e[1] << " added patch " << patchI
1150 << " " << mesh.boundaryMesh()[patchI].name()
1153 interfaceToPatch.insert(e, patchI);
1156 return interfaceToPatch;
1160 // Find region that covers most of cell zone
1161 label findCorrespondingRegion
1163 const labelList& existingZoneID, // per cell the (unique) zoneID
1164 const labelList& cellRegion,
1165 const label nCellRegions,
1167 const label minOverlapSize
1170 // Per region the number of cells in zoneI
1171 labelList cellsInZone(nCellRegions, 0);
1173 forAll(cellRegion, cellI)
1175 if (existingZoneID[cellI] == zoneI)
1177 cellsInZone[cellRegion[cellI]]++;
1181 Pstream::listCombineGather(cellsInZone, plusEqOp<label>());
1182 Pstream::listCombineScatter(cellsInZone);
1184 // Pick region with largest overlap of zoneI
1185 label regionI = findMax(cellsInZone);
1188 if (cellsInZone[regionI] < minOverlapSize)
1190 // Region covers too little of zone. Not good enough.
1195 // Check that region contains no cells that aren't in cellZone.
1196 forAll(cellRegion, cellI)
1198 if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1200 // cellI in regionI but not in zoneI
1205 // If one in error, all should be in error. Note that branch gets taken
1207 reduce(regionI, minOp<label>());
1214 //// Checks if cellZone has corresponding cellRegion.
1215 //label findCorrespondingRegion
1217 // const cellZoneMesh& cellZones,
1218 // const labelList& existingZoneID, // per cell the (unique) zoneID
1219 // const labelList& cellRegion,
1220 // const label nCellRegions,
1221 // const label zoneI
1224 // // Region corresponding to zone. Start off with special value: no
1225 // // corresponding region.
1226 // label regionI = labelMax;
1228 // const cellZone& cz = cellZones[zoneI];
1232 // // My local portion is empty. Maps to any empty cellZone. Mark with
1233 // // special value which can get overwritten by other processors.
1238 // regionI = cellRegion[cz[0]];
1242 // if (cellRegion[cz[i]] != regionI)
1244 // regionI = labelMax;
1250 // // Determine same zone over all processors.
1251 // reduce(regionI, maxOp<label>());
1254 // // 2. All of region present?
1256 // if (regionI == labelMax)
1260 // else if (regionI != -1)
1262 // forAll(cellRegion, cellI)
1264 // if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1266 // // cellI in regionI but not in zoneI
1271 // // If one in error, all should be in error. Note that branch gets taken
1273 // reduce(regionI, minOp<label>());
1280 // Get zone per cell
1281 // - non-unique zoning
1285 const polyMesh& mesh,
1286 const cellZoneMesh& cellZones,
1288 labelList& neiZoneID
1292 zoneID.setSize(mesh.nCells());
1295 forAll(cellZones, zoneI)
1297 const cellZone& cz = cellZones[zoneI];
1301 label cellI = cz[i];
1302 if (zoneID[cellI] == -1)
1304 zoneID[cellI] = zoneI;
1308 FatalErrorIn("getZoneID(..)")
1309 << "Cell " << cellI << " with cell centre "
1310 << mesh.cellCentres()[cellI]
1311 << " is multiple zones. This is not allowed." << endl
1312 << "It is in zone " << cellZones[zoneID[cellI]].name()
1313 << " and in zone " << cellZones[zoneI].name()
1314 << exit(FatalError);
1319 // Neighbour zoneID.
1320 neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
1322 forAll(neiZoneID, i)
1324 neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
1326 syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
1332 const bool sloppyCellZones,
1333 const polyMesh& mesh,
1335 const label nCellRegions,
1336 const labelList& cellRegion,
1338 labelList& regionToZone,
1339 wordList& regionNames,
1340 labelList& zoneToRegion
1343 const cellZoneMesh& cellZones = mesh.cellZones();
1345 regionToZone.setSize(nCellRegions, -1);
1346 regionNames.setSize(nCellRegions);
1347 zoneToRegion.setSize(cellZones.size(), -1);
1349 // Get current per cell zoneID
1350 labelList zoneID(mesh.nCells(), -1);
1351 labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1352 getZoneID(mesh, cellZones, zoneID, neiZoneID);
1354 // Sizes per cellzone
1355 labelList zoneSizes(cellZones.size(), 0);
1357 List<wordList> zoneNames(Pstream::nProcs());
1358 zoneNames[Pstream::myProcNo()] = cellZones.names();
1359 Pstream::gatherList(zoneNames);
1360 Pstream::scatterList(zoneNames);
1362 forAll(zoneNames, procI)
1364 if (zoneNames[procI] != zoneNames[0])
1366 FatalErrorIn("matchRegions(..)")
1367 << "cellZones not synchronised across processors." << endl
1368 << "Master has cellZones " << zoneNames[0] << endl
1369 << "Processor " << procI
1370 << " has cellZones " << zoneNames[procI]
1371 << exit(FatalError);
1375 forAll(cellZones, zoneI)
1377 zoneSizes[zoneI] = returnReduce
1379 cellZones[zoneI].size(),
1386 if (sloppyCellZones)
1388 Info<< "Trying to match regions to existing cell zones;"
1389 << " region can be subset of cell zone." << nl << endl;
1391 forAll(cellZones, zoneI)
1393 label regionI = findCorrespondingRegion
1399 label(0.5*zoneSizes[zoneI]) // minimum overlap
1404 Info<< "Sloppily matched region " << regionI
1405 //<< " size " << regionSizes[regionI]
1406 << " to zone " << zoneI << " size " << zoneSizes[zoneI]
1408 zoneToRegion[zoneI] = regionI;
1409 regionToZone[regionI] = zoneI;
1410 regionNames[regionI] = cellZones[zoneI].name();
1416 Info<< "Trying to match regions to existing cell zones." << nl << endl;
1418 forAll(cellZones, zoneI)
1420 label regionI = findCorrespondingRegion
1426 1 // minimum overlap
1431 zoneToRegion[zoneI] = regionI;
1432 regionToZone[regionI] = zoneI;
1433 regionNames[regionI] = cellZones[zoneI].name();
1437 // Allocate region names for unmatched regions.
1438 forAll(regionToZone, regionI)
1440 if (regionToZone[regionI] == -1)
1442 regionNames[regionI] = "domain" + Foam::name(regionI);
1448 void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
1450 // Write to manual decomposition option
1452 labelIOList cellToRegion
1457 mesh.facesInstance(),
1465 cellToRegion.write();
1467 Info<< "Writing region per cell file (for manual decomposition) to "
1468 << cellToRegion.objectPath() << nl << endl;
1470 // Write for postprocessing
1472 volScalarField cellToRegion
1477 mesh.time().timeName(),
1484 dimensionedScalar("zero", dimless, 0),
1485 zeroGradientFvPatchScalarField::typeName
1487 forAll(cellRegion, cellI)
1489 cellToRegion[cellI] = cellRegion[cellI];
1491 cellToRegion.write();
1493 Info<< "Writing region per cell as volScalarField to "
1494 << cellToRegion.objectPath() << nl << endl;
1502 int main(int argc, char *argv[])
1504 argList::validOptions.insert("cellZones", "");
1505 argList::validOptions.insert("cellZonesOnly", "");
1506 argList::validOptions.insert("cellZonesFileOnly", "cellZonesName");
1507 argList::validOptions.insert("blockedFaces", "faceSet");
1508 argList::validOptions.insert("makeCellZones", "");
1509 argList::validOptions.insert("largestOnly", "");
1510 argList::validOptions.insert("insidePoint", "point");
1511 argList::validOptions.insert("overwrite", "");
1512 argList::validOptions.insert("detectOnly", "");
1513 argList::validOptions.insert("sloppyCellZones", "");
1515 # include "setRootCase.H"
1516 # include "createTime.H"
1517 runTime.functionObjects().off();
1518 # include "createMesh.H"
1519 const word oldInstance = mesh.pointsInstance();
1521 word blockedFacesName;
1522 if (args.optionFound("blockedFaces"))
1524 blockedFacesName = args.option("blockedFaces");
1525 Info<< "Reading blocked internal faces from faceSet "
1526 << blockedFacesName << nl << endl;
1529 bool makeCellZones = args.optionFound("makeCellZones");
1530 bool largestOnly = args.optionFound("largestOnly");
1531 bool insidePoint = args.optionFound("insidePoint");
1532 bool useCellZones = args.optionFound("cellZones");
1533 bool useCellZonesOnly = args.optionFound("cellZonesOnly");
1534 bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
1535 bool overwrite = args.optionFound("overwrite");
1536 bool detectOnly = args.optionFound("detectOnly");
1537 bool sloppyCellZones = args.optionFound("sloppyCellZones");
1541 (useCellZonesOnly || useCellZonesFile)
1543 blockedFacesName != word::null
1548 FatalErrorIn(args.executable())
1549 << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
1550 << " (which specify complete split)"
1551 << " in combination with -blockedFaces or -cellZones"
1552 << " (which imply a split based on topology)"
1553 << exit(FatalError);
1558 if (insidePoint && largestOnly)
1560 FatalErrorIn(args.executable())
1561 << "You cannot specify both -largestOnly"
1562 << " (keep region with most cells)"
1563 << " and -insidePoint (keep region containing point)"
1564 << exit(FatalError);
1568 const cellZoneMesh& cellZones = mesh.cellZones();
1571 labelList zoneID(mesh.nCells(), -1);
1572 // Neighbour zoneID.
1573 labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1574 getZoneID(mesh, cellZones, zoneID, neiZoneID);
1578 // Determine per cell the region it belongs to
1579 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1581 // cellRegion is the labelList with the region per cell.
1582 labelList cellRegion;
1584 labelList regionToZone;
1586 wordList regionNames;
1588 labelList zoneToRegion;
1590 label nCellRegions = 0;
1591 if (useCellZonesOnly)
1593 Info<< "Using current cellZones to split mesh into regions."
1594 << " This requires all"
1595 << " cells to be in one and only one cellZone." << nl << endl;
1597 label unzonedCellI = findIndex(zoneID, -1);
1598 if (unzonedCellI != -1)
1600 FatalErrorIn(args.executable())
1601 << "For the cellZonesOnly option all cells "
1602 << "have to be in a cellZone." << endl
1603 << "Cell " << unzonedCellI
1604 << " at" << mesh.cellCentres()[unzonedCellI]
1605 << " is not in a cellZone. There might be more unzoned cells."
1606 << exit(FatalError);
1608 cellRegion = zoneID;
1609 nCellRegions = gMax(cellRegion)+1;
1610 regionToZone.setSize(nCellRegions);
1611 regionNames.setSize(nCellRegions);
1612 zoneToRegion.setSize(cellZones.size(), -1);
1613 for (label regionI = 0; regionI < nCellRegions; regionI++)
1615 regionToZone[regionI] = regionI;
1616 zoneToRegion[regionI] = regionI;
1617 regionNames[regionI] = cellZones[regionI].name();
1620 else if (useCellZonesFile)
1622 const word zoneFile = args.option("cellZonesFileOnly");
1623 Info<< "Reading split from cellZones file " << zoneFile << endl
1624 << "This requires all"
1625 << " cells to be in one and only one cellZone." << nl << endl;
1627 cellZoneMesh newCellZones
1632 mesh.facesInstance(),
1633 polyMesh::meshSubDir,
1635 IOobject::MUST_READ,
1642 labelList newZoneID(mesh.nCells(), -1);
1643 labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1644 getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
1646 label unzonedCellI = findIndex(newZoneID, -1);
1647 if (unzonedCellI != -1)
1649 FatalErrorIn(args.executable())
1650 << "For the cellZonesFileOnly option all cells "
1651 << "have to be in a cellZone." << endl
1652 << "Cell " << unzonedCellI
1653 << " at" << mesh.cellCentres()[unzonedCellI]
1654 << " is not in a cellZone. There might be more unzoned cells."
1655 << exit(FatalError);
1657 cellRegion = newZoneID;
1658 nCellRegions = gMax(cellRegion)+1;
1659 zoneToRegion.setSize(newCellZones.size(), -1);
1660 regionToZone.setSize(nCellRegions);
1661 regionNames.setSize(nCellRegions);
1662 for (label regionI = 0; regionI < nCellRegions; regionI++)
1664 regionToZone[regionI] = regionI;
1665 zoneToRegion[regionI] = regionI;
1666 regionNames[regionI] = newCellZones[regionI].name();
1671 // Determine connected regions
1672 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1674 // Mark additional faces that are blocked
1675 boolList blockedFace;
1677 // Read from faceSet
1678 if (blockedFacesName.size())
1680 faceSet blockedFaceSet(mesh, blockedFacesName);
1682 << returnReduce(blockedFaceSet.size(), sumOp<label>())
1683 << " blocked faces from set " << blockedFacesName << nl << endl;
1685 blockedFace.setSize(mesh.nFaces(), false);
1687 forAllConstIter(faceSet, blockedFaceSet, iter)
1689 blockedFace[iter.key()] = true;
1693 // Imply from differing cellZones
1696 blockedFace.setSize(mesh.nFaces(), false);
1698 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
1700 label own = mesh.faceOwner()[faceI];
1701 label nei = mesh.faceNeighbour()[faceI];
1703 if (zoneID[own] != zoneID[nei])
1705 blockedFace[faceI] = true;
1709 // Different cellZones on either side of processor patch.
1710 forAll(neiZoneID, i)
1712 label faceI = i+mesh.nInternalFaces();
1714 if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
1716 blockedFace[faceI] = true;
1721 // Do a topological walk to determine regions
1722 regionSplit regions(mesh, blockedFace);
1723 nCellRegions = regions.nRegions();
1724 cellRegion.transfer(regions);
1726 // Make up region names. If possible match them to existing zones.
1740 Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
1743 // Write decomposition to file
1744 writeCellToRegion(mesh, cellRegion);
1751 labelList regionSizes(nCellRegions, 0);
1753 forAll(cellRegion, cellI)
1755 regionSizes[cellRegion[cellI]]++;
1757 forAll(regionSizes, regionI)
1759 reduce(regionSizes[regionI], sumOp<label>());
1762 Info<< "Region\tCells" << nl
1763 << "------\t-----" << endl;
1765 forAll(regionSizes, regionI)
1767 Info<< regionI << '\t' << regionSizes[regionI] << nl;
1773 // Print region to zone
1774 Info<< "Region\tZone\tName" << nl
1775 << "------\t----\t----" << endl;
1776 forAll(regionToZone, regionI)
1778 Info<< regionI << '\t' << regionToZone[regionI] << '\t'
1779 << regionNames[regionI] << nl;
1785 // Since we're going to mess with patches make sure all non-processor ones
1786 // are on all processors.
1787 mesh.boundaryMesh().checkParallelSync(true);
1790 // Sizes of interface between regions. From pair of regions to number of
1792 edgeList interfaces;
1793 EdgeMap<label> interfaceSizes;
1798 true, // sum in parallel?
1804 Info<< "Sizes inbetween regions:" << nl << nl
1805 << "Region\tRegion\tFaces" << nl
1806 << "------\t------\t-----" << endl;
1808 forAll(interfaces, interI)
1810 const edge& e = interfaces[interI];
1812 Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
1823 // Read objects in time directory
1824 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1826 IOobjectList objects(mesh, runTime.timeName());
1830 PtrList<volScalarField> vsFlds;
1831 ReadFields(mesh, objects, vsFlds);
1833 PtrList<volVectorField> vvFlds;
1834 ReadFields(mesh, objects, vvFlds);
1836 PtrList<volSphericalTensorField> vstFlds;
1837 ReadFields(mesh, objects, vstFlds);
1839 PtrList<volSymmTensorField> vsymtFlds;
1840 ReadFields(mesh, objects, vsymtFlds);
1842 PtrList<volTensorField> vtFlds;
1843 ReadFields(mesh, objects, vtFlds);
1845 // Read surface fields.
1847 PtrList<surfaceScalarField> ssFlds;
1848 ReadFields(mesh, objects, ssFlds);
1850 PtrList<surfaceVectorField> svFlds;
1851 ReadFields(mesh, objects, svFlds);
1853 PtrList<surfaceSphericalTensorField> sstFlds;
1854 ReadFields(mesh, objects, sstFlds);
1856 PtrList<surfaceSymmTensorField> ssymtFlds;
1857 ReadFields(mesh, objects, ssymtFlds);
1859 PtrList<surfaceTensorField> stFlds;
1860 ReadFields(mesh, objects, stFlds);
1865 // Remove any demand-driven fields ('S', 'V' etc)
1869 if (nCellRegions == 1)
1871 Info<< "Only one region. Doing nothing." << endl;
1873 else if (makeCellZones)
1875 Info<< "Putting cells into cellZones instead of splitting mesh."
1878 // Check if region overlaps with existing zone. If so keep.
1880 for (label regionI = 0; regionI < nCellRegions; regionI++)
1882 label zoneI = regionToZone[regionI];
1886 Info<< " Region " << regionI << " : corresponds to existing"
1888 << zoneI << ' ' << cellZones[zoneI].name() << endl;
1892 // Create new cellZone.
1893 labelList regionCells = findIndices(cellRegion, regionI);
1895 word zoneName = "region" + Foam::name(regionI);
1897 zoneI = cellZones.findZoneID(zoneName);
1901 zoneI = cellZones.size();
1902 mesh.cellZones().setSize(zoneI+1);
1903 mesh.cellZones().set
1909 regionCells, //addressing
1911 cellZones //cellZoneMesh
1917 mesh.cellZones()[zoneI].clearAddressing();
1918 mesh.cellZones()[zoneI] = regionCells;
1920 Info<< " Region " << regionI << " : created new cellZone "
1921 << zoneI << ' ' << cellZones[zoneI].name() << endl;
1924 mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
1929 mesh.setInstance(runTime.timeName());
1933 mesh.setInstance(oldInstance);
1936 Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
1942 // Write cellSets for convenience
1943 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1945 Info<< "Writing cellSets corresponding to cellZones." << nl << endl;
1947 forAll(cellZones, zoneI)
1949 const cellZone& cz = cellZones[zoneI];
1951 cellSet(mesh, cz.name(), cz).write();
1956 // Add patches for interfaces
1957 // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1959 // Add all possible patches. Empty ones get filtered later on.
1960 Info<< nl << "Adding patches" << nl << endl;
1962 EdgeMap<label> interfaceToPatch
1987 point insidePoint(args.optionLookup("insidePoint")());
1991 label cellI = mesh.findCell(insidePoint);
1993 Info<< nl << "Found point " << insidePoint << " in cell " << cellI
1998 regionI = cellRegion[cellI];
2001 reduce(regionI, maxOp<label>());
2004 << "Subsetting region " << regionI
2005 << " containing point " << insidePoint << endl;
2009 FatalErrorIn(args.executable())
2010 << "Point " << insidePoint
2011 << " is not inside the mesh." << nl
2012 << "Bounding box of the mesh:" << mesh.bounds()
2013 << exit(FatalError);
2016 createAndWriteRegion
2023 (overwrite ? oldInstance : runTime.timeName())
2026 else if (largestOnly)
2028 label regionI = findMax(regionSizes);
2031 << "Subsetting region " << regionI
2032 << " of size " << regionSizes[regionI] << endl;
2034 createAndWriteRegion
2041 (overwrite ? oldInstance : runTime.timeName())
2047 for (label regionI = 0; regionI < nCellRegions; regionI++)
2050 << "Region " << regionI << nl
2051 << "-------- " << endl;
2053 createAndWriteRegion
2060 (overwrite ? oldInstance : runTime.timeName())
2066 Info<< "End\n" << endl;
2072 // ************************************************************************* //