FreeFOAM The Cross-Platform CFD Toolkit
snappyHexMesh.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  snappyHexMesh
26 
27 Description
28  Automatic split hex mesher. Refines and snaps to surface.
29 
30 Usage
31 
32  - snappyHexMesh [OPTIONS]
33 
34  @param -overwrite
35  Overwrite existing data.
36 
37  @param -case <dir>\n
38  Case directory.
39 
40  @param -parallel \n
41  Run in parallel.
42 
43  @param -help \n
44  Display help message.
45 
46  @param -doc \n
47  Display Doxygen API documentation page for this application.
48 
49  @param -srcDoc \n
50  Display Doxygen source documentation page for this application.
51 
52 \*---------------------------------------------------------------------------*/
53 
54 #include <OpenFOAM/argList.H>
55 #include <OpenFOAM/Time.H>
56 #include <finiteVolume/fvMesh.H>
62 #include <autoMesh/shellSurfaces.H>
65 #include <OpenFOAM/wallPolyPatch.H>
69 
70 
71 using namespace Foam;
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 // Check writing tolerance before doing any serious work
76 scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
77 {
78  const boundBox& meshBb = mesh.bounds();
79  scalar mergeDist = mergeTol * meshBb.mag();
80  scalar writeTol = std::pow
81  (
82  scalar(10.0),
84  );
85 
86  Info<< nl
87  << "Overall mesh bounding box : " << meshBb << nl
88  << "Relative tolerance : " << mergeTol << nl
89  << "Absolute matching distance : " << mergeDist << nl
90  << endl;
91 
92  if (mesh.time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
93  {
94  FatalErrorIn("getMergeDistance(const polyMesh&, const scalar)")
95  << "Your current settings specify ASCII writing with "
96  << IOstream::defaultPrecision() << " digits precision." << endl
97  << "Your merging tolerance (" << mergeTol << ") is finer than this."
98  << endl
99  << "Please change your writeFormat to binary"
100  << " or increase the writePrecision" << endl
101  << "or adjust the merge tolerance (-mergeTol)."
102  << exit(FatalError);
103  }
104 
105  return mergeDist;
106 }
107 
108 
109 // Write mesh and additional information
110 void writeMesh
111 (
112  const string& msg,
113  const meshRefinement& meshRefiner,
114  const label debug
115 )
116 {
117  const fvMesh& mesh = meshRefiner.mesh();
118 
119  meshRefiner.printMeshInfo(debug, msg);
120  Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
121 
124  {
125  meshRefiner.write
126  (
127  meshRefinement::OBJINTERSECTIONS,
128  mesh.time().path()/meshRefiner.timeName()
129  );
130  }
131  Info<< "Written mesh in = "
132  << mesh.time().cpuTimeIncrement() << " s." << endl;
133 }
134 
135 
136 
137 int main(int argc, char *argv[])
138 {
139  argList::validOptions.insert("overwrite", "");
140 # include <OpenFOAM/setRootCase.H>
141 # include <OpenFOAM/createTime.H>
142  runTime.functionObjects().off();
143 # include <OpenFOAM/createMesh.H>
144 
145  Info<< "Read mesh in = "
146  << runTime.cpuTimeIncrement() << " s" << endl;
147 
148  const bool overwrite = args.optionFound("overwrite");
149 
150 
151  // Check patches and faceZones are synchronised
152  mesh.boundaryMesh().checkParallelSync(true);
154 
155 
156  // Read decomposePar dictionary
157  IOdictionary decomposeDict
158  (
159  IOobject
160  (
161  "decomposeParDict",
162  runTime.system(),
163  mesh,
166  )
167  );
168 
169  // Read meshing dictionary
170  IOdictionary meshDict
171  (
172  IOobject
173  (
174  "snappyHexMeshDict",
175  runTime.system(),
176  mesh,
179  )
180  );
181 
182  // all surface geometry
183  const dictionary& geometryDict = meshDict.subDict("geometry");
184 
185  // refinement parameters
186  const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
187 
188  // mesh motion and mesh quality parameters
189  const dictionary& motionDict = meshDict.subDict("meshQualityControls");
190 
191  // snap-to-surface parameters
192  const dictionary& snapDict = meshDict.subDict("snapControls");
193 
194  // layer addition parameters
195  const dictionary& layerDict = meshDict.subDict("addLayersControls");
196 
197 
198  const scalar mergeDist = getMergeDistance
199  (
200  mesh,
201  readScalar(meshDict.lookup("mergeTolerance"))
202  );
203 
204 
205 
206  // Debug
207  // ~~~~~
208 
209  const label debug(readLabel(meshDict.lookup("debug")));
210  if (debug > 0)
211  {
212  meshRefinement::debug = debug;
213  autoRefineDriver::debug = debug;
214  autoSnapDriver::debug = debug;
215  autoLayerDriver::debug = debug;
216  }
217 
218 
219  // Read geometry
220  // ~~~~~~~~~~~~~
221 
222  searchableSurfaces allGeometry
223  (
224  IOobject
225  (
226  "abc", // dummy name
227  mesh.time().constant(), // instance
228  //mesh.time().findInstance("triSurface", word::null),// instance
229  "triSurface", // local
230  mesh.time(), // registry
233  ),
234  geometryDict
235  );
236 
237 
238  // Read refinement surfaces
239  // ~~~~~~~~~~~~~~~~~~~~~~~~
240 
241  Info<< "Reading refinement surfaces." << endl;
242  refinementSurfaces surfaces
243  (
244  allGeometry,
245  refineDict.subDict("refinementSurfaces")
246  );
247  Info<< "Read refinement surfaces in = "
248  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
249 
250 
251  // Read refinement shells
252  // ~~~~~~~~~~~~~~~~~~~~~~
253 
254  Info<< "Reading refinement shells." << endl;
255  shellSurfaces shells
256  (
257  allGeometry,
258  refineDict.subDict("refinementRegions")
259  );
260  Info<< "Read refinement shells in = "
261  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
262 
263 
264  Info<< "Setting refinement level of surface to be consistent"
265  << " with shells." << endl;
266  surfaces.setMinLevelFields(shells);
267  Info<< "Checked shell refinement in = "
268  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
269 
270 
271  // Refinement engine
272  // ~~~~~~~~~~~~~~~~~
273 
274  Info<< nl
275  << "Determining initial surface intersections" << nl
276  << "-----------------------------------------" << nl
277  << endl;
278 
279  // Main refinement engine
280  meshRefinement meshRefiner
281  (
282  mesh,
283  mergeDist, // tolerance used in sorting coordinates
284  overwrite, // overwrite mesh files?
285  surfaces, // for surface intersection refinement
286  shells // for volume (inside/outside) refinement
287  );
288  Info<< "Calculated surface intersections in = "
289  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
290 
291  // Some stats
292  meshRefiner.printMeshInfo(debug, "Initial mesh");
293 
294  meshRefiner.write
295  (
296  debug&meshRefinement::OBJINTERSECTIONS,
297  mesh.time().path()/meshRefiner.timeName()
298  );
299 
300 
301  // Add all the surface regions as patches
302  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
303 
304  labelList globalToPatch;
305  {
306  Info<< nl
307  << "Adding patches for surface regions" << nl
308  << "----------------------------------" << nl
309  << endl;
310 
311  // From global region number to mesh patch.
312  globalToPatch.setSize(surfaces.nRegions(), -1);
313 
314  Info<< "Patch\tRegion" << nl
315  << "-----\t------"
316  << endl;
317 
318  const labelList& surfaceGeometry = surfaces.surfaces();
319  forAll(surfaceGeometry, surfI)
320  {
321  label geomI = surfaceGeometry[surfI];
322 
323  const wordList& regNames = allGeometry.regionNames()[geomI];
324 
325  Info<< surfaces.names()[surfI] << ':' << nl << nl;
326 
327  forAll(regNames, i)
328  {
329  label patchI = meshRefiner.addMeshedPatch
330  (
331  regNames[i],
332  wallPolyPatch::typeName
333  );
334 
335  Info<< patchI << '\t' << regNames[i] << nl;
336 
337  globalToPatch[surfaces.globalRegion(surfI, i)] = patchI;
338  }
339 
340  Info<< nl;
341  }
342  Info<< "Added patches in = "
343  << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
344  }
345 
346 
347  // Parallel
348  // ~~~~~~~~
349 
350  // Decomposition
351  autoPtr<decompositionMethod> decomposerPtr
352  (
354  (
355  decomposeDict,
356  mesh
357  )
358  );
359  decompositionMethod& decomposer = decomposerPtr();
360 
361  if (Pstream::parRun() && !decomposer.parallelAware())
362  {
364  << "You have selected decomposition method "
365  << decomposer.typeName
366  << " which is not parallel aware." << endl
367  << "Please select one that is (hierarchical, parMetis)"
368  << exit(FatalError);
369  }
370 
371  // Mesh distribution engine (uses tolerance to reconstruct meshes)
372  fvMeshDistribute distributor(mesh, mergeDist);
373 
374 
375 
376 
377 
378  // Now do the real work -refinement -snapping -layers
379  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
380 
381  Switch wantRefine(meshDict.lookup("castellatedMesh"));
382  Switch wantSnap(meshDict.lookup("snap"));
383  Switch wantLayers(meshDict.lookup("addLayers"));
384 
385  if (wantRefine)
386  {
387  cpuTime timer;
388 
389  autoRefineDriver refineDriver
390  (
391  meshRefiner,
392  decomposer,
393  distributor,
394  globalToPatch
395  );
396 
397  // Refinement parameters
398  refinementParameters refineParams(refineDict);
399 
400  if (!overwrite)
401  {
402  const_cast<Time&>(mesh.time())++;
403  }
404 
405  refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
406 
407  writeMesh
408  (
409  "Refined mesh",
410  meshRefiner,
411  debug
412  );
413 
414  Info<< "Mesh refined in = "
415  << timer.cpuTimeIncrement() << " s." << endl;
416  }
417 
418  if (wantSnap)
419  {
420  cpuTime timer;
421 
422  autoSnapDriver snapDriver
423  (
424  meshRefiner,
425  globalToPatch
426  );
427 
428  // Snap parameters
429  snapParameters snapParams(snapDict);
430 
431  if (!overwrite)
432  {
433  const_cast<Time&>(mesh.time())++;
434  }
435 
436  snapDriver.doSnap(snapDict, motionDict, snapParams);
437 
438  writeMesh
439  (
440  "Snapped mesh",
441  meshRefiner,
442  debug
443  );
444 
445  Info<< "Mesh snapped in = "
446  << timer.cpuTimeIncrement() << " s." << endl;
447  }
448 
449  if (wantLayers)
450  {
451  cpuTime timer;
452 
453  autoLayerDriver layerDriver(meshRefiner);
454 
455  // Layer addition parameters
456  layerParameters layerParams(layerDict, mesh.boundaryMesh());
457 
459  bool preBalance;
460  {
461  refinementParameters refineParams(refineDict);
462 
463  preBalance = returnReduce
464  (
465  (mesh.nCells() >= refineParams.maxLocalCells()),
466  orOp<bool>()
467  );
468  }
469 
470 
471  if (!overwrite)
472  {
473  const_cast<Time&>(mesh.time())++;
474  }
475 
476  layerDriver.doLayers
477  (
478  layerDict,
479  motionDict,
480  layerParams,
481  preBalance,
482  decomposer,
483  distributor
484  );
485 
486  writeMesh
487  (
488  "Layer mesh",
489  meshRefiner,
490  debug
491  );
492 
493  Info<< "Layers added in = "
494  << timer.cpuTimeIncrement() << " s." << endl;
495  }
496 
497 
498  Info<< "Finished meshing in = "
499  << runTime.elapsedCpuTime() << " s." << endl;
500 
501  Info<< "End\n" << endl;
502 
503  return(0);
504 }
505 
506 
507 // ************************ vim: set sw=4 sts=4 et: ************************ //