FreeFOAM The Cross-Platform CFD Toolkit
polyDualMeshApp.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  polyDualMesh
26 
27 Description
28  Calculate the dual of a polyMesh. Adheres to all the feature and patch
29  edges.
30 
31  Detects any boundary edge > angle and creates multiple boundary faces
32  for it. Normal behaviour is to have each point become a cell
33  (1.5 behaviour)
34 
35 Usage
36 
37  - polyDualMesh [OPTIONS] <feature angle [0-180]>
38 
39  @param -concaveMultiCells
40  Creates multiple cells for each point on a concave edge. Might limit
41  the amount of distortion on some meshes.
42 
43  @param -splitAllFaces
44  Normally only constructs a single face between two cells. This single face
45  might be too distorted. splitAllFaces will create a single face for every
46  original cell the face passes through. The mesh will thus have
47  multiple faces inbetween two cells! (so is not strictly upper-triangular
48  anymore - checkMesh will complain)
49 
50  @param -doNotPreserveFaceZones:
51  By default all faceZones are preserved by marking all faces, edges and
52  points on them as features. The -doNotPreserveFaceZones disables this
53  behaviour.
54 
55  @param -overwrite \n
56  Overwrite existing data.
57 
58  @param -case <dir>\n
59  Case directory.
60 
61  @param -help \n
62  Display help message.
63 
64  @param -doc \n
65  Display Doxygen API documentation page for this application.
66 
67  @param -srcDoc \n
68  Display Doxygen source documentation page for this application.
69 
70 Note
71  This is just a driver for meshDualiser. Substitute your own
72  simpleMarkFeatures to have different behaviour.
73 
74 \*---------------------------------------------------------------------------*/
75 
76 #include <OpenFOAM/argList.H>
77 #include <OpenFOAM/Time.H>
78 #include <OpenFOAM/timeSelector.H>
79 #include <finiteVolume/fvMesh.H>
82 #include <OpenFOAM/mapPolyMesh.H>
84 #include <meshTools/meshTools.H>
85 #include <OpenFOAM/OFstream.H>
86 #include "meshDualiser.H"
87 #include <OpenFOAM/ReadFields.H>
88 #include <finiteVolume/volFields.H>
90 
91 using namespace Foam;
92 
93 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
94 
95 // Naive feature detection. All boundary edges with angle > featureAngle become
96 // feature edges. All points on feature edges become feature points. All
97 // boundary faces become feature faces.
98 void simpleMarkFeatures
99 (
100  const polyMesh& mesh,
101  const PackedBoolList& isBoundaryEdge,
102  const scalar featureAngle,
103  const bool concaveMultiCells,
104  const bool doNotPreserveFaceZones,
105 
106  labelList& featureFaces,
107  labelList& featureEdges,
108  labelList& singleCellFeaturePoints,
109  labelList& multiCellFeaturePoints
110 )
111 {
112  scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
113 
114  const polyBoundaryMesh& patches = mesh.boundaryMesh();
115 
116  // Working sets
117  labelHashSet featureEdgeSet;
118  labelHashSet singleCellFeaturePointSet;
119  labelHashSet multiCellFeaturePointSet;
120 
121 
122  // 1. Mark all edges between patches
123  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124 
125  forAll(patches, patchI)
126  {
127  const polyPatch& pp = patches[patchI];
128  const labelList& meshEdges = pp.meshEdges();
129 
130  // All patch corner edges. These need to be feature points & edges!
131  for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
132  {
133  label meshEdgeI = meshEdges[edgeI];
134  featureEdgeSet.insert(meshEdgeI);
135  singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]);
136  singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]);
137  }
138  }
139 
140 
141 
142  // 2. Mark all geometric feature edges
143  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144  // Make distinction between convex features where the boundary point becomes
145  // a single cell and concave features where the boundary point becomes
146  // multiple 'half' cells.
147 
148  // Addressing for all outside faces
149  primitivePatch allBoundary
150  (
152  (
153  mesh.faces(),
154  mesh.nFaces()-mesh.nInternalFaces(),
155  mesh.nInternalFaces()
156  ),
157  mesh.points()
158  );
159 
160  // Check for non-manifold points (surface pinched at point)
161  allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);
162 
163  // Check for non-manifold edges (surface pinched at edge)
164  const labelListList& edgeFaces = allBoundary.edgeFaces();
165  const labelList& meshPoints = allBoundary.meshPoints();
166 
167  forAll(edgeFaces, edgeI)
168  {
169  const labelList& eFaces = edgeFaces[edgeI];
170 
171  if (eFaces.size() > 2)
172  {
173  const edge& e = allBoundary.edges()[edgeI];
174 
175  //Info<< "Detected non-manifold boundary edge:" << edgeI
176  // << " coords:"
177  // << allBoundary.points()[meshPoints[e[0]]]
178  // << allBoundary.points()[meshPoints[e[1]]] << endl;
179 
180  singleCellFeaturePointSet.insert(meshPoints[e[0]]);
181  singleCellFeaturePointSet.insert(meshPoints[e[1]]);
182  }
183  }
184 
185  // Check for features.
186  forAll(edgeFaces, edgeI)
187  {
188  const labelList& eFaces = edgeFaces[edgeI];
189 
190  if (eFaces.size() == 2)
191  {
192  label f0 = eFaces[0];
193  label f1 = eFaces[1];
194 
195  // check angle
196  const vector& n0 = allBoundary.faceNormals()[f0];
197  const vector& n1 = allBoundary.faceNormals()[f1];
198 
199  if ((n0 & n1) < minCos)
200  {
201  const edge& e = allBoundary.edges()[edgeI];
202  label v0 = meshPoints[e[0]];
203  label v1 = meshPoints[e[1]];
204 
205  label meshEdgeI = meshTools::findEdge(mesh, v0, v1);
206  featureEdgeSet.insert(meshEdgeI);
207 
208  // Check if convex or concave by looking at angle
209  // between face centres and normal
210  vector c1c0
211  (
212  allBoundary[f1].centre(allBoundary.points())
213  - allBoundary[f0].centre(allBoundary.points())
214  );
215 
216  if (concaveMultiCells && (c1c0 & n0) > SMALL)
217  {
218  // Found concave edge. Make into multiCell features
219  Info<< "Detected concave feature edge:" << edgeI
220  << " cos:" << (c1c0 & n0)
221  << " coords:"
222  << allBoundary.points()[v0]
223  << allBoundary.points()[v1]
224  << endl;
225 
226  singleCellFeaturePointSet.erase(v0);
227  multiCellFeaturePointSet.insert(v0);
228  singleCellFeaturePointSet.erase(v1);
229  multiCellFeaturePointSet.insert(v1);
230  }
231  else
232  {
233  // Convex. singleCell feature.
234  if (!multiCellFeaturePointSet.found(v0))
235  {
236  singleCellFeaturePointSet.insert(v0);
237  }
238  if (!multiCellFeaturePointSet.found(v1))
239  {
240  singleCellFeaturePointSet.insert(v1);
241  }
242  }
243  }
244  }
245  }
246 
247 
248  // 3. Mark all feature faces
249  // ~~~~~~~~~~~~~~~~~~~~~~~~~
250 
251  // Face centres that need inclusion in the dual mesh
252  labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces());
253  // A. boundary faces.
254  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
255  {
256  featureFaceSet.insert(faceI);
257  }
258 
259  // B. face zones.
260  const faceZoneMesh& faceZones = mesh.faceZones();
261 
262  if (doNotPreserveFaceZones)
263  {
264  if (faceZones.size() > 0)
265  {
266  WarningIn("simpleMarkFeatures(..)")
267  << "Detected " << faceZones.size()
268  << " faceZones. These will not be preserved."
269  << endl;
270  }
271  }
272  else
273  {
274  if (faceZones.size() > 0)
275  {
276  Info<< "Detected " << faceZones.size()
277  << " faceZones. Preserving these by marking their"
278  << " points, edges and faces as features." << endl;
279  }
280 
281  forAll(faceZones, zoneI)
282  {
283  const faceZone& fz = faceZones[zoneI];
284 
285  Info<< "Inserting all faces in faceZone " << fz.name()
286  << " as features." << endl;
287 
288  forAll(fz, i)
289  {
290  label faceI = fz[i];
291  const face& f = mesh.faces()[faceI];
292  const labelList& fEdges = mesh.faceEdges()[faceI];
293 
294  featureFaceSet.insert(faceI);
295  forAll(f, fp)
296  {
297  // Mark point as multi cell point (since both sides of
298  // face should have different cells)
299  singleCellFeaturePointSet.erase(f[fp]);
300  multiCellFeaturePointSet.insert(f[fp]);
301 
302  // Make sure there are points on the edges.
303  featureEdgeSet.insert(fEdges[fp]);
304  }
305  }
306  }
307  }
308 
309  // Transfer to arguments
310  featureFaces = featureFaceSet.toc();
311  featureEdges = featureEdgeSet.toc();
312  singleCellFeaturePoints = singleCellFeaturePointSet.toc();
313  multiCellFeaturePoints = multiCellFeaturePointSet.toc();
314 }
315 
316 
317 // Dump features to .obj files
318 void dumpFeatures
319 (
320  const polyMesh& mesh,
321  const labelList& featureFaces,
322  const labelList& featureEdges,
323  const labelList& singleCellFeaturePoints,
324  const labelList& multiCellFeaturePoints
325 )
326 {
327  {
328  OFstream str("featureFaces.obj");
329  Info<< "Dumping centres of featureFaces to obj file " << str.name()
330  << endl;
331  forAll(featureFaces, i)
332  {
333  meshTools::writeOBJ(str, mesh.faceCentres()[featureFaces[i]]);
334  }
335  }
336  {
337  OFstream str("featureEdges.obj");
338  Info<< "Dumping featureEdges to obj file " << str.name() << endl;
339  label vertI = 0;
340 
341  forAll(featureEdges, i)
342  {
343  const edge& e = mesh.edges()[featureEdges[i]];
344  meshTools::writeOBJ(str, mesh.points()[e[0]]);
345  vertI++;
346  meshTools::writeOBJ(str, mesh.points()[e[1]]);
347  vertI++;
348  str<< "l " << vertI-1 << ' ' << vertI << nl;
349  }
350  }
351  {
352  OFstream str("singleCellFeaturePoints.obj");
353  Info<< "Dumping featurePoints that become a single cell to obj file "
354  << str.name() << endl;
355  forAll(singleCellFeaturePoints, i)
356  {
357  meshTools::writeOBJ(str, mesh.points()[singleCellFeaturePoints[i]]);
358  }
359  }
360  {
361  OFstream str("multiCellFeaturePoints.obj");
362  Info<< "Dumping featurePoints that become multiple cells to obj file "
363  << str.name() << endl;
364  forAll(multiCellFeaturePoints, i)
365  {
366  meshTools::writeOBJ(str, mesh.points()[multiCellFeaturePoints[i]]);
367  }
368  }
369 }
370 
371 
372 int main(int argc, char *argv[])
373 {
375  timeSelector::addOptions(true, false);
376 
377  argList::validArgs.append("feature angle[0-180]");
378  argList::validOptions.insert("splitAllFaces", "");
379  argList::validOptions.insert("concaveMultiCells", "");
380  argList::validOptions.insert("doNotPreserveFaceZones", "");
381  argList::validOptions.insert("overwrite", "");
382 
383 # include <OpenFOAM/setRootCase.H>
384 # include <OpenFOAM/createTime.H>
385 
387 
388 # include <OpenFOAM/createMesh.H>
389 
390  const word oldInstance = mesh.pointsInstance();
391 
392  // Mark boundary edges and points.
393  // (Note: in 1.4.2 we can use the built-in mesh point ordering
394  // facility instead)
395  PackedBoolList isBoundaryEdge(mesh.nEdges());
396  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
397  {
398  const labelList& fEdges = mesh.faceEdges()[faceI];
399 
400  forAll(fEdges, i)
401  {
402  isBoundaryEdge.set(fEdges[i], 1);
403  }
404  }
405 
406  scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
407 
408  scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
409 
410  Info<< "Feature:" << featureAngle << endl
411  << "minCos :" << minCos << endl
412  << endl;
413 
414 
415  const bool splitAllFaces = args.optionFound("splitAllFaces");
416  if (splitAllFaces)
417  {
418  Info<< "Splitting all internal faces to create multiple faces"
419  << " between two cells." << nl
420  << endl;
421  }
422 
423  const bool overwrite = args.optionFound("overwrite");
424  const bool doNotPreserveFaceZones = args.optionFound
425  (
426  "doNotPreserveFaceZones"
427  );
428  const bool concaveMultiCells = args.optionFound("concaveMultiCells");
429  if (concaveMultiCells)
430  {
431  Info<< "Generating multiple cells for points on concave feature edges."
432  << nl << endl;
433  }
434 
435 
436  // Face(centre)s that need inclusion in the dual mesh
437  labelList featureFaces;
438  // Edge(centre)s ,,
439  labelList featureEdges;
440  // Points (that become a single cell) that need inclusion in the dual mesh
441  labelList singleCellFeaturePoints;
442  // Points (that become a multiple cells) ,,
443  labelList multiCellFeaturePoints;
444 
445  // Sample implementation of feature detection.
446  simpleMarkFeatures
447  (
448  mesh,
449  isBoundaryEdge,
450  featureAngle,
451  concaveMultiCells,
452  doNotPreserveFaceZones,
453 
454  featureFaces,
455  featureEdges,
456  singleCellFeaturePoints,
457  multiCellFeaturePoints
458  );
459 
460  // If we want to split all polyMesh faces into one dualface per cell
461  // we are passing through we also need a point
462  // at the polyMesh facecentre and edgemid of the faces we want to
463  // split.
464  if (splitAllFaces)
465  {
466  featureEdges = identity(mesh.nEdges());
467  featureFaces = identity(mesh.nFaces());
468  }
469 
470  // Write obj files for debugging
471  dumpFeatures
472  (
473  mesh,
474  featureFaces,
475  featureEdges,
476  singleCellFeaturePoints,
477  multiCellFeaturePoints
478  );
479 
480 
481 
482  // Read objects in time directory
483  IOobjectList objects(mesh, runTime.timeName());
484 
485  // Read vol fields.
487  ReadFields(mesh, objects, vsFlds);
488 
490  ReadFields(mesh, objects, vvFlds);
491 
493  ReadFields(mesh, objects, vstFlds);
494 
495  PtrList<volSymmTensorField> vsymtFlds;
496  ReadFields(mesh, objects, vsymtFlds);
497 
499  ReadFields(mesh, objects, vtFlds);
500 
501  // Read surface fields.
503  ReadFields(mesh, objects, ssFlds);
504 
506  ReadFields(mesh, objects, svFlds);
507 
509  ReadFields(mesh, objects, sstFlds);
510 
512  ReadFields(mesh, objects, ssymtFlds);
513 
515  ReadFields(mesh, objects, stFlds);
516 
517 
518  // Topo change container
519  polyTopoChange meshMod(mesh.boundaryMesh().size());
520 
521  // Mesh dualiser engine
522  meshDualiser dualMaker(mesh);
523 
524  // Insert all commands into polyTopoChange to create dual of mesh. This does
525  // all the hard work.
526  dualMaker.setRefinement
527  (
528  splitAllFaces,
529  featureFaces,
530  featureEdges,
531  singleCellFeaturePoints,
532  multiCellFeaturePoints,
533  meshMod
534  );
535 
536  // Create mesh, return map from old to new mesh.
537  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
538 
539  // Update fields
540  mesh.updateMesh(map);
541 
542  // Optionally inflate mesh
543  if (map().hasMotionPoints())
544  {
545  mesh.movePoints(map().preMotionPoints());
546  }
547 
548  if (!overwrite)
549  {
550  runTime++;
551  }
552  else
553  {
554  mesh.setInstance(oldInstance);
555  }
556 
557  Info<< "Writing dual mesh to " << runTime.timeName() << endl;
558 
559  mesh.write();
560 
561  Info<< "End\n" << endl;
562 
563  return 0;
564 }
565 
566 
567 // ************************ vim: set sw=4 sts=4 et: ************************ //