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