FreeFOAM The Cross-Platform CFD Toolkit
autoRefineMesh.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 Application
25  autoRefineMesh
26 
27 Description
28  Utility to refine cells near to a surface.
29 
30 Usage
31 
32  - autoRefineMesh [OPTIONS]
33 
34  @param -case <dir>\n
35  Case directory.
36 
37  @param -help \n
38  Display help message.
39 
40  @param -doc \n
41  Display Doxygen API documentation page for this application.
42 
43  @param -srcDoc \n
44  Display Doxygen source documentation page for this application.
45 
46 \*---------------------------------------------------------------------------*/
47 
48 #include <OpenFOAM/argList.H>
49 #include <OpenFOAM/Time.H>
50 #include <OpenFOAM/polyMesh.H>
52 #include <OpenFOAM/OFstream.H>
54 
55 #include <triSurface/triSurface.H>
57 
58 #include <meshTools/cellSet.H>
59 #include <meshTools/pointSet.H>
62 #include <meshTools/cellToPoint.H>
63 #include <meshTools/pointToCell.H>
64 #include <meshTools/cellToCell.H>
65 #include <meshTools/surfaceSets.H>
68 #include <OpenFOAM/mapPolyMesh.H>
69 #include <OpenFOAM/labelIOList.H>
72 #include <meshTools/meshSearch.H>
73 
74 using namespace Foam;
75 
76 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
77 
78 
79 // Max cos angle for edges to be considered aligned with axis.
80 static const scalar edgeTol = 1E-3;
81 
82 
83 void writeSet(const cellSet& cells, const string& msg)
84 {
85  Info<< "Writing " << msg << " (" << cells.size() << ") to cellSet "
86  << cells.instance()/cells.local()/cells.name()
87  << endl;
88  cells.write();
89 }
90 
91 
92 direction getNormalDir(const twoDPointCorrector* correct2DPtr)
93 {
94  direction dir = 3;
95 
96  if (correct2DPtr)
97  {
98  const vector& normal = correct2DPtr->planeNormal();
99 
100  if (mag(normal & vector(1, 0, 0)) > 1-edgeTol)
101  {
102  dir = 0;
103  }
104  else if (mag(normal & vector(0, 1, 0)) > 1-edgeTol)
105  {
106  dir = 1;
107  }
108  else if (mag(normal & vector(0, 0, 1)) > 1-edgeTol)
109  {
110  dir = 2;
111  }
112  }
113  return dir;
114 }
115 
116 
117 
118 // Calculate some edge statistics on mesh. Return min. edge length over all
119 // directions but exclude component (0=x, 1=y, 2=z, other=none)
120 scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
121 {
122  label nX = 0;
123  label nY = 0;
124  label nZ = 0;
125 
126  scalar minX = GREAT;
127  scalar maxX = -GREAT;
128  vector x(1, 0, 0);
129 
130  scalar minY = GREAT;
131  scalar maxY = -GREAT;
132  vector y(0, 1, 0);
133 
134  scalar minZ = GREAT;
135  scalar maxZ = -GREAT;
136  vector z(0, 0, 1);
137 
138  scalar minOther = GREAT;
139  scalar maxOther = -GREAT;
140 
141  const edgeList& edges = mesh.edges();
142 
143  forAll(edges, edgeI)
144  {
145  const edge& e = edges[edgeI];
146 
147  vector eVec(e.vec(mesh.points()));
148 
149  scalar eMag = mag(eVec);
150 
151  eVec /= eMag;
152 
153  if (mag(eVec & x) > 1-edgeTol)
154  {
155  minX = min(minX, eMag);
156  maxX = max(maxX, eMag);
157  nX++;
158  }
159  else if (mag(eVec & y) > 1-edgeTol)
160  {
161  minY = min(minY, eMag);
162  maxY = max(maxY, eMag);
163  nY++;
164  }
165  else if (mag(eVec & z) > 1-edgeTol)
166  {
167  minZ = min(minZ, eMag);
168  maxZ = max(maxZ, eMag);
169  nZ++;
170  }
171  else
172  {
173  minOther = min(minOther, eMag);
174  maxOther = max(maxOther, eMag);
175  }
176  }
177 
178  Info<< "Mesh bounding box:" << boundBox(mesh.points()) << nl << nl
179  << "Mesh edge statistics:" << nl
180  << " x aligned : number:" << nX << "\tminLen:" << minX
181  << "\tmaxLen:" << maxX << nl
182  << " y aligned : number:" << nY << "\tminLen:" << minY
183  << "\tmaxLen:" << maxY << nl
184  << " z aligned : number:" << nZ << "\tminLen:" << minZ
185  << "\tmaxLen:" << maxZ << nl
186  << " other : number:" << mesh.nEdges() - nX - nY - nZ
187  << "\tminLen:" << minOther
188  << "\tmaxLen:" << maxOther << nl << endl;
189 
190  if (excludeCmpt == 0)
191  {
192  return min(minY, min(minZ, minOther));
193  }
194  else if (excludeCmpt == 1)
195  {
196  return min(minX, min(minZ, minOther));
197  }
198  else if (excludeCmpt == 2)
199  {
200  return min(minX, min(minY, minOther));
201  }
202  else
203  {
204  return min(minX, min(minY, min(minZ, minOther)));
205  }
206 }
207 
208 
209 // Adds empty patch if not yet there. Returns patchID.
210 label addPatch(polyMesh& mesh, const word& patchName)
211 {
212  label patchI = mesh.boundaryMesh().findPatchID(patchName);
213 
214  if (patchI == -1)
215  {
216  const polyBoundaryMesh& patches = mesh.boundaryMesh();
217 
218  List<polyPatch*> newPatches(patches.size() + 1);
219 
220  // Add empty patch as 0th entry (Note: only since subsetMesh wants this)
221  patchI = 0;
222 
223  newPatches[patchI] =
224  new emptyPolyPatch
225  (
226  Foam::word(patchName),
227  0,
228  mesh.nInternalFaces(),
229  patchI,
230  patches
231  );
232 
233  forAll(patches, i)
234  {
235  const polyPatch& pp = patches[i];
236 
237  newPatches[i+1] =
238  pp.clone
239  (
240  patches,
241  i+1,
242  pp.size(),
243  pp.start()
244  ).ptr();
245  }
246 
247  mesh.removeBoundary();
248  mesh.addPatches(newPatches);
249 
250  Info<< "Created patch oldInternalFaces at " << patchI << endl;
251  }
252  else
253  {
254  Info<< "Reusing patch oldInternalFaces at " << patchI << endl;
255  }
256 
257 
258  return patchI;
259 }
260 
261 
262 // Take surface and select cells based on surface curvature.
263 void selectCurvatureCells
264 (
265  const polyMesh& mesh,
266  const fileName& surfName,
267  const triSurfaceSearch& querySurf,
268  const scalar nearDist,
269  const scalar curvature,
270  cellSet& cells
271 )
272 {
273  // Use surfaceToCell to do actual calculation.
274 
275  // Since we're adding make sure set is on disk.
276  cells.write();
277 
278  // Take centre of cell 0 as outside point since info not used.
279 
280  surfaceToCell cutSource
281  (
282  mesh,
283  surfName,
284  querySurf.surface(),
285  querySurf,
286  pointField(1, mesh.cellCentres()[0]),
287  false,
288  false,
289  false,
290  nearDist,
291  curvature
292  );
293 
294  cutSource.applyToSet(topoSetSource::ADD, cells);
295 }
296 
297 
298 // cutCells contains currently selected cells to be refined. Add neighbours
299 // on the inside or outside to them.
300 void addCutNeighbours
301 (
302  const polyMesh& mesh,
303  const bool selectInside,
304  const bool selectOutside,
305  const labelHashSet& inside,
306  const labelHashSet& outside,
307  labelHashSet& cutCells
308 )
309 {
310  // Pick up face neighbours of cutCells
311 
312  labelHashSet addCutFaces(cutCells.size());
313 
314  for
315  (
316  labelHashSet::const_iterator iter = cutCells.begin();
317  iter != cutCells.end();
318  ++iter
319  )
320  {
321  label cellI = iter.key();
322 
323  const labelList& cFaces = mesh.cells()[cellI];
324 
325  forAll(cFaces, i)
326  {
327  label faceI = cFaces[i];
328 
329  if (mesh.isInternalFace(faceI))
330  {
331  label nbr = mesh.faceOwner()[faceI];
332 
333  if (nbr == cellI)
334  {
335  nbr = mesh.faceNeighbour()[faceI];
336  }
337 
338  if (selectInside && inside.found(nbr))
339  {
340  addCutFaces.insert(nbr);
341  }
342  else if (selectOutside && outside.found(nbr))
343  {
344  addCutFaces.insert(nbr);
345  }
346  }
347  }
348  }
349 
350  Info<< " Selected an additional " << addCutFaces.size()
351  << " neighbours of cutCells to refine" << endl;
352 
353  for
354  (
355  labelHashSet::const_iterator iter = addCutFaces.begin();
356  iter != addCutFaces.end();
357  ++iter
358  )
359  {
360  cutCells.insert(iter.key());
361  }
362 }
363 
364 
365 // Return true if any cells had to be split to keep a difference between
366 // neighbouring refinement levels < limitDiff.
367 // Gets cells which will be refined (so increase the refinement level) and
368 // updates it.
369 bool limitRefinementLevel
370 (
371  const primitiveMesh& mesh,
372  const label limitDiff,
373  const labelHashSet& excludeCells,
374  const labelList& refLevel,
375  labelHashSet& cutCells
376 )
377 {
378  // Do simple check on validity of refinement level.
379  forAll(refLevel, cellI)
380  {
381  if (!excludeCells.found(cellI))
382  {
383  const labelList& cCells = mesh.cellCells()[cellI];
384 
385  forAll(cCells, i)
386  {
387  label nbr = cCells[i];
388 
389  if (!excludeCells.found(nbr))
390  {
391  if (refLevel[cellI] - refLevel[nbr] >= limitDiff)
392  {
393  FatalErrorIn("limitRefinementLevel")
394  << "Level difference between neighbouring cells "
395  << cellI << " and " << nbr
396  << " greater than or equal to " << limitDiff << endl
397  << "refLevels:" << refLevel[cellI] << ' '
398  << refLevel[nbr] << abort(FatalError);
399  }
400  }
401  }
402  }
403  }
404 
405 
406  labelHashSet addCutCells(cutCells.size());
407 
408  for
409  (
410  labelHashSet::const_iterator iter = cutCells.begin();
411  iter != cutCells.end();
412  ++iter
413  )
414  {
415  // cellI will be refined.
416  label cellI = iter.key();
417 
418  const labelList& cCells = mesh.cellCells()[cellI];
419 
420  forAll(cCells, i)
421  {
422  label nbr = cCells[i];
423 
424  if (!excludeCells.found(nbr) && !cutCells.found(nbr))
425  {
426  if (refLevel[cellI] + 1 - refLevel[nbr] >= limitDiff)
427  {
428  addCutCells.insert(nbr);
429  }
430  }
431  }
432  }
433 
434  if (addCutCells.size())
435  {
436  // Add cells to cutCells.
437 
438  Info<< "Added an additional " << addCutCells.size() << " cells"
439  << " to satisfy 1:" << limitDiff << " refinement level"
440  << endl;
441  for
442  (
443  labelHashSet::const_iterator iter = addCutCells.begin();
444  iter != addCutCells.end();
445  ++iter
446  )
447  {
448  cutCells.insert(iter.key());
449  }
450  return true;
451  }
452  else
453  {
454  Info<< "Added no additional cells"
455  << " to satisfy 1:" << limitDiff << " refinement level"
456  << endl;
457 
458  return false;
459  }
460 }
461 
462 
463 // Do all refinement (i.e. refCells) according to refineDict and update
464 // refLevel afterwards for added cells
465 void doRefinement
466 (
467  polyMesh& mesh,
468  const dictionary& refineDict,
469  const labelHashSet& refCells,
470  labelList& refLevel
471 )
472 {
473  label oldCells = mesh.nCells();
474 
475  // Multi-iteration, multi-direction topology modifier.
476  multiDirRefinement multiRef
477  (
478  mesh,
479  refCells.toc(),
480  refineDict
481  );
482 
483  //
484  // Update refLevel for split cells
485  //
486 
487  refLevel.setSize(mesh.nCells());
488 
489  for (label cellI = oldCells; cellI < mesh.nCells(); cellI++)
490  {
491  refLevel[cellI] = 0;
492  }
493 
494  const labelListList& addedCells = multiRef.addedCells();
495 
496  forAll(addedCells, oldCellI)
497  {
498  const labelList& added = addedCells[oldCellI];
499 
500  if (added.size())
501  {
502  // Give all cells resulting from split the refinement level
503  // of the master.
504  label masterLevel = ++refLevel[oldCellI];
505 
506  forAll(added, i)
507  {
508  refLevel[added[i]] = masterLevel;
509  }
510  }
511  }
512 }
513 
514 
515 // Subset mesh and update refLevel and cutCells
516 void subsetMesh
517 (
518  polyMesh& mesh,
519  const label writeMesh,
520  const label patchI, // patchID for exposed faces
521  const labelHashSet& cellsToRemove,
522  cellSet& cutCells,
523  labelIOList& refLevel
524 )
525 {
526  removeCells cellRemover(mesh);
527 
528  labelList cellLabels(cellsToRemove.toc());
529 
530  Info<< "Mesh has:" << mesh.nCells() << " cells."
531  << " Removing:" << cellLabels.size() << " cells" << endl;
532 
533  // exposed faces
534  labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
535 
536  polyTopoChange meshMod(mesh);
537  cellRemover.setRefinement
538  (
539  cellLabels,
540  exposedFaces,
541  labelList(exposedFaces.size(), patchI),
542  meshMod
543  );
544 
545  // Do all changes
546  Info<< "Morphing ..." << endl;
547 
548  const Time& runTime = mesh.time();
549 
550  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
551 
552  if (morphMap().hasMotionPoints())
553  {
554  mesh.movePoints(morphMap().preMotionPoints());
555  }
556 
557  // Update topology on cellRemover
558  cellRemover.updateMesh(morphMap());
559 
560  // Update refLevel for removed cells.
561  const labelList& cellMap = morphMap().cellMap();
562 
563  labelList newRefLevel(cellMap.size());
564 
565  forAll(cellMap, i)
566  {
567  newRefLevel[i] = refLevel[cellMap[i]];
568  }
569 
570  // Transfer back to refLevel
571  refLevel.transfer(newRefLevel);
572 
573  if (writeMesh)
574  {
575  Info<< "Writing refined mesh to time " << runTime.timeName() << nl
576  << endl;
577 
579  mesh.write();
580  refLevel.write();
581  }
582 
583  // Update cutCells for removed cells.
584  cutCells.updateMesh(morphMap());
585 }
586 
587 
588 // Divide the cells according to status compared to surface. Constructs sets:
589 // - cutCells : all cells cut by surface
590 // - inside : all cells inside surface
591 // - outside : ,, outside ,,
592 // and a combined set:
593 // - selected : sum of above according to selectCut/Inside/Outside flags.
594 void classifyCells
595 (
596  const polyMesh& mesh,
597  const fileName& surfName,
598  const triSurfaceSearch& querySurf,
599  const pointField& outsidePts,
600 
601  const bool selectCut,
602  const bool selectInside,
603  const bool selectOutside,
604 
605  const label nCutLayers,
606 
607  cellSet& inside,
608  cellSet& outside,
609  cellSet& cutCells,
610  cellSet& selected
611 )
612 {
613  // Cut faces with surface and classify cells
615  (
616  mesh,
617  surfName,
618  querySurf.surface(),
619  querySurf,
620  outsidePts,
621 
622  nCutLayers,
623 
624  inside,
625  outside,
626  cutCells
627  );
628 
629  // Combine wanted parts into selected
630  if (selectCut)
631  {
632  selected.addSet(cutCells);
633  }
634  if (selectInside)
635  {
636  selected.addSet(inside);
637  }
638  if (selectOutside)
639  {
640  selected.addSet(outside);
641  }
642 
643  Info<< "Determined cell status:" << endl
644  << " inside :" << inside.size() << endl
645  << " outside :" << outside.size() << endl
646  << " cutCells:" << cutCells.size() << endl
647  << " selected:" << selected.size() << endl
648  << endl;
649 
650  writeSet(inside, "inside cells");
651  writeSet(outside, "outside cells");
652  writeSet(cutCells, "cut cells");
653  writeSet(selected, "selected cells");
654 }
655 
656 
657 // Main program:
658 
659 int main(int argc, char *argv[])
660 {
662 
663 # include <OpenFOAM/setRootCase.H>
664 # include <OpenFOAM/createTime.H>
665 # include <OpenFOAM/createPolyMesh.H>
666 
667  // If nessecary add oldInternalFaces patch
668  label newPatchI = addPatch(mesh, "oldInternalFaces");
669 
670 
671  //
672  // Read motionProperties dictionary
673  //
674 
675  Info<< "Checking for motionProperties\n" << endl;
676 
677  IOobject motionObj
678  (
679  "motionProperties",
680  runTime.constant(),
681  mesh,
684  );
685 
686  // corrector for mesh motion
687  twoDPointCorrector* correct2DPtr = NULL;
688 
689  if (motionObj.headerOk())
690  {
691  Info<< "Reading " << runTime.constant() / "motionProperties"
692  << endl << endl;
693 
694  IOdictionary motionProperties(motionObj);
695 
696  Switch twoDMotion(motionProperties.lookup("twoDMotion"));
697 
698  if (twoDMotion)
699  {
700  Info<< "Correcting for 2D motion" << endl << endl;
701  correct2DPtr = new twoDPointCorrector(mesh);
702  }
703  }
704 
705  IOdictionary refineDict
706  (
707  IOobject
708  (
709  "autoRefineMeshDict",
710  runTime.system(),
711  mesh,
714  )
715  );
716 
717  fileName surfName(refineDict.lookup("surface"));
718  surfName.expand();
719  label nCutLayers(readLabel(refineDict.lookup("nCutLayers")));
720  label cellLimit(readLabel(refineDict.lookup("cellLimit")));
721  bool selectCut(readBool(refineDict.lookup("selectCut")));
722  bool selectInside(readBool(refineDict.lookup("selectInside")));
723  bool selectOutside(readBool(refineDict.lookup("selectOutside")));
724  bool selectHanging(readBool(refineDict.lookup("selectHanging")));
725 
726  scalar minEdgeLen(readScalar(refineDict.lookup("minEdgeLen")));
727  scalar maxEdgeLen(readScalar(refineDict.lookup("maxEdgeLen")));
728  scalar curvature(readScalar(refineDict.lookup("curvature")));
729  scalar curvDist(readScalar(refineDict.lookup("curvatureDistance")));
730  pointField outsidePts(refineDict.lookup("outsidePoints"));
731  label refinementLimit(readLabel(refineDict.lookup("splitLevel")));
732  bool writeMesh(readBool(refineDict.lookup("writeMesh")));
733 
734  Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
735  << " cells cut by surface : " << selectCut << nl
736  << " cells inside of surface : " << selectInside << nl
737  << " cells outside of surface : " << selectOutside << nl
738  << " hanging cells : " << selectHanging << nl
739  << endl;
740 
741 
742  if (nCutLayers > 0 && selectInside)
743  {
745  << "Illogical settings : Both nCutLayers>0 and selectInside true."
746  << endl
747  << "This would mean that inside cells get removed but should be"
748  << " included in final mesh" << endl;
749  }
750 
751  // Surface.
752  triSurface surf(surfName);
753 
754  // Search engine on surface
755  triSurfaceSearch querySurf(surf);
756 
757  // Search engine on mesh. No face decomposition since mesh unwarped.
758  meshSearch queryMesh(mesh, false);
759 
760  // Check all 'outside' points
761  forAll(outsidePts, outsideI)
762  {
763  const point& outsidePoint = outsidePts[outsideI];
764 
765  if (queryMesh.findCell(outsidePoint, -1, false) == -1)
766  {
768  << "outsidePoint " << outsidePoint
769  << " is not inside any cell"
770  << exit(FatalError);
771  }
772  }
773 
774 
775 
776  // Current refinement level. Read if present.
777  labelIOList refLevel
778  (
779  IOobject
780  (
781  "refinementLevel",
782  runTime.timeName(),
784  mesh,
787  ),
788  labelList(mesh.nCells(), 0)
789  );
790 
791  label maxLevel = max(refLevel);
792 
793  if (maxLevel > 0)
794  {
795  Info<< "Read existing refinement level from file "
796  << refLevel.objectPath() << nl
797  << " min level : " << min(refLevel) << nl
798  << " max level : " << maxLevel << nl
799  << endl;
800  }
801  else
802  {
803  Info<< "Created zero refinement level in file "
804  << refLevel.objectPath() << nl
805  << endl;
806  }
807 
808 
809 
810 
811  // Print edge stats on original mesh. Leave out 2d-normal direction
812  direction normalDir(getNormalDir(correct2DPtr));
813  scalar meshMinEdgeLen = getEdgeStats(mesh, normalDir);
814 
815  while (meshMinEdgeLen > minEdgeLen)
816  {
817  // Get inside/outside/cutCells cellSets.
818  cellSet inside(mesh, "inside", mesh.nCells()/10);
819  cellSet outside(mesh, "outside", mesh.nCells()/10);
820  cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
821  cellSet selected(mesh, "selected", mesh.nCells()/10);
822 
823  classifyCells
824  (
825  mesh,
826  surfName,
827  querySurf,
828  outsidePts,
829 
830  selectCut, // for determination of selected
831  selectInside, // for determination of selected
832  selectOutside, // for determination of selected
833 
834  nCutLayers, // How many layers of cutCells to include
835 
836  inside,
837  outside,
838  cutCells,
839  selected // not used but determined anyway.
840  );
841 
842  Info<< " Selected " << cutCells.name() << " with "
843  << cutCells.size() << " cells" << endl;
844 
845  if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
846  {
847  // Done refining enough close to surface. Now switch to more
848  // intelligent refinement where only cutCells on surface curvature
849  // are refined.
850  cellSet curveCells(mesh, "curveCells", mesh.nCells()/10);
851 
852  selectCurvatureCells
853  (
854  mesh,
855  surfName,
856  querySurf,
857  maxEdgeLen,
858  curvature,
859  curveCells
860  );
861 
862  Info<< " Selected " << curveCells.name() << " with "
863  << curveCells.size() << " cells" << endl;
864 
865  // Add neighbours to cutCells. This is if selectCut is not
866  // true and so outside and/or inside cells get exposed there is
867  // also refinement in them.
868  if (!selectCut)
869  {
870  addCutNeighbours
871  (
872  mesh,
873  selectInside,
874  selectOutside,
875  inside,
876  outside,
877  cutCells
878  );
879  }
880 
881  // Subset cutCells to only curveCells
882  cutCells.subset(curveCells);
883 
884  Info<< " Removed cells not on surface curvature. Selected "
885  << cutCells.size() << endl;
886  }
887 
888 
889  if (nCutLayers > 0)
890  {
891  // Subset mesh to remove inside cells altogether. Updates cutCells,
892  // refLevel.
893  subsetMesh(mesh, writeMesh, newPatchI, inside, cutCells, refLevel);
894  }
895 
896 
897  // Added cells from 2:1 refinement level restriction.
898  while
899  (
900  limitRefinementLevel
901  (
902  mesh,
903  refinementLimit,
904  labelHashSet(),
905  refLevel,
906  cutCells
907  )
908  )
909  {}
910 
911 
912  Info<< " Current cells : " << mesh.nCells() << nl
913  << " Selected for refinement :" << cutCells.size() << nl
914  << endl;
915 
916  if (cutCells.empty())
917  {
918  Info<< "Stopping refining since 0 cells selected to be refined ..."
919  << nl << endl;
920  break;
921  }
922 
923  if ((mesh.nCells() + 8*cutCells.size()) > cellLimit)
924  {
925  Info<< "Stopping refining since cell limit reached ..." << nl
926  << "Would refine from " << mesh.nCells()
927  << " to " << mesh.nCells() + 8*cutCells.size() << " cells."
928  << nl << endl;
929  break;
930  }
931 
932  doRefinement
933  (
934  mesh,
935  refineDict,
936  cutCells,
937  refLevel
938  );
939 
940  Info<< " After refinement:" << mesh.nCells() << nl
941  << endl;
942 
943  if (writeMesh)
944  {
945  Info<< " Writing refined mesh to time " << runTime.timeName()
946  << nl << endl;
947 
949  mesh.write();
950  refLevel.write();
951  }
952 
953  // Update mesh edge stats.
954  meshMinEdgeLen = getEdgeStats(mesh, normalDir);
955  }
956 
957 
958  if (selectHanging)
959  {
960  // Get inside/outside/cutCells cellSets.
961  cellSet inside(mesh, "inside", mesh.nCells()/10);
962  cellSet outside(mesh, "outside", mesh.nCells()/10);
963  cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
964  cellSet selected(mesh, "selected", mesh.nCells()/10);
965 
966  classifyCells
967  (
968  mesh,
969  surfName,
970  querySurf,
971  outsidePts,
972 
973  selectCut,
974  selectInside,
975  selectOutside,
976 
977  nCutLayers,
978 
979  inside,
980  outside,
981  cutCells,
982  selected
983  );
984 
985 
986  // Find any cells which have all their points on the outside of the
987  // selected set and refine them
988  labelHashSet hanging = surfaceSets::getHangingCells(mesh, selected);
989 
990  Info<< "Detected " << hanging.size() << " hanging cells"
991  << " (cells with all points on"
992  << " outside of cellSet 'selected').\nThese would probably result"
993  << " in flattened cells when snapping the mesh to the surface"
994  << endl;
995 
996  Info<< "Refining " << hanging.size() << " hanging cells" << nl
997  << endl;
998 
999  // Added cells from 2:1 refinement level restriction.
1000  while
1001  (
1002  limitRefinementLevel
1003  (
1004  mesh,
1005  refinementLimit,
1006  labelHashSet(),
1007  refLevel,
1008  hanging
1009  )
1010  )
1011  {}
1012 
1013  doRefinement(mesh, refineDict, hanging, refLevel);
1014 
1015  Info<< "Writing refined mesh to time " << runTime.timeName() << nl
1016  << endl;
1017 
1018  // Write final mesh
1020  mesh.write();
1021  refLevel.write();
1022 
1023  }
1024  else if (!writeMesh)
1025  {
1026  Info<< "Writing refined mesh to time " << runTime.timeName() << nl
1027  << endl;
1028 
1029  // Write final mesh. (will have been written already if writeMesh=true)
1031  mesh.write();
1032  refLevel.write();
1033  }
1034 
1035 
1036  Info << "End\n" << endl;
1037 
1038  return 0;
1039 }
1040 
1041 
1042 // ************************ vim: set sw=4 sts=4 et: ************************ //