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-2018 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 "shellSurfaces.H"
36 #include "mapDistributePolyMesh.H"
37 #include "unitConversion.H"
38 #include "snapParameters.H"
39 #include "localPointRegion.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45 
46 defineTypeNameAndDebug(snappyRefineDriver, 0);
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.keepPoints(),
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.keepPoints(),
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.keepPoints(),
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_.timeName();
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_.timeName();
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_.timeName();
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.keepPoints()[0]
696  );
697 
698  if (debug&meshRefinement::MESH)
699  {
700  Pout<< "Writing subsetted mesh to time "
701  << meshRefiner_.timeName() << '.' << endl;
702  meshRefiner_.write
703  (
706  (
709  ),
710  mesh.time().path()/meshRefiner_.timeName()
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.keepPoints(),
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_.timeName();
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.keepPoints()[0]
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.keepPoints()[0],
960  refineParams.allowFreeStandingZoneFaces()
961  );
962 
963  if (debug&meshRefinement::MESH)
964  {
965  Pout<< "Writing zoned mesh to time "
966  << meshRefiner_.timeName() << '.' << endl;
967  meshRefiner_.write
968  (
971  (
974  ),
975  mesh.time().path()/meshRefiner_.timeName()
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.keepPoints()[0]
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.
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 parallellized. 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.keepPoints()[0]
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_.timeName() << '.' << endl;
1081  meshRefiner_.write
1082  (
1085  (
1088  ),
1089  mesh.time().path()/meshRefiner_.timeName()
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 
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 keepPoint.
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 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
fileName path() const
Return path.
Definition: Time.H:266
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
label nFaces() const
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const labelList & patchID() const
Per boundary face label the patch index.
label nTotalCells() const
Return total number of cells in decomposed mesh.
const pointField & keepPoints() const
Areas to keep.
static writeType writeLevel()
Get/set write level.
const cellList & cells() const
set value to -1 any face that was refined
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
Simple container to keep together refinement specific information.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
labelList findCells(const polyMesh &) const
Checks that cells are in mesh. Returns cells they are in.
label minRefineCells() const
When to stop refining.
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
dynamicFvMesh & mesh
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
bool allowFreeStandingZoneFaces() const
Are zone faces allowed only in between different cell zones.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Simple container to keep together snap specific information.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:74
void doRefine(const dictionary &refineDict, const refinementParameters &refineParams, const snapParameters &snapParams, const bool prepareForSnapping, const dictionary &motionDict)
Do all the refinement.
bool useTopologicalSnapDetection() const
Use old topology based problem-cell removal.
label maxGlobalCells() const
Total number of cells.
scalar maxLoadUnbalance() const
Allowed load unbalance.
static void checkCoupledFaceZones(const polyMesh &)
Helper function: check that face zones are synced.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1394
List< label > labelList
A List of labels.
Definition: labelList.H:56
Abstract base class for decomposition.
Foam::polyBoundaryMesh.
static const char nl
Definition: Ostream.H:265
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
scalar planarAngle() const
Angle when two intersections are considered to be planar.
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
snappyRefineDriver(meshRefinement &meshRefiner, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch)
Construct from components.
void setSize(const label)
Reset size of List.
Definition: List.C:281
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
const fileName & instance() const
Definition: IOobject.H:378
label patchi
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
A collection of cell labels.
Definition: cellSet.H:48
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const dimensionedScalar c
Speed of light in a vacuum.
void clear()
Remove all regIOobject owned by the registry.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
label n
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual Ostream & write(const token &)=0
Write next token to stream.
virtual bool write(const bool write=true) const
Write using setting from DB.
T * first()
Return the first entry.
Definition: UILList.H:109
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
scalar curvature() const
Curvature.
Namespace for OpenFOAM.
label maxLocalCells() const
Per processor max number of cells.
label nBufferLayers() const
Number of layers between different refinement levels.