snappyRefineDriver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "snappyRefineDriver.H"
27 #include "meshRefinement.H"
28 #include "fvMesh.H"
29 #include "Time.H"
30 #include "cellSet.H"
31 #include "syncTools.H"
32 #include "refinementParameters.H"
33 #include "refinementSurfaces.H"
34 #include "refinementFeatures.H"
35 #include "refinementRegions.H"
36 #include "polyDistributionMap.H"
37 #include "unitConversion.H"
38 #include "snapParameters.H"
39 #include "localPointRegion.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45 
47 
48 } // End namespace Foam
49 
50 
51 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
52 
53 // Construct from components
55 (
56  meshRefinement& meshRefiner,
57  decompositionMethod& decomposer,
58  fvMeshDistribute& distributor,
59  const labelList& globalToMasterPatch,
60  const labelList& globalToSlavePatch
61 )
62 :
63  meshRefiner_(meshRefiner),
64  decomposer_(decomposer),
65  distributor_(distributor),
66  globalToMasterPatch_(globalToMasterPatch),
67  globalToSlavePatch_(globalToSlavePatch)
68 {}
69 
70 
71 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
72 
73 Foam::label Foam::snappyRefineDriver::featureEdgeRefine
74 (
75  const refinementParameters& refineParams,
76  const label maxIter,
77  const label minRefine
78 )
79 {
80  const fvMesh& mesh = meshRefiner_.mesh();
81 
82  label iter = 0;
83 
84  if (meshRefiner_.features().size() && maxIter > 0)
85  {
86  for (; iter < maxIter; iter++)
87  {
88  Info<< nl
89  << "Feature refinement iteration " << iter << nl
90  << "------------------------------" << nl
91  << endl;
92 
93  labelList candidateCells
94  (
95  meshRefiner_.refineCandidates
96  (
97  refineParams.selectionPoints().inside(),
98  refineParams.curvature(),
99  refineParams.planarAngle(),
100 
101  true, // featureRefinement
102  false, // featureDistanceRefinement
103  false, // internalRefinement
104  false, // surfaceRefinement
105  false, // curvatureRefinement
106  false, // gapRefinement
107  refineParams.maxGlobalCells(),
108  refineParams.maxLocalCells()
109  )
110  );
111  labelList cellsToRefine
112  (
113  meshRefiner_.meshCutter().consistentRefinement
114  (
115  candidateCells,
116  true
117  )
118  );
119  Info<< "Determined cells to refine in = "
120  << mesh.time().cpuTimeIncrement() << " s" << endl;
121 
122 
123 
124  label nCellsToRefine = cellsToRefine.size();
125  reduce(nCellsToRefine, sumOp<label>());
126 
127  Info<< "Selected for feature refinement : " << nCellsToRefine
128  << " cells (out of " << mesh.globalData().nTotalCells()
129  << ')' << endl;
130 
131  if (nCellsToRefine <= minRefine)
132  {
133  Info<< "Stopping refining since too few cells selected."
134  << nl << endl;
135  break;
136  }
137 
138 
139  if (debug > 0)
140  {
141  const_cast<Time&>(mesh.time())++;
142  }
143 
144 
145  if
146  (
148  (
149  (mesh.nCells() >= refineParams.maxLocalCells()),
150  orOp<bool>()
151  )
152  )
153  {
154  meshRefiner_.balanceAndRefine
155  (
156  "feature refinement iteration " + name(iter),
157  decomposer_,
158  distributor_,
159  cellsToRefine,
160  refineParams.maxLoadUnbalance()
161  );
162  }
163  else
164  {
165  meshRefiner_.refineAndBalance
166  (
167  "feature refinement iteration " + name(iter),
168  decomposer_,
169  distributor_,
170  cellsToRefine,
171  refineParams.maxLoadUnbalance()
172  );
173  }
174  }
175  }
176  return iter;
177 }
178 
179 
180 Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
181 (
182  const refinementParameters& refineParams,
183  const label maxIter
184 )
185 {
186  const fvMesh& mesh = meshRefiner_.mesh();
187 
188  // Determine the maximum refinement level over all surfaces. This
189  // determines the minimum number of surface refinement iterations.
190  label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
191 
192  label iter;
193  for (iter = 0; iter < maxIter; iter++)
194  {
195  Info<< nl
196  << "Surface refinement iteration " << iter << nl
197  << "------------------------------" << nl
198  << endl;
199 
200 
201  // Determine cells to refine
202  // ~~~~~~~~~~~~~~~~~~~~~~~~~
203  // Only look at surface intersections (minLevel and surface curvature),
204  // do not do internal refinement (refinementShells)
205 
206  labelList candidateCells
207  (
208  meshRefiner_.refineCandidates
209  (
210  refineParams.selectionPoints().inside(),
211  refineParams.curvature(),
212  refineParams.planarAngle(),
213 
214  false, // featureRefinement
215  false, // featureDistanceRefinement
216  false, // internalRefinement
217  true, // surfaceRefinement
218  true, // curvatureRefinement
219  false, // gapRefinement
220  refineParams.maxGlobalCells(),
221  refineParams.maxLocalCells()
222  )
223  );
224  labelList cellsToRefine
225  (
226  meshRefiner_.meshCutter().consistentRefinement
227  (
228  candidateCells,
229  true
230  )
231  );
232  Info<< "Determined cells to refine in = "
233  << mesh.time().cpuTimeIncrement() << " s" << endl;
234 
235 
236  label nCellsToRefine = cellsToRefine.size();
237  reduce(nCellsToRefine, sumOp<label>());
238 
239  Info<< "Selected for refinement : " << nCellsToRefine
240  << " cells (out of " << mesh.globalData().nTotalCells()
241  << ')' << endl;
242 
243  // Stop when no cells to refine or have done minimum necessary
244  // iterations and not enough cells to refine.
245  if
246  (
247  nCellsToRefine == 0
248  || (
249  iter >= overallMaxLevel
250  && nCellsToRefine <= refineParams.minRefineCells()
251  )
252  )
253  {
254  Info<< "Stopping refining since too few cells selected."
255  << nl << endl;
256  break;
257  }
258 
259 
260  if (debug)
261  {
262  const_cast<Time&>(mesh.time())++;
263  }
264 
265 
266  if
267  (
269  (
270  (mesh.nCells() >= refineParams.maxLocalCells()),
271  orOp<bool>()
272  )
273  )
274  {
275  meshRefiner_.balanceAndRefine
276  (
277  "surface refinement iteration " + name(iter),
278  decomposer_,
279  distributor_,
280  cellsToRefine,
281  refineParams.maxLoadUnbalance()
282  );
283  }
284  else
285  {
286  meshRefiner_.refineAndBalance
287  (
288  "surface refinement iteration " + name(iter),
289  decomposer_,
290  distributor_,
291  cellsToRefine,
292  refineParams.maxLoadUnbalance()
293  );
294  }
295  }
296  return iter;
297 }
298 
299 
300 Foam::label Foam::snappyRefineDriver::gapOnlyRefine
301 (
302  const refinementParameters& refineParams,
303  const label maxIter
304 )
305 {
306  const fvMesh& mesh = meshRefiner_.mesh();
307 
308  // Determine the maximum refinement level over all surfaces. This
309  // determines the minimum number of surface refinement iterations.
310 
311  label maxIncrement = 0;
312  const labelList& maxLevel = meshRefiner_.surfaces().maxLevel();
313  const labelList& gapLevel = meshRefiner_.surfaces().gapLevel();
314 
315  forAll(maxLevel, i)
316  {
317  maxIncrement = max(maxIncrement, gapLevel[i]-maxLevel[i]);
318  }
319 
320  label iter = 0;
321 
322  if (maxIncrement == 0)
323  {
324  return iter;
325  }
326 
327  for (iter = 0; iter < maxIter; iter++)
328  {
329  Info<< nl
330  << "Gap refinement iteration " << iter << nl
331  << "--------------------------" << nl
332  << endl;
333 
334 
335  // Determine cells to refine
336  // ~~~~~~~~~~~~~~~~~~~~~~~~~
337  // Only look at surface intersections (minLevel and surface curvature),
338  // do not do internal refinement (refinementShells)
339 
340  labelList candidateCells
341  (
342  meshRefiner_.refineCandidates
343  (
344  refineParams.selectionPoints().inside(),
345  refineParams.curvature(),
346  refineParams.planarAngle(),
347 
348  false, // featureRefinement
349  false, // featureDistanceRefinement
350  false, // internalRefinement
351  false, // surfaceRefinement
352  false, // curvatureRefinement
353  true, // gapRefinement
354  refineParams.maxGlobalCells(),
355  refineParams.maxLocalCells()
356  )
357  );
358 
359  if (debug&meshRefinement::MESH)
360  {
361  Pout<< "Dumping " << candidateCells.size()
362  << " cells to cellSet candidateCellsFromGap." << endl;
363  cellSet c(mesh, "candidateCellsFromGap", candidateCells);
364  c.instance() = meshRefiner_.name();
365  c.write();
366  }
367 
368  // Grow by one layer to make sure we're covering the gap
369  {
370  boolList isCandidateCell(mesh.nCells(), false);
371  forAll(candidateCells, i)
372  {
373  isCandidateCell[candidateCells[i]] = true;
374  }
375 
376  for (label i=0; i<1; i++)
377  {
378  boolList newIsCandidateCell(isCandidateCell);
379 
380  // Internal faces
381  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
382  {
383  label own = mesh.faceOwner()[facei];
384  label nei = mesh.faceNeighbour()[facei];
385 
386  if (isCandidateCell[own] != isCandidateCell[nei])
387  {
388  newIsCandidateCell[own] = true;
389  newIsCandidateCell[nei] = true;
390  }
391  }
392 
393  // Get coupled boundary condition values
394  boolList neiisCandidateCell;
396  (
397  mesh,
398  isCandidateCell,
399  neiisCandidateCell
400  );
401 
402  // Boundary faces
403  for
404  (
405  label facei = mesh.nInternalFaces();
406  facei < mesh.nFaces();
407  facei++
408  )
409  {
410  label own = mesh.faceOwner()[facei];
411  label bFacei = facei-mesh.nInternalFaces();
412 
413  if (isCandidateCell[own] != neiisCandidateCell[bFacei])
414  {
415  newIsCandidateCell[own] = true;
416  }
417  }
418 
419  isCandidateCell.transfer(newIsCandidateCell);
420  }
421 
422  label n = 0;
423  forAll(isCandidateCell, celli)
424  {
425  if (isCandidateCell[celli])
426  {
427  n++;
428  }
429  }
430  candidateCells.setSize(n);
431  n = 0;
432  forAll(isCandidateCell, celli)
433  {
434  if (isCandidateCell[celli])
435  {
436  candidateCells[n++] = celli;
437  }
438  }
439  }
440 
441 
442  if (debug&meshRefinement::MESH)
443  {
444  Pout<< "Dumping " << candidateCells.size()
445  << " cells to cellSet candidateCellsFromGapPlusBuffer." << endl;
446  cellSet c(mesh, "candidateCellsFromGapPlusBuffer", candidateCells);
447  c.instance() = meshRefiner_.name();
448  c.write();
449  }
450 
451 
452  labelList cellsToRefine
453  (
454  meshRefiner_.meshCutter().consistentRefinement
455  (
456  candidateCells,
457  true
458  )
459  );
460  Info<< "Determined cells to refine in = "
461  << mesh.time().cpuTimeIncrement() << " s" << endl;
462 
463 
464  label nCellsToRefine = cellsToRefine.size();
465  reduce(nCellsToRefine, sumOp<label>());
466 
467  Info<< "Selected for refinement : " << nCellsToRefine
468  << " cells (out of " << mesh.globalData().nTotalCells()
469  << ')' << endl;
470 
471  // Stop when no cells to refine or have done minimum necessary
472  // iterations and not enough cells to refine.
473  if
474  (
475  nCellsToRefine == 0
476  || (
477  iter >= maxIncrement
478  && nCellsToRefine <= refineParams.minRefineCells()
479  )
480  )
481  {
482  Info<< "Stopping refining since too few cells selected."
483  << nl << endl;
484  break;
485  }
486 
487 
488  if (debug)
489  {
490  const_cast<Time&>(mesh.time())++;
491  }
492 
493 
494  if
495  (
497  (
498  (mesh.nCells() >= refineParams.maxLocalCells()),
499  orOp<bool>()
500  )
501  )
502  {
503  meshRefiner_.balanceAndRefine
504  (
505  "gap refinement iteration " + name(iter),
506  decomposer_,
507  distributor_,
508  cellsToRefine,
509  refineParams.maxLoadUnbalance()
510  );
511  }
512  else
513  {
514  meshRefiner_.refineAndBalance
515  (
516  "gap refinement iteration " + name(iter),
517  decomposer_,
518  distributor_,
519  cellsToRefine,
520  refineParams.maxLoadUnbalance()
521  );
522  }
523  }
524  return iter;
525 }
526 
527 
528 Foam::label Foam::snappyRefineDriver::danglingCellRefine
529 (
530  const refinementParameters& refineParams,
531  const label nFaces,
532  const label maxIter
533 )
534 {
535  const fvMesh& mesh = meshRefiner_.mesh();
536 
537  label iter;
538  for (iter = 0; iter < maxIter; iter++)
539  {
540  Info<< nl
541  << "Dangling coarse cells refinement iteration " << iter << nl
542  << "--------------------------------------------" << nl
543  << endl;
544 
545 
546  // Determine cells to refine
547  // ~~~~~~~~~~~~~~~~~~~~~~~~~
548 
549  const cellList& cells = mesh.cells();
550  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
551 
552  labelList candidateCells;
553  {
554  cellSet candidateCellSet(mesh, "candidateCells", cells.size()/1000);
555 
556  forAll(cells, celli)
557  {
558  const cell& cFaces = cells[celli];
559 
560  label nIntFaces = 0;
561  forAll(cFaces, i)
562  {
563  label bFacei = cFaces[i]-mesh.nInternalFaces();
564  if (bFacei < 0)
565  {
566  nIntFaces++;
567  }
568  else
569  {
570  label patchi = pbm.patchID()[bFacei];
571  if (pbm[patchi].coupled())
572  {
573  nIntFaces++;
574  }
575  }
576  }
577 
578  if (nIntFaces == nFaces)
579  {
580  candidateCellSet.insert(celli);
581  }
582  }
583 
584  if (debug&meshRefinement::MESH)
585  {
586  Pout<< "Dumping " << candidateCellSet.size()
587  << " cells to cellSet candidateCellSet." << endl;
588  candidateCellSet.instance() = meshRefiner_.name();
589  candidateCellSet.write();
590  }
591  candidateCells = candidateCellSet.toc();
592  }
593 
594 
595 
596  labelList cellsToRefine
597  (
598  meshRefiner_.meshCutter().consistentRefinement
599  (
600  candidateCells,
601  true
602  )
603  );
604  Info<< "Determined cells to refine in = "
605  << mesh.time().cpuTimeIncrement() << " s" << endl;
606 
607 
608  label nCellsToRefine = cellsToRefine.size();
609  reduce(nCellsToRefine, sumOp<label>());
610 
611  Info<< "Selected for refinement : " << nCellsToRefine
612  << " cells (out of " << mesh.globalData().nTotalCells()
613  << ')' << endl;
614 
615  // Stop when no cells to refine. After a few iterations check if too
616  // few cells
617  if
618  (
619  nCellsToRefine == 0
620  || (
621  iter >= 1
622  && nCellsToRefine <= refineParams.minRefineCells()
623  )
624  )
625  {
626  Info<< "Stopping refining since too few cells selected."
627  << nl << endl;
628  break;
629  }
630 
631 
632  if (debug)
633  {
634  const_cast<Time&>(mesh.time())++;
635  }
636 
637 
638  if
639  (
641  (
642  (mesh.nCells() >= refineParams.maxLocalCells()),
643  orOp<bool>()
644  )
645  )
646  {
647  meshRefiner_.balanceAndRefine
648  (
649  "coarse cell refinement iteration " + name(iter),
650  decomposer_,
651  distributor_,
652  cellsToRefine,
653  refineParams.maxLoadUnbalance()
654  );
655  }
656  else
657  {
658  meshRefiner_.refineAndBalance
659  (
660  "coarse cell refinement iteration " + name(iter),
661  decomposer_,
662  distributor_,
663  cellsToRefine,
664  refineParams.maxLoadUnbalance()
665  );
666  }
667  }
668  return iter;
669 }
670 
671 
672 void Foam::snappyRefineDriver::removeInsideCells
673 (
674  const refinementParameters& refineParams,
675  const label nBufferLayers
676 )
677 {
678  Info<< nl
679  << "Removing mesh beyond surface intersections" << nl
680  << "------------------------------------------" << nl
681  << endl;
682 
683  const fvMesh& mesh = meshRefiner_.mesh();
684 
685  if (debug)
686  {
687  const_cast<Time&>(mesh.time())++;
688  }
689 
690  meshRefiner_.splitMesh
691  (
692  nBufferLayers, // nBufferLayers
693  globalToMasterPatch_,
694  globalToSlavePatch_,
695  refineParams.selectionPoints()
696  );
697 
698  if (debug&meshRefinement::MESH)
699  {
700  Pout<< "Writing subsetted mesh to time "
701  << meshRefiner_.name() << '.' << endl;
702  meshRefiner_.write
703  (
706  (
709  ),
710  mesh.time().path()/meshRefiner_.name()
711  );
712  Pout<< "Dumped mesh in = "
713  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
714  }
715 }
716 
717 
718 Foam::label Foam::snappyRefineDriver::shellRefine
719 (
720  const refinementParameters& refineParams,
721  const label maxIter
722 )
723 {
724  const fvMesh& mesh = meshRefiner_.mesh();
725 
726  // Mark current boundary faces with 0. Have meshRefiner maintain them.
727  meshRefiner_.userFaceData().setSize(1);
728 
729  // mark list to remove any refined faces
730  meshRefiner_.userFaceData()[0].first() = meshRefinement::REMOVE;
731  meshRefiner_.userFaceData()[0].second() = createWithValues<labelList>
732  (
733  mesh.nFaces(),
734  -1,
735  meshRefiner_.intersectedFaces(),
736  0
737  );
738 
739  // Determine the maximum refinement level over all volume refinement
740  // regions. This determines the minimum number of shell refinement
741  // iterations.
742  label overallMaxShellLevel = meshRefiner_.shells().maxLevel();
743 
744  label iter;
745  for (iter = 0; iter < maxIter; iter++)
746  {
747  Info<< nl
748  << "Shell refinement iteration " << iter << nl
749  << "----------------------------" << nl
750  << endl;
751 
752  labelList candidateCells
753  (
754  meshRefiner_.refineCandidates
755  (
756  refineParams.selectionPoints().inside(),
757  refineParams.curvature(),
758  refineParams.planarAngle(),
759 
760  false, // featureRefinement
761  true, // featureDistanceRefinement
762  true, // internalRefinement
763  false, // surfaceRefinement
764  false, // curvatureRefinement
765  false, // gapRefinement
766  refineParams.maxGlobalCells(),
767  refineParams.maxLocalCells()
768  )
769  );
770 
771  if (debug&meshRefinement::MESH)
772  {
773  Pout<< "Dumping " << candidateCells.size()
774  << " cells to cellSet candidateCellsFromShells." << endl;
775 
776  cellSet c(mesh, "candidateCellsFromShells", candidateCells);
777  c.instance() = meshRefiner_.name();
778  c.write();
779  }
780 
781  // Problem choosing starting faces for bufferlayers (bFaces)
782  // - we can't use the current intersected boundary faces
783  // (intersectedFaces) since this grows indefinitely
784  // - if we use 0 faces we don't satisfy bufferLayers from the
785  // surface.
786  // - possibly we want to have bFaces only the initial set of faces
787  // and maintain the list while doing the refinement.
788  labelList bFaces
789  (
790  findIndices(meshRefiner_.userFaceData()[0].second(), 0)
791  );
792 
793  // Info<< "Collected boundary faces : "
794  // << returnReduce(bFaces.size(), sumOp<label>()) << endl;
795 
796  labelList cellsToRefine;
797 
798  if (refineParams.nBufferLayers() <= 2)
799  {
800  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement
801  (
802  refineParams.nBufferLayers(),
803  candidateCells, // cells to refine
804  bFaces, // faces for nBufferLayers
805  1, // point difference
806  meshRefiner_.intersectedPoints() // points to check
807  );
808  }
809  else
810  {
811  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement2
812  (
813  refineParams.nBufferLayers(),
814  candidateCells, // cells to refine
815  bFaces // faces for nBufferLayers
816  );
817  }
818 
819  Info<< "Determined cells to refine in = "
820  << mesh.time().cpuTimeIncrement() << " s" << endl;
821 
822 
823  label nCellsToRefine = cellsToRefine.size();
824  reduce(nCellsToRefine, sumOp<label>());
825 
826  Info<< "Selected for internal refinement : " << nCellsToRefine
827  << " cells (out of " << mesh.globalData().nTotalCells()
828  << ')' << endl;
829 
830  // Stop when no cells to refine or have done minimum necessary
831  // iterations and not enough cells to refine.
832  if
833  (
834  nCellsToRefine == 0
835  || (
836  iter >= overallMaxShellLevel
837  && nCellsToRefine <= refineParams.minRefineCells()
838  )
839  )
840  {
841  Info<< "Stopping refining since too few cells selected."
842  << nl << endl;
843  break;
844  }
845 
846 
847  if (debug)
848  {
849  const_cast<Time&>(mesh.time())++;
850  }
851 
852  if
853  (
855  (
856  (mesh.nCells() >= refineParams.maxLocalCells()),
857  orOp<bool>()
858  )
859  )
860  {
861  meshRefiner_.balanceAndRefine
862  (
863  "shell refinement iteration " + name(iter),
864  decomposer_,
865  distributor_,
866  cellsToRefine,
867  refineParams.maxLoadUnbalance()
868  );
869  }
870  else
871  {
872  meshRefiner_.refineAndBalance
873  (
874  "shell refinement iteration " + name(iter),
875  decomposer_,
876  distributor_,
877  cellsToRefine,
878  refineParams.maxLoadUnbalance()
879  );
880  }
881  }
882  meshRefiner_.userFaceData().clear();
883 
884  return iter;
885 }
886 
887 
888 void Foam::snappyRefineDriver::baffleAndSplitMesh
889 (
890  const refinementParameters& refineParams,
891  const snapParameters& snapParams,
892  const bool handleSnapProblems,
893  const dictionary& motionDict
894 )
895 {
896  Info<< nl
897  << "Splitting mesh at surface intersections" << nl
898  << "---------------------------------------" << nl
899  << endl;
900 
901  const fvMesh& mesh = meshRefiner_.mesh();
902 
903  // Introduce baffles at surface intersections. Note:
904  // meshRefiment::surfaceIndex() will
905  // be like boundary face from now on so not coupled anymore.
906  meshRefiner_.baffleAndSplitMesh
907  (
908  handleSnapProblems, // detect&remove potential snap problem
909 
910  // Snap problem cell detection
911  snapParams,
912  refineParams.useTopologicalSnapDetection(),
913  false, // perpendicular edge connected cells
914  scalarField(0), // per region perpendicular angle
915 
916  // Free standing baffles
917  !handleSnapProblems, // merge free standing baffles?
918  refineParams.planarAngle(),
919 
920  motionDict,
921  const_cast<Time&>(mesh.time()),
922  globalToMasterPatch_,
923  globalToSlavePatch_,
924  refineParams.selectionPoints()
925  );
926 }
927 
928 
929 void Foam::snappyRefineDriver::zonify
930 (
931  const refinementParameters& refineParams
932 )
933 {
934  // Mesh is at its finest. Do zoning
935  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
936  // This puts all faces with intersection across a zoneable surface
937  // into that surface's faceZone. All cells inside faceZone get given the
938  // same cellZone.
939 
940  const labelList namedSurfaces =
941  surfaceZonesInfo::getNamedSurfaces(meshRefiner_.surfaces().surfZones());
942 
943  if (namedSurfaces.size())
944  {
945  Info<< nl
946  << "Introducing zones for interfaces" << nl
947  << "--------------------------------" << nl
948  << endl;
949 
950  const fvMesh& mesh = meshRefiner_.mesh();
951 
952  if (debug)
953  {
954  const_cast<Time&>(mesh.time())++;
955  }
956 
957  meshRefiner_.zonify
958  (
959  refineParams.selectionPoints().inside(),
960  refineParams.allowFreeStandingZoneFaces()
961  );
962 
963  if (debug&meshRefinement::MESH)
964  {
965  Pout<< "Writing zoned mesh to time "
966  << meshRefiner_.name() << '.' << endl;
967  meshRefiner_.write
968  (
971  (
974  ),
975  mesh.time().path()/meshRefiner_.name()
976  );
977  }
978 
979  // Check that all faces are synced
981  }
982 }
983 
984 
985 void Foam::snappyRefineDriver::splitAndMergeBaffles
986 (
987  const refinementParameters& refineParams,
988  const snapParameters& snapParams,
989  const bool handleSnapProblems,
990  const dictionary& motionDict
991 )
992 {
993  Info<< nl
994  << "Handling cells with snap problems" << nl
995  << "---------------------------------" << nl
996  << endl;
997 
998  const fvMesh& mesh = meshRefiner_.mesh();
999 
1000  // Introduce baffles and split mesh
1001  if (debug)
1002  {
1003  const_cast<Time&>(mesh.time())++;
1004  }
1005 
1006  const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
1007 
1008  meshRefiner_.baffleAndSplitMesh
1009  (
1010  handleSnapProblems,
1011 
1012  // Snap problem cell detection
1013  snapParams,
1014  refineParams.useTopologicalSnapDetection(),
1015  handleSnapProblems, // remove perp edge connected cells
1016  perpAngle, // perp angle
1017 
1018  // Free standing baffles
1019  true, // merge free standing baffles?
1020  refineParams.planarAngle(), // planar angle
1021 
1022  motionDict,
1023  const_cast<Time&>(mesh.time()),
1024  globalToMasterPatch_,
1025  globalToSlavePatch_,
1026  refineParams.selectionPoints()
1027  );
1028 
1029  if (debug)
1030  {
1031  const_cast<Time&>(mesh.time())++;
1032  }
1033 
1034  // Duplicate points on baffles that are on more than one cell
1035  // region. This will help snapping pull them to separate surfaces.
1036  meshRefiner_.dupNonManifoldPoints();
1037 
1038 
1039  // Merge all baffles that are still remaining after duplicating points.
1040  List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
1041 
1042  label nCouples = returnReduce(couples.size(), sumOp<label>());
1043 
1044  Info<< "Detected unsplittable baffles : " << nCouples << endl;
1045 
1046  if (nCouples > 0)
1047  {
1048  // Actually merge baffles. Note: not exactly parallelised. Should
1049  // convert baffle faces into processor faces if they resulted
1050  // from them.
1051  meshRefiner_.mergeBaffles(couples);
1052 
1053  if (debug)
1054  {
1055  // Debug:test all is still synced across proc patches
1056  meshRefiner_.checkData();
1057  }
1058 
1059  // Remove any now dangling parts
1060  meshRefiner_.splitMeshRegions
1061  (
1062  globalToMasterPatch_,
1063  globalToSlavePatch_,
1064  refineParams.selectionPoints()
1065  );
1066 
1067  if (debug)
1068  {
1069  // Debug:test all is still synced across proc patches
1070  meshRefiner_.checkData();
1071  }
1072 
1073  Info<< "Merged free-standing baffles in = "
1074  << mesh.time().cpuTimeIncrement() << " s." << endl;
1075  }
1076 
1077  if (debug&meshRefinement::MESH)
1078  {
1079  Pout<< "Writing handleProblemCells mesh to time "
1080  << meshRefiner_.name() << '.' << endl;
1081  meshRefiner_.write
1082  (
1085  (
1088  ),
1089  mesh.time().path()/meshRefiner_.name()
1090  );
1091  }
1092 }
1093 
1094 
1095 void Foam::snappyRefineDriver::mergePatchFaces
1096 (
1097  const refinementParameters& refineParams,
1098  const dictionary& motionDict
1099 )
1100 {
1101  Info<< nl
1102  << "Merge refined boundary faces" << nl
1103  << "----------------------------" << nl
1104  << endl;
1105 
1106  const fvMesh& mesh = meshRefiner_.mesh();
1107 
1108  meshRefiner_.mergePatchFacesUndo
1109  (
1110  Foam::cos(degToRad(45.0)),
1111  Foam::cos(degToRad(45.0)),
1112  meshRefiner_.meshedPatches(),
1113  motionDict,
1114  labelList(mesh.nFaces(), -1)
1115  );
1116 
1117  if (debug)
1118  {
1119  meshRefiner_.checkData();
1120  }
1121 
1122  meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
1123 
1124  if (debug)
1125  {
1126  meshRefiner_.checkData();
1127  }
1128 }
1129 
1130 
1132 (
1133  const dictionary& refineDict,
1134  const refinementParameters& refineParams,
1135  const snapParameters& snapParams,
1136  const bool prepareForSnapping,
1137  const dictionary& motionDict
1138 )
1139 {
1140  Info<< nl
1141  << "Refinement phase" << nl
1142  << "----------------" << nl
1143  << endl;
1144 
1145  const fvMesh& mesh = meshRefiner_.mesh();
1146 
1147  // Check that all the keep points are inside the mesh.
1148  refineParams.findCells(mesh);
1149 
1150  // Refine around feature edges
1151  featureEdgeRefine
1152  (
1153  refineParams,
1154  100, // maxIter
1155  0 // min cells to refine
1156  );
1157 
1158  // Refine based on surface
1159  surfaceOnlyRefine
1160  (
1161  refineParams,
1162  100 // maxIter
1163  );
1164 
1165  gapOnlyRefine
1166  (
1167  refineParams,
1168  100 // maxIter
1169  );
1170 
1171  // Remove cells (a certain distance) beyond surface intersections
1172  removeInsideCells
1173  (
1174  refineParams,
1175  1 // nBufferLayers
1176  );
1177 
1178  // Internal mesh refinement
1179  shellRefine
1180  (
1181  refineParams,
1182  100 // maxIter
1183  );
1184 
1185  // Refine any hexes with 5 or 6 faces refined to make smooth edges
1186  danglingCellRefine
1187  (
1188  refineParams,
1189  21, // 1 coarse face + 5 refined faces
1190  100 // maxIter
1191  );
1192  danglingCellRefine
1193  (
1194  refineParams,
1195  24, // 0 coarse faces + 6 refined faces
1196  100 // maxIter
1197  );
1198 
1199  // Introduce baffles at surface intersections. Remove sections unreachable
1200  // from insidePoints.
1201  baffleAndSplitMesh
1202  (
1203  refineParams,
1204  snapParams,
1205  prepareForSnapping,
1206  motionDict
1207  );
1208 
1209  // Mesh is at its finest. Do optional zoning.
1210  zonify(refineParams);
1211 
1212  // Pull baffles apart
1213  splitAndMergeBaffles
1214  (
1215  refineParams,
1216  snapParams,
1217  prepareForSnapping,
1218  motionDict
1219  );
1220 
1221  // Do something about cells with refined faces on the boundary
1222  if (prepareForSnapping)
1223  {
1224  mergePatchFaces(refineParams, motionDict);
1225  }
1226 
1227 
1228  if (Pstream::parRun())
1229  {
1230  Info<< nl
1231  << "Doing final balancing" << nl
1232  << "---------------------" << nl
1233  << endl;
1234 
1235  // if (debug)
1236  //{
1237  // const_cast<Time&>(mesh.time())++;
1238  //}
1239 
1240  // Do final balancing. Keep zoned faces on one processor since the
1241  // snap phase will convert them to baffles and this only works for
1242  // internal faces.
1243  meshRefiner_.balance
1244  (
1245  true,
1246  false,
1247  scalarField(mesh.nCells(), 1), // dummy weights
1248  decomposer_,
1249  distributor_
1250  );
1251 
1252 
1253  if (debug)
1254  {
1255  meshRefiner_.checkZoneFaces();
1256  }
1257  }
1258 }
1259 
1260 
1261 // ************************************************************************* //
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.H:294
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:85
virtual Ostream & write(const char)=0
Write character.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:60
Abstract base class for decomposition.
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
const word & name() const
Return const reference to name.
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:406
label nTotalCells() const
Return total number of cells in decomposed mesh.
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
@ REMOVE
set value to -1 any face that was refined
static void checkCoupledFaceZones(const polyMesh &)
Helper function: check that face zones are synced.
static writeType writeLevel()
Get/set write level.
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1563
label nCells() const
const List< point > & inside() const
Return the points inside the surface regions to selected cells.
Simple container to keep together refinement specific information.
scalar planarAngle() const
Angle when two intersections are considered to be planar.
label maxLocalCells() const
Per processor max number of cells.
scalar curvature() const
Curvature.
scalar maxLoadUnbalance() const
Allowed load unbalance.
label maxGlobalCells() const
Total number of cells.
const cellSelectionPoints & selectionPoints() const
Return the points to select cells inside and outside.
labelList findCells(const polyMesh &) const
Checks that cells are in mesh. Returns cells they are in.
Simple container to keep together snap specific information.
void doRefine(const dictionary &refineDict, const refinementParameters &refineParams, const snapParameters &snapParams, const bool prepareForSnapping, const dictionary &motionDict)
Do all the refinement.
snappyRefineDriver(meshRefinement &meshRefiner, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch)
Construct from components.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
label patchi
const cellShapeList & cells
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: labelList.H:56
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
List< cell > cellList
list of cells
Definition: cellList.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
messageStream Info
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
static const char nl
Definition: Ostream.H:260
dimensionedScalar cos(const dimensionedScalar &ds)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Unit conversion functions.