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