meshRefinementBaffles.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-2014 OpenFOAM Foundation
9  Copyright (C) 2015-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
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 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "meshRefinement.H"
30 #include "fvMesh.H"
31 #include "syncTools.H"
32 #include "Time.H"
33 #include "refinementSurfaces.H"
34 #include "faceSet.H"
35 #include "indirectPrimitivePatch.H"
36 #include "polyTopoChange.H"
37 #include "meshTools.H"
38 #include "polyModifyFace.H"
39 #include "polyModifyCell.H"
40 #include "polyAddFace.H"
41 #include "polyRemoveFace.H"
42 #include "localPointRegion.H"
43 #include "duplicatePoints.H"
44 #include "regionSplit.H"
45 #include "removeCells.H"
46 #include "unitConversion.H"
47 #include "OBJstream.H"
48 #include "patchFaceOrientation.H"
49 #include "PatchEdgeFaceWave.H"
50 #include "edgeTopoDistanceData.H"
51 #include "polyMeshAdder.H"
52 #include "IOmanip.H"
53 #include "refinementParameters.H"
54 #include "shellSurfaces.H"
55 #include "zeroGradientFvPatchFields.H"
56 #include "volFields.H"
57 
58 #include "FaceCellWave.H"
59 #include "wallPoints.H"
60 
61 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
62 
63 Foam::label Foam::meshRefinement::createBaffle
64 (
65  const label faceI,
66  const label ownPatch,
67  const label neiPatch,
68  polyTopoChange& meshMod
69 ) const
70 {
71  const face& f = mesh_.faces()[faceI];
72  label zoneID = mesh_.faceZones().whichZone(faceI);
73  bool zoneFlip = false;
74 
75  if (zoneID >= 0)
76  {
77  const faceZone& fZone = mesh_.faceZones()[zoneID];
78  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
79  }
80 
81  meshMod.setAction
82  (
83  polyModifyFace
84  (
85  f, // modified face
86  faceI, // label of face
87  mesh_.faceOwner()[faceI], // owner
88  -1, // neighbour
89  false, // face flip
90  ownPatch, // patch for face
91  false, // remove from zone
92  zoneID, // zone for face
93  zoneFlip // face flip in zone
94  )
95  );
96 
97 
98  label dupFaceI = -1;
99 
100  if (mesh_.isInternalFace(faceI))
101  {
102  if (neiPatch == -1)
103  {
104  FatalErrorInFunction
105  << "No neighbour patch for internal face " << faceI
106  << " fc:" << mesh_.faceCentres()[faceI]
107  << " ownPatch:" << ownPatch << abort(FatalError);
108  }
109 
110  bool reverseFlip = false;
111  if (zoneID >= 0)
112  {
113  reverseFlip = !zoneFlip;
114  }
115 
116  dupFaceI = meshMod.setAction
117  (
118  polyAddFace
119  (
120  f.reverseFace(), // modified face
121  mesh_.faceNeighbour()[faceI],// owner
122  -1, // neighbour
123  -1, // masterPointID
124  -1, // masterEdgeID
125  faceI, // masterFaceID,
126  true, // face flip
127  neiPatch, // patch for face
128  zoneID, // zone for face
129  reverseFlip // face flip in zone
130  )
131  );
132  }
133  return dupFaceI;
134 }
135 
136 
143 //bool Foam::meshRefinement::validBaffleTopology
144 //(
145 // const label faceI,
146 // const vector& n1,
147 // const vector& n2,
148 // const vector& testDir
149 //) const
150 //{
151 //
152 // label patchI = mesh_.boundaryMesh().whichPatch(faceI);
153 // if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
154 // {
155 // return true;
156 // }
157 // else if (mag(n1&n2) > cos(degToRad(30.0)))
158 // {
159 // // Both normals aligned. Check that test vector perpendicularish to
160 // // surface normal
161 // scalar magTestDir = mag(testDir);
162 // if (magTestDir > VSMALL)
163 // {
164 // if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45.0)))
165 // {
166 // //Pout<< "** disabling baffling face "
167 // // << mesh_.faceCentres()[faceI] << endl;
168 // return false;
169 // }
170 // }
171 // }
172 // return true;
173 //}
174 
175 
176 void Foam::meshRefinement::getIntersections
177 (
178  const labelList& surfacesToTest,
179  const pointField& neiCc,
180  const labelList& testFaces,
181 
182  labelList& globalRegion1,
183  labelList& globalRegion2
184 ) const
185 {
186  autoPtr<OBJstream> str;
187  if (debug&OBJINTERSECTIONS)
188  {
189  mkDir(mesh_.time().path()/timeName());
190  str.reset
191  (
192  new OBJstream
193  (
194  mesh_.time().path()/timeName()/"intersections.obj"
195  )
196  );
197 
198  Pout<< "getIntersections : Writing surface intersections to file "
199  << str().name() << nl << endl;
200  }
201 
202 
203  globalRegion1.setSize(mesh_.nFaces());
204  globalRegion1 = -1;
205  globalRegion2.setSize(mesh_.nFaces());
206  globalRegion2 = -1;
207 
208 
209  // Collect segments
210  // ~~~~~~~~~~~~~~~~
211 
212  pointField start(testFaces.size());
213  pointField end(testFaces.size());
214 
215  {
216  labelList minLevel;
217  calcCellCellRays
218  (
219  neiCc,
220  labelList(neiCc.size(), -1),
221  testFaces,
222  start,
223  end,
224  minLevel
225  );
226  }
227 
228 
229  // Do test for intersections
230  // ~~~~~~~~~~~~~~~~~~~~~~~~~
231 
232  labelList surface1;
233  List<pointIndexHit> hit1;
234  labelList region1;
235  labelList surface2;
236  List<pointIndexHit> hit2;
237  labelList region2;
238  surfaces_.findNearestIntersection
239  (
240  surfacesToTest,
241  start,
242  end,
243 
244  surface1,
245  hit1,
246  region1,
247 
248  surface2,
249  hit2,
250  region2
251  );
252 
253 
254  forAll(testFaces, i)
255  {
256  label faceI = testFaces[i];
257 
258  if (hit1[i].hit() && hit2[i].hit())
259  {
260  if (str)
261  {
262  str().write(linePointRef(start[i], hit1[i].rawPoint()));
263  str().write
264  (
265  linePointRef(hit1[i].rawPoint(), hit2[i].rawPoint())
266  );
267  str().write(linePointRef(hit2[i].rawPoint(), end[i]));
268  }
269 
270  // Pick up the patches
271  globalRegion1[faceI] =
272  surfaces_.globalRegion(surface1[i], region1[i]);
273  globalRegion2[faceI] =
274  surfaces_.globalRegion(surface2[i], region2[i]);
275 
276  if (globalRegion1[faceI] == -1 || globalRegion2[faceI] == -1)
277  {
278  FatalErrorInFunction
279  << "problem." << abort(FatalError);
280  }
281  }
282  }
283 }
284 
285 
286 void Foam::meshRefinement::getBafflePatches
287 (
288  const label nErodeCellZones,
289  const labelList& globalToMasterPatch,
290  const pointField& locationsInMesh,
291  const wordList& zonesInMesh,
292 
293  const labelList& neiLevel,
294  const pointField& neiCc,
295 
296  labelList& ownPatch,
297  labelList& neiPatch
298 ) const
299 {
300  // This determines the patches for the intersected faces to
301  // - remove the outside of the mesh
302  // - introduce baffles for (non-faceZone) intersections
303  // Any baffles for faceZones (faceType 'baffle'/'boundary') get introduced
304  // later
305 
306 
307  // 1. Determine intersections with unnamed surfaces and cell zones
308  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
309 
310  labelList cellToZone;
311  labelList unnamedRegion1;
312  labelList unnamedRegion2;
313  labelList namedSurfaceRegion;
314  {
315  bitSet posOrientation;
316  zonify
317  (
318  true, // allowFreeStandingZoneFaces
319  nErodeCellZones,
320  -2, // zone to put unreached cells into
321  locationsInMesh,
322  zonesInMesh,
323 
324  cellToZone,
325  unnamedRegion1,
326  unnamedRegion2,
327  namedSurfaceRegion,
328  posOrientation
329  );
330  }
331 
332 
333  // 2. Baffle all boundary faces
334  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
335  // The logic is that all intersections with unnamed surfaces become
336  // baffles except where we're inbetween a cellZone and background
337  // or inbetween two different cellZones.
338 
339  labelList neiCellZone;
340  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
341 
342  ownPatch.setSize(mesh_.nFaces());
343  ownPatch = -1;
344  neiPatch.setSize(mesh_.nFaces());
345  neiPatch = -1;
346 
347  forAll(ownPatch, faceI)
348  {
349  if (unnamedRegion1[faceI] != -1 || unnamedRegion2[faceI] != -1)
350  {
351  label ownMasterPatch = -1;
352  if (unnamedRegion1[faceI] != -1)
353  {
354  ownMasterPatch = globalToMasterPatch[unnamedRegion1[faceI]];
355  }
356  label neiMasterPatch = -1;
357  if (unnamedRegion2[faceI] != -1)
358  {
359  neiMasterPatch = globalToMasterPatch[unnamedRegion2[faceI]];
360  }
361 
362 
363  // Now we always want to produce a baffle except if
364  // one side is a valid cellZone
365 
366  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
367  label neiZone = -1;
368 
369  if (mesh_.isInternalFace(faceI))
370  {
371  neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
372  }
373  else
374  {
375  neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
376  }
377 
378 
379  if
380  (
381  (ownZone != neiZone)
382  && (
383  (ownZone >= 0 && neiZone != -2)
384  || (neiZone >= 0 && ownZone != -2)
385  )
386  && (
387  namedSurfaceRegion.size() == 0
388  || namedSurfaceRegion[faceI] == -1
389  )
390  )
391  {
392  // One side is a valid cellZone and the neighbour is a different
393  // one (or -1 but not -2). Do not baffle unless the user has
394  // put both an unnamed and a named surface there. In that
395  // case assume that the user wants both: baffle and faceZone.
396  }
397  else
398  {
399  ownPatch[faceI] = ownMasterPatch;
400  neiPatch[faceI] = neiMasterPatch;
401  }
402  }
403  }
404 
405  // No need to parallel sync since intersection data (surfaceIndex_ etc.)
406  // already guaranteed to be synced...
407  // However:
408  // - owncc and neicc are reversed on different procs so might pick
409  // up different regions reversed? No problem. Neighbour on one processor
410  // might not be owner on the other processor but the neighbour is
411  // not used when creating baffles from proc faces.
412  // - tolerances issues occasionally crop up.
413  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
414  syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>());
415 }
416 
417 
418 Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches
419 (
420  const bool allowBoundary,
421  const labelList& globalToMasterPatch,
422  const labelList& globalToSlavePatch
423 ) const
424 {
425  Map<labelPair> bafflePatch(mesh_.nFaces()/1000);
426 
427  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
428  const faceZoneMesh& fZones = mesh_.faceZones();
429 
430  forAll(surfZones, surfI)
431  {
432  const wordList& faceZoneNames = surfZones[surfI].faceZoneNames();
433 
434  forAll(faceZoneNames, fzi)
435  {
436  // Get zone
437  const word& faceZoneName = faceZoneNames[fzi];
438  label zoneI = fZones.findZoneID(faceZoneName);
439  const faceZone& fZone = fZones[zoneI];
440 
441  // Get patch allocated for zone
442  label globalRegionI = surfaces_.globalRegion(surfI, fzi);
443  labelPair zPatches
444  (
445  globalToMasterPatch[globalRegionI],
446  globalToSlavePatch[globalRegionI]
447  );
448 
449  Info<< "For zone " << fZone.name() << " found patches "
450  << mesh_.boundaryMesh()[zPatches[0]].name() << " and "
451  << mesh_.boundaryMesh()[zPatches[1]].name()
452  << endl;
453 
454  forAll(fZone, i)
455  {
456  label faceI = fZone[i];
457 
458  if (allowBoundary || mesh_.isInternalFace(faceI))
459  {
460  labelPair patches = zPatches;
461  if (fZone.flipMap()[i])
462  {
463  patches = reverse(patches);
464  }
465 
466  if (!bafflePatch.insert(faceI, patches))
467  {
468  FatalErrorInFunction
469  << "Face " << faceI
470  << " fc:" << mesh_.faceCentres()[faceI]
471  << " in zone " << fZone.name()
472  << " is in multiple zones!"
473  << abort(FatalError);
474  }
475  }
476  }
477  }
478  }
479  return bafflePatch;
480 }
481 
482 
483 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
484 (
485  const labelList& ownPatch,
486  const labelList& neiPatch
487 )
488 {
489  if
490  (
491  ownPatch.size() != mesh_.nFaces()
492  || neiPatch.size() != mesh_.nFaces()
493  )
494  {
495  FatalErrorInFunction
496  << "Illegal size :"
497  << " ownPatch:" << ownPatch.size()
498  << " neiPatch:" << neiPatch.size()
499  << ". Should be number of faces:" << mesh_.nFaces()
500  << abort(FatalError);
501  }
502 
503  if (debug)
504  {
505  labelList syncedOwnPatch(ownPatch);
506  syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>());
507  labelList syncedNeiPatch(neiPatch);
508  syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>());
509 
510  forAll(syncedOwnPatch, faceI)
511  {
512  if
513  (
514  (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
515  || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
516  )
517  {
518  FatalErrorInFunction
519  << "Non synchronised at face:" << faceI
520  << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
521  << " fc:" << mesh_.faceCentres()[faceI] << endl
522  << "ownPatch:" << ownPatch[faceI]
523  << " syncedOwnPatch:" << syncedOwnPatch[faceI]
524  << " neiPatch:" << neiPatch[faceI]
525  << " syncedNeiPatch:" << syncedNeiPatch[faceI]
526  << abort(FatalError);
527  }
528  }
529  }
530 
531  // Topochange container
532  polyTopoChange meshMod(mesh_);
533 
534  label nBaffles = 0;
535 
536  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
537  {
538  if (ownPatch[faceI] != -1)
539  {
540  // Create baffle or repatch face. Return label of inserted baffle
541  // face.
542  createBaffle
543  (
544  faceI,
545  ownPatch[faceI], // owner side patch
546  neiPatch[faceI], // neighbour side patch
547  meshMod
548  );
549  nBaffles++;
550  }
551  }
552  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
553 
554  forAll(pbm, patchI)
555  {
556  const polyPatch& pp = pbm[patchI];
557 
558  label coupledPatchI = -1;
559  if
560  (
561  pp.coupled()
562  && !refCast<const coupledPolyPatch>(pp).separated()
563  )
564  {
565  coupledPatchI = patchI;
566  }
567 
568  forAll(pp, i)
569  {
570  label faceI = pp.start()+i;
571 
572  if (ownPatch[faceI] != -1)
573  {
574  createBaffle
575  (
576  faceI,
577  ownPatch[faceI], // owner side patch
578  neiPatch[faceI], // neighbour side patch
579  meshMod
580  );
581 
582  if (coupledPatchI != -1)
583  {
584  faceToCoupledPatch_.insert(faceI, coupledPatchI);
585  }
586 
587  nBaffles++;
588  }
589  }
590  }
591 
592 
593  autoPtr<mapPolyMesh> mapPtr;
594  if (returnReduce(nBaffles, sumOp<label>()))
595  {
596  // Change the mesh (no inflation, parallel sync)
597  mapPtr = meshMod.changeMesh(mesh_, false, true);
598  mapPolyMesh& map = *mapPtr;
599 
600  // Update fields
601  mesh_.updateMesh(map);
602 
603  // Move mesh if in inflation mode
604  if (map.hasMotionPoints())
605  {
606  mesh_.movePoints(map.preMotionPoints());
607  }
608  else
609  {
610  // Delete mesh volumes.
611  mesh_.clearOut();
612  }
613 
614 
615  // Reset the instance for if in overwrite mode
616  mesh_.setInstance(timeName());
617 
618  //- Redo the intersections on the newly create baffle faces. Note that
619  // this changes also the cell centre positions.
620  faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
621 
622  const labelList& reverseFaceMap = map.reverseFaceMap();
623  const labelList& faceMap = map.faceMap();
624 
625  // Pick up owner side of baffle
626  forAll(ownPatch, oldFaceI)
627  {
628  label faceI = reverseFaceMap[oldFaceI];
629 
630  if (ownPatch[oldFaceI] != -1 && faceI >= 0)
631  {
632  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
633 
634  baffledFacesSet.insert(ownFaces);
635  }
636  }
637  // Pick up neighbour side of baffle (added faces)
638  forAll(faceMap, faceI)
639  {
640  label oldFaceI = faceMap[faceI];
641 
642  if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
643  {
644  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
645 
646  baffledFacesSet.insert(ownFaces);
647  }
648  }
649  baffledFacesSet.sync(mesh_);
650 
651  updateMesh(map, baffledFacesSet.toc());
652  }
653 
654  return mapPtr;
655 }
656 
657 
658 Foam::labelList Foam::meshRefinement::getZones
659 (
660  const List<surfaceZonesInfo::faceZoneType>& fzTypes
661 ) const
662 {
663  const faceZoneMesh& faceZones = mesh_.faceZones();
664 
665  DynamicList<label> zoneIDs(faceZones.size());
666 
667  forAll(faceZones, zoneI)
668  {
669  const faceZone& fZone = faceZones[zoneI];
670 
671  label mpI, spI;
672  surfaceZonesInfo::faceZoneType fzType;
673  bool hasInfo = getFaceZoneInfo(fZone.name(), mpI, spI, fzType);
674 
675  if (hasInfo && fzTypes.found(fzType))
676  {
677  zoneIDs.append(zoneI);
678  }
679  }
680  return zoneIDs;
681 }
682 
683 
684 // Subset those baffles where both faces are on the same zone
685 Foam::List<Foam::labelPair> Foam::meshRefinement::subsetBaffles
686 (
687  const polyMesh& mesh,
688  const labelList& zoneIDs,
689  const List<labelPair>& baffles
690 )
691 {
692  const faceZoneMesh& faceZones = mesh.faceZones();
693 
694  // Mark zone per face
695  labelList faceToZone(mesh.nFaces(), -1);
696 
697  for (const label zoneID : zoneIDs)
698  {
699  labelUIndList(faceToZone, faceZones[zoneID]) = zoneID;
700  }
701 
702 
703  // Subset baffles
704  DynamicList<labelPair> newBaffles(baffles.size());
705  forAll(baffles, i)
706  {
707  const labelPair& p = baffles[i];
708  if (faceToZone[p[0]] != -1 && (faceToZone[p[0]] == faceToZone[p[1]]))
709  {
710  newBaffles.append(p);
711  }
712  }
713 
714  return newBaffles;
715 }
716 
717 
718 void Foam::meshRefinement::getZoneFaces
719 (
720  const labelList& zoneIDs,
721  labelList& faceZoneID,
722  labelList& ownPatch,
723  labelList& neiPatch,
724  labelList& nBaffles
725 ) const
726 {
727  const faceZoneMesh& faceZones = mesh_.faceZones();
728 
729  // Per (internal) face the patch related to the faceZone
730  ownPatch.setSize(mesh_.nFaces());
731  ownPatch= -1;
732  neiPatch.setSize(mesh_.nFaces());
733  neiPatch = -1;
734  faceZoneID.setSize(mesh_.nFaces());
735  faceZoneID = -1;
736  nBaffles.setSize(zoneIDs.size());
737  nBaffles = Zero;
738 
739  forAll(zoneIDs, j)
740  {
741  label zoneI = zoneIDs[j];
742  const faceZone& fz = faceZones[zoneI];
743  const word& masterName = faceZoneToMasterPatch_[fz.name()];
744  label masterPatchI = mesh_.boundaryMesh().findPatchID(masterName);
745  const word& slaveName = faceZoneToSlavePatch_[fz.name()];
746  label slavePatchI = mesh_.boundaryMesh().findPatchID(slaveName);
747 
748  if (masterPatchI == -1 || slavePatchI == -1)
749  {
750  FatalErrorInFunction
751  << "Problem: masterPatchI:" << masterPatchI
752  << " slavePatchI:" << slavePatchI << exit(FatalError);
753  }
754 
755  forAll(fz, i)
756  {
757  label faceI = fz[i];
758  if (mesh_.isInternalFace(faceI))
759  {
760  if (fz.flipMap()[i])
761  {
762  ownPatch[faceI] = slavePatchI;
763  neiPatch[faceI] = masterPatchI;
764  }
765  else
766  {
767  ownPatch[faceI] = masterPatchI;
768  neiPatch[faceI] = slavePatchI;
769  }
770  faceZoneID[faceI] = zoneI;
771 
772  nBaffles[j]++;
773  }
774  }
775  }
776 }
777 
778 
779 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
780 (
781  const labelList& zoneIDs,
782  List<labelPair>& baffles,
783  labelList& originatingFaceZone
784 )
785 {
786  autoPtr<mapPolyMesh> map;
787 
788  if (zoneIDs.size() > 0)
789  {
790  const faceZoneMesh& faceZones = mesh_.faceZones();
791 
792  // Split internal faces on interface surfaces
793  Info<< "Converting zoned faces into baffles ..." << endl;
794 
795  // Get faceZone and patch(es) per face (or -1 if face not on faceZone)
796  labelList faceZoneID;
797  labelList ownPatch;
798  labelList neiPatch;
799  labelList nBaffles;
800  getZoneFaces(zoneIDs, faceZoneID, ownPatch, neiPatch, nBaffles);
801 
802  label nLocalBaffles = sum(nBaffles);
803 
804 
805  label nTotalBaffles = returnReduce(nLocalBaffles, sumOp<label>());
806 
807  if (nTotalBaffles > 0)
808  {
809  Pstream::listCombineGather(nBaffles, plusEqOp<label>());
810  Pstream::listCombineScatter(nBaffles);
811 
812  Info<< nl
813  << setf(ios_base::left)
814  << setw(30) << "FaceZone"
815  << setw(10) << "FaceType"
816  << setw(10) << "nBaffles"
817  << nl
818  << setw(30) << "--------"
819  << setw(10) << "--------"
820  << setw(10) << "--------"
821  << endl;
822 
823  forAll(zoneIDs, j)
824  {
825  label zoneI = zoneIDs[j];
826  const faceZone& fz = faceZones[zoneI];
827 
828  label mpI, spI;
829  surfaceZonesInfo::faceZoneType fzType;
830  bool hasInfo = getFaceZoneInfo(fz.name(), mpI, spI, fzType);
831 
832  if (hasInfo)
833  {
834  Info<< setf(ios_base::left)
835  << setw(30) << fz.name()
836  << setw(10)
837  << surfaceZonesInfo::faceZoneTypeNames[fzType]
838  << setw(10) << nBaffles[j]
839  << nl;
840  }
841  }
842  Info<< endl;
843 
844  // Create baffles.
845  map = createBaffles(ownPatch, neiPatch);
846 
847  // Get pairs of faces created.
848  // Just loop over faceMap and store baffle if we encounter a slave
849  // face.
850 
851  baffles.setSize(nLocalBaffles);
852  originatingFaceZone.setSize(nLocalBaffles);
853  label baffleI = 0;
854 
855  const labelList& faceMap = map().faceMap();
856  const labelList& reverseFaceMap = map().reverseFaceMap();
857 
858  for
859  (
860  label faceI = mesh_.nInternalFaces();
861  faceI < mesh_.nFaces();
862  faceI++
863  )
864  {
865  label oldFaceI = faceMap[faceI];
866  label masterFaceI = reverseFaceMap[oldFaceI];
867  if (masterFaceI != faceI && ownPatch[oldFaceI] != -1)
868  {
869  baffles[baffleI] = labelPair(masterFaceI, faceI);
870  originatingFaceZone[baffleI] = faceZoneID[oldFaceI];
871  baffleI++;
872  }
873  }
874 
875  if (baffleI != baffles.size())
876  {
877  FatalErrorInFunction
878  << "Had " << baffles.size() << " baffles to create "
879  << " but encountered " << baffleI
880  << " slave faces originating from patcheable faces."
881  << abort(FatalError);
882  }
883 
884  if (debug&MESH)
885  {
886  const_cast<Time&>(mesh_.time())++;
887  Pout<< "Writing zone-baffled mesh to time " << timeName()
888  << endl;
889  write
890  (
891  debugType(debug),
892  writeType(writeLevel() | WRITEMESH),
893  mesh_.time().path()/"baffles"
894  );
895  }
896  }
897  Info<< "Created " << nTotalBaffles << " baffles in = "
898  << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
899  }
900  else
901  {
902  baffles.clear();
903  originatingFaceZone.clear();
904  }
905 
906  return map;
907 }
908 
909 
910 Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
911 (
912  const List<labelPair>& couples,
913  const scalar planarAngle
914 ) const
915 {
916  // Done by counting the number of baffles faces per mesh edge. If edge
917  // has 2 boundary faces and both are baffle faces it is the edge of a baffle
918  // region.
919 
920  // All duplicate faces on edge of the patch are to be merged.
921  // So we count for all edges of duplicate faces how many duplicate
922  // faces use them.
923  labelList nBafflesPerEdge(mesh_.nEdges(), Zero);
924 
925 
926  // This algorithm is quite tricky. We don't want to use edgeFaces and
927  // also want it to run in parallel so it is now an algorithm over
928  // all (boundary) faces instead.
929  // We want to pick up any edges that are only used by the baffle
930  // or internal faces but not by any other boundary faces. So
931  // - increment count on an edge by 1 if it is used by any (uncoupled)
932  // boundary face.
933  // - increment count on an edge by 1000000 if it is used by a baffle face
934  // - sum in parallel
935  //
936  // So now any edge that is used by baffle faces only will have the
937  // value 2*1000000+2*1.
938 
939 
940  const label baffleValue = 1000000;
941 
942 
943  // Count number of boundary faces per edge
944  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
945 
946  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
947 
948  forAll(patches, patchI)
949  {
950  const polyPatch& pp = patches[patchI];
951 
952  // Count number of boundary faces. Discard coupled boundary faces.
953  if (!pp.coupled())
954  {
955  label faceI = pp.start();
956 
957  forAll(pp, i)
958  {
959  const labelList& fEdges = mesh_.faceEdges(faceI);
960 
961  forAll(fEdges, fEdgeI)
962  {
963  nBafflesPerEdge[fEdges[fEdgeI]]++;
964  }
965  faceI++;
966  }
967  }
968  }
969 
970 
971  DynamicList<label> fe0;
972  DynamicList<label> fe1;
973 
974 
975  // Count number of duplicate boundary faces per edge
976  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
977 
978  forAll(couples, i)
979  {
980  {
981  label f0 = couples[i].first();
982  const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
983  forAll(fEdges0, fEdgeI)
984  {
985  nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
986  }
987  }
988 
989  {
990  label f1 = couples[i].second();
991  const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
992  forAll(fEdges1, fEdgeI)
993  {
994  nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
995  }
996  }
997  }
998 
999  // Add nBaffles on shared edges
1000  syncTools::syncEdgeList
1001  (
1002  mesh_,
1003  nBafflesPerEdge,
1004  plusEqOp<label>(), // in-place add
1005  label(0) // initial value
1006  );
1007 
1008 
1009  // Baffles which are not next to other boundaries and baffles will have
1010  // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which
1011  // are both baffle faces)
1012 
1013  List<labelPair> filteredCouples(couples.size());
1014  label filterI = 0;
1015 
1016  forAll(couples, i)
1017  {
1018  const labelPair& couple = couples[i];
1019 
1020  if
1021  (
1022  patches.whichPatch(couple.first())
1023  == patches.whichPatch(couple.second())
1024  )
1025  {
1026  const labelList& fEdges = mesh_.faceEdges(couple.first());
1027 
1028  forAll(fEdges, fEdgeI)
1029  {
1030  label edgeI = fEdges[fEdgeI];
1031 
1032  if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
1033  {
1034  filteredCouples[filterI++] = couple;
1035  break;
1036  }
1037  }
1038  }
1039  }
1040  filteredCouples.setSize(filterI);
1041 
1042 
1043  label nFiltered = returnReduce(filteredCouples.size(), sumOp<label>());
1044 
1045  Info<< "freeStandingBaffles : detected "
1046  << nFiltered
1047  << " free-standing baffles out of "
1048  << returnReduce(couples.size(), sumOp<label>())
1049  << nl << endl;
1050 
1051 
1052  if (nFiltered > 0)
1053  {
1054  // Collect segments
1055  // ~~~~~~~~~~~~~~~~
1056 
1057  pointField start(filteredCouples.size());
1058  pointField end(filteredCouples.size());
1059 
1060  const pointField& cellCentres = mesh_.cellCentres();
1061 
1062  forAll(filteredCouples, i)
1063  {
1064  const labelPair& couple = filteredCouples[i];
1065  start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
1066  end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
1067  }
1068 
1069  // Extend segments a bit
1070  {
1071  const vectorField smallVec(ROOTSMALL*(end-start));
1072  start -= smallVec;
1073  end += smallVec;
1074  }
1075 
1076 
1077  // Do test for intersections
1078  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1079  labelList surface1;
1080  List<pointIndexHit> hit1;
1081  labelList region1;
1082  vectorField normal1;
1083 
1084  labelList surface2;
1085  List<pointIndexHit> hit2;
1086  labelList region2;
1087  vectorField normal2;
1088 
1089  surfaces_.findNearestIntersection
1090  (
1091  identity(surfaces_.surfaces().size()),
1092  start,
1093  end,
1094 
1095  surface1,
1096  hit1,
1097  region1,
1098  normal1,
1099 
1100  surface2,
1101  hit2,
1102  region2,
1103  normal2
1104  );
1105 
1106  //mkDir(mesh_.time().path()/timeName());
1107  //OBJstream str
1108  //(
1109  // mesh_.time().path()/timeName()/"flatBaffles.obj"
1110  //);
1111 
1112  const scalar planarAngleCos = Foam::cos(degToRad(planarAngle));
1113 
1114  label filterI = 0;
1115  forAll(filteredCouples, i)
1116  {
1117  const labelPair& couple = filteredCouples[i];
1118 
1119  // Note: for a baffle-surface we do not want to merge the baffle.
1120  // We could either check for hitting the same triangle (but you
1121  // might hit same point on neighbouring triangles due to tolerance
1122  // issues) or better just to compare the hit point.
1123  // This might still go wrong for a ray in the plane of the triangle
1124  // which might hit two different points on the same triangle due
1125  // to tolerances...
1126 
1127  if
1128  (
1129  hit1[i].hit()
1130  && hit2[i].hit()
1131  && mag(hit1[i].hitPoint()-hit2[i].hitPoint()) > mergeDistance_
1132  )
1133  {
1134  // Two different hits. Check angle.
1135  //str.write
1136  //(
1137  // linePointRef(hit1[i].hitPoint(), hit2[i].hitPoint()),
1138  // normal1[i],
1139  // normal2[i]
1140  //);
1141 
1142  if ((normal1[i]&normal2[i]) > planarAngleCos)
1143  {
1144  // Both normals aligned
1145  vector n = end[i]-start[i];
1146  scalar magN = mag(n);
1147  if (magN > VSMALL)
1148  {
1149  filteredCouples[filterI++] = couple;
1150  }
1151  }
1152  }
1153  else if (hit1[i].hit() || hit2[i].hit())
1154  {
1155  // Single hit. Do not include in freestanding baffles.
1156  }
1157  }
1158 
1159  filteredCouples.setSize(filterI);
1160 
1161  Info<< "freeStandingBaffles : detected "
1162  << returnReduce(filterI, sumOp<label>())
1163  << " planar (within " << planarAngle
1164  << " degrees) free-standing baffles out of "
1165  << nFiltered
1166  << nl << endl;
1167  }
1168 
1169  return filteredCouples;
1170 }
1171 
1172 
1173 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
1174 (
1175  const List<labelPair>& couples,
1176  const Map<label>& faceToPatch
1177 )
1178 {
1179  autoPtr<mapPolyMesh> mapPtr;
1180 
1181  if (returnReduce(couples.size()+faceToPatch.size(), sumOp<label>()))
1182  {
1183  // Mesh change engine
1184  polyTopoChange meshMod(mesh_);
1185 
1186  const faceList& faces = mesh_.faces();
1187  const labelList& faceOwner = mesh_.faceOwner();
1188  const faceZoneMesh& faceZones = mesh_.faceZones();
1189 
1190  forAll(couples, i)
1191  {
1192  label face0 = couples[i].first();
1193  label face1 = couples[i].second();
1194 
1195  // face1 < 0 signals a coupled face that has been converted to
1196  // baffle
1197 
1198  label own0 = faceOwner[face0];
1199  label own1 = faceOwner[face1];
1200 
1201  if (face1 < 0 || own0 < own1)
1202  {
1203  // Use face0 as the new internal face.
1204  label zoneID = faceZones.whichZone(face0);
1205  bool zoneFlip = false;
1206 
1207  if (zoneID >= 0)
1208  {
1209  const faceZone& fZone = faceZones[zoneID];
1210  zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
1211  }
1212 
1213  label nei = (face1 < 0 ? -1 : own1);
1214 
1215  meshMod.setAction(polyRemoveFace(face1));
1216  meshMod.setAction
1217  (
1218  polyModifyFace
1219  (
1220  faces[face0], // modified face
1221  face0, // label of face being modified
1222  own0, // owner
1223  nei, // neighbour
1224  false, // face flip
1225  -1, // patch for face
1226  false, // remove from zone
1227  zoneID, // zone for face
1228  zoneFlip // face flip in zone
1229  )
1230  );
1231  }
1232  else
1233  {
1234  // Use face1 as the new internal face.
1235  label zoneID = faceZones.whichZone(face1);
1236  bool zoneFlip = false;
1237 
1238  if (zoneID >= 0)
1239  {
1240  const faceZone& fZone = faceZones[zoneID];
1241  zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
1242  }
1243 
1244  meshMod.setAction(polyRemoveFace(face0));
1245  meshMod.setAction
1246  (
1247  polyModifyFace
1248  (
1249  faces[face1], // modified face
1250  face1, // label of face being modified
1251  own1, // owner
1252  own0, // neighbour
1253  false, // face flip
1254  -1, // patch for face
1255  false, // remove from zone
1256  zoneID, // zone for face
1257  zoneFlip // face flip in zone
1258  )
1259  );
1260  }
1261  }
1262 
1263  forAllConstIters(faceToPatch, iter)
1264  {
1265  const label faceI = iter.key();
1266  const label patchI = iter.val();
1267 
1268  if (!mesh_.isInternalFace(faceI))
1269  {
1270  FatalErrorInFunction
1271  << "problem: face:" << faceI
1272  << " at:" << mesh_.faceCentres()[faceI]
1273  << "(wanted patch:" << patchI
1274  << ") is an internal face" << exit(FatalError);
1275  }
1276 
1277  label zoneID = faceZones.whichZone(faceI);
1278  bool zoneFlip = false;
1279 
1280  if (zoneID >= 0)
1281  {
1282  const faceZone& fZone = faceZones[zoneID];
1283  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
1284  }
1285 
1286  meshMod.setAction
1287  (
1288  polyModifyFace
1289  (
1290  faces[faceI], // modified face
1291  faceI, // label of face being modified
1292  faceOwner[faceI], // owner
1293  -1, // neighbour
1294  false, // face flip
1295  patchI, // patch for face
1296  false, // remove from zone
1297  zoneID, // zone for face
1298  zoneFlip // face flip in zone
1299  )
1300  );
1301  }
1302 
1303 
1304  // Change the mesh (no inflation)
1305  mapPtr = meshMod.changeMesh(mesh_, false, true);
1306  mapPolyMesh& map = *mapPtr;
1307 
1308  // Update fields
1309  mesh_.updateMesh(map);
1310 
1311  // Move mesh (since morphing does not do this)
1312  if (map.hasMotionPoints())
1313  {
1314  mesh_.movePoints(map.preMotionPoints());
1315  }
1316  else
1317  {
1318  // Delete mesh volumes.
1319  mesh_.clearOut();
1320  }
1321 
1322  // Reset the instance for if in overwrite mode
1323  mesh_.setInstance(timeName());
1324 
1325  // Update intersections. Recalculate intersections on merged faces since
1326  // this seems to give problems? Note: should not be necessary since
1327  // baffles preserve intersections from when they were created.
1328  labelList newExposedFaces(2*couples.size());
1329  label newI = 0;
1330 
1331  forAll(couples, i)
1332  {
1333  const label newFace0 = map.reverseFaceMap()[couples[i].first()];
1334  if (newFace0 != -1)
1335  {
1336  newExposedFaces[newI++] = newFace0;
1337  }
1338 
1339  const label newFace1 = map.reverseFaceMap()[couples[i].second()];
1340  if (newFace1 != -1)
1341  {
1342  newExposedFaces[newI++] = newFace1;
1343  }
1344  }
1345  newExposedFaces.setSize(newI);
1346  updateMesh(map, newExposedFaces);
1347  }
1348 
1349  return mapPtr;
1350 }
1351 
1352 
1353 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeZoneBaffles
1354 (
1355  const bool doInternalZones,
1356  const bool doBaffleZones
1357 )
1358 {
1359  labelList zoneIDs;
1360  {
1361  DynamicList<surfaceZonesInfo::faceZoneType> fzTypes;
1362  if (doInternalZones)
1363  {
1364  fzTypes.append(surfaceZonesInfo::INTERNAL);
1365  }
1366  if (doBaffleZones)
1367  {
1368  fzTypes.append(surfaceZonesInfo::BAFFLE);
1369  }
1370  zoneIDs = getZones(fzTypes);
1371  }
1372 
1373  List<labelPair> zoneBaffles
1374  (
1375  subsetBaffles
1376  (
1377  mesh_,
1378  zoneIDs,
1379  localPointRegion::findDuplicateFacePairs(mesh_)
1380  )
1381  );
1382 
1383  autoPtr<mapPolyMesh> mapPtr;
1384  if (returnReduce(zoneBaffles.size(), sumOp<label>()))
1385  {
1386  mapPtr = mergeBaffles(zoneBaffles, Map<label>(0));
1387  }
1388  return mapPtr;
1389 }
1390 
1391 
1392 // Finds region per cell for cells inside closed named surfaces
1393 void Foam::meshRefinement::findCellZoneGeometric
1394 (
1395  const pointField& neiCc,
1396  const labelList& closedNamedSurfaces, // indices of closed surfaces
1397  labelList& namedSurfaceRegion, // per face: named surface region
1398  const labelList& surfaceToCellZone, // cell zone index per surface
1399 
1400  labelList& cellToZone
1401 ) const
1402 {
1403  const pointField& cellCentres = mesh_.cellCentres();
1404  const labelList& faceOwner = mesh_.faceOwner();
1405  const labelList& faceNeighbour = mesh_.faceNeighbour();
1406 
1407  // Check if cell centre is inside
1408  labelList insideSurfaces;
1409  surfaces_.findInside
1410  (
1411  closedNamedSurfaces,
1412  cellCentres,
1413  insideSurfaces
1414  );
1415 
1416  forAll(insideSurfaces, cellI)
1417  {
1418  label surfI = insideSurfaces[cellI];
1419 
1420  if (surfI != -1)
1421  {
1422  if (cellToZone[cellI] == -2)
1423  {
1424  cellToZone[cellI] = surfaceToCellZone[surfI];
1425  }
1426  else if (cellToZone[cellI] == -1)
1427  {
1428  // ? Allow named surface to override background zone (-1)
1429  // This is used in the multiRegionHeater tutorial where the
1430  // locationInMesh is inside a named surface.
1431  cellToZone[cellI] = surfaceToCellZone[surfI];
1432  }
1433  }
1434  }
1435 
1436 
1437  // Some cells with cell centres close to surface might have
1438  // had been put into wrong surface. Recheck with perturbed cell centre.
1439 
1440 
1441  // 1. Collect points
1442 
1443  // Count points to test.
1444  label nCandidates = 0;
1445  forAll(namedSurfaceRegion, faceI)
1446  {
1447  if (namedSurfaceRegion[faceI] != -1)
1448  {
1449  if (mesh_.isInternalFace(faceI))
1450  {
1451  nCandidates += 2;
1452  }
1453  else
1454  {
1455  nCandidates += 1;
1456  }
1457  }
1458  }
1459 
1460  // Collect points.
1461  pointField candidatePoints(nCandidates);
1462  nCandidates = 0;
1463  forAll(namedSurfaceRegion, faceI)
1464  {
1465  if (namedSurfaceRegion[faceI] != -1)
1466  {
1467  label own = faceOwner[faceI];
1468  const point& ownCc = cellCentres[own];
1469 
1470  if (mesh_.isInternalFace(faceI))
1471  {
1472  label nei = faceNeighbour[faceI];
1473  const point& neiCc = cellCentres[nei];
1474  // Perturbed cc
1475  const vector d = 1e-4*(neiCc - ownCc);
1476  candidatePoints[nCandidates++] = ownCc-d;
1477  candidatePoints[nCandidates++] = neiCc+d;
1478  }
1479  else
1480  {
1481  //const point& neiFc = mesh_.faceCentres()[faceI];
1482  const point& neiFc = neiCc[faceI-mesh_.nInternalFaces()];
1483 
1484  // Perturbed cc
1485  const vector d = 1e-4*(neiFc - ownCc);
1486  candidatePoints[nCandidates++] = ownCc-d;
1487  }
1488  }
1489  }
1490 
1491 
1492  // 2. Test points for inside
1493 
1494  surfaces_.findInside
1495  (
1496  closedNamedSurfaces,
1497  candidatePoints,
1498  insideSurfaces
1499  );
1500 
1501 
1502  // 3. Update zone information
1503 
1504  nCandidates = 0;
1505  forAll(namedSurfaceRegion, faceI)
1506  {
1507  if (namedSurfaceRegion[faceI] != -1)
1508  {
1509  label own = faceOwner[faceI];
1510 
1511  if (mesh_.isInternalFace(faceI))
1512  {
1513  label ownSurfI = insideSurfaces[nCandidates++];
1514  if (ownSurfI != -1)
1515  {
1516  cellToZone[own] = surfaceToCellZone[ownSurfI];
1517  }
1518 
1519  label neiSurfI = insideSurfaces[nCandidates++];
1520  if (neiSurfI != -1)
1521  {
1522  label nei = faceNeighbour[faceI];
1523 
1524  cellToZone[nei] = surfaceToCellZone[neiSurfI];
1525  }
1526  }
1527  else
1528  {
1529  label ownSurfI = insideSurfaces[nCandidates++];
1530  if (ownSurfI != -1)
1531  {
1532  cellToZone[own] = surfaceToCellZone[ownSurfI];
1533  }
1534  }
1535  }
1536  }
1537 
1538 
1539  // Adapt the namedSurfaceRegion
1540  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1541  // for if any cells were not completely covered.
1542 
1543  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1544  {
1545  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1546  label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
1547 
1548  if (namedSurfaceRegion[faceI] == -1 && (ownZone != neiZone))
1549  {
1550  // Give face the zone of min cell zone (but only if the
1551  // cellZone originated from a closed, named surface)
1552 
1553  label minZone;
1554  if (ownZone == -1)
1555  {
1556  minZone = neiZone;
1557  }
1558  else if (neiZone == -1)
1559  {
1560  minZone = ownZone;
1561  }
1562  else
1563  {
1564  minZone = min(ownZone, neiZone);
1565  }
1566 
1567  // Make sure the cellZone originated from a closed surface. Use
1568  // hardcoded region 0 inside named surface.
1569  label geomSurfI = surfaceToCellZone.find(minZone);
1570 
1571  if (geomSurfI != -1)
1572  {
1573  namedSurfaceRegion[faceI] =
1574  surfaces_.globalRegion(geomSurfI, 0);
1575  }
1576  }
1577  }
1578 
1579  labelList neiCellZone;
1580  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
1581 
1582  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1583 
1584  forAll(patches, patchI)
1585  {
1586  const polyPatch& pp = patches[patchI];
1587 
1588  if (pp.coupled())
1589  {
1590  forAll(pp, i)
1591  {
1592  label faceI = pp.start()+i;
1593  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1594  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1595 
1596  if (namedSurfaceRegion[faceI] == -1 && (ownZone != neiZone))
1597  {
1598  // Give face the min cell zone
1599  label minZone;
1600  if (ownZone == -1)
1601  {
1602  minZone = neiZone;
1603  }
1604  else if (neiZone == -1)
1605  {
1606  minZone = ownZone;
1607  }
1608  else
1609  {
1610  minZone = min(ownZone, neiZone);
1611  }
1612 
1613  // Make sure the cellZone originated from a closed surface.
1614  // Use hardcoded region 0 inside named surface.
1615  label geomSurfI = surfaceToCellZone.find(minZone);
1616 
1617  if (geomSurfI != -1)
1618  {
1619  namedSurfaceRegion[faceI] =
1620  surfaces_.globalRegion(geomSurfI, 0);
1621  }
1622  }
1623  }
1624  }
1625  }
1626 
1627  // Sync
1628  syncTools::syncFaceList(mesh_, namedSurfaceRegion, maxEqOp<label>());
1629 }
1630 
1631 
1632 void Foam::meshRefinement::findCellZoneInsideWalk
1633 (
1634  const pointField& locationsInMesh,
1635  const labelList& zonesInMesh,
1636  const labelList& faceToZone, // per face -1 or some index >= 0
1637 
1638  labelList& cellToZone
1639 ) const
1640 {
1641  // Analyse regions. Reuse regionsplit
1642  boolList blockedFace(mesh_.nFaces());
1643  //selectSeparatedCoupledFaces(blockedFace);
1644 
1645  forAll(blockedFace, faceI)
1646  {
1647  if (faceToZone[faceI] == -1)
1648  {
1649  blockedFace[faceI] = false;
1650  }
1651  else
1652  {
1653  blockedFace[faceI] = true;
1654  }
1655  }
1656  // No need to sync since faceToZone already is synced
1657 
1658  // Set region per cell based on walking
1659  regionSplit cellRegion(mesh_, blockedFace);
1660  blockedFace.clear();
1661 
1662  // Force calculation of face decomposition (used in findCell)
1663  (void)mesh_.tetBasePtIs();
1664 
1665  // For all locationsInMesh find the cell
1666  forAll(locationsInMesh, i)
1667  {
1668  // Get location and index of zone ("none" for cellZone -1)
1669  const point& insidePoint = locationsInMesh[i];
1670  label zoneID = zonesInMesh[i];
1671 
1672  // Find the region containing the insidePoint
1673  label keepRegionI = findRegion
1674  (
1675  mesh_,
1676  cellRegion,
1677  mergeDistance_ * vector::one,
1678  insidePoint
1679  );
1680 
1681  Info<< "For cellZone "
1682  << (
1683  zoneID == -1
1684  ? "none"
1685  : mesh_.cellZones()[zoneID].name()
1686  )
1687  << " found point " << insidePoint
1688  << " in global region " << keepRegionI
1689  << " out of " << cellRegion.nRegions() << " regions." << endl;
1690 
1691  if (keepRegionI == -1)
1692  {
1693  FatalErrorInFunction
1694  << "Point " << insidePoint
1695  << " is not inside the mesh." << nl
1696  << "Bounding box of the mesh:" << mesh_.bounds()
1697  << exit(FatalError);
1698  }
1699 
1700 
1701 
1702  // Set all cells with this region to the zoneID
1703  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1704 
1705  label nWarnings = 0;
1706 
1707  forAll(cellRegion, cellI)
1708  {
1709  if (cellRegion[cellI] == keepRegionI)
1710  {
1711  if (cellToZone[cellI] == -2)
1712  {
1713  // First visit of cell
1714  cellToZone[cellI] = zoneID;
1715  }
1716  else if (cellToZone[cellI] != zoneID)
1717  {
1718  if (cellToZone[cellI] != -1 && nWarnings < 10)
1719  {
1720  WarningInFunction
1721  << "Cell " << cellI
1722  << " at " << mesh_.cellCentres()[cellI]
1723  << " is inside cellZone "
1724  << (
1725  zoneID == -1
1726  ? "none"
1727  : mesh_.cellZones()[zoneID].name()
1728  )
1729  << " from locationInMesh " << insidePoint
1730  << " but already marked as being in zone "
1731  << mesh_.cellZones()[cellToZone[cellI]].name()
1732  << endl
1733  << "This can happen if your surfaces are not"
1734  << " (sufficiently) closed."
1735  << endl;
1736  nWarnings++;
1737  }
1738 
1739  // Override the zone
1740  cellToZone[cellI] = zoneID;
1741  }
1742  }
1743  }
1744  }
1745 }
1746 
1747 
1748 void Foam::meshRefinement::findCellZoneInsideWalk
1749 (
1750  const pointField& locationsInMesh,
1751  const wordList& zoneNamesInMesh,
1752  const labelList& faceToZone, // per face -1 or some index >= 0
1753 
1754  labelList& cellToZone
1755 ) const
1756 {
1757  const cellZoneMesh& czs = mesh_.cellZones();
1758 
1759  labelList zoneIDs(zoneNamesInMesh.size());
1760  forAll(zoneNamesInMesh, i)
1761  {
1762  zoneIDs[i] = czs.findZoneID(zoneNamesInMesh[i]);
1763  }
1764  findCellZoneInsideWalk
1765  (
1766  locationsInMesh,
1767  zoneIDs,
1768  faceToZone,
1769  cellToZone
1770  );
1771 }
1772 
1773 
1774 bool Foam::meshRefinement::calcRegionToZone
1775 (
1776  const label backgroundZoneID,
1777  const label surfZoneI,
1778  const label ownRegion,
1779  const label neiRegion,
1780 
1781  labelList& regionToCellZone
1782 ) const
1783 {
1784  bool changed = false;
1785 
1786  // Check whether inbetween different regions
1787  if (ownRegion != neiRegion)
1788  {
1789  // Jump. Change one of the sides to my type.
1790 
1791  // 1. Interface between my type and unset region.
1792  // Set region to keepRegion
1793 
1794  if (regionToCellZone[ownRegion] == -2)
1795  {
1796  if (surfZoneI == -1)
1797  {
1798  // Special: face is -on faceZone -not real boundary
1799  // -not on cellZone
1800  // so make regions same on either side
1801  if (regionToCellZone[neiRegion] != -2)
1802  {
1803  regionToCellZone[ownRegion] = regionToCellZone[neiRegion];
1804  changed = true;
1805  }
1806  }
1807  else if (regionToCellZone[neiRegion] == surfZoneI)
1808  {
1809  // Face between unset and my region. Put unset
1810  // region into background region
1811  //MEJ: see comment in findCellZoneTopo
1812  if (backgroundZoneID != -2)
1813  {
1814  regionToCellZone[ownRegion] = backgroundZoneID;
1815  changed = true;
1816  }
1817  }
1818  else if (regionToCellZone[neiRegion] != -2)
1819  {
1820  // Face between unset and other region.
1821  // Put unset region into my region
1822  regionToCellZone[ownRegion] = surfZoneI;
1823  changed = true;
1824  }
1825  }
1826  else if (regionToCellZone[neiRegion] == -2)
1827  {
1828  if (surfZoneI == -1)
1829  {
1830  // Special: face is -on faceZone -not real boundary
1831  // -not on cellZone
1832  // so make regions same on either side
1833  regionToCellZone[neiRegion] = regionToCellZone[ownRegion];
1834  changed = true;
1835  }
1836  else if (regionToCellZone[ownRegion] == surfZoneI)
1837  {
1838  // Face between unset and my region. Put unset
1839  // region into background region
1840  if (backgroundZoneID != -2)
1841  {
1842  regionToCellZone[neiRegion] = backgroundZoneID;
1843  changed = true;
1844  }
1845  }
1846  else if (regionToCellZone[ownRegion] != -2)
1847  {
1848  // Face between unset and other region.
1849  // Put unset region into my region
1850  regionToCellZone[neiRegion] = surfZoneI;
1851  changed = true;
1852  }
1853  }
1854  }
1855  return changed;
1856 }
1857 
1858 
1859 void Foam::meshRefinement::findCellZoneTopo
1860 (
1861  const label backgroundZoneID,
1862  const pointField& locationsInMesh,
1863  const labelList& unnamedSurfaceRegion,
1864  const labelList& namedSurfaceRegion,
1865  const labelList& surfaceToCellZone,
1866  labelList& cellToZone
1867 ) const
1868 {
1869  // This routine fixes small problems with left over unassigned regions
1870  // (after all off the unreachable bits of the mesh have been removed).
1871  // This routine splits the mesh into regions, based on the intersection
1872  // with a surface. The problem is that we know the surface which
1873  // (intersected) face belongs to (in namedSurfaceRegion) but we don't
1874  // know which side of the face it relates to. So all we are doing here
1875  // is get the correspondence between surface/cellZone and regionSplit
1876  // region. See the logic in calcRegionToZone.
1877  // Basically it looks at the neighbours of a face on a named surface.
1878  // If one side is in the cellZone corresponding to the surface and
1879  // the other side is unassigned (-2) it sets this to the background zone.
1880  // So the zone to set these unmarked cells to is provided as argument:
1881  // - backgroundZoneID = -2 : do not change so remove cells
1882  // - backgroundZoneID = -1 : put into background
1883 
1884  // Assumes:
1885  // - region containing keepPoint does not go into a cellZone
1886  // - all other regions can be found by crossing faces marked in
1887  // namedSurfaceRegion.
1888 
1889  // Analyse regions. Reuse regionsplit
1890  boolList blockedFace(mesh_.nFaces());
1891 
1892  forAll(unnamedSurfaceRegion, faceI)
1893  {
1894  if
1895  (
1896  unnamedSurfaceRegion[faceI] == -1
1897  && namedSurfaceRegion[faceI] == -1
1898  )
1899  {
1900  blockedFace[faceI] = false;
1901  }
1902  else
1903  {
1904  blockedFace[faceI] = true;
1905  }
1906  }
1907  // No need to sync since namedSurfaceRegion already is synced
1908 
1909  // Set region per cell based on walking
1910  regionSplit cellRegion(mesh_, blockedFace);
1911  blockedFace.clear();
1912 
1913  // Per mesh region the zone the cell should be put in.
1914  // -2 : not analysed yet
1915  // -1 : keepPoint region. Not put into any cellzone.
1916  // >= 0 : index of cellZone
1917  labelList regionToCellZone(cellRegion.nRegions(), -2);
1918 
1919  // See which cells already are set in the cellToZone (from geometric
1920  // searching) and use these to take over their zones.
1921  // Note: could be improved to count number of cells per region.
1922  forAll(cellToZone, cellI)
1923  {
1924  label regionI = cellRegion[cellI];
1925  if (cellToZone[cellI] != -2 && regionToCellZone[regionI] == -2)
1926  {
1927  regionToCellZone[regionI] = cellToZone[cellI];
1928  }
1929  }
1930 
1931  // Synchronise regionToCellZone.
1932  // Note:
1933  // - region numbers are identical on all processors
1934  // - keepRegion is identical ,,
1935  // - cellZones are identical ,,
1936  Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
1937  Pstream::listCombineScatter(regionToCellZone);
1938 
1939 
1940  // Find the region containing the keepPoint
1941  forAll(locationsInMesh, i)
1942  {
1943  const point& keepPoint = locationsInMesh[i];
1944  label keepRegionI = findRegion
1945  (
1946  mesh_,
1947  cellRegion,
1948  mergeDistance_ * vector::one,
1949  keepPoint
1950  );
1951 
1952  Info<< "Found point " << keepPoint
1953  << " in global region " << keepRegionI
1954  << " out of " << cellRegion.nRegions() << " regions." << endl;
1955 
1956  if (keepRegionI == -1)
1957  {
1958  FatalErrorInFunction
1959  << "Point " << keepPoint
1960  << " is not inside the mesh." << nl
1961  << "Bounding box of the mesh:" << mesh_.bounds()
1962  << exit(FatalError);
1963  }
1964 
1965  // Mark default region with zone -1. Note that all regions should
1966  // already be matched to a cellZone through the loop over cellToZone.
1967  // This is just to mop up any remaining regions. Not sure whether
1968  // this is needed?
1969  if (regionToCellZone[keepRegionI] == -2)
1970  {
1971  regionToCellZone[keepRegionI] = -1;
1972  }
1973  }
1974 
1975 
1976  // Find correspondence between cell zone and surface
1977  // by changing cell zone every time we cross a surface.
1978  while (true)
1979  {
1980  // Synchronise regionToCellZone.
1981  // Note:
1982  // - region numbers are identical on all processors
1983  // - keepRegion is identical ,,
1984  // - cellZones are identical ,,
1985  // This done at top of loop to account for geometric matching
1986  // not being synchronised.
1987  Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
1988  Pstream::listCombineScatter(regionToCellZone);
1989 
1990 
1991  bool changed = false;
1992 
1993  // Internal faces
1994 
1995  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1996  {
1997  label regionI = namedSurfaceRegion[faceI];
1998 
1999  // Connected even if no cellZone defined for surface
2000  if (unnamedSurfaceRegion[faceI] == -1 && regionI != -1)
2001  {
2002  // Calculate region to zone from cellRegions on either side
2003  // of internal face.
2004 
2005  label surfI = surfaces_.whichSurface(regionI).first();
2006 
2007  bool changedCell = calcRegionToZone
2008  (
2009  backgroundZoneID,
2010  surfaceToCellZone[surfI],
2011  cellRegion[mesh_.faceOwner()[faceI]],
2012  cellRegion[mesh_.faceNeighbour()[faceI]],
2013  regionToCellZone
2014  );
2015 
2016  changed = changed | changedCell;
2017  }
2018  }
2019 
2020  // Boundary faces
2021 
2022  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2023 
2024  // Get coupled neighbour cellRegion
2025  labelList neiCellRegion;
2026  syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
2027 
2028  // Calculate region to zone from cellRegions on either side of coupled
2029  // face.
2030  forAll(patches, patchI)
2031  {
2032  const polyPatch& pp = patches[patchI];
2033 
2034  if (pp.coupled())
2035  {
2036  forAll(pp, i)
2037  {
2038  label faceI = pp.start()+i;
2039 
2040  label regionI = namedSurfaceRegion[faceI];
2041 
2042  // Connected even if no cellZone defined for surface
2043  if (unnamedSurfaceRegion[faceI] == -1 && regionI != -1)
2044  {
2045  label surfI = surfaces_.whichSurface(regionI).first();
2046 
2047  bool changedCell = calcRegionToZone
2048  (
2049  backgroundZoneID,
2050  surfaceToCellZone[surfI],
2051  cellRegion[mesh_.faceOwner()[faceI]],
2052  neiCellRegion[faceI-mesh_.nInternalFaces()],
2053  regionToCellZone
2054  );
2055 
2056  changed = changed | changedCell;
2057  }
2058  }
2059  }
2060  }
2061 
2062  if (!returnReduce(changed, orOp<bool>()))
2063  {
2064  break;
2065  }
2066  }
2067 
2068 
2069  if (debug)
2070  {
2071  Pout<< "meshRefinement::findCellZoneTopo :"
2072  << " nRegions:" << regionToCellZone.size()
2073  << " of which visited (-1 = background, >= 0 : cellZone) :"
2074  << endl;
2075 
2076  forAll(regionToCellZone, regionI)
2077  {
2078  if (regionToCellZone[regionI] != -2)
2079  {
2080  Pout<< "Region " << regionI
2081  << " becomes cellZone:" << regionToCellZone[regionI]
2082  << endl;
2083  }
2084  }
2085  }
2086 
2087  // Rework into cellToZone
2088  forAll(cellToZone, cellI)
2089  {
2090  label regionI = cellRegion[cellI];
2091  if (cellToZone[cellI] == -2 && regionToCellZone[regionI] != -2)
2092  {
2093  cellToZone[cellI] = regionToCellZone[regionI];
2094  }
2095  }
2096 }
2097 
2098 
2099 void Foam::meshRefinement::erodeCellZone
2100 (
2101  const label nErodeCellZones,
2102  const label backgroundZoneID,
2103  const labelList& unnamedSurfaceRegion,
2104  const labelList& namedSurfaceRegion,
2105  labelList& cellToZone
2106 ) const
2107 {
2108  // This routine fixes small problems with left over unassigned regions
2109  // (after all off the unreachable bits of the mesh have been removed).
2110  // The problem is that the cell zone information might be inconsistent
2111  // with the face zone information. So what we do here is to erode
2112  // any cell zones until we hit a named face.
2113  // - backgroundZoneID = -2 : do not change so remove cells
2114  // - backgroundZoneID = -1 : put into background
2115  // Note that is the opposite of findCellZoneTopo which moves unassigned
2116  // regions into a neighbouring region(=cellZone) unless there is an
2117  // intersected faces inbetween the two.
2118 
2119  for (label iter = 0; iter < nErodeCellZones; iter++)
2120  {
2121  label nChanged = 0;
2122 
2123  labelList erodedCellToZone(cellToZone);
2124 
2125  // Do internal faces
2126  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2127  {
2128  if
2129  (
2130  unnamedSurfaceRegion[facei] == -1
2131  && namedSurfaceRegion[facei] == -1
2132  )
2133  {
2134  label own = mesh_.faceOwner()[facei];
2135  label nei = mesh_.faceNeighbour()[facei];
2136  if (cellToZone[own] == -2 && cellToZone[nei] >= -1)
2137  {
2138  erodedCellToZone[nei] = backgroundZoneID;
2139  nChanged++;
2140  }
2141  else if (cellToZone[nei] == -2 && cellToZone[own] >= -1)
2142  {
2143  erodedCellToZone[own] = backgroundZoneID;
2144  nChanged++;
2145  }
2146  }
2147  }
2148 
2149  // Do boundary faces
2150 
2151  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2152 
2153  // Get coupled neighbour cellRegion
2154  labelList neiCellZone;
2155  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2156 
2157  // Calculate region to zone from cellRegions on either side of coupled
2158  // face.
2159  forAll(patches, patchi)
2160  {
2161  const polyPatch& pp = patches[patchi];
2162 
2163  if (pp.coupled())
2164  {
2165  forAll(pp, i)
2166  {
2167  label facei = pp.start()+i;
2168  if
2169  (
2170  unnamedSurfaceRegion[facei] == -1
2171  && namedSurfaceRegion[facei] == -1
2172  )
2173  {
2174  label own = mesh_.faceOwner()[facei];
2175  label bFacei = facei-mesh_.nInternalFaces();
2176  if (neiCellZone[bFacei] == -2 && cellToZone[own] >= -1)
2177  {
2178  erodedCellToZone[own] = backgroundZoneID;
2179  nChanged++;
2180  }
2181  }
2182  }
2183  }
2184  }
2185 
2186  cellToZone.transfer(erodedCellToZone);
2187 
2188  reduce(nChanged, sumOp<label>());
2189  if (debug)
2190  {
2191  Pout<< "erodeCellZone : eroded " << nChanged
2192  << " cells (moved from cellZone to background zone "
2193  << backgroundZoneID << endl;
2194  }
2195 
2196  if (nChanged == 0)
2197  {
2198  break;
2199  }
2200  }
2201 }
2202 
2203 
2204 void Foam::meshRefinement::growCellZone
2205 (
2206  const label nGrowCellZones,
2207  const label backgroundZoneID,
2208  labelList& unnamedSurfaceRegion1,
2209  labelList& unnamedSurfaceRegion2,
2210  labelList& namedSurfaceRegion, // potentially zero size if no faceZones
2211  labelList& cellToZone
2212 ) const
2213 {
2214  if (nGrowCellZones != 1)
2215  {
2216  WarningInFunction
2217  << "Growing currently only supported for single layer."
2218  << " Exiting zone growing" << endl;
2219  return;
2220  }
2221 
2222 
2223  // See meshRefinement::markProximityRefinementWave:
2224  // - walk out nGrowCellZones iterations from outside of cellZone
2225  // and wall into unassigned cells
2226  // - detect any cells inbetween
2227  // - multiple zones
2228  // - zone and wall
2229  // and
2230  // - assign cells to the cellZone
2231  // - unblock faces (along the path) inbetween
2232 
2233 
2234  // Special tag for walls
2235  const FixedList<label, 3> wallTag
2236  ({
2237  labelMax, // Special value for wall face. Loses out to cellZone tag
2238  labelMax,
2239  labelMax
2240  });
2241 
2242  // Work arrays
2243  pointField origins(1);
2244  scalarField distSqrs(1);
2245  List<FixedList<label, 3>> surfaces(1);
2246 
2247 
2248  // Set up blockage. Marked with 0 distance (so always wins)
2249 
2250  //List<wallPoints> allFaceInfo(mesh_.nFaces());
2251  //for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2252  //{
2253  // if
2254  // (
2255  // unnamedSurfaceRegion1[facei] != -1
2256  // || unnamedSurfaceRegion2[facei] != -1
2257  // )
2258  // {
2259  // origins[0] = mesh_.faceCentres()[facei];
2260  // distSqrs[0] = 0.0; // Initial distance
2261  // surfaces[0] = wallTag;
2262  // allFaceInfo[facei] = wallPoints(origins, distSqrs, surfaces);
2263  // }
2264  //}
2265  List<wallPoints> allCellInfo(mesh_.nCells());
2266  forAll(cellToZone, celli)
2267  {
2268  if (cellToZone[celli] >= 0)
2269  {
2270  const FixedList<label, 3> zoneTag
2271  ({
2272  cellToZone[celli], // zone
2273  0, // 'region'
2274  labelMax // 'sub-region'
2275  });
2276 
2277  origins[0] = mesh_.cellCentres()[celli];
2278  distSqrs[0] = 0.0; // Initial distance
2279  surfaces[0] = zoneTag;
2280  allCellInfo[celli] = wallPoints(origins, distSqrs, surfaces);
2281  }
2282  }
2283 
2284 
2285  DynamicList<wallPoints> faceDist(mesh_.nFaces()/128);
2286  DynamicList<label> changedFaces(mesh_.nFaces()/128);
2287 
2288  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
2289  {
2290  const label own = mesh_.faceOwner()[facei];
2291  const label nei = mesh_.faceNeighbour()[facei];
2292  if (cellToZone[own] >= 0 && cellToZone[nei] < 0)
2293  {
2294  // boundary between cellZone (own) and background/unvisited (nei)
2295 
2296  origins[0] = mesh_.faceCentres()[facei];
2297  distSqrs[0] = 0.0; // Initial distance
2298  surfaces[0] = FixedList<label, 3>
2299  ({
2300  cellToZone[own], // zone
2301  0, // 'region'
2302  labelMax // 'sub-region'
2303  });
2304  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2305  changedFaces.append(facei);
2306  }
2307  else if (cellToZone[own] < 0 && cellToZone[nei] >= 0)
2308  {
2309  // boundary between cellZone (nei) and background/unvisited (own)
2310 
2311  origins[0] = mesh_.faceCentres()[facei];
2312  distSqrs[0] = 0.0; // Initial distance
2313  surfaces[0] = FixedList<label, 3>
2314  ({
2315  cellToZone[nei], // zone
2316  0, // 'region'
2317  labelMax // 'sub-region'
2318  });
2319  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2320  changedFaces.append(facei);
2321  }
2322  else if
2323  (
2324  unnamedSurfaceRegion1[facei] != -1
2325  || unnamedSurfaceRegion2[facei] != -1
2326  )
2327  {
2328  // Seed (yet unpatched) wall faces
2329 
2330  origins[0] = mesh_.faceCentres()[facei];
2331  distSqrs[0] = 0.0; // Initial distance
2332  surfaces[0] = wallTag;
2333  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2334  changedFaces.append(facei);
2335  }
2336  }
2337 
2338  // Get coupled neighbour cellRegion
2339  labelList neiCellZone;
2340  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2341 
2342  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2343 
2344  // Calculate region to zone from cellRegions on either side of coupled
2345  // face.
2346  forAll(patches, patchi)
2347  {
2348  const polyPatch& pp = patches[patchi];
2349 
2350  if (pp.coupled())
2351  {
2352  // Like internal faces
2353  forAll(pp, i)
2354  {
2355  label facei = pp.start()+i;
2356  label own = mesh_.faceOwner()[facei];
2357  label bFacei = facei-mesh_.nInternalFaces();
2358  if (cellToZone[own] >= 0 && neiCellZone[bFacei] < 0)
2359  {
2360  origins[0] = mesh_.faceCentres()[facei];
2361  distSqrs[0] = 0.0;
2362  surfaces[0] = FixedList<label, 3>
2363  ({
2364  cellToZone[own], // zone
2365  0, // 'region'
2366  labelMax // 'sub-region'
2367  });
2368  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2369  changedFaces.append(facei);
2370  }
2371  else if (cellToZone[own] < 0 && neiCellZone[bFacei] >= 0)
2372  {
2373  // Handled on nbr processor
2374  }
2375  else if
2376  (
2377  unnamedSurfaceRegion1[facei] != -1
2378  || unnamedSurfaceRegion2[facei] != -1
2379  )
2380  {
2381  // Seed (yet unpatched) wall faces
2382  origins[0] = mesh_.faceCentres()[facei];
2383  distSqrs[0] = 0.0; // Initial distance
2384  surfaces[0] = wallTag;
2385  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2386  changedFaces.append(facei);
2387  }
2388  }
2389  }
2390  else
2391  {
2392  // Seed wall faces
2393  forAll(pp, i)
2394  {
2395  label facei = pp.start()+i;
2396  label own = mesh_.faceOwner()[facei];
2397  if (cellToZone[own] < 0)
2398  {
2399  origins[0] = mesh_.faceCentres()[facei];
2400  distSqrs[0] = 0.0; // Initial distance
2401  surfaces[0] = wallTag;
2402  faceDist.append(wallPoints(origins, distSqrs, surfaces));
2403  changedFaces.append(facei);
2404  }
2405  }
2406  }
2407  }
2408 
2409 
2410  List<wallPoints> allFaceInfo(mesh_.nFaces());
2411  FaceCellWave<wallPoints> wallDistCalc
2412  (
2413  mesh_,
2414  changedFaces,
2415  faceDist,
2416  allFaceInfo,
2417  allCellInfo,
2418  0 // max iterations
2419  );
2420  wallDistCalc.iterate(nGrowCellZones);
2421 
2422 
2423  // Check for cells which are within nGrowCellZones of two cellZones or
2424  // one cellZone and a wall
2425  // TBD. Currently only one layer of growth handled.
2426 
2427  bitSet isChangedCell(mesh_.nCells());
2428 
2429  forAll(allCellInfo, celli)
2430  {
2431  if (allCellInfo[celli].valid(wallDistCalc.data()))
2432  {
2433  const List<FixedList<label, 3>>& surfaces =
2434  allCellInfo[celli].surface();
2435 
2436  if (surfaces.size())
2437  {
2438  // Cell close to cellZone. Remove any free-standing baffles.
2439  // Done by marking as changed cell. Wip.
2440  isChangedCell.set(celli);
2441  }
2442 
2443  if (surfaces.size() > 1)
2444  {
2445  // Check if inbetween two cellZones or cellZone and wall
2446  // Find 'nearest' other cellZone
2447  scalar minDistSqr = Foam::sqr(GREAT);
2448  label minZone = -1;
2449  for (label i = 0; i < surfaces.size(); i++)
2450  {
2451  const label zonei = surfaces[i][0];
2452  const scalar d2 = allCellInfo[celli].distSqr()[i];
2453  if (zonei != cellToZone[celli] && d2 < minDistSqr)
2454  {
2455  minDistSqr = d2;
2456  minZone = zonei; // zoneID
2457  }
2458  }
2459 
2460  if (minDistSqr < Foam::sqr(GREAT))
2461  {
2462  if (minZone != cellToZone[celli] && minZone != wallTag[0])
2463  {
2464  cellToZone[celli] = minZone;
2465  isChangedCell.set(celli);
2466  }
2467  }
2468  }
2469  }
2470  }
2471 
2472  // Fix any left-over unvisited cells
2473  if (backgroundZoneID != -2)
2474  {
2475  forAll(cellToZone, celli)
2476  {
2477  if (cellToZone[celli] == -2)
2478  {
2479  cellToZone[celli] = backgroundZoneID;
2480  }
2481  }
2482  }
2483 
2484 
2485  // Make sure to unset faces of changed cell
2486 
2487  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2488 
2489 
2490  label nUnnamed = 0;
2491  label nNamed = 0;
2492  for (const label celli : isChangedCell)
2493  {
2494  const cell& cFaces = mesh_.cells()[celli];
2495  for (const label facei : cFaces)
2496  {
2497  const label own = mesh_.faceOwner()[facei];
2498  const label ownZone = cellToZone[own];
2499 
2500  label nbrZone = -1;
2501  if (mesh_.isInternalFace(facei))
2502  {
2503  const label neiZone = cellToZone[mesh_.faceNeighbour()[facei]];
2504  nbrZone = (own != celli ? ownZone : neiZone);
2505  }
2506  else
2507  {
2508  nbrZone = neiCellZone[facei-mesh_.nInternalFaces()];
2509  }
2510 
2511  if (nbrZone == cellToZone[celli])
2512  {
2513  if
2514  (
2515  unnamedSurfaceRegion1[facei] != -1
2516  || unnamedSurfaceRegion2[facei] != -1
2517  )
2518  {
2519  unnamedSurfaceRegion1[facei] = -1;
2520  unnamedSurfaceRegion2[facei] = -1;
2521  nUnnamed++;
2522  }
2523  if
2524  (
2525  namedSurfaceRegion.size()
2526  && namedSurfaceRegion[facei] != -1
2527  )
2528  {
2529  namedSurfaceRegion[facei] = -1;
2530  nNamed++;
2531  }
2532  }
2533  }
2534  }
2535 
2536  reduce(nUnnamed, sumOp<label>());
2537  reduce(nNamed, sumOp<label>());
2538 
2539  // Do always; might bypass if nNamed,nUnnamed zero
2540  syncTools::syncFaceList
2541  (
2542  mesh_,
2543  unnamedSurfaceRegion1,
2544  maxEqOp<label>()
2545  );
2546  syncTools::syncFaceList
2547  (
2548  mesh_,
2549  unnamedSurfaceRegion2,
2550  maxEqOp<label>()
2551  );
2552  if (namedSurfaceRegion.size())
2553  {
2554  syncTools::syncFaceList
2555  (
2556  mesh_,
2557  namedSurfaceRegion,
2558  maxEqOp<label>()
2559  );
2560  }
2561 
2562  if (debug)
2563  {
2564  Pout<< "growCellZone : grown cellZones by "
2565  << returnReduce(isChangedCell.count(), sumOp<label>())
2566  << " cells (moved from background to nearest cellZone)"
2567  << endl;
2568  Pout<< "growCellZone : unmarked " << nUnnamed
2569  << " unzoned intersections; " << nNamed << " zoned intersections; "
2570  << endl;
2571  }
2572 }
2573 
2574 
2575 void Foam::meshRefinement::makeConsistentFaceIndex
2576 (
2577  const labelList& surfaceMap,
2578  const labelList& cellToZone,
2579  labelList& namedSurfaceRegion
2580 ) const
2581 {
2582  // Make namedSurfaceRegion consistent with cellToZone - clear out any
2583  // blocked faces inbetween same cell zone (or background (=-1))
2584  // Do not do this for surfaces relating to 'pure' faceZones i.e.
2585  // faceZones without a cellZone. Note that we cannot check here
2586  // for different cellZones on either side but no namedSurfaceRegion
2587  // since cellZones can now originate from locationsInMesh as well
2588  // (instead of only through named surfaces)
2589 
2590  const labelList& faceOwner = mesh_.faceOwner();
2591  const labelList& faceNeighbour = mesh_.faceNeighbour();
2592 
2593  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2594  {
2595  label ownZone = cellToZone[faceOwner[faceI]];
2596  label neiZone = cellToZone[faceNeighbour[faceI]];
2597  label globalI = namedSurfaceRegion[faceI];
2598 
2599  if
2600  (
2601  ownZone == neiZone
2602  && globalI != -1
2603  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2604  )
2605  {
2606  namedSurfaceRegion[faceI] = -1;
2607  }
2608  }
2609 
2610  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2611 
2612  // Get coupled neighbour cellZone
2613  labelList neiCellZone;
2614  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2615 
2616  // Use coupled cellZone to do check
2617  forAll(patches, patchI)
2618  {
2619  const polyPatch& pp = patches[patchI];
2620 
2621  if (pp.coupled())
2622  {
2623  forAll(pp, i)
2624  {
2625  label faceI = pp.start()+i;
2626 
2627  label ownZone = cellToZone[faceOwner[faceI]];
2628  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
2629  label globalI = namedSurfaceRegion[faceI];
2630 
2631  if
2632  (
2633  ownZone == neiZone
2634  && globalI != -1
2635  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2636  )
2637  {
2638  namedSurfaceRegion[faceI] = -1;
2639  }
2640  }
2641  }
2642  else
2643  {
2644  // Unzonify boundary faces
2645  forAll(pp, i)
2646  {
2647  label faceI = pp.start()+i;
2648  label globalI = namedSurfaceRegion[faceI];
2649 
2650  if
2651  (
2652  globalI != -1
2653  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2654  )
2655  {
2656  namedSurfaceRegion[faceI] = -1;
2657  }
2658  }
2659  }
2660  }
2661 }
2662 
2663 
2664 void Foam::meshRefinement::getIntersections
2665 (
2666  const labelList& surfacesToTest,
2667  const pointField& neiCc,
2668  const labelList& testFaces,
2669 
2670  labelList& namedSurfaceRegion,
2671  bitSet& posOrientation
2672 ) const
2673 {
2674  namedSurfaceRegion.setSize(mesh_.nFaces());
2675  namedSurfaceRegion = -1;
2676 
2677  posOrientation.setSize(mesh_.nFaces());
2678  posOrientation = false;
2679 
2680  // Statistics: number of faces per faceZone
2681  labelList nSurfFaces(surfaces_.surfZones().size(), Zero);
2682 
2683  // Collect segments
2684  // ~~~~~~~~~~~~~~~~
2685 
2686  pointField start(testFaces.size());
2687  pointField end(testFaces.size());
2688 
2689  {
2690  labelList minLevel;
2691  calcCellCellRays
2692  (
2693  neiCc,
2694  labelList(neiCc.size(), -1),
2695  testFaces,
2696  start,
2697  end,
2698  minLevel
2699  );
2700  }
2701 
2702  // Do test for intersections
2703  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2704  // Note that we intersect all intersected faces again. Could reuse
2705  // the information already in surfaceIndex_.
2706 
2707  labelList surface1;
2708  labelList region1;
2709  List<pointIndexHit> hit1;
2710  vectorField normal1;
2711  labelList surface2;
2712  labelList region2;
2713  List<pointIndexHit> hit2;
2714  vectorField normal2;
2715 
2716  surfaces_.findNearestIntersection
2717  (
2718  surfacesToTest,
2719  start,
2720  end,
2721 
2722  surface1,
2723  hit1,
2724  region1,
2725  normal1,
2726 
2727  surface2,
2728  hit2,
2729  region2,
2730  normal2
2731  );
2732 
2733  forAll(testFaces, i)
2734  {
2735  label faceI = testFaces[i];
2736  const vector& area = mesh_.faceAreas()[faceI];
2737 
2738  if (surface1[i] != -1)
2739  {
2740  // If both hit should probably choose 'nearest'
2741  if
2742  (
2743  surface2[i] != -1
2744  && (
2745  magSqr(hit2[i].hitPoint())
2746  < magSqr(hit1[i].hitPoint())
2747  )
2748  )
2749  {
2750  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2751  (
2752  surface2[i],
2753  region2[i]
2754  );
2755  posOrientation.set(faceI, ((area&normal2[i]) > 0));
2756  nSurfFaces[surface2[i]]++;
2757  }
2758  else
2759  {
2760  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2761  (
2762  surface1[i],
2763  region1[i]
2764  );
2765  posOrientation.set(faceI, ((area&normal1[i]) > 0));
2766  nSurfFaces[surface1[i]]++;
2767  }
2768  }
2769  else if (surface2[i] != -1)
2770  {
2771  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2772  (
2773  surface2[i],
2774  region2[i]
2775  );
2776  posOrientation.set(faceI, ((area&normal2[i]) > 0));
2777  nSurfFaces[surface2[i]]++;
2778  }
2779  }
2780 
2781 
2782  // surfaceIndex might have different surfaces on both sides if
2783  // there happen to be a (obviously thin) surface with different
2784  // regions between the cell centres. If one is on a named surface
2785  // and the other is not this might give problems so sync.
2786  syncTools::syncFaceList
2787  (
2788  mesh_,
2789  namedSurfaceRegion,
2790  maxEqOp<label>()
2791  );
2792 
2793  // Print a bit
2794  if (debug)
2795  {
2796  forAll(nSurfFaces, surfI)
2797  {
2798  Pout<< "Surface:"
2799  << surfaces_.names()[surfI]
2800  << " nZoneFaces:" << nSurfFaces[surfI] << nl;
2801  }
2802  Pout<< endl;
2803  }
2804 }
2805 
2806 
2807 void Foam::meshRefinement::zonify
2808 (
2809  const bool allowFreeStandingZoneFaces,
2810  const label nErodeCellZones,
2811  const label backgroundZoneID,
2812  const pointField& locationsInMesh,
2813  const wordList& zonesInMesh,
2814 
2815  labelList& cellToZone,
2816  labelList& unnamedRegion1,
2817  labelList& unnamedRegion2,
2818  labelList& namedSurfaceRegion,
2819  bitSet& posOrientation
2820 ) const
2821 {
2822  // Determine zones for cells and faces
2823  // cellToZone:
2824  // -2 : unset
2825  // -1 : not in any zone (zone 'none' or background zone)
2826  // >=0 : zoneID
2827  // namedSurfaceRegion, posOrientation:
2828  // -1 : face not intersected by named surface
2829  // >=0 : globalRegion (surface+region)
2830  // (and posOrientation: surface normal v.s. face normal)
2831 
2832  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
2833 
2834  // Swap neighbouring cell centres and cell level
2835  labelList neiLevel(mesh_.nBoundaryFaces());
2836  pointField neiCc(mesh_.nBoundaryFaces());
2837  calcNeighbourData(neiLevel, neiCc);
2838 
2839 
2840  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
2841  labelList unnamedSurfaces(surfaceZonesInfo::getUnnamedSurfaces(surfZones));
2842 
2843  // Get map from surface to cellZone (or -1)
2844  labelList surfaceToCellZone;
2845  if (namedSurfaces.size())
2846  {
2847  // Get/add cellZones corresponding to surface names
2848  surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
2849  (
2850  surfZones,
2851  namedSurfaces,
2852  mesh_
2853  );
2854  }
2855 
2856 
2857  cellToZone.setSize(mesh_.nCells());
2858  cellToZone = -2;
2859 
2860  namedSurfaceRegion.clear();
2861  posOrientation.clear();
2862 
2863 
2864 
2865  // 1. Test all (unnamed & named) surfaces
2866 
2867  // Unnamed surfaces. Global surface region of intersection (or -1)
2868  getIntersections
2869  (
2870  unnamedSurfaces,
2871  neiCc,
2872  intersectedFaces(),
2873  unnamedRegion1,
2874  unnamedRegion2
2875  );
2876 
2877 
2878  if (namedSurfaces.size())
2879  {
2880  getIntersections
2881  (
2882  namedSurfaces,
2883  neiCc,
2884  intersectedFaces(),
2885  namedSurfaceRegion,
2886  posOrientation
2887  );
2888  }
2889 
2890 
2891  // 2. Walk from locationsInMesh. Hard set cellZones.
2892  // Note: walk through faceZones! (these might get overridden later)
2893 
2894  if (locationsInMesh.size())
2895  {
2896  Info<< "Setting cellZones according to locationsInMesh:" << endl;
2897 
2898  labelList locationsZoneIDs(zonesInMesh.size(), -1);
2899  forAll(locationsInMesh, i)
2900  {
2901  const word& name = zonesInMesh[i];
2902 
2903  Info<< "Location : " << locationsInMesh[i] << nl
2904  << " cellZone : " << name << endl;
2905 
2906  if (name != "none")
2907  {
2908  label zoneID = mesh_.cellZones().findZoneID(name);
2909  if (zoneID == -1)
2910  {
2911  FatalErrorInFunction << "problem" << abort(FatalError);
2912  }
2913  locationsZoneIDs[i] = zoneID;
2914  }
2915  }
2916  Info<< endl;
2917 
2918 
2919  // Assign cellZone according to seed points. Walk through named surfaces
2920  findCellZoneInsideWalk
2921  (
2922  locationsInMesh, // locations
2923  locationsZoneIDs, // index of cellZone (or -1)
2924  unnamedRegion1, // per face -1 (unblocked) or >= 0 (blocked)
2925  cellToZone
2926  );
2927  }
2928 
2929 
2930  // 3. Mark named-surfaces-with-insidePoint. Hard set cellZones.
2931 
2932  labelList locationSurfaces
2933  (
2934  surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
2935  );
2936 
2937  if (locationSurfaces.size())
2938  {
2939  Info<< "Found " << locationSurfaces.size()
2940  << " named surfaces with a provided inside point."
2941  << " Assigning cells inside these surfaces"
2942  << " to the corresponding cellZone."
2943  << nl << endl;
2944 
2945  // Collect per surface the -insidePoint -cellZone
2946  pointField insidePoints(locationSurfaces.size());
2947  labelList insidePointCellZoneIDs(locationSurfaces.size(), -1);
2948  forAll(locationSurfaces, i)
2949  {
2950  label surfI = locationSurfaces[i];
2951  insidePoints[i] = surfZones[surfI].zoneInsidePoint();
2952 
2953  const word& name = surfZones[surfI].cellZoneName();
2954  if (name != "none")
2955  {
2956  label zoneID = mesh_.cellZones().findZoneID(name);
2957  if (zoneID == -1)
2958  {
2959  FatalErrorInFunction
2960  << "problem"
2961  << abort(FatalError);
2962  }
2963  insidePointCellZoneIDs[i] = zoneID;
2964  }
2965  }
2966 
2967 
2968  // Stop at unnamed or named surface
2969  labelList allRegion1(mesh_.nFaces(), -1);
2970  forAll(namedSurfaceRegion, faceI)
2971  {
2972  allRegion1[faceI] = max
2973  (
2974  unnamedRegion1[faceI],
2975  namedSurfaceRegion[faceI]
2976  );
2977  }
2978 
2979  findCellZoneInsideWalk
2980  (
2981  insidePoints, // locations
2982  insidePointCellZoneIDs, // index of cellZone
2983  allRegion1, // per face -1 (unblocked) or >= 0 (blocked)
2984  cellToZone
2985  );
2986  }
2987 
2988 
2989  // 4. Mark named-surfaces-with-geometric faces. Do geometric test. Soft set
2990  // cellZones. Correct through making consistent.
2991 
2992  // Closed surfaces with cellZone specified.
2993  labelList closedNamedSurfaces
2994  (
2995  surfaceZonesInfo::getClosedNamedSurfaces
2996  (
2997  surfZones,
2998  surfaces_.geometry(),
2999  surfaces_.surfaces()
3000  )
3001  );
3002 
3003  if (closedNamedSurfaces.size())
3004  {
3005  Info<< "Found " << closedNamedSurfaces.size()
3006  << " closed, named surfaces. Assigning cells in/outside"
3007  << " these surfaces to the corresponding cellZone."
3008  << nl << endl;
3009 
3010  findCellZoneGeometric
3011  (
3012  neiCc,
3013  closedNamedSurfaces, // indices of closed surfaces
3014  namedSurfaceRegion, // per face index of named surface + region
3015  surfaceToCellZone, // cell zone index per surface
3016 
3017  cellToZone
3018  );
3019  }
3020 
3021 
3022  // 5. Find any unassigned regions (from regionSplit)
3023 
3024  if (namedSurfaces.size())
3025  {
3026  if (nErodeCellZones == 0)
3027  {
3028  Info<< "Walking from known cellZones; crossing a faceZone "
3029  << "face changes cellZone" << nl << endl;
3030 
3031  // Put unassigned regions into any connected cellZone
3032  findCellZoneTopo
3033  (
3034  backgroundZoneID,
3035  pointField(0),
3036  unnamedRegion1, // Intersections with unnamed surfaces
3037  namedSurfaceRegion, // Intersections with named surfaces
3038  surfaceToCellZone,
3039  cellToZone
3040  );
3041  }
3042  else if (nErodeCellZones < 0)
3043  {
3044  Info<< "Growing cellZones by " << -nErodeCellZones
3045  << " layers" << nl << endl;
3046 
3047  growCellZone
3048  (
3049  -nErodeCellZones,
3050  backgroundZoneID,
3051  unnamedRegion1,
3052  unnamedRegion2,
3053  namedSurfaceRegion,
3054  cellToZone
3055  );
3056  }
3057  else
3058  {
3059  Info<< "Eroding cellZone cells to make these consistent with"
3060  << " faceZone faces" << nl << endl;
3061 
3062  // Erode cell zone cells (connected to an unassigned cell)
3063  // and put them into backgroundZone
3064  erodeCellZone
3065  (
3066  nErodeCellZones,
3067  backgroundZoneID,
3068  unnamedRegion1,
3069  namedSurfaceRegion,
3070  cellToZone
3071  );
3072  }
3073 
3074 
3075  // Make sure namedSurfaceRegion is unset inbetween same cell zones.
3076  if (!allowFreeStandingZoneFaces)
3077  {
3078  Info<< "Only keeping zone faces inbetween different cellZones."
3079  << nl << endl;
3080 
3081  // Surfaces with faceZone but no cellZone
3082  labelList standaloneNamedSurfaces
3083  (
3084  surfaceZonesInfo::getStandaloneNamedSurfaces
3085  (
3086  surfZones
3087  )
3088  );
3089 
3090  // Construct map from surface index to index in
3091  // standaloneNamedSurfaces (or -1)
3092  labelList surfaceMap(surfZones.size(), -1);
3093  forAll(standaloneNamedSurfaces, i)
3094  {
3095  surfaceMap[standaloneNamedSurfaces[i]] = i;
3096  }
3097 
3098  makeConsistentFaceIndex
3099  (
3100  surfaceMap,
3101  cellToZone,
3102  namedSurfaceRegion
3103  );
3104  }
3105  }
3106  else if (nErodeCellZones < 0 && gMax(cellToZone) >= 0)
3107  {
3108  // With multiple locationsInMesh and non-trivial cellZones we allow
3109  // some growing of the cellZones to avoid any background cells
3110 
3111  Info<< "Growing cellZones by " << -nErodeCellZones
3112  << " layers" << nl << endl;
3113 
3114  growCellZone
3115  (
3116  -nErodeCellZones,
3117  backgroundZoneID,
3118  unnamedRegion1,
3119  unnamedRegion2,
3120  namedSurfaceRegion, // note: potentially zero sized
3121  cellToZone
3122  );
3123 
3124  // Make sure namedSurfaceRegion is unset inbetween same cell zones.
3125  if (!allowFreeStandingZoneFaces && namedSurfaceRegion.size())
3126  {
3127  Info<< "Only keeping zone faces inbetween different cellZones."
3128  << nl << endl;
3129 
3130  // Surfaces with faceZone but no cellZone
3131  labelList standaloneNamedSurfaces
3132  (
3133  surfaceZonesInfo::getStandaloneNamedSurfaces
3134  (
3135  surfZones
3136  )
3137  );
3138 
3139  // Construct map from surface index to index in
3140  // standaloneNamedSurfaces (or -1)
3141  labelList surfaceMap(surfZones.size(), -1);
3142  forAll(standaloneNamedSurfaces, i)
3143  {
3144  surfaceMap[standaloneNamedSurfaces[i]] = i;
3145  }
3146 
3147  makeConsistentFaceIndex
3148  (
3149  surfaceMap,
3150  cellToZone,
3151  namedSurfaceRegion
3152  );
3153  }
3154  }
3155 
3156 
3157  // Some stats
3158  if (debug)
3159  {
3160  label nZones = gMax(cellToZone)+1;
3161 
3162  label nUnvisited = 0;
3163  label nBackgroundCells = 0;
3164  labelList nZoneCells(nZones, Zero);
3165  forAll(cellToZone, cellI)
3166  {
3167  label zoneI = cellToZone[cellI];
3168  if (zoneI >= 0)
3169  {
3170  nZoneCells[zoneI]++;
3171  }
3172  else if (zoneI == -1)
3173  {
3174  nBackgroundCells++;
3175  }
3176  else if (zoneI == -2)
3177  {
3178  nUnvisited++;
3179  }
3180  else
3181  {
3182  FatalErrorInFunction
3183  << "problem" << exit(FatalError);
3184  }
3185  }
3186  reduce(nUnvisited, sumOp<label>());
3187  reduce(nBackgroundCells, sumOp<label>());
3188  forAll(nZoneCells, zoneI)
3189  {
3190  reduce(nZoneCells[zoneI], sumOp<label>());
3191  }
3192  Info<< "nUnvisited :" << nUnvisited << endl;
3193  Info<< "nBackgroundCells:" << nBackgroundCells << endl;
3194  Info<< "nZoneCells :" << nZoneCells << endl;
3195  }
3196  if (debug&MESH)
3197  {
3198  const_cast<Time&>(mesh_.time())++;
3199  Pout<< "Writing cell zone allocation on mesh to time "
3200  << timeName() << endl;
3201  write
3202  (
3203  debugType(debug),
3204  writeType(writeLevel() | WRITEMESH),
3205  mesh_.time().path()/"cell2Zone"
3206  );
3207  volScalarField volCellToZone
3208  (
3209  IOobject
3210  (
3211  "cellToZone",
3212  timeName(),
3213  mesh_,
3214  IOobject::NO_READ,
3215  IOobject::AUTO_WRITE,
3216  false
3217  ),
3218  mesh_,
3219  dimensionedScalar(dimless, Zero),
3220  zeroGradientFvPatchScalarField::typeName
3221  );
3222 
3223  forAll(cellToZone, cellI)
3224  {
3225  volCellToZone[cellI] = cellToZone[cellI];
3226  }
3227  volCellToZone.write();
3228 
3229 
3230  //mkDir(mesh_.time().path()/timeName());
3231  //OBJstream str
3232  //(
3233  // mesh_.time().path()/timeName()/"zoneBoundaryFaces.obj"
3234  //);
3235  //Pout<< "Writing zone boundaries to " << str.name() << endl;
3236  //for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3237  //{
3238  // const label ownZone = cellToZone[mesh_.faceOwner()[facei]];
3239  // const label neiZone = cellToZone[mesh_.faceNeighbour()[facei]];
3240  // if (ownZone != neiZone)
3241  // {
3242  // str.write(mesh_.faces()[facei], mesh_.points(), false);
3243  // }
3244  //}
3245  //labelList neiCellZone;
3246  //syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
3247  //for
3248  //(
3249  // label facei = mesh_.nInternalFaces();
3250  // facei < mesh_.nFaces();
3251  // facei++
3252  //)
3253  //{
3254  // const label ownZone = cellToZone[mesh_.faceOwner()[facei]];
3255  // const label bFacei = facei-mesh_.nInternalFaces();
3256  // const label neiZone = neiCellZone[bFacei];
3257  // if (ownZone != neiZone)
3258  // {
3259  // str.write(mesh_.faces()[facei], mesh_.points(), false);
3260  // }
3261  //}
3262 
3263  //mkDir(mesh_.time().path()/timeName());
3264  //OBJstream str1
3265  //(
3266  // mesh_.time().path()/timeName()/"unnamedRegion1.obj"
3267  //);
3268  //OBJstream str2
3269  //(
3270  // mesh_.time().path()/timeName()/"unnamedRegion2.obj"
3271  //);
3272  //Pout<< "Writing unnamed1 to " << str1.name() << endl;
3273  //Pout<< "Writing unnamed2 to " << str2.name() << endl;
3274  //for (label facei = 0; facei < mesh_.nFaces(); facei++)
3275  //{
3276  // if
3277  // (
3278  // unnamedRegion1[facei] < -1
3279  // || unnamedRegion2[facei] < -1
3280  // )
3281  // {
3282  // FatalErrorInFunction << "face:"
3283  // << mesh_.faceCentres()[facei]
3284  // << " unnamed1:" << unnamedRegion1[facei]
3285  // << " unnamed2:" << unnamedRegion2[facei]
3286  // << exit(FatalError);
3287  // }
3288  //
3289  // if (unnamedRegion1[facei] >= 0)
3290  // {
3291  // str1.write(mesh_.faces()[facei], mesh_.points(), false);
3292  // }
3293  //
3294  // if (unnamedRegion2[facei] >= 0)
3295  // {
3296  // str2.write(mesh_.faces()[facei], mesh_.points(), false);
3297  // }
3298  //}
3299 
3300  //if (namedSurfaceRegion.size())
3301  //{
3302  // OBJstream strNamed
3303  // (
3304  // mesh_.time().path()/timeName()/"namedSurfaceRegion.obj"
3305  // );
3306  // Pout<< "Writing named to " << strNamed.name() << endl;
3307  // for (label facei = 0; facei < mesh_.nFaces(); facei++)
3308  // {
3309  // const face& f = mesh_.faces()[facei];
3310  // if (namedSurfaceRegion[facei] < -1)
3311  // {
3312  // FatalErrorInFunction << "face:"
3313  // << mesh_.faceCentres()[facei]
3314  // << " unnamed1:" << unnamedRegion1[facei]
3315  // << " unnamed2:" << unnamedRegion2[facei]
3316  // << " named:" << namedSurfaceRegion[facei]
3317  // << exit(FatalError);
3318  // }
3319  // if (namedSurfaceRegion[facei] >= 0)
3320  // {
3321  // strNamed.write(f, mesh_.points(), false);
3322  // }
3323  // }
3324  //}
3325  }
3326 }
3327 
3328 
3329 void Foam::meshRefinement::handleSnapProblems
3330 (
3331  const snapParameters& snapParams,
3332  const bool useTopologicalSnapDetection,
3333  const bool removeEdgeConnectedCells,
3334  const scalarField& perpendicularAngle,
3335  const dictionary& motionDict,
3336  Time& runTime,
3337  const labelList& globalToMasterPatch,
3338  const labelList& globalToSlavePatch
3339 )
3340 {
3341  Info<< nl
3342  << "Introducing baffles to block off problem cells" << nl
3343  << "----------------------------------------------" << nl
3344  << endl;
3345 
3346  labelList facePatch;
3347  if (useTopologicalSnapDetection)
3348  {
3349  facePatch = markFacesOnProblemCells
3350  (
3351  motionDict,
3352  removeEdgeConnectedCells,
3353  perpendicularAngle,
3354  globalToMasterPatch
3355  );
3356  }
3357  else
3358  {
3359  facePatch = markFacesOnProblemCellsGeometric
3360  (
3361  snapParams,
3362  motionDict,
3363  globalToMasterPatch,
3364  globalToSlavePatch
3365  );
3366  }
3367  Info<< "Analyzed problem cells in = "
3368  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3369 
3370  if (debug&MESH)
3371  {
3372  faceSet problemFaces(mesh_, "problemFaces", mesh_.nFaces()/100);
3373 
3374  forAll(facePatch, faceI)
3375  {
3376  if (facePatch[faceI] != -1)
3377  {
3378  problemFaces.insert(faceI);
3379  }
3380  }
3381  problemFaces.instance() = timeName();
3382  Pout<< "Dumping " << problemFaces.size()
3383  << " problem faces to " << problemFaces.objectPath() << endl;
3384  problemFaces.write();
3385  }
3386 
3387  Info<< "Introducing baffles to delete problem cells." << nl << endl;
3388 
3389  if (debug)
3390  {
3391  ++runTime;
3392  }
3393 
3394  // Create baffles with same owner and neighbour for now.
3395  createBaffles(facePatch, facePatch);
3396 
3397  if (debug)
3398  {
3399  // Debug:test all is still synced across proc patches
3400  checkData();
3401  }
3402  Info<< "Created baffles in = "
3403  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3404 
3405  printMeshInfo(debug, "After introducing baffles");
3406 
3407  if (debug&MESH)
3408  {
3409  const_cast<Time&>(mesh_.time())++;
3410  Pout<< "Writing extra baffled mesh to time "
3411  << timeName() << endl;
3412  write
3413  (
3414  debugType(debug),
3415  writeType(writeLevel() | WRITEMESH),
3416  runTime.path()/"extraBaffles"
3417  );
3418  Pout<< "Dumped debug data in = "
3419  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3420  }
3421 }
3422 
3423 
3424 Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
3425 (
3426  const labelList& faceToZone,
3427  const labelList& cellToZone,
3428  const labelList& neiCellZone
3429 ) const
3430 {
3431  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
3432  const labelList& faceOwner = mesh_.faceOwner();
3433  const labelList& faceNeighbour = mesh_.faceNeighbour();
3434 
3435 
3436  // We want to pick up the faces to orient. These faces come in
3437  // two variants:
3438  // - faces originating from stand-alone faceZones
3439  // (these will most likely have no cellZone on either side so
3440  // ownZone and neiZone both -1)
3441  // - sticky-up faces originating from a 'bulge' in a outside of
3442  // a cellZone. These will have the same cellZone on either side.
3443  // How to orient these is not really clearly defined so do them
3444  // same as stand-alone faceZone faces for now. (Normally these will
3445  // already have been removed by the 'allowFreeStandingZoneFaces=false'
3446  // default setting)
3447 
3448  // Note that argument neiCellZone will have -1 on uncoupled boundaries.
3449 
3450  DynamicList<label> faceLabels(mesh_.nFaces()/100);
3451 
3452  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3453  {
3454  if (faceToZone[faceI] != -1)
3455  {
3456  // Free standing baffle?
3457  label ownZone = cellToZone[faceOwner[faceI]];
3458  label neiZone = cellToZone[faceNeighbour[faceI]];
3459  if (ownZone == neiZone)
3460  {
3461  faceLabels.append(faceI);
3462  }
3463  }
3464  }
3465  forAll(patches, patchI)
3466  {
3467  const polyPatch& pp = patches[patchI];
3468 
3469  forAll(pp, i)
3470  {
3471  label faceI = pp.start()+i;
3472  if (faceToZone[faceI] != -1)
3473  {
3474  // Free standing baffle?
3475  label ownZone = cellToZone[faceOwner[faceI]];
3476  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
3477  if (ownZone == neiZone)
3478  {
3479  faceLabels.append(faceI);
3480  }
3481  }
3482  }
3483  }
3484  return faceLabels.shrink();
3485 }
3486 
3487 
3488 void Foam::meshRefinement::calcPatchNumMasterFaces
3489 (
3490  const bitSet& isMasterFace,
3491  const indirectPrimitivePatch& patch,
3492  labelList& nMasterFacesPerEdge
3493 ) const
3494 {
3495  // Number of (master)faces per edge
3496  nMasterFacesPerEdge.setSize(patch.nEdges());
3497  nMasterFacesPerEdge = 0;
3498 
3499  forAll(patch.addressing(), faceI)
3500  {
3501  const label meshFaceI = patch.addressing()[faceI];
3502 
3503  if (isMasterFace.test(meshFaceI))
3504  {
3505  const labelList& fEdges = patch.faceEdges()[faceI];
3506  forAll(fEdges, fEdgeI)
3507  {
3508  nMasterFacesPerEdge[fEdges[fEdgeI]]++;
3509  }
3510  }
3511  }
3512 
3513  syncTools::syncEdgeList
3514  (
3515  mesh_,
3516  patch.meshEdges(mesh_.edges(), mesh_.pointEdges()),
3517  nMasterFacesPerEdge,
3518  plusEqOp<label>(),
3519  label(0)
3520  );
3521 }
3522 
3523 
3524 Foam::label Foam::meshRefinement::markPatchZones
3525 (
3526  const indirectPrimitivePatch& patch,
3527  const labelList& nMasterFacesPerEdge,
3528  labelList& faceToZone
3529 ) const
3530 {
3531  List<edgeTopoDistanceData<label>> allEdgeInfo(patch.nEdges());
3532  List<edgeTopoDistanceData<label>> allFaceInfo(patch.size());
3533 
3534 
3535  // Protect all non-manifold edges
3536  {
3537  label nProtected = 0;
3538 
3539  forAll(nMasterFacesPerEdge, edgeI)
3540  {
3541  if (nMasterFacesPerEdge[edgeI] > 2)
3542  {
3543  allEdgeInfo[edgeI] = edgeTopoDistanceData<label>(0, -2);
3544  nProtected++;
3545  }
3546  }
3547  //Info<< "Protected from visiting "
3548  // << returnReduce(nProtected, sumOp<label>())
3549  // << " non-manifold edges" << nl << endl;
3550  }
3551 
3552 
3553  // Hand out zones
3554 
3555  DynamicList<label> changedEdges;
3556  DynamicList<edgeTopoDistanceData<label>> changedInfo;
3557 
3558  const scalar tol = PatchEdgeFaceWave
3559  <
3560  indirectPrimitivePatch,
3561  edgeTopoDistanceData<label>
3562  >::propagationTol();
3563 
3564  int dummyTrackData;
3565 
3566  const globalIndex globalFaces(patch.size());
3567 
3568  label faceI = 0;
3569 
3570  label currentZoneI = 0;
3571 
3572  while (true)
3573  {
3574  // Pick an unset face
3575  label globalSeed = labelMax;
3576  for (; faceI < allFaceInfo.size(); faceI++)
3577  {
3578  if (!allFaceInfo[faceI].valid(dummyTrackData))
3579  {
3580  globalSeed = globalFaces.toGlobal(faceI);
3581  break;
3582  }
3583  }
3584 
3585  reduce(globalSeed, minOp<label>());
3586 
3587  if (globalSeed == labelMax)
3588  {
3589  break;
3590  }
3591 
3592  label procI = globalFaces.whichProcID(globalSeed);
3593  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
3594 
3595  //Info<< "Seeding zone " << currentZoneI
3596  // << " from processor " << procI << " face " << seedFaceI
3597  // << endl;
3598 
3599  if (procI == Pstream::myProcNo())
3600  {
3601  edgeTopoDistanceData<label>& faceInfo = allFaceInfo[seedFaceI];
3602 
3603 
3604  // Set face
3605  faceInfo = edgeTopoDistanceData<label>(0, currentZoneI);
3606 
3607  // .. and seed its edges
3608  const labelList& fEdges = patch.faceEdges()[seedFaceI];
3609  forAll(fEdges, fEdgeI)
3610  {
3611  label edgeI = fEdges[fEdgeI];
3612 
3613  edgeTopoDistanceData<label>& edgeInfo = allEdgeInfo[edgeI];
3614 
3615  if
3616  (
3617  edgeInfo.updateEdge<int>
3618  (
3619  mesh_,
3620  patch,
3621  edgeI,
3622  seedFaceI,
3623  faceInfo,
3624  tol,
3625  dummyTrackData
3626  )
3627  )
3628  {
3629  changedEdges.append(edgeI);
3630  changedInfo.append(edgeInfo);
3631  }
3632  }
3633  }
3634 
3635 
3636  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
3637  {
3638  break;
3639  }
3640 
3641 
3642  // Walk
3643  PatchEdgeFaceWave
3644  <
3645  indirectPrimitivePatch,
3646  edgeTopoDistanceData<label>
3647  > calc
3648  (
3649  mesh_,
3650  patch,
3651  changedEdges,
3652  changedInfo,
3653  allEdgeInfo,
3654  allFaceInfo,
3655  returnReduce(patch.nEdges(), sumOp<label>())
3656  );
3657 
3658  currentZoneI++;
3659  }
3660 
3661 
3662  faceToZone.setSize(patch.size());
3663  forAll(allFaceInfo, faceI)
3664  {
3665  if (!allFaceInfo[faceI].valid(dummyTrackData))
3666  {
3667  FatalErrorInFunction
3668  << "Problem: unvisited face " << faceI
3669  << " at " << patch.faceCentres()[faceI]
3670  << exit(FatalError);
3671  }
3672  faceToZone[faceI] = allFaceInfo[faceI].data();
3673  }
3674 
3675  return currentZoneI;
3676 }
3677 
3678 
3679 void Foam::meshRefinement::consistentOrientation
3680 (
3681  const bitSet& isMasterFace,
3682  const indirectPrimitivePatch& patch,
3683  const labelList& nMasterFacesPerEdge,
3684  const labelList& faceToZone,
3685  const Map<label>& zoneToOrientation,
3686  bitSet& meshFlipMap
3687 ) const
3688 {
3689  const polyBoundaryMesh& bm = mesh_.boundaryMesh();
3690 
3691  // Data on all edges and faces
3692  List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
3693  List<patchFaceOrientation> allFaceInfo(patch.size());
3694 
3695  // Make sure we don't walk through
3696  // - slaves of coupled faces
3697  // - non-manifold edges
3698  {
3699  label nProtected = 0;
3700 
3701  forAll(patch.addressing(), faceI)
3702  {
3703  const label meshFaceI = patch.addressing()[faceI];
3704  const label patchI = bm.whichPatch(meshFaceI);
3705 
3706  if
3707  (
3708  patchI != -1
3709  && bm[patchI].coupled()
3710  && !isMasterFace.test(meshFaceI)
3711  )
3712  {
3713  // Slave side. Mark so doesn't get visited.
3714  allFaceInfo[faceI] = orientedSurface::NOFLIP;
3715  nProtected++;
3716  }
3717  }
3718  //Info<< "Protected from visiting "
3719  // << returnReduce(nProtected, sumOp<label>())
3720  // << " slaves of coupled faces" << nl << endl;
3721  }
3722  {
3723  label nProtected = 0;
3724 
3725  forAll(nMasterFacesPerEdge, edgeI)
3726  {
3727  if (nMasterFacesPerEdge[edgeI] > 2)
3728  {
3729  allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
3730  nProtected++;
3731  }
3732  }
3733 
3734  Info<< "Protected from visiting "
3735  << returnReduce(nProtected, sumOp<label>())
3736  << " non-manifold edges" << nl << endl;
3737  }
3738 
3739 
3740 
3741  DynamicList<label> changedEdges;
3742  DynamicList<patchFaceOrientation> changedInfo;
3743 
3744  const scalar tol = PatchEdgeFaceWave
3745  <
3746  indirectPrimitivePatch,
3747  patchFaceOrientation
3748  >::propagationTol();
3749 
3750  int dummyTrackData;
3751 
3752  globalIndex globalFaces(patch.size());
3753 
3754  while (true)
3755  {
3756  // Pick an unset face
3757  label globalSeed = labelMax;
3758  forAll(allFaceInfo, faceI)
3759  {
3760  if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
3761  {
3762  globalSeed = globalFaces.toGlobal(faceI);
3763  break;
3764  }
3765  }
3766 
3767  reduce(globalSeed, minOp<label>());
3768 
3769  if (globalSeed == labelMax)
3770  {
3771  break;
3772  }
3773 
3774  label procI = globalFaces.whichProcID(globalSeed);
3775  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
3776 
3777  //Info<< "Seeding from processor " << procI << " face " << seedFaceI
3778  // << endl;
3779 
3780  if (procI == Pstream::myProcNo())
3781  {
3782  // Determine orientation of seedFace
3783 
3784  patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];
3785 
3786  // Start off with correct orientation
3787  faceInfo = orientedSurface::NOFLIP;
3788 
3789  if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
3790  {
3791  faceInfo.flip();
3792  }
3793 
3794 
3795  const labelList& fEdges = patch.faceEdges()[seedFaceI];
3796  forAll(fEdges, fEdgeI)
3797  {
3798  label edgeI = fEdges[fEdgeI];
3799 
3800  patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
3801 
3802  if
3803  (
3804  edgeInfo.updateEdge<int>
3805  (
3806  mesh_,
3807  patch,
3808  edgeI,
3809  seedFaceI,
3810  faceInfo,
3811  tol,
3812  dummyTrackData
3813  )
3814  )
3815  {
3816  changedEdges.append(edgeI);
3817  changedInfo.append(edgeInfo);
3818  }
3819  }
3820  }
3821 
3822 
3823  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
3824  {
3825  break;
3826  }
3827 
3828 
3829 
3830  // Walk
3831  PatchEdgeFaceWave
3832  <
3833  indirectPrimitivePatch,
3834  patchFaceOrientation
3835  > calc
3836  (
3837  mesh_,
3838  patch,
3839  changedEdges,
3840  changedInfo,
3841  allEdgeInfo,
3842  allFaceInfo,
3843  returnReduce(patch.nEdges(), sumOp<label>())
3844  );
3845  }
3846 
3847 
3848  // Push master zone info over to slave (since slave faces never visited)
3849  {
3850  labelList neiStatus
3851  (
3852  mesh_.nBoundaryFaces(),
3853  orientedSurface::UNVISITED
3854  );
3855 
3856  forAll(patch.addressing(), i)
3857  {
3858  const label meshFaceI = patch.addressing()[i];
3859  if (!mesh_.isInternalFace(meshFaceI))
3860  {
3861  neiStatus[meshFaceI-mesh_.nInternalFaces()] =
3862  allFaceInfo[i].flipStatus();
3863  }
3864  }
3865  syncTools::swapBoundaryFaceList(mesh_, neiStatus);
3866 
3867  forAll(patch.addressing(), i)
3868  {
3869  const label meshFaceI = patch.addressing()[i];
3870  const label patchI = bm.whichPatch(meshFaceI);
3871 
3872  if
3873  (
3874  patchI != -1
3875  && bm[patchI].coupled()
3876  && !isMasterFace.test(meshFaceI)
3877  )
3878  {
3879  // Slave side. Take flipped from neighbour
3880  label bFaceI = meshFaceI-mesh_.nInternalFaces();
3881 
3882  if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
3883  {
3884  allFaceInfo[i] = orientedSurface::FLIP;
3885  }
3886  else if (neiStatus[bFaceI] == orientedSurface::FLIP)
3887  {
3888  allFaceInfo[i] = orientedSurface::NOFLIP;
3889  }
3890  else
3891  {
3892  FatalErrorInFunction
3893  << "Incorrect status for face " << meshFaceI
3894  << abort(FatalError);
3895  }
3896  }
3897  }
3898  }
3899 
3900 
3901  // Convert to meshFlipMap and adapt faceZones
3902 
3903  meshFlipMap.setSize(mesh_.nFaces());
3904  meshFlipMap = false;
3905 
3906  forAll(allFaceInfo, faceI)
3907  {
3908  label meshFaceI = patch.addressing()[faceI];
3909 
3910  if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
3911  {
3912  meshFlipMap.unset(meshFaceI);
3913  }
3914  else if (allFaceInfo[faceI] == orientedSurface::FLIP)
3915  {
3916  meshFlipMap.set(meshFaceI);
3917  }
3918  else
3919  {
3920  FatalErrorInFunction
3921  << "Problem : unvisited face " << faceI
3922  << " centre:" << mesh_.faceCentres()[meshFaceI]
3923  << abort(FatalError);
3924  }
3925  }
3926 }
3927 
3928 
3929 void Foam::meshRefinement::zonify
3930 (
3931  // Get per face whether is it master (of a coupled set of faces)
3932  const bitSet& isMasterFace,
3933  const labelList& cellToZone,
3934  const labelList& neiCellZone,
3935  const labelList& faceToZone,
3936  const bitSet& meshFlipMap,
3937  polyTopoChange& meshMod
3938 ) const
3939 {
3940  const labelList& faceOwner = mesh_.faceOwner();
3941  const labelList& faceNeighbour = mesh_.faceNeighbour();
3942 
3943  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3944  {
3945  label faceZoneI = faceToZone[faceI];
3946 
3947  if (faceZoneI != -1)
3948  {
3949  // Orient face zone to have slave cells in min cell zone.
3950  // Note: logic to use flipMap should be consistent with logic
3951  // to pick up the freeStandingBaffleFaces!
3952 
3953  label ownZone = cellToZone[faceOwner[faceI]];
3954  label neiZone = cellToZone[faceNeighbour[faceI]];
3955 
3956  bool flip;
3957 
3958  if (ownZone == neiZone)
3959  {
3960  // free-standing face. Use geometrically derived orientation
3961  flip = meshFlipMap.test(faceI);
3962  }
3963  else
3964  {
3965  flip =
3966  (
3967  ownZone == -1
3968  || (neiZone != -1 && ownZone > neiZone)
3969  );
3970  }
3971 
3972  meshMod.setAction
3973  (
3974  polyModifyFace
3975  (
3976  mesh_.faces()[faceI], // modified face
3977  faceI, // label of face
3978  faceOwner[faceI], // owner
3979  faceNeighbour[faceI], // neighbour
3980  false, // face flip
3981  -1, // patch for face
3982  false, // remove from zone
3983  faceZoneI, // zone for face
3984  flip // face flip in zone
3985  )
3986  );
3987  }
3988  }
3989 
3990 
3991  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
3992 
3993  // Set owner as no-flip
3994  forAll(patches, patchI)
3995  {
3996  const polyPatch& pp = patches[patchI];
3997 
3998  label faceI = pp.start();
3999 
4000  forAll(pp, i)
4001  {
4002  label faceZoneI = faceToZone[faceI];
4003 
4004  if (faceZoneI != -1)
4005  {
4006  label ownZone = cellToZone[faceOwner[faceI]];
4007  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
4008 
4009  bool flip;
4010 
4011  if (ownZone == neiZone)
4012  {
4013  // free-standing face. Use geometrically derived orientation
4014  flip = meshFlipMap.test(faceI);
4015  }
4016  else
4017  {
4018  flip =
4019  (
4020  ownZone == -1
4021  || (neiZone != -1 && ownZone > neiZone)
4022  );
4023  }
4024 
4025  meshMod.setAction
4026  (
4027  polyModifyFace
4028  (
4029  mesh_.faces()[faceI], // modified face
4030  faceI, // label of face
4031  faceOwner[faceI], // owner
4032  -1, // neighbour
4033  false, // face flip
4034  patchI, // patch for face
4035  false, // remove from zone
4036  faceZoneI, // zone for face
4037  flip // face flip in zone
4038  )
4039  );
4040  }
4041  faceI++;
4042  }
4043  }
4044 
4045 
4046  // Put the cells into the correct zone
4047  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4048 
4049  forAll(cellToZone, cellI)
4050  {
4051  label zoneI = cellToZone[cellI];
4052 
4053  if (zoneI >= 0)
4054  {
4055  meshMod.setAction
4056  (
4057  polyModifyCell
4058  (
4059  cellI,
4060  false, // removeFromZone
4061  zoneI
4062  )
4063  );
4064  }
4065  }
4066 }
4067 
4068 
4069 void Foam::meshRefinement::allocateInterRegionFaceZone
4070 (
4071  const label ownZone,
4072  const label neiZone,
4073  wordPairHashTable& zonesToFaceZone,
4074  LabelPairMap<word>& zoneIDsToFaceZone
4075 ) const
4076 {
4077  const cellZoneMesh& cellZones = mesh_.cellZones();
4078 
4079  if (ownZone != neiZone)
4080  {
4081  // Make sure lowest number cellZone is master. Non-cellZone
4082  // areas are slave
4083  const bool swap =
4084  (
4085  ownZone == -1
4086  || (neiZone != -1 && ownZone > neiZone)
4087  );
4088 
4089  // Quick check whether we already have pair of zones
4090  labelPair key(ownZone, neiZone);
4091  if (swap)
4092  {
4093  key.flip();
4094  }
4095 
4096  if (!zoneIDsToFaceZone.found(key))
4097  {
4098  // Not found. Allocate.
4099  const word ownZoneName =
4100  (
4101  ownZone != -1
4102  ? cellZones[ownZone].name()
4103  : "none"
4104  );
4105  const word neiZoneName =
4106  (
4107  neiZone != -1
4108  ? cellZones[neiZone].name()
4109  : "none"
4110  );
4111 
4112  // Get lowest zone first
4113  Pair<word> wordKey(ownZoneName, neiZoneName);
4114  if (swap)
4115  {
4116  wordKey.flip();
4117  }
4118 
4119  word fzName = wordKey.first() + "_to_" + wordKey.second();
4120 
4121  zoneIDsToFaceZone.insert(key, fzName);
4122  zonesToFaceZone.insert(wordKey, fzName);
4123  }
4124  }
4125 }
4126 
4127 
4128 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
4129 
4130 void Foam::meshRefinement::baffleAndSplitMesh
4131 (
4132  const bool doHandleSnapProblems,
4133  const snapParameters& snapParams,
4134  const bool useTopologicalSnapDetection,
4135  const bool removeEdgeConnectedCells,
4136  const scalarField& perpendicularAngle,
4137  const label nErodeCellZones,
4138  const dictionary& motionDict,
4139  Time& runTime,
4140  const labelList& globalToMasterPatch,
4141  const labelList& globalToSlavePatch,
4142 
4143  const pointField& locationsInMesh,
4144  const wordList& zonesInMesh,
4145  const pointField& locationsOutsideMesh,
4146  const writer<scalar>& leakPathFormatter
4147 )
4148 {
4149  // Introduce baffles
4150  // ~~~~~~~~~~~~~~~~~
4151 
4152  // Split the mesh along internal faces wherever there is a pierce between
4153  // two cell centres.
4154 
4155  Info<< "Introducing baffles for "
4156  << returnReduce(countHits(), sumOp<label>())
4157  << " faces that are intersected by the surface." << nl << endl;
4158 
4159  // Swap neighbouring cell centres and cell level
4160  labelList neiLevel(mesh_.nBoundaryFaces());
4161  pointField neiCc(mesh_.nBoundaryFaces());
4162  calcNeighbourData(neiLevel, neiCc);
4163 
4164  labelList ownPatch, neiPatch;
4165  getBafflePatches
4166  (
4167  nErodeCellZones,
4168  globalToMasterPatch,
4169 
4170  locationsInMesh,
4171  zonesInMesh,
4172 
4173  neiLevel,
4174  neiCc,
4175 
4176  ownPatch,
4177  neiPatch
4178  );
4179 
4180  createBaffles(ownPatch, neiPatch);
4181 
4182  if (debug)
4183  {
4184  // Debug:test all is still synced across proc patches
4185  checkData();
4186  }
4187 
4188  Info<< "Created baffles in = "
4189  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4190 
4191  printMeshInfo(debug, "After introducing baffles");
4192 
4193  if (debug&MESH)
4194  {
4195  const_cast<Time&>(mesh_.time())++;
4196  Pout<< "Writing baffled mesh to time " << timeName()
4197  << endl;
4198  write
4199  (
4200  debugType(debug),
4201  writeType(writeLevel() | WRITEMESH),
4202  runTime.path()/"baffles"
4203  );
4204  Pout<< "Dumped debug data in = "
4205  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4206  }
4207 
4208 
4209  // Introduce baffles to delete problem cells
4210  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4211 
4212  // Create some additional baffles where we want surface cells removed.
4213 
4214  if (doHandleSnapProblems)
4215  {
4216  handleSnapProblems
4217  (
4218  snapParams,
4219  useTopologicalSnapDetection,
4220  removeEdgeConnectedCells,
4221  perpendicularAngle,
4222  motionDict,
4223  runTime,
4224  globalToMasterPatch,
4225  globalToSlavePatch
4226  );
4227 
4228  // Removing additional cells might have created disconnected bits
4229  // so re-do the surface intersections
4230  {
4231  // Swap neighbouring cell centres and cell level
4232  neiLevel.setSize(mesh_.nBoundaryFaces());
4233  neiCc.setSize(mesh_.nBoundaryFaces());
4234  calcNeighbourData(neiLevel, neiCc);
4235 
4236  labelList ownPatch, neiPatch;
4237  getBafflePatches
4238  (
4239  nErodeCellZones,
4240  globalToMasterPatch,
4241 
4242  locationsInMesh,
4243  zonesInMesh,
4244 
4245  neiLevel,
4246  neiCc,
4247 
4248  ownPatch,
4249  neiPatch
4250  );
4251 
4252  createBaffles(ownPatch, neiPatch);
4253  }
4254 
4255  if (debug)
4256  {
4257  // Debug:test all is still synced across proc patches
4258  checkData();
4259  }
4260  }
4261 
4262 
4263  // Select part of mesh
4264  // ~~~~~~~~~~~~~~~~~~~
4265 
4266  Info<< nl
4267  << "Remove unreachable sections of mesh" << nl
4268  << "-----------------------------------" << nl
4269  << endl;
4270 
4271  if (debug)
4272  {
4273  ++runTime;
4274  }
4275 
4276  splitMeshRegions
4277  (
4278  globalToMasterPatch,
4279  globalToSlavePatch,
4280  locationsInMesh,
4281  locationsOutsideMesh,
4282  true, // Exit if any connection between inside and outside
4283  leakPathFormatter
4284  );
4285 
4286  if (debug)
4287  {
4288  // Debug:test all is still synced across proc patches
4289  checkData();
4290  }
4291  Info<< "Split mesh in = "
4292  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4293 
4294  printMeshInfo(debug, "After subsetting");
4295 
4296  if (debug&MESH)
4297  {
4298  ++runTime;
4299  Pout<< "Writing subsetted mesh to time " << timeName()
4300  << endl;
4301  write
4302  (
4303  debugType(debug),
4304  writeType(writeLevel() | WRITEMESH),
4305  runTime.path()/timeName()
4306  );
4307  Pout<< "Dumped debug data in = "
4308  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4309  }
4310 }
4311 
4312 
4313 void Foam::meshRefinement::mergeFreeStandingBaffles
4314 (
4315  const snapParameters& snapParams,
4316  const bool useTopologicalSnapDetection,
4317  const bool removeEdgeConnectedCells,
4318  const scalarField& perpendicularAngle,
4319  const scalar planarAngle,
4320  const dictionary& motionDict,
4321  Time& runTime,
4322  const labelList& globalToMasterPatch,
4323  const labelList& globalToSlavePatch,
4324  const pointField& locationsInMesh,
4325  const pointField& locationsOutsideMesh,
4326  const writer<scalar>& leakPathFormatter
4327 )
4328 {
4329  // Merge baffles
4330  // ~~~~~~~~~~~~~
4331 
4332  Info<< nl
4333  << "Merge free-standing baffles" << nl
4334  << "---------------------------" << nl
4335  << endl;
4336 
4337 
4338  // List of pairs of freestanding baffle faces.
4339  List<labelPair> couples
4340  (
4341  freeStandingBaffles // filter out freestanding baffles
4342  (
4343  localPointRegion::findDuplicateFacePairs(mesh_),
4344  planarAngle
4345  )
4346  );
4347 
4348  label nCouples = couples.size();
4349  reduce(nCouples, sumOp<label>());
4350 
4351  Info<< "Detected free-standing baffles : " << nCouples << endl;
4352 
4353  if (nCouples > 0)
4354  {
4355  // Actually merge baffles. Note: not exactly parallellized. Should
4356  // convert baffle faces into processor faces if they resulted
4357  // from them.
4358  mergeBaffles(couples, Map<label>(0));
4359 
4360  // Detect any problem cells resulting from merging of baffles
4361  // and delete them
4362  handleSnapProblems
4363  (
4364  snapParams,
4365  useTopologicalSnapDetection,
4366  removeEdgeConnectedCells,
4367  perpendicularAngle,
4368  motionDict,
4369  runTime,
4370  globalToMasterPatch,
4371  globalToSlavePatch
4372  );
4373 
4374  // Very occasionally removing a problem cell might create a disconnected
4375  // region so re-check
4376 
4377  Info<< nl
4378  << "Remove unreachable sections of mesh" << nl
4379  << "-----------------------------------" << nl
4380  << endl;
4381 
4382  if (debug)
4383  {
4384  ++runTime;
4385  }
4386 
4387  splitMeshRegions
4388  (
4389  globalToMasterPatch,
4390  globalToSlavePatch,
4391  locationsInMesh,
4392  locationsOutsideMesh,
4393  true, // Exit if any connection between inside and outside
4394  leakPathFormatter
4395  );
4396 
4397 
4398  if (debug)
4399  {
4400  // Debug:test all is still synced across proc patches
4401  checkData();
4402  }
4403  }
4404  Info<< "Merged free-standing baffles in = "
4405  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4406 }
4407 
4408 
4409 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
4410 (
4411  const label nBufferLayers,
4412  const label nErodeCellZones,
4413  const labelList& globalToMasterPatch,
4414  const labelList& globalToSlavePatch,
4415 
4416  const pointField& locationsInMesh,
4417  const wordList& zonesInMesh,
4418  const pointField& locationsOutsideMesh,
4419  const writer<scalar>& leakPathFormatter
4420 )
4421 {
4422  // Determine patches to put intersections into
4423  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4424 
4425  // Swap neighbouring cell centres and cell level
4426  labelList neiLevel(mesh_.nBoundaryFaces());
4427  pointField neiCc(mesh_.nBoundaryFaces());
4428  calcNeighbourData(neiLevel, neiCc);
4429 
4430  // Find intersections with all unnamed(!) surfaces
4431  labelList ownPatch, neiPatch;
4432  getBafflePatches
4433  (
4434  nErodeCellZones,
4435  globalToMasterPatch,
4436 
4437  locationsInMesh,
4438  zonesInMesh,
4439 
4440  neiLevel,
4441  neiCc,
4442 
4443  ownPatch,
4444  neiPatch
4445  );
4446 
4447  // Analyse regions. Reuse regionsplit
4448  boolList blockedFace(mesh_.nFaces(), false);
4449 
4450  forAll(ownPatch, faceI)
4451  {
4452  if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
4453  {
4454  blockedFace[faceI] = true;
4455  }
4456  }
4457  syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>());
4458 
4459 
4460  regionSplit cellRegion(mesh_, blockedFace);
4461 
4462  // Set unreachable cells to -1
4463  findRegions
4464  (
4465  mesh_,
4466  mergeDistance_ * vector::one, // perturbVec
4467  locationsInMesh,
4468  locationsOutsideMesh,
4469  false, // do not exit if outside location found
4470  leakPathFormatter,
4471  cellRegion.nRegions(),
4472  cellRegion,
4473  blockedFace
4474  );
4475 
4476  return splitMesh
4477  (
4478  nBufferLayers,
4479  globalToMasterPatch,
4480  globalToSlavePatch,
4481  cellRegion,
4482  ownPatch,
4483  neiPatch
4484  );
4485 }
4486 
4487 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
4488 (
4489  const label nBufferLayers,
4490  const labelList& globalToMasterPatch,
4491  const labelList& globalToSlavePatch,
4492 
4493  labelList& cellRegion,
4494  labelList& ownPatch,
4495  labelList& neiPatch
4496 )
4497 {
4498  // Walk out nBufferlayers from region boundary
4499  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4500  // (modifies cellRegion, ownPatch)
4501  // Takes over face patch onto points and then back to faces and cells
4502  // (so cell-face-point walk)
4503 
4504  const labelList& faceOwner = mesh_.faceOwner();
4505  const labelList& faceNeighbour = mesh_.faceNeighbour();
4506 
4507  // Patch for exposed faces for lack of anything sensible.
4508  label defaultPatch = 0;
4509  if (globalToMasterPatch.size())
4510  {
4511  defaultPatch = globalToMasterPatch[0];
4512  }
4513 
4514  for (label i = 0; i < nBufferLayers; i++)
4515  {
4516  // 1. From cells (via faces) to points
4517 
4518  labelList pointBaffle(mesh_.nPoints(), -1);
4519 
4520  forAll(faceNeighbour, faceI)
4521  {
4522  const face& f = mesh_.faces()[faceI];
4523 
4524  label ownRegion = cellRegion[faceOwner[faceI]];
4525  label neiRegion = cellRegion[faceNeighbour[faceI]];
4526 
4527  if (ownRegion == -1 && neiRegion != -1)
4528  {
4529  // Note max(..) since possibly regionSplit might have split
4530  // off extra unreachable parts of mesh. Note: or can this only
4531  // happen for boundary faces?
4532  forAll(f, fp)
4533  {
4534  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
4535  }
4536  }
4537  else if (ownRegion != -1 && neiRegion == -1)
4538  {
4539  label newPatchI = neiPatch[faceI];
4540  if (newPatchI == -1)
4541  {
4542  newPatchI = max(defaultPatch, ownPatch[faceI]);
4543  }
4544  forAll(f, fp)
4545  {
4546  pointBaffle[f[fp]] = newPatchI;
4547  }
4548  }
4549  }
4550  for
4551  (
4552  label faceI = mesh_.nInternalFaces();
4553  faceI < mesh_.nFaces();
4554  faceI++
4555  )
4556  {
4557  const face& f = mesh_.faces()[faceI];
4558 
4559  label ownRegion = cellRegion[faceOwner[faceI]];
4560 
4561  if (ownRegion == -1)
4562  {
4563  forAll(f, fp)
4564  {
4565  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
4566  }
4567  }
4568  }
4569 
4570  // Sync
4571  syncTools::syncPointList
4572  (
4573  mesh_,
4574  pointBaffle,
4575  maxEqOp<label>(),
4576  label(-1) // null value
4577  );
4578 
4579 
4580  // 2. From points back to faces
4581 
4582  const labelListList& pointFaces = mesh_.pointFaces();
4583 
4584  forAll(pointFaces, pointI)
4585  {
4586  if (pointBaffle[pointI] != -1)
4587  {
4588  const labelList& pFaces = pointFaces[pointI];
4589 
4590  forAll(pFaces, pFaceI)
4591  {
4592  label faceI = pFaces[pFaceI];
4593 
4594  if (ownPatch[faceI] == -1)
4595  {
4596  ownPatch[faceI] = pointBaffle[pointI];
4597  }
4598  }
4599  }
4600  }
4601  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
4602 
4603 
4604  // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
4605 
4606  labelList newOwnPatch(ownPatch);
4607 
4608  forAll(ownPatch, faceI)
4609  {
4610  if (ownPatch[faceI] != -1)
4611  {
4612  label own = faceOwner[faceI];
4613 
4614  if (cellRegion[own] == -1)
4615  {
4616  cellRegion[own] = labelMax;
4617 
4618  const cell& ownFaces = mesh_.cells()[own];
4619  forAll(ownFaces, j)
4620  {
4621  if (ownPatch[ownFaces[j]] == -1)
4622  {
4623  newOwnPatch[ownFaces[j]] = ownPatch[faceI];
4624  }
4625  }
4626  }
4627  if (mesh_.isInternalFace(faceI))
4628  {
4629  label nei = faceNeighbour[faceI];
4630 
4631  if (cellRegion[nei] == -1)
4632  {
4633  cellRegion[nei] = labelMax;
4634 
4635  const cell& neiFaces = mesh_.cells()[nei];
4636  forAll(neiFaces, j)
4637  {
4638  if (ownPatch[neiFaces[j]] == -1)
4639  {
4640  newOwnPatch[neiFaces[j]] = ownPatch[faceI];
4641  }
4642  }
4643  }
4644  }
4645  }
4646  }
4647 
4648  ownPatch.transfer(newOwnPatch);
4649 
4650  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
4651  }
4652 
4653 
4654 
4655  // Subset
4656  // ~~~~~~
4657 
4658  // Get cells to remove
4659  DynamicList<label> cellsToRemove(mesh_.nCells());
4660  forAll(cellRegion, cellI)
4661  {
4662  if (cellRegion[cellI] == -1)
4663  {
4664  cellsToRemove.append(cellI);
4665  }
4666  }
4667  cellsToRemove.shrink();
4668 
4669  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
4670  reduce(nCellsToKeep, sumOp<label>());
4671 
4672  Info<< "Keeping all cells containing inside points" << endl
4673  << "Selected for keeping : " << nCellsToKeep << " cells." << endl;
4674 
4675 
4676  // Remove cells
4677  removeCells cellRemover(mesh_);
4678 
4679  // Pick up patches for exposed faces
4680  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
4681  labelList exposedPatches(exposedFaces.size());
4682 
4683  forAll(exposedFaces, i)
4684  {
4685  label faceI = exposedFaces[i];
4686 
4687  if (ownPatch[faceI] != -1)
4688  {
4689  exposedPatches[i] = ownPatch[faceI];
4690  }
4691  else
4692  {
4693  WarningInFunction
4694  << "For exposed face " << faceI
4695  << " fc:" << mesh_.faceCentres()[faceI]
4696  << " found no patch." << endl
4697  << " Taking patch " << defaultPatch
4698  << " instead." << endl;
4699  exposedPatches[i] = defaultPatch;
4700  }
4701  }
4702 
4703  return doRemoveCells
4704  (
4705  cellsToRemove,
4706  exposedFaces,
4707  exposedPatches,
4708  cellRemover
4709  );
4710 }
4711 
4712 
4713 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLimitShells
4714 (
4715  const label nBufferLayers,
4716  const label nErodeCellZones,
4717  const labelList& globalToMasterPatch,
4718  const labelList& globalToSlavePatch,
4719  const pointField& locationsInMesh,
4720  const wordList& zonesInMesh
4721 )
4722 {
4723  // Determine patches to put intersections into
4724  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4725 
4726  // Swap neighbouring cell centres and cell level
4727  labelList neiLevel(mesh_.nBoundaryFaces());
4728  pointField neiCc(mesh_.nBoundaryFaces());
4729  calcNeighbourData(neiLevel, neiCc);
4730 
4731  // Find intersections with all unnamed(!) surfaces
4732  labelList ownPatch, neiPatch;
4733  getBafflePatches
4734  (
4735  nErodeCellZones,
4736  globalToMasterPatch,
4737 
4738  locationsInMesh,
4739  zonesInMesh,
4740 
4741  neiLevel,
4742  neiCc,
4743 
4744  ownPatch,
4745  neiPatch
4746  );
4747 
4748 
4749  labelList cellRegion(mesh_.nCells(), Zero);
4750  // Find any cells inside a limit shell with minLevel -1
4751  labelList levelShell;
4752  limitShells_.findLevel
4753  (
4754  mesh_.cellCentres(),
4755  labelList(mesh_.nCells(), -1), // pick up only shells with -1
4756  levelShell
4757  );
4758  forAll(levelShell, celli)
4759  {
4760  if (levelShell[celli] != -1)
4761  {
4762  // Mark cell region so it gets deletec
4763  cellRegion[celli] = -1;
4764  }
4765  }
4766 
4767  autoPtr<mapPolyMesh> mapPtr = splitMesh
4768  (
4769  nBufferLayers,
4770  globalToMasterPatch,
4771  globalToSlavePatch,
4772  cellRegion,
4773  ownPatch,
4774  neiPatch
4775  );
4776 
4777  if (debug&meshRefinement::MESH)
4778  {
4779  const_cast<Time&>(mesh_.time())++;
4780  Pout<< "Writing mesh after removing limitShells"
4781  << " to time " << timeName() << endl;
4782  write
4783  (
4784  debugType(debug),
4785  writeType
4786  (
4787  writeLevel()
4788  | WRITEMESH
4789  ),
4790  mesh_.time().path()/timeName()
4791  );
4792  }
4793  return mapPtr;
4794 }
4795 
4796 
4797 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
4798 (
4799  const localPointRegion& regionSide
4800 )
4801 {
4802  // Topochange container
4803  polyTopoChange meshMod(mesh_);
4804 
4805  label nNonManifPoints = returnReduce
4806  (
4807  regionSide.meshPointMap().size(),
4808  sumOp<label>()
4809  );
4810 
4811  Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
4812  << " non-manifold points (out of "
4813  << mesh_.globalData().nTotalPoints()
4814  << ')' << endl;
4815 
4816 
4817  autoPtr<mapPolyMesh> mapPtr;
4818 
4819  if (nNonManifPoints)
4820  {
4821  // Topo change engine
4822  duplicatePoints pointDuplicator(mesh_);
4823 
4824  // Insert changes into meshMod
4825  pointDuplicator.setRefinement(regionSide, meshMod);
4826 
4827  // Change the mesh (no inflation, parallel sync)
4828  mapPtr = meshMod.changeMesh(mesh_, false, true);
4829  mapPolyMesh& map = *mapPtr;
4830 
4831  // Update fields
4832  mesh_.updateMesh(map);
4833 
4834  // Move mesh if in inflation mode
4835  if (map.hasMotionPoints())
4836  {
4837  mesh_.movePoints(map.preMotionPoints());
4838  }
4839  else
4840  {
4841  // Delete mesh volumes.
4842  mesh_.clearOut();
4843  }
4844 
4845  // Reset the instance for if in overwrite mode
4846  mesh_.setInstance(timeName());
4847 
4848  // Update intersections. Is mapping only (no faces created, positions
4849  // stay same) so no need to recalculate intersections.
4850  updateMesh(map, labelList());
4851  }
4852 
4853  return mapPtr;
4854 }
4855 
4856 
4857 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
4858 {
4859  // Analyse which points need to be duplicated
4860  localPointRegion regionSide(mesh_);
4861 
4862  return dupNonManifoldPoints(regionSide);
4863 }
4864 
4865 
4866 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergePoints
4867 (
4868  const labelList& pointToDuplicate
4869 )
4870 {
4871  label nPointPairs = 0;
4872  forAll(pointToDuplicate, pointI)
4873  {
4874  label otherPointI = pointToDuplicate[pointI];
4875  if (otherPointI != -1)
4876  {
4877  nPointPairs++;
4878  }
4879  }
4880 
4881  autoPtr<mapPolyMesh> mapPtr;
4882 
4883  if (returnReduce(nPointPairs, sumOp<label>()))
4884  {
4885  Map<label> pointToMaster(2*nPointPairs);
4886  forAll(pointToDuplicate, pointI)
4887  {
4888  label otherPointI = pointToDuplicate[pointI];
4889  if (otherPointI != -1)
4890  {
4891  // Slave point
4892  pointToMaster.insert(pointI, otherPointI);
4893  }
4894  }
4895 
4896  // Topochange container
4897  polyTopoChange meshMod(mesh_);
4898 
4899  // Insert changes
4900  polyMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
4901 
4902  // Change the mesh (no inflation, parallel sync)
4903  mapPtr = meshMod.changeMesh(mesh_, false, true);
4904  mapPolyMesh& map = *mapPtr;
4905 
4906  // Update fields
4907  mesh_.updateMesh(map);
4908 
4909  // Move mesh if in inflation mode
4910  if (map.hasMotionPoints())
4911  {
4912  mesh_.movePoints(map.preMotionPoints());
4913  }
4914  else
4915  {
4916  // Delete mesh volumes.
4917  mesh_.clearOut();
4918  }
4919 
4920  // Reset the instance for if in overwrite mode
4921  mesh_.setInstance(timeName());
4922 
4923  // Update intersections. Is mapping only (no faces created, positions
4924  // stay same) so no need to recalculate intersections.
4925  updateMesh(map, labelList());
4926  }
4927 
4928  return mapPtr;
4929 }
4930 
4931 
4932 // Duplicate points on 'boundary' zones. Do not duplicate points on
4933 // 'internal' or 'baffle' zone. Whether points are on normal patches does
4934 // not matter
4935 Foam::autoPtr<Foam::mapPolyMesh>
4936 Foam::meshRefinement::dupNonManifoldBoundaryPoints()
4937 {
4938  const labelList boundaryFaceZones
4939  (
4940  getZones
4941  (
4942  List<surfaceZonesInfo::faceZoneType>
4943  (
4944  1,
4945  surfaceZonesInfo::BOUNDARY
4946  )
4947  )
4948  );
4949  labelList internalOrBaffleFaceZones;
4950  {
4951  List<surfaceZonesInfo::faceZoneType> fzTypes(2);
4952  fzTypes[0] = surfaceZonesInfo::INTERNAL;
4953  fzTypes[1] = surfaceZonesInfo::BAFFLE;
4954  internalOrBaffleFaceZones = getZones(fzTypes);
4955  }
4956 
4957 
4958 
4959  // 0 : point used by normal, unzoned boundary faces
4960  // 1 : point used by 'boundary' zone
4961  // 2 : point used by internal/baffle zone
4962  PackedList<2> pointStatus(mesh_.nPoints(), 0u);
4963 
4964  forAll(boundaryFaceZones, j)
4965  {
4966  const faceZone& fZone = mesh_.faceZones()[boundaryFaceZones[j]];
4967  forAll(fZone, i)
4968  {
4969  const face& f = mesh_.faces()[fZone[i]];
4970  forAll(f, fp)
4971  {
4972  pointStatus[f[fp]] = max(pointStatus[f[fp]], 1u);
4973  }
4974  }
4975  }
4976  forAll(internalOrBaffleFaceZones, j)
4977  {
4978  const faceZone& fZone = mesh_.faceZones()[internalOrBaffleFaceZones[j]];
4979  forAll(fZone, i)
4980  {
4981  const face& f = mesh_.faces()[fZone[i]];
4982  forAll(f, fp)
4983  {
4984  pointStatus[f[fp]] = max(pointStatus[f[fp]], 2u);
4985  }
4986  }
4987  }
4988 
4989  syncTools::syncPointList
4990  (
4991  mesh_,
4992  pointStatus,
4993  maxEqOp<unsigned int>(), // combine op
4994  0u // null value
4995  );
4996 
4997  // Pick up points on boundary zones that are not on internal/baffle zones
4998  label n = 0;
4999  forAll(pointStatus, pointI)
5000  {
5001  if (pointStatus[pointI] == 1u)
5002  {
5003  n++;
5004  }
5005  }
5006 
5007  label globalNPoints = returnReduce(n, sumOp<label>());
5008  Info<< "Duplicating " << globalNPoints << " points on"
5009  << " faceZones of type "
5010  << surfaceZonesInfo::faceZoneTypeNames[surfaceZonesInfo::BOUNDARY]
5011  << endl;
5012 
5013  autoPtr<mapPolyMesh> map;
5014 
5015  if (globalNPoints)
5016  {
5017  labelList candidatePoints(n);
5018  n = 0;
5019  forAll(pointStatus, pointI)
5020  {
5021  if (pointStatus[pointI] == 1u)
5022  {
5023  candidatePoints[n++] = pointI;
5024  }
5025  }
5026  localPointRegion regionSide(mesh_, candidatePoints);
5027  map = dupNonManifoldPoints(regionSide);
5028  }
5029  return map;
5030 }
5031 
5032 
5033 // Zoning
5034 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
5035 (
5036  const bool allowFreeStandingZoneFaces,
5037  const label nErodeCellZones,
5038  const pointField& locationsInMesh,
5039  const wordList& zonesInMesh,
5040  wordPairHashTable& zonesToFaceZone
5041 )
5042 {
5043  if (locationsInMesh.size() != zonesInMesh.size())
5044  {
5045  FatalErrorInFunction << "problem" << abort(FatalError);
5046  }
5047 
5048  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
5049  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
5050 
5051 
5052  // Swap neighbouring cell centres and cell level
5053  labelList neiLevel(mesh_.nBoundaryFaces());
5054  pointField neiCc(mesh_.nBoundaryFaces());
5055  calcNeighbourData(neiLevel, neiCc);
5056 
5057 
5058  // Add any faceZones, cellZones originating from surface to the mesh
5059  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5060  labelList surfaceToCellZone;
5061  labelListList surfaceToFaceZones;
5062 
5063  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
5064  if (namedSurfaces.size())
5065  {
5066  Info<< "Setting cellZones according to named surfaces:" << endl;
5067  forAll(namedSurfaces, i)
5068  {
5069  label surfI = namedSurfaces[i];
5070 
5071  Info<< "Surface : " << surfaces_.names()[surfI] << nl
5072  << " faceZones : " << surfZones[surfI].faceZoneNames() << nl
5073  << " cellZone : " << surfZones[surfI].cellZoneName()
5074  << endl;
5075  }
5076  Info<< endl;
5077 
5078  // Add zones to mesh
5079  surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
5080  (
5081  surfZones,
5082  namedSurfaces,
5083  mesh_
5084  );
5085  surfaceToFaceZones = surfaceZonesInfo::addFaceZonesToMesh
5086  (
5087  surfZones,
5088  namedSurfaces,
5089  mesh_
5090  );
5091  }
5092 
5093 
5094 
5095  // Put the cells into the correct zone
5096  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5097 
5098  // Zone per cell:
5099  // -2 : unset : not allowed!
5100  // -1 : not in any zone (zone 'none')
5101  // >=0: zoneID
5102  // namedSurfaceRegion: zero sized or:
5103  // -1 : face not intersecting a named surface
5104  // >=0 : index of named surface
5105  labelList cellToZone;
5106  labelList namedSurfaceRegion;
5107  bitSet posOrientation;
5108  {
5109  labelList unnamedRegion1;
5110  labelList unnamedRegion2;
5111 
5112  zonify
5113  (
5114  allowFreeStandingZoneFaces,
5115  nErodeCellZones,// Use erosion (>0) or regionSplit to clean up
5116  -1, // Set all cells with cellToZone -2 to -1
5117  locationsInMesh,
5118  zonesInMesh,
5119 
5120  cellToZone,
5121  unnamedRegion1,
5122  unnamedRegion2,
5123  namedSurfaceRegion,
5124  posOrientation
5125  );
5126  }
5127 
5128 
5129  // Convert namedSurfaceRegion (index of named surfaces) to
5130  // actual faceZone index
5131 
5132  //- Per face index of faceZone or -1
5133  labelList faceToZone(mesh_.nFaces(), -1);
5134 
5135  forAll(namedSurfaceRegion, faceI)
5136  {
5137  //label surfI = namedSurfaceIndex[faceI];
5138  label globalI = namedSurfaceRegion[faceI];
5139  if (globalI != -1)
5140  {
5141  const labelPair spr = surfaces_.whichSurface(globalI);
5142  faceToZone[faceI] = surfaceToFaceZones[spr.first()][spr.second()];
5143  }
5144  }
5145 
5146 
5147 
5148  // Allocate and assign faceZones from cellZones
5149  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5150 
5151  {
5152  // 1. Detect inter-region face and allocate names
5153 
5154  HashTable<word, labelPair, labelPair::Hash<>> zoneIDsToFaceZone;
5155 
5156  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5157  {
5158  if (faceToZone[faceI] == -1) // Not named surface
5159  {
5160  // Face not yet in a faceZone. (it might already have been
5161  // done so by a 'named' surface). Check if inbetween different
5162  // cellZones
5163  allocateInterRegionFaceZone
5164  (
5165  cellToZone[mesh_.faceOwner()[faceI]],
5166  cellToZone[mesh_.faceNeighbour()[faceI]],
5167  zonesToFaceZone,
5168  zoneIDsToFaceZone
5169  );
5170  }
5171  }
5172 
5173  labelList neiCellZone;
5174  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
5175 
5176  forAll(neiCellZone, bFaceI)
5177  {
5178  label faceI = bFaceI + mesh_.nInternalFaces();
5179  if (faceToZone[faceI] == -1)
5180  {
5181  allocateInterRegionFaceZone
5182  (
5183  cellToZone[mesh_.faceOwner()[faceI]],
5184  neiCellZone[bFaceI],
5185  zonesToFaceZone,
5186  zoneIDsToFaceZone
5187  );
5188  }
5189  }
5190 
5191 
5192  // 2.Combine faceZoneNames allocated on different processors
5193 
5194  Pstream::mapCombineGather(zonesToFaceZone, eqOp<word>());
5195  Pstream::mapCombineScatter(zonesToFaceZone);
5196 
5197 
5198  // 3. Allocate faceZones from (now synchronised) faceZoneNames
5199  // Note: the faceZoneNames contain the same data but in different
5200  // order. We could sort the contents but instead just loop
5201  // in sortedToc order.
5202 
5203  Info<< "Setting faceZones according to neighbouring cellZones:"
5204  << endl;
5205 
5206  // From cellZone indices to faceZone index
5207  HashTable<label, labelPair, labelPair::Hash<>> fZoneLookup
5208  (
5209  zonesToFaceZone.size()
5210  );
5211 
5212  const cellZoneMesh& cellZones = mesh_.cellZones();
5213 
5214  {
5215  List<Pair<word>> czs(zonesToFaceZone.sortedToc());
5216 
5217  forAll(czs, i)
5218  {
5219  const Pair<word>& cz = czs[i];
5220  const word& fzName = zonesToFaceZone[cz];
5221 
5222  Info<< indent<< "cellZones : "
5223  << cz[0] << ' ' << cz[1] << nl
5224  << " faceZone : " << fzName << endl;
5225 
5226  label faceZoneI = surfaceZonesInfo::addFaceZone
5227  (
5228  fzName, // name
5229  labelList(0), // addressing
5230  boolList(0), // flipMap
5231  mesh_
5232  );
5233 
5234  label cz0 = cellZones.findZoneID(cz[0]);
5235  label cz1 = cellZones.findZoneID(cz[1]);
5236 
5237  fZoneLookup.insert(labelPair(cz0, cz1), faceZoneI);
5238  }
5239  }
5240 
5241 
5242  // 4. Set faceToZone with new faceZones
5243 
5244 
5245  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5246  {
5247  if (faceToZone[faceI] == -1)
5248  {
5249  // Face not yet in a faceZone. (it might already have been
5250  // done so by a 'named' surface). Check if inbetween different
5251  // cellZones
5252 
5253  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
5254  label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
5255  if (ownZone != neiZone)
5256  {
5257  const bool swap =
5258  (
5259  ownZone == -1
5260  || (neiZone != -1 && ownZone > neiZone)
5261  );
5262  labelPair key(ownZone, neiZone);
5263  if (swap)
5264  {
5265  key.flip();
5266  }
5267  faceToZone[faceI] = fZoneLookup[key];
5268  }
5269  }
5270  }
5271  forAll(neiCellZone, bFaceI)
5272  {
5273  label faceI = bFaceI + mesh_.nInternalFaces();
5274  if (faceToZone[faceI] == -1)
5275  {
5276  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
5277  label neiZone = neiCellZone[bFaceI];
5278  if (ownZone != neiZone)
5279  {
5280  const bool swap =
5281  (
5282  ownZone == -1
5283  || (neiZone != -1 && ownZone > neiZone)
5284  );
5285  labelPair key(ownZone, neiZone);
5286  if (swap)
5287  {
5288  key.flip();
5289  }
5290  faceToZone[faceI] = fZoneLookup[key];
5291  }
5292  }
5293  }
5294  Info<< endl;
5295  }
5296 
5297 
5298 
5299 
5300  // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
5301  labelList neiCellZone;
5302  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
5303  forAll(patches, patchI)
5304  {
5305  const polyPatch& pp = patches[patchI];
5306 
5307  if (!pp.coupled())
5308  {
5309  label bFaceI = pp.start()-mesh_.nInternalFaces();
5310  forAll(pp, i)
5311  {
5312  neiCellZone[bFaceI++] = -1;
5313  }
5314  }
5315  }
5316 
5317 
5318 
5319  // Get per face whether is it master (of a coupled set of faces)
5320  const bitSet isMasterFace(syncTools::getMasterFaces(mesh_));
5321 
5322 
5323  // faceZones
5324  // ~~~~~~~~~
5325  // Faces on faceZones come in two variants:
5326  // - faces on the outside of a cellZone. They will be oriented to
5327  // point out of the maximum cellZone.
5328  // - free-standing faces. These will be oriented according to the
5329  // local surface normal. We do this in a two step algorithm:
5330  // - do a consistent orientation
5331  // - check number of faces with consistent orientation
5332  // - if <0 flip the whole patch
5333  bitSet meshFlipMap(mesh_.nFaces(), false);
5334  {
5335  // Collect all data on zone faces without cellZones on either side.
5336  const indirectPrimitivePatch patch
5337  (
5338  IndirectList<face>
5339  (
5340  mesh_.faces(),
5341  freeStandingBaffleFaces
5342  (
5343  faceToZone,
5344  cellToZone,
5345  neiCellZone
5346  )
5347  ),
5348  mesh_.points()
5349  );
5350 
5351  label nFreeStanding = returnReduce(patch.size(), sumOp<label>());
5352  if (nFreeStanding > 0)
5353  {
5354  Info<< "Detected " << nFreeStanding << " free-standing zone faces"
5355  << endl;
5356 
5357  if (debug)
5358  {
5359  OBJstream str(mesh_.time().path()/"freeStanding.obj");
5360  str.write(patch.localFaces(), patch.localPoints(), false);
5361  }
5362 
5363 
5364  // Detect non-manifold edges
5365  labelList nMasterFacesPerEdge;
5366  calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
5367 
5368  // Mark zones. Even a single original surface might create multiple
5369  // disconnected/non-manifold-connected zones
5370  labelList faceToConnectedZone;
5371  const label nZones = markPatchZones
5372  (
5373  patch,
5374  nMasterFacesPerEdge,
5375  faceToConnectedZone
5376  );
5377 
5378  Map<label> nPosOrientation(2*nZones);
5379  for (label zoneI = 0; zoneI < nZones; zoneI++)
5380  {
5381  nPosOrientation.insert(zoneI, 0);
5382  }
5383 
5384  // Make orientations consistent in a topological way. This just
5385  // checks the first face per zone for whether nPosOrientation
5386  // is negative (which it never is at this point)
5387  consistentOrientation
5388  (
5389  isMasterFace,
5390  patch,
5391  nMasterFacesPerEdge,
5392  faceToConnectedZone,
5393  nPosOrientation,
5394 
5395  meshFlipMap
5396  );
5397 
5398  // Count per region the number of orientations (taking the new
5399  // flipMap into account)
5400  forAll(patch.addressing(), faceI)
5401  {
5402  label meshFaceI = patch.addressing()[faceI];
5403 
5404  if (isMasterFace.test(meshFaceI))
5405  {
5406  label n = 1;
5407  if
5408  (
5409  posOrientation.test(meshFaceI)
5410  == meshFlipMap.test(meshFaceI)
5411  )
5412  {
5413  n = -1;
5414  }
5415 
5416  nPosOrientation.find(faceToConnectedZone[faceI])() += n;
5417  }
5418  }
5419  Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
5420  Pstream::mapCombineScatter(nPosOrientation);
5421 
5422 
5423  Info<< "Split " << nFreeStanding << " free-standing zone faces"
5424  << " into " << nZones << " disconnected regions with size"
5425  << " (negative denotes wrong orientation) :"
5426  << endl;
5427 
5428  for (label zoneI = 0; zoneI < nZones; zoneI++)
5429  {
5430  Info<< " " << zoneI << "\t" << nPosOrientation[zoneI]
5431  << endl;
5432  }
5433  Info<< endl;
5434 
5435 
5436  // Re-apply with new counts (in nPosOrientation). This will cause
5437  // zones with a negative count to be flipped.
5438  consistentOrientation
5439  (
5440  isMasterFace,
5441  patch,
5442  nMasterFacesPerEdge,
5443  faceToConnectedZone,
5444  nPosOrientation,
5445 
5446  meshFlipMap
5447  );
5448  }
5449  }
5450 
5451 
5452 
5453 
5454  // Topochange container
5455  polyTopoChange meshMod(mesh_);
5456 
5457  // Insert changes to put cells and faces into zone
5458  zonify
5459  (
5460  isMasterFace,
5461  cellToZone,
5462  neiCellZone,
5463  faceToZone,
5464  meshFlipMap,
5465  meshMod
5466  );
5467 
5468  // Change the mesh (no inflation, parallel sync)
5469  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
5470 
5471  // Update fields
5472  mesh_.updateMesh(map());
5473 
5474  // Move mesh if in inflation mode
5475  if (map().hasMotionPoints())
5476  {
5477  mesh_.movePoints(map().preMotionPoints());
5478  }
5479  else
5480  {
5481  // Delete mesh volumes.
5482  mesh_.clearOut();
5483  }
5484 
5485  // Reset the instance for if in overwrite mode
5486  mesh_.setInstance(timeName());
5487 
5488  // Print some stats (note: zones are synchronised)
5489  if (mesh_.cellZones().size() > 0)
5490  {
5491  Info<< "CellZones:" << endl;
5492  forAll(mesh_.cellZones(), zoneI)
5493  {
5494  const cellZone& cz = mesh_.cellZones()[zoneI];
5495  Info<< " " << cz.name()
5496  << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
5497  << endl;
5498  }
5499  Info<< endl;
5500  }
5501  if (mesh_.faceZones().size() > 0)
5502  {
5503  Info<< "FaceZones:" << endl;
5504  forAll(mesh_.faceZones(), zoneI)
5505  {
5506  const faceZone& fz = mesh_.faceZones()[zoneI];
5507  Info<< " " << fz.name()
5508  << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
5509  << endl;
5510  }
5511  Info<< endl;
5512  }
5513 
5514  // None of the faces has changed, only the zones. Still...
5515  updateMesh(map(), labelList());
5516 
5517  return map;
5518 }
5519 
5520 
5521 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::Pstream::mapCombineGather
static void mapCombineGather(const List< commsStruct > &comms, Container &Values, const CombineOp &cop, const int tag, const label comm)
Definition: combineGatherScatter.C:551
nZones
label nZones
Definition: interpolatedFaces.H:24
shellSurfaces.H
Foam::Pair::second
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:122
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
Foam::setf
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:169
volFields.H
Foam::meshRefinement::mergePoints
autoPtr< mapPolyMesh > mergePoints(const labelList &pointToDuplicate)
Merge duplicate points.
Definition: meshRefinementBaffles.C:4867
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::HashTable::size
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::reverse
void reverse(UList< T > &list, const label n)
Definition: UListI.H:396
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
meshTools.H
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::IOobject::AUTO_WRITE
Definition: IOobject.H:129
Foam::localPointRegion
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
Definition: localPointRegion.H:69
Foam::dimless
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Dimensionless.
Definition: dimensionSets.H:50
Foam::VectorSpace< Vector< Cmpt >, Cmpt, 3 >::one
static const Vector< Cmpt > one
Definition: VectorSpace.H:116
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::surfaceZonesInfo::faceZoneType
faceZoneType
What to do with faceZone faces.
Definition: surfaceZonesInfo.H:86
Foam::faceMap
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Definition: blockMeshMergeTopological.C:94
Foam::IOobject::name
const word & name() const
Return name.
Definition: IOobjectI.H:70
Foam::meshRefinement::createBaffles
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
Definition: meshRefinementBaffles.C:484
Foam::labelMax
constexpr label labelMax
Definition: label.H:61
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:62
Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have 'insidePoint'.
Definition: surfaceZonesInfo.C:401
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:63
Foam::OBJstream
OFstream that keeps track of vertices.
Definition: OBJstream.H:58
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Definition: polyBoundaryMesh.H:62
Foam::HashTable::toc
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:121
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::polyTopoChange::changeMesh
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Definition: polyTopoChange.C:2959
Foam::DynamicList< label >
polyAddFace.H
Foam::checkData
label checkData(const fvMesh &mesh, const instantList &timeDirs, wordList &objectNames)
Check if fields are good to use (available at all times)
Foam::primitiveMesh::nFaces
label nFaces() const
Number of mesh faces.
Definition: primitiveMeshI.H:90
polyRemoveFace.H
Foam::polyPatch::coupled
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:329
polyTopoChange.H
Foam::meshRefinement::getZones
labelList getZones(const List< surfaceZonesInfo::faceZoneType > &fzTypes) const
Get zones of given type.
Definition: meshRefinementBaffles.C:659
localPointRegion.H
Foam::polyTopoChange
Direct mesh changes based on v1.3 polyTopoChange syntax.
Definition: polyTopoChange.H:99
Foam::HashTable::insert
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
Foam::Map
A HashTable to objects of type <T> with a label key.
Definition: lumpedPointController.H:69
unitConversion.H
Unit conversion functions.
Foam::syncTools::swapBoundaryFaceList
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:439
Foam::faceSet
A list of face labels.
Definition: faceSet.H:51
Foam::orEqOp
Definition: ops.H:86
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:69
Foam::surfaceZonesInfo::getNamedSurfaces
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
Definition: surfaceZonesInfo.C:263
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:182
Foam::polyModifyFace
Class describing modification of a face.
Definition: polyModifyFace.H:50
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
Foam::mapPolyMesh::preMotionPoints
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:613
Foam::syncTools::getMasterFaces
static bitSet getMasterFaces(const polyMesh &mesh)
Definition: syncTools.C:126
removeCells.H
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::cellZoneMesh
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
Definition: cellZoneMeshFwd.H:44
Foam::bitSet::test
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:520
Foam::DynamicID< faceZoneMesh >
Foam::OBJstream::write
virtual Ostream & write(const char c)
Write character.
Definition: OBJstream.C:78
Foam::meshRefinement::removeLimitShells
autoPtr< mapPolyMesh > removeLimitShells(const label nBufferLayers, const label nErodeCellZones, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList &regionsInMesh)
Remove cells from limitRegions if level -1.
Definition: meshRefinementBaffles.C:4714
syncTools.H
Foam::polyBoundaryMesh::start
label start() const
The start label of the boundary faces in the polyMesh face list.
Definition: polyBoundaryMesh.C:638
zoneIDs
const labelIOList & zoneIDs
Definition: correctPhi.H:59
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::duplicatePoints
Duplicate points.
Definition: duplicatePoints.H:59
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition: syncToolsTemplates.C:724
Foam::sumOp
Definition: ops.H:213
Foam::orientedSurface::NOFLIP
Definition: orientedSurface.H:68
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::meshRefinement::subsetBaffles
static List< labelPair > subsetBaffles(const polyMesh &mesh, const labelList &zoneIDs, const List< labelPair > &baffles)
Subset baffles according to zones.
Definition: meshRefinementBaffles.C:686
Foam::cellZone
A subset of mesh cells.
Definition: cellZone.H:62
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::polyMesh::faceZones
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:486
duplicatePoints.H
Foam::labelPair
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:54
Foam::wordList
List< word > wordList
A List of words.
Definition: fileName.H:59
Foam::faceSet::sync
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
Definition: faceSet.C:130
Foam::surfaceZonesInfo::getStandaloneNamedSurfaces
static labelList getStandaloneNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces without a cellZone.
Definition: surfaceZonesInfo.C:288
Foam::faceZoneMesh
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
Definition: faceZoneMeshFwd.H:44
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::regionSide
Determines the 'side' for every face and connected to a singly-connected (through edges) region of fa...
Definition: regionSide.H:63
Foam::snapParameters
Simple container to keep together snap specific information.
Definition: snapParameters.H:52
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::syncTools::syncFaceList
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:390
PatchEdgeFaceWave.H
Foam::surfaceZonesInfo::BOUNDARY
Definition: surfaceZonesInfo.H:90
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
Foam::Field< vector >
Foam::faceZone
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:65
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:68
refinementSurfaces.H
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
Foam::OSstream::name
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:107
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:474
Foam::syncTools::syncEdgeList
static void syncEdgeList(const polyMesh &mesh, List< T > &edgeValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
Definition: syncToolsTemplates.C:803
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1107
faceSet.H
Foam::List::transfer
void transfer(List< T > &list)
Definition: List.C:459
regionSplit.H
Foam::IndirectList
A List with indirect addressing.
Definition: IndirectList.H:56
Foam::ZoneMesh< faceZone, polyMesh >
IOmanip.H
Istream and Ostream manipulators taking arguments.
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:43
Foam::polyBoundaryMesh::whichPatch
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Definition: polyBoundaryMesh.C:810
Foam::regionSplit
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:140
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: List.H:62
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::regionSplit::nRegions
label nRegions() const
Return total number of regions.
Definition: regionSplit.H:294
Foam::orientedSurface::UNVISITED
Definition: orientedSurface.H:66
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::syncTools::swapBoundaryCellList
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
Definition: syncToolsTemplates.C:1152
indirectPrimitivePatch.H
Foam::ZoneMesh::whichZone
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:315
timeName
word timeName
Definition: getTimeIndex.H:3
Foam::meshRefinement::mergeZoneBaffles
autoPtr< mapPolyMesh > mergeZoneBaffles(const bool doInternalZones, const bool doBaffleZones)
Merge all baffles on faceZones.
Definition: meshRefinementBaffles.C:1354
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:121
Foam::eqOp
Definition: ops.H:71
Foam::Pair::flip
void flip()
Flip the Pair in-place.
Definition: PairI.H:160
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
meshRefinement.H
wallPoints.H
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam::writer< scalar >
zoneID
const labelIOList & zoneID
Definition: interpolatedFaces.H:22
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
fvMesh.H
edgeTopoDistanceData.H
Foam::faceZone::whichFace
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:379
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::zone::name
const word & name() const
Return name.
Definition: zone.H:158
Foam::setw
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
polyMeshAdder.H
Foam::indent
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:320
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::HashTable::sortedToc
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition: HashTable.C:136
Foam::polyPatch::start
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:313
Foam::maxEqOp
Definition: ops.H:80
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:48
polyModifyFace.H
Foam::surfaceZonesInfo::getClosedNamedSurfaces
static labelList getClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed and.
Definition: surfaceZonesInfo.C:314
Foam::meshRefinement::getZoneFaces
void getZoneFaces(const labelList &zoneIDs, labelList &faceZoneID, labelList &ownPatch, labelList &neiPatch, labelList &nBaffles) const
Get per-face information (faceZone, master/slave patch)
Definition: meshRefinementBaffles.C:719
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::duplicatePoints::setRefinement
void setRefinement(const localPointRegion &regionSide, polyTopoChange &)
Play commands into polyTopoChange to duplicate points. Gets.
Definition: duplicatePoints.C:59
Foam::surfaceZonesInfo::getUnnamedSurfaces
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Definition: surfaceZonesInfo.C:242
Foam::HashTable
A HashTable similar to std::unordered_map.
Definition: HashTable.H:105
Foam::mapPolyMesh::reverseFaceMap
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:501
Foam::indirectPrimitivePatch
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
Definition: indirectPrimitivePatch.H:49
Foam::Pstream::listCombineGather
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
Definition: combineGatherScatter.C:290
Foam::surfaceZonesInfo::BAFFLE
Definition: surfaceZonesInfo.H:89
Foam::surfaceZonesInfo::faceZoneTypeNames
static const Enum< faceZoneType > faceZoneTypeNames
Definition: surfaceZonesInfo.H:93
Time.H
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::meshRefinement::baffleAndSplitMesh
void baffleAndSplitMesh(const bool handleSnapProblems, const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const label nErodeCellZones, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList &regionsInMesh, const pointField &locationsOutsideMesh, const writer< scalar > &leakPathFormatter)
Split off unreachable areas of mesh.
Definition: meshRefinementBaffles.C:4131
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::wordPairHashTable
HashTable< word, Pair< word >, FixedList< word, 2 >::Hash<> > wordPairHashTable
HashTable of Pair<word>
Definition: wordPairHashTable.H:51
Foam::meshRefinement::dupNonManifoldBoundaryPoints
autoPtr< mapPolyMesh > dupNonManifoldBoundaryPoints()
Find boundary points that are on faceZones of type boundary.
Definition: meshRefinementBaffles.C:4936
Foam::meshRefinement::mergeFreeStandingBaffles
void mergeFreeStandingBaffles(const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const scalar planarAngle, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const pointField &locationsOutsideMesh, const writer< scalar > &leakPathFormatter)
Merge free-standing baffles.
Definition: meshRefinementBaffles.C:4314
Foam::linePointRef
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1094
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:381
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::polyTopoChange::setAction
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
Definition: polyTopoChange.C:2460
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:464
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::Pair< label >
Foam::Time::path
fileName path() const
Return path.
Definition: Time.H:358
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
Foam::Pstream::mapCombineScatter
static void mapCombineScatter(const List< commsStruct > &comms, Container &Values, const int tag, const label comm)
Scatter data. Reverse of combineGather.
Definition: combineGatherScatter.C:666
Foam::meshRefinement::mergeBaffles
autoPtr< mapPolyMesh > mergeBaffles(const List< labelPair > &, const Map< label > &faceToPatch)
Merge baffles. Gets pairs of faces and boundary faces to move.
Definition: meshRefinementBaffles.C:1174
f
labelList f(nPoints)
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
refinementParameters.H
Foam::polyRemoveFace
Class containing data for face removal.
Definition: polyRemoveFace.H:48
Foam::meshRefinement::writeType
writeType
Enumeration for what to write. Used as a bit-pattern.
Definition: meshRefinement.H:112
Foam::List< label >
Foam::PackedList< 2 >
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:102
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:77
FaceCellWave.H
Foam::meshRefinement::createZoneBaffles
autoPtr< mapPolyMesh > createZoneBaffles(const labelList &zoneIDs, List< labelPair > &baffles, labelList &originatingFaceZone)
Create baffles for faces on faceZones. Return created baffles.
Definition: meshRefinementBaffles.C:780
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:115
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:181
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
Foam::mapPolyMesh::faceMap
const labelList & faceMap() const
Old face map.
Definition: mapPolyMesh.H:410
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:325
Foam::vtk::write
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
Definition: foamVtkOutputTemplates.C:35
Foam::fieldTypes::area
const wordList area
Standard area field types (scalar, vector, tensor, etc)
Foam::mapPolyMesh
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:161
Foam::mapPolyMesh::hasMotionPoints
bool hasMotionPoints() const
Has valid preMotionPoints?
Definition: mapPolyMesh.H:619
patchFaceOrientation.H
insidePoints
insidePoints((1 2 3))
Foam::orientedSurface::FLIP
Definition: orientedSurface.H:67
Foam::cpuTimeCxx::cpuTimeIncrement
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTimeCxx.C:75
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
Foam::surfaceZonesInfo::addFaceZonesToMesh
static labelListList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
Definition: surfaceZonesInfo.C:543
Foam::plusEqOp
Definition: ops.H:72
Foam::Pstream::listCombineScatter
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
Definition: combineGatherScatter.C:432
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::List::setSize
void setSize(const label newSize)
Alias for resize(const label)
Definition: ListI.H:146
Foam::meshRefinement::dupNonManifoldPoints
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
Definition: meshRefinementBaffles.C:4857
Foam::IOobject::NO_READ
Definition: IOobject.H:123
Foam::mkDir
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
Definition: MSwindows.C:507
Foam::faceZone::flipMap
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:271
zeroGradientFvPatchFields.H
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:303
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
Foam::polyMeshAdder::mergePoints
static void mergePoints(const polyMesh &, const Map< label > &pointToMaster, polyTopoChange &meshMod)
Helper: Merge points.
Definition: polyMeshAdder.C:2212
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
OBJstream.H
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1113
Foam::meshRefinement::MESH
Definition: meshRefinement.H:94
Foam::localPointRegion::findDuplicateFacePairs
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
Definition: localPointRegion.C:625
polyModifyCell.H
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:85
Foam::surfaceZonesInfo::addCellZonesToMesh
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
Definition: surfaceZonesInfo.C:458
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58
Foam::meshRefinement::debugType
debugType
Enumeration for what to debug. Used as a bit-pattern.
Definition: meshRefinement.H:92
Foam::surfaceZonesInfo::addFaceZone
static label addFaceZone(const word &name, const labelList &addressing, const boolList &flipMap, polyMesh &mesh)
Definition: surfaceZonesInfo.C:510
Foam::surfaceZonesInfo::INTERNAL
Definition: surfaceZonesInfo.H:88
Foam::faceZoneID
DynamicID< faceZoneMesh > faceZoneID
Foam::faceZoneID.
Definition: ZoneIDs.H:46