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