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

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

Generated by   doxygen 1.8.14