FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
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
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  splitMeshRegions
26 
27 Description
28  Splits mesh into multiple regions.
29 
30  Each region is defined as a domain whose cells can all be reached by
31  cell-face-cell walking without crossing
32  - boundary faces
33  - additional faces from faceset (-blockedFaces faceSet).
34  - any face inbetween differing cellZones (-cellZones)
35 
36  Output is:
37  - volScalarField with regions as different scalars (-detectOnly) or
38  - mesh with multiple regions or
39  - mesh with cells put into cellZones (-makeCellZones)
40 
41 Note
42  - cellZonesOnly does not do a walk and uses the cellZones only. Use
43  this if you don't mind having disconnected domains in a single region.
44  This option requires all cells to be in one (and one only) cellZone.
45 
46  - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
47  from the specified file. This allows one to explicitly specify the region
48  distribution and still have multiple cellZones per region.
49 
50  - useCellZonesOnly does not do a walk and uses the cellZones only. Use
51  this if you don't mind having disconnected domains in a single region.
52  This option requires all cells to be in one (and one only) cellZone.
53 
54 
55  - Should work in parallel.
56  cellZones can differ on either side of processor boundaries in which case
57  the faces get moved from processor patch to directMapped patch. Not
58  very well tested.
59 
60  - If a cell zone gets split into more than one region it can detect
61  the largest matching region (-sloppyCellZones). This will accept any
62  region that covers more than 50% of the zone. It has to be a subset
63  so cannot have any cells in any other zone.
64 
65  - writes maps like decomposePar back to original mesh:
66  - pointRegionAddressing : for every point in this region the point in
67  the original mesh
68  - cellRegionAddressing : ,, cell ,, cell ,,
69  - faceRegionAddressing : ,, face ,, face in
70  the original mesh + 'turning index'. For a face in the same orientation
71  this is the original facelabel+1, for a turned face this is -facelabel-1
72 
73 Usage
74 
75  - splitMeshRegions [OPTIONS]
76 
77  @param -cellZones \n
78  Split different cell zones.
79 
80  @param -detectOnly \n
81  Do no processing.
82 
83  @param -blockedFaces <faceSet name>\n
84  Split at give face set.
85 
86  @param -sloppyCellZones \n
87  Try to match regions to existing cell zones.
88 
89  @param -makeCellZones \n
90  Create mesh with cells in different cell zones.
91 
92  @param -overwrite \n
93  Overwrite existing data.
94 
95  @param -case <dir>\n
96  Case directory.
97 
98  @param -parallel \n
99  Run in parallel.
100 
101  @param -help \n
102  Display help message.
103 
104  @param -doc \n
105  Display Doxygen API documentation page for this application.
106 
107  @param -srcDoc \n
108  Display Doxygen source documentation page for this application.
109 
110  @param -insidePoint <point>\n
111  Only keep region containing specified point.
112 
113  @param -largestOnly \n
114  Only keep largest region.
115 
116 \*---------------------------------------------------------------------------*/
117 
118 #include <OpenFOAM/SortableList.H>
119 #include <OpenFOAM/argList.H>
120 #include <meshTools/regionSplit.H>
122 #include <OpenFOAM/IOobjectList.H>
123 #include <finiteVolume/volFields.H>
124 #include <meshTools/faceSet.H>
125 #include <meshTools/cellSet.H>
127 #include <dynamicMesh/removeCells.H>
128 #include <OpenFOAM/EdgeMap.H>
129 #include <OpenFOAM/syncTools.H>
130 #include <OpenFOAM/ReadFields.H>
133 
134 using namespace Foam;
135 
136 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
137 
138 template<class GeoField>
139 void addPatchFields(fvMesh& mesh, const word& patchFieldType)
140 {
142  (
143  mesh.objectRegistry::lookupClass<GeoField>()
144  );
145 
146  for
147  (
148  typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
149  iter != flds.end();
150  ++iter
151  )
152  {
153  const GeoField& fld = *iter();
154 
155  typename GeoField::GeometricBoundaryField& bfld =
156  const_cast<typename GeoField::GeometricBoundaryField&>
157  (
158  fld.boundaryField()
159  );
160 
161  label sz = bfld.size();
162  bfld.setSize(sz+1);
163  bfld.set
164  (
165  sz,
166  GeoField::PatchFieldType::New
167  (
168  patchFieldType,
169  mesh.boundary()[sz],
170  fld.dimensionedInternalField()
171  )
172  );
173  }
174 }
175 
176 
177 // Remove last patch field
178 template<class GeoField>
179 void trimPatchFields(fvMesh& mesh, const label nPatches)
180 {
182  (
183  mesh.objectRegistry::lookupClass<GeoField>()
184  );
185 
186  for
187  (
188  typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
189  iter != flds.end();
190  ++iter
191  )
192  {
193  const GeoField& fld = *iter();
194 
195  const_cast<typename GeoField::GeometricBoundaryField&>
196  (
197  fld.boundaryField()
198  ).setSize(nPatches);
199  }
200 }
201 
202 
203 // Reorder patch field
204 template<class GeoField>
205 void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
206 {
208  (
209  mesh.objectRegistry::lookupClass<GeoField>()
210  );
211 
212  for
213  (
214  typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
215  iter != flds.end();
216  ++iter
217  )
218  {
219  const GeoField& fld = *iter();
220 
221  typename GeoField::GeometricBoundaryField& bfld =
222  const_cast<typename GeoField::GeometricBoundaryField&>
223  (
224  fld.boundaryField()
225  );
226 
227  bfld.reorder(oldToNew);
228  }
229 }
230 
231 
232 // Adds patch if not yet there. Returns patchID.
233 label addPatch(fvMesh& mesh, const polyPatch& patch)
234 {
235  polyBoundaryMesh& polyPatches =
236  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
237 
238  label patchI = polyPatches.findPatchID(patch.name());
239  if (patchI != -1)
240  {
241  if (polyPatches[patchI].type() == patch.type())
242  {
243  // Already there
244  return patchI;
245  }
246  else
247  {
248  FatalErrorIn("addPatch(fvMesh&, const polyPatch*)")
249  << "Already have patch " << patch.name()
250  << " but of type " << patch.type()
251  << exit(FatalError);
252  }
253  }
254 
255 
256  label insertPatchI = polyPatches.size();
257  label startFaceI = mesh.nFaces();
258 
259  forAll(polyPatches, patchI)
260  {
261  const polyPatch& pp = polyPatches[patchI];
262 
263  if (isA<processorPolyPatch>(pp))
264  {
265  insertPatchI = patchI;
266  startFaceI = pp.start();
267  break;
268  }
269  }
270 
271 
272  // Below is all quite a hack. Feel free to change once there is a better
273  // mechanism to insert and reorder patches.
274 
275  // Clear local fields and e.g. polyMesh parallelInfo.
276  mesh.clearOut();
277 
278  label sz = polyPatches.size();
279 
280  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
281 
282  // Add polyPatch at the end
283  polyPatches.setSize(sz+1);
284  polyPatches.set
285  (
286  sz,
287  patch.clone
288  (
289  polyPatches,
290  insertPatchI, //index
291  0, //size
292  startFaceI //start
293  )
294  );
295  fvPatches.setSize(sz+1);
296  fvPatches.set
297  (
298  sz,
300  (
301  polyPatches[sz], // point to newly added polyPatch
302  mesh.boundary()
303  )
304  );
305 
306  addPatchFields<volScalarField>
307  (
308  mesh,
310  );
311  addPatchFields<volVectorField>
312  (
313  mesh,
315  );
316  addPatchFields<volSphericalTensorField>
317  (
318  mesh,
320  );
321  addPatchFields<volSymmTensorField>
322  (
323  mesh,
325  );
326  addPatchFields<volTensorField>
327  (
328  mesh,
330  );
331 
332  // Surface fields
333 
334  addPatchFields<surfaceScalarField>
335  (
336  mesh,
338  );
339  addPatchFields<surfaceVectorField>
340  (
341  mesh,
343  );
344  addPatchFields<surfaceSphericalTensorField>
345  (
346  mesh,
348  );
349  addPatchFields<surfaceSymmTensorField>
350  (
351  mesh,
353  );
354  addPatchFields<surfaceTensorField>
355  (
356  mesh,
358  );
359 
360  // Create reordering list
361  // patches before insert position stay as is
362  labelList oldToNew(sz+1);
363  for (label i = 0; i < insertPatchI; i++)
364  {
365  oldToNew[i] = i;
366  }
367  // patches after insert position move one up
368  for (label i = insertPatchI; i < sz; i++)
369  {
370  oldToNew[i] = i+1;
371  }
372  // appended patch gets moved to insert position
373  oldToNew[sz] = insertPatchI;
374 
375  // Shuffle into place
376  polyPatches.reorder(oldToNew);
377  fvPatches.reorder(oldToNew);
378 
379  reorderPatchFields<volScalarField>(mesh, oldToNew);
380  reorderPatchFields<volVectorField>(mesh, oldToNew);
381  reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
382  reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
383  reorderPatchFields<volTensorField>(mesh, oldToNew);
384  reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
385  reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
386  reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
387  reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
388  reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
389 
390  return insertPatchI;
391 }
392 
393 
394 // Reorder and delete patches.
395 void reorderPatches
396 (
397  fvMesh& mesh,
398  const labelList& oldToNew,
399  const label nNewPatches
400 )
401 {
402  polyBoundaryMesh& polyPatches =
403  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
404  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
405 
406  // Shuffle into place
407  polyPatches.reorder(oldToNew);
408  fvPatches.reorder(oldToNew);
409 
410  reorderPatchFields<volScalarField>(mesh, oldToNew);
411  reorderPatchFields<volVectorField>(mesh, oldToNew);
412  reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
413  reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
414  reorderPatchFields<volTensorField>(mesh, oldToNew);
415  reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
416  reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
417  reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
418  reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
419  reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
420 
421  // Remove last.
422  polyPatches.setSize(nNewPatches);
423  fvPatches.setSize(nNewPatches);
424  trimPatchFields<volScalarField>(mesh, nNewPatches);
425  trimPatchFields<volVectorField>(mesh, nNewPatches);
426  trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
427  trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
428  trimPatchFields<volTensorField>(mesh, nNewPatches);
429  trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
430  trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
431  trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
432  trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
433  trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
434 }
435 
436 
437 template<class GeoField>
438 void subsetVolFields
439 (
440  const fvMesh& mesh,
441  const fvMesh& subMesh,
442  const labelList& cellMap,
443  const labelList& faceMap,
444  const labelHashSet& addedPatches
445 )
446 {
447  const labelList patchMap(identity(mesh.boundaryMesh().size()));
448 
450  (
451  mesh.objectRegistry::lookupClass<GeoField>()
452  );
454  {
455  const GeoField& fld = *iter();
456 
457  Info<< "Mapping field " << fld.name() << endl;
458 
459  tmp<GeoField> tSubFld
460  (
462  (
463  fld,
464  subMesh,
465  patchMap,
466  cellMap,
467  faceMap
468  )
469  );
470 
471  // Hack: set value to 0 for introduced patches (since don't
472  // get initialised.
473  forAll(tSubFld().boundaryField(), patchI)
474  {
475  if (addedPatches.found(patchI))
476  {
477  tSubFld().boundaryField()[patchI] ==
479  }
480  }
481 
482  // Store on subMesh
483  GeoField* subFld = tSubFld.ptr();
484  subFld->rename(fld.name());
485  subFld->writeOpt() = IOobject::AUTO_WRITE;
486  subFld->store();
487  }
488 }
489 
490 
491 template<class GeoField>
492 void subsetSurfaceFields
493 (
494  const fvMesh& mesh,
495  const fvMesh& subMesh,
496  const labelList& faceMap,
497  const labelHashSet& addedPatches
498 )
499 {
500  const labelList patchMap(identity(mesh.boundaryMesh().size()));
501 
503  (
504  mesh.objectRegistry::lookupClass<GeoField>()
505  );
507  {
508  const GeoField& fld = *iter();
509 
510  Info<< "Mapping field " << fld.name() << endl;
511 
512  tmp<GeoField> tSubFld
513  (
515  (
516  fld,
517  subMesh,
518  patchMap,
519  faceMap
520  )
521  );
522 
523  // Hack: set value to 0 for introduced patches (since don't
524  // get initialised.
525  forAll(tSubFld().boundaryField(), patchI)
526  {
527  if (addedPatches.found(patchI))
528  {
529  tSubFld().boundaryField()[patchI] ==
531  }
532  }
533 
534  // Store on subMesh
535  GeoField* subFld = tSubFld.ptr();
536  subFld->rename(fld.name());
537  subFld->writeOpt() = IOobject::AUTO_WRITE;
538  subFld->store();
539  }
540 }
541 
542 // Select all cells not in the region
543 labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
544 {
545  DynamicList<label> nonRegionCells(cellRegion.size());
546  forAll(cellRegion, cellI)
547  {
548  if (cellRegion[cellI] != regionI)
549  {
550  nonRegionCells.append(cellI);
551  }
552  }
553  return nonRegionCells.shrink();
554 }
555 
556 
557 // Get per region-region interface the sizes. If sumParallel sums sizes.
558 // Returns interfaces as straight list for looping in identical order.
559 void getInterfaceSizes
560 (
561  const polyMesh& mesh,
562  const labelList& cellRegion,
563  const bool sumParallel,
564 
565  edgeList& interfaces,
566  EdgeMap<label>& interfaceSizes
567 )
568 {
569 
570  // Internal faces
571  // ~~~~~~~~~~~~~~
572 
573  forAll(mesh.faceNeighbour(), faceI)
574  {
575  label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
576  label neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
577 
578  if (ownRegion != neiRegion)
579  {
581  (
582  min(ownRegion, neiRegion),
583  max(ownRegion, neiRegion)
584  );
585 
586  EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
587 
588  if (iter != interfaceSizes.end())
589  {
590  iter()++;
591  }
592  else
593  {
594  interfaceSizes.insert(interface, 1);
595  }
596  }
597  }
598 
599  // Boundary faces
600  // ~~~~~~~~~~~~~~
601 
602  // Neighbour cellRegion.
603  labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
604 
605  forAll(coupledRegion, i)
606  {
607  label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
608  coupledRegion[i] = cellRegion[cellI];
609  }
610  syncTools::swapBoundaryFaceList(mesh, coupledRegion, false);
611 
612  forAll(coupledRegion, i)
613  {
614  label faceI = i+mesh.nInternalFaces();
615  label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
616  label neiRegion = coupledRegion[i];
617 
618  if (ownRegion != neiRegion)
619  {
621  (
622  min(ownRegion, neiRegion),
623  max(ownRegion, neiRegion)
624  );
625 
626  EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
627 
628  if (iter != interfaceSizes.end())
629  {
630  iter()++;
631  }
632  else
633  {
634  interfaceSizes.insert(interface, 1);
635  }
636  }
637  }
638 
639 
640  if (sumParallel && Pstream::parRun())
641  {
642  if (Pstream::master())
643  {
644  // Receive and add to my sizes
645  for
646  (
647  int slave=Pstream::firstSlave();
648  slave<=Pstream::lastSlave();
649  slave++
650  )
651  {
652  IPstream fromSlave(Pstream::blocking, slave);
653 
654  EdgeMap<label> slaveSizes(fromSlave);
655 
656  forAllConstIter(EdgeMap<label>, slaveSizes, slaveIter)
657  {
658  EdgeMap<label>::iterator masterIter =
659  interfaceSizes.find(slaveIter.key());
660 
661  if (masterIter != interfaceSizes.end())
662  {
663  masterIter() += slaveIter();
664  }
665  else
666  {
667  interfaceSizes.insert(slaveIter.key(), slaveIter());
668  }
669  }
670  }
671 
672  // Distribute
673  for
674  (
675  int slave=Pstream::firstSlave();
676  slave<=Pstream::lastSlave();
677  slave++
678  )
679  {
680  // Receive the edges using shared points from the slave.
681  OPstream toSlave(Pstream::blocking, slave);
682  toSlave << interfaceSizes;
683  }
684  }
685  else
686  {
687  // Send to master
688  {
690  toMaster << interfaceSizes;
691  }
692  // Receive from master
693  {
695  fromMaster >> interfaceSizes;
696  }
697  }
698  }
699 
700  // Make sure all processors have interfaces in same order
701  interfaces = interfaceSizes.toc();
702  if (sumParallel)
703  {
704  Pstream::scatter(interfaces);
705  }
706 }
707 
708 
709 // Create mesh for region.
710 autoPtr<mapPolyMesh> createRegionMesh
711 (
712  const labelList& cellRegion,
713  const EdgeMap<label>& interfaceToPatch,
714  const fvMesh& mesh,
715  const label regionI,
716  const word& regionName,
717  autoPtr<fvMesh>& newMesh
718 )
719 {
720  // Create dummy system/fv*
721  {
722  IOobject io
723  (
724  "fvSchemes",
725  mesh.time().system(),
726  regionName,
727  mesh,
730  false
731  );
732 
733  Info<< "Testing:" << io.objectPath() << endl;
734 
735  if (!io.headerOk())
736  // if (!exists(io.objectPath()))
737  {
738  Info<< "Writing dummy " << regionName/io.name() << endl;
739  dictionary dummyDict;
740  dictionary divDict;
741  dummyDict.add("divSchemes", divDict);
742  dictionary gradDict;
743  dummyDict.add("gradSchemes", gradDict);
744  dictionary laplDict;
745  dummyDict.add("laplacianSchemes", laplDict);
746 
747  IOdictionary(io, dummyDict).regIOobject::write();
748  }
749  }
750  {
751  IOobject io
752  (
753  "fvSolution",
754  mesh.time().system(),
755  regionName,
756  mesh,
759  false
760  );
761 
762  if (!io.headerOk())
763  //if (!exists(io.objectPath()))
764  {
765  Info<< "Writing dummy " << regionName/io.name() << endl;
766  dictionary dummyDict;
767  IOdictionary(io, dummyDict).regIOobject::write();
768  }
769  }
770 
771 
772  // Neighbour cellRegion.
773  labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
774 
775  forAll(coupledRegion, i)
776  {
777  label cellI = mesh.faceOwner()[i+mesh.nInternalFaces()];
778  coupledRegion[i] = cellRegion[cellI];
779  }
780  syncTools::swapBoundaryFaceList(mesh, coupledRegion, false);
781 
782 
783  // Topology change container. Start off from existing mesh.
784  polyTopoChange meshMod(mesh);
785 
786  // Cell remover engine
787  removeCells cellRemover(mesh);
788 
789  // Select all but region cells
790  labelList cellsToRemove(getNonRegionCells(cellRegion, regionI));
791 
792  // Find out which faces will get exposed. Note that this
793  // gets faces in mesh face order. So both regions will get same
794  // face in same order (important!)
795  labelList exposedFaces = cellRemover.getExposedFaces(cellsToRemove);
796 
797  labelList exposedPatchIDs(exposedFaces.size());
798  forAll(exposedFaces, i)
799  {
800  label faceI = exposedFaces[i];
801 
802  label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
803  label neiRegion = -1;
804 
805  if (mesh.isInternalFace(faceI))
806  {
807  neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
808  }
809  else
810  {
811  neiRegion = coupledRegion[faceI-mesh.nInternalFaces()];
812  }
813 
814  label otherRegion = -1;
815 
816  if (ownRegion == regionI && neiRegion != regionI)
817  {
818  otherRegion = neiRegion;
819  }
820  else if (ownRegion != regionI && neiRegion == regionI)
821  {
822  otherRegion = ownRegion;
823  }
824  else
825  {
826  FatalErrorIn("createRegionMesh(..)")
827  << "Exposed face:" << faceI
828  << " fc:" << mesh.faceCentres()[faceI]
829  << " has owner region " << ownRegion
830  << " and neighbour region " << neiRegion
831  << " when handling region:" << regionI
832  << exit(FatalError);
833  }
834 
835  if (otherRegion != -1)
836  {
837  edge interface(regionI, otherRegion);
838 
839  // Find the patch.
840  if (regionI < otherRegion)
841  {
842  exposedPatchIDs[i] = interfaceToPatch[interface];
843  }
844  else
845  {
846  exposedPatchIDs[i] = interfaceToPatch[interface]+1;
847  }
848  }
849  }
850 
851  // Remove faces
852  cellRemover.setRefinement
853  (
854  cellsToRemove,
855  exposedFaces,
856  exposedPatchIDs,
857  meshMod
858  );
859 
860  autoPtr<mapPolyMesh> map = meshMod.makeMesh
861  (
862  newMesh,
863  IOobject
864  (
865  regionName,
866  mesh.time().timeName(),
867  mesh.time(),
870  ),
871  mesh
872  );
873 
874  return map;
875 }
876 
877 
878 void createAndWriteRegion
879 (
880  const fvMesh& mesh,
881  const labelList& cellRegion,
882  const wordList& regionNames,
883  const EdgeMap<label>& interfaceToPatch,
884  const label regionI,
885  const word& newMeshInstance
886 )
887 {
888  Info<< "Creating mesh for region " << regionI
889  << ' ' << regionNames[regionI] << endl;
890 
891  autoPtr<fvMesh> newMesh;
892  autoPtr<mapPolyMesh> map = createRegionMesh
893  (
894  cellRegion,
895  interfaceToPatch,
896  mesh,
897  regionI,
898  regionNames[regionI],
899  newMesh
900  );
901 
902 
903  // Make map of all added patches
904  labelHashSet addedPatches(2*interfaceToPatch.size());
905  forAllConstIter(EdgeMap<label>, interfaceToPatch, iter)
906  {
907  addedPatches.insert(iter());
908  addedPatches.insert(iter()+1);
909  }
910 
911  Info<< "Mapping fields" << endl;
912 
913  // Map existing fields
914  newMesh().updateMesh(map());
915 
916  // Add subsetted fields
917  subsetVolFields<volScalarField>
918  (
919  mesh,
920  newMesh(),
921  map().cellMap(),
922  map().faceMap(),
923  addedPatches
924  );
925  subsetVolFields<volVectorField>
926  (
927  mesh,
928  newMesh(),
929  map().cellMap(),
930  map().faceMap(),
931  addedPatches
932  );
933  subsetVolFields<volSphericalTensorField>
934  (
935  mesh,
936  newMesh(),
937  map().cellMap(),
938  map().faceMap(),
939  addedPatches
940  );
941  subsetVolFields<volSymmTensorField>
942  (
943  mesh,
944  newMesh(),
945  map().cellMap(),
946  map().faceMap(),
947  addedPatches
948  );
949  subsetVolFields<volTensorField>
950  (
951  mesh,
952  newMesh(),
953  map().cellMap(),
954  map().faceMap(),
955  addedPatches
956  );
957 
958  subsetSurfaceFields<surfaceScalarField>
959  (
960  mesh,
961  newMesh(),
962  map().faceMap(),
963  addedPatches
964  );
965  subsetSurfaceFields<surfaceVectorField>
966  (
967  mesh,
968  newMesh(),
969  map().faceMap(),
970  addedPatches
971  );
972  subsetSurfaceFields<surfaceSphericalTensorField>
973  (
974  mesh,
975  newMesh(),
976  map().faceMap(),
977  addedPatches
978  );
979  subsetSurfaceFields<surfaceSymmTensorField>
980  (
981  mesh,
982  newMesh(),
983  map().faceMap(),
984  addedPatches
985  );
986  subsetSurfaceFields<surfaceTensorField>
987  (
988  mesh,
989  newMesh(),
990  map().faceMap(),
991  addedPatches
992  );
993 
994 
995  const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
996  newPatches.checkParallelSync(true);
997 
998  // Delete empty patches
999  // ~~~~~~~~~~~~~~~~~~~~
1000 
1001  // Create reordering list to move patches-to-be-deleted to end
1002  labelList oldToNew(newPatches.size(), -1);
1003  label newI = 0;
1004 
1005  Info<< "Deleting empty patches" << endl;
1006 
1007  // Assumes all non-proc boundaries are on all processors!
1008  forAll(newPatches, patchI)
1009  {
1010  const polyPatch& pp = newPatches[patchI];
1011 
1012  if (!isA<processorPolyPatch>(pp))
1013  {
1014  if (returnReduce(pp.size(), sumOp<label>()) > 0)
1015  {
1016  oldToNew[patchI] = newI++;
1017  }
1018  }
1019  }
1020 
1021  // Same for processor patches (but need no reduction)
1022  forAll(newPatches, patchI)
1023  {
1024  const polyPatch& pp = newPatches[patchI];
1025 
1026  if (isA<processorPolyPatch>(pp) && pp.size())
1027  {
1028  oldToNew[patchI] = newI++;
1029  }
1030  }
1031 
1032  const label nNewPatches = newI;
1033 
1034  // Move all deleteable patches to the end
1035  forAll(oldToNew, patchI)
1036  {
1037  if (oldToNew[patchI] == -1)
1038  {
1039  oldToNew[patchI] = newI++;
1040  }
1041  }
1042 
1043  reorderPatches(newMesh(), oldToNew, nNewPatches);
1044 
1045 
1046  Info<< "Writing new mesh" << endl;
1047 
1048  newMesh().setInstance(newMeshInstance);
1049  newMesh().write();
1050 
1051  // Write addressing files like decomposePar
1052  Info<< "Writing addressing to base mesh" << endl;
1053 
1054  labelIOList pointProcAddressing
1055  (
1056  IOobject
1057  (
1058  "pointRegionAddressing",
1059  newMesh().facesInstance(),
1060  newMesh().meshSubDir,
1061  newMesh(),
1062  IOobject::NO_READ,
1063  IOobject::NO_WRITE,
1064  false
1065  ),
1066  map().pointMap()
1067  );
1068  Info<< "Writing map " << pointProcAddressing.name()
1069  << " from region" << regionI
1070  << " points back to base mesh." << endl;
1071  pointProcAddressing.write();
1072 
1073  labelIOList faceProcAddressing
1074  (
1075  IOobject
1076  (
1077  "faceRegionAddressing",
1078  newMesh().facesInstance(),
1079  newMesh().meshSubDir,
1080  newMesh(),
1081  IOobject::NO_READ,
1082  IOobject::NO_WRITE,
1083  false
1084  ),
1085  newMesh().nFaces()
1086  );
1087  forAll(faceProcAddressing, faceI)
1088  {
1089  // face + turning index. (see decomposePar)
1090  // Is the face pointing in the same direction?
1091  label oldFaceI = map().faceMap()[faceI];
1092 
1093  if
1094  (
1095  map().cellMap()[newMesh().faceOwner()[faceI]]
1096  == mesh.faceOwner()[oldFaceI]
1097  )
1098  {
1099  faceProcAddressing[faceI] = oldFaceI+1;
1100  }
1101  else
1102  {
1103  faceProcAddressing[faceI] = -(oldFaceI+1);
1104  }
1105  }
1106  Info<< "Writing map " << faceProcAddressing.name()
1107  << " from region" << regionI
1108  << " faces back to base mesh." << endl;
1109  faceProcAddressing.write();
1110 
1111  labelIOList cellProcAddressing
1112  (
1113  IOobject
1114  (
1115  "cellRegionAddressing",
1116  newMesh().facesInstance(),
1117  newMesh().meshSubDir,
1118  newMesh(),
1119  IOobject::NO_READ,
1120  IOobject::NO_WRITE,
1121  false
1122  ),
1123  map().cellMap()
1124  );
1125  Info<< "Writing map " <<cellProcAddressing.name()
1126  << " from region" << regionI
1127  << " cells back to base mesh." << endl;
1128  cellProcAddressing.write();
1129 }
1130 
1131 
1132 // Create for every region-region interface with non-zero size two patches.
1133 // First one is for minimumregion to maximumregion.
1134 // Note that patches get created in same order on all processors (if parallel)
1135 // since looping over synchronised 'interfaces'.
1136 EdgeMap<label> addRegionPatches
1137 (
1138  fvMesh& mesh,
1139  const labelList& cellRegion,
1140  const label nCellRegions,
1141  const edgeList& interfaces,
1142  const EdgeMap<label>& interfaceSizes,
1143  const wordList& regionNames
1144 )
1145 {
1146  // Check that all patches are present in same order.
1147  mesh.boundaryMesh().checkParallelSync(true);
1148 
1149  Info<< nl << "Adding patches" << nl << endl;
1150 
1151  EdgeMap<label> interfaceToPatch(nCellRegions);
1152 
1153  forAll(interfaces, interI)
1154  {
1155  const edge& e = interfaces[interI];
1156 
1157  if (interfaceSizes[e] > 0)
1158  {
1159  const word inter1 = regionNames[e[0]] + "_to_" + regionNames[e[1]];
1160  const word inter2 = regionNames[e[1]] + "_to_" + regionNames[e[0]];
1161 
1163  (
1164  inter1,
1165  0, // overridden
1166  0, // overridden
1167  0, // overridden
1168  regionNames[e[1]], // sampleRegion
1170  inter2, // samplePatch
1171  point::zero, // offset
1172  mesh.boundaryMesh()
1173  );
1174 
1175  label patchI = addPatch(mesh, patch1);
1176 
1178  (
1179  inter2,
1180  0,
1181  0,
1182  0,
1183  regionNames[e[0]], // sampleRegion
1185  inter1,
1186  point::zero, // offset
1187  mesh.boundaryMesh()
1188  );
1189  addPatch(mesh, patch2);
1190 
1191  Info<< "For interface between region " << e[0]
1192  << " and " << e[1] << " added patch " << patchI
1193  << " " << mesh.boundaryMesh()[patchI].name()
1194  << endl;
1195 
1196  interfaceToPatch.insert(e, patchI);
1197  }
1198  }
1199  return interfaceToPatch;
1200 }
1201 
1202 
1203 // Find region that covers most of cell zone
1204 label findCorrespondingRegion
1205 (
1206  const labelList& existingZoneID, // per cell the (unique) zoneID
1207  const labelList& cellRegion,
1208  const label nCellRegions,
1209  const label zoneI,
1210  const label minOverlapSize
1211 )
1212 {
1213  // Per region the number of cells in zoneI
1214  labelList cellsInZone(nCellRegions, 0);
1215 
1216  forAll(cellRegion, cellI)
1217  {
1218  if (existingZoneID[cellI] == zoneI)
1219  {
1220  cellsInZone[cellRegion[cellI]]++;
1221  }
1222  }
1223 
1225  Pstream::listCombineScatter(cellsInZone);
1226 
1227  // Pick region with largest overlap of zoneI
1228  label regionI = findMax(cellsInZone);
1229 
1230 
1231  if (cellsInZone[regionI] < minOverlapSize)
1232  {
1233  // Region covers too little of zone. Not good enough.
1234  regionI = -1;
1235  }
1236  else
1237  {
1238  // Check that region contains no cells that aren't in cellZone.
1239  forAll(cellRegion, cellI)
1240  {
1241  if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1242  {
1243  // cellI in regionI but not in zoneI
1244  regionI = -1;
1245  break;
1246  }
1247  }
1248  // If one in error, all should be in error. Note that branch gets taken
1249  // on all procs.
1250  reduce(regionI, minOp<label>());
1251  }
1252 
1253  return regionI;
1254 }
1255 
1256 
1258 //label findCorrespondingRegion
1259 //(
1260 // const cellZoneMesh& cellZones,
1261 // const labelList& existingZoneID, // per cell the (unique) zoneID
1262 // const labelList& cellRegion,
1263 // const label nCellRegions,
1264 // const label zoneI
1265 //)
1266 //{
1267 // // Region corresponding to zone. Start off with special value: no
1268 // // corresponding region.
1269 // label regionI = labelMax;
1270 //
1271 // const cellZone& cz = cellZones[zoneI];
1272 //
1273 // if (cz.empty())
1274 // {
1275 // // My local portion is empty. Maps to any empty cellZone. Mark with
1276 // // special value which can get overwritten by other processors.
1277 // regionI = -1;
1278 // }
1279 // else
1280 // {
1281 // regionI = cellRegion[cz[0]];
1282 //
1283 // forAll(cz, i)
1284 // {
1285 // if (cellRegion[cz[i]] != regionI)
1286 // {
1287 // regionI = labelMax;
1288 // break;
1289 // }
1290 // }
1291 // }
1292 //
1293 // // Determine same zone over all processors.
1294 // reduce(regionI, maxOp<label>());
1295 //
1296 //
1297 // // 2. All of region present?
1298 //
1299 // if (regionI == labelMax)
1300 // {
1301 // regionI = -1;
1302 // }
1303 // else if (regionI != -1)
1304 // {
1305 // forAll(cellRegion, cellI)
1306 // {
1307 // if (cellRegion[cellI] == regionI && existingZoneID[cellI] != zoneI)
1308 // {
1309 // // cellI in regionI but not in zoneI
1310 // regionI = -1;
1311 // break;
1312 // }
1313 // }
1314 // // If one in error, all should be in error. Note that branch gets taken
1315 // // on all procs.
1316 // reduce(regionI, minOp<label>());
1317 // }
1318 //
1319 // return regionI;
1320 //}
1321 
1322 
1323 // Get zone per cell
1324 // - non-unique zoning
1325 // - coupled zones
1326 void getZoneID
1327 (
1328  const polyMesh& mesh,
1329  const cellZoneMesh& cellZones,
1330  labelList& zoneID,
1331  labelList& neiZoneID
1332 )
1333 {
1334  // Existing zoneID
1335  zoneID.setSize(mesh.nCells());
1336  zoneID = -1;
1337 
1338  forAll(cellZones, zoneI)
1339  {
1340  const cellZone& cz = cellZones[zoneI];
1341 
1342  forAll(cz, i)
1343  {
1344  label cellI = cz[i];
1345  if (zoneID[cellI] == -1)
1346  {
1347  zoneID[cellI] = zoneI;
1348  }
1349  else
1350  {
1351  FatalErrorIn("getZoneID(..)")
1352  << "Cell " << cellI << " with cell centre "
1353  << mesh.cellCentres()[cellI]
1354  << " is multiple zones. This is not allowed." << endl
1355  << "It is in zone " << cellZones[zoneID[cellI]].name()
1356  << " and in zone " << cellZones[zoneI].name()
1357  << exit(FatalError);
1358  }
1359  }
1360  }
1361 
1362  // Neighbour zoneID.
1363  neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
1364 
1365  forAll(neiZoneID, i)
1366  {
1367  neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
1368  }
1369  syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
1370 }
1371 
1372 
1373 void matchRegions
1374 (
1375  const bool sloppyCellZones,
1376  const polyMesh& mesh,
1377 
1378  const label nCellRegions,
1379  const labelList& cellRegion,
1380 
1381  labelList& regionToZone,
1382  wordList& regionNames,
1383  labelList& zoneToRegion
1384 )
1385 {
1386  const cellZoneMesh& cellZones = mesh.cellZones();
1387 
1388  regionToZone.setSize(nCellRegions, -1);
1389  regionNames.setSize(nCellRegions);
1390  zoneToRegion.setSize(cellZones.size(), -1);
1391 
1392  // Get current per cell zoneID
1393  labelList zoneID(mesh.nCells(), -1);
1394  labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1395  getZoneID(mesh, cellZones, zoneID, neiZoneID);
1396 
1397  // Sizes per cellzone
1398  labelList zoneSizes(cellZones.size(), 0);
1399  {
1400  List<wordList> zoneNames(Pstream::nProcs());
1401  zoneNames[Pstream::myProcNo()] = cellZones.names();
1402  Pstream::gatherList(zoneNames);
1403  Pstream::scatterList(zoneNames);
1404 
1405  forAll(zoneNames, procI)
1406  {
1407  if (zoneNames[procI] != zoneNames[0])
1408  {
1409  FatalErrorIn("matchRegions(..)")
1410  << "cellZones not synchronised across processors." << endl
1411  << "Master has cellZones " << zoneNames[0] << endl
1412  << "Processor " << procI
1413  << " has cellZones " << zoneNames[procI]
1414  << exit(FatalError);
1415  }
1416  }
1417 
1418  forAll(cellZones, zoneI)
1419  {
1420  zoneSizes[zoneI] = returnReduce
1421  (
1422  cellZones[zoneI].size(),
1423  sumOp<label>()
1424  );
1425  }
1426  }
1427 
1428 
1429  if (sloppyCellZones)
1430  {
1431  Info<< "Trying to match regions to existing cell zones;"
1432  << " region can be subset of cell zone." << nl << endl;
1433 
1434  forAll(cellZones, zoneI)
1435  {
1436  label regionI = findCorrespondingRegion
1437  (
1438  zoneID,
1439  cellRegion,
1440  nCellRegions,
1441  zoneI,
1442  label(0.5*zoneSizes[zoneI]) // minimum overlap
1443  );
1444 
1445  if (regionI != -1)
1446  {
1447  Info<< "Sloppily matched region " << regionI
1448  //<< " size " << regionSizes[regionI]
1449  << " to zone " << zoneI << " size " << zoneSizes[zoneI]
1450  << endl;
1451  zoneToRegion[zoneI] = regionI;
1452  regionToZone[regionI] = zoneI;
1453  regionNames[regionI] = cellZones[zoneI].name();
1454  }
1455  }
1456  }
1457  else
1458  {
1459  Info<< "Trying to match regions to existing cell zones." << nl << endl;
1460 
1461  forAll(cellZones, zoneI)
1462  {
1463  label regionI = findCorrespondingRegion
1464  (
1465  zoneID,
1466  cellRegion,
1467  nCellRegions,
1468  zoneI,
1469  1 // minimum overlap
1470  );
1471 
1472  if (regionI != -1)
1473  {
1474  zoneToRegion[zoneI] = regionI;
1475  regionToZone[regionI] = zoneI;
1476  regionNames[regionI] = cellZones[zoneI].name();
1477  }
1478  }
1479  }
1480  // Allocate region names for unmatched regions.
1481  forAll(regionToZone, regionI)
1482  {
1483  if (regionToZone[regionI] == -1)
1484  {
1485  regionNames[regionI] = "domain" + Foam::name(regionI);
1486  }
1487  }
1488 }
1489 
1490 
1491 void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
1492 {
1493  // Write to manual decomposition option
1494  {
1495  labelIOList cellToRegion
1496  (
1497  IOobject
1498  (
1499  "cellToRegion",
1500  mesh.facesInstance(),
1501  mesh,
1504  false
1505  ),
1506  cellRegion
1507  );
1508  cellToRegion.write();
1509 
1510  Info<< "Writing region per cell file (for manual decomposition) to "
1511  << cellToRegion.objectPath() << nl << endl;
1512  }
1513  // Write for postprocessing
1514  {
1515  volScalarField cellToRegion
1516  (
1517  IOobject
1518  (
1519  "cellToRegion",
1520  mesh.time().timeName(),
1521  mesh,
1524  false
1525  ),
1526  mesh,
1527  dimensionedScalar("zero", dimless, 0),
1528  zeroGradientFvPatchScalarField::typeName
1529  );
1530  forAll(cellRegion, cellI)
1531  {
1532  cellToRegion[cellI] = cellRegion[cellI];
1533  }
1534  cellToRegion.write();
1535 
1536  Info<< "Writing region per cell as volScalarField to "
1537  << cellToRegion.objectPath() << nl << endl;
1538  }
1539 }
1540 
1541 
1542 
1543 // Main program:
1544 
1545 int main(int argc, char *argv[])
1546 {
1547  argList::validOptions.insert("cellZones", "");
1548  argList::validOptions.insert("cellZonesOnly", "");
1549  argList::validOptions.insert("cellZonesFileOnly", "cellZonesName");
1550  argList::validOptions.insert("blockedFaces", "faceSet");
1551  argList::validOptions.insert("makeCellZones", "");
1552  argList::validOptions.insert("largestOnly", "");
1553  argList::validOptions.insert("insidePoint", "point");
1554  argList::validOptions.insert("overwrite", "");
1555  argList::validOptions.insert("detectOnly", "");
1556  argList::validOptions.insert("sloppyCellZones", "");
1557 
1558 # include <OpenFOAM/setRootCase.H>
1559 # include <OpenFOAM/createTime.H>
1560  runTime.functionObjects().off();
1561 # include <OpenFOAM/createMesh.H>
1562  const word oldInstance = mesh.pointsInstance();
1563 
1564  word blockedFacesName;
1565  if (args.optionFound("blockedFaces"))
1566  {
1567  blockedFacesName = args.option("blockedFaces");
1568  Info<< "Reading blocked internal faces from faceSet "
1569  << blockedFacesName << nl << endl;
1570  }
1571 
1572  bool makeCellZones = args.optionFound("makeCellZones");
1573  bool largestOnly = args.optionFound("largestOnly");
1574  bool insidePoint = args.optionFound("insidePoint");
1575  bool useCellZones = args.optionFound("cellZones");
1576  bool useCellZonesOnly = args.optionFound("cellZonesOnly");
1577  bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
1578  bool overwrite = args.optionFound("overwrite");
1579  bool detectOnly = args.optionFound("detectOnly");
1580  bool sloppyCellZones = args.optionFound("sloppyCellZones");
1581 
1582  if
1583  (
1584  (useCellZonesOnly || useCellZonesFile)
1585  && (
1586  blockedFacesName != word::null
1587  || useCellZones
1588  )
1589  )
1590  {
1592  << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
1593  << " (which specify complete split)"
1594  << " in combination with -blockedFaces or -cellZones"
1595  << " (which imply a split based on topology)"
1596  << exit(FatalError);
1597  }
1598 
1599 
1600 
1601  if (insidePoint && largestOnly)
1602  {
1604  << "You cannot specify both -largestOnly"
1605  << " (keep region with most cells)"
1606  << " and -insidePoint (keep region containing point)"
1607  << exit(FatalError);
1608  }
1609 
1610 
1611  const cellZoneMesh& cellZones = mesh.cellZones();
1612 
1613  // Existing zoneID
1614  labelList zoneID(mesh.nCells(), -1);
1615  // Neighbour zoneID.
1616  labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1617  getZoneID(mesh, cellZones, zoneID, neiZoneID);
1618 
1619 
1620 
1621  // Determine per cell the region it belongs to
1622  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1623 
1624  // cellRegion is the labelList with the region per cell.
1625  labelList cellRegion;
1626  // Region per zone
1627  labelList regionToZone;
1628  // Name of region
1629  wordList regionNames;
1630  // Zone to region
1631  labelList zoneToRegion;
1632 
1633  label nCellRegions = 0;
1634  if (useCellZonesOnly)
1635  {
1636  Info<< "Using current cellZones to split mesh into regions."
1637  << " This requires all"
1638  << " cells to be in one and only one cellZone." << nl << endl;
1639 
1640  label unzonedCellI = findIndex(zoneID, -1);
1641  if (unzonedCellI != -1)
1642  {
1644  << "For the cellZonesOnly option all cells "
1645  << "have to be in a cellZone." << endl
1646  << "Cell " << unzonedCellI
1647  << " at" << mesh.cellCentres()[unzonedCellI]
1648  << " is not in a cellZone. There might be more unzoned cells."
1649  << exit(FatalError);
1650  }
1651  cellRegion = zoneID;
1652  nCellRegions = gMax(cellRegion)+1;
1653  regionToZone.setSize(nCellRegions);
1654  regionNames.setSize(nCellRegions);
1655  zoneToRegion.setSize(cellZones.size(), -1);
1656  for (label regionI = 0; regionI < nCellRegions; regionI++)
1657  {
1658  regionToZone[regionI] = regionI;
1659  zoneToRegion[regionI] = regionI;
1660  regionNames[regionI] = cellZones[regionI].name();
1661  }
1662  }
1663  else if (useCellZonesFile)
1664  {
1665  const word zoneFile = args.option("cellZonesFileOnly");
1666  Info<< "Reading split from cellZones file " << zoneFile << endl
1667  << "This requires all"
1668  << " cells to be in one and only one cellZone." << nl << endl;
1669 
1670  cellZoneMesh newCellZones
1671  (
1672  IOobject
1673  (
1674  zoneFile,
1675  mesh.facesInstance(),
1677  mesh,
1680  false
1681  ),
1682  mesh
1683  );
1684 
1685  labelList newZoneID(mesh.nCells(), -1);
1686  labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
1687  getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
1688 
1689  label unzonedCellI = findIndex(newZoneID, -1);
1690  if (unzonedCellI != -1)
1691  {
1693  << "For the cellZonesFileOnly option all cells "
1694  << "have to be in a cellZone." << endl
1695  << "Cell " << unzonedCellI
1696  << " at" << mesh.cellCentres()[unzonedCellI]
1697  << " is not in a cellZone. There might be more unzoned cells."
1698  << exit(FatalError);
1699  }
1700  cellRegion = newZoneID;
1701  nCellRegions = gMax(cellRegion)+1;
1702  zoneToRegion.setSize(newCellZones.size(), -1);
1703  regionToZone.setSize(nCellRegions);
1704  regionNames.setSize(nCellRegions);
1705  for (label regionI = 0; regionI < nCellRegions; regionI++)
1706  {
1707  regionToZone[regionI] = regionI;
1708  zoneToRegion[regionI] = regionI;
1709  regionNames[regionI] = newCellZones[regionI].name();
1710  }
1711  }
1712  else
1713  {
1714  // Determine connected regions
1715  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1716 
1717  // Mark additional faces that are blocked
1718  boolList blockedFace;
1719 
1720  // Read from faceSet
1721  if (blockedFacesName.size())
1722  {
1723  faceSet blockedFaceSet(mesh, blockedFacesName);
1724  Info<< "Read "
1725  << returnReduce(blockedFaceSet.size(), sumOp<label>())
1726  << " blocked faces from set " << blockedFacesName << nl << endl;
1727 
1728  blockedFace.setSize(mesh.nFaces(), false);
1729 
1730  forAllConstIter(faceSet, blockedFaceSet, iter)
1731  {
1732  blockedFace[iter.key()] = true;
1733  }
1734  }
1735 
1736  // Imply from differing cellZones
1737  if (useCellZones)
1738  {
1739  blockedFace.setSize(mesh.nFaces(), false);
1740 
1741  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
1742  {
1743  label own = mesh.faceOwner()[faceI];
1744  label nei = mesh.faceNeighbour()[faceI];
1745 
1746  if (zoneID[own] != zoneID[nei])
1747  {
1748  blockedFace[faceI] = true;
1749  }
1750  }
1751 
1752  // Different cellZones on either side of processor patch.
1753  forAll(neiZoneID, i)
1754  {
1755  label faceI = i+mesh.nInternalFaces();
1756 
1757  if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
1758  {
1759  blockedFace[faceI] = true;
1760  }
1761  }
1762  }
1763 
1764  // Do a topological walk to determine regions
1765  regionSplit regions(mesh, blockedFace);
1766  nCellRegions = regions.nRegions();
1767  cellRegion.transfer(regions);
1768 
1769  // Make up region names. If possible match them to existing zones.
1770  matchRegions
1771  (
1772  sloppyCellZones,
1773  mesh,
1774  nCellRegions,
1775  cellRegion,
1776 
1777  regionToZone,
1778  regionNames,
1779  zoneToRegion
1780  );
1781  }
1782 
1783  Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
1784 
1785 
1786  // Write decomposition to file
1787  writeCellToRegion(mesh, cellRegion);
1788 
1789 
1790 
1791  // Sizes per region
1792  // ~~~~~~~~~~~~~~~~
1793 
1794  labelList regionSizes(nCellRegions, 0);
1795 
1796  forAll(cellRegion, cellI)
1797  {
1798  regionSizes[cellRegion[cellI]]++;
1799  }
1800  forAll(regionSizes, regionI)
1801  {
1802  reduce(regionSizes[regionI], sumOp<label>());
1803  }
1804 
1805  Info<< "Region\tCells" << nl
1806  << "------\t-----" << endl;
1807 
1808  forAll(regionSizes, regionI)
1809  {
1810  Info<< regionI << '\t' << regionSizes[regionI] << nl;
1811  }
1812  Info<< endl;
1813 
1814 
1815 
1816  // Print region to zone
1817  Info<< "Region\tZone\tName" << nl
1818  << "------\t----\t----" << endl;
1819  forAll(regionToZone, regionI)
1820  {
1821  Info<< regionI << '\t' << regionToZone[regionI] << '\t'
1822  << regionNames[regionI] << nl;
1823  }
1824  Info<< endl;
1825 
1826 
1827 
1828  // Since we're going to mess with patches make sure all non-processor ones
1829  // are on all processors.
1830  mesh.boundaryMesh().checkParallelSync(true);
1831 
1832 
1833  // Sizes of interface between regions. From pair of regions to number of
1834  // faces.
1835  edgeList interfaces;
1836  EdgeMap<label> interfaceSizes;
1837  getInterfaceSizes
1838  (
1839  mesh,
1840  cellRegion,
1841  true, // sum in parallel?
1842 
1843  interfaces,
1844  interfaceSizes
1845  );
1846 
1847  Info<< "Sizes inbetween regions:" << nl << nl
1848  << "Region\tRegion\tFaces" << nl
1849  << "------\t------\t-----" << endl;
1850 
1851  forAll(interfaces, interI)
1852  {
1853  const edge& e = interfaces[interI];
1854 
1855  Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
1856  }
1857  Info<< endl;
1858 
1859 
1860  if (detectOnly)
1861  {
1862  return 0;
1863  }
1864 
1865 
1866  // Read objects in time directory
1867  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1868 
1869  IOobjectList objects(mesh, runTime.timeName());
1870 
1871  // Read vol fields.
1872 
1873  PtrList<volScalarField> vsFlds;
1874  ReadFields(mesh, objects, vsFlds);
1875 
1876  PtrList<volVectorField> vvFlds;
1877  ReadFields(mesh, objects, vvFlds);
1878 
1880  ReadFields(mesh, objects, vstFlds);
1881 
1882  PtrList<volSymmTensorField> vsymtFlds;
1883  ReadFields(mesh, objects, vsymtFlds);
1884 
1885  PtrList<volTensorField> vtFlds;
1886  ReadFields(mesh, objects, vtFlds);
1887 
1888  // Read surface fields.
1889 
1891  ReadFields(mesh, objects, ssFlds);
1892 
1894  ReadFields(mesh, objects, svFlds);
1895 
1897  ReadFields(mesh, objects, sstFlds);
1898 
1900  ReadFields(mesh, objects, ssymtFlds);
1901 
1903  ReadFields(mesh, objects, stFlds);
1904 
1905  Info<< endl;
1906 
1907 
1908  // Remove any demand-driven fields ('S', 'V' etc)
1909  mesh.clearOut();
1910 
1911 
1912  if (nCellRegions == 1)
1913  {
1914  Info<< "Only one region. Doing nothing." << endl;
1915  }
1916  else if (makeCellZones)
1917  {
1918  Info<< "Putting cells into cellZones instead of splitting mesh."
1919  << endl;
1920 
1921  // Check if region overlaps with existing zone. If so keep.
1922 
1923  for (label regionI = 0; regionI < nCellRegions; regionI++)
1924  {
1925  label zoneI = regionToZone[regionI];
1926 
1927  if (zoneI != -1)
1928  {
1929  Info<< " Region " << regionI << " : corresponds to existing"
1930  << " cellZone "
1931  << zoneI << ' ' << cellZones[zoneI].name() << endl;
1932  }
1933  else
1934  {
1935  // Create new cellZone.
1936  labelList regionCells = findIndices(cellRegion, regionI);
1937 
1938  word zoneName = "region" + Foam::name(regionI);
1939 
1940  zoneI = cellZones.findZoneID(zoneName);
1941 
1942  if (zoneI == -1)
1943  {
1944  zoneI = cellZones.size();
1945  mesh.cellZones().setSize(zoneI+1);
1946  mesh.cellZones().set
1947  (
1948  zoneI,
1949  new cellZone
1950  (
1951  zoneName, //name
1952  regionCells, //addressing
1953  zoneI, //index
1954  cellZones //cellZoneMesh
1955  )
1956  );
1957  }
1958  else
1959  {
1960  mesh.cellZones()[zoneI].clearAddressing();
1961  mesh.cellZones()[zoneI] = regionCells;
1962  }
1963  Info<< " Region " << regionI << " : created new cellZone "
1964  << zoneI << ' ' << cellZones[zoneI].name() << endl;
1965  }
1966  }
1968 
1969  if (!overwrite)
1970  {
1971  runTime++;
1972  mesh.setInstance(runTime.timeName());
1973  }
1974  else
1975  {
1976  mesh.setInstance(oldInstance);
1977  }
1978 
1979  Info<< "Writing cellZones as new mesh to time " << runTime.timeName()
1980  << nl << endl;
1981 
1982  mesh.write();
1983 
1984 
1985  // Write cellSets for convenience
1986  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1987 
1988  Info<< "Writing cellSets corresponding to cellZones." << nl << endl;
1989 
1990  forAll(cellZones, zoneI)
1991  {
1992  const cellZone& cz = cellZones[zoneI];
1993 
1994  cellSet(mesh, cz.name(), cz).write();
1995  }
1996  }
1997  else
1998  {
1999  // Add patches for interfaces
2000  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
2001 
2002  // Add all possible patches. Empty ones get filtered later on.
2003  Info<< nl << "Adding patches" << nl << endl;
2004 
2005  EdgeMap<label> interfaceToPatch
2006  (
2007  addRegionPatches
2008  (
2009  mesh,
2010  cellRegion,
2011  nCellRegions,
2012  interfaces,
2013  interfaceSizes,
2014  regionNames
2015  )
2016  );
2017 
2018 
2019  if (!overwrite)
2020  {
2021  runTime++;
2022  }
2023 
2024 
2025  // Create regions
2026  // ~~~~~~~~~~~~~~
2027 
2028  if (insidePoint)
2029  {
2030  point insidePoint(args.optionLookup("insidePoint")());
2031 
2032  label regionI = -1;
2033 
2034  label cellI = mesh.findCell(insidePoint);
2035 
2036  Info<< nl << "Found point " << insidePoint << " in cell " << cellI
2037  << endl;
2038 
2039  if (cellI != -1)
2040  {
2041  regionI = cellRegion[cellI];
2042  }
2043 
2044  reduce(regionI, maxOp<label>());
2045 
2046  Info<< nl
2047  << "Subsetting region " << regionI
2048  << " containing point " << insidePoint << endl;
2049 
2050  if (regionI == -1)
2051  {
2053  << "Point " << insidePoint
2054  << " is not inside the mesh." << nl
2055  << "Bounding box of the mesh:" << mesh.bounds()
2056  << exit(FatalError);
2057  }
2058 
2059  createAndWriteRegion
2060  (
2061  mesh,
2062  cellRegion,
2063  regionNames,
2064  interfaceToPatch,
2065  regionI,
2066  (overwrite ? oldInstance : runTime.timeName())
2067  );
2068  }
2069  else if (largestOnly)
2070  {
2071  label regionI = findMax(regionSizes);
2072 
2073  Info<< nl
2074  << "Subsetting region " << regionI
2075  << " of size " << regionSizes[regionI] << endl;
2076 
2077  createAndWriteRegion
2078  (
2079  mesh,
2080  cellRegion,
2081  regionNames,
2082  interfaceToPatch,
2083  regionI,
2084  (overwrite ? oldInstance : runTime.timeName())
2085  );
2086  }
2087  else
2088  {
2089  // Split all
2090  for (label regionI = 0; regionI < nCellRegions; regionI++)
2091  {
2092  Info<< nl
2093  << "Region " << regionI << nl
2094  << "-------- " << endl;
2095 
2096  createAndWriteRegion
2097  (
2098  mesh,
2099  cellRegion,
2100  regionNames,
2101  interfaceToPatch,
2102  regionI,
2103  (overwrite ? oldInstance : runTime.timeName())
2104  );
2105  }
2106  }
2107  }
2108 
2109  Info<< "End\n" << endl;
2110 
2111  return 0;
2112 }
2113 
2114 
2115 // ************************ vim: set sw=4 sts=4 et: ************************ //