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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2018-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 refineMesh
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 Utility to refine cells in multiple directions.
35
36 Command-line option handling:
37 - If -all specified or no refineMeshDict exists or, refine all cells
38 - If -dict <file> specified refine according to <file>
39 - If refineMeshDict exists refine according to refineMeshDict
40
41 When the refinement or all cells is selected apply 3D refinement for 3D
42 cases and 2D refinement for 2D cases.
43
44\*---------------------------------------------------------------------------*/
45
46#include "argList.H"
47#include "polyMesh.H"
48#include "Time.H"
49#include "cellSet.H"
50#include "multiDirRefinement.H"
51#include "labelIOList.H"
52#include "IOdictionary.H"
53#include "syncTools.H"
54
55using namespace Foam;
56
57// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58
59// Max cos angle for edges to be considered aligned with axis.
60static const scalar edgeTol = 1e-3;
61
62
63// Print edge statistics on mesh.
64void printEdgeStats(const polyMesh& mesh)
65{
66 label nX = 0;
67 label nY = 0;
68 label nZ = 0;
69
70 scalar minX = GREAT;
71 scalar maxX = -GREAT;
72 static const vector x(1, 0, 0);
73
74 scalar minY = GREAT;
75 scalar maxY = -GREAT;
76 static const vector y(0, 1, 0);
77
78 scalar minZ = GREAT;
79 scalar maxZ = -GREAT;
80 static const vector z(0, 0, 1);
81
82 scalar minOther = GREAT;
83 scalar maxOther = -GREAT;
84
85 bitSet isMasterEdge(syncTools::getMasterEdges(mesh));
86
87 const edgeList& edges = mesh.edges();
88
89 forAll(edges, edgeI)
90 {
91 if (isMasterEdge.test(edgeI))
92 {
93 const edge& e = edges[edgeI];
94
95 vector eVec(e.vec(mesh.points()));
96
97 scalar eMag = mag(eVec);
98
99 eVec /= eMag;
100
101 if (mag(eVec & x) > 1-edgeTol)
102 {
103 minX = min(minX, eMag);
104 maxX = max(maxX, eMag);
105 nX++;
106 }
107 else if (mag(eVec & y) > 1-edgeTol)
108 {
109 minY = min(minY, eMag);
110 maxY = max(maxY, eMag);
111 nY++;
112 }
113 else if (mag(eVec & z) > 1-edgeTol)
114 {
115 minZ = min(minZ, eMag);
116 maxZ = max(maxZ, eMag);
117 nZ++;
118 }
119 else
120 {
121 minOther = min(minOther, eMag);
122 maxOther = max(maxOther, eMag);
123 }
124 }
125 }
126
127 label nEdges = mesh.nEdges();
129 reduce(nX, sumOp<label>());
130 reduce(nY, sumOp<label>());
131 reduce(nZ, sumOp<label>());
132
133 reduce(minX, minOp<scalar>());
134 reduce(maxX, maxOp<scalar>());
135
136 reduce(minY, minOp<scalar>());
137 reduce(maxY, maxOp<scalar>());
138
139 reduce(minZ, minOp<scalar>());
140 reduce(maxZ, maxOp<scalar>());
141
142 reduce(minOther, minOp<scalar>());
143 reduce(maxOther, maxOp<scalar>());
144
145
146 Info<< "Mesh edge statistics:" << nl
147 << " x aligned : number:" << nX << "\tminLen:" << minX
148 << "\tmaxLen:" << maxX << nl
149 << " y aligned : number:" << nY << "\tminLen:" << minY
150 << "\tmaxLen:" << maxY << nl
151 << " z aligned : number:" << nZ << "\tminLen:" << minZ
152 << "\tmaxLen:" << maxZ << nl
153 << " other : number:" << nEdges - nX - nY - nZ
154 << "\tminLen:" << minOther
155 << "\tmaxLen:" << maxOther << nl << endl;
156}
157
158
159int main(int argc, char *argv[])
160{
161 argList::addNote
162 (
163 "Refine cells in multiple directions"
164 );
165
166 #include "addOverwriteOption.H"
167 #include "addRegionOption.H"
168
169 argList::addOption("dict", "file", "Alternative refineMeshDict");
170
171 argList::addBoolOption
172 (
173 "all",
174 "Refine all cells"
175 );
176
177 argList::noFunctionObjects(); // Never use function objects
178
179 #include "setRootCase.H"
180 #include "createTime.H"
181 #include "createNamedPolyMesh.H"
182
183 const word oldInstance = mesh.pointsInstance();
184
185 printEdgeStats(mesh);
186
187 //
188 // Read/construct control dictionary
189 //
190
191 const bool refineAllCells = args.found("all");
192 const bool overwrite = args.found("overwrite");
193
194 // List of cells to refine
195 labelList refCells;
196
197 // Dictionary to control refinement
198 dictionary refineDict;
199
200 // The -all option has precedence over -dict, or anything else
201 if (!refineAllCells)
202 {
203 const word dictName("refineMeshDict");
204
205 // Obtain dictPath here for messages
206 fileName dictPath = args.getOrDefault<fileName>("dict", "");
207
208 IOobject dictIO = IOobject::selectIO
209 (
211 (
212 dictName,
213 runTime.system(),
214 mesh,
215 IOobject::MUST_READ
216 ),
217 dictPath
218 );
219
220 // The reported dictPath will not be full resolved for the output
221 // (it will just be the -dict value) but this is purely cosmetic.
222
223 if (dictIO.typeHeaderOk<IOdictionary>(true))
224 {
225 refineDict = IOdictionary(dictIO);
226
227 Info<< "Refining according to ";
228
229 if (dictPath.empty())
230 {
231 Info<< dictName;
232 }
233 else
234 {
235 Info<< dictPath;
236 }
237 Info<< nl << endl;
238 }
239 else if (dictPath.empty())
240 {
241 Info<< "Refinement dictionary " << dictName << " not found" << nl;
242 }
243 else
244 {
246 << "Cannot open specified refinement dictionary "
247 << dictPath
248 << exit(FatalError);
249 }
250 }
251
252 if (refineDict.size())
253 {
254 const word setName(refineDict.get<word>("set"));
255
256 cellSet cells(mesh, setName);
257
258 Info<< "Read " << returnReduce(cells.size(), sumOp<label>())
259 << " cells from cellSet "
260 << cells.instance()/cells.local()/cells.name()
261 << endl << endl;
262
263 refCells = cells.toc();
264 }
265 else
266 {
267 Info<< "Refining all cells" << nl << endl;
268
269 // Select all cells
270 refCells = identity(mesh.nCells());
271
272 if (mesh.nGeometricD() == 3)
273 {
274 Info<< "3D case; refining all directions" << nl << endl;
275
277 directions[0] = "tan1";
278 directions[1] = "tan2";
279 directions[2] = "normal";
280 refineDict.add("directions", directions);
281
282 // Use hex cutter
283 refineDict.add("useHexTopology", "true");
284 }
285 else
286 {
287 const Vector<label> dirs(mesh.geometricD());
288
290
291 if (dirs.x() == -1)
292 {
293 Info<< "2D case; refining in directions y,z\n" << endl;
294 directions[0] = "tan2";
295 directions[1] = "normal";
296 }
297 else if (dirs.y() == -1)
298 {
299 Info<< "2D case; refining in directions x,z\n" << endl;
300 directions[0] = "tan1";
301 directions[1] = "normal";
302 }
303 else
304 {
305 Info<< "2D case; refining in directions x,y\n" << endl;
306 directions[0] = "tan1";
307 directions[1] = "tan2";
308 }
309
310 refineDict.add("directions", directions);
311
312 // Use standard cutter
313 refineDict.add("useHexTopology", "false");
314 }
315
316 refineDict.add("coordinateSystem", "global");
317
318 dictionary coeffsDict;
319 coeffsDict.add("tan1", vector(1, 0, 0));
320 coeffsDict.add("tan2", vector(0, 1, 0));
321 refineDict.add("globalCoeffs", coeffsDict);
322
323 refineDict.add("geometricCut", "false");
324 refineDict.add("writeMesh", "false");
325 }
326
327
328 string oldTimeName(runTime.timeName());
329
330 if (!overwrite)
331 {
332 ++runTime;
333 }
334
335
336 // Multi-directional refinement (does multiple iterations)
337 multiDirRefinement multiRef(mesh, refCells, refineDict);
338
339
340 // Write resulting mesh
341 if (overwrite)
342 {
343 mesh.setInstance(oldInstance);
344 }
345 mesh.write();
346
347
348 // Get list of cell splits.
349 // (is for every cell in old mesh the cells they have been split into)
350 const labelListList& oldToNew = multiRef.addedCells();
351
352
353 // Create cellSet with added cells for easy inspection
354 cellSet newCells(mesh, "refinedCells", refCells.size());
355
356 for (const labelList& added : oldToNew)
357 {
358 newCells.insert(added);
359 }
360
361 Info<< "Writing refined cells ("
362 << returnReduce(newCells.size(), sumOp<label>())
363 << ") to cellSet "
364 << newCells.instance()/newCells.local()/newCells.name()
365 << endl << endl;
366
367 newCells.write();
368
369
370 //
371 // Invert cell split to construct map from new to old
372 //
373
374 labelIOList newToOld
375 (
377 (
378 "cellMap",
379 runTime.timeName(),
380 polyMesh::meshSubDir,
381 mesh,
382 IOobject::NO_READ,
383 IOobject::AUTO_WRITE
384 ),
385 mesh.nCells()
386 );
387 newToOld.note() =
388 "From cells in mesh at "
389 + runTime.timeName()
390 + " to cells in mesh at "
391 + oldTimeName;
392
393
394 forAll(oldToNew, oldCelli)
395 {
396 const labelList& added = oldToNew[oldCelli];
397
398 if (added.size())
399 {
400 for (const label celli : added)
401 {
402 newToOld[celli] = oldCelli;
403 }
404 }
405 else
406 {
407 // Unrefined cell
408 newToOld[oldCelli] = oldCelli;
409 }
410 }
411
412 Info<< "Writing map from new to old cell to "
413 << newToOld.objectPath() << nl << endl;
414
415 newToOld.write();
416
417 printEdgeStats(mesh);
418
419 Info<< "End\n" << endl;
420
421 return 0;
422}
423
424
425// ************************************************************************* //
scalar y
Y[inertIndex] max(0.0)
reduce(hasMovingMesh, orOp< bool >())
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:65
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
A collection of cell labels.
Definition: cellSet.H:54
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
entry * add(entry *entryPtr, bool mergeEntry=false)
Add a new entry.
Definition: dictionary.C:640
Set of directions for each cell in the mesh. Either uniform and size=1 or one set of directions per c...
Definition: directions.H:70
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A class for handling file names.
Definition: fileName.H:76
Does multiple pass refinement to refine cells in multiple directions.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const word dictName("faMeshDefinition")
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const cellShapeList & cells
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
IOobject dictIO
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333