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

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17