refineMesh.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) 2011-2016 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 Application
25  refineMesh
26 
27 Description
28  Utility to refine cells in multiple directions.
29 
30  Command-line option handling:
31  - If -all specified or no refineMeshDict exists or, refine all cells
32  - If -dict <file> specified refine according to <file>
33  - If refineMeshDict exists refine according to refineMeshDict
34 
35  When the refinement or all cells is selected apply 3D refinement for 3D
36  cases and 2D refinement for 2D cases.
37 
38 \*---------------------------------------------------------------------------*/
39 
40 #include "argList.H"
41 #include "polyMesh.H"
42 #include "Time.H"
43 #include "undoableMeshCutter.H"
44 #include "hexCellLooper.H"
45 #include "cellSet.H"
46 #include "twoDPointCorrector.H"
47 #include "directions.H"
48 #include "OFstream.H"
49 #include "multiDirRefinement.H"
50 #include "labelIOList.H"
51 #include "wedgePolyPatch.H"
52 #include "plane.H"
53 #include "SubField.H"
54 
55 using namespace Foam;
56 
57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58 
59 // Max cos angle for edges to be considered aligned with axis.
60 static const scalar edgeTol = 1e-3;
61 
62 
63 // Print edge statistics on mesh.
64 void printEdgeStats(const primitiveMesh& mesh)
65 {
66  label nX = 0;
67  label nY = 0;
68  label nZ = 0;
69 
70  scalar minX = GREAT;
71  scalar maxX = -GREAT;
72  vector x(1, 0, 0);
73 
74  scalar minY = GREAT;
75  scalar maxY = -GREAT;
76  vector y(0, 1, 0);
77 
78  scalar minZ = GREAT;
79  scalar maxZ = -GREAT;
80  vector z(0, 0, 1);
81 
82  scalar minOther = GREAT;
83  scalar maxOther = -GREAT;
84 
85  const edgeList& edges = mesh.edges();
86 
87  forAll(edges, edgeI)
88  {
89  const edge& e = edges[edgeI];
90 
91  vector eVec(e.vec(mesh.points()));
92 
93  scalar eMag = mag(eVec);
94 
95  eVec /= eMag;
96 
97  if (mag(eVec & x) > 1-edgeTol)
98  {
99  minX = min(minX, eMag);
100  maxX = max(maxX, eMag);
101  nX++;
102  }
103  else if (mag(eVec & y) > 1-edgeTol)
104  {
105  minY = min(minY, eMag);
106  maxY = max(maxY, eMag);
107  nY++;
108  }
109  else if (mag(eVec & z) > 1-edgeTol)
110  {
111  minZ = min(minZ, eMag);
112  maxZ = max(maxZ, eMag);
113  nZ++;
114  }
115  else
116  {
117  minOther = min(minOther, eMag);
118  maxOther = max(maxOther, eMag);
119  }
120  }
121 
122  Pout<< "Mesh edge statistics:" << endl
123  << " x aligned : number:" << nX << "\tminLen:" << minX
124  << "\tmaxLen:" << maxX << endl
125  << " y aligned : number:" << nY << "\tminLen:" << minY
126  << "\tmaxLen:" << maxY << endl
127  << " z aligned : number:" << nZ << "\tminLen:" << minZ
128  << "\tmaxLen:" << maxZ << endl
129  << " other : number:" << mesh.nEdges() - nX - nY - nZ
130  << "\tminLen:" << minOther
131  << "\tmaxLen:" << maxOther << endl << endl;
132 }
133 
134 
135 int main(int argc, char *argv[])
136 {
138  (
139  "refine cells in multiple directions"
140  );
141 
142  #include "addOverwriteOption.H"
143  #include "addRegionOption.H"
144  #include "addDictOption.H"
145 
147  (
148  "all",
149  "Refine all cells"
150  );
151 
152  #include "setRootCase.H"
153  #include "createTime.H"
154  runTime.functionObjects().off();
155  #include "createNamedPolyMesh.H"
156  const word oldInstance = mesh.pointsInstance();
157 
158  printEdgeStats(mesh);
159 
160  //
161  // Read/construct control dictionary
162  //
163 
164  const bool readDict = args.optionFound("dict");
165  const bool refineAllCells = args.optionFound("all");
166  const bool overwrite = args.optionFound("overwrite");
167 
168  // List of cells to refine
169  labelList refCells;
170 
171  // Dictionary to control refinement
172  dictionary refineDict;
173  const word dictName("refineMeshDict");
174 
175  if (readDict)
176  {
177  fileName dictPath = args["dict"];
178  if (isDir(dictPath))
179  {
180  dictPath = dictPath/dictName;
181  }
182 
184  (
185  dictPath,
186  mesh,
188  );
189 
190  if (!dictIO.headerOk())
191  {
193  << "Cannot open specified refinement dictionary "
194  << dictPath
195  << exit(FatalError);
196  }
197 
198  Info<< "Refining according to " << dictPath << nl << endl;
199 
200  refineDict = IOdictionary(dictIO);
201  }
202  else if (!refineAllCells)
203  {
205  (
206  dictName,
207  runTime.system(),
208  mesh,
210  );
211 
212  if (dictIO.headerOk())
213  {
214  Info<< "Refining according to " << dictName << nl << endl;
215 
216  refineDict = IOdictionary(dictIO);
217  }
218  else
219  {
220  Info<< "Refinement dictionary " << dictName << " not found" << endl;
221  }
222  }
223 
224  if (refineDict.size())
225  {
226  const word setName(refineDict.lookup("set"));
227 
228  cellSet cells(mesh, setName);
229 
230  Pout<< "Read " << cells.size() << " cells from cellSet "
231  << cells.instance()/cells.local()/cells.name()
232  << endl << endl;
233 
234  refCells = cells.toc();
235  }
236  else
237  {
238  Info<< "Refining all cells" << nl << endl;
239 
240  // Select all cells
241  refCells.setSize(mesh.nCells());
242 
243  forAll(mesh.cells(), celli)
244  {
245  refCells[celli] = celli;
246  }
247 
248  if (mesh.nGeometricD() == 3)
249  {
250  Info<< "3D case; refining all directions" << nl << endl;
251 
252  wordList directions(3);
253  directions[0] = "tan1";
254  directions[1] = "tan2";
255  directions[2] = "normal";
256  refineDict.add("directions", directions);
257 
258  // Use hex cutter
259  refineDict.add("useHexTopology", "true");
260  }
261  else
262  {
263  const Vector<label> dirs(mesh.geometricD());
264 
265  wordList directions(2);
266 
267  if (dirs.x() == -1)
268  {
269  Info<< "2D case; refining in directions y,z\n" << endl;
270  directions[0] = "tan2";
271  directions[1] = "normal";
272  }
273  else if (dirs.y() == -1)
274  {
275  Info<< "2D case; refining in directions x,z\n" << endl;
276  directions[0] = "tan1";
277  directions[1] = "normal";
278  }
279  else
280  {
281  Info<< "2D case; refining in directions x,y\n" << endl;
282  directions[0] = "tan1";
283  directions[1] = "tan2";
284  }
285 
286  refineDict.add("directions", directions);
287 
288  // Use standard cutter
289  refineDict.add("useHexTopology", "false");
290  }
291 
292  refineDict.add("coordinateSystem", "global");
293 
294  dictionary coeffsDict;
295  coeffsDict.add("tan1", vector(1, 0, 0));
296  coeffsDict.add("tan2", vector(0, 1, 0));
297  refineDict.add("globalCoeffs", coeffsDict);
298 
299  refineDict.add("geometricCut", "false");
300  refineDict.add("writeMesh", "false");
301  }
302 
303 
304  string oldTimeName(runTime.timeName());
305 
306  if (!overwrite)
307  {
308  runTime++;
309  }
310 
311 
312  // Multi-directional refinement (does multiple iterations)
313  multiDirRefinement multiRef(mesh, refCells, refineDict);
314 
315 
316  // Write resulting mesh
317  if (overwrite)
318  {
319  mesh.setInstance(oldInstance);
320  }
321  mesh.write();
322 
323 
324  // Get list of cell splits.
325  // (is for every cell in old mesh the cells they have been split into)
326  const labelListList& oldToNew = multiRef.addedCells();
327 
328 
329  // Create cellSet with added cells for easy inspection
330  cellSet newCells(mesh, "refinedCells", refCells.size());
331 
332  forAll(oldToNew, oldCelli)
333  {
334  const labelList& added = oldToNew[oldCelli];
335 
336  forAll(added, i)
337  {
338  newCells.insert(added[i]);
339  }
340  }
341 
342  Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
343  << newCells.instance()/newCells.local()/newCells.name()
344  << endl << endl;
345 
346  newCells.write();
347 
348 
349 
350 
351  //
352  // Invert cell split to construct map from new to old
353  //
354 
355  labelIOList newToOld
356  (
357  IOobject
358  (
359  "cellMap",
360  runTime.timeName(),
362  mesh,
365  ),
366  mesh.nCells()
367  );
368  newToOld.note() =
369  "From cells in mesh at "
370  + runTime.timeName()
371  + " to cells in mesh at "
372  + oldTimeName;
373 
374 
375  forAll(oldToNew, oldCelli)
376  {
377  const labelList& added = oldToNew[oldCelli];
378 
379  if (added.size())
380  {
381  forAll(added, i)
382  {
383  newToOld[added[i]] = oldCelli;
384  }
385  }
386  else
387  {
388  // Unrefined cell
389  newToOld[oldCelli] = oldCelli;
390  }
391  }
392 
393  Info<< "Writing map from new to old cell to "
394  << newToOld.objectPath() << nl << endl;
395 
396  newToOld.write();
397 
398  printEdgeStats(mesh);
399 
400  Info<< "End\n" << endl;
401 
402  return 0;
403 }
404 
405 
406 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 class for handling file names.
Definition: fileName.H:69
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:309
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Does multiple pass refinement to refine cells in multiple directions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Set of directions for each cell in the mesh. Either uniform and size=1 or one set of directions per c...
Definition: directions.H:64
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual const pointField & points() const =0
Return mesh points.
bool isDir(const fileName &)
Does the name exist as a DIRECTORY in the file system?
Definition: POSIX.C:486
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
const cellList & cells() const
bool add(entry *, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:737
label nCells() const
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
string & note()
Return non-constant access to the optional note.
Definition: IOobject.H:272
scalar y
dynamicFvMesh & mesh
const cellShapeList & cells
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A class for handling words, derived from string.
Definition: word.H:59
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
static const char nl
Definition: Ostream.H:262
label nEdges() const
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
IOobject dictIO(dictName, runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE)
void setSize(const label)
Reset size of List.
Definition: List.C:295
word dictName("noiseDict")
A collection of cell labels.
Definition: cellSet.H:48
vector vec(const pointField &) const
Return the vector (end - start)
Definition: edgeI.H:157
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
bool headerOk()
Read and check header info.
Definition: IOobject.C:400
virtual Ostream & write(const token &)=0
Write next token to stream.
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
Definition: argList.C:83
static void addNote(const string &)
Add extra notes for the usage information.
Definition: argList.C:124
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451