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

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

Generated by   doxygen 1.8.14