FreeFOAM The Cross-Platform CFD Toolkit
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) 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  refineMesh
26 
27 Description
28  Utility to refine cells in multiple directions.
29 
30  Either supply -all option to refine all cells (3D refinement for 3D
31  cases; 2D for 2D cases) or reads a refineMeshDict with
32  - cellSet to refine
33  - directions to refine
34 
35 Usage
36 
37  - refineMesh [OPTIONS]
38 
39  @param -dict <specify refineMeshDict>\n
40  Refine according to refineMeshDict.
41 
42  @param -overwrite \n
43  Overwrite existing data.
44 
45  @param -case <dir>\n
46  Case directory.
47 
48  @param -parallel \n
49  Run in parallel.
50 
51  @param -help \n
52  Display help message.
53 
54  @param -doc \n
55  Display Doxygen API documentation page for this application.
56 
57  @param -srcDoc \n
58  Display Doxygen source documentation page for this application.
59 
60 \*---------------------------------------------------------------------------*/
61 
62 #include <OpenFOAM/argList.H>
63 #include <OpenFOAM/polyMesh.H>
64 #include <OpenFOAM/Time.H>
67 #include <meshTools/cellSet.H>
69 #include <dynamicMesh/directions.H>
70 #include <OpenFOAM/OFstream.H>
72 #include <OpenFOAM/labelIOList.H>
74 #include <OpenFOAM/plane.H>
75 
76 using namespace Foam;
77 
78 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
79 
80 
81 // Max cos angle for edges to be considered aligned with axis.
82 static const scalar edgeTol = 1E-3;
83 
84 
85 // Calculate some edge statistics on mesh.
86 void printEdgeStats(const primitiveMesh& mesh)
87 {
88  label nX = 0;
89  label nY = 0;
90  label nZ = 0;
91 
92  scalar minX = GREAT;
93  scalar maxX = -GREAT;
94  vector x(1, 0, 0);
95 
96  scalar minY = GREAT;
97  scalar maxY = -GREAT;
98  vector y(0, 1, 0);
99 
100  scalar minZ = GREAT;
101  scalar maxZ = -GREAT;
102  vector z(0, 0, 1);
103 
104  scalar minOther = GREAT;
105  scalar maxOther = -GREAT;
106 
107  const edgeList& edges = mesh.edges();
108 
109  forAll(edges, edgeI)
110  {
111  const edge& e = edges[edgeI];
112 
113  vector eVec(e.vec(mesh.points()));
114 
115  scalar eMag = mag(eVec);
116 
117  eVec /= eMag;
118 
119  if (mag(eVec & x) > 1-edgeTol)
120  {
121  minX = min(minX, eMag);
122  maxX = max(maxX, eMag);
123  nX++;
124  }
125  else if (mag(eVec & y) > 1-edgeTol)
126  {
127  minY = min(minY, eMag);
128  maxY = max(maxY, eMag);
129  nY++;
130  }
131  else if (mag(eVec & z) > 1-edgeTol)
132  {
133  minZ = min(minZ, eMag);
134  maxZ = max(maxZ, eMag);
135  nZ++;
136  }
137  else
138  {
139  minOther = min(minOther, eMag);
140  maxOther = max(maxOther, eMag);
141  }
142  }
143 
144  Pout<< "Mesh edge statistics:" << endl
145  << " x aligned : number:" << nX << "\tminLen:" << minX
146  << "\tmaxLen:" << maxX << endl
147  << " y aligned : number:" << nY << "\tminLen:" << minY
148  << "\tmaxLen:" << maxY << endl
149  << " z aligned : number:" << nZ << "\tminLen:" << minZ
150  << "\tmaxLen:" << maxZ << endl
151  << " other : number:" << mesh.nEdges() - nX - nY - nZ
152  << "\tminLen:" << minOther
153  << "\tmaxLen:" << maxOther << endl << endl;
154 }
155 
156 
157 // Return index of coordinate axis.
158 label axis(const vector& normal)
159 {
160  label axisIndex = -1;
161 
162  if (mag(normal & point(1, 0, 0)) > (1-edgeTol))
163  {
164  axisIndex = 0;
165  }
166  else if (mag(normal & point(0, 1, 0)) > (1-edgeTol))
167  {
168  axisIndex = 1;
169  }
170  else if (mag(normal & point(0, 0, 1)) > (1-edgeTol))
171  {
172  axisIndex = 2;
173  }
174 
175  return axisIndex;
176 }
177 
178 
179 //- Returns -1 or cartesian coordinate component (0=x, 1=y, 2=z) of normal
180 // in case of 2D mesh
181 label twoDNess(const polyMesh& mesh)
182 {
183  const pointField& ctrs = mesh.cellCentres();
184 
185  if (ctrs.size() < 2)
186  {
187  return -1;
188  }
189 
190 
191  //
192  // 1. All cell centres on single plane aligned with x, y or z
193  //
194 
195  // Determine 3 points to base plane on.
196  vector vec10 = ctrs[1] - ctrs[0];
197  vec10 /= mag(vec10);
198 
199  label otherCellI = -1;
200 
201  for (label cellI = 2; cellI < ctrs.size(); cellI++)
202  {
203  vector vec(ctrs[cellI] - ctrs[0]);
204  vec /= mag(vec);
205 
206  if (mag(vec & vec10) < 0.9)
207  {
208  // ctrs[cellI] not in line with n
209  otherCellI = cellI;
210 
211  break;
212  }
213  }
214 
215  if (otherCellI == -1)
216  {
217  // Cannot find cell to make decent angle with cell0-cell1 vector.
218  // Note: what to do here? All cells (almost) in one line. Maybe 1D case?
219  return -1;
220  }
221 
222  plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);
223 
224 
225  forAll(ctrs, cellI)
226  {
227  const labelList& cEdges = mesh.cellEdges()[cellI];
228 
229  scalar minLen = GREAT;
230 
231  forAll(cEdges, i)
232  {
233  minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
234  }
235 
236  if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
237  {
238  // Centres not in plane
239  return -1;
240  }
241  }
242 
243  label axisIndex = axis(cellPlane.normal());
244 
245  if (axisIndex == -1)
246  {
247  return axisIndex;
248  }
249 
250 
251  const polyBoundaryMesh& patches = mesh.boundaryMesh();
252 
253 
254  //
255  // 2. No edges without points on boundary
256  //
257 
258  // Mark boundary points
259  boolList boundaryPoint(mesh.points().size(), false);
260 
261  forAll(patches, patchI)
262  {
263  const polyPatch& patch = patches[patchI];
264 
265  forAll(patch, patchFaceI)
266  {
267  const face& f = patch[patchFaceI];
268 
269  forAll(f, fp)
270  {
271  boundaryPoint[f[fp]] = true;
272  }
273  }
274  }
275 
276 
277  const edgeList& edges = mesh.edges();
278 
279  forAll(edges, edgeI)
280  {
281  const edge& e = edges[edgeI];
282 
283  if (!boundaryPoint[e.start()] && !boundaryPoint[e.end()])
284  {
285  // Edge has no point on boundary.
286  return -1;
287  }
288  }
289 
290 
291  // 3. For all non-wedge patches: all faces either perp or aligned with
292  // cell-plane normal. (wedge patches already checked upon construction)
293 
294  forAll(patches, patchI)
295  {
296  const polyPatch& patch = patches[patchI];
297 
298  if (!isA<wedgePolyPatch>(patch))
299  {
300  const vectorField& n = patch.faceAreas();
301 
302  scalarField cosAngle = mag(n/mag(n) & cellPlane.normal());
303 
304  if (mag(min(cosAngle) - max(cosAngle)) > 1E-6)
305  {
306  // cosAngle should be either ~1 over all faces (2D front and
307  // back) or ~0 (all other patches perp to 2D)
308  return -1;
309  }
310  }
311  }
312 
313  return axisIndex;
314 }
315 
316 
317 // Main program:
318 
319 int main(int argc, char *argv[])
320 {
322  Foam::argList::validOptions.insert("overwrite", "");
323 
324 # include <OpenFOAM/setRootCase.H>
325 # include <OpenFOAM/createTime.H>
326  runTime.functionObjects().off();
327 # include <OpenFOAM/createPolyMesh.H>
328  const word oldInstance = mesh.pointsInstance();
329 
330  printEdgeStats(mesh);
331 
332 
333  //
334  // Read/construct control dictionary
335  //
336 
337  bool readDict = args.optionFound("dict");
338  bool overwrite = args.optionFound("overwrite");
339 
340  // List of cells to refine
341  labelList refCells;
342 
343  // Dictionary to control refinement
344  dictionary refineDict;
345 
346  if (readDict)
347  {
348  Info<< "Refining according to refineMeshDict" << nl << endl;
349 
350  refineDict =
352  (
353  IOobject
354  (
355  "refineMeshDict",
356  runTime.system(),
357  mesh,
360  )
361  );
362 
363  word setName(refineDict.lookup("set"));
364 
365  cellSet cells(mesh, setName);
366 
367  Pout<< "Read " << cells.size() << " cells from cellSet "
368  << cells.instance()/cells.local()/cells.name()
369  << endl << endl;
370 
371  refCells = cells.toc();
372  }
373  else
374  {
375  Info<< "Refining all cells" << nl << endl;
376 
377  // Select all cells
378  refCells.setSize(mesh.nCells());
379 
380  forAll(mesh.cells(), cellI)
381  {
382  refCells[cellI] = cellI;
383  }
384 
385 
386  // Set refinement directions based on 2D/3D
387  label axisIndex = twoDNess(mesh);
388 
389  if (axisIndex == -1)
390  {
391  Info<< "3D case; refining all directions" << nl << endl;
392 
393  wordList directions(3);
394  directions[0] = "tan1";
395  directions[1] = "tan2";
396  directions[2] = "normal";
397  refineDict.add("directions", directions);
398 
399  // Use hex cutter
400  refineDict.add("useHexTopology", "true");
401  }
402  else
403  {
404  wordList directions(2);
405 
406  if (axisIndex == 0)
407  {
408  Info<< "2D case; refining in directions y,z\n" << endl;
409  directions[0] = "tan2";
410  directions[1] = "normal";
411  }
412  else if (axisIndex == 1)
413  {
414  Info<< "2D case; refining in directions x,z\n" << endl;
415  directions[0] = "tan1";
416  directions[1] = "normal";
417  }
418  else
419  {
420  Info<< "2D case; refining in directions x,y\n" << endl;
421  directions[0] = "tan1";
422  directions[1] = "tan2";
423  }
424 
425  refineDict.add("directions", directions);
426 
427  // Use standard cutter
428  refineDict.add("useHexTopology", "false");
429  }
430 
431  refineDict.add("coordinateSystem", "global");
432 
433  dictionary coeffsDict;
434  coeffsDict.add("tan1", vector(1, 0, 0));
435  coeffsDict.add("tan2", vector(0, 1, 0));
436  refineDict.add("globalCoeffs", coeffsDict);
437 
438  refineDict.add("geometricCut", "false");
439  refineDict.add("writeMesh", "false");
440  }
441 
442 
443  string oldTimeName(runTime.timeName());
444 
445  if (!overwrite)
446  {
447  runTime++;
448  }
449 
450 
451  // Multi-directional refinement (does multiple iterations)
452  multiDirRefinement multiRef(mesh, refCells, refineDict);
453 
454 
455  // Write resulting mesh
456  if (overwrite)
457  {
458  mesh.setInstance(oldInstance);
459  }
460  mesh.write();
461 
462 
463  // Get list of cell splits.
464  // (is for every cell in old mesh the cells they have been split into)
465  const labelListList& oldToNew = multiRef.addedCells();
466 
467 
468  // Create cellSet with added cells for easy inspection
469  cellSet newCells(mesh, "refinedCells", refCells.size());
470 
471  forAll(oldToNew, oldCellI)
472  {
473  const labelList& added = oldToNew[oldCellI];
474 
475  forAll(added, i)
476  {
477  newCells.insert(added[i]);
478  }
479  }
480 
481  Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
482  << newCells.instance()/newCells.local()/newCells.name()
483  << endl << endl;
484 
485  newCells.write();
486 
487 
488 
489 
490  //
491  // Invert cell split to construct map from new to old
492  //
493 
494  labelIOList newToOld
495  (
496  IOobject
497  (
498  "cellMap",
499  runTime.timeName(),
501  mesh,
504  ),
505  mesh.nCells()
506  );
507  newToOld.note() =
508  "From cells in mesh at "
509  + runTime.timeName()
510  + " to cells in mesh at "
511  + oldTimeName;
512 
513 
514  forAll(oldToNew, oldCellI)
515  {
516  const labelList& added = oldToNew[oldCellI];
517 
518  if (added.size())
519  {
520  forAll(added, i)
521  {
522  newToOld[added[i]] = oldCellI;
523  }
524  }
525  else
526  {
527  // Unrefined cell
528  newToOld[oldCellI] = oldCellI;
529  }
530  }
531 
532  Info<< "Writing map from new to old cell to "
533  << newToOld.objectPath() << nl << endl;
534 
535  newToOld.write();
536 
537 
538  // Some statistics.
539 
540  printEdgeStats(mesh);
541 
542  Info<< "End\n" << endl;
543 
544  return 0;
545 }
546 
547 
548 // ************************ vim: set sw=4 sts=4 et: ************************ //