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