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-2021 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 
2412  const bitSet isBlockedFace(mesh_.nFaces());
2413  wallPoints::trackData td(isBlockedFace);
2414  FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
2415  (
2416  mesh_,
2417  changedFaces,
2418  faceDist,
2419  allFaceInfo,
2420  allCellInfo,
2421  0, // max iterations
2422  td
2423  );
2424  wallDistCalc.iterate(nGrowCellZones);
2425 
2426 
2427  // Check for cells which are within nGrowCellZones of two cellZones or
2428  // one cellZone and a wall
2429  // TBD. Currently only one layer of growth handled.
2430 
2431  bitSet isChangedCell(mesh_.nCells());
2432 
2433  forAll(allCellInfo, celli)
2434  {
2435  if (allCellInfo[celli].valid(wallDistCalc.data()))
2436  {
2437  const List<FixedList<label, 3>>& surfaces =
2438  allCellInfo[celli].surface();
2439 
2440  if (surfaces.size())
2441  {
2442  // Cell close to cellZone. Remove any free-standing baffles.
2443  // Done by marking as changed cell. Wip.
2444  isChangedCell.set(celli);
2445  }
2446 
2447  if (surfaces.size() > 1)
2448  {
2449  // Check if inbetween two cellZones or cellZone and wall
2450  // Find 'nearest' other cellZone
2451  scalar minDistSqr = Foam::sqr(GREAT);
2452  label minZone = -1;
2453  for (label i = 0; i < surfaces.size(); i++)
2454  {
2455  const label zonei = surfaces[i][0];
2456  const scalar d2 = allCellInfo[celli].distSqr()[i];
2457  if (zonei != cellToZone[celli] && d2 < minDistSqr)
2458  {
2459  minDistSqr = d2;
2460  minZone = zonei; // zoneID
2461  }
2462  }
2463 
2464  if (minDistSqr < Foam::sqr(GREAT))
2465  {
2466  if (minZone != cellToZone[celli] && minZone != wallTag[0])
2467  {
2468  cellToZone[celli] = minZone;
2469  isChangedCell.set(celli);
2470  }
2471  }
2472  }
2473  }
2474  }
2475 
2476  // Fix any left-over unvisited cells
2477  if (backgroundZoneID != -2)
2478  {
2479  forAll(cellToZone, celli)
2480  {
2481  if (cellToZone[celli] == -2)
2482  {
2483  cellToZone[celli] = backgroundZoneID;
2484  }
2485  }
2486  }
2487 
2488 
2489  // Make sure to unset faces of changed cell
2490 
2491  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2492 
2493 
2494  label nUnnamed = 0;
2495  label nNamed = 0;
2496  for (const label celli : isChangedCell)
2497  {
2498  const cell& cFaces = mesh_.cells()[celli];
2499  for (const label facei : cFaces)
2500  {
2501  const label own = mesh_.faceOwner()[facei];
2502  const label ownZone = cellToZone[own];
2503 
2504  label nbrZone = -1;
2505  if (mesh_.isInternalFace(facei))
2506  {
2507  const label neiZone = cellToZone[mesh_.faceNeighbour()[facei]];
2508  nbrZone = (own != celli ? ownZone : neiZone);
2509  }
2510  else
2511  {
2512  nbrZone = neiCellZone[facei-mesh_.nInternalFaces()];
2513  }
2514 
2515  if (nbrZone == cellToZone[celli])
2516  {
2517  if
2518  (
2519  unnamedSurfaceRegion1[facei] != -1
2520  || unnamedSurfaceRegion2[facei] != -1
2521  )
2522  {
2523  unnamedSurfaceRegion1[facei] = -1;
2524  unnamedSurfaceRegion2[facei] = -1;
2525  nUnnamed++;
2526  }
2527  if
2528  (
2529  namedSurfaceRegion.size()
2530  && namedSurfaceRegion[facei] != -1
2531  )
2532  {
2533  namedSurfaceRegion[facei] = -1;
2534  nNamed++;
2535  }
2536  }
2537  }
2538  }
2539 
2540  reduce(nUnnamed, sumOp<label>());
2541  reduce(nNamed, sumOp<label>());
2542 
2543  // Do always; might bypass if nNamed,nUnnamed zero
2544  syncTools::syncFaceList
2545  (
2546  mesh_,
2547  unnamedSurfaceRegion1,
2548  maxEqOp<label>()
2549  );
2550  syncTools::syncFaceList
2551  (
2552  mesh_,
2553  unnamedSurfaceRegion2,
2554  maxEqOp<label>()
2555  );
2556  if (namedSurfaceRegion.size())
2557  {
2558  syncTools::syncFaceList
2559  (
2560  mesh_,
2561  namedSurfaceRegion,
2562  maxEqOp<label>()
2563  );
2564  }
2565 
2566  if (debug)
2567  {
2568  Pout<< "growCellZone : grown cellZones by "
2569  << returnReduce(isChangedCell.count(), sumOp<label>())
2570  << " cells (moved from background to nearest cellZone)"
2571  << endl;
2572  Pout<< "growCellZone : unmarked " << nUnnamed
2573  << " unzoned intersections; " << nNamed << " zoned intersections; "
2574  << endl;
2575  }
2576 }
2577 
2578 
2579 void Foam::meshRefinement::makeConsistentFaceIndex
2580 (
2581  const labelList& surfaceMap,
2582  const labelList& cellToZone,
2583  labelList& namedSurfaceRegion
2584 ) const
2585 {
2586  // Make namedSurfaceRegion consistent with cellToZone - clear out any
2587  // blocked faces inbetween same cell zone (or background (=-1))
2588  // Do not do this for surfaces relating to 'pure' faceZones i.e.
2589  // faceZones without a cellZone. Note that we cannot check here
2590  // for different cellZones on either side but no namedSurfaceRegion
2591  // since cellZones can now originate from locationsInMesh as well
2592  // (instead of only through named surfaces)
2593 
2594  const labelList& faceOwner = mesh_.faceOwner();
2595  const labelList& faceNeighbour = mesh_.faceNeighbour();
2596 
2597  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2598  {
2599  label ownZone = cellToZone[faceOwner[faceI]];
2600  label neiZone = cellToZone[faceNeighbour[faceI]];
2601  label globalI = namedSurfaceRegion[faceI];
2602 
2603  if
2604  (
2605  ownZone == neiZone
2606  && globalI != -1
2607  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2608  )
2609  {
2610  namedSurfaceRegion[faceI] = -1;
2611  }
2612  }
2613 
2614  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2615 
2616  // Get coupled neighbour cellZone
2617  labelList neiCellZone;
2618  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
2619 
2620  // Use coupled cellZone to do check
2621  forAll(patches, patchI)
2622  {
2623  const polyPatch& pp = patches[patchI];
2624 
2625  if (pp.coupled())
2626  {
2627  forAll(pp, i)
2628  {
2629  label faceI = pp.start()+i;
2630 
2631  label ownZone = cellToZone[faceOwner[faceI]];
2632  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
2633  label globalI = namedSurfaceRegion[faceI];
2634 
2635  if
2636  (
2637  ownZone == neiZone
2638  && globalI != -1
2639  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2640  )
2641  {
2642  namedSurfaceRegion[faceI] = -1;
2643  }
2644  }
2645  }
2646  else
2647  {
2648  // Unzonify boundary faces
2649  forAll(pp, i)
2650  {
2651  label faceI = pp.start()+i;
2652  label globalI = namedSurfaceRegion[faceI];
2653 
2654  if
2655  (
2656  globalI != -1
2657  && surfaceMap[surfaces_.whichSurface(globalI).first()] == -1
2658  )
2659  {
2660  namedSurfaceRegion[faceI] = -1;
2661  }
2662  }
2663  }
2664  }
2665 }
2666 
2667 
2668 void Foam::meshRefinement::getIntersections
2669 (
2670  const labelList& surfacesToTest,
2671  const pointField& neiCc,
2672  const labelList& testFaces,
2673 
2674  labelList& namedSurfaceRegion,
2675  bitSet& posOrientation
2676 ) const
2677 {
2678  namedSurfaceRegion.setSize(mesh_.nFaces());
2679  namedSurfaceRegion = -1;
2680 
2681  posOrientation.setSize(mesh_.nFaces());
2682  posOrientation = false;
2683 
2684  // Statistics: number of faces per faceZone
2685  labelList nSurfFaces(surfaces_.surfZones().size(), Zero);
2686 
2687  // Collect segments
2688  // ~~~~~~~~~~~~~~~~
2689 
2690  pointField start(testFaces.size());
2691  pointField end(testFaces.size());
2692 
2693  {
2694  labelList minLevel;
2695  calcCellCellRays
2696  (
2697  neiCc,
2698  labelList(neiCc.size(), -1),
2699  testFaces,
2700  start,
2701  end,
2702  minLevel
2703  );
2704  }
2705 
2706  // Do test for intersections
2707  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2708  // Note that we intersect all intersected faces again. Could reuse
2709  // the information already in surfaceIndex_.
2710 
2711  labelList surface1;
2712  labelList region1;
2713  List<pointIndexHit> hit1;
2714  vectorField normal1;
2715  labelList surface2;
2716  labelList region2;
2717  List<pointIndexHit> hit2;
2718  vectorField normal2;
2719 
2720  surfaces_.findNearestIntersection
2721  (
2722  surfacesToTest,
2723  start,
2724  end,
2725 
2726  surface1,
2727  hit1,
2728  region1,
2729  normal1,
2730 
2731  surface2,
2732  hit2,
2733  region2,
2734  normal2
2735  );
2736 
2737  forAll(testFaces, i)
2738  {
2739  label faceI = testFaces[i];
2740  const vector& area = mesh_.faceAreas()[faceI];
2741 
2742  if (surface1[i] != -1)
2743  {
2744  // If both hit should probably choose 'nearest'
2745  if
2746  (
2747  surface2[i] != -1
2748  && (
2749  magSqr(hit2[i].hitPoint())
2750  < magSqr(hit1[i].hitPoint())
2751  )
2752  )
2753  {
2754  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2755  (
2756  surface2[i],
2757  region2[i]
2758  );
2759  posOrientation.set(faceI, ((area&normal2[i]) > 0));
2760  nSurfFaces[surface2[i]]++;
2761  }
2762  else
2763  {
2764  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2765  (
2766  surface1[i],
2767  region1[i]
2768  );
2769  posOrientation.set(faceI, ((area&normal1[i]) > 0));
2770  nSurfFaces[surface1[i]]++;
2771  }
2772  }
2773  else if (surface2[i] != -1)
2774  {
2775  namedSurfaceRegion[faceI] = surfaces_.globalRegion
2776  (
2777  surface2[i],
2778  region2[i]
2779  );
2780  posOrientation.set(faceI, ((area&normal2[i]) > 0));
2781  nSurfFaces[surface2[i]]++;
2782  }
2783  }
2784 
2785 
2786  // surfaceIndex might have different surfaces on both sides if
2787  // there happen to be a (obviously thin) surface with different
2788  // regions between the cell centres. If one is on a named surface
2789  // and the other is not this might give problems so sync.
2790  syncTools::syncFaceList
2791  (
2792  mesh_,
2793  namedSurfaceRegion,
2794  maxEqOp<label>()
2795  );
2796 
2797  // Print a bit
2798  if (debug)
2799  {
2800  forAll(nSurfFaces, surfI)
2801  {
2802  Pout<< "Surface:"
2803  << surfaces_.names()[surfI]
2804  << " nZoneFaces:" << nSurfFaces[surfI] << nl;
2805  }
2806  Pout<< endl;
2807  }
2808 }
2809 
2810 
2811 void Foam::meshRefinement::zonify
2812 (
2813  const bool allowFreeStandingZoneFaces,
2814  const label nErodeCellZones,
2815  const label backgroundZoneID,
2816  const pointField& locationsInMesh,
2817  const wordList& zonesInMesh,
2818 
2819  labelList& cellToZone,
2820  labelList& unnamedRegion1,
2821  labelList& unnamedRegion2,
2822  labelList& namedSurfaceRegion,
2823  bitSet& posOrientation
2824 ) const
2825 {
2826  // Determine zones for cells and faces
2827  // cellToZone:
2828  // -2 : unset
2829  // -1 : not in any zone (zone 'none' or background zone)
2830  // >=0 : zoneID
2831  // namedSurfaceRegion, posOrientation:
2832  // -1 : face not intersected by named surface
2833  // >=0 : globalRegion (surface+region)
2834  // (and posOrientation: surface normal v.s. face normal)
2835 
2836  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
2837 
2838  // Swap neighbouring cell centres and cell level
2839  labelList neiLevel(mesh_.nBoundaryFaces());
2840  pointField neiCc(mesh_.nBoundaryFaces());
2841  calcNeighbourData(neiLevel, neiCc);
2842 
2843 
2844  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
2845  labelList unnamedSurfaces(surfaceZonesInfo::getUnnamedSurfaces(surfZones));
2846 
2847  // Get map from surface to cellZone (or -1)
2848  labelList surfaceToCellZone;
2849  if (namedSurfaces.size())
2850  {
2851  // Get/add cellZones corresponding to surface names
2852  surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
2853  (
2854  surfZones,
2855  namedSurfaces,
2856  mesh_
2857  );
2858  }
2859 
2860 
2861  cellToZone.setSize(mesh_.nCells());
2862  cellToZone = -2;
2863 
2864  namedSurfaceRegion.clear();
2865  posOrientation.clear();
2866 
2867 
2868 
2869  // 1. Test all (unnamed & named) surfaces
2870 
2871  // Unnamed surfaces. Global surface region of intersection (or -1)
2872  getIntersections
2873  (
2874  unnamedSurfaces,
2875  neiCc,
2876  intersectedFaces(),
2877  unnamedRegion1,
2878  unnamedRegion2
2879  );
2880 
2881 
2882  if (namedSurfaces.size())
2883  {
2884  getIntersections
2885  (
2886  namedSurfaces,
2887  neiCc,
2888  intersectedFaces(),
2889  namedSurfaceRegion,
2890  posOrientation
2891  );
2892  }
2893 
2894 
2895  // 2. Walk from locationsInMesh. Hard set cellZones.
2896  // Note: walk through faceZones! (these might get overridden later)
2897 
2898  if (locationsInMesh.size())
2899  {
2900  Info<< "Setting cellZones according to locationsInMesh:" << endl;
2901 
2902  labelList locationsZoneIDs(zonesInMesh.size(), -1);
2903  forAll(locationsInMesh, i)
2904  {
2905  const word& name = zonesInMesh[i];
2906 
2907  Info<< "Location : " << locationsInMesh[i] << nl
2908  << " cellZone : " << name << endl;
2909 
2910  if (name != "none")
2911  {
2912  label zoneID = mesh_.cellZones().findZoneID(name);
2913  if (zoneID == -1)
2914  {
2915  FatalErrorInFunction << "problem" << abort(FatalError);
2916  }
2917  locationsZoneIDs[i] = zoneID;
2918  }
2919  }
2920  Info<< endl;
2921 
2922 
2923  // Assign cellZone according to seed points. Walk through named surfaces
2924  findCellZoneInsideWalk
2925  (
2926  locationsInMesh, // locations
2927  locationsZoneIDs, // index of cellZone (or -1)
2928  unnamedRegion1, // per face -1 (unblocked) or >= 0 (blocked)
2929  cellToZone
2930  );
2931  }
2932 
2933 
2934  // 3. Mark named-surfaces-with-insidePoint. Hard set cellZones.
2935 
2936  labelList locationSurfaces
2937  (
2938  surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
2939  );
2940 
2941  if (locationSurfaces.size())
2942  {
2943  Info<< "Found " << locationSurfaces.size()
2944  << " named surfaces with a provided inside point."
2945  << " Assigning cells inside these surfaces"
2946  << " to the corresponding cellZone."
2947  << nl << endl;
2948 
2949  // Collect per surface the -insidePoint -cellZone
2950  pointField insidePoints(locationSurfaces.size());
2951  labelList insidePointCellZoneIDs(locationSurfaces.size(), -1);
2952  forAll(locationSurfaces, i)
2953  {
2954  label surfI = locationSurfaces[i];
2955  insidePoints[i] = surfZones[surfI].zoneInsidePoint();
2956 
2957  const word& name = surfZones[surfI].cellZoneName();
2958  if (name != "none")
2959  {
2960  label zoneID = mesh_.cellZones().findZoneID(name);
2961  if (zoneID == -1)
2962  {
2963  FatalErrorInFunction
2964  << "problem"
2965  << abort(FatalError);
2966  }
2967  insidePointCellZoneIDs[i] = zoneID;
2968  }
2969  }
2970 
2971 
2972  // Stop at unnamed or named surface
2973  labelList allRegion1(mesh_.nFaces(), -1);
2974  forAll(namedSurfaceRegion, faceI)
2975  {
2976  allRegion1[faceI] = max
2977  (
2978  unnamedRegion1[faceI],
2979  namedSurfaceRegion[faceI]
2980  );
2981  }
2982 
2983  findCellZoneInsideWalk
2984  (
2985  insidePoints, // locations
2986  insidePointCellZoneIDs, // index of cellZone
2987  allRegion1, // per face -1 (unblocked) or >= 0 (blocked)
2988  cellToZone
2989  );
2990  }
2991 
2992 
2993  // 4. Mark named-surfaces-with-geometric faces. Do geometric test. Soft set
2994  // cellZones. Correct through making consistent.
2995 
2996  // Closed surfaces with cellZone specified.
2997  labelList closedNamedSurfaces
2998  (
2999  surfaceZonesInfo::getClosedNamedSurfaces
3000  (
3001  surfZones,
3002  surfaces_.geometry(),
3003  surfaces_.surfaces()
3004  )
3005  );
3006 
3007  if (closedNamedSurfaces.size())
3008  {
3009  Info<< "Found " << closedNamedSurfaces.size()
3010  << " closed, named surfaces. Assigning cells in/outside"
3011  << " these surfaces to the corresponding cellZone."
3012  << nl << endl;
3013 
3014  findCellZoneGeometric
3015  (
3016  neiCc,
3017  closedNamedSurfaces, // indices of closed surfaces
3018  namedSurfaceRegion, // per face index of named surface + region
3019  surfaceToCellZone, // cell zone index per surface
3020 
3021  cellToZone
3022  );
3023  }
3024 
3025 
3026  // 5. Find any unassigned regions (from regionSplit)
3027 
3028  if (namedSurfaces.size())
3029  {
3030  if (nErodeCellZones == 0)
3031  {
3032  Info<< "Walking from known cellZones; crossing a faceZone "
3033  << "face changes cellZone" << nl << endl;
3034 
3035  // Put unassigned regions into any connected cellZone
3036  findCellZoneTopo
3037  (
3038  backgroundZoneID,
3039  pointField(0),
3040  unnamedRegion1, // Intersections with unnamed surfaces
3041  namedSurfaceRegion, // Intersections with named surfaces
3042  surfaceToCellZone,
3043  cellToZone
3044  );
3045  }
3046  else if (nErodeCellZones < 0)
3047  {
3048  Info<< "Growing cellZones by " << -nErodeCellZones
3049  << " layers" << nl << endl;
3050 
3051  growCellZone
3052  (
3053  -nErodeCellZones,
3054  backgroundZoneID,
3055  unnamedRegion1,
3056  unnamedRegion2,
3057  namedSurfaceRegion,
3058  cellToZone
3059  );
3060  }
3061  else
3062  {
3063  Info<< "Eroding cellZone cells to make these consistent with"
3064  << " faceZone faces" << nl << endl;
3065 
3066  // Erode cell zone cells (connected to an unassigned cell)
3067  // and put them into backgroundZone
3068  erodeCellZone
3069  (
3070  nErodeCellZones,
3071  backgroundZoneID,
3072  unnamedRegion1,
3073  namedSurfaceRegion,
3074  cellToZone
3075  );
3076  }
3077 
3078 
3079  // Make sure namedSurfaceRegion is unset inbetween same cell zones.
3080  if (!allowFreeStandingZoneFaces)
3081  {
3082  Info<< "Only keeping zone faces inbetween different cellZones."
3083  << nl << endl;
3084 
3085  // Surfaces with faceZone but no cellZone
3086  labelList standaloneNamedSurfaces
3087  (
3088  surfaceZonesInfo::getStandaloneNamedSurfaces
3089  (
3090  surfZones
3091  )
3092  );
3093 
3094  // Construct map from surface index to index in
3095  // standaloneNamedSurfaces (or -1)
3096  labelList surfaceMap(surfZones.size(), -1);
3097  forAll(standaloneNamedSurfaces, i)
3098  {
3099  surfaceMap[standaloneNamedSurfaces[i]] = i;
3100  }
3101 
3102  makeConsistentFaceIndex
3103  (
3104  surfaceMap,
3105  cellToZone,
3106  namedSurfaceRegion
3107  );
3108  }
3109  }
3110  else if (nErodeCellZones < 0 && gMax(cellToZone) >= 0)
3111  {
3112  // With multiple locationsInMesh and non-trivial cellZones we allow
3113  // some growing of the cellZones to avoid any background cells
3114 
3115  Info<< "Growing cellZones by " << -nErodeCellZones
3116  << " layers" << nl << endl;
3117 
3118  growCellZone
3119  (
3120  -nErodeCellZones,
3121  backgroundZoneID,
3122  unnamedRegion1,
3123  unnamedRegion2,
3124  namedSurfaceRegion, // note: potentially zero sized
3125  cellToZone
3126  );
3127 
3128  // Make sure namedSurfaceRegion is unset inbetween same cell zones.
3129  if (!allowFreeStandingZoneFaces && namedSurfaceRegion.size())
3130  {
3131  Info<< "Only keeping zone faces inbetween different cellZones."
3132  << nl << endl;
3133 
3134  // Surfaces with faceZone but no cellZone
3135  labelList standaloneNamedSurfaces
3136  (
3137  surfaceZonesInfo::getStandaloneNamedSurfaces
3138  (
3139  surfZones
3140  )
3141  );
3142 
3143  // Construct map from surface index to index in
3144  // standaloneNamedSurfaces (or -1)
3145  labelList surfaceMap(surfZones.size(), -1);
3146  forAll(standaloneNamedSurfaces, i)
3147  {
3148  surfaceMap[standaloneNamedSurfaces[i]] = i;
3149  }
3150 
3151  makeConsistentFaceIndex
3152  (
3153  surfaceMap,
3154  cellToZone,
3155  namedSurfaceRegion
3156  );
3157  }
3158  }
3159 
3160 
3161  // Some stats
3162  if (debug)
3163  {
3164  label nZones = gMax(cellToZone)+1;
3165 
3166  label nUnvisited = 0;
3167  label nBackgroundCells = 0;
3168  labelList nZoneCells(nZones, Zero);
3169  forAll(cellToZone, cellI)
3170  {
3171  label zoneI = cellToZone[cellI];
3172  if (zoneI >= 0)
3173  {
3174  nZoneCells[zoneI]++;
3175  }
3176  else if (zoneI == -1)
3177  {
3178  nBackgroundCells++;
3179  }
3180  else if (zoneI == -2)
3181  {
3182  nUnvisited++;
3183  }
3184  else
3185  {
3186  FatalErrorInFunction
3187  << "problem" << exit(FatalError);
3188  }
3189  }
3190  reduce(nUnvisited, sumOp<label>());
3191  reduce(nBackgroundCells, sumOp<label>());
3192  forAll(nZoneCells, zoneI)
3193  {
3194  reduce(nZoneCells[zoneI], sumOp<label>());
3195  }
3196  Info<< "nUnvisited :" << nUnvisited << endl;
3197  Info<< "nBackgroundCells:" << nBackgroundCells << endl;
3198  Info<< "nZoneCells :" << nZoneCells << endl;
3199  }
3200  if (debug&MESH)
3201  {
3202  const_cast<Time&>(mesh_.time())++;
3203  Pout<< "Writing cell zone allocation on mesh to time "
3204  << timeName() << endl;
3205  write
3206  (
3207  debugType(debug),
3208  writeType(writeLevel() | WRITEMESH),
3209  mesh_.time().path()/"cell2Zone"
3210  );
3211  volScalarField volCellToZone
3212  (
3213  IOobject
3214  (
3215  "cellToZone",
3216  timeName(),
3217  mesh_,
3218  IOobject::NO_READ,
3219  IOobject::AUTO_WRITE,
3220  false
3221  ),
3222  mesh_,
3223  dimensionedScalar(dimless, Zero),
3224  zeroGradientFvPatchScalarField::typeName
3225  );
3226 
3227  forAll(cellToZone, cellI)
3228  {
3229  volCellToZone[cellI] = cellToZone[cellI];
3230  }
3231  volCellToZone.write();
3232 
3233 
3234  //mkDir(mesh_.time().path()/timeName());
3235  //OBJstream str
3236  //(
3237  // mesh_.time().path()/timeName()/"zoneBoundaryFaces.obj"
3238  //);
3239  //Pout<< "Writing zone boundaries to " << str.name() << endl;
3240  //for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
3241  //{
3242  // const label ownZone = cellToZone[mesh_.faceOwner()[facei]];
3243  // const label neiZone = cellToZone[mesh_.faceNeighbour()[facei]];
3244  // if (ownZone != neiZone)
3245  // {
3246  // str.write(mesh_.faces()[facei], mesh_.points(), false);
3247  // }
3248  //}
3249  //labelList neiCellZone;
3250  //syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
3251  //for
3252  //(
3253  // label facei = mesh_.nInternalFaces();
3254  // facei < mesh_.nFaces();
3255  // facei++
3256  //)
3257  //{
3258  // const label ownZone = cellToZone[mesh_.faceOwner()[facei]];
3259  // const label bFacei = facei-mesh_.nInternalFaces();
3260  // const label neiZone = neiCellZone[bFacei];
3261  // if (ownZone != neiZone)
3262  // {
3263  // str.write(mesh_.faces()[facei], mesh_.points(), false);
3264  // }
3265  //}
3266 
3267  //mkDir(mesh_.time().path()/timeName());
3268  //OBJstream str1
3269  //(
3270  // mesh_.time().path()/timeName()/"unnamedRegion1.obj"
3271  //);
3272  //OBJstream str2
3273  //(
3274  // mesh_.time().path()/timeName()/"unnamedRegion2.obj"
3275  //);
3276  //Pout<< "Writing unnamed1 to " << str1.name() << endl;
3277  //Pout<< "Writing unnamed2 to " << str2.name() << endl;
3278  //for (label facei = 0; facei < mesh_.nFaces(); facei++)
3279  //{
3280  // if
3281  // (
3282  // unnamedRegion1[facei] < -1
3283  // || unnamedRegion2[facei] < -1
3284  // )
3285  // {
3286  // FatalErrorInFunction << "face:"
3287  // << mesh_.faceCentres()[facei]
3288  // << " unnamed1:" << unnamedRegion1[facei]
3289  // << " unnamed2:" << unnamedRegion2[facei]
3290  // << exit(FatalError);
3291  // }
3292  //
3293  // if (unnamedRegion1[facei] >= 0)
3294  // {
3295  // str1.write(mesh_.faces()[facei], mesh_.points(), false);
3296  // }
3297  //
3298  // if (unnamedRegion2[facei] >= 0)
3299  // {
3300  // str2.write(mesh_.faces()[facei], mesh_.points(), false);
3301  // }
3302  //}
3303 
3304  //if (namedSurfaceRegion.size())
3305  //{
3306  // OBJstream strNamed
3307  // (
3308  // mesh_.time().path()/timeName()/"namedSurfaceRegion.obj"
3309  // );
3310  // Pout<< "Writing named to " << strNamed.name() << endl;
3311  // for (label facei = 0; facei < mesh_.nFaces(); facei++)
3312  // {
3313  // const face& f = mesh_.faces()[facei];
3314  // if (namedSurfaceRegion[facei] < -1)
3315  // {
3316  // FatalErrorInFunction << "face:"
3317  // << mesh_.faceCentres()[facei]
3318  // << " unnamed1:" << unnamedRegion1[facei]
3319  // << " unnamed2:" << unnamedRegion2[facei]
3320  // << " named:" << namedSurfaceRegion[facei]
3321  // << exit(FatalError);
3322  // }
3323  // if (namedSurfaceRegion[facei] >= 0)
3324  // {
3325  // strNamed.write(f, mesh_.points(), false);
3326  // }
3327  // }
3328  //}
3329  }
3330 }
3331 
3332 
3333 void Foam::meshRefinement::handleSnapProblems
3334 (
3335  const snapParameters& snapParams,
3336  const bool useTopologicalSnapDetection,
3337  const bool removeEdgeConnectedCells,
3338  const scalarField& perpendicularAngle,
3339  const dictionary& motionDict,
3340  Time& runTime,
3341  const labelList& globalToMasterPatch,
3342  const labelList& globalToSlavePatch
3343 )
3344 {
3345  Info<< nl
3346  << "Introducing baffles to block off problem cells" << nl
3347  << "----------------------------------------------" << nl
3348  << endl;
3349 
3350  labelList facePatch;
3351  if (useTopologicalSnapDetection)
3352  {
3353  facePatch = markFacesOnProblemCells
3354  (
3355  motionDict,
3356  removeEdgeConnectedCells,
3357  perpendicularAngle,
3358  globalToMasterPatch
3359  );
3360  }
3361  else
3362  {
3363  facePatch = markFacesOnProblemCellsGeometric
3364  (
3365  snapParams,
3366  motionDict,
3367  globalToMasterPatch,
3368  globalToSlavePatch
3369  );
3370  }
3371  Info<< "Analyzed problem cells in = "
3372  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3373 
3374  if (debug&MESH)
3375  {
3376  faceSet problemFaces(mesh_, "problemFaces", mesh_.nFaces()/100);
3377 
3378  forAll(facePatch, faceI)
3379  {
3380  if (facePatch[faceI] != -1)
3381  {
3382  problemFaces.insert(faceI);
3383  }
3384  }
3385  problemFaces.instance() = timeName();
3386  Pout<< "Dumping " << problemFaces.size()
3387  << " problem faces to " << problemFaces.objectPath() << endl;
3388  problemFaces.write();
3389  }
3390 
3391  Info<< "Introducing baffles to delete problem cells." << nl << endl;
3392 
3393  if (debug)
3394  {
3395  ++runTime;
3396  }
3397 
3398  // Create baffles with same owner and neighbour for now.
3399  createBaffles(facePatch, facePatch);
3400 
3401  if (debug)
3402  {
3403  // Debug:test all is still synced across proc patches
3404  checkData();
3405  }
3406  Info<< "Created baffles in = "
3407  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3408 
3409  printMeshInfo(debug, "After introducing baffles");
3410 
3411  if (debug&MESH)
3412  {
3413  const_cast<Time&>(mesh_.time())++;
3414  Pout<< "Writing extra baffled mesh to time "
3415  << timeName() << endl;
3416  write
3417  (
3418  debugType(debug),
3419  writeType(writeLevel() | WRITEMESH),
3420  runTime.path()/"extraBaffles"
3421  );
3422  Pout<< "Dumped debug data in = "
3423  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
3424  }
3425 }
3426 
3427 
3428 Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
3429 (
3430  const labelList& faceToZone,
3431  const labelList& cellToZone,
3432  const labelList& neiCellZone
3433 ) const
3434 {
3435  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
3436  const labelList& faceOwner = mesh_.faceOwner();
3437  const labelList& faceNeighbour = mesh_.faceNeighbour();
3438 
3439 
3440  // We want to pick up the faces to orient. These faces come in
3441  // two variants:
3442  // - faces originating from stand-alone faceZones
3443  // (these will most likely have no cellZone on either side so
3444  // ownZone and neiZone both -1)
3445  // - sticky-up faces originating from a 'bulge' in a outside of
3446  // a cellZone. These will have the same cellZone on either side.
3447  // How to orient these is not really clearly defined so do them
3448  // same as stand-alone faceZone faces for now. (Normally these will
3449  // already have been removed by the 'allowFreeStandingZoneFaces=false'
3450  // default setting)
3451 
3452  // Note that argument neiCellZone will have -1 on uncoupled boundaries.
3453 
3454  DynamicList<label> faceLabels(mesh_.nFaces()/100);
3455 
3456  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3457  {
3458  if (faceToZone[faceI] != -1)
3459  {
3460  // Free standing baffle?
3461  label ownZone = cellToZone[faceOwner[faceI]];
3462  label neiZone = cellToZone[faceNeighbour[faceI]];
3463  if (ownZone == neiZone)
3464  {
3465  faceLabels.append(faceI);
3466  }
3467  }
3468  }
3469  forAll(patches, patchI)
3470  {
3471  const polyPatch& pp = patches[patchI];
3472 
3473  forAll(pp, i)
3474  {
3475  label faceI = pp.start()+i;
3476  if (faceToZone[faceI] != -1)
3477  {
3478  // Free standing baffle?
3479  label ownZone = cellToZone[faceOwner[faceI]];
3480  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
3481  if (ownZone == neiZone)
3482  {
3483  faceLabels.append(faceI);
3484  }
3485  }
3486  }
3487  }
3488  return faceLabels.shrink();
3489 }
3490 
3491 
3492 void Foam::meshRefinement::calcPatchNumMasterFaces
3493 (
3494  const bitSet& isMasterFace,
3495  const indirectPrimitivePatch& patch,
3496  labelList& nMasterFacesPerEdge
3497 ) const
3498 {
3499  // Number of (master)faces per edge
3500  nMasterFacesPerEdge.setSize(patch.nEdges());
3501  nMasterFacesPerEdge = 0;
3502 
3503  forAll(patch.addressing(), faceI)
3504  {
3505  const label meshFaceI = patch.addressing()[faceI];
3506 
3507  if (isMasterFace.test(meshFaceI))
3508  {
3509  const labelList& fEdges = patch.faceEdges()[faceI];
3510  forAll(fEdges, fEdgeI)
3511  {
3512  nMasterFacesPerEdge[fEdges[fEdgeI]]++;
3513  }
3514  }
3515  }
3516 
3517  syncTools::syncEdgeList
3518  (
3519  mesh_,
3520  patch.meshEdges(mesh_.edges(), mesh_.pointEdges()),
3521  nMasterFacesPerEdge,
3522  plusEqOp<label>(),
3523  label(0)
3524  );
3525 }
3526 
3527 
3528 Foam::label Foam::meshRefinement::markPatchZones
3529 (
3530  const indirectPrimitivePatch& patch,
3531  const labelList& nMasterFacesPerEdge,
3532  labelList& faceToZone
3533 ) const
3534 {
3535  List<edgeTopoDistanceData<label>> allEdgeInfo(patch.nEdges());
3536  List<edgeTopoDistanceData<label>> allFaceInfo(patch.size());
3537 
3538 
3539  // Protect all non-manifold edges
3540  {
3541  label nProtected = 0;
3542 
3543  forAll(nMasterFacesPerEdge, edgeI)
3544  {
3545  if (nMasterFacesPerEdge[edgeI] > 2)
3546  {
3547  allEdgeInfo[edgeI] = edgeTopoDistanceData<label>(0, -2);
3548  nProtected++;
3549  }
3550  }
3551  //Info<< "Protected from visiting "
3552  // << returnReduce(nProtected, sumOp<label>())
3553  // << " non-manifold edges" << nl << endl;
3554  }
3555 
3556 
3557  // Hand out zones
3558 
3559  DynamicList<label> changedEdges;
3560  DynamicList<edgeTopoDistanceData<label>> changedInfo;
3561 
3562  const scalar tol = PatchEdgeFaceWave
3563  <
3564  indirectPrimitivePatch,
3565  edgeTopoDistanceData<label>
3566  >::propagationTol();
3567 
3568  int dummyTrackData;
3569 
3570  const globalIndex globalFaces(patch.size());
3571 
3572  label faceI = 0;
3573 
3574  label currentZoneI = 0;
3575 
3576  while (true)
3577  {
3578  // Pick an unset face
3579  label globalSeed = labelMax;
3580  for (; faceI < allFaceInfo.size(); faceI++)
3581  {
3582  if (!allFaceInfo[faceI].valid(dummyTrackData))
3583  {
3584  globalSeed = globalFaces.toGlobal(faceI);
3585  break;
3586  }
3587  }
3588 
3589  reduce(globalSeed, minOp<label>());
3590 
3591  if (globalSeed == labelMax)
3592  {
3593  break;
3594  }
3595 
3596  label procI = globalFaces.whichProcID(globalSeed);
3597  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
3598 
3599  //Info<< "Seeding zone " << currentZoneI
3600  // << " from processor " << procI << " face " << seedFaceI
3601  // << endl;
3602 
3603  if (procI == Pstream::myProcNo())
3604  {
3605  edgeTopoDistanceData<label>& faceInfo = allFaceInfo[seedFaceI];
3606 
3607 
3608  // Set face
3609  faceInfo = edgeTopoDistanceData<label>(0, currentZoneI);
3610 
3611  // .. and seed its edges
3612  const labelList& fEdges = patch.faceEdges()[seedFaceI];
3613  forAll(fEdges, fEdgeI)
3614  {
3615  label edgeI = fEdges[fEdgeI];
3616 
3617  edgeTopoDistanceData<label>& edgeInfo = allEdgeInfo[edgeI];
3618 
3619  if
3620  (
3621  edgeInfo.updateEdge<int>
3622  (
3623  mesh_,
3624  patch,
3625  edgeI,
3626  seedFaceI,
3627  faceInfo,
3628  tol,
3629  dummyTrackData
3630  )
3631  )
3632  {
3633  changedEdges.append(edgeI);
3634  changedInfo.append(edgeInfo);
3635  }
3636  }
3637  }
3638 
3639 
3640  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
3641  {
3642  break;
3643  }
3644 
3645 
3646  // Walk
3647  PatchEdgeFaceWave
3648  <
3649  indirectPrimitivePatch,
3650  edgeTopoDistanceData<label>
3651  > calc
3652  (
3653  mesh_,
3654  patch,
3655  changedEdges,
3656  changedInfo,
3657  allEdgeInfo,
3658  allFaceInfo,
3659  returnReduce(patch.nEdges(), sumOp<label>())
3660  );
3661 
3662  currentZoneI++;
3663  }
3664 
3665 
3666  faceToZone.setSize(patch.size());
3667  forAll(allFaceInfo, faceI)
3668  {
3669  if (!allFaceInfo[faceI].valid(dummyTrackData))
3670  {
3671  FatalErrorInFunction
3672  << "Problem: unvisited face " << faceI
3673  << " at " << patch.faceCentres()[faceI]
3674  << exit(FatalError);
3675  }
3676  faceToZone[faceI] = allFaceInfo[faceI].data();
3677  }
3678 
3679  return currentZoneI;
3680 }
3681 
3682 
3683 void Foam::meshRefinement::consistentOrientation
3684 (
3685  const bitSet& isMasterFace,
3686  const indirectPrimitivePatch& patch,
3687  const labelList& nMasterFacesPerEdge,
3688  const labelList& faceToZone,
3689  const Map<label>& zoneToOrientation,
3690  bitSet& meshFlipMap
3691 ) const
3692 {
3693  const polyBoundaryMesh& bm = mesh_.boundaryMesh();
3694 
3695  // Data on all edges and faces
3696  List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
3697  List<patchFaceOrientation> allFaceInfo(patch.size());
3698 
3699  // Make sure we don't walk through
3700  // - slaves of coupled faces
3701  // - non-manifold edges
3702  {
3703  label nProtected = 0;
3704 
3705  forAll(patch.addressing(), faceI)
3706  {
3707  const label meshFaceI = patch.addressing()[faceI];
3708  const label patchI = bm.whichPatch(meshFaceI);
3709 
3710  if
3711  (
3712  patchI != -1
3713  && bm[patchI].coupled()
3714  && !isMasterFace.test(meshFaceI)
3715  )
3716  {
3717  // Slave side. Mark so doesn't get visited.
3718  allFaceInfo[faceI] = orientedSurface::NOFLIP;
3719  nProtected++;
3720  }
3721  }
3722  //Info<< "Protected from visiting "
3723  // << returnReduce(nProtected, sumOp<label>())
3724  // << " slaves of coupled faces" << nl << endl;
3725  }
3726  {
3727  label nProtected = 0;
3728 
3729  forAll(nMasterFacesPerEdge, edgeI)
3730  {
3731  if (nMasterFacesPerEdge[edgeI] > 2)
3732  {
3733  allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
3734  nProtected++;
3735  }
3736  }
3737 
3738  Info<< "Protected from visiting "
3739  << returnReduce(nProtected, sumOp<label>())
3740  << " non-manifold edges" << nl << endl;
3741  }
3742 
3743 
3744 
3745  DynamicList<label> changedEdges;
3746  DynamicList<patchFaceOrientation> changedInfo;
3747 
3748  const scalar tol = PatchEdgeFaceWave
3749  <
3750  indirectPrimitivePatch,
3751  patchFaceOrientation
3752  >::propagationTol();
3753 
3754  int dummyTrackData;
3755 
3756  globalIndex globalFaces(patch.size());
3757 
3758  while (true)
3759  {
3760  // Pick an unset face
3761  label globalSeed = labelMax;
3762  forAll(allFaceInfo, faceI)
3763  {
3764  if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
3765  {
3766  globalSeed = globalFaces.toGlobal(faceI);
3767  break;
3768  }
3769  }
3770 
3771  reduce(globalSeed, minOp<label>());
3772 
3773  if (globalSeed == labelMax)
3774  {
3775  break;
3776  }
3777 
3778  label procI = globalFaces.whichProcID(globalSeed);
3779  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
3780 
3781  //Info<< "Seeding from processor " << procI << " face " << seedFaceI
3782  // << endl;
3783 
3784  if (procI == Pstream::myProcNo())
3785  {
3786  // Determine orientation of seedFace
3787 
3788  patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];
3789 
3790  // Start off with correct orientation
3791  faceInfo = orientedSurface::NOFLIP;
3792 
3793  if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
3794  {
3795  faceInfo.flip();
3796  }
3797 
3798 
3799  const labelList& fEdges = patch.faceEdges()[seedFaceI];
3800  forAll(fEdges, fEdgeI)
3801  {
3802  label edgeI = fEdges[fEdgeI];
3803 
3804  patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
3805 
3806  if
3807  (
3808  edgeInfo.updateEdge<int>
3809  (
3810  mesh_,
3811  patch,
3812  edgeI,
3813  seedFaceI,
3814  faceInfo,
3815  tol,
3816  dummyTrackData
3817  )
3818  )
3819  {
3820  changedEdges.append(edgeI);
3821  changedInfo.append(edgeInfo);
3822  }
3823  }
3824  }
3825 
3826 
3827  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
3828  {
3829  break;
3830  }
3831 
3832 
3833 
3834  // Walk
3835  PatchEdgeFaceWave
3836  <
3837  indirectPrimitivePatch,
3838  patchFaceOrientation
3839  > calc
3840  (
3841  mesh_,
3842  patch,
3843  changedEdges,
3844  changedInfo,
3845  allEdgeInfo,
3846  allFaceInfo,
3847  returnReduce(patch.nEdges(), sumOp<label>())
3848  );
3849  }
3850 
3851 
3852  // Push master zone info over to slave (since slave faces never visited)
3853  {
3854  labelList neiStatus
3855  (
3856  mesh_.nBoundaryFaces(),
3857  orientedSurface::UNVISITED
3858  );
3859 
3860  forAll(patch.addressing(), i)
3861  {
3862  const label meshFaceI = patch.addressing()[i];
3863  if (!mesh_.isInternalFace(meshFaceI))
3864  {
3865  neiStatus[meshFaceI-mesh_.nInternalFaces()] =
3866  allFaceInfo[i].flipStatus();
3867  }
3868  }
3869  syncTools::swapBoundaryFaceList(mesh_, neiStatus);
3870 
3871  forAll(patch.addressing(), i)
3872  {
3873  const label meshFaceI = patch.addressing()[i];
3874  const label patchI = bm.whichPatch(meshFaceI);
3875 
3876  if
3877  (
3878  patchI != -1
3879  && bm[patchI].coupled()
3880  && !isMasterFace.test(meshFaceI)
3881  )
3882  {
3883  // Slave side. Take flipped from neighbour
3884  label bFaceI = meshFaceI-mesh_.nInternalFaces();
3885 
3886  if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
3887  {
3888  allFaceInfo[i] = orientedSurface::FLIP;
3889  }
3890  else if (neiStatus[bFaceI] == orientedSurface::FLIP)
3891  {
3892  allFaceInfo[i] = orientedSurface::NOFLIP;
3893  }
3894  else
3895  {
3896  FatalErrorInFunction
3897  << "Incorrect status for face " << meshFaceI
3898  << abort(FatalError);
3899  }
3900  }
3901  }
3902  }
3903 
3904 
3905  // Convert to meshFlipMap and adapt faceZones
3906 
3907  meshFlipMap.setSize(mesh_.nFaces());
3908  meshFlipMap = false;
3909 
3910  forAll(allFaceInfo, faceI)
3911  {
3912  label meshFaceI = patch.addressing()[faceI];
3913 
3914  if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
3915  {
3916  meshFlipMap.unset(meshFaceI);
3917  }
3918  else if (allFaceInfo[faceI] == orientedSurface::FLIP)
3919  {
3920  meshFlipMap.set(meshFaceI);
3921  }
3922  else
3923  {
3924  FatalErrorInFunction
3925  << "Problem : unvisited face " << faceI
3926  << " centre:" << mesh_.faceCentres()[meshFaceI]
3927  << abort(FatalError);
3928  }
3929  }
3930 }
3931 
3932 
3933 void Foam::meshRefinement::zonify
3934 (
3935  // Get per face whether is it master (of a coupled set of faces)
3936  const bitSet& isMasterFace,
3937  const labelList& cellToZone,
3938  const labelList& neiCellZone,
3939  const labelList& faceToZone,
3940  const bitSet& meshFlipMap,
3941  polyTopoChange& meshMod
3942 ) const
3943 {
3944  const labelList& faceOwner = mesh_.faceOwner();
3945  const labelList& faceNeighbour = mesh_.faceNeighbour();
3946 
3947  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3948  {
3949  label faceZoneI = faceToZone[faceI];
3950 
3951  if (faceZoneI != -1)
3952  {
3953  // Orient face zone to have slave cells in min cell zone.
3954  // Note: logic to use flipMap should be consistent with logic
3955  // to pick up the freeStandingBaffleFaces!
3956 
3957  label ownZone = cellToZone[faceOwner[faceI]];
3958  label neiZone = cellToZone[faceNeighbour[faceI]];
3959 
3960  bool flip;
3961 
3962  if (ownZone == neiZone)
3963  {
3964  // free-standing face. Use geometrically derived orientation
3965  flip = meshFlipMap.test(faceI);
3966  }
3967  else
3968  {
3969  flip =
3970  (
3971  ownZone == -1
3972  || (neiZone != -1 && ownZone > neiZone)
3973  );
3974  }
3975 
3976  meshMod.setAction
3977  (
3978  polyModifyFace
3979  (
3980  mesh_.faces()[faceI], // modified face
3981  faceI, // label of face
3982  faceOwner[faceI], // owner
3983  faceNeighbour[faceI], // neighbour
3984  false, // face flip
3985  -1, // patch for face
3986  false, // remove from zone
3987  faceZoneI, // zone for face
3988  flip // face flip in zone
3989  )
3990  );
3991  }
3992  }
3993 
3994 
3995  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
3996 
3997  // Set owner as no-flip
3998  forAll(patches, patchI)
3999  {
4000  const polyPatch& pp = patches[patchI];
4001 
4002  label faceI = pp.start();
4003 
4004  forAll(pp, i)
4005  {
4006  label faceZoneI = faceToZone[faceI];
4007 
4008  if (faceZoneI != -1)
4009  {
4010  label ownZone = cellToZone[faceOwner[faceI]];
4011  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
4012 
4013  bool flip;
4014 
4015  if (ownZone == neiZone)
4016  {
4017  // free-standing face. Use geometrically derived orientation
4018  flip = meshFlipMap.test(faceI);
4019  }
4020  else
4021  {
4022  flip =
4023  (
4024  ownZone == -1
4025  || (neiZone != -1 && ownZone > neiZone)
4026  );
4027  }
4028 
4029  meshMod.setAction
4030  (
4031  polyModifyFace
4032  (
4033  mesh_.faces()[faceI], // modified face
4034  faceI, // label of face
4035  faceOwner[faceI], // owner
4036  -1, // neighbour
4037  false, // face flip
4038  patchI, // patch for face
4039  false, // remove from zone
4040  faceZoneI, // zone for face
4041  flip // face flip in zone
4042  )
4043  );
4044  }
4045  faceI++;
4046  }
4047  }
4048 
4049 
4050  // Put the cells into the correct zone
4051  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4052 
4053  forAll(cellToZone, cellI)
4054  {
4055  label zoneI = cellToZone[cellI];
4056 
4057  if (zoneI >= 0)
4058  {
4059  meshMod.setAction
4060  (
4061  polyModifyCell
4062  (
4063  cellI,
4064  false, // removeFromZone
4065  zoneI
4066  )
4067  );
4068  }
4069  }
4070 }
4071 
4072 
4073 void Foam::meshRefinement::allocateInterRegionFaceZone
4074 (
4075  const label ownZone,
4076  const label neiZone,
4077  wordPairHashTable& zonesToFaceZone,
4078  LabelPairMap<word>& zoneIDsToFaceZone
4079 ) const
4080 {
4081  const cellZoneMesh& cellZones = mesh_.cellZones();
4082 
4083  if (ownZone != neiZone)
4084  {
4085  // Make sure lowest number cellZone is master. Non-cellZone
4086  // areas are slave
4087  const bool swap =
4088  (
4089  ownZone == -1
4090  || (neiZone != -1 && ownZone > neiZone)
4091  );
4092 
4093  // Quick check whether we already have pair of zones
4094  labelPair key(ownZone, neiZone);
4095  if (swap)
4096  {
4097  key.flip();
4098  }
4099 
4100  if (!zoneIDsToFaceZone.found(key))
4101  {
4102  // Not found. Allocate.
4103  const word ownZoneName =
4104  (
4105  ownZone != -1
4106  ? cellZones[ownZone].name()
4107  : "none"
4108  );
4109  const word neiZoneName =
4110  (
4111  neiZone != -1
4112  ? cellZones[neiZone].name()
4113  : "none"
4114  );
4115 
4116  // Get lowest zone first
4117  Pair<word> wordKey(ownZoneName, neiZoneName);
4118  if (swap)
4119  {
4120  wordKey.flip();
4121  }
4122 
4123  word fzName = wordKey.first() + "_to_" + wordKey.second();
4124 
4125  zoneIDsToFaceZone.insert(key, fzName);
4126  zonesToFaceZone.insert(wordKey, fzName);
4127  }
4128  }
4129 }
4130 
4131 
4132 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
4133 
4134 void Foam::meshRefinement::baffleAndSplitMesh
4135 (
4136  const bool doHandleSnapProblems,
4137  const snapParameters& snapParams,
4138  const bool useTopologicalSnapDetection,
4139  const bool removeEdgeConnectedCells,
4140  const scalarField& perpendicularAngle,
4141  const label nErodeCellZones,
4142  const dictionary& motionDict,
4143  Time& runTime,
4144  const labelList& globalToMasterPatch,
4145  const labelList& globalToSlavePatch,
4146 
4147  const pointField& locationsInMesh,
4148  const wordList& zonesInMesh,
4149  const pointField& locationsOutsideMesh,
4150  const writer<scalar>& leakPathFormatter
4151 )
4152 {
4153  // Introduce baffles
4154  // ~~~~~~~~~~~~~~~~~
4155 
4156  // Split the mesh along internal faces wherever there is a pierce between
4157  // two cell centres.
4158 
4159  Info<< "Introducing baffles for "
4160  << returnReduce(countHits(), sumOp<label>())
4161  << " faces that are intersected by the surface." << nl << endl;
4162 
4163  // Swap neighbouring cell centres and cell level
4164  labelList neiLevel(mesh_.nBoundaryFaces());
4165  pointField neiCc(mesh_.nBoundaryFaces());
4166  calcNeighbourData(neiLevel, neiCc);
4167 
4168  labelList ownPatch, neiPatch;
4169  getBafflePatches
4170  (
4171  nErodeCellZones,
4172  globalToMasterPatch,
4173 
4174  locationsInMesh,
4175  zonesInMesh,
4176 
4177  neiLevel,
4178  neiCc,
4179 
4180  ownPatch,
4181  neiPatch
4182  );
4183 
4184  createBaffles(ownPatch, neiPatch);
4185 
4186  if (debug)
4187  {
4188  // Debug:test all is still synced across proc patches
4189  checkData();
4190  }
4191 
4192  Info<< "Created baffles in = "
4193  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4194 
4195  printMeshInfo(debug, "After introducing baffles");
4196 
4197  if (debug&MESH)
4198  {
4199  const_cast<Time&>(mesh_.time())++;
4200  Pout<< "Writing baffled mesh to time " << timeName()
4201  << endl;
4202  write
4203  (
4204  debugType(debug),
4205  writeType(writeLevel() | WRITEMESH),
4206  runTime.path()/"baffles"
4207  );
4208  Pout<< "Dumped debug data in = "
4209  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4210  }
4211 
4212 
4213  // Introduce baffles to delete problem cells
4214  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4215 
4216  // Create some additional baffles where we want surface cells removed.
4217 
4218  if (doHandleSnapProblems)
4219  {
4220  handleSnapProblems
4221  (
4222  snapParams,
4223  useTopologicalSnapDetection,
4224  removeEdgeConnectedCells,
4225  perpendicularAngle,
4226  motionDict,
4227  runTime,
4228  globalToMasterPatch,
4229  globalToSlavePatch
4230  );
4231 
4232  // Removing additional cells might have created disconnected bits
4233  // so re-do the surface intersections
4234  {
4235  // Swap neighbouring cell centres and cell level
4236  neiLevel.setSize(mesh_.nBoundaryFaces());
4237  neiCc.setSize(mesh_.nBoundaryFaces());
4238  calcNeighbourData(neiLevel, neiCc);
4239 
4240  labelList ownPatch, neiPatch;
4241  getBafflePatches
4242  (
4243  nErodeCellZones,
4244  globalToMasterPatch,
4245 
4246  locationsInMesh,
4247  zonesInMesh,
4248 
4249  neiLevel,
4250  neiCc,
4251 
4252  ownPatch,
4253  neiPatch
4254  );
4255 
4256  createBaffles(ownPatch, neiPatch);
4257  }
4258 
4259  if (debug)
4260  {
4261  // Debug:test all is still synced across proc patches
4262  checkData();
4263  }
4264  }
4265 
4266 
4267  // Select part of mesh
4268  // ~~~~~~~~~~~~~~~~~~~
4269 
4270  Info<< nl
4271  << "Remove unreachable sections of mesh" << nl
4272  << "-----------------------------------" << nl
4273  << endl;
4274 
4275  if (debug)
4276  {
4277  ++runTime;
4278  }
4279 
4280  splitMeshRegions
4281  (
4282  globalToMasterPatch,
4283  globalToSlavePatch,
4284  locationsInMesh,
4285  locationsOutsideMesh,
4286  true, // Exit if any connection between inside and outside
4287  leakPathFormatter
4288  );
4289 
4290  if (debug)
4291  {
4292  // Debug:test all is still synced across proc patches
4293  checkData();
4294  }
4295  Info<< "Split mesh in = "
4296  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4297 
4298  printMeshInfo(debug, "After subsetting");
4299 
4300  if (debug&MESH)
4301  {
4302  ++runTime;
4303  Pout<< "Writing subsetted mesh to time " << timeName()
4304  << endl;
4305  write
4306  (
4307  debugType(debug),
4308  writeType(writeLevel() | WRITEMESH),
4309  runTime.path()/timeName()
4310  );
4311  Pout<< "Dumped debug data in = "
4312  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4313  }
4314 }
4315 
4316 
4317 void Foam::meshRefinement::mergeFreeStandingBaffles
4318 (
4319  const snapParameters& snapParams,
4320  const bool useTopologicalSnapDetection,
4321  const bool removeEdgeConnectedCells,
4322  const scalarField& perpendicularAngle,
4323  const scalar planarAngle,
4324  const dictionary& motionDict,
4325  Time& runTime,
4326  const labelList& globalToMasterPatch,
4327  const labelList& globalToSlavePatch,
4328  const pointField& locationsInMesh,
4329  const pointField& locationsOutsideMesh,
4330  const writer<scalar>& leakPathFormatter
4331 )
4332 {
4333  // Merge baffles
4334  // ~~~~~~~~~~~~~
4335 
4336  Info<< nl
4337  << "Merge free-standing baffles" << nl
4338  << "---------------------------" << nl
4339  << endl;
4340 
4341 
4342  // List of pairs of freestanding baffle faces.
4343  List<labelPair> couples
4344  (
4345  freeStandingBaffles // filter out freestanding baffles
4346  (
4347  localPointRegion::findDuplicateFacePairs(mesh_),
4348  planarAngle
4349  )
4350  );
4351 
4352  label nCouples = couples.size();
4353  reduce(nCouples, sumOp<label>());
4354 
4355  Info<< "Detected free-standing baffles : " << nCouples << endl;
4356 
4357  if (nCouples > 0)
4358  {
4359  // Actually merge baffles. Note: not exactly parallellized. Should
4360  // convert baffle faces into processor faces if they resulted
4361  // from them.
4362  mergeBaffles(couples, Map<label>(0));
4363 
4364  // Detect any problem cells resulting from merging of baffles
4365  // and delete them
4366  handleSnapProblems
4367  (
4368  snapParams,
4369  useTopologicalSnapDetection,
4370  removeEdgeConnectedCells,
4371  perpendicularAngle,
4372  motionDict,
4373  runTime,
4374  globalToMasterPatch,
4375  globalToSlavePatch
4376  );
4377 
4378  // Very occasionally removing a problem cell might create a disconnected
4379  // region so re-check
4380 
4381  Info<< nl
4382  << "Remove unreachable sections of mesh" << nl
4383  << "-----------------------------------" << nl
4384  << endl;
4385 
4386  if (debug)
4387  {
4388  ++runTime;
4389  }
4390 
4391  splitMeshRegions
4392  (
4393  globalToMasterPatch,
4394  globalToSlavePatch,
4395  locationsInMesh,
4396  locationsOutsideMesh,
4397  true, // Exit if any connection between inside and outside
4398  leakPathFormatter
4399  );
4400 
4401 
4402  if (debug)
4403  {
4404  // Debug:test all is still synced across proc patches
4405  checkData();
4406  }
4407  }
4408  Info<< "Merged free-standing baffles in = "
4409  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
4410 }
4411 
4412 
4413 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
4414 (
4415  const label nBufferLayers,
4416  const label nErodeCellZones,
4417  const labelList& globalToMasterPatch,
4418  const labelList& globalToSlavePatch,
4419 
4420  const pointField& locationsInMesh,
4421  const wordList& zonesInMesh,
4422  const pointField& locationsOutsideMesh,
4423  const writer<scalar>& leakPathFormatter
4424 )
4425 {
4426  // Determine patches to put intersections into
4427  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4428 
4429  // Swap neighbouring cell centres and cell level
4430  labelList neiLevel(mesh_.nBoundaryFaces());
4431  pointField neiCc(mesh_.nBoundaryFaces());
4432  calcNeighbourData(neiLevel, neiCc);
4433 
4434  // Find intersections with all unnamed(!) surfaces
4435  labelList ownPatch, neiPatch;
4436  getBafflePatches
4437  (
4438  nErodeCellZones,
4439  globalToMasterPatch,
4440 
4441  locationsInMesh,
4442  zonesInMesh,
4443 
4444  neiLevel,
4445  neiCc,
4446 
4447  ownPatch,
4448  neiPatch
4449  );
4450 
4451  // Analyse regions. Reuse regionsplit
4452  boolList blockedFace(mesh_.nFaces(), false);
4453 
4454  forAll(ownPatch, faceI)
4455  {
4456  if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
4457  {
4458  blockedFace[faceI] = true;
4459  }
4460  }
4461  syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>());
4462 
4463 
4464  regionSplit cellRegion(mesh_, blockedFace);
4465 
4466  // Set unreachable cells to -1
4467  findRegions
4468  (
4469  mesh_,
4470  mergeDistance_ * vector::one, // perturbVec
4471  locationsInMesh,
4472  locationsOutsideMesh,
4473  false, // do not exit if outside location found
4474  leakPathFormatter,
4475  cellRegion.nRegions(),
4476  cellRegion,
4477  blockedFace
4478  );
4479 
4480  return splitMesh
4481  (
4482  nBufferLayers,
4483  globalToMasterPatch,
4484  globalToSlavePatch,
4485  cellRegion,
4486  ownPatch,
4487  neiPatch
4488  );
4489 }
4490 
4491 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
4492 (
4493  const label nBufferLayers,
4494  const labelList& globalToMasterPatch,
4495  const labelList& globalToSlavePatch,
4496 
4497  labelList& cellRegion,
4498  labelList& ownPatch,
4499  labelList& neiPatch
4500 )
4501 {
4502  // Walk out nBufferlayers from region boundary
4503  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4504  // (modifies cellRegion, ownPatch)
4505  // Takes over face patch onto points and then back to faces and cells
4506  // (so cell-face-point walk)
4507 
4508  const labelList& faceOwner = mesh_.faceOwner();
4509  const labelList& faceNeighbour = mesh_.faceNeighbour();
4510 
4511  // Patch for exposed faces for lack of anything sensible.
4512  label defaultPatch = 0;
4513  if (globalToMasterPatch.size())
4514  {
4515  defaultPatch = globalToMasterPatch[0];
4516  }
4517 
4518  for (label i = 0; i < nBufferLayers; i++)
4519  {
4520  // 1. From cells (via faces) to points
4521 
4522  labelList pointBaffle(mesh_.nPoints(), -1);
4523 
4524  forAll(faceNeighbour, faceI)
4525  {
4526  const face& f = mesh_.faces()[faceI];
4527 
4528  const label ownRegion = cellRegion[faceOwner[faceI]];
4529  const label neiRegion = cellRegion[faceNeighbour[faceI]];
4530 
4531  if (ownRegion == -1 && neiRegion != -1)
4532  {
4533  // Note max(..) since possibly regionSplit might have split
4534  // off extra unreachable parts of mesh. Note: or can this only
4535  // happen for boundary faces?
4536  forAll(f, fp)
4537  {
4538  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
4539  }
4540  }
4541  else if (ownRegion != -1 && neiRegion == -1)
4542  {
4543  label newPatchI = neiPatch[faceI];
4544  if (newPatchI == -1)
4545  {
4546  newPatchI = max(defaultPatch, ownPatch[faceI]);
4547  }
4548  forAll(f, fp)
4549  {
4550  pointBaffle[f[fp]] = newPatchI;
4551  }
4552  }
4553  }
4554 
4555  labelList neiCellRegion;
4556  syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
4557  for
4558  (
4559  label faceI = mesh_.nInternalFaces();
4560  faceI < mesh_.nFaces();
4561  faceI++
4562  )
4563  {
4564  const face& f = mesh_.faces()[faceI];
4565 
4566  const label ownRegion = cellRegion[faceOwner[faceI]];
4567  const label neiRegion = neiCellRegion[faceI-mesh_.nInternalFaces()];
4568 
4569  if (ownRegion == -1 && neiRegion != -1)
4570  {
4571  forAll(f, fp)
4572  {
4573  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
4574  }
4575  }
4576  }
4577 
4578  // Sync
4579  syncTools::syncPointList
4580  (
4581  mesh_,
4582  pointBaffle,
4583  maxEqOp<label>(),
4584  label(-1) // null value
4585  );
4586 
4587 
4588  // 2. From points back to faces
4589 
4590  const labelListList& pointFaces = mesh_.pointFaces();
4591 
4592  forAll(pointFaces, pointI)
4593  {
4594  if (pointBaffle[pointI] != -1)
4595  {
4596  const labelList& pFaces = pointFaces[pointI];
4597 
4598  forAll(pFaces, pFaceI)
4599  {
4600  label faceI = pFaces[pFaceI];
4601 
4602  if (ownPatch[faceI] == -1)
4603  {
4604  ownPatch[faceI] = pointBaffle[pointI];
4605  }
4606  }
4607  }
4608  }
4609  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
4610 
4611 
4612  // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
4613 
4614  labelList newOwnPatch(ownPatch);
4615 
4616  forAll(ownPatch, faceI)
4617  {
4618  if (ownPatch[faceI] != -1)
4619  {
4620  label own = faceOwner[faceI];
4621 
4622  if (cellRegion[own] == -1)
4623  {
4624  cellRegion[own] = labelMax;
4625 
4626  const cell& ownFaces = mesh_.cells()[own];
4627  forAll(ownFaces, j)
4628  {
4629  if (ownPatch[ownFaces[j]] == -1)
4630  {
4631  newOwnPatch[ownFaces[j]] = ownPatch[faceI];
4632  }
4633  }
4634  }
4635  if (mesh_.isInternalFace(faceI))
4636  {
4637  label nei = faceNeighbour[faceI];
4638 
4639  if (cellRegion[nei] == -1)
4640  {
4641  cellRegion[nei] = labelMax;
4642 
4643  const cell& neiFaces = mesh_.cells()[nei];
4644  forAll(neiFaces, j)
4645  {
4646  if (ownPatch[neiFaces[j]] == -1)
4647  {
4648  newOwnPatch[neiFaces[j]] = ownPatch[faceI];
4649  }
4650  }
4651  }
4652  }
4653  }
4654  }
4655 
4656  ownPatch.transfer(newOwnPatch);
4657 
4658  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
4659  }
4660 
4661 
4662 
4663  // Subset
4664  // ~~~~~~
4665 
4666  // Get cells to remove
4667  DynamicList<label> cellsToRemove(mesh_.nCells());
4668  forAll(cellRegion, cellI)
4669  {
4670  if (cellRegion[cellI] == -1)
4671  {
4672  cellsToRemove.append(cellI);
4673  }
4674  }
4675  cellsToRemove.shrink();
4676 
4677  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
4678  reduce(nCellsToKeep, sumOp<label>());
4679 
4680  Info<< "Keeping all cells containing inside points" << endl
4681  << "Selected for keeping : " << nCellsToKeep << " cells." << endl;
4682 
4683 
4684  // Remove cells
4685  removeCells cellRemover(mesh_);
4686 
4687  // Pick up patches for exposed faces
4688  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
4689  labelList exposedPatches(exposedFaces.size());
4690 
4691  forAll(exposedFaces, i)
4692  {
4693  label faceI = exposedFaces[i];
4694 
4695  if (ownPatch[faceI] != -1)
4696  {
4697  exposedPatches[i] = ownPatch[faceI];
4698  }
4699  else
4700  {
4701  WarningInFunction
4702  << "For exposed face " << faceI
4703  << " fc:" << mesh_.faceCentres()[faceI]
4704  << " found no patch." << endl
4705  << " Taking patch " << defaultPatch
4706  << " instead." << endl;
4707  exposedPatches[i] = defaultPatch;
4708  }
4709  }
4710 
4711  return doRemoveCells
4712  (
4713  cellsToRemove,
4714  exposedFaces,
4715  exposedPatches,
4716  cellRemover
4717  );
4718 }
4719 
4720 
4721 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLimitShells
4722 (
4723  const label nBufferLayers,
4724  const label nErodeCellZones,
4725  const labelList& globalToMasterPatch,
4726  const labelList& globalToSlavePatch,
4727  const pointField& locationsInMesh,
4728  const wordList& zonesInMesh
4729 )
4730 {
4731  // Determine patches to put intersections into
4732  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4733 
4734  // Swap neighbouring cell centres and cell level
4735  labelList neiLevel(mesh_.nBoundaryFaces());
4736  pointField neiCc(mesh_.nBoundaryFaces());
4737  calcNeighbourData(neiLevel, neiCc);
4738 
4739  // Find intersections with all unnamed(!) surfaces
4740  labelList ownPatch, neiPatch;
4741  getBafflePatches
4742  (
4743  nErodeCellZones,
4744  globalToMasterPatch,
4745 
4746  locationsInMesh,
4747  zonesInMesh,
4748 
4749  neiLevel,
4750  neiCc,
4751 
4752  ownPatch,
4753  neiPatch
4754  );
4755 
4756 
4757  labelList cellRegion(mesh_.nCells(), Zero);
4758  // Find any cells inside a limit shell with minLevel -1
4759  labelList levelShell;
4760  limitShells_.findLevel
4761  (
4762  mesh_.cellCentres(),
4763  labelList(mesh_.nCells(), -1), // pick up only shells with -1
4764  levelShell
4765  );
4766  forAll(levelShell, celli)
4767  {
4768  if (levelShell[celli] != -1)
4769  {
4770  // Mark cell region so it gets deletec
4771  cellRegion[celli] = -1;
4772  }
4773  }
4774 
4775  autoPtr<mapPolyMesh> mapPtr = splitMesh
4776  (
4777  nBufferLayers,
4778  globalToMasterPatch,
4779  globalToSlavePatch,
4780  cellRegion,
4781  ownPatch,
4782  neiPatch
4783  );
4784 
4785  if (debug&meshRefinement::MESH)
4786  {
4787  const_cast<Time&>(mesh_.time())++;
4788  Pout<< "Writing mesh after removing limitShells"
4789  << " to time " << timeName() << endl;
4790  write
4791  (
4792  debugType(debug),
4793  writeType
4794  (
4795  writeLevel()
4796  | WRITEMESH
4797  ),
4798  mesh_.time().path()/timeName()
4799  );
4800  }
4801  return mapPtr;
4802 }
4803 
4804 
4805 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
4806 (
4807  const localPointRegion& regionSide
4808 )
4809 {
4810  // Topochange container
4811  polyTopoChange meshMod(mesh_);
4812 
4813  label nNonManifPoints = returnReduce
4814  (
4815  regionSide.meshPointMap().size(),
4816  sumOp<label>()
4817  );
4818 
4819  Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
4820  << " non-manifold points (out of "
4821  << mesh_.globalData().nTotalPoints()
4822  << ')' << endl;
4823 
4824 
4825  autoPtr<mapPolyMesh> mapPtr;
4826 
4827  if (nNonManifPoints)
4828  {
4829  // Topo change engine
4830  duplicatePoints pointDuplicator(mesh_);
4831 
4832  // Insert changes into meshMod
4833  pointDuplicator.setRefinement(regionSide, meshMod);
4834 
4835  // Change the mesh (no inflation, parallel sync)
4836  mapPtr = meshMod.changeMesh(mesh_, false, true);
4837  mapPolyMesh& map = *mapPtr;
4838 
4839  // Update fields
4840  mesh_.updateMesh(map);
4841 
4842  // Move mesh if in inflation mode
4843  if (map.hasMotionPoints())
4844  {
4845  mesh_.movePoints(map.preMotionPoints());
4846  }
4847  else
4848  {
4849  // Delete mesh volumes.
4850  mesh_.clearOut();
4851  }
4852 
4853  // Reset the instance for if in overwrite mode
4854  mesh_.setInstance(timeName());
4855 
4856  // Update intersections. Is mapping only (no faces created, positions
4857  // stay same) so no need to recalculate intersections.
4858  updateMesh(map, labelList());
4859  }
4860 
4861  return mapPtr;
4862 }
4863 
4864 
4865 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
4866 {
4867  // Analyse which points need to be duplicated
4868  localPointRegion regionSide(mesh_);
4869 
4870  return dupNonManifoldPoints(regionSide);
4871 }
4872 
4873 
4874 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergePoints
4875 (
4876  const labelList& pointToDuplicate
4877 )
4878 {
4879  label nPointPairs = 0;
4880  forAll(pointToDuplicate, pointI)
4881  {
4882  label otherPointI = pointToDuplicate[pointI];
4883  if (otherPointI != -1)
4884  {
4885  nPointPairs++;
4886  }
4887  }
4888 
4889  autoPtr<mapPolyMesh> mapPtr;
4890 
4891  if (returnReduce(nPointPairs, sumOp<label>()))
4892  {
4893  Map<label> pointToMaster(2*nPointPairs);
4894  forAll(pointToDuplicate, pointI)
4895  {
4896  label otherPointI = pointToDuplicate[pointI];
4897  if (otherPointI != -1)
4898  {
4899  // Slave point
4900  pointToMaster.insert(pointI, otherPointI);
4901  }
4902  }
4903 
4904  // Topochange container
4905  polyTopoChange meshMod(mesh_);
4906 
4907  // Insert changes
4908  polyMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
4909 
4910  // Change the mesh (no inflation, parallel sync)
4911  mapPtr = meshMod.changeMesh(mesh_, false, true);
4912  mapPolyMesh& map = *mapPtr;
4913 
4914  // Update fields
4915  mesh_.updateMesh(map);
4916 
4917  // Move mesh if in inflation mode
4918  if (map.hasMotionPoints())
4919  {
4920  mesh_.movePoints(map.preMotionPoints());
4921  }
4922  else
4923  {
4924  // Delete mesh volumes.
4925  mesh_.clearOut();
4926  }
4927 
4928  // Reset the instance for if in overwrite mode
4929  mesh_.setInstance(timeName());
4930 
4931  // Update intersections. Is mapping only (no faces created, positions
4932  // stay same) so no need to recalculate intersections.
4933  updateMesh(map, labelList());
4934  }
4935 
4936  return mapPtr;
4937 }
4938 
4939 
4940 // Duplicate points on 'boundary' zones. Do not duplicate points on
4941 // 'internal' or 'baffle' zone. Whether points are on normal patches does
4942 // not matter
4943 Foam::autoPtr<Foam::mapPolyMesh>
4944 Foam::meshRefinement::dupNonManifoldBoundaryPoints()
4945 {
4946  const labelList boundaryFaceZones
4947  (
4948  getZones
4949  (
4950  List<surfaceZonesInfo::faceZoneType>
4951  (
4952  1,
4953  surfaceZonesInfo::BOUNDARY
4954  )
4955  )
4956  );
4957  labelList internalOrBaffleFaceZones;
4958  {
4959  List<surfaceZonesInfo::faceZoneType> fzTypes(2);
4960  fzTypes[0] = surfaceZonesInfo::INTERNAL;
4961  fzTypes[1] = surfaceZonesInfo::BAFFLE;
4962  internalOrBaffleFaceZones = getZones(fzTypes);
4963  }
4964 
4965 
4966 
4967  // 0 : point used by normal, unzoned boundary faces
4968  // 1 : point used by 'boundary' zone
4969  // 2 : point used by internal/baffle zone
4970  PackedList<2> pointStatus(mesh_.nPoints(), 0u);
4971 
4972  forAll(boundaryFaceZones, j)
4973  {
4974  const faceZone& fZone = mesh_.faceZones()[boundaryFaceZones[j]];
4975  forAll(fZone, i)
4976  {
4977  const face& f = mesh_.faces()[fZone[i]];
4978  forAll(f, fp)
4979  {
4980  pointStatus[f[fp]] = max(pointStatus[f[fp]], 1u);
4981  }
4982  }
4983  }
4984  forAll(internalOrBaffleFaceZones, j)
4985  {
4986  const faceZone& fZone = mesh_.faceZones()[internalOrBaffleFaceZones[j]];
4987  forAll(fZone, i)
4988  {
4989  const face& f = mesh_.faces()[fZone[i]];
4990  forAll(f, fp)
4991  {
4992  pointStatus[f[fp]] = max(pointStatus[f[fp]], 2u);
4993  }
4994  }
4995  }
4996 
4997  syncTools::syncPointList
4998  (
4999  mesh_,
5000  pointStatus,
5001  maxEqOp<unsigned int>(), // combine op
5002  0u // null value
5003  );
5004 
5005  // Pick up points on boundary zones that are not on internal/baffle zones
5006  label n = 0;
5007  forAll(pointStatus, pointI)
5008  {
5009  if (pointStatus[pointI] == 1u)
5010  {
5011  n++;
5012  }
5013  }
5014 
5015  label globalNPoints = returnReduce(n, sumOp<label>());
5016  Info<< "Duplicating " << globalNPoints << " points on"
5017  << " faceZones of type "
5018  << surfaceZonesInfo::faceZoneTypeNames[surfaceZonesInfo::BOUNDARY]
5019  << endl;
5020 
5021  autoPtr<mapPolyMesh> map;
5022 
5023  if (globalNPoints)
5024  {
5025  labelList candidatePoints(n);
5026  n = 0;
5027  forAll(pointStatus, pointI)
5028  {
5029  if (pointStatus[pointI] == 1u)
5030  {
5031  candidatePoints[n++] = pointI;
5032  }
5033  }
5034  localPointRegion regionSide(mesh_, candidatePoints);
5035  map = dupNonManifoldPoints(regionSide);
5036  }
5037  return map;
5038 }
5039 
5040 
5041 // Zoning
5042 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
5043 (
5044  const bool allowFreeStandingZoneFaces,
5045  const label nErodeCellZones,
5046  const pointField& locationsInMesh,
5047  const wordList& zonesInMesh,
5048  wordPairHashTable& zonesToFaceZone
5049 )
5050 {
5051  if (locationsInMesh.size() != zonesInMesh.size())
5052  {
5053  FatalErrorInFunction << "problem" << abort(FatalError);
5054  }
5055 
5056  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
5057  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
5058 
5059 
5060  // Swap neighbouring cell centres and cell level
5061  labelList neiLevel(mesh_.nBoundaryFaces());
5062  pointField neiCc(mesh_.nBoundaryFaces());
5063  calcNeighbourData(neiLevel, neiCc);
5064 
5065 
5066  // Add any faceZones, cellZones originating from surface to the mesh
5067  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5068  labelList surfaceToCellZone;
5069  labelListList surfaceToFaceZones;
5070 
5071  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
5072  if (namedSurfaces.size())
5073  {
5074  Info<< "Setting cellZones according to named surfaces:" << endl;
5075  forAll(namedSurfaces, i)
5076  {
5077  label surfI = namedSurfaces[i];
5078 
5079  Info<< "Surface : " << surfaces_.names()[surfI] << nl
5080  << " faceZones : " << surfZones[surfI].faceZoneNames() << nl
5081  << " cellZone : " << surfZones[surfI].cellZoneName()
5082  << endl;
5083  }
5084  Info<< endl;
5085 
5086  // Add zones to mesh
5087  surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
5088  (
5089  surfZones,
5090  namedSurfaces,
5091  mesh_
5092  );
5093  surfaceToFaceZones = surfaceZonesInfo::addFaceZonesToMesh
5094  (
5095  surfZones,
5096  namedSurfaces,
5097  mesh_
5098  );
5099  }
5100 
5101 
5102 
5103  // Put the cells into the correct zone
5104  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5105 
5106  // Zone per cell:
5107  // -2 : unset : not allowed!
5108  // -1 : not in any zone (zone 'none')
5109  // >=0: zoneID
5110  // namedSurfaceRegion: zero sized or:
5111  // -1 : face not intersecting a named surface
5112  // >=0 : index of named surface
5113  labelList cellToZone;
5114  labelList namedSurfaceRegion;
5115  bitSet posOrientation;
5116  {
5117  labelList unnamedRegion1;
5118  labelList unnamedRegion2;
5119 
5120  zonify
5121  (
5122  allowFreeStandingZoneFaces,
5123  nErodeCellZones,// Use erosion (>0) or regionSplit to clean up
5124  -1, // Set all cells with cellToZone -2 to -1
5125  locationsInMesh,
5126  zonesInMesh,
5127 
5128  cellToZone,
5129  unnamedRegion1,
5130  unnamedRegion2,
5131  namedSurfaceRegion,
5132  posOrientation
5133  );
5134  }
5135 
5136 
5137  // Convert namedSurfaceRegion (index of named surfaces) to
5138  // actual faceZone index
5139 
5140  //- Per face index of faceZone or -1
5141  labelList faceToZone(mesh_.nFaces(), -1);
5142 
5143  forAll(namedSurfaceRegion, faceI)
5144  {
5145  //label surfI = namedSurfaceIndex[faceI];
5146  label globalI = namedSurfaceRegion[faceI];
5147  if (globalI != -1)
5148  {
5149  const labelPair spr = surfaces_.whichSurface(globalI);
5150  faceToZone[faceI] = surfaceToFaceZones[spr.first()][spr.second()];
5151  }
5152  }
5153 
5154 
5155 
5156  // Allocate and assign faceZones from cellZones
5157  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5158 
5159  {
5160  // 1. Detect inter-region face and allocate names
5161 
5162  LabelPairMap<word> zoneIDsToFaceZone;
5163 
5164  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5165  {
5166  if (faceToZone[faceI] == -1) // Not named surface
5167  {
5168  // Face not yet in a faceZone. (it might already have been
5169  // done so by a 'named' surface). Check if inbetween different
5170  // cellZones
5171  allocateInterRegionFaceZone
5172  (
5173  cellToZone[mesh_.faceOwner()[faceI]],
5174  cellToZone[mesh_.faceNeighbour()[faceI]],
5175  zonesToFaceZone,
5176  zoneIDsToFaceZone
5177  );
5178  }
5179  }
5180 
5181  labelList neiCellZone;
5182  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
5183 
5184  forAll(neiCellZone, bFaceI)
5185  {
5186  label faceI = bFaceI + mesh_.nInternalFaces();
5187  if (faceToZone[faceI] == -1)
5188  {
5189  allocateInterRegionFaceZone
5190  (
5191  cellToZone[mesh_.faceOwner()[faceI]],
5192  neiCellZone[bFaceI],
5193  zonesToFaceZone,
5194  zoneIDsToFaceZone
5195  );
5196  }
5197  }
5198 
5199 
5200  // 2.Combine faceZoneNames allocated on different processors
5201 
5202  Pstream::mapCombineGather(zonesToFaceZone, eqOp<word>());
5203  Pstream::mapCombineScatter(zonesToFaceZone);
5204 
5205 
5206  // 3. Allocate faceZones from (now synchronised) faceZoneNames
5207  // Note: the faceZoneNames contain the same data but in different
5208  // order. We could sort the contents but instead just loop
5209  // in sortedToc order.
5210 
5211  Info<< "Setting faceZones according to neighbouring cellZones:"
5212  << endl;
5213 
5214  // From cellZone indices to faceZone index
5215  LabelPairMap<label> fZoneLookup(zonesToFaceZone.size());
5216 
5217  const cellZoneMesh& cellZones = mesh_.cellZones();
5218 
5219  {
5220  List<Pair<word>> czs(zonesToFaceZone.sortedToc());
5221 
5222  forAll(czs, i)
5223  {
5224  const Pair<word>& cz = czs[i];
5225  const word& fzName = zonesToFaceZone[cz];
5226 
5227  Info<< indent<< "cellZones : "
5228  << cz[0] << ' ' << cz[1] << nl
5229  << " faceZone : " << fzName << endl;
5230 
5231  label faceZoneI = surfaceZonesInfo::addFaceZone
5232  (
5233  fzName, // name
5234  labelList(0), // addressing
5235  boolList(0), // flipMap
5236  mesh_
5237  );
5238 
5239  label cz0 = cellZones.findZoneID(cz[0]);
5240  label cz1 = cellZones.findZoneID(cz[1]);
5241 
5242  fZoneLookup.insert(labelPair(cz0, cz1), faceZoneI);
5243  }
5244  }
5245 
5246 
5247  // 4. Set faceToZone with new faceZones
5248 
5249 
5250  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
5251  {
5252  if (faceToZone[faceI] == -1)
5253  {
5254  // Face not yet in a faceZone. (it might already have been
5255  // done so by a 'named' surface). Check if inbetween different
5256  // cellZones
5257 
5258  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
5259  label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
5260  if (ownZone != neiZone)
5261  {
5262  const bool swap =
5263  (
5264  ownZone == -1
5265  || (neiZone != -1 && ownZone > neiZone)
5266  );
5267  labelPair key(ownZone, neiZone);
5268  if (swap)
5269  {
5270  key.flip();
5271  }
5272  faceToZone[faceI] = fZoneLookup[key];
5273  }
5274  }
5275  }
5276  forAll(neiCellZone, bFaceI)
5277  {
5278  label faceI = bFaceI + mesh_.nInternalFaces();
5279  if (faceToZone[faceI] == -1)
5280  {
5281  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
5282  label neiZone = neiCellZone[bFaceI];
5283  if (ownZone != neiZone)
5284  {
5285  const bool swap =
5286  (
5287  ownZone == -1
5288  || (neiZone != -1 && ownZone > neiZone)
5289  );
5290  labelPair key(ownZone, neiZone);
5291  if (swap)
5292  {
5293  key.flip();
5294  }
5295  faceToZone[faceI] = fZoneLookup[key];
5296  }
5297  }
5298  }
5299  Info<< endl;
5300  }
5301 
5302 
5303 
5304 
5305  // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
5306  labelList neiCellZone;
5307  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
5308  forAll(patches, patchI)
5309  {
5310  const polyPatch& pp = patches[patchI];
5311 
5312  if (!pp.coupled())
5313  {
5314  label bFaceI = pp.start()-mesh_.nInternalFaces();
5315  forAll(pp, i)
5316  {
5317  neiCellZone[bFaceI++] = -1;
5318  }
5319  }
5320  }
5321 
5322 
5323 
5324  // Get per face whether is it master (of a coupled set of faces)
5325  const bitSet isMasterFace(syncTools::getMasterFaces(mesh_));
5326 
5327 
5328  // faceZones
5329  // ~~~~~~~~~
5330  // Faces on faceZones come in two variants:
5331  // - faces on the outside of a cellZone. They will be oriented to
5332  // point out of the maximum cellZone.
5333  // - free-standing faces. These will be oriented according to the
5334  // local surface normal. We do this in a two step algorithm:
5335  // - do a consistent orientation
5336  // - check number of faces with consistent orientation
5337  // - if <0 flip the whole patch
5338  bitSet meshFlipMap(mesh_.nFaces(), false);
5339  {
5340  // Collect all data on zone faces without cellZones on either side.
5341  const indirectPrimitivePatch patch
5342  (
5343  IndirectList<face>
5344  (
5345  mesh_.faces(),
5346  freeStandingBaffleFaces
5347  (
5348  faceToZone,
5349  cellToZone,
5350  neiCellZone
5351  )
5352  ),
5353  mesh_.points()
5354  );
5355 
5356  label nFreeStanding = returnReduce(patch.size(), sumOp<label>());
5357  if (nFreeStanding > 0)
5358  {
5359  Info<< "Detected " << nFreeStanding << " free-standing zone faces"
5360  << endl;
5361 
5362  if (debug)
5363  {
5364  OBJstream str(mesh_.time().path()/"freeStanding.obj");
5365  str.write(patch.localFaces(), patch.localPoints(), false);
5366  }
5367 
5368 
5369  // Detect non-manifold edges
5370  labelList nMasterFacesPerEdge;
5371  calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
5372 
5373  // Mark zones. Even a single original surface might create multiple
5374  // disconnected/non-manifold-connected zones
5375  labelList faceToConnectedZone;
5376  const label nZones = markPatchZones
5377  (
5378  patch,
5379  nMasterFacesPerEdge,
5380  faceToConnectedZone
5381  );
5382 
5383  Map<label> nPosOrientation(2*nZones);
5384  for (label zoneI = 0; zoneI < nZones; zoneI++)
5385  {
5386  nPosOrientation.insert(zoneI, 0);
5387  }
5388 
5389  // Make orientations consistent in a topological way. This just
5390  // checks the first face per zone for whether nPosOrientation
5391  // is negative (which it never is at this point)
5392  consistentOrientation
5393  (
5394  isMasterFace,
5395  patch,
5396  nMasterFacesPerEdge,
5397  faceToConnectedZone,
5398  nPosOrientation,
5399 
5400  meshFlipMap
5401  );
5402 
5403  // Count per region the number of orientations (taking the new
5404  // flipMap into account)
5405  forAll(patch.addressing(), faceI)
5406  {
5407  label meshFaceI = patch.addressing()[faceI];
5408 
5409  if (isMasterFace.test(meshFaceI))
5410  {
5411  label n = 1;
5412  if
5413  (
5414  posOrientation.test(meshFaceI)
5415  == meshFlipMap.test(meshFaceI)
5416  )
5417  {
5418  n = -1;
5419  }
5420 
5421  nPosOrientation.find(faceToConnectedZone[faceI])() += n;
5422  }
5423  }
5424  Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
5425  Pstream::mapCombineScatter(nPosOrientation);
5426 
5427 
5428  Info<< "Split " << nFreeStanding << " free-standing zone faces"
5429  << " into " << nZones << " disconnected regions with size"
5430  << " (negative denotes wrong orientation) :"
5431  << endl;
5432 
5433  for (label zoneI = 0; zoneI < nZones; zoneI++)
5434  {
5435  Info<< " " << zoneI << "\t" << nPosOrientation[zoneI]
5436  << endl;
5437  }
5438  Info<< endl;
5439 
5440 
5441  // Re-apply with new counts (in nPosOrientation). This will cause
5442  // zones with a negative count to be flipped.
5443  consistentOrientation
5444  (
5445  isMasterFace,
5446  patch,
5447  nMasterFacesPerEdge,
5448  faceToConnectedZone,
5449  nPosOrientation,
5450 
5451  meshFlipMap
5452  );
5453  }
5454  }
5455 
5456 
5457 
5458 
5459  // Topochange container
5460  polyTopoChange meshMod(mesh_);
5461 
5462  // Insert changes to put cells and faces into zone
5463  zonify
5464  (
5465  isMasterFace,
5466  cellToZone,
5467  neiCellZone,
5468  faceToZone,
5469  meshFlipMap,
5470  meshMod
5471  );
5472 
5473  // Change the mesh (no inflation, parallel sync)
5474  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
5475 
5476  // Update fields
5477  mesh_.updateMesh(map());
5478 
5479  // Move mesh if in inflation mode
5480  if (map().hasMotionPoints())
5481  {
5482  mesh_.movePoints(map().preMotionPoints());
5483  }
5484  else
5485  {
5486  // Delete mesh volumes.
5487  mesh_.clearOut();
5488  }
5489 
5490  // Reset the instance for if in overwrite mode
5491  mesh_.setInstance(timeName());
5492 
5493  // Print some stats (note: zones are synchronised)
5494  if (mesh_.cellZones().size() > 0)
5495  {
5496  Info<< "CellZones:" << endl;
5497  forAll(mesh_.cellZones(), zoneI)
5498  {
5499  const cellZone& cz = mesh_.cellZones()[zoneI];
5500  Info<< " " << cz.name()
5501  << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
5502  << endl;
5503  }
5504  Info<< endl;
5505  }
5506  if (mesh_.faceZones().size() > 0)
5507  {
5508  Info<< "FaceZones:" << endl;
5509  forAll(mesh_.faceZones(), zoneI)
5510  {
5511  const faceZone& fz = mesh_.faceZones()[zoneI];
5512  Info<< " " << fz.name()
5513  << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
5514  << endl;
5515  }
5516  Info<< endl;
5517  }
5518 
5519  // None of the faces has changed, only the zones. Still...
5520  updateMesh(map(), labelList());
5521 
5522  return map;
5523 }
5524 
5525 
5526 // ************************************************************************* //
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:67
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:4875
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::syncTools::syncEdgeList
static void syncEdgeList(const polyMesh &mesh, List< T > &edgeValues, const CombineOp &cop, const T &nullValue, const TransformOp &top, const FlipOp &fop)
Synchronize values on all mesh edges.
Definition: syncToolsTemplates.C:800
Foam::reverse
void reverse(UList< T > &list, const label n)
Definition: UListI.H:449
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:194
Foam::localPointRegion
Takes mesh with 'baffles' (= boundary faces sharing points). Determines for selected points on bounda...
Definition: localPointRegion.H:69
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::one
static const Vector< scalar > 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::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:65
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::faceZone::flipMap
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
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:63
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
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)
polyRemoveFace.H
Foam::polyPatch::coupled
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:377
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::glTF::key
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
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:445
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:65
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:175
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:369
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:4722
syncTools.H
Foam::polyBoundaryMesh::start
label start() const
The start label of the boundary faces in the polyMesh face list.
Definition: polyBoundaryMesh.C:611
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:721
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)
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:62
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:396
PatchEdgeFaceWave.H
Foam::surfaceZonesInfo::BOUNDARY
Definition: surfaceZonesInfo.H:90
Foam::Field< vector >
Foam::faceZone
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
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::OSstream::name
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:511
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:456
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::polyMesh::faceZones
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:486
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
Foam::polyBoundaryMesh::whichPatch
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Definition: polyBoundaryMesh.C:812
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:59
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::polyTopoChange::changeMesh
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Definition: polyTopoChange.C:2980
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:1392
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:289
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:123
Foam::eqOp
Definition: ops.H:71
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:372
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
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:339
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:361
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::IOobject::name
const word & name() const noexcept
Return name.
Definition: IOobjectI.H:65
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:4135
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::wordPairHashTable
HashTable< word, wordPair, Foam::Hash< wordPair > > wordPairHashTable
HashTable of wordPair.
Definition: wordPairHashes.H:60
Foam::meshRefinement::dupNonManifoldBoundaryPoints
autoPtr< mapPolyMesh > dupNonManifoldBoundaryPoints()
Find boundary points that are on faceZones of type boundary.
Definition: meshRefinementBaffles.C:4944
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:4318
Foam::linePointRef
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:103
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:453
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:2481
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:463
Foam::nl
constexpr char nl
Definition: Ostream.H:404
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
Foam::zoneIdentifier::name
const word & name() const noexcept
The zone name.
Definition: zoneIdentifier.H:123
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::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:116
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
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:341
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:36
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
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
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::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:90
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
coupled
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
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::meshRefinement::dupNonManifoldPoints
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
Definition: meshRefinementBaffles.C:4865
Foam::IOobject::NO_READ
Definition: IOobject.H:188
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
zeroGradientFvPatchFields.H
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:54
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].reset(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::polyMeshAdder::mergePoints
static void mergePoints(const polyMesh &, const Map< label > &pointToMaster, polyTopoChange &meshMod)
Helper: Merge points.
Definition: polyMeshAdder.C:2198
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
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189
polyModifyCell.H
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:79
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