FreeFOAM The Cross-Platform CFD Toolkit
autoHexMeshDriver.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 "autoHexMeshDriver.H"
27 #include <finiteVolume/fvMesh.H>
28 #include <OpenFOAM/Time.H>
29 #include <OpenFOAM/boundBox.H>
30 #include <OpenFOAM/wallPolyPatch.H>
31 #include <meshTools/cellSet.H>
32 #include <OpenFOAM/syncTools.H>
36 #include "autoRefineDriver.H"
37 #include "autoSnapDriver.H"
38 #include "autoLayerDriver.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
46 }
47 
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 // Check writing tolerance before doing any serious work
52 Foam::scalar Foam::autoHexMeshDriver::getMergeDistance(const scalar mergeTol)
53  const
54 {
55  const boundBox& meshBb = mesh_.bounds();
56  scalar mergeDist = mergeTol * meshBb.mag();
57  scalar writeTol = std::pow
58  (
59  scalar(10.0),
61  );
62 
63  Info<< nl
64  << "Overall mesh bounding box : " << meshBb << nl
65  << "Relative tolerance : " << mergeTol << nl
66  << "Absolute matching distance : " << mergeDist << nl
67  << endl;
68 
69  if (mesh_.time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
70  {
71  FatalErrorIn("autoHexMeshDriver::getMergeDistance(const scalar) const")
72  << "Your current settings specify ASCII writing with "
73  << IOstream::defaultPrecision() << " digits precision." << endl
74  << "Your merging tolerance (" << mergeTol << ") is finer than this."
75  << endl
76  << "Please change your writeFormat to binary"
77  << " or increase the writePrecision" << endl
78  << "or adjust the merge tolerance (-mergeTol)."
79  << exit(FatalError);
80  }
81 
82  return mergeDist;
83 }
84 
85 
87 //void Foam::autoHexMeshDriver::orientOutside
88 //(
89 // PtrList<searchableSurface>& shells
90 //)
91 //{
92 // // Determine outside point.
93 // boundBox overallBb = boundBox::invertedBox;
94 //
95 // bool hasSurface = false;
96 //
97 // forAll(shells, shellI)
98 // {
99 // if (isA<triSurfaceMesh>(shells[shellI]))
100 // {
101 // const triSurfaceMesh& shell =
102 // refCast<const triSurfaceMesh>(shells[shellI]);
103 //
104 // hasSurface = true;
105 //
106 // boundBox shellBb(shell.localPoints(), false);
107 //
108 // overallBb.min() = min(overallBb.min(), shellBb.min());
109 // overallBb.max() = max(overallBb.max(), shellBb.max());
110 // }
111 // }
112 //
113 // if (hasSurface)
114 // {
115 // const point outsidePt = 2 * overallBb.span();
116 //
117 // //Info<< "Using point " << outsidePt << " to orient shells" << endl;
118 //
119 // forAll(shells, shellI)
120 // {
121 // if (isA<triSurfaceMesh>(shells[shellI]))
122 // {
123 // triSurfaceMesh& shell =
124 // refCast<triSurfaceMesh>(shells[shellI]);
125 //
126 // if (!refinementSurfaces::isSurfaceClosed(shell))
127 // {
128 // FatalErrorIn("orientOutside(PtrList<searchableSurface>&)")
129 // << "Refinement shell "
130 // << shell.searchableSurface::name()
131 // << " is not closed." << exit(FatalError);
132 // }
133 //
134 // refinementSurfaces::orientSurface(outsidePt, shell);
135 // }
136 // }
137 // }
138 //}
139 
140 
141 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
142 
143 // Construct from components
144 Foam::autoHexMeshDriver::autoHexMeshDriver
145 (
146  fvMesh& mesh,
147  const bool overwrite,
148  const dictionary& dict,
149  const dictionary& decomposeDict
150 )
151 :
152  mesh_(mesh),
153  dict_(dict),
154  debug_(readLabel(dict_.lookup("debug"))),
155  mergeDist_(getMergeDistance(readScalar(dict_.lookup("mergeTolerance"))))
156 {
157  if (debug_ > 0)
158  {
159  meshRefinement::debug = debug_;
160  autoHexMeshDriver::debug = debug_;
161  autoRefineDriver::debug = debug;
162  autoSnapDriver::debug = debug;
163  autoLayerDriver::debug = debug;
164  }
165 
166  refinementParameters refineParams(dict, 1);
167 
168  Info<< "Overall cell limit : "
169  << refineParams.maxGlobalCells() << endl;
170  Info<< "Per processor cell limit : "
171  << refineParams.maxLocalCells() << endl;
172  Info<< "Minimum number of cells to refine : "
173  << refineParams.minRefineCells() << endl;
174  Info<< "Curvature : "
175  << refineParams.curvature() << nl << endl;
176  Info<< "Layers between different refinement levels : "
177  << refineParams.nBufferLayers() << endl;
178 
179  PtrList<dictionary> shellDicts(dict_.lookup("refinementShells"));
180 
181  PtrList<dictionary> surfaceDicts(dict_.lookup("surfaces"));
182 
183 
184  // Read geometry
185  // ~~~~~~~~~~~~~
186 
187  {
188  Info<< "Reading all geometry." << endl;
189 
190  // Construct dictionary with all shells and all refinement surfaces
191  dictionary geometryDict;
192 
193  forAll(shellDicts, shellI)
194  {
195  dictionary shellDict = shellDicts[shellI];
196  const word name(shellDict.lookup("name"));
197  shellDict.remove("name");
198  shellDict.remove("level");
199  shellDict.remove("refineInside");
200  geometryDict.add(name, shellDict);
201  }
202 
203  forAll(surfaceDicts, surfI)
204  {
205  dictionary surfDict = surfaceDicts[surfI];
206  const word name(string::validate<word>(surfDict.lookup("file")));
207  surfDict.remove("file");
208  surfDict.remove("regions");
209  if (!surfDict.found("name"))
210  {
211  surfDict.add("name", name);
212  }
213  surfDict.add("type", triSurfaceMesh::typeName);
214  geometryDict.add(name, surfDict);
215  }
216 
217  allGeometryPtr_.reset
218  (
220  (
221  IOobject
222  (
223  "abc", // dummy name
224  //mesh_.time().findInstance("triSurface", word::null),
225  // instance
226  mesh_.time().constant(), // instance
227  "triSurface", // local
228  mesh_.time(), // registry
231  ),
232  geometryDict
233  )
234  );
235 
236  Info<< "Read geometry in = "
237  << mesh_.time().cpuTimeIncrement() << " s" << endl;
238  }
239 
240 
241  // Read refinement surfaces
242  // ~~~~~~~~~~~~~~~~~~~~~~~~
243 
244  {
245  Info<< "Reading surfaces and constructing search trees." << endl;
246 
247  surfacesPtr_.reset
248  (
250  (
251  allGeometryPtr_(),
252  surfaceDicts
253  )
254  );
255  Info<< "Read surfaces in = "
256  << mesh_.time().cpuTimeIncrement() << " s" << endl;
257  }
258 
259  // Read refinement shells
260  // ~~~~~~~~~~~~~~~~~~~~~~
261 
262  {
263  Info<< "Reading refinement shells." << endl;
264  shellsPtr_.reset
265  (
266  new shellSurfaces
267  (
268  allGeometryPtr_(),
269  shellDicts
270  )
271  );
272  Info<< "Read refinement shells in = "
273  << mesh_.time().cpuTimeIncrement() << " s" << endl;
274 
276  //Info<< "Orienting triSurface shells so point far away is outside."
277  // << endl;
278  //orientOutside(shells_);
279  //Info<< "Oriented shells in = "
280  // << mesh_.time().cpuTimeIncrement() << " s" << endl;
281 
282  Info<< "Setting refinement level of surface to be consistent"
283  << " with shells." << endl;
284  surfacesPtr_().setMinLevelFields(shells());
285  Info<< "Checked shell refinement in = "
286  << mesh_.time().cpuTimeIncrement() << " s" << endl;
287  }
288 
289 
290  // Check faceZones are synchronised
292 
293 
294  // Refinement engine
295  // ~~~~~~~~~~~~~~~~~
296 
297  {
298  Info<< nl
299  << "Determining initial surface intersections" << nl
300  << "-----------------------------------------" << nl
301  << endl;
302 
303  // Main refinement engine
304  meshRefinerPtr_.reset
305  (
306  new meshRefinement
307  (
308  mesh,
309  mergeDist_, // tolerance used in sorting coordinates
310  overwrite,
311  surfaces(),
312  shells()
313  )
314  );
315  Info<< "Calculated surface intersections in = "
316  << mesh_.time().cpuTimeIncrement() << " s" << endl;
317 
318  // Some stats
319  meshRefinerPtr_().printMeshInfo(debug_, "Initial mesh");
320 
321  meshRefinerPtr_().write
322  (
324  mesh_.time().path()/meshRefinerPtr_().timeName()
325  );
326  }
327 
328 
329  // Add all the surface regions as patches
330  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
331 
332  {
333  Info<< nl
334  << "Adding patches for surface regions" << nl
335  << "----------------------------------" << nl
336  << endl;
337 
338  // From global region number to mesh patch.
339  globalToPatch_.setSize(surfaces().nRegions(), -1);
340 
341  Info<< "Patch\tRegion" << nl
342  << "-----\t------"
343  << endl;
344 
345  const labelList& surfaceGeometry = surfaces().surfaces();
346  forAll(surfaceGeometry, surfI)
347  {
348  label geomI = surfaceGeometry[surfI];
349 
350  const wordList& regNames = allGeometryPtr_().regionNames()[geomI];
351 
352  Info<< surfaces().names()[surfI] << ':' << nl << nl;
353 
354  forAll(regNames, i)
355  {
356  label patchI = meshRefinerPtr_().addMeshedPatch
357  (
358  regNames[i],
359  wallPolyPatch::typeName
360  );
361 
362  Info<< patchI << '\t' << regNames[i] << nl;
363 
364  globalToPatch_[surfaces().globalRegion(surfI, i)] = patchI;
365  }
366 
367  Info<< nl;
368  }
369  Info<< "Added patches in = "
370  << mesh_.time().cpuTimeIncrement() << " s" << nl << endl;
371  }
372 
373 
378  //
379  //labelList namedSurfaces(surfaces().getNamedSurfaces());
380  //if (namedSurfaces.size())
381  //{
382  // Info<< nl
383  // << "Introducing cyclics for faceZones" << nl
384  // << "---------------------------------" << nl
385  // << endl;
386  //
387  // // From surface to cyclic patch
388  // surfaceToCyclicPatch_.setSize(surfaces().size(), -1);
389  //
390  // Info<< "Patch\tZone" << nl
391  // << "----\t-----"
392  // << endl;
393  //
394  // forAll(namedSurfaces, i)
395  // {
396  // label surfI = namedSurfaces[i];
397  //
398  // surfaceToCyclicPatch_[surfI] = meshRefinement::addPatch
399  // (
400  // mesh,
401  // surfaces().faceZoneNames()[surfI],
402  // cyclicPolyPatch::typeName
403  // );
404  //
405  // Info<< surfaceToCyclicPatch_[surfI] << '\t'
406  // << surfaces().faceZoneNames()[surfI] << nl << endl;
407  // }
408  // Info<< "Added cyclic patches in = "
409  // << mesh_.time().cpuTimeIncrement() << " s" << endl;
410  //}
411 
412 
413  // Parallel
414  // ~~~~~~~~
415 
416  {
417  // Decomposition
418  decomposerPtr_ = decompositionMethod::New
419  (
420  decomposeDict,
421  mesh_
422  );
423  decompositionMethod& decomposer = decomposerPtr_();
424 
425 
426  if (Pstream::parRun() && !decomposer.parallelAware())
427  {
428  FatalErrorIn("autoHexMeshDriver::autoHexMeshDriver"
429  "(const IOobject&, fvMesh&)")
430  << "You have selected decomposition method "
431  << decomposer.typeName
432  << " which is not parallel aware." << endl
433  << "Please select one that is (parMetis, hierarchical)"
434  << exit(FatalError);
435  }
436 
437  // Mesh distribution engine (uses tolerance to reconstruct meshes)
438  distributorPtr_.reset(new fvMeshDistribute(mesh_, mergeDist_));
439  }
440 }
441 
442 
443 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
444 
445 void Foam::autoHexMeshDriver::writeMesh(const string& msg) const
446 {
447  const meshRefinement& meshRefiner = meshRefinerPtr_();
448 
449  meshRefiner.printMeshInfo(debug_, msg);
450  Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
451 
454  {
455  meshRefiner.write
456  (
457  meshRefinement::OBJINTERSECTIONS,
458  mesh_.time().path()/meshRefiner.timeName()
459  );
460  }
461  Info<< "Written mesh in = "
462  << mesh_.time().cpuTimeIncrement() << " s." << endl;
463 }
464 
465 
467 {
468  Switch wantRefine(dict_.lookup("doRefine"));
469  Switch wantSnap(dict_.lookup("doSnap"));
470  Switch wantLayers(dict_.lookup("doLayers"));
471 
472  Info<< "Do refinement : " << wantRefine << nl
473  << "Do snapping : " << wantSnap << nl
474  << "Do layers : " << wantLayers << nl
475  << endl;
476 
477  if (wantRefine)
478  {
479  const dictionary& motionDict = dict_.subDict("motionDict");
480 
481  autoRefineDriver refineDriver
482  (
483  meshRefinerPtr_(),
484  decomposerPtr_(),
485  distributorPtr_(),
486  globalToPatch_
487  );
488 
489  // Get all the refinement specific params
490  refinementParameters refineParams(dict_, 1);
491 
492  refineDriver.doRefine(dict_, refineParams, wantSnap, motionDict);
493 
494  // Write mesh
495  writeMesh("Refined mesh");
496  }
497 
498  if (wantSnap)
499  {
500  const dictionary& snapDict = dict_.subDict("snapDict");
501  const dictionary& motionDict = dict_.subDict("motionDict");
502 
503  autoSnapDriver snapDriver
504  (
505  meshRefinerPtr_(),
506  globalToPatch_
507  );
508 
509  // Get all the snapping specific params
510  snapParameters snapParams(snapDict, 1);
511 
512  snapDriver.doSnap(snapDict, motionDict, snapParams);
513 
514  // Write mesh.
515  writeMesh("Snapped mesh");
516  }
517 
518  if (wantLayers)
519  {
520  const dictionary& motionDict = dict_.subDict("motionDict");
521  const dictionary& shrinkDict = dict_.subDict("shrinkDict");
522  PtrList<dictionary> surfaceDicts(dict_.lookup("surfaces"));
523 
524  autoLayerDriver layerDriver(meshRefinerPtr_());
525 
526  // Get all the layer specific params
527  layerParameters layerParams
528  (
529  surfaceDicts,
530  surfacesPtr_(),
531  globalToPatch_,
532  shrinkDict,
533  mesh_.boundaryMesh()
534  );
535 
536  layerDriver.doLayers
537  (
538  shrinkDict,
539  motionDict,
540  layerParams,
541  true, // pre-balance
542  decomposerPtr_(),
543  distributorPtr_()
544  );
545 
546  // Write mesh.
547  writeMesh("Layer mesh");
548  }
549 }
550 
551 
552 // ************************ vim: set sw=4 sts=4 et: ************************ //