snappyRefineDriver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2015 OpenFOAM Foundation
9  Copyright (C) 2015-2023 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "snappyRefineDriver.H"
30 #include "meshRefinement.H"
31 #include "fvMesh.H"
32 #include "Time.H"
33 #include "cellSet.H"
34 #include "syncTools.H"
35 #include "refinementParameters.H"
36 #include "refinementSurfaces.H"
37 #include "refinementFeatures.H"
38 #include "shellSurfaces.H"
39 #include "mapDistributePolyMesh.H"
40 #include "unitConversion.H"
41 #include "snapParameters.H"
42 #include "localPointRegion.H"
43 #include "IOmanip.H"
44 #include "labelVector.H"
45 #include "profiling.H"
46 #include "searchableSurfaces.H"
47 #include "fvMeshSubset.H"
48 #include "interpolationTable.H"
49 #include "snappyVoxelMeshDriver.H"
50 #include "regionSplit.H"
51 #include "removeCells.H"
52 
53 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
54 
55 namespace Foam
56 {
57  defineTypeNameAndDebug(snappyRefineDriver, 0);
58 } // End namespace Foam
59 
60 
61 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
62 
63 Foam::snappyRefineDriver::snappyRefineDriver
64 (
65  meshRefinement& meshRefiner,
66  decompositionMethod& decomposer,
67  fvMeshDistribute& distributor,
68  const labelUList& globalToMasterPatch,
69  const labelUList& globalToSlavePatch,
70  coordSetWriter& setFormatter,
71  refPtr<surfaceWriter>& surfFormatter,
72  const bool dryRun
73 )
74 :
75  meshRefiner_(meshRefiner),
76  decomposer_(decomposer),
77  distributor_(distributor),
78  globalToMasterPatch_(globalToMasterPatch),
79  globalToSlavePatch_(globalToSlavePatch),
80  setFormatter_(setFormatter),
81  surfFormatter_(surfFormatter),
82  dryRun_(dryRun)
83 {}
84 
85 
86 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
87 
88 Foam::label Foam::snappyRefineDriver::featureEdgeRefine
89 (
90  const refinementParameters& refineParams,
91  const label maxIter,
92  const label minRefine
93 )
94 {
95  if (dryRun_)
96  {
97  return 0;
98  }
99 
100  if (refineParams.minRefineCells() == -1)
101  {
102  // Special setting to be able to restart shm on meshes with inconsistent
103  // cellLevel/pointLevel
104  return 0;
105  }
106 
107  addProfiling(edge, "snappyHexMesh::refine::edge");
108  const fvMesh& mesh = meshRefiner_.mesh();
109 
110  label iter = 0;
111 
112  if (meshRefiner_.features().size() && maxIter > 0)
113  {
114  for (; iter < maxIter; iter++)
115  {
116  Info<< nl
117  << "Feature refinement iteration " << iter << nl
118  << "------------------------------" << nl
119  << endl;
120 
121  labelList candidateCells
122  (
123  meshRefiner_.refineCandidates
124  (
125  refineParams.locationsInMesh(),
126  refineParams.curvature(),
127  refineParams.planarAngle(),
128 
129  true, // featureRefinement
130  false, // featureDistanceRefinement
131  false, // internalRefinement
132  false, // surfaceRefinement
133  false, // curvatureRefinement
134  false, // smallFeatureRefinement
135  false, // gapRefinement
136  false, // bigGapRefinement
137  false, // spreadGapSize
138  refineParams.maxGlobalCells(),
139  refineParams.maxLocalCells()
140  )
141  );
142  labelList cellsToRefine
143  (
144  meshRefiner_.meshCutter().consistentRefinement
145  (
146  candidateCells,
147  true
148  )
149  );
150  Info<< "Determined cells to refine in = "
151  << mesh.time().cpuTimeIncrement() << " s" << endl;
152 
153 
154 
155  const label nCellsToRefine =
156  returnReduce(cellsToRefine.size(), sumOp<label>());
157 
158  Info<< "Selected for feature refinement : " << nCellsToRefine
159  << " cells (out of " << mesh.globalData().nTotalCells()
160  << ')' << endl;
161 
162  if (nCellsToRefine <= minRefine)
163  {
164  Info<< "Stopping refining since too few cells selected."
165  << nl << endl;
166  break;
167  }
168 
169 
170  if (debug > 0)
171  {
172  const_cast<Time&>(mesh.time())++;
173  }
174 
175 
176  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
177  {
178  meshRefiner_.balanceAndRefine
179  (
180  "feature refinement iteration " + name(iter),
181  decomposer_,
182  distributor_,
183  cellsToRefine,
184  refineParams.maxLoadUnbalance(),
185  refineParams.maxCellUnbalance()
186  );
187  }
188  else
189  {
190  meshRefiner_.refineAndBalance
191  (
192  "feature refinement iteration " + name(iter),
193  decomposer_,
194  distributor_,
195  cellsToRefine,
196  refineParams.maxLoadUnbalance(),
197  refineParams.maxCellUnbalance()
198  );
199  }
200  }
201  }
202  return iter;
203 }
204 
205 
206 Foam::label Foam::snappyRefineDriver::smallFeatureRefine
207 (
208  const refinementParameters& refineParams,
209  const label maxIter
210 )
211 {
212  if (dryRun_)
213  {
214  return 0;
215  }
216 
217  if (refineParams.minRefineCells() == -1)
218  {
219  // Special setting to be able to restart shm on meshes with inconsistent
220  // cellLevel/pointLevel
221  return 0;
222  }
223 
224  addProfiling(feature, "snappyHexMesh::refine::smallFeature");
225  const fvMesh& mesh = meshRefiner_.mesh();
226 
227  label iter = 0;
228 
229  // See if any surface has an extendedGapLevel
230  const labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
231  const labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
232  const labelList curvMaxLevel(meshRefiner_.surfaces().maxCurvatureLevel());
233 
234  if
235  (
236  max(surfaceMaxLevel) == 0
237  && max(shellMaxLevel) == 0
238  && max(curvMaxLevel) == 0
239  )
240  {
241  return iter;
242  }
243 
244  for (; iter < maxIter; iter++)
245  {
246  Info<< nl
247  << "Small surface feature refinement iteration " << iter << nl
248  << "--------------------------------------------" << nl
249  << endl;
250 
251 
252  // Determine cells to refine
253  // ~~~~~~~~~~~~~~~~~~~~~~~~~
254 
255  labelList candidateCells
256  (
257  meshRefiner_.refineCandidates
258  (
259  refineParams.locationsInMesh(),
260  refineParams.curvature(),
261  refineParams.planarAngle(),
262 
263  false, // featureRefinement
264  false, // featureDistanceRefinement
265  false, // internalRefinement
266  false, // surfaceRefinement
267  false, // curvatureRefinement
268  true, // smallFeatureRefinement
269  false, // gapRefinement
270  false, // bigGapRefinement
271  false, // spreadGapSize
272  refineParams.maxGlobalCells(),
273  refineParams.maxLocalCells()
274  )
275  );
276 
277  labelList cellsToRefine
278  (
279  meshRefiner_.meshCutter().consistentRefinement
280  (
281  candidateCells,
282  true
283  )
284  );
285  Info<< "Determined cells to refine in = "
286  << mesh.time().cpuTimeIncrement() << " s" << endl;
287 
288 
289  const label nCellsToRefine =
290  returnReduce(cellsToRefine.size(), sumOp<label>());
291 
292  Info<< "Selected for refinement : " << nCellsToRefine
293  << " cells (out of " << mesh.globalData().nTotalCells()
294  << ')' << endl;
295 
296  // Stop when no cells to refine or have done minimum necessary
297  // iterations and not enough cells to refine.
298  if (nCellsToRefine == 0)
299  {
300  Info<< "Stopping refining since too few cells selected."
301  << nl << endl;
302  break;
303  }
304 
305 
306  if (debug)
307  {
308  const_cast<Time&>(mesh.time())++;
309  }
310 
311 
312  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
313  {
314  meshRefiner_.balanceAndRefine
315  (
316  "small feature refinement iteration " + name(iter),
317  decomposer_,
318  distributor_,
319  cellsToRefine,
320  refineParams.maxLoadUnbalance(),
321  refineParams.maxCellUnbalance()
322  );
323  }
324  else
325  {
326  meshRefiner_.refineAndBalance
327  (
328  "small feature refinement iteration " + name(iter),
329  decomposer_,
330  distributor_,
331  cellsToRefine,
332  refineParams.maxLoadUnbalance(),
333  refineParams.maxCellUnbalance()
334  );
335  }
336  }
337  return iter;
338 }
339 
340 
341 Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
342 (
343  const refinementParameters& refineParams,
344  const label maxIter,
345  const label leakBlockageIter
346 )
347 {
348  if (dryRun_)
349  {
350  return 0;
351  }
352 
353  if (refineParams.minRefineCells() == -1)
354  {
355  // Special setting to be able to restart shm on meshes with inconsistent
356  // cellLevel/pointLevel
357  return 0;
358  }
359 
360  addProfiling(surface, "snappyHexMesh::refine::surface");
361  const fvMesh& mesh = meshRefiner_.mesh();
362  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
363 
364  // Determine the maximum refinement level over all surfaces. This
365  // determines the minimum number of surface refinement iterations.
366  label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
367 
368  label iter;
369  for (iter = 0; iter < maxIter; iter++)
370  {
371  Info<< nl
372  << "Surface refinement iteration " << iter << nl
373  << "------------------------------" << nl
374  << endl;
375 
376 
377  // Do optional leak closing (by removing cells)
378  if (iter >= leakBlockageIter)
379  {
380  // Block off intersections with unzoned surfaces with specified
381  // leakLevel < iter
382  const labelList unnamedSurfaces
383  (
384  surfaceZonesInfo::getUnnamedSurfaces
385  (
386  surfaces.surfZones()
387  )
388  );
389 
390  DynamicList<label> selectedSurfaces(unnamedSurfaces.size());
391  for (const label surfi : unnamedSurfaces)
392  {
393  const label regioni = surfaces.globalRegion(surfi, 0);
394 
395  // Take shortcut: assume all cells on surface are refined to
396  // its refinement level at iteration iter. So just use the
397  // iteration to see if the surface is a candidate.
398  if (iter > surfaces.leakLevel()[regioni])
399  {
400  selectedSurfaces.append(surfi);
401  }
402  }
403 
404  if
405  (
406  selectedSurfaces.size()
407  && refineParams.locationsOutsideMesh().size()
408  )
409  {
410  meshRefiner_.blockLeakFaces
411  (
412  globalToMasterPatch_,
413  globalToSlavePatch_,
414  refineParams.locationsInMesh(),
415  refineParams.zonesInMesh(),
416  refineParams.locationsOutsideMesh(),
417  selectedSurfaces
418  );
419  }
420  }
421 
422 
423  // Determine cells to refine
424  // ~~~~~~~~~~~~~~~~~~~~~~~~~
425  // Only look at surface intersections (minLevel and surface curvature),
426  // do not do internal refinement (refinementShells)
427 
428  labelList candidateCells
429  (
430  meshRefiner_.refineCandidates
431  (
432  refineParams.locationsInMesh(),
433  refineParams.curvature(),
434  refineParams.planarAngle(),
435 
436  false, // featureRefinement
437  false, // featureDistanceRefinement
438  false, // internalRefinement
439  true, // surfaceRefinement
440  true, // curvatureRefinement
441  false, // smallFeatureRefinement
442  false, // gapRefinement
443  false, // bigGapRefinement
444  false, // spreadGapSize
445  refineParams.maxGlobalCells(),
446  refineParams.maxLocalCells()
447  )
448  );
449  labelList cellsToRefine
450  (
451  meshRefiner_.meshCutter().consistentRefinement
452  (
453  candidateCells,
454  true
455  )
456  );
457  Info<< "Determined cells to refine in = "
458  << mesh.time().cpuTimeIncrement() << " s" << endl;
459 
460 
461  const label nCellsToRefine =
462  returnReduce(cellsToRefine.size(), sumOp<label>());
463 
464  Info<< "Selected for refinement : " << nCellsToRefine
465  << " cells (out of " << mesh.globalData().nTotalCells()
466  << ')' << endl;
467 
468  // Stop when no cells to refine or have done minimum necessary
469  // iterations and not enough cells to refine.
470  if
471  (
472  nCellsToRefine == 0
473  || (
474  iter >= overallMaxLevel
475  && nCellsToRefine <= refineParams.minRefineCells()
476  )
477  )
478  {
479  Info<< "Stopping refining since too few cells selected."
480  << nl << endl;
481  break;
482  }
483 
484 
485  if (debug)
486  {
487  const_cast<Time&>(mesh.time())++;
488  }
489 
490 
491  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
492  {
493  meshRefiner_.balanceAndRefine
494  (
495  "surface refinement iteration " + name(iter),
496  decomposer_,
497  distributor_,
498  cellsToRefine,
499  refineParams.maxLoadUnbalance(),
500  refineParams.maxCellUnbalance()
501  );
502  }
503  else
504  {
505  meshRefiner_.refineAndBalance
506  (
507  "surface refinement iteration " + name(iter),
508  decomposer_,
509  distributor_,
510  cellsToRefine,
511  refineParams.maxLoadUnbalance(),
512  refineParams.maxCellUnbalance()
513  );
514  }
515  }
516  return iter;
517 }
518 
519 
520 Foam::label Foam::snappyRefineDriver::gapOnlyRefine
521 (
522  const refinementParameters& refineParams,
523  const label maxIter
524 )
525 {
526  if (dryRun_)
527  {
528  return 0;
529  }
530 
531  if (refineParams.minRefineCells() == -1)
532  {
533  // Special setting to be able to restart shm on meshes with inconsistent
534  // cellLevel/pointLevel
535  return 0;
536  }
537 
538  const fvMesh& mesh = meshRefiner_.mesh();
539 
540  // Determine the maximum refinement level over all surfaces. This
541  // determines the minimum number of surface refinement iterations.
542 
543  label maxIncrement = 0;
544  const labelList& maxLevel = meshRefiner_.surfaces().maxLevel();
545  const labelList& gapLevel = meshRefiner_.surfaces().gapLevel();
546 
547  forAll(maxLevel, i)
548  {
549  maxIncrement = max(maxIncrement, gapLevel[i]-maxLevel[i]);
550  }
551 
552  label iter = 0;
553 
554  if (maxIncrement == 0)
555  {
556  return iter;
557  }
558 
559  for (iter = 0; iter < maxIter; iter++)
560  {
561  Info<< nl
562  << "Gap refinement iteration " << iter << nl
563  << "--------------------------" << nl
564  << endl;
565 
566 
567  // Determine cells to refine
568  // ~~~~~~~~~~~~~~~~~~~~~~~~~
569  // Only look at surface intersections (minLevel and surface curvature),
570  // do not do internal refinement (refinementShells)
571 
572  labelList candidateCells
573  (
574  meshRefiner_.refineCandidates
575  (
576  refineParams.locationsInMesh(),
577  refineParams.curvature(),
578  refineParams.planarAngle(),
579 
580  false, // featureRefinement
581  false, // featureDistanceRefinement
582  false, // internalRefinement
583  false, // surfaceRefinement
584  false, // curvatureRefinement
585  false, // smallFeatureRefinement
586  true, // gapRefinement
587  false, // bigGapRefinement
588  false, // spreadGapSize
589  refineParams.maxGlobalCells(),
590  refineParams.maxLocalCells()
591  )
592  );
593 
594  if (debug&meshRefinement::MESH)
595  {
596  Pout<< "Writing current mesh to time "
597  << meshRefiner_.timeName() << endl;
598  meshRefiner_.write
599  (
600  meshRefinement::debugType(debug),
601  meshRefinement::writeType
602  (
603  meshRefinement::writeLevel()
604  | meshRefinement::WRITEMESH
605  ),
606  mesh.time().path()/meshRefiner_.timeName()
607  );
608  Pout<< "Dumped mesh in = "
609  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
610 
611 
612  Pout<< "Dumping " << candidateCells.size()
613  << " cells to cellSet candidateCellsFromGap." << endl;
614  cellSet c(mesh, "candidateCellsFromGap", candidateCells);
615  c.instance() = meshRefiner_.timeName();
616  c.write();
617  }
618 
619  // Grow by one layer to make sure we're covering the gap
620  {
621  boolList isCandidateCell(mesh.nCells(), false);
622  forAll(candidateCells, i)
623  {
624  isCandidateCell[candidateCells[i]] = true;
625  }
626 
627  for (label i=0; i<1; i++)
628  {
629  boolList newIsCandidateCell(isCandidateCell);
630 
631  // Internal faces
632  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
633  {
634  label own = mesh.faceOwner()[facei];
635  label nei = mesh.faceNeighbour()[facei];
636 
637  if (isCandidateCell[own] != isCandidateCell[nei])
638  {
639  newIsCandidateCell[own] = true;
640  newIsCandidateCell[nei] = true;
641  }
642  }
643 
644  // Get coupled boundary condition values
645  boolList neiIsCandidateCell;
646  syncTools::swapBoundaryCellList
647  (
648  mesh,
649  isCandidateCell,
650  neiIsCandidateCell
651  );
652 
653  // Boundary faces
654  for
655  (
656  label facei = mesh.nInternalFaces();
657  facei < mesh.nFaces();
658  facei++
659  )
660  {
661  label own = mesh.faceOwner()[facei];
662  label bFacei = facei-mesh.nInternalFaces();
663 
664  if (isCandidateCell[own] != neiIsCandidateCell[bFacei])
665  {
666  newIsCandidateCell[own] = true;
667  }
668  }
669 
670  isCandidateCell.transfer(newIsCandidateCell);
671  }
672 
673  label n = 0;
674  forAll(isCandidateCell, celli)
675  {
676  if (isCandidateCell[celli])
677  {
678  n++;
679  }
680  }
681  candidateCells.setSize(n);
682  n = 0;
683  forAll(isCandidateCell, celli)
684  {
685  if (isCandidateCell[celli])
686  {
687  candidateCells[n++] = celli;
688  }
689  }
690  }
691 
692 
693  if (debug&meshRefinement::MESH)
694  {
695  Pout<< "Dumping " << candidateCells.size()
696  << " cells to cellSet candidateCellsFromGapPlusBuffer." << endl;
697  cellSet c(mesh, "candidateCellsFromGapPlusBuffer", candidateCells);
698  c.instance() = meshRefiner_.timeName();
699  c.write();
700  }
701 
702 
703  labelList cellsToRefine
704  (
705  meshRefiner_.meshCutter().consistentRefinement
706  (
707  candidateCells,
708  true
709  )
710  );
711  Info<< "Determined cells to refine in = "
712  << mesh.time().cpuTimeIncrement() << " s" << endl;
713 
714 
715  const label nCellsToRefine =
716  returnReduce(cellsToRefine.size(), sumOp<label>());
717 
718  Info<< "Selected for refinement : " << nCellsToRefine
719  << " cells (out of " << mesh.globalData().nTotalCells()
720  << ')' << endl;
721 
722  // Stop when no cells to refine or have done minimum necessary
723  // iterations and not enough cells to refine.
724  if
725  (
726  nCellsToRefine == 0
727  || (
728  iter >= maxIncrement
729  && nCellsToRefine <= refineParams.minRefineCells()
730  )
731  )
732  {
733  Info<< "Stopping refining since too few cells selected."
734  << nl << endl;
735  break;
736  }
737 
738 
739  if (debug)
740  {
741  const_cast<Time&>(mesh.time())++;
742  }
743 
744 
745  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
746  {
747  meshRefiner_.balanceAndRefine
748  (
749  "gap refinement iteration " + name(iter),
750  decomposer_,
751  distributor_,
752  cellsToRefine,
753  refineParams.maxLoadUnbalance(),
754  refineParams.maxCellUnbalance()
755  );
756  }
757  else
758  {
759  meshRefiner_.refineAndBalance
760  (
761  "gap refinement iteration " + name(iter),
762  decomposer_,
763  distributor_,
764  cellsToRefine,
765  refineParams.maxLoadUnbalance(),
766  refineParams.maxCellUnbalance()
767  );
768  }
769  }
770  return iter;
771 }
772 
773 
774 Foam::label Foam::snappyRefineDriver::surfaceProximityBlock
775 (
776  const refinementParameters& refineParams,
777  const label maxIter
778 )
779 {
780  if (refineParams.minRefineCells() == -1)
781  {
782  // Special setting to be able to restart shm on meshes with inconsistent
783  // cellLevel/pointLevel
784  return 0;
785  }
786 
787  fvMesh& mesh = meshRefiner_.mesh();
788 
789  if (min(meshRefiner_.surfaces().blockLevel()) == labelMax)
790  {
791  return 0;
792  }
793 
794  label iter = 0;
795 
796  for (iter = 0; iter < maxIter; iter++)
797  {
798  Info<< nl
799  << "Gap blocking iteration " << iter << nl
800  << "------------------------" << nl
801  << endl;
802 
803 
804  // Determine cells to block
805  // ~~~~~~~~~~~~~~~~~~~~~~~~
806 
807  meshRefiner_.removeGapCells
808  (
809  refineParams.planarAngle(),
810  meshRefiner_.surfaces().blockLevel(),
811  globalToMasterPatch_,
812  refineParams.nFilterIter()
813  );
814 
815  if (debug)
816  {
817  const_cast<Time&>(mesh.time())++;
818  Pout<< "Writing gap blocking iteration "
819  << iter << " mesh to time " << meshRefiner_.timeName()
820  << endl;
821  meshRefiner_.write
822  (
823  meshRefinement::debugType(debug),
824  meshRefinement::writeType
825  (
826  meshRefinement::writeLevel()
827  | meshRefinement::WRITEMESH
828  ),
829  mesh.time().path()/meshRefiner_.timeName()
830  );
831  }
832  }
833  return iter;
834 }
835 
836 
837 Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
838 (
839  const refinementParameters& refineParams,
840  const bool spreadGapSize,
841  const label maxIter
842 )
843 {
844  if (refineParams.minRefineCells() == -1)
845  {
846  // Special setting to be able to restart shm on meshes with inconsistent
847  // cellLevel/pointLevel
848  return 0;
849  }
850 
851  if (dryRun_)
852  {
853  return 0;
854  }
855 
856  const fvMesh& mesh = meshRefiner_.mesh();
857 
858  label iter = 0;
859 
860  // See if any surface has an extendedGapLevel
861  labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
862  labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
863 
864  label overallMaxLevel(max(max(surfaceMaxLevel), max(shellMaxLevel)));
865 
866  if (overallMaxLevel == 0)
867  {
868  return iter;
869  }
870 
871 
872  for (; iter < maxIter; iter++)
873  {
874  Info<< nl
875  << "Big gap refinement iteration " << iter << nl
876  << "------------------------------" << nl
877  << endl;
878 
879 
880  // Determine cells to refine
881  // ~~~~~~~~~~~~~~~~~~~~~~~~~
882 
883  labelList candidateCells
884  (
885  meshRefiner_.refineCandidates
886  (
887  refineParams.locationsInMesh(),
888  refineParams.curvature(),
889  refineParams.planarAngle(),
890 
891  false, // featureRefinement
892  false, // featureDistanceRefinement
893  false, // internalRefinement
894  false, // surfaceRefinement
895  false, // curvatureRefinement
896  false, // smallFeatureRefinement
897  false, // gapRefinement
898  true, // bigGapRefinement
899  spreadGapSize, // spreadGapSize
900  refineParams.maxGlobalCells(),
901  refineParams.maxLocalCells()
902  )
903  );
904 
905 
906  if (debug&meshRefinement::MESH)
907  {
908  Pout<< "Writing current mesh to time "
909  << meshRefiner_.timeName() << endl;
910  meshRefiner_.write
911  (
912  meshRefinement::debugType(debug),
913  meshRefinement::writeType
914  (
915  meshRefinement::writeLevel()
916  | meshRefinement::WRITEMESH
917  ),
918  mesh.time().path()/meshRefiner_.timeName()
919  );
920  Pout<< "Dumped mesh in = "
921  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
922 
923  Pout<< "Dumping " << candidateCells.size()
924  << " cells to cellSet candidateCellsFromBigGap." << endl;
925  cellSet c(mesh, "candidateCellsFromBigGap", candidateCells);
926  c.instance() = meshRefiner_.timeName();
927  c.write();
928  }
929 
930  labelList cellsToRefine
931  (
932  meshRefiner_.meshCutter().consistentRefinement
933  (
934  candidateCells,
935  true
936  )
937  );
938  Info<< "Determined cells to refine in = "
939  << mesh.time().cpuTimeIncrement() << " s" << endl;
940 
941 
942  const label nCellsToRefine =
943  returnReduce(cellsToRefine.size(), sumOp<label>());
944 
945  Info<< "Selected for refinement : " << nCellsToRefine
946  << " cells (out of " << mesh.globalData().nTotalCells()
947  << ')' << endl;
948 
949  // Stop when no cells to refine or have done minimum necessary
950  // iterations and not enough cells to refine.
951  if
952  (
953  nCellsToRefine == 0
954  || (
955  iter >= overallMaxLevel
956  && nCellsToRefine <= refineParams.minRefineCells()
957  )
958  )
959  {
960  Info<< "Stopping refining since too few cells selected."
961  << nl << endl;
962  break;
963  }
964 
965 
966  if (debug)
967  {
968  const_cast<Time&>(mesh.time())++;
969  }
970 
971 
972  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
973  {
974  meshRefiner_.balanceAndRefine
975  (
976  "big gap refinement iteration " + name(iter),
977  decomposer_,
978  distributor_,
979  cellsToRefine,
980  refineParams.maxLoadUnbalance(),
981  refineParams.maxCellUnbalance()
982  );
983  }
984  else
985  {
986  meshRefiner_.refineAndBalance
987  (
988  "big gap refinement iteration " + name(iter),
989  decomposer_,
990  distributor_,
991  cellsToRefine,
992  refineParams.maxLoadUnbalance(),
993  refineParams.maxCellUnbalance()
994  );
995  }
996  }
997  return iter;
998 }
999 
1000 
1001 Foam::label Foam::snappyRefineDriver::danglingCellRefine
1002 (
1003  const refinementParameters& refineParams,
1004  const label nFaces,
1005  const label maxIter
1006 )
1007 {
1008  if (refineParams.minRefineCells() == -1)
1009  {
1010  // Special setting to be able to restart shm on meshes with inconsistent
1011  // cellLevel/pointLevel
1012  return 0;
1013  }
1014 
1015  if (dryRun_)
1016  {
1017  return 0;
1018  }
1019 
1020  addProfiling(dangling, "snappyHexMesh::refine::danglingCell");
1021  const fvMesh& mesh = meshRefiner_.mesh();
1022 
1023  label iter;
1024  for (iter = 0; iter < maxIter; iter++)
1025  {
1026  Info<< nl
1027  << "Dangling coarse cells refinement iteration " << iter << nl
1028  << "--------------------------------------------" << nl
1029  << endl;
1030 
1031 
1032  // Determine cells to refine
1033  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1034 
1035  const cellList& cells = mesh.cells();
1036  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
1037 
1038  labelList candidateCells;
1039  {
1040  cellSet candidateCellSet(mesh, "candidateCells", cells.size()/1000);
1041 
1042  forAll(cells, celli)
1043  {
1044  label nIntFaces = 0;
1045  for (const label meshFacei : cells[celli])
1046  {
1047  const label patchi = pbm.patchID(meshFacei);
1048 
1049  if (patchi < 0 || pbm[patchi].coupled())
1050  {
1051  ++nIntFaces;
1052  }
1053  }
1054 
1055  if (nIntFaces == nFaces)
1056  {
1057  candidateCellSet.insert(celli);
1058  }
1059  }
1060 
1061  if (debug&meshRefinement::MESH)
1062  {
1063  Pout<< "Dumping " << candidateCellSet.size()
1064  << " cells to cellSet candidateCellSet." << endl;
1065  candidateCellSet.instance() = meshRefiner_.timeName();
1066  candidateCellSet.write();
1067  }
1068  candidateCells = candidateCellSet.toc();
1069  }
1070 
1071 
1072 
1073  labelList cellsToRefine
1074  (
1075  meshRefiner_.meshCutter().consistentRefinement
1076  (
1077  candidateCells,
1078  true
1079  )
1080  );
1081  Info<< "Determined cells to refine in = "
1082  << mesh.time().cpuTimeIncrement() << " s" << endl;
1083 
1084 
1085  const label nCellsToRefine =
1086  returnReduce(cellsToRefine.size(), sumOp<label>());
1087 
1088  Info<< "Selected for refinement : " << nCellsToRefine
1089  << " cells (out of " << mesh.globalData().nTotalCells()
1090  << ')' << endl;
1091 
1092  // Stop when no cells to refine. After a few iterations check if too
1093  // few cells
1094  if
1095  (
1096  nCellsToRefine == 0
1097  || (
1098  iter >= 1
1099  && nCellsToRefine <= refineParams.minRefineCells()
1100  )
1101  )
1102  {
1103  Info<< "Stopping refining since too few cells selected."
1104  << nl << endl;
1105 
1106  if (refineParams.balanceAtEnd())
1107  {
1108  Info<< "Final mesh balancing" << endl;
1109 
1110  meshRefiner_.balance
1111  (
1112  "",
1113  decomposer_,
1114  distributor_,
1115  cellsToRefine,
1116  0,
1117  -1
1118  );
1119  }
1120 
1121  break;
1122  }
1123 
1124 
1125  if (debug)
1126  {
1127  const_cast<Time&>(mesh.time())++;
1128  }
1129 
1130 
1131  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
1132  {
1133  meshRefiner_.balanceAndRefine
1134  (
1135  "coarse cell refinement iteration " + name(iter),
1136  decomposer_,
1137  distributor_,
1138  cellsToRefine,
1139  refineParams.maxLoadUnbalance(),
1140  refineParams.maxCellUnbalance()
1141  );
1142  }
1143  else
1144  {
1145  meshRefiner_.refineAndBalance
1146  (
1147  "coarse cell refinement iteration " + name(iter),
1148  decomposer_,
1149  distributor_,
1150  cellsToRefine,
1151  refineParams.maxLoadUnbalance(),
1152  refineParams.maxCellUnbalance()
1153  );
1154  }
1155  }
1156  return iter;
1157 }
1158 
1159 
1160 // Detect cells with opposing intersected faces of differing refinement
1161 // level and refine them.
1162 Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
1163 (
1164  const refinementParameters& refineParams,
1165  const label maxIter
1166 )
1167 {
1168  if (refineParams.minRefineCells() == -1)
1169  {
1170  // Special setting to be able to restart shm on meshes with inconsistent
1171  // cellLevel/pointLevel
1172  return 0;
1173  }
1174 
1175  if (dryRun_)
1176  {
1177  return 0;
1178  }
1179 
1180  addProfiling(interface, "snappyHexMesh::refine::transition");
1181  const fvMesh& mesh = meshRefiner_.mesh();
1182 
1183  label iter = 0;
1184 
1185  if (refineParams.interfaceRefine())
1186  {
1187  for (;iter < maxIter; iter++)
1188  {
1189  Info<< nl
1190  << "Refinement transition refinement iteration " << iter << nl
1191  << "--------------------------------------------" << nl
1192  << endl;
1193 
1194  const labelList& surfaceIndex = meshRefiner_.surfaceIndex();
1195  const hexRef8& cutter = meshRefiner_.meshCutter();
1196  const vectorField& fA = mesh.faceAreas();
1197  const labelList& faceOwner = mesh.faceOwner();
1198 
1199 
1200  // Determine cells to refine
1201  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1202 
1203  const cellList& cells = mesh.cells();
1204 
1205  labelList candidateCells;
1206  {
1207  // Pass1: pick up cells with differing face level
1208 
1209  cellSet transitionCells
1210  (
1211  mesh,
1212  "transitionCells",
1213  cells.size()/100
1214  );
1215 
1216  forAll(cells, celli)
1217  {
1218  const cell& cFaces = cells[celli];
1219  label cLevel = cutter.cellLevel()[celli];
1220 
1221  forAll(cFaces, cFacei)
1222  {
1223  label facei = cFaces[cFacei];
1224 
1225  if (surfaceIndex[facei] != -1)
1226  {
1227  label fLevel = cutter.faceLevel(facei);
1228  if (fLevel != cLevel)
1229  {
1230  transitionCells.insert(celli);
1231  }
1232  }
1233  }
1234  }
1235 
1236 
1237  cellSet candidateCellSet
1238  (
1239  mesh,
1240  "candidateCells",
1241  cells.size()/1000
1242  );
1243 
1244  // Pass2: check for oppositeness
1245 
1246  //for (const label celli : transitionCells)
1247  //{
1248  // const cell& cFaces = cells[celli];
1249  // const point& cc = cellCentres[celli];
1250  // const scalar rCVol = pow(cellVolumes[celli], -5.0/3.0);
1251  //
1252  // // Determine principal axes of cell
1253  // symmTensor R(Zero);
1254  //
1255  // forAll(cFaces, i)
1256  // {
1257  // label facei = cFaces[i];
1258  //
1259  // const point& fc = faceCentres[facei];
1260  //
1261  // // Calculate face-pyramid volume
1262  // scalar pyrVol = 1.0/3.0 * fA[facei] & (fc-cc);
1263  //
1264  // if (faceOwner[facei] != celli)
1265  // {
1266  // pyrVol = -pyrVol;
1267  // }
1268  //
1269  // // Calculate face-pyramid centre
1270  // vector pc = (3.0/4.0)*fc + (1.0/4.0)*cc;
1271  //
1272  // R += pyrVol*sqr(pc-cc)*rCVol;
1273  // }
1274  //
1275  // //- MEJ: Problem: truncation errors cause complex evs
1276  // vector lambdas(eigenValues(R));
1277  // const tensor axes(eigenVectors(R, lambdas));
1278  //
1279  //
1280  // // Check if this cell has
1281  // // - opposing sides intersected
1282  // // - which are of different refinement level
1283  // // - plus the inbetween face
1284  //
1285  // labelVector plusFaceLevel(labelVector(-1, -1, -1));
1286  // labelVector minFaceLevel(labelVector(-1, -1, -1));
1287  //
1288  // forAll(cFaces, cFacei)
1289  // {
1290  // label facei = cFaces[cFacei];
1291  //
1292  // if (surfaceIndex[facei] != -1)
1293  // {
1294  // label fLevel = cutter.faceLevel(facei);
1295  //
1296  // // Get outwards pointing normal
1297  // vector n = fA[facei]/mag(fA[facei]);
1298  // if (faceOwner[facei] != celli)
1299  // {
1300  // n = -n;
1301  // }
1302  //
1303  // // What is major direction and sign
1304  // direction cmpt = vector::X;
1305  // scalar maxComp = (n&axes.x());
1306  //
1307  // scalar yComp = (n&axes.y());
1308  // scalar zComp = (n&axes.z());
1309  //
1310  // if (mag(yComp) > mag(maxComp))
1311  // {
1312  // maxComp = yComp;
1313  // cmpt = vector::Y;
1314  // }
1315  //
1316  // if (mag(zComp) > mag(maxComp))
1317  // {
1318  // maxComp = zComp;
1319  // cmpt = vector::Z;
1320  // }
1321  //
1322  // if (maxComp > 0)
1323  // {
1324  // plusFaceLevel[cmpt] = max
1325  // (
1326  // plusFaceLevel[cmpt],
1327  // fLevel
1328  // );
1329  // }
1330  // else
1331  // {
1332  // minFaceLevel[cmpt] = max
1333  // (
1334  // minFaceLevel[cmpt],
1335  // fLevel
1336  // );
1337  // }
1338  // }
1339  // }
1340  //
1341  // // Check if we picked up any opposite differing level
1342  // for (direction dir = 0; dir < vector::nComponents; dir++)
1343  // {
1344  // if
1345  // (
1346  // plusFaceLevel[dir] != -1
1347  // && minFaceLevel[dir] != -1
1348  // && plusFaceLevel[dir] != minFaceLevel[dir]
1349  // )
1350  // {
1351  // candidateCellSet.insert(celli);
1352  // }
1353  // }
1354  //}
1355 
1356  const scalar oppositeCos = Foam::cos(degToRad(135.0));
1357 
1358  for (const label celli : transitionCells)
1359  {
1360  const cell& cFaces = cells[celli];
1361  label cLevel = cutter.cellLevel()[celli];
1362 
1363  // Detect opposite intersection
1364  bool foundOpposite = false;
1365 
1366  forAll(cFaces, cFacei)
1367  {
1368  label facei = cFaces[cFacei];
1369 
1370  if
1371  (
1372  surfaceIndex[facei] != -1
1373  && cutter.faceLevel(facei) > cLevel
1374  )
1375  {
1376  // Get outwards pointing normal
1377  vector n = fA[facei]/mag(fA[facei]);
1378  if (faceOwner[facei] != celli)
1379  {
1380  n = -n;
1381  }
1382 
1383  // Check for any opposite intersection
1384  forAll(cFaces, cFaceI2)
1385  {
1386  label face2i = cFaces[cFaceI2];
1387 
1388  if
1389  (
1390  face2i != facei
1391  && surfaceIndex[face2i] != -1
1392  && cutter.faceLevel(face2i) > cLevel
1393  )
1394  {
1395  // Get outwards pointing normal
1396  vector n2 = fA[face2i]/mag(fA[face2i]);
1397  if (faceOwner[face2i] != celli)
1398  {
1399  n2 = -n2;
1400  }
1401 
1402 
1403  if ((n&n2) < oppositeCos)
1404  {
1405  foundOpposite = true;
1406  break;
1407  }
1408  }
1409  }
1410 
1411  if (foundOpposite)
1412  {
1413  break;
1414  }
1415  }
1416  }
1417 
1418 
1419  if (foundOpposite)
1420  {
1421  candidateCellSet.insert(celli);
1422  }
1423  }
1424 
1425  if (debug&meshRefinement::MESH)
1426  {
1427  Pout<< "Dumping " << candidateCellSet.size()
1428  << " cells to cellSet candidateCellSet." << endl;
1429  candidateCellSet.instance() = meshRefiner_.timeName();
1430  candidateCellSet.write();
1431  }
1432  candidateCells = candidateCellSet.toc();
1433  }
1434 
1435 
1436 
1437  labelList cellsToRefine
1438  (
1439  meshRefiner_.meshCutter().consistentRefinement
1440  (
1441  candidateCells,
1442  true
1443  )
1444  );
1445  Info<< "Determined cells to refine in = "
1446  << mesh.time().cpuTimeIncrement() << " s" << endl;
1447 
1448 
1449  const label nCellsToRefine =
1450  returnReduce(cellsToRefine.size(), sumOp<label>());
1451 
1452  Info<< "Selected for refinement : " << nCellsToRefine
1453  << " cells (out of " << mesh.globalData().nTotalCells()
1454  << ')' << endl;
1455 
1456  // Stop when no cells to refine. After a few iterations check if too
1457  // few cells
1458  if
1459  (
1460  nCellsToRefine == 0
1461  || (
1462  iter >= 1
1463  && nCellsToRefine <= refineParams.minRefineCells()
1464  )
1465  )
1466  {
1467  Info<< "Stopping refining since too few cells selected."
1468  << nl << endl;
1469  break;
1470  }
1471 
1472 
1473  if (debug)
1474  {
1475  const_cast<Time&>(mesh.time())++;
1476  }
1477 
1478 
1479  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
1480  {
1481  meshRefiner_.balanceAndRefine
1482  (
1483  "interface cell refinement iteration " + name(iter),
1484  decomposer_,
1485  distributor_,
1486  cellsToRefine,
1487  refineParams.maxLoadUnbalance(),
1488  refineParams.maxCellUnbalance()
1489  );
1490  }
1491  else
1492  {
1493  meshRefiner_.refineAndBalance
1494  (
1495  "interface cell refinement iteration " + name(iter),
1496  decomposer_,
1497  distributor_,
1498  cellsToRefine,
1499  refineParams.maxLoadUnbalance(),
1500  refineParams.maxCellUnbalance()
1501  );
1502  }
1503  }
1504  }
1505  return iter;
1506 }
1507 
1508 bool Foam::snappyRefineDriver::usesHigherLevel
1509 (
1510  const labelUList& boundaryPointLevel,
1511  const labelUList& f,
1512  const label cLevel
1513 ) const
1514 {
1515  for (const label pointi : f)
1516  {
1517  if (boundaryPointLevel[pointi] > cLevel)
1518  {
1519  return true;
1520  }
1521  }
1522  return false;
1523 }
1524 
1525 
1526 Foam::label Foam::snappyRefineDriver::boundaryRefinementInterfaceRefine
1527 (
1528  const refinementParameters& refineParams,
1529  const label maxIter
1530 )
1531 {
1532  if (refineParams.minRefineCells() == -1)
1533  {
1534  // Special setting to be able to restart shm on meshes with inconsistent
1535  // cellLevel/pointLevel
1536  return 0;
1537  }
1538 
1539  if (dryRun_)
1540  {
1541  return 0;
1542  }
1543 
1544  addProfiling(interface, "snappyHexMesh::refine::transition");
1545  const fvMesh& mesh = meshRefiner_.mesh();
1546 
1547  label iter = 0;
1548 
1549  if (refineParams.interfaceRefine())
1550  {
1551  for (;iter < maxIter; iter++)
1552  {
1553  Info<< nl
1554  << "Boundary refinement iteration " << iter << nl
1555  << "-------------------------------" << nl
1556  << endl;
1557 
1558  const labelList& surfaceIndex = meshRefiner_.surfaceIndex();
1559  const hexRef8& cutter = meshRefiner_.meshCutter();
1560  const labelList& cellLevel = cutter.cellLevel();
1561  const faceList& faces = mesh.faces();
1562  const cellList& cells = mesh.cells();
1563 
1564 
1565  // Determine cells to refine
1566  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1567 
1568  // Point/face on boundary
1569  bitSet isBoundaryPoint(mesh.nPoints());
1570  bitSet isBoundaryFace(mesh.nFaces());
1571  {
1572  forAll(surfaceIndex, facei)
1573  {
1574  if (surfaceIndex[facei] != -1)
1575  {
1576  isBoundaryFace.set(facei);
1577  isBoundaryPoint.set(faces[facei]);
1578  }
1579  }
1580  const labelList meshPatchIDs(meshRefiner_.meshedPatches());
1581  for (const label patchi : meshPatchIDs)
1582  {
1583  const polyPatch& pp = mesh.boundaryMesh()[patchi];
1584  forAll(pp, i)
1585  {
1586  isBoundaryFace.set(pp.start()+i);
1587  isBoundaryPoint.set(pp[i]);
1588  }
1589  }
1590 
1591  syncTools::syncPointList
1592  (
1593  mesh,
1594  isBoundaryPoint,
1595  orEqOp<unsigned int>(),
1596  0
1597  );
1598  }
1599 
1600  // Mark max boundary face level onto boundary points. All points
1601  // not on a boundary face stay 0.
1602  labelList boundaryPointLevel(mesh.nPoints(), 0);
1603  {
1604  forAll(cells, celli)
1605  {
1606  const cell& cFaces = cells[celli];
1607  const label cLevel = cellLevel[celli];
1608 
1609  for (const label facei : cFaces)
1610  {
1611  if (isBoundaryFace(facei))
1612  {
1613  const face& f = faces[facei];
1614  for (const label pointi : f)
1615  {
1616  boundaryPointLevel[pointi] =
1617  max
1618  (
1619  boundaryPointLevel[pointi],
1620  cLevel
1621  );
1622  }
1623  }
1624  }
1625  }
1626 
1627  syncTools::syncPointList
1628  (
1629  mesh,
1630  boundaryPointLevel,
1631  maxEqOp<label>(),
1632  label(0)
1633  );
1634  }
1635 
1636 
1637  // Detect cells with a point but not face on the boundary
1638  labelList candidateCells;
1639  {
1640  const cellList& cells = mesh.cells();
1641 
1642  cellSet candidateCellSet
1643  (
1644  mesh,
1645  "candidateCells",
1646  cells.size()/100
1647  );
1648 
1649  forAll(cells, celli)
1650  {
1651  const cell& cFaces = cells[celli];
1652  const label cLevel = cellLevel[celli];
1653 
1654  bool isBoundaryCell = false;
1655  for (const label facei : cFaces)
1656  {
1657  if (isBoundaryFace(facei))
1658  {
1659  isBoundaryCell = true;
1660  break;
1661  }
1662  }
1663 
1664  if (!isBoundaryCell)
1665  {
1666  for (const label facei : cFaces)
1667  {
1668  const face& f = mesh.faces()[facei];
1669  if (usesHigherLevel(boundaryPointLevel, f, cLevel))
1670  {
1671  candidateCellSet.insert(celli);
1672  break;
1673  }
1674  }
1675  }
1676  }
1677  if (debug&meshRefinement::MESH)
1678  {
1679  Pout<< "Dumping " << candidateCellSet.size()
1680  << " cells to cellSet candidateCellSet." << endl;
1681  candidateCellSet.instance() = meshRefiner_.timeName();
1682  candidateCellSet.write();
1683  }
1684  candidateCells = candidateCellSet.toc();
1685  }
1686 
1687  labelList cellsToRefine
1688  (
1689  meshRefiner_.meshCutter().consistentRefinement
1690  (
1691  candidateCells,
1692  true
1693  )
1694  );
1695  Info<< "Determined cells to refine in = "
1696  << mesh.time().cpuTimeIncrement() << " s" << endl;
1697 
1698 
1699  const label nCellsToRefine =
1700  returnReduce(cellsToRefine.size(), sumOp<label>());
1701 
1702  Info<< "Selected for refinement : " << nCellsToRefine
1703  << " cells (out of " << mesh.globalData().nTotalCells()
1704  << ')' << endl;
1705 
1706  // Stop when no cells to refine. After a few iterations check if too
1707  // few cells
1708  if
1709  (
1710  nCellsToRefine == 0
1711  // || (
1712  // iter >= 1
1713  // && nCellsToRefine <= refineParams.minRefineCells()
1714  // )
1715  )
1716  {
1717  Info<< "Stopping refining since too few cells selected."
1718  << nl << endl;
1719  break;
1720  }
1721 
1722 
1723  if (debug)
1724  {
1725  const_cast<Time&>(mesh.time())++;
1726  }
1727 
1728 
1729  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
1730  {
1731  meshRefiner_.balanceAndRefine
1732  (
1733  "boundary cell refinement iteration " + name(iter),
1734  decomposer_,
1735  distributor_,
1736  cellsToRefine,
1737  refineParams.maxLoadUnbalance(),
1738  refineParams.maxCellUnbalance()
1739  );
1740  }
1741  else
1742  {
1743  meshRefiner_.refineAndBalance
1744  (
1745  "boundary cell refinement iteration " + name(iter),
1746  decomposer_,
1747  distributor_,
1748  cellsToRefine,
1749  refineParams.maxLoadUnbalance(),
1750  refineParams.maxCellUnbalance()
1751  );
1752  }
1753  }
1754  }
1755  return iter;
1756 }
1757 
1758 
1759 void Foam::snappyRefineDriver::removeInsideCells
1760 (
1761  const refinementParameters& refineParams,
1762  const label nBufferLayers
1763 )
1764 {
1765  // Skip if no limitRegion and zero bufferLayers
1766  if (meshRefiner_.limitShells().shells().size() == 0 && nBufferLayers == 0)
1767  {
1768  return;
1769  }
1770 
1771  if (dryRun_)
1772  {
1773  return;
1774  }
1775 
1776  Info<< nl
1777  << "Removing mesh beyond surface intersections" << nl
1778  << "------------------------------------------" << nl
1779  << endl;
1780 
1781  const fvMesh& mesh = meshRefiner_.mesh();
1782 
1783  if (debug)
1784  {
1785  const_cast<Time&>(mesh.time())++;
1786  }
1787 
1788  // Remove any cells inside limitShells with level -1
1789  if (meshRefiner_.limitShells().shells().size())
1790  {
1791  meshRefiner_.removeLimitShells
1792  (
1793  nBufferLayers,
1794  1,
1795  globalToMasterPatch_,
1796  globalToSlavePatch_,
1797  refineParams.locationsInMesh(),
1798  refineParams.zonesInMesh(),
1799  refineParams.locationsOutsideMesh()
1800  );
1801  }
1802 
1803 
1804  // Fix any additional (e.g. locationsOutsideMesh). Note: probably not
1805  // necessary.
1806  meshRefiner_.splitMesh
1807  (
1808  nBufferLayers, // nBufferLayers
1809  refineParams.nErodeCellZone(),
1810  globalToMasterPatch_,
1811  globalToSlavePatch_,
1812  refineParams.locationsInMesh(),
1813  refineParams.zonesInMesh(),
1814  refineParams.locationsOutsideMesh(),
1815  !refineParams.useLeakClosure(),
1816  setFormatter_,
1817  surfFormatter_
1818  );
1819 
1820  if (debug&meshRefinement::MESH)
1821  {
1822  const_cast<Time&>(mesh.time())++;
1823 
1824  Pout<< "Writing subsetted mesh to time "
1825  << meshRefiner_.timeName() << endl;
1826  meshRefiner_.write
1827  (
1828  meshRefinement::debugType(debug),
1829  meshRefinement::writeType
1830  (
1831  meshRefinement::writeLevel()
1832  | meshRefinement::WRITEMESH
1833  ),
1834  mesh.time().path()/meshRefiner_.timeName()
1835  );
1836  Pout<< "Dumped mesh in = "
1837  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1838  }
1839 }
1840 
1841 
1842 Foam::label Foam::snappyRefineDriver::shellRefine
1843 (
1844  const refinementParameters& refineParams,
1845  const label maxIter
1846 )
1847 {
1848  if (dryRun_)
1849  {
1850  return 0;
1851  }
1852 
1853  if (refineParams.minRefineCells() == -1)
1854  {
1855  // Special setting to be able to restart shm on meshes with inconsistent
1856  // cellLevel/pointLevel
1857  return 0;
1858  }
1859 
1860  addProfiling(shell, "snappyHexMesh::refine::shell");
1861  const fvMesh& mesh = meshRefiner_.mesh();
1862 
1863  // Mark current boundary faces with 0. Have meshRefiner maintain them.
1864  meshRefiner_.userFaceData().setSize(1);
1865 
1866  // mark list to remove any refined faces
1867  meshRefiner_.userFaceData()[0].first() = meshRefinement::REMOVE;
1868  meshRefiner_.userFaceData()[0].second() = ListOps::createWithValue<label>
1869  (
1870  mesh.nFaces(),
1871  meshRefiner_.intersectedFaces(),
1872  0, // set value
1873  -1 // default value
1874  );
1875 
1876  // Determine the maximum refinement level over all volume refinement
1877  // regions. This determines the minimum number of shell refinement
1878  // iterations.
1879  label overallMaxShellLevel = meshRefiner_.shells().maxLevel();
1880 
1881  label iter;
1882  for (iter = 0; iter < maxIter; iter++)
1883  {
1884  Info<< nl
1885  << "Shell refinement iteration " << iter << nl
1886  << "----------------------------" << nl
1887  << endl;
1888 
1889  labelList candidateCells
1890  (
1891  meshRefiner_.refineCandidates
1892  (
1893  refineParams.locationsInMesh(),
1894  refineParams.curvature(),
1895  refineParams.planarAngle(),
1896 
1897  false, // featureRefinement
1898  true, // featureDistanceRefinement
1899  true, // internalRefinement
1900  false, // surfaceRefinement
1901  false, // curvatureRefinement
1902  false, // smallFeatureRefinement
1903  false, // gapRefinement
1904  false, // bigGapRefinement
1905  false, // spreadGapSize
1906  refineParams.maxGlobalCells(),
1907  refineParams.maxLocalCells()
1908  )
1909  );
1910 
1911  if (debug&meshRefinement::MESH)
1912  {
1913  Pout<< "Dumping " << candidateCells.size()
1914  << " cells to cellSet candidateCellsFromShells." << endl;
1915 
1916  cellSet c(mesh, "candidateCellsFromShells", candidateCells);
1917  c.instance() = meshRefiner_.timeName();
1918  c.write();
1919  }
1920 
1921  // Problem choosing starting faces for bufferlayers (bFaces)
1922  // - we can't use the current intersected boundary faces
1923  // (intersectedFaces) since this grows indefinitely
1924  // - if we use 0 faces we don't satisfy bufferLayers from the
1925  // surface.
1926  // - possibly we want to have bFaces only the initial set of faces
1927  // and maintain the list while doing the refinement.
1928  labelList bFaces
1929  (
1930  findIndices(meshRefiner_.userFaceData()[0].second(), 0)
1931  );
1932 
1933  //Info<< "Collected boundary faces : "
1934  // << returnReduce(bFaces.size(), sumOp<label>()) << endl;
1935 
1936  labelList cellsToRefine;
1937 
1938  if (refineParams.nBufferLayers() <= 2)
1939  {
1940  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement
1941  (
1942  refineParams.nBufferLayers(),
1943  candidateCells, // cells to refine
1944  bFaces, // faces for nBufferLayers
1945  1, // point difference
1946  meshRefiner_.intersectedPoints() // points to check
1947  );
1948  }
1949  else
1950  {
1951  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement2
1952  (
1953  refineParams.nBufferLayers(),
1954  candidateCells, // cells to refine
1955  bFaces // faces for nBufferLayers
1956  );
1957  }
1958 
1959  Info<< "Determined cells to refine in = "
1960  << mesh.time().cpuTimeIncrement() << " s" << endl;
1961 
1962 
1963  const label nCellsToRefine =
1964  returnReduce(cellsToRefine.size(), sumOp<label>());
1965 
1966  Info<< "Selected for internal refinement : " << nCellsToRefine
1967  << " cells (out of " << mesh.globalData().nTotalCells()
1968  << ')' << endl;
1969 
1970  // Stop when no cells to refine or have done minimum necessary
1971  // iterations and not enough cells to refine.
1972  if
1973  (
1974  nCellsToRefine == 0
1975  || (
1976  iter >= overallMaxShellLevel
1977  && nCellsToRefine <= refineParams.minRefineCells()
1978  )
1979  )
1980  {
1981  Info<< "Stopping refining since too few cells selected."
1982  << nl << endl;
1983  break;
1984  }
1985 
1986 
1987  if (debug)
1988  {
1989  const_cast<Time&>(mesh.time())++;
1990  }
1991 
1992  if (returnReduceOr(mesh.nCells() >= refineParams.maxLocalCells()))
1993  {
1994  meshRefiner_.balanceAndRefine
1995  (
1996  "shell refinement iteration " + name(iter),
1997  decomposer_,
1998  distributor_,
1999  cellsToRefine,
2000  refineParams.maxLoadUnbalance(),
2001  refineParams.maxCellUnbalance()
2002  );
2003  }
2004  else
2005  {
2006  meshRefiner_.refineAndBalance
2007  (
2008  "shell refinement iteration " + name(iter),
2009  decomposer_,
2010  distributor_,
2011  cellsToRefine,
2012  refineParams.maxLoadUnbalance(),
2013  refineParams.maxCellUnbalance()
2014  );
2015  }
2016  }
2017  meshRefiner_.userFaceData().clear();
2018 
2019  return iter;
2020 }
2021 
2022 
2023 Foam::label Foam::snappyRefineDriver::directionalShellRefine
2024 (
2025  const refinementParameters& refineParams,
2026  const label maxIter
2027 )
2028 {
2029  if (dryRun_)
2030  {
2031  return 0;
2032  }
2033 
2034  addProfiling(shell, "snappyHexMesh::refine::directionalShell");
2035  const fvMesh& mesh = meshRefiner_.mesh();
2036  const shellSurfaces& shells = meshRefiner_.shells();
2037 
2038  labelList& cellLevel =
2039  const_cast<labelIOList&>(meshRefiner_.meshCutter().cellLevel());
2040  labelList& pointLevel =
2041  const_cast<labelIOList&>(meshRefiner_.meshCutter().pointLevel());
2042 
2043 
2044  // Determine the minimum and maximum cell levels that are candidates for
2045  // directional refinement
2046  const labelPairList dirSelect(shells.directionalSelectLevel());
2047  label overallMinLevel = labelMax;
2048  label overallMaxLevel = labelMin;
2049  forAll(dirSelect, shelli)
2050  {
2051  overallMinLevel = min(dirSelect[shelli].first(), overallMinLevel);
2052  overallMaxLevel = max(dirSelect[shelli].second(), overallMaxLevel);
2053  }
2054 
2055  if (overallMinLevel > overallMaxLevel)
2056  {
2057  return 0;
2058  }
2059 
2060  // Maintain directional refinement levels
2061  List<labelVector> dirCellLevel(cellLevel.size());
2062  forAll(cellLevel, celli)
2063  {
2064  dirCellLevel[celli] = labelVector::uniform(cellLevel[celli]);
2065  }
2066 
2067  label iter;
2068  for (iter = 0; iter < maxIter; iter++)
2069  {
2070  Info<< nl
2071  << "Directional shell refinement iteration " << iter << nl
2072  << "----------------------------------------" << nl
2073  << endl;
2074 
2075  label nAllRefine = 0;
2076 
2077  for (direction dir = 0; dir < vector::nComponents; dir++)
2078  {
2079  // Select the cells that need to be refined in certain direction:
2080  // - cell inside/outside shell
2081  // - original cellLevel (using mapping) mentioned in levelIncrement
2082  // - dirCellLevel not yet up to cellLevel+levelIncrement
2083 
2084 
2085  // Extract component of directional level
2086  labelList currentLevel(dirCellLevel.size());
2087  forAll(dirCellLevel, celli)
2088  {
2089  currentLevel[celli] = dirCellLevel[celli][dir];
2090  }
2091 
2092  labelList candidateCells
2093  (
2094  meshRefiner_.directionalRefineCandidates
2095  (
2096  refineParams.maxGlobalCells(),
2097  refineParams.maxLocalCells(),
2098  currentLevel,
2099  dir
2100  )
2101  );
2102 
2103  // Extend to keep 2:1 ratio
2104  labelList cellsToRefine
2105  (
2106  meshRefiner_.meshCutter().consistentRefinement
2107  (
2108  currentLevel,
2109  candidateCells,
2110  true
2111  )
2112  );
2113 
2114  Info<< "Determined cells to refine in = "
2115  << mesh.time().cpuTimeIncrement() << " s" << endl;
2116 
2117  const label nCellsToRefine =
2118  returnReduce(cellsToRefine.size(), sumOp<label>());
2119 
2120  Info<< "Selected for direction " << vector::componentNames[dir]
2121  << " refinement : " << nCellsToRefine
2122  << " cells (out of " << mesh.globalData().nTotalCells()
2123  << ')' << endl;
2124 
2125  nAllRefine += nCellsToRefine;
2126 
2127  // Stop when no cells to refine or have done minimum necessary
2128  // iterations and not enough cells to refine.
2129  if (nCellsToRefine > 0)
2130  {
2131  if (debug)
2132  {
2133  const_cast<Time&>(mesh.time())++;
2134  }
2135 
2136  const bitSet isRefineCell(mesh.nCells(), cellsToRefine);
2137 
2138  autoPtr<mapPolyMesh> map
2139  (
2140  meshRefiner_.directionalRefine
2141  (
2142  "directional refinement iteration " + name(iter),
2143  dir,
2144  cellsToRefine
2145  )
2146  );
2147 
2148  Info<< "Refined mesh in = "
2149  << mesh.time().cpuTimeIncrement() << " s" << endl;
2150 
2151  meshRefinement::updateList
2152  (
2153  map().cellMap(),
2154  labelVector(0, 0, 0),
2155  dirCellLevel
2156  );
2157 
2158  // Note: edges will have been split. The points might have
2159  // inherited pointLevel from either side of the edge which
2160  // might not be the same for coupled edges so sync
2161  syncTools::syncPointList
2162  (
2163  mesh,
2164  pointLevel,
2165  maxEqOp<label>(),
2166  labelMin
2167  );
2168 
2169  forAll(map().cellMap(), celli)
2170  {
2171  if (isRefineCell[map().cellMap()[celli]])
2172  {
2173  dirCellLevel[celli][dir]++;
2174  }
2175  }
2176 
2177  // Do something with the pointLevel. See discussion about the
2178  // cellLevel. Do we keep min/max ?
2179  forAll(map().pointMap(), pointi)
2180  {
2181  label oldPointi = map().pointMap()[pointi];
2182  if (map().reversePointMap()[oldPointi] != pointi)
2183  {
2184  // Is added point (splitting an edge)
2185  pointLevel[pointi]++;
2186  }
2187  }
2188  }
2189  }
2190 
2191 
2192  if (nAllRefine == 0)
2193  {
2194  Info<< "Stopping refining since no cells selected."
2195  << nl << endl;
2196  break;
2197  }
2198 
2199  meshRefiner_.printMeshInfo
2200  (
2201  debug,
2202  "After directional refinement iteration " + name(iter),
2203  true
2204  );
2205 
2206  if (debug&meshRefinement::MESH)
2207  {
2208  Pout<< "Writing directional refinement iteration "
2209  << iter << " mesh to time " << meshRefiner_.timeName() << endl;
2210  meshRefiner_.write
2211  (
2212  meshRefinement::debugType(debug),
2213  meshRefinement::writeType
2214  (
2215  meshRefinement::writeLevel()
2216  | meshRefinement::WRITEMESH
2217  ),
2218  mesh.time().path()/meshRefiner_.timeName()
2219  );
2220  }
2221  }
2222 
2223  // Adjust cellLevel from dirLevel? As max? Or the min?
2224  // For now: use max. The idea is that if there is a wall
2225  // any directional refinement is likely to be aligned with
2226  // the wall (wall layers) so any snapping/layering would probably
2227  // want to use this highest refinement level.
2228 
2229  forAll(cellLevel, celli)
2230  {
2231  cellLevel[celli] = cmptMax(dirCellLevel[celli]);
2232  }
2233 
2234  return iter;
2235 }
2236 
2237 
2238 void Foam::snappyRefineDriver::mergeAndSmoothRatio
2239 (
2240  const scalarList& allSeedPointDist,
2241  const label nSmoothExpansion,
2242  List<Tuple2<scalar, scalar>>& keyAndValue
2243 )
2244 {
2245  // Merge duplicate distance from coupled locations to get unique
2246  // distances to operate on, do on master
2247  SortableList<scalar> unmergedDist(allSeedPointDist);
2248  DynamicList<scalar> mergedDist;
2249 
2250  scalar prevDist = GREAT;
2251  forAll(unmergedDist, i)
2252  {
2253  scalar curDist = unmergedDist[i];
2254  scalar difference = mag(curDist - prevDist);
2255  if (difference > meshRefiner_.mergeDistance())
2256  //if (difference > 0.01)
2257  {
2258  mergedDist.append(curDist);
2259  prevDist = curDist;
2260  }
2261  }
2262 
2263  // Sort the unique distances
2264  SortableList<scalar> sortedDist(mergedDist);
2265  labelList indexSet = sortedDist.indices();
2266 
2267  // Get updated position starting from original (undistorted) mesh
2268  scalarList seedPointsNewLocation = sortedDist;
2269 
2270  scalar initResidual = 0.0;
2271  scalar prevIterResidual = GREAT;
2272 
2273  for (label iter = 0; iter < nSmoothExpansion; iter++)
2274  {
2275 
2276  // Position based edge averaging algorithm operated on
2277  // all seed plane locations in normalized form.
2278  //
2279  // 0 1 2 3 4 5 6 (edge numbers)
2280  // ---x---x---x---x---x---x---
2281  // 0 1 2 3 4 5 (point numbers)
2282  //
2283  // Average of edge 1-3 in terms of position
2284  // = (point3 - point0)/3
2285  // Keeping points 0-1 frozen, new position of point 2
2286  // = position2 + (average of edge 1-3 as above)
2287  for(label i = 2; i<mergedDist.size()-1; i++)
2288  {
2289  scalar oldX00 = sortedDist[i-2];
2290  scalar oldX1 = sortedDist[i+1];
2291  scalar curX0 = seedPointsNewLocation[i-1];
2292  seedPointsNewLocation[i] = curX0 + (oldX1 - oldX00)/3;
2293  }
2294 
2295  const scalarField residual(seedPointsNewLocation-sortedDist);
2296  {
2297  scalar res(sumMag(residual));
2298 
2299  if (iter == 0)
2300  {
2301  initResidual = res;
2302  }
2303  res /= initResidual;
2304 
2305  if (mag(prevIterResidual - res) < SMALL)
2306  {
2307  if (debug)
2308  {
2309  Pout<< "Converged with iteration " << iter
2310  << " initResidual: " << initResidual
2311  << " final residual : " << res << endl;
2312  }
2313  break;
2314  }
2315  else
2316  {
2317  prevIterResidual = res;
2318  }
2319  }
2320 
2321  // Update the field for next iteration, avoid moving points
2322  sortedDist = seedPointsNewLocation;
2323 
2324  }
2325 
2326  keyAndValue.setSize(mergedDist.size());
2327 
2328  forAll(mergedDist, i)
2329  {
2330  keyAndValue[i].first() = mergedDist[i];
2331  label index = indexSet[i];
2332  keyAndValue[i].second() = seedPointsNewLocation[index];
2333  }
2334 }
2335 
2336 
2337 Foam::label Foam::snappyRefineDriver::directionalSmooth
2338 (
2339  const refinementParameters& refineParams
2340 )
2341 {
2342  addProfiling(split, "snappyHexMesh::refine::smooth");
2343  Info<< nl
2344  << "Directional expansion ratio smoothing" << nl
2345  << "-------------------------------------" << nl
2346  << endl;
2347 
2348  fvMesh& baseMesh = meshRefiner_.mesh();
2349  const searchableSurfaces& geometry = meshRefiner_.surfaces().geometry();
2350  const shellSurfaces& shells = meshRefiner_.shells();
2351 
2352  label iter = 0;
2353 
2354  forAll(shells.nSmoothExpansion(), shellI)
2355  {
2356  if
2357  (
2358  shells.nSmoothExpansion()[shellI] > 0
2359  || shells.nSmoothPosition()[shellI] > 0
2360  )
2361  {
2362  label surfi = shells.shells()[shellI];
2363  const vector& userDirection = shells.smoothDirection()[shellI];
2364 
2365 
2366  // Extract inside points
2367  labelList pointLabels;
2368  {
2369  // Get inside points
2370  List<volumeType> volType;
2371  geometry[surfi].getVolumeType(baseMesh.points(), volType);
2372 
2373  label nInside = 0;
2374  forAll(volType, pointi)
2375  {
2376  if (volType[pointi] == volumeType::INSIDE)
2377  {
2378  nInside++;
2379  }
2380  }
2381  pointLabels.setSize(nInside);
2382  nInside = 0;
2383  forAll(volType, pointi)
2384  {
2385  if (volType[pointi] == volumeType::INSIDE)
2386  {
2387  pointLabels[nInside++] = pointi;
2388  }
2389  }
2390 
2391  //bitSet isInsidePoint(baseMesh.nPoints());
2392  //forAll(volType, pointi)
2393  //{
2394  // if (volType[pointi] == volumeType::INSIDE)
2395  // {
2396  // isInsidePoint.set(pointi);
2397  // }
2398  //}
2399  //pointLabels = isInsidePoint.used();
2400  }
2401 
2402  // Mark all directed edges
2403  bitSet isXEdge(baseMesh.edges().size());
2404  {
2405  const edgeList& edges = baseMesh.edges();
2406  forAll(edges, edgei)
2407  {
2408  const edge& e = edges[edgei];
2409  vector eVec(e.vec(baseMesh.points()));
2410  eVec /= mag(eVec);
2411  if (mag(eVec&userDirection) > 0.9)
2412  {
2413  isXEdge.set(edgei);
2414  }
2415  }
2416  }
2417 
2418  // Get the extreme of smoothing region and
2419  // normalize all points within
2420  const scalar totalLength =
2421  geometry[surfi].bounds().span()
2422  & userDirection;
2423  const scalar startPosition =
2424  geometry[surfi].bounds().min()
2425  & userDirection;
2426 
2427  scalarField normalizedPosition(pointLabels.size(), Zero);
2428  forAll(pointLabels, i)
2429  {
2430  label pointi = pointLabels[i];
2431  normalizedPosition[i] =
2432  (
2433  ((baseMesh.points()[pointi]&userDirection) - startPosition)
2434  / totalLength
2435  );
2436  }
2437 
2438  // Sort the normalized position
2439  labelList order(sortedOrder(normalizedPosition));
2440 
2441  DynamicList<scalar> seedPointDist;
2442 
2443  // Select points from finest refinement (one point-per plane)
2444  scalar prevDist = GREAT;
2445  forAll(order, i)
2446  {
2447  label pointi = order[i];
2448  scalar curDist = normalizedPosition[pointi];
2449  if (mag(curDist - prevDist) > meshRefiner_.mergeDistance())
2450  {
2451  seedPointDist.append(curDist);
2452  prevDist = curDist;
2453  }
2454  }
2455 
2456  // Collect data from all processors
2457  scalarList allSeedPointDist;
2458  {
2459  List<scalarList> gatheredDist(Pstream::nProcs());
2460  gatheredDist[Pstream::myProcNo()] = seedPointDist;
2461  Pstream::gatherList(gatheredDist);
2462 
2463  // Combine processor lists into one big list.
2464  allSeedPointDist =
2465  ListListOps::combine<scalarList>
2466  (
2467  gatheredDist, accessOp<scalarList>()
2468  );
2469  }
2470 
2471  // Pre-set the points not to smooth (after expansion)
2472  bitSet isFrozenPoint(baseMesh.nPoints(), true);
2473 
2474  {
2475  scalar minSeed = min(allSeedPointDist);
2476  scalar maxSeed = max(allSeedPointDist);
2477  Pstream::broadcasts(UPstream::worldComm, minSeed, maxSeed);
2478 
2479  forAll(normalizedPosition, posI)
2480  {
2481  const scalar pos = normalizedPosition[posI];
2482  if
2483  (
2484  (mag(pos-minSeed) < meshRefiner_.mergeDistance())
2485  || (mag(pos-maxSeed) < meshRefiner_.mergeDistance())
2486  )
2487  {
2488  // Boundary point: freeze
2489  isFrozenPoint.set(pointLabels[posI]);
2490  }
2491  else
2492  {
2493  // Internal to moving region
2494  isFrozenPoint.unset(pointLabels[posI]);
2495  }
2496  }
2497  }
2498 
2499  Info<< "Smoothing " << geometry[surfi].name() << ':' << nl
2500  << " Direction : " << userDirection << nl
2501  << " Number of points : "
2502  << returnReduce(pointLabels.size(), sumOp<label>())
2503  << " (out of " << baseMesh.globalData().nTotalPoints()
2504  << ")" << nl
2505  << " Smooth expansion iterations : "
2506  << shells.nSmoothExpansion()[shellI] << nl
2507  << " Smooth position iterations : "
2508  << shells.nSmoothPosition()[shellI] << nl
2509  << " Number of planes : "
2510  << allSeedPointDist.size()
2511  << endl;
2512 
2513  // Make lookup from original normalized distance to new value
2514  List<Tuple2<scalar, scalar>> keyAndValue(allSeedPointDist.size());
2515 
2516  // Filter unique seed distances and iterate for user given steps
2517  // or until convergence. Then get back map from old to new distance
2518  if (Pstream::master())
2519  {
2520  mergeAndSmoothRatio
2521  (
2522  allSeedPointDist,
2523  shells.nSmoothExpansion()[shellI],
2524  keyAndValue
2525  );
2526  }
2527 
2528  Pstream::broadcast(keyAndValue);
2529 
2530  // Construct an iterpolation table for further queries
2531  // - although normalized values are used for query,
2532  // it might flow out of bounds due to precision, hence clamped
2533  const interpolationTable<scalar> table
2534  (
2535  keyAndValue,
2536  bounds::repeatableBounding::CLAMP,
2537  "undefined"
2538  );
2539 
2540  // Now move the points directly on the baseMesh.
2541  pointField baseNewPoints(baseMesh.points());
2542  forAll(pointLabels, i)
2543  {
2544  label pointi = pointLabels[i];
2545  const point& curPoint = baseMesh.points()[pointi];
2546  scalar curDist = normalizedPosition[i];
2547  scalar newDist = table(curDist);
2548  scalar newPosition = startPosition + newDist*totalLength;
2549  baseNewPoints[pointi] +=
2550  userDirection * (newPosition - (curPoint &userDirection));
2551  }
2552 
2553  // Moving base mesh with expansion ratio smoothing
2554  vectorField disp(baseNewPoints-baseMesh.points());
2555  syncTools::syncPointList
2556  (
2557  baseMesh,
2558  disp,
2559  maxMagSqrEqOp<vector>(),
2560  point::zero
2561  );
2562  baseMesh.movePoints(baseMesh.points()+disp);
2563 
2564  // Reset any moving state
2565  baseMesh.moving(false);
2566 
2567  if (debug&meshRefinement::MESH)
2568  {
2569  const_cast<Time&>(baseMesh.time())++;
2570 
2571  Pout<< "Writing directional expansion ratio smoothed"
2572  << " mesh to time " << meshRefiner_.timeName() << endl;
2573 
2574  meshRefiner_.write
2575  (
2576  meshRefinement::debugType(debug),
2577  meshRefinement::writeType
2578  (
2579  meshRefinement::writeLevel()
2580  | meshRefinement::WRITEMESH
2581  ),
2582  baseMesh.time().path()/meshRefiner_.timeName()
2583  );
2584  }
2585 
2586  // Now we have moved the points in user specified region. Smooth
2587  // them with neighbour points to avoid skewed cells
2588  // Instead of moving actual mesh, operate on copy
2589  pointField baseMeshPoints(baseMesh.points());
2590  scalar initResidual = 0.0;
2591  scalar prevIterResidual = GREAT;
2592  for (iter = 0; iter < shells.nSmoothPosition()[shellI]; iter++)
2593  {
2594  {
2595  const edgeList& edges = baseMesh.edges();
2596  const labelListList& pointEdges = baseMesh.pointEdges();
2597 
2598  pointField unsmoothedPoints(baseMeshPoints);
2599 
2600  scalarField sumOther(baseMesh.nPoints(), Zero);
2601  labelList nSumOther(baseMesh.nPoints(), Zero);
2602  labelList nSumXEdges(baseMesh.nPoints(), Zero);
2603  forAll(edges, edgei)
2604  {
2605  const edge& e = edges[edgei];
2606  sumOther[e[0]] +=
2607  (unsmoothedPoints[e[1]]&userDirection);
2608  nSumOther[e[0]]++;
2609  sumOther[e[1]] +=
2610  (unsmoothedPoints[e[0]]&userDirection);
2611  nSumOther[e[1]]++;
2612  if (isXEdge[edgei])
2613  {
2614  nSumXEdges[e[0]]++;
2615  nSumXEdges[e[1]]++;
2616  }
2617  }
2618 
2619  syncTools::syncPointList
2620  (
2621  baseMesh,
2622  nSumXEdges,
2623  plusEqOp<label>(),
2624  label(0)
2625  );
2626 
2627  forAll(pointLabels, i)
2628  {
2629  label pointi = pointLabels[i];
2630 
2631  if (nSumXEdges[pointi] < 2)
2632  {
2633  // Hanging node. Remove the (single!) X edge so it
2634  // will follow points above or below instead
2635  const labelList& pEdges = pointEdges[pointi];
2636  forAll(pEdges, pE)
2637  {
2638  label edgei = pEdges[pE];
2639  if (isXEdge[edgei])
2640  {
2641  const edge& e = edges[edgei];
2642  label otherPt = e.otherVertex(pointi);
2643  nSumOther[pointi]--;
2644  sumOther[pointi] -=
2645  (
2646  unsmoothedPoints[otherPt]
2647  & userDirection
2648  );
2649  }
2650  }
2651  }
2652  }
2653 
2654  syncTools::syncPointList
2655  (
2656  baseMesh,
2657  sumOther,
2658  plusEqOp<scalar>(),
2659  scalar(0)
2660  );
2661  syncTools::syncPointList
2662  (
2663  baseMesh,
2664  nSumOther,
2665  plusEqOp<label>(),
2666  label(0)
2667  );
2668 
2669  forAll(pointLabels, i)
2670  {
2671  label pointi = pointLabels[i];
2672 
2673  if ((nSumOther[pointi] >= 2) && !isFrozenPoint[pointi])
2674  {
2675  scalar smoothPos =
2676  0.5
2677  *(
2678  (unsmoothedPoints[pointi]&userDirection)
2679  +sumOther[pointi]/nSumOther[pointi]
2680  );
2681 
2682  vector& v = baseNewPoints[pointi];
2683  v += (smoothPos-(v&userDirection))*userDirection;
2684  }
2685  }
2686 
2687  const vectorField residual(baseNewPoints - baseMeshPoints);
2688  {
2689  scalar res(gSum(mag(residual)));
2690 
2691  if (iter == 0)
2692  {
2693  initResidual = res;
2694  }
2695  res /= initResidual;
2696 
2697  if (mag(prevIterResidual - res) < SMALL)
2698  {
2699  Info<< "Converged smoothing in iteration " << iter
2700  << " initResidual: " << initResidual
2701  << " final residual : " << res << endl;
2702  break;
2703  }
2704  else
2705  {
2706  prevIterResidual = res;
2707  }
2708  }
2709 
2710  // Just copy new location instead of moving base mesh
2711  baseMeshPoints = baseNewPoints;
2712  }
2713  }
2714 
2715  // Move mesh to new location
2716  vectorField dispSmooth(baseMeshPoints-baseMesh.points());
2717  syncTools::syncPointList
2718  (
2719  baseMesh,
2720  dispSmooth,
2721  maxMagSqrEqOp<vector>(),
2722  point::zero
2723  );
2724  baseMesh.movePoints(baseMesh.points()+dispSmooth);
2725 
2726  // Reset any moving state
2727  baseMesh.moving(false);
2728 
2729  if (debug&meshRefinement::MESH)
2730  {
2731  const_cast<Time&>(baseMesh.time())++;
2732 
2733  Pout<< "Writing positional smoothing iteration "
2734  << iter << " mesh to time " << meshRefiner_.timeName()
2735  << endl;
2736  meshRefiner_.write
2737  (
2738  meshRefinement::debugType(debug),
2739  meshRefinement::writeType
2740  (
2741  meshRefinement::writeLevel()
2742  | meshRefinement::WRITEMESH
2743  ),
2744  baseMesh.time().path()/meshRefiner_.timeName()
2745  );
2746  }
2747  }
2748  }
2749  return iter;
2750 }
2751 
2752 
2753 void Foam::snappyRefineDriver::baffleAndSplitMesh
2754 (
2755  const refinementParameters& refineParams,
2756  const snapParameters& snapParams,
2757  const bool handleSnapProblems,
2758  const dictionary& motionDict
2759 )
2760 {
2761  if (dryRun_)
2762  {
2763  return;
2764  }
2765 
2766  addProfiling(split, "snappyHexMesh::refine::splitting");
2767  Info<< nl
2768  << "Splitting mesh at surface intersections" << nl
2769  << "---------------------------------------" << nl
2770  << endl;
2771 
2772  const fvMesh& mesh = meshRefiner_.mesh();
2773 
2774  if (debug)
2775  {
2776  const_cast<Time&>(mesh.time())++;
2777  }
2778 
2779  // Introduce baffles at surface intersections. Note:
2780  // meshRefinement::surfaceIndex() will
2781  // be like boundary face from now on so not coupled anymore.
2782  meshRefiner_.baffleAndSplitMesh
2783  (
2784  handleSnapProblems, // detect&remove potential snap problem
2785 
2786  // Snap problem cell detection
2787  snapParams,
2788  refineParams.useTopologicalSnapDetection(),
2789  false, // perpendicular edge connected cells
2790  scalarField(0), // per region perpendicular angle
2791  refineParams.nErodeCellZone(),
2792 
2793  motionDict,
2794  const_cast<Time&>(mesh.time()),
2795  globalToMasterPatch_,
2796  globalToSlavePatch_,
2797  refineParams.locationsInMesh(),
2798  refineParams.zonesInMesh(),
2799  refineParams.locationsOutsideMesh(),
2800  !refineParams.useLeakClosure(),
2801  setFormatter_,
2802  surfFormatter_
2803  );
2804 
2805 
2806  if (!handleSnapProblems) // merge free standing baffles?
2807  {
2808  meshRefiner_.mergeFreeStandingBaffles
2809  (
2810  snapParams,
2811  refineParams.useTopologicalSnapDetection(),
2812  false, // perpendicular edge connected cells
2813  scalarField(0), // per region perpendicular angle
2814  refineParams.planarAngle(),
2815  motionDict,
2816  const_cast<Time&>(mesh.time()),
2817  globalToMasterPatch_,
2818  globalToSlavePatch_,
2819  refineParams.locationsInMesh(),
2820  refineParams.locationsOutsideMesh()
2821  );
2822  }
2823 }
2824 
2825 
2826 void Foam::snappyRefineDriver::zonify
2827 (
2828  const refinementParameters& refineParams,
2829  wordPairHashTable& zonesToFaceZone
2830 )
2831 {
2832  // Mesh is at its finest. Do zoning
2833  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2834  // This puts all faces with intersection across a zoneable surface
2835  // into that surface's faceZone. All cells inside faceZone get given the
2836  // same cellZone.
2837 
2838  const labelList namedSurfaces =
2839  surfaceZonesInfo::getNamedSurfaces(meshRefiner_.surfaces().surfZones());
2840 
2841  if
2842  (
2843  namedSurfaces.size()
2844  || refineParams.zonesInMesh().size()
2845  )
2846  {
2847  Info<< nl
2848  << "Introducing zones for interfaces" << nl
2849  << "--------------------------------" << nl
2850  << endl;
2851 
2852  const fvMesh& mesh = meshRefiner_.mesh();
2853 
2854  if (debug)
2855  {
2856  const_cast<Time&>(mesh.time())++;
2857  }
2858 
2859  meshRefiner_.zonify
2860  (
2861  refineParams.allowFreeStandingZoneFaces(),
2862  refineParams.nErodeCellZone(),
2863  refineParams.locationsInMesh(),
2864  refineParams.zonesInMesh(),
2865  refineParams.locationsOutsideMesh(),
2866  !refineParams.useLeakClosure(),
2867  setFormatter_,
2868  surfFormatter_,
2869  zonesToFaceZone
2870  );
2871 
2872  if (debug&meshRefinement::MESH)
2873  {
2874  Pout<< "Writing zoned mesh to time "
2875  << meshRefiner_.timeName() << endl;
2876  meshRefiner_.write
2877  (
2878  meshRefinement::debugType(debug),
2879  meshRefinement::writeType
2880  (
2881  meshRefinement::writeLevel()
2882  | meshRefinement::WRITEMESH
2883  ),
2884  mesh.time().path()/meshRefiner_.timeName()
2885  );
2886  }
2887 
2888  // Check that all faces are synced
2889  meshRefinement::checkCoupledFaceZones(mesh);
2890  }
2891 }
2892 
2893 
2894 void Foam::snappyRefineDriver::splitAndMergeBaffles
2895 (
2896  const refinementParameters& refineParams,
2897  const snapParameters& snapParams,
2898  const bool handleSnapProblems,
2899  const dictionary& motionDict
2900 )
2901 {
2902  if (dryRun_)
2903  {
2904  return;
2905  }
2906 
2907  Info<< nl
2908  << "Handling cells with snap problems" << nl
2909  << "---------------------------------" << nl
2910  << endl;
2911 
2912  const fvMesh& mesh = meshRefiner_.mesh();
2913 
2914  // Introduce baffles and split mesh
2915  if (debug)
2916  {
2917  const_cast<Time&>(mesh.time())++;
2918  }
2919 
2920  const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
2921 
2922  meshRefiner_.baffleAndSplitMesh
2923  (
2924  handleSnapProblems,
2925 
2926  // Snap problem cell detection
2927  snapParams,
2928  refineParams.useTopologicalSnapDetection(),
2929  handleSnapProblems, // remove perp edge connected cells
2930  perpAngle, // perp angle
2931  refineParams.nErodeCellZone(),
2932 
2933  motionDict,
2934  const_cast<Time&>(mesh.time()),
2935  globalToMasterPatch_,
2936  globalToSlavePatch_,
2937  refineParams.locationsInMesh(),
2938  refineParams.zonesInMesh(),
2939  refineParams.locationsOutsideMesh(),
2940  !refineParams.useLeakClosure(),
2941  setFormatter_,
2942  surfFormatter_
2943  );
2944 
2945  // Merge free-standing baffles always
2946  meshRefiner_.mergeFreeStandingBaffles
2947  (
2948  snapParams,
2949  refineParams.useTopologicalSnapDetection(),
2950  handleSnapProblems,
2951  perpAngle,
2952  refineParams.planarAngle(),
2953  motionDict,
2954  const_cast<Time&>(mesh.time()),
2955  globalToMasterPatch_,
2956  globalToSlavePatch_,
2957  refineParams.locationsInMesh(),
2958  refineParams.locationsOutsideMesh()
2959  );
2960 
2961  if (debug)
2962  {
2963  const_cast<Time&>(mesh.time())++;
2964  }
2965 
2966  // Duplicate points on baffles that are on more than one cell
2967  // region. This will help snapping pull them to separate surfaces.
2968  meshRefiner_.dupNonManifoldPoints();
2969 
2970 
2971  // Merge all baffles that are still remaining after duplicating points.
2972  List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
2973 
2974  const label nCouples = returnReduce(couples.size(), sumOp<label>());
2975 
2976  Info<< "Detected unsplittable baffles : " << nCouples << endl;
2977 
2978  if (nCouples > 0)
2979  {
2980  // Actually merge baffles. Note: not exactly parallellized. Should
2981  // convert baffle faces into processor faces if they resulted
2982  // from them.
2983  meshRefiner_.mergeBaffles(couples, Map<label>(0));
2984 
2985  if (debug)
2986  {
2987  // Debug:test all is still synced across proc patches
2988  meshRefiner_.checkData();
2989  }
2990 
2991  // Remove any now dangling parts
2992  meshRefiner_.splitMeshRegions
2993  (
2994  globalToMasterPatch_,
2995  globalToSlavePatch_,
2996  refineParams.locationsInMesh(),
2997  refineParams.locationsOutsideMesh(),
2998  true, // exit if any path to outside points
2999  setFormatter_
3000  );
3001 
3002  if (debug)
3003  {
3004  // Debug:test all is still synced across proc patches
3005  meshRefiner_.checkData();
3006  }
3007 
3008  Info<< "Merged free-standing baffles in = "
3009  << mesh.time().cpuTimeIncrement() << " s." << endl;
3010  }
3011 
3012  if (debug&meshRefinement::MESH)
3013  {
3014  Pout<< "Writing handleProblemCells mesh to time "
3015  << meshRefiner_.timeName() << endl;
3016  meshRefiner_.write
3017  (
3018  meshRefinement::debugType(debug),
3019  meshRefinement::writeType
3020  (
3021  meshRefinement::writeLevel()
3022  | meshRefinement::WRITEMESH
3023  ),
3024  mesh.time().path()/meshRefiner_.timeName()
3025  );
3026  }
3027 }
3028 
3029 
3030 void Foam::snappyRefineDriver::addFaceZones
3031 (
3032  meshRefinement& meshRefiner,
3033  const refinementParameters& refineParams,
3034  const HashTable<Pair<word>>& faceZoneToPatches
3035 )
3036 {
3037  if (faceZoneToPatches.size())
3038  {
3039  Info<< nl
3040  << "Adding patches for face zones" << nl
3041  << "-----------------------------" << nl
3042  << endl;
3043 
3044  Info<< setf(ios_base::left)
3045  << setw(6) << "Patch"
3046  << setw(20) << "Type"
3047  << setw(30) << "Name"
3048  << setw(30) << "FaceZone"
3049  << setw(10) << "FaceType"
3050  << nl
3051  << setw(6) << "-----"
3052  << setw(20) << "----"
3053  << setw(30) << "----"
3054  << setw(30) << "--------"
3055  << setw(10) << "--------"
3056  << endl;
3057 
3058  const polyMesh& mesh = meshRefiner.mesh();
3059 
3060  // Add patches for added inter-region faceZones
3061  forAllConstIters(faceZoneToPatches, iter)
3062  {
3063  const word& fzName = iter.key();
3064  const Pair<word>& patchNames = iter.val();
3065 
3066  // Get any user-defined faceZone data
3067  surfaceZonesInfo::faceZoneType fzType;
3068  dictionary patchInfo = refineParams.getZoneInfo(fzName, fzType);
3069 
3070  const word& masterName = fzName;
3071  //const word slaveName = fzName + "_slave";
3072  //const word slaveName = czNames.second()+"_to_"+czNames.first();
3073  const word& slaveName = patchNames.second();
3074 
3075  label mpi = meshRefiner.addMeshedPatch(masterName, patchInfo);
3076 
3077  Info<< setf(ios_base::left)
3078  << setw(6) << mpi
3079  << setw(20) << mesh.boundaryMesh()[mpi].type()
3080  << setw(30) << masterName
3081  << setw(30) << fzName
3082  << setw(10) << surfaceZonesInfo::faceZoneTypeNames[fzType]
3083  << nl;
3084 
3085 
3086  label sli = meshRefiner.addMeshedPatch(slaveName, patchInfo);
3087 
3088  Info<< setf(ios_base::left)
3089  << setw(6) << sli
3090  << setw(20) << mesh.boundaryMesh()[sli].type()
3091  << setw(30) << slaveName
3092  << setw(30) << fzName
3093  << setw(10) << surfaceZonesInfo::faceZoneTypeNames[fzType]
3094  << nl;
3095 
3096  meshRefiner.addFaceZone(fzName, masterName, slaveName, fzType);
3097  }
3098 
3099  Info<< endl;
3100  }
3101 }
3102 
3103 
3104 void Foam::snappyRefineDriver::mergePatchFaces
3105 (
3106  const meshRefinement::FaceMergeType mergeType,
3107  const refinementParameters& refineParams,
3108  const dictionary& motionDict
3109 )
3110 {
3111  if (dryRun_)
3112  {
3113  return;
3114  }
3115 
3116  addProfiling(merge, "snappyHexMesh::refine::merge");
3117  Info<< nl
3118  << "Merge refined boundary faces" << nl
3119  << "----------------------------" << nl
3120  << endl;
3121 
3122  const fvMesh& mesh = meshRefiner_.mesh();
3123 
3124  if
3125  (
3126  mergeType == meshRefinement::FaceMergeType::GEOMETRIC
3127  || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
3128  )
3129  {
3130  meshRefiner_.mergePatchFacesUndo
3131  (
3132  Foam::cos(degToRad(45.0)),
3133  Foam::cos(degToRad(45.0)),
3134  meshRefiner_.meshedPatches(),
3135  motionDict,
3136  labelList(mesh.nFaces(), -1),
3137  mergeType
3138  );
3139  }
3140  else
3141  {
3142  // Still merge refined boundary faces if all four are on same patch
3143  meshRefiner_.mergePatchFaces
3144  (
3145  Foam::cos(degToRad(45.0)),
3146  Foam::cos(degToRad(45.0)),
3147  4, // only merge faces split into 4
3148  meshRefiner_.meshedPatches(),
3149  meshRefinement::FaceMergeType::GEOMETRIC // no merge across patches
3150  );
3151  }
3152 
3153  if (debug)
3154  {
3155  meshRefiner_.checkData();
3156  }
3157 
3158  meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
3159 
3160  if (debug)
3161  {
3162  meshRefiner_.checkData();
3163  }
3164 }
3165 
3166 
3167 void Foam::snappyRefineDriver::deleteSmallRegions
3168 (
3169  const refinementParameters& refineParams
3170 )
3171 {
3172  const fvMesh& mesh = meshRefiner_.mesh();
3173  const cellZoneMesh& czm = mesh.cellZones();
3174 
3175  //const labelList zoneIDs
3176  //(
3177  // meshRefiner_.getZones
3178  // (
3179  // List<surfaceZonesInfo::faceZoneType> fzTypes
3180  // ({
3181  // surfaceZonesInfo::BAFFLE,
3182  // surfaceZonesInfo::BOUNDARY,
3183  // });
3184  // )
3185  //);
3186  const labelList zoneIDs(identity(mesh.faceZones().size()));
3187 
3188  // Get faceZone and patch(es) per face (or -1 if face not on faceZone)
3189  labelList faceZoneID;
3190  labelList ownPatch;
3191  labelList neiPatch;
3192  labelList nB; // local count of faces per faceZone
3193  meshRefiner_.getZoneFaces(zoneIDs, faceZoneID, ownPatch, neiPatch, nB);
3194 
3195 
3196  // Mark all faces on outside of zones. Note: assumes that faceZones
3197  // are consistent with the outside of cellZones ...
3198 
3199  boolList isBlockedFace(mesh.nFaces(), false);
3200  meshRefiner_.selectSeparatedCoupledFaces(isBlockedFace);
3201 
3202  forAll(ownPatch, facei)
3203  {
3204  if (ownPatch[facei] != -1)
3205  {
3206  isBlockedFace[facei] = true;
3207  }
3208  }
3209 
3210  // Map from cell to zone. Not that background cells are not -1 but
3211  // cellZones.size()
3212  labelList cellToZone(mesh.nCells(), czm.size());
3213  for (const auto& cz : czm)
3214  {
3215  UIndirectList<label>(cellToZone, cz) = cz.index();
3216  }
3217 
3218  // Walk to split into regions
3219  const regionSplit cellRegion(mesh, isBlockedFace);
3220 
3221  // Count number of cells per zone and per region
3222  labelList nCellsPerRegion(cellRegion.nRegions(), 0);
3223  labelList regionToZone(cellRegion.nRegions(), -2);
3224  labelList nCellsPerZone(czm.size()+1, 0);
3225  forAll(cellRegion, celli)
3226  {
3227  const label regioni = cellRegion[celli];
3228  const label zonei = cellToZone[celli];
3229 
3230  // Zone for this region
3231  regionToZone[regioni] = zonei;
3232 
3233  nCellsPerRegion[regioni]++;
3234  nCellsPerZone[zonei]++;
3235  }
3236  Pstream::listCombineReduce(nCellsPerRegion, plusEqOp<label>());
3237  Pstream::listCombineReduce(regionToZone, maxEqOp<label>());
3238  Pstream::listCombineReduce(nCellsPerZone, plusEqOp<label>());
3239 
3240 
3241  // Mark small regions. Note that all processors have the same information
3242  // so should take the same decision. Alternative is to do master only
3243  // and scatter the decision.
3244  forAll(nCellsPerRegion, regioni)
3245  {
3246  const label zonei = regionToZone[regioni];
3247  label& nRegionCells = nCellsPerRegion[regioni];
3248 
3249  if
3250  (
3251  nRegionCells < refineParams.minCellFraction()*nCellsPerZone[zonei]
3252  || nRegionCells < refineParams.nMinCells()
3253  )
3254  {
3255  Info<< "Deleting region " << regioni
3256  << " (size " << nRegionCells
3257  << ") of zone size " << nCellsPerZone[zonei]
3258  << endl;
3259 
3260  // Mark region to be deleted. 0 size (= global) should never
3261  // occur.
3262  nRegionCells = 0;
3263  }
3264  }
3265 
3266  DynamicList<label> cellsToRemove(mesh.nCells()/128);
3267  forAll(cellRegion, celli)
3268  {
3269  if (nCellsPerRegion[cellRegion[celli]] == 0)
3270  {
3271  cellsToRemove.append(celli);
3272  }
3273  }
3274  const label nTotCellsToRemove = returnReduce
3275  (
3276  cellsToRemove.size(),
3277  sumOp<label>()
3278  );
3279  if (nTotCellsToRemove > 0)
3280  {
3281  Info<< "Deleting " << nTotCellsToRemove
3282  << " cells in small regions" << endl;
3283 
3284  removeCells cellRemover(mesh);
3285 
3286  cellsToRemove.shrink();
3287  const labelList exposedFaces
3288  (
3289  cellRemover.getExposedFaces(cellsToRemove)
3290  );
3291  const labelList exposedPatch
3292  (
3293  UIndirectList<label>(ownPatch, exposedFaces)
3294  );
3295  (void)meshRefiner_.doRemoveCells
3296  (
3297  cellsToRemove,
3298  exposedFaces,
3299  exposedPatch,
3300  cellRemover
3301  );
3302  }
3303 }
3304 
3305 
3306 void Foam::snappyRefineDriver::doRefine
3307 (
3308  const dictionary& refineDict,
3309  const refinementParameters& refineParams,
3310  const snapParameters& snapParams,
3311  const bool prepareForSnapping,
3312  const meshRefinement::FaceMergeType mergeType,
3313  const dictionary& motionDict
3314 )
3315 {
3316  addProfiling(refine, "snappyHexMesh::refine");
3317  Info<< nl
3318  << "Refinement phase" << nl
3319  << "----------------" << nl
3320  << endl;
3321 
3322  const fvMesh& mesh = meshRefiner_.mesh();
3323 
3324 
3325  // Check that all the keep points are inside the mesh.
3326  refineParams.findCells(true, mesh, refineParams.locationsInMesh());
3327 
3328  // Check that all the keep points are inside the mesh.
3329  if (dryRun_)
3330  {
3331  refineParams.findCells(true, mesh, refineParams.locationsOutsideMesh());
3332  }
3333 
3334  // Estimate cell sizes
3335  if (dryRun_)
3336  {
3337  snappyVoxelMeshDriver voxelDriver
3338  (
3339  meshRefiner_,
3340  globalToMasterPatch_,
3341  globalToSlavePatch_
3342  );
3343  voxelDriver.doRefine(refineParams);
3344  }
3345 
3346 
3347  // Refine around feature edges
3348  featureEdgeRefine
3349  (
3350  refineParams,
3351  100, // maxIter
3352  0 // min cells to refine
3353  );
3354 
3355 
3356  // Initial automatic gap-level refinement: gaps smaller than the cell size
3357  if
3358  (
3359  max(meshRefiner_.surfaces().maxGapLevel()) > 0
3360  || max(meshRefiner_.shells().maxGapLevel()) > 0
3361  || max(meshRefiner_.surfaces().maxCurvatureLevel()) > 0
3362  )
3363  {
3364  // In case we use automatic gap level refinement do some pre-refinement
3365  // (fast) since it is so slow.
3366 
3367  // Refine based on surface
3368  surfaceOnlyRefine
3369  (
3370  refineParams,
3371  20, // maxIter
3372  labelMax // no leak detection yet since all in same cell
3373  );
3374 
3375  // Refine cells that contain a gap
3376  smallFeatureRefine
3377  (
3378  refineParams,
3379  100 // maxIter
3380  );
3381  }
3382 
3383 
3384  // Refine based on surface
3385  surfaceOnlyRefine
3386  (
3387  refineParams,
3388  100, // maxIter
3389  0 // Iteration to start leak detection and closure
3390  );
3391 
3392  // Pass1 of automatic gap-level refinement: surface-intersected cells
3393  // in narrow gaps. Done early so we can remove the inside
3394  gapOnlyRefine
3395  (
3396  refineParams,
3397  100 // maxIter
3398  );
3399 
3400  // Remove cells inbetween two surfaces
3401  surfaceProximityBlock
3402  (
3403  refineParams,
3404  1 //100 // maxIter
3405  );
3406 
3407  // Remove cells (a certain distance) beyond surface intersections
3408  removeInsideCells
3409  (
3410  refineParams,
3411  1 // nBufferLayers
3412  );
3413 
3414  // Pass2 of automatic gap-level refinement: all cells in gaps
3415  bigGapOnlyRefine
3416  (
3417  refineParams,
3418  true, // spreadGapSize
3419  100 // maxIter
3420  );
3421 
3422  // Internal mesh refinement
3423  shellRefine
3424  (
3425  refineParams,
3426  100 // maxIter
3427  );
3428 
3429  // Remove any extra cells from limitRegion with level -1, without
3430  // adding any buffer layer this time
3431  removeInsideCells
3432  (
3433  refineParams,
3434  0 // nBufferLayers
3435  );
3436 
3437  // Refine any hexes with 5 or 6 faces refined to make smooth edges
3438  danglingCellRefine
3439  (
3440  refineParams,
3441  21, // 1 coarse face + 5 refined faces
3442  100 // maxIter
3443  );
3444  danglingCellRefine
3445  (
3446  refineParams,
3447  24, // 0 coarse faces + 6 refined faces
3448  100 // maxIter
3449  );
3450 
3451  // Refine any cells with differing refinement level on either side
3452  refinementInterfaceRefine
3453  (
3454  refineParams,
3455  10 // maxIter
3456  );
3457 
3458  // Directional shell refinement
3459  directionalShellRefine
3460  (
3461  refineParams,
3462  100 // maxIter
3463  );
3464 
3465  // Block gaps (always, ignore surface leakLevel)
3466  if (refineParams.locationsOutsideMesh().size())
3467  {
3468  // For now: only check leaks on meshed surfaces. The problem is that
3469  // blockLeakFaces always generates baffles any not just faceZones ...
3470  const labelList unnamedSurfaces
3471  (
3472  surfaceZonesInfo::getUnnamedSurfaces
3473  (
3474  meshRefiner_.surfaces().surfZones()
3475  )
3476  );
3477  if (unnamedSurfaces.size() && !dryRun_)
3478  {
3479  meshRefiner_.blockLeakFaces
3480  (
3481  globalToMasterPatch_,
3482  globalToSlavePatch_,
3483  refineParams.locationsInMesh(),
3484  refineParams.zonesInMesh(),
3485  refineParams.locationsOutsideMesh(),
3486  unnamedSurfaces
3487  );
3488  }
3489  }
3490 
3491  if
3492  (
3493  max(meshRefiner_.shells().nSmoothExpansion()) > 0
3494  || max(meshRefiner_.shells().nSmoothPosition()) > 0
3495  )
3496  {
3497  directionalSmooth(refineParams);
3498  }
3499 
3500 
3501  // Introduce baffles at surface intersections. Remove sections unreachable
3502  // from keepPoint.
3503  baffleAndSplitMesh
3504  (
3505  refineParams,
3506  snapParams,
3507  prepareForSnapping,
3508  motionDict
3509  );
3510 
3511 
3512  //- Commented out for now since causes zoning errors (sigsegv) on
3513  // case with faceZones. TBD.
3516  //boundaryRefinementInterfaceRefine
3517  //(
3518  // refineParams,
3519  // 10 // maxIter
3520  //);
3521 
3522 
3523  // Mesh is at its finest. Do optional zoning (cellZones and faceZones)
3524  wordPairHashTable zonesToFaceZone;
3525  zonify(refineParams, zonesToFaceZone);
3526 
3527  // Create pairs of patches for faceZones
3528  {
3529  HashTable<Pair<word>> faceZoneToPatches(zonesToFaceZone.size());
3530 
3531  // Note: zonesToFaceZone contains the same data on different
3532  // processors but in different order. We could sort the
3533  // contents but instead just loop in sortedToc order.
3534  List<Pair<word>> czs(zonesToFaceZone.sortedToc());
3535 
3536  forAll(czs, i)
3537  {
3538  const Pair<word>& czNames = czs[i];
3539  const word& fzName = zonesToFaceZone[czNames];
3540 
3541  const word& masterName = fzName;
3542  const word slaveName = czNames.second() + "_to_" + czNames.first();
3543  Pair<word> patches(masterName, slaveName);
3544  faceZoneToPatches.insert(fzName, patches);
3545  }
3546  addFaceZones(meshRefiner_, refineParams, faceZoneToPatches);
3547  }
3548 
3549  // Pull baffles apart
3550  splitAndMergeBaffles
3551  (
3552  refineParams,
3553  snapParams,
3554  prepareForSnapping,
3555  motionDict
3556  );
3557 
3558  // Do something about cells with refined faces on the boundary
3559  if (prepareForSnapping)
3560  {
3561  mergePatchFaces(mergeType, refineParams, motionDict);
3562  }
3563 
3564 
3565  if (refineParams.minCellFraction() > 0 || refineParams.nMinCells() > 0)
3566  {
3567  // Some small disconnected bits of mesh might remain since at
3568  // this point faceZones have not been converted into e.g. baffles.
3569  // We don't know whether e.g. the baffles are reset to be cyclicAMI
3570  // thus reconnecting. For now check if there are any particularly
3571  // small regions.
3572  deleteSmallRegions(refineParams);
3573  }
3574 
3575 
3576  if (!dryRun_ && Pstream::parRun())
3577  {
3578  Info<< nl
3579  << "Doing final balancing" << nl
3580  << "---------------------" << nl
3581  << endl;
3582 
3583  // Do final balancing. Keep zoned faces on one processor since the
3584  // snap phase will convert them to baffles and this only works for
3585  // internal faces.
3586  meshRefiner_.balance
3587  (
3588  true, // keepZoneFaces
3589  false, // keepBaffles
3590  scalarField(mesh.nCells(), 1), // cellWeights
3591  decomposer_,
3592  distributor_
3593  );
3594  }
3595 }
3596 
3597 
3598 // ************************************************************************* //
Foam::scalarList
List< scalar > scalarList
List of scalar.
Definition: scalarList.H:32
pbm
const polyBoundaryMesh & pbm
Definition: checkPatchTopology.H:27
Foam::sumMag
dimensioned< typename typeOfMag< Type >::type > sumMag(const DimensionedField< Type, GeoMesh > &f1)
Definition: DimensionedFieldFunctions.C:300
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
Foam::cmptMax
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:236
zoneIDs
const labelIOList & zoneIDs
Definition: correctPhi.H:59
Foam::Time::path
fileName path() const
Return path = rootPath/caseName. Same as TimePaths::path()
Definition: Time.H:503
Foam::direction
uint8_t direction
Definition: direction.H:46
Foam::cellList
List< cell > cellList
List of cell.
Definition: cellListFwd.H:39
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
pointLabels
labelList pointLabels(nPoints, -1)
Foam::primitiveMesh::nPoints
label nPoints() const noexcept
Number of mesh points.
Definition: primitiveMeshI.H:30
Foam::bounds::normalBounding::CLAMP
Clamp value to the start/end value.
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1122
Foam::edgeList
List< edge > edgeList
List of edge.
Definition: edgeList.H:32
Foam::sortedOrder
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:40
unitConversion.H
Unit conversion functions.
Foam::nl
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:50
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531
Foam::polyBoundaryMesh::patchID
const labelList & patchID() const
Per boundary face label the patch index.
Definition: polyBoundaryMesh.C:514
Foam::UPstream::parRun
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:1049
Foam::meshRefinement::writeLevel
static writeType writeLevel()
Get/set write level.
Definition: meshRefinement.C:3777
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
interface
interfaceProperties interface(alpha1, U, thermo->transportPropertiesDict())
Foam::labelPairList
List< labelPair > labelPairList
List of labelPair.
Definition: labelPair.H:33
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:131
Foam::surfaceZonesInfo::faceZoneTypeNames
static const Enum< faceZoneType > faceZoneTypeNames
Definition: surfaceZonesInfo.H:95
Foam::expressions::Detail::nComponents
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
Definition: exprTraits.C:40
Foam::labelMin
constexpr label labelMin
Definition: label.H:54
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
Definition: UPstream.H:1074
Foam::meshRefinement::REMOVE
set value to -1 any face that was refined
Definition: meshRefinement.H:132
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:360
Foam::UPstream::worldComm
static label worldComm
Communicator for all ranks. May differ from commGlobal() if local worlds are in use.
Definition: UPstream.H:409
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:83
Foam::returnReduce
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Definition: PstreamReduceOps.H:463
Foam::labelListList
List< labelList > labelListList
List of labelList.
Definition: labelList.H:38
Foam::faceZoneID
DynamicID< faceZoneMesh > faceZoneID
Foam::faceZoneID.
Definition: ZoneIDs.H:43
Foam::refPtr
A class for managing references or pointers (no reference counting)
Definition: HashPtrTable.H:49
Foam::labelUList
UList< label > labelUList
A UList of labels.
Definition: UList.H:78
Foam::Pstream::broadcast
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all communicator ranks. Does nothing in non-paral...
Definition: PstreamBroadcast.C:28
Foam::localPointRegion::findDuplicateFacePairs
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
Definition: localPointRegion.C:618
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:421
Foam::faceList
List< face > faceList
List of faces.
Definition: faceListFwd.H:39
Foam::meshRefinement::updateList
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
Definition: meshRefinementTemplates.C:31
Foam::pos
dimensionedScalar pos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:170
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
Definition: UPstream.H:1065
Foam::fvMeshDistribute
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Definition: fvMeshDistribute.H:66
addProfiling
#define addProfiling(Name, Descr)
Define profiling trigger with specified name and description string.
Definition: profilingTrigger.H:131
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
Foam::Pstream::gatherList
static void gatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
Gather data, but keep individual values separate. Uses the specified communication schedule...
Definition: PstreamGatherList.C:41
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::labelIOList
IOList< label > labelIOList
IO for a List of label.
Definition: labelIOList.H:32
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:316
Foam::meshRefinement::FaceMergeType
FaceMergeType
Enumeration for what to do with co-planar patch faces on a single.
Definition: meshRefinement.H:140
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:575
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:258
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Definition: exprTraits.C:127
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Definition: labelLists.C:44
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Definition: polyMesh.H:608
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
Foam::UPtrList::size
label size() const noexcept
The number of entries in the list.
Definition: UPtrListI.H:106
patchNames
wordList patchNames(nPatches)
Foam::coordSetWriter
Base class for writing coordSet(s) and tracks with fields.
Definition: coordSetWriter.H:82
Foam::meshRefinement::checkCoupledFaceZones
static void checkCoupledFaceZones(const polyMesh &)
Helper function: check that face zones are synced.
Definition: meshRefinement.C:2222
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1116
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam::setf
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:169
Foam::polyMesh::globalData
const globalMeshData & globalData() const
Return parallel info (demand-driven)
Definition: polyMesh.C:1311
Foam::decompositionMethod
Abstract base class for domain decomposition.
Definition: decompositionMethod.H:47
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:71
Foam::vector
Vector< scalar > vector
Definition: vector.H:57
Foam::volumeType::INSIDE
A location inside the volume.
Definition: volumeType.H:65
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1103
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:26
Foam::surfaceZonesInfo::getUnnamedSurfaces
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Definition: surfaceZonesInfo.C:257
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition: syncToolsTemplates.C:753
Foam::fieldTypes::surface
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
IOmanip.H
Istream and Ostream manipulators taking arguments.
Foam::ensightOutput::debug
int debug
Static debugging option.
Foam::polyMesh::faceZones
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:670
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
f
labelList f(nPoints)
Foam::globalMeshData::nTotalCells
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
Definition: globalMeshData.H:463
Foam::cpuTimePosix::cpuTimeIncrement
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTimePosix.C:80
Foam::snappyRefineDriver::addFaceZones
static void addFaceZones(meshRefinement &meshRefiner, const refinementParameters &refineParams, const HashTable< Pair< word >> &faceZoneToPatches)
Helper: add faceZones and patches.
Definition: snappyRefineDriver.C:3024
Foam::Pstream::broadcasts
static void broadcasts(const label comm, Type &arg1, Args &&... args)
Broadcast multiple items to all communicator ranks. Does nothing in non-parallel. ...
Definition: PstreamBroadcast.C:58
Foam::meshRefinement
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
Definition: meshRefinement.H:83
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:1433
split
static bool split(const std::string &line, std::string &key, std::string &val)
Definition: cpuInfo.C:32
Foam::snappyRefineDriver::doRefine
void doRefine(const dictionary &refineDict, const refinementParameters &refineParams, const snapParameters &snapParams, const bool prepareForSnapping, const meshRefinement::FaceMergeType mergeType, const dictionary &motionDict)
Do all the refinement.
Definition: snappyRefineDriver.C:3300
Foam::point
vector point
Point is a vector.
Definition: point.H:37
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:89
Foam::primitiveMesh::faceAreas
const vectorField & faceAreas() const
Definition: primitiveMeshFaceCentresAndAreas.C:83
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::UPstream::master
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
Definition: UPstream.H:1082
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:54
Foam::polyMesh::cellZones
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
Definition: polyMesh.H:678
Foam::surfaceZonesInfo::getNamedSurfaces
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
Definition: surfaceZonesInfo.C:278
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::labelMax
constexpr label labelMax
Definition: label.H:55
Foam::VectorSpace< Vector< label >, label, 3 >::uniform
static Vector< label > uniform(const label &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:157
Foam::meshRefinement::writeType
writeType
Enumeration for what to write. Used as a bit-pattern.
Definition: meshRefinement.H:114
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:48
Foam::surfaceZonesInfo::faceZoneType
faceZoneType
What to do with faceZone faces.
Definition: surfaceZonesInfo.H:88
Foam::meshRefinement::debugType
debugType
Enumeration for what to debug. Used as a bit-pattern.
Definition: meshRefinement.H:92
Foam::labelVector
Vector< label > labelVector
Vector of labels.
Definition: labelVector.H:47
Foam::setw
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:62
Foam::UILList::first
T * first()
The first entry in the list.
Definition: UILList.H:542
coupled
bool coupled
Definition: createCoupledRegions.H:1
Foam::returnReduceOr
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
Definition: PstreamReduceOps.H:499
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:44
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:60
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::wordPairHashTable
HashTable< word, wordPair, Foam::Hash< wordPair > > wordPairHashTable
HashTable of wordPair.
Definition: wordPairHashes.H:55
Foam::cellZoneMesh
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
Definition: cellZoneMeshFwd.H:38
pp
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:26
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
Foam::Pstream::listCombineReduce
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
Definition: PstreamCombineGather.C:530
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127