FreeFOAM The Cross-Platform CFD Toolkit
autoRefineDriver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "autoRefineDriver.H"
28 #include <finiteVolume/fvMesh.H>
29 #include <OpenFOAM/Time.H>
30 #include <meshTools/cellSet.H>
31 #include <OpenFOAM/syncTools.H>
35 #include <autoMesh/shellSurfaces.H>
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43 
44 defineTypeNameAndDebug(autoRefineDriver, 0);
45 
46 } // End namespace Foam
47 
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 // Read explicit feature edges
52 Foam::label Foam::autoRefineDriver::readFeatureEdges
53 (
54  const PtrList<dictionary>& featDicts,
55  PtrList<featureEdgeMesh>& featureMeshes,
56  labelList& featureLevels
57 ) const
58 {
59  Info<< "Reading external feature lines." << endl;
60 
61  const fvMesh& mesh = meshRefiner_.mesh();
62 
63  featureMeshes.setSize(featDicts.size());
64  featureLevels.setSize(featDicts.size());
65 
66  forAll(featDicts, i)
67  {
68  const dictionary& dict = featDicts[i];
69 
70  fileName featFileName(dict.lookup("file"));
71 
72  featureMeshes.set
73  (
74  i,
75  new featureEdgeMesh
76  (
77  IOobject
78  (
79  featFileName, // name
80  //mesh.time().findInstance("triSurface", featFileName),
81  // // instance
82  mesh.time().constant(), // instance
83  "triSurface", // local
84  mesh.time(), // registry
87  false
88  )
89  )
90  );
91 
92  featureMeshes[i].mergePoints(meshRefiner_.mergeDistance());
93  featureLevels[i] = readLabel(dict.lookup("level"));
94 
95  Info<< "Refinement level " << featureLevels[i]
96  << " for all cells crossed by feature " << featFileName
97  << " (" << featureMeshes[i].points().size() << " points, "
98  << featureMeshes[i].edges().size() << " edges)." << endl;
99  }
100 
101  Info<< "Read feature lines in = "
102  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
103 
104  return featureMeshes.size();
105 }
106 
107 
108 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
109 
110 // Construct from components
111 Foam::autoRefineDriver::autoRefineDriver
112 (
113  meshRefinement& meshRefiner,
114  decompositionMethod& decomposer,
115  fvMeshDistribute& distributor,
116  const labelList& globalToPatch
117 )
118 :
119  meshRefiner_(meshRefiner),
120  decomposer_(decomposer),
121  distributor_(distributor),
122  globalToPatch_(globalToPatch)
123 {}
124 
125 
126 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
127 
128 Foam::label Foam::autoRefineDriver::featureEdgeRefine
129 (
130  const refinementParameters& refineParams,
131  const PtrList<dictionary>& featDicts,
132  const label maxIter,
133  const label minRefine
134 )
135 {
136  const fvMesh& mesh = meshRefiner_.mesh();
137 
138  // Read explicit feature edges
139  PtrList<featureEdgeMesh> featureMeshes;
140  // Per feature the refinement level
141  labelList featureLevels;
142  readFeatureEdges(featDicts, featureMeshes, featureLevels);
143 
144 
145  label iter = 0;
146 
147  if (featureMeshes.size() && maxIter > 0)
148  {
149  for (; iter < maxIter; iter++)
150  {
151  Info<< nl
152  << "Feature refinement iteration " << iter << nl
153  << "------------------------------" << nl
154  << endl;
155 
156  labelList candidateCells
157  (
158  meshRefiner_.refineCandidates
159  (
160  refineParams.keepPoints()[0], // For now only use one.
161  refineParams.curvature(),
162 
163  featureMeshes,
164  featureLevels,
165 
166  true, // featureRefinement
167  false, // internalRefinement
168  false, // surfaceRefinement
169  false, // curvatureRefinement
170  refineParams.maxGlobalCells(),
171  refineParams.maxLocalCells()
172  )
173  );
174  labelList cellsToRefine
175  (
176  meshRefiner_.meshCutter().consistentRefinement
177  (
178  candidateCells,
179  true
180  )
181  );
182  Info<< "Determined cells to refine in = "
183  << mesh.time().cpuTimeIncrement() << " s" << endl;
184 
185 
186 
187  label nCellsToRefine = cellsToRefine.size();
188  reduce(nCellsToRefine, sumOp<label>());
189 
190  Info<< "Selected for feature refinement : " << nCellsToRefine
191  << " cells (out of " << mesh.globalData().nTotalCells()
192  << ')' << endl;
193 
194  if (nCellsToRefine <= minRefine)
195  {
196  Info<< "Stopping refining since too few cells selected."
197  << nl << endl;
198  break;
199  }
200 
201 
202  if (debug > 0)
203  {
204  const_cast<Time&>(mesh.time())++;
205  }
206 
207 
208  if
209  (
211  (
212  (mesh.nCells() >= refineParams.maxLocalCells()),
213  orOp<bool>()
214  )
215  )
216  {
217  meshRefiner_.balanceAndRefine
218  (
219  "feature refinement iteration " + name(iter),
220  decomposer_,
221  distributor_,
222  cellsToRefine,
223  refineParams.maxLoadUnbalance()
224  );
225  }
226  else
227  {
228  meshRefiner_.refineAndBalance
229  (
230  "feature refinement iteration " + name(iter),
231  decomposer_,
232  distributor_,
233  cellsToRefine,
234  refineParams.maxLoadUnbalance()
235  );
236  }
237  }
238  }
239  return iter;
240 }
241 
242 
243 Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
244 (
245  const refinementParameters& refineParams,
246  const label maxIter
247 )
248 {
249  const fvMesh& mesh = meshRefiner_.mesh();
250 
251  // Determine the maximum refinement level over all surfaces. This
252  // determines the minumum number of surface refinement iterations.
253  label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
254 
255  label iter;
256  for (iter = 0; iter < maxIter; iter++)
257  {
258  Info<< nl
259  << "Surface refinement iteration " << iter << nl
260  << "------------------------------" << nl
261  << endl;
262 
263 
264  // Determine cells to refine
265  // ~~~~~~~~~~~~~~~~~~~~~~~~~
266  // Only look at surface intersections (minLevel and surface curvature),
267  // do not do internal refinement (refinementShells)
268 
269  const PtrList<featureEdgeMesh> dummyFeatures;
270 
271  labelList candidateCells
272  (
273  meshRefiner_.refineCandidates
274  (
275  refineParams.keepPoints()[0],
276  refineParams.curvature(),
277 
278  dummyFeatures, // dummy featureMeshes;
279  labelList(0), // dummy featureLevels;
280 
281  false, // featureRefinement
282  false, // internalRefinement
283  true, // surfaceRefinement
284  true, // curvatureRefinement
285  refineParams.maxGlobalCells(),
286  refineParams.maxLocalCells()
287  )
288  );
289  labelList cellsToRefine
290  (
291  meshRefiner_.meshCutter().consistentRefinement
292  (
293  candidateCells,
294  true
295  )
296  );
297  Info<< "Determined cells to refine in = "
298  << mesh.time().cpuTimeIncrement() << " s" << endl;
299 
300 
301  label nCellsToRefine = cellsToRefine.size();
302  reduce(nCellsToRefine, sumOp<label>());
303 
304  Info<< "Selected for refinement : " << nCellsToRefine
305  << " cells (out of " << mesh.globalData().nTotalCells()
306  << ')' << endl;
307 
308  // Stop when no cells to refine or have done minimum nessecary
309  // iterations and not enough cells to refine.
310  if
311  (
312  nCellsToRefine == 0
313  || (
314  iter >= overallMaxLevel
315  && nCellsToRefine <= refineParams.minRefineCells()
316  )
317  )
318  {
319  Info<< "Stopping refining since too few cells selected."
320  << nl << endl;
321  break;
322  }
323 
324 
325  if (debug)
326  {
327  const_cast<Time&>(mesh.time())++;
328  }
329 
330 
331  if
332  (
334  (
335  (mesh.nCells() >= refineParams.maxLocalCells()),
336  orOp<bool>()
337  )
338  )
339  {
340  meshRefiner_.balanceAndRefine
341  (
342  "surface refinement iteration " + name(iter),
343  decomposer_,
344  distributor_,
345  cellsToRefine,
346  refineParams.maxLoadUnbalance()
347  );
348  }
349  else
350  {
351  meshRefiner_.refineAndBalance
352  (
353  "surface refinement iteration " + name(iter),
354  decomposer_,
355  distributor_,
356  cellsToRefine,
357  refineParams.maxLoadUnbalance()
358  );
359  }
360  }
361  return iter;
362 }
363 
364 
365 void Foam::autoRefineDriver::removeInsideCells
366 (
367  const refinementParameters& refineParams,
368  const label nBufferLayers
369 )
370 {
371  Info<< nl
372  << "Removing mesh beyond surface intersections" << nl
373  << "------------------------------------------" << nl
374  << endl;
375 
376  const fvMesh& mesh = meshRefiner_.mesh();
377 
378  if (debug)
379  {
380  const_cast<Time&>(mesh.time())++;
381  }
382 
383  meshRefiner_.splitMesh
384  (
385  nBufferLayers, // nBufferLayers
386  globalToPatch_,
387  refineParams.keepPoints()[0]
388  );
389 
390  if (debug)
391  {
392  Pout<< "Writing subsetted mesh to time "
393  << meshRefiner_.timeName() << '.' << endl;
394  meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
395  Pout<< "Dumped mesh in = "
396  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
397  }
398 }
399 
400 
401 Foam::label Foam::autoRefineDriver::shellRefine
402 (
403  const refinementParameters& refineParams,
404  const label maxIter
405 )
406 {
407  const fvMesh& mesh = meshRefiner_.mesh();
408 
409  // Mark current boundary faces with 0. Have meshRefiner maintain them.
410  meshRefiner_.userFaceData().setSize(1);
411 
412  // mark list to remove any refined faces
413  meshRefiner_.userFaceData()[0].first() = meshRefinement::REMOVE;
414  meshRefiner_.userFaceData()[0].second() = createWithValues<labelList>
415  (
416  mesh.nFaces(),
417  -1,
418  meshRefiner_.intersectedFaces(),
419  0
420  );
421 
422  // Determine the maximum refinement level over all volume refinement
423  // regions. This determines the minumum number of shell refinement
424  // iterations.
425  label overallMaxShellLevel = meshRefiner_.shells().maxLevel();
426 
427  label iter;
428  for (iter = 0; iter < maxIter; iter++)
429  {
430  Info<< nl
431  << "Shell refinement iteration " << iter << nl
432  << "----------------------------" << nl
433  << endl;
434 
435  const PtrList<featureEdgeMesh> dummyFeatures;
436 
437  labelList candidateCells
438  (
439  meshRefiner_.refineCandidates
440  (
441  refineParams.keepPoints()[0],
442  refineParams.curvature(),
443 
444  dummyFeatures, // dummy featureMeshes;
445  labelList(0), // dummy featureLevels;
446 
447  false, // featureRefinement
448  true, // internalRefinement
449  false, // surfaceRefinement
450  false, // curvatureRefinement
451  refineParams.maxGlobalCells(),
452  refineParams.maxLocalCells()
453  )
454  );
455 
456  if (debug)
457  {
458  Pout<< "Dumping " << candidateCells.size()
459  << " cells to cellSet candidateCellsFromShells." << endl;
460 
461  cellSet(mesh, "candidateCellsFromShells", candidateCells).write();
462  }
463 
464  // Problem choosing starting faces for bufferlayers (bFaces)
465  // - we can't use the current intersected boundary faces
466  // (intersectedFaces) since this grows indefinitely
467  // - if we use 0 faces we don't satisfy bufferLayers from the
468  // surface.
469  // - possibly we want to have bFaces only the initial set of faces
470  // and maintain the list while doing the refinement.
471  labelList bFaces
472  (
473  findIndices(meshRefiner_.userFaceData()[0].second(), 0)
474  );
475 
476  //Info<< "Collected boundary faces : "
477  // << returnReduce(bFaces.size(), sumOp<label>()) << endl;
478 
479  labelList cellsToRefine;
480 
481  if (refineParams.nBufferLayers() <= 2)
482  {
483  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement
484  (
485  refineParams.nBufferLayers(),
486  candidateCells, // cells to refine
487  bFaces, // faces for nBufferLayers
488  1, // point difference
489  meshRefiner_.intersectedPoints() // points to check
490  );
491  }
492  else
493  {
494  cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement2
495  (
496  refineParams.nBufferLayers(),
497  candidateCells, // cells to refine
498  bFaces // faces for nBufferLayers
499  );
500  }
501 
502  Info<< "Determined cells to refine in = "
503  << mesh.time().cpuTimeIncrement() << " s" << endl;
504 
505 
506  label nCellsToRefine = cellsToRefine.size();
507  reduce(nCellsToRefine, sumOp<label>());
508 
509  Info<< "Selected for internal refinement : " << nCellsToRefine
510  << " cells (out of " << mesh.globalData().nTotalCells()
511  << ')' << endl;
512 
513  // Stop when no cells to refine or have done minimum nessecary
514  // iterations and not enough cells to refine.
515  if
516  (
517  nCellsToRefine == 0
518  || (
519  iter >= overallMaxShellLevel
520  && nCellsToRefine <= refineParams.minRefineCells()
521  )
522  )
523  {
524  Info<< "Stopping refining since too few cells selected."
525  << nl << endl;
526  break;
527  }
528 
529 
530  if (debug)
531  {
532  const_cast<Time&>(mesh.time())++;
533  }
534 
535  if
536  (
538  (
539  (mesh.nCells() >= refineParams.maxLocalCells()),
540  orOp<bool>()
541  )
542  )
543  {
544  meshRefiner_.balanceAndRefine
545  (
546  "shell refinement iteration " + name(iter),
547  decomposer_,
548  distributor_,
549  cellsToRefine,
550  refineParams.maxLoadUnbalance()
551  );
552  }
553  else
554  {
555  meshRefiner_.refineAndBalance
556  (
557  "shell refinement iteration " + name(iter),
558  decomposer_,
559  distributor_,
560  cellsToRefine,
561  refineParams.maxLoadUnbalance()
562  );
563  }
564  }
565  meshRefiner_.userFaceData().clear();
566 
567  return iter;
568 }
569 
570 
571 void Foam::autoRefineDriver::baffleAndSplitMesh
572 (
573  const refinementParameters& refineParams,
574  const bool handleSnapProblems,
575  const dictionary& motionDict
576 )
577 {
578  Info<< nl
579  << "Splitting mesh at surface intersections" << nl
580  << "---------------------------------------" << nl
581  << endl;
582 
583  const fvMesh& mesh = meshRefiner_.mesh();
584 
585  // Introduce baffles at surface intersections. Note:
586  // meshRefiment::surfaceIndex() will
587  // be like boundary face from now on so not coupled anymore.
588  meshRefiner_.baffleAndSplitMesh
589  (
590  handleSnapProblems, // detect&remove potential snap problem
591  false, // perpendicular edge connected cells
592  scalarField(0), // per region perpendicular angle
593  !handleSnapProblems, // merge free standing baffles?
594  motionDict,
595  const_cast<Time&>(mesh.time()),
596  globalToPatch_,
597  refineParams.keepPoints()[0]
598  );
599 }
600 
601 
602 void Foam::autoRefineDriver::zonify
603 (
604  const refinementParameters& refineParams
605 )
606 {
607  // Mesh is at its finest. Do zoning
608  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
609  // This puts all faces with intersection across a zoneable surface
610  // into that surface's faceZone. All cells inside faceZone get given the
611  // same cellZone.
612 
613  if (meshRefiner_.surfaces().getNamedSurfaces().size())
614  {
615  Info<< nl
616  << "Introducing zones for interfaces" << nl
617  << "--------------------------------" << nl
618  << endl;
619 
620  const fvMesh& mesh = meshRefiner_.mesh();
621 
622  if (debug)
623  {
624  const_cast<Time&>(mesh.time())++;
625  }
626 
627  meshRefiner_.zonify
628  (
629  refineParams.keepPoints()[0],
630  refineParams.allowFreeStandingZoneFaces()
631  );
632 
633  if (debug)
634  {
635  Pout<< "Writing zoned mesh to time "
636  << meshRefiner_.timeName() << '.' << endl;
637  meshRefiner_.write
638  (
639  debug,
640  mesh.time().path()/meshRefiner_.timeName()
641  );
642  }
643 
644  // Check that all faces are synced
646  }
647 }
648 
649 
650 void Foam::autoRefineDriver::splitAndMergeBaffles
651 (
652  const refinementParameters& refineParams,
653  const bool handleSnapProblems,
654  const dictionary& motionDict
655 )
656 {
657  Info<< nl
658  << "Handling cells with snap problems" << nl
659  << "---------------------------------" << nl
660  << endl;
661 
662  const fvMesh& mesh = meshRefiner_.mesh();
663 
664  // Introduce baffles and split mesh
665  if (debug)
666  {
667  const_cast<Time&>(mesh.time())++;
668  }
669 
670  const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
671 
672  meshRefiner_.baffleAndSplitMesh
673  (
674  handleSnapProblems,
675  handleSnapProblems, // remove perp edge connected cells
676  perpAngle, // perp angle
677  false, // merge free standing baffles?
678  motionDict,
679  const_cast<Time&>(mesh.time()),
680  globalToPatch_,
681  refineParams.keepPoints()[0]
682  );
683 
684  if (debug)
685  {
686  const_cast<Time&>(mesh.time())++;
687  }
688 
689  // Duplicate points on baffles that are on more than one cell
690  // region. This will help snapping pull them to separate surfaces.
691  meshRefiner_.dupNonManifoldPoints();
692 
693 
694  // Merge all baffles that are still remaining after duplicating points.
695  List<labelPair> couples
696  (
697  meshRefiner_.getDuplicateFaces // get all baffles
698  (
699  identity(mesh.nFaces()-mesh.nInternalFaces())
700  + mesh.nInternalFaces()
701  )
702  );
703 
704  label nCouples = returnReduce(couples.size(), sumOp<label>());
705 
706  Info<< "Detected unsplittable baffles : "
707  << nCouples << endl;
708 
709  if (nCouples > 0)
710  {
711  // Actually merge baffles. Note: not exactly parallellized. Should
712  // convert baffle faces into processor faces if they resulted
713  // from them.
714  meshRefiner_.mergeBaffles(couples);
715 
716  if (debug)
717  {
718  // Debug:test all is still synced across proc patches
719  meshRefiner_.checkData();
720  }
721  Info<< "Merged free-standing baffles in = "
722  << mesh.time().cpuTimeIncrement() << " s." << endl;
723  }
724 
725  if (debug)
726  {
727  Pout<< "Writing handleProblemCells mesh to time "
728  << meshRefiner_.timeName() << '.' << endl;
729  meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
730  }
731 }
732 
733 
734 void Foam::autoRefineDriver::mergePatchFaces
735 (
736  const refinementParameters& refineParams
737 )
738 {
739  const fvMesh& mesh = meshRefiner_.mesh();
740 
741  Info<< nl
742  << "Merge refined boundary faces" << nl
743  << "----------------------------" << nl
744  << endl;
745 
746  if (debug)
747  {
748  const_cast<Time&>(mesh.time())++;
749  }
750 
751  meshRefiner_.mergePatchFaces
752  (
755  meshRefiner_.meshedPatches()
756  );
757 
758  if (debug)
759  {
760  meshRefiner_.checkData();
761  }
762 
763  meshRefiner_.mergeEdges(Foam::cos(45*mathematicalConstant::pi/180.0));
764 
765  if (debug)
766  {
767  meshRefiner_.checkData();
768  }
769 }
770 
771 
773 (
774  const dictionary& refineDict,
775  const refinementParameters& refineParams,
776  const bool prepareForSnapping,
777  const dictionary& motionDict
778 )
779 {
780  Info<< nl
781  << "Refinement phase" << nl
782  << "----------------" << nl
783  << endl;
784 
785  const fvMesh& mesh = meshRefiner_.mesh();
786 
787  // Check that all the keep points are inside the mesh.
788  refineParams.findCells(mesh);
789 
790  PtrList<dictionary> featDicts(refineDict.lookup("features"));
791 
792  // Refine around feature edges
793  featureEdgeRefine
794  (
795  refineParams,
796  featDicts,
797  100, // maxIter
798  0 // min cells to refine
799  );
800 
801  // Refine based on surface
802  surfaceOnlyRefine
803  (
804  refineParams,
805  100 // maxIter
806  );
807 
808  // Remove cells (a certain distance) beyond surface intersections
809  removeInsideCells
810  (
811  refineParams,
812  1 // nBufferLayers
813  );
814 
815  // Internal mesh refinement
816  shellRefine
817  (
818  refineParams,
819  100 // maxIter
820  );
821 
822  // Introduce baffles at surface intersections
823  baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
824 
825  // Mesh is at its finest. Do optional zoning.
826  zonify(refineParams);
827 
828  // Pull baffles apart
829  splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
830 
831  // Do something about cells with refined faces on the boundary
832  if (prepareForSnapping)
833  {
834  mergePatchFaces(refineParams);
835  }
836 
837 
838  if (Pstream::parRun())
839  {
840  Info<< nl
841  << "Doing final balancing" << nl
842  << "---------------------" << nl
843  << endl;
844 
845  if (debug)
846  {
847  const_cast<Time&>(mesh.time())++;
848  }
849 
850  // Do final balancing. Keep zoned faces on one processor since the
851  // snap phase will convert them to baffles and this only works for
852  // internal faces.
853  meshRefiner_.balance
854  (
855  true,
856  false,
857  scalarField(mesh.nCells(), 1), // dummy weights
858  decomposer_,
859  distributor_
860  );
861  }
862 }
863 
864 
865 // ************************ vim: set sw=4 sts=4 et: ************************ //