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) 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  polyDualMesh
26 
27 Description
28  Calculates the dual of a polyMesh. Adheres to all the feature and patch
29  edges.
30 
31 Usage
32  \b polyDualMesh featureAngle
33 
34  Detects any boundary edge > angle and creates multiple boundary faces
35  for it. Normal behaviour is to have each point become a cell
36  (1.5 behaviour)
37 
38  Options:
39  - \par -concaveMultiCells
40  Creates multiple cells for each point on a concave edge. Might limit
41  the amount of distortion on some meshes.
42 
43  - \par -splitAllFaces
44  Normally only constructs a single face between two cells. This single
45  face might be too distorted. splitAllFaces will create a single face for
46  every 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  - \par -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 Note
56  It is just a driver for meshDualiser. Substitute your own simpleMarkFeatures
57  to have different behaviour.
58 
59 \*---------------------------------------------------------------------------*/
60 
61 #include "argList.H"
62 #include "Time.H"
63 #include "fvMesh.H"
64 #include "unitConversion.H"
65 #include "polyTopoChange.H"
66 #include "mapPolyMesh.H"
67 #include "PackedBoolList.H"
68 #include "meshTools.H"
69 #include "OFstream.H"
70 #include "meshDualiser.H"
71 #include "ReadFields.H"
72 #include "volFields.H"
73 #include "surfaceFields.H"
74 
75 using namespace Foam;
76 
77 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
78 
79 // Naive feature detection. All boundary edges with angle > featureAngle become
80 // feature edges. All points on feature edges become feature points. All
81 // boundary faces become feature faces.
82 void simpleMarkFeatures
83 (
84  const polyMesh& mesh,
85  const PackedBoolList& isBoundaryEdge,
86  const scalar featureAngle,
87  const bool concaveMultiCells,
88  const bool doNotPreserveFaceZones,
89 
90  labelList& featureFaces,
91  labelList& featureEdges,
92  labelList& singleCellFeaturePoints,
93  labelList& multiCellFeaturePoints
94 )
95 {
96  scalar minCos = Foam::cos(degToRad(featureAngle));
97 
98  const polyBoundaryMesh& patches = mesh.boundaryMesh();
99 
100  // Working sets
101  labelHashSet featureEdgeSet;
102  labelHashSet singleCellFeaturePointSet;
103  labelHashSet multiCellFeaturePointSet;
104 
105 
106  // 1. Mark all edges between patches
107  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
108 
109  forAll(patches, patchi)
110  {
111  const polyPatch& pp = patches[patchi];
112  const labelList& meshEdges = pp.meshEdges();
113 
114  // All patch corner edges. These need to be feature points & edges!
115  for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
116  {
117  label meshEdgeI = meshEdges[edgeI];
118  featureEdgeSet.insert(meshEdgeI);
119  singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]);
120  singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]);
121  }
122  }
123 
124 
125 
126  // 2. Mark all geometric feature edges
127  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
128  // Make distinction between convex features where the boundary point becomes
129  // a single cell and concave features where the boundary point becomes
130  // multiple 'half' cells.
131 
132  // Addressing for all outside faces
133  primitivePatch allBoundary
134  (
136  (
137  mesh.faces(),
138  mesh.nFaces()-mesh.nInternalFaces(),
139  mesh.nInternalFaces()
140  ),
141  mesh.points()
142  );
143 
144  // Check for non-manifold points (surface pinched at point)
145  allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);
146 
147  // Check for non-manifold edges (surface pinched at edge)
148  const labelListList& edgeFaces = allBoundary.edgeFaces();
149  const labelList& meshPoints = allBoundary.meshPoints();
150 
151  forAll(edgeFaces, edgeI)
152  {
153  const labelList& eFaces = edgeFaces[edgeI];
154 
155  if (eFaces.size() > 2)
156  {
157  const edge& e = allBoundary.edges()[edgeI];
158 
159  //Info<< "Detected non-manifold boundary edge:" << edgeI
160  // << " coords:"
161  // << allBoundary.points()[meshPoints[e[0]]]
162  // << allBoundary.points()[meshPoints[e[1]]] << endl;
163 
164  singleCellFeaturePointSet.insert(meshPoints[e[0]]);
165  singleCellFeaturePointSet.insert(meshPoints[e[1]]);
166  }
167  }
168 
169  // Check for features.
170  forAll(edgeFaces, edgeI)
171  {
172  const labelList& eFaces = edgeFaces[edgeI];
173 
174  if (eFaces.size() == 2)
175  {
176  label f0 = eFaces[0];
177  label f1 = eFaces[1];
178 
179  // check angle
180  const vector& n0 = allBoundary.faceNormals()[f0];
181  const vector& n1 = allBoundary.faceNormals()[f1];
182 
183  if ((n0 & n1) < minCos)
184  {
185  const edge& e = allBoundary.edges()[edgeI];
186  label v0 = meshPoints[e[0]];
187  label v1 = meshPoints[e[1]];
188 
189  label meshEdgeI = meshTools::findEdge(mesh, v0, v1);
190  featureEdgeSet.insert(meshEdgeI);
191 
192  // Check if convex or concave by looking at angle
193  // between face centres and normal
194  vector c1c0
195  (
196  allBoundary[f1].centre(allBoundary.points())
197  - allBoundary[f0].centre(allBoundary.points())
198  );
199 
200  if (concaveMultiCells && (c1c0 & n0) > SMALL)
201  {
202  // Found concave edge. Make into multiCell features
203  Info<< "Detected concave feature edge:" << edgeI
204  << " cos:" << (c1c0 & n0)
205  << " coords:"
206  << allBoundary.points()[v0]
207  << allBoundary.points()[v1]
208  << endl;
209 
210  singleCellFeaturePointSet.erase(v0);
211  multiCellFeaturePointSet.insert(v0);
212  singleCellFeaturePointSet.erase(v1);
213  multiCellFeaturePointSet.insert(v1);
214  }
215  else
216  {
217  // Convex. singleCell feature.
218  if (!multiCellFeaturePointSet.found(v0))
219  {
220  singleCellFeaturePointSet.insert(v0);
221  }
222  if (!multiCellFeaturePointSet.found(v1))
223  {
224  singleCellFeaturePointSet.insert(v1);
225  }
226  }
227  }
228  }
229  }
230 
231 
232  // 3. Mark all feature faces
233  // ~~~~~~~~~~~~~~~~~~~~~~~~~
234 
235  // Face centres that need inclusion in the dual mesh
236  labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces());
237  // A. boundary faces.
238  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
239  {
240  featureFaceSet.insert(facei);
241  }
242 
243  // B. face zones.
244  const faceZoneMesh& faceZones = mesh.faceZones();
245 
246  if (doNotPreserveFaceZones)
247  {
248  if (faceZones.size() > 0)
249  {
251  << "Detected " << faceZones.size()
252  << " faceZones. These will not be preserved."
253  << endl;
254  }
255  }
256  else
257  {
258  if (faceZones.size() > 0)
259  {
260  Info<< "Detected " << faceZones.size()
261  << " faceZones. Preserving these by marking their"
262  << " points, edges and faces as features." << endl;
263  }
264 
265  forAll(faceZones, zoneI)
266  {
267  const faceZone& fz = faceZones[zoneI];
268 
269  Info<< "Inserting all faces in faceZone " << fz.name()
270  << " as features." << endl;
271 
272  forAll(fz, i)
273  {
274  label facei = fz[i];
275  const face& f = mesh.faces()[facei];
276  const labelList& fEdges = mesh.faceEdges()[facei];
277 
278  featureFaceSet.insert(facei);
279  forAll(f, fp)
280  {
281  // Mark point as multi cell point (since both sides of
282  // face should have different cells)
283  singleCellFeaturePointSet.erase(f[fp]);
284  multiCellFeaturePointSet.insert(f[fp]);
285 
286  // Make sure there are points on the edges.
287  featureEdgeSet.insert(fEdges[fp]);
288  }
289  }
290  }
291  }
292 
293  // Transfer to arguments
294  featureFaces = featureFaceSet.toc();
295  featureEdges = featureEdgeSet.toc();
296  singleCellFeaturePoints = singleCellFeaturePointSet.toc();
297  multiCellFeaturePoints = multiCellFeaturePointSet.toc();
298 }
299 
300 
301 // Dump features to .obj files
302 void dumpFeatures
303 (
304  const polyMesh& mesh,
305  const labelList& featureFaces,
306  const labelList& featureEdges,
307  const labelList& singleCellFeaturePoints,
308  const labelList& multiCellFeaturePoints
309 )
310 {
311  {
312  OFstream str("featureFaces.obj");
313  Info<< "Dumping centres of featureFaces to obj file " << str.name()
314  << endl;
315  forAll(featureFaces, i)
316  {
317  meshTools::writeOBJ(str, mesh.faceCentres()[featureFaces[i]]);
318  }
319  }
320  {
321  OFstream str("featureEdges.obj");
322  Info<< "Dumping featureEdges to obj file " << str.name() << endl;
323  label vertI = 0;
324 
325  forAll(featureEdges, i)
326  {
327  const edge& e = mesh.edges()[featureEdges[i]];
328  meshTools::writeOBJ(str, mesh.points()[e[0]]);
329  vertI++;
330  meshTools::writeOBJ(str, mesh.points()[e[1]]);
331  vertI++;
332  str<< "l " << vertI-1 << ' ' << vertI << nl;
333  }
334  }
335  {
336  OFstream str("singleCellFeaturePoints.obj");
337  Info<< "Dumping featurePoints that become a single cell to obj file "
338  << str.name() << endl;
339  forAll(singleCellFeaturePoints, i)
340  {
341  meshTools::writeOBJ(str, mesh.points()[singleCellFeaturePoints[i]]);
342  }
343  }
344  {
345  OFstream str("multiCellFeaturePoints.obj");
346  Info<< "Dumping featurePoints that become multiple cells to obj file "
347  << str.name() << endl;
348  forAll(multiCellFeaturePoints, i)
349  {
350  meshTools::writeOBJ(str, mesh.points()[multiCellFeaturePoints[i]]);
351  }
352  }
353 }
354 
355 
356 int main(int argc, char *argv[])
357 {
358  #include "addOverwriteOption.H"
360 
361  argList::validArgs.append("featureAngle [0-180]");
363  (
364  "splitAllFaces",
365  "have multiple faces inbetween cells"
366  );
368  (
369  "concaveMultiCells",
370  "split cells on concave boundary edges into multiple cells"
371  );
373  (
374  "doNotPreserveFaceZones",
375  "disable the default behaviour of preserving faceZones by having"
376  " multiple faces inbetween cells"
377  );
378 
379  #include "setRootCase.H"
380  #include "createTime.H"
381  #include "createMesh.H"
382 
383  const word oldInstance = mesh.pointsInstance();
384 
385  // Mark boundary edges and points.
386  // (Note: in 1.4.2 we can use the built-in mesh point ordering
387  // facility instead)
388  PackedBoolList isBoundaryEdge(mesh.nEdges());
389  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
390  {
391  const labelList& fEdges = mesh.faceEdges()[facei];
392 
393  forAll(fEdges, i)
394  {
395  isBoundaryEdge.set(fEdges[i], 1);
396  }
397  }
398 
399  const scalar featureAngle = args.argRead<scalar>(1);
400  const scalar minCos = Foam::cos(degToRad(featureAngle));
401 
402  Info<< "Feature:" << featureAngle << endl
403  << "minCos :" << minCos << endl
404  << endl;
405 
406 
407  const bool splitAllFaces = args.optionFound("splitAllFaces");
408  if (splitAllFaces)
409  {
410  Info<< "Splitting all internal faces to create multiple faces"
411  << " between two cells." << nl
412  << endl;
413  }
414 
415  const bool overwrite = args.optionFound("overwrite");
416  const bool doNotPreserveFaceZones = args.optionFound
417  (
418  "doNotPreserveFaceZones"
419  );
420  const bool concaveMultiCells = args.optionFound("concaveMultiCells");
421  if (concaveMultiCells)
422  {
423  Info<< "Generating multiple cells for points on concave feature edges."
424  << nl << endl;
425  }
426 
427 
428  // Face(centre)s that need inclusion in the dual mesh
429  labelList featureFaces;
430  // Edge(centre)s ,,
431  labelList featureEdges;
432  // Points (that become a single cell) that need inclusion in the dual mesh
433  labelList singleCellFeaturePoints;
434  // Points (that become a multiple cells) ,,
435  labelList multiCellFeaturePoints;
436 
437  // Sample implementation of feature detection.
438  simpleMarkFeatures
439  (
440  mesh,
441  isBoundaryEdge,
442  featureAngle,
443  concaveMultiCells,
444  doNotPreserveFaceZones,
445 
446  featureFaces,
447  featureEdges,
448  singleCellFeaturePoints,
449  multiCellFeaturePoints
450  );
451 
452  // If we want to split all polyMesh faces into one dualface per cell
453  // we are passing through we also need a point
454  // at the polyMesh facecentre and edgemid of the faces we want to
455  // split.
456  if (splitAllFaces)
457  {
458  featureEdges = identity(mesh.nEdges());
459  featureFaces = identity(mesh.nFaces());
460  }
461 
462  // Write obj files for debugging
463  dumpFeatures
464  (
465  mesh,
466  featureFaces,
467  featureEdges,
468  singleCellFeaturePoints,
469  multiCellFeaturePoints
470  );
471 
472 
473 
474  // Read objects in time directory
475  IOobjectList objects(mesh, runTime.timeName());
476 
477  // Read vol fields.
479  ReadFields(mesh, objects, vsFlds);
480 
482  ReadFields(mesh, objects, vvFlds);
483 
485  ReadFields(mesh, objects, vstFlds);
486 
487  PtrList<volSymmTensorField> vsymtFlds;
488  ReadFields(mesh, objects, vsymtFlds);
489 
491  ReadFields(mesh, objects, vtFlds);
492 
493  // Read surface fields.
495  ReadFields(mesh, objects, ssFlds);
496 
498  ReadFields(mesh, objects, svFlds);
499 
501  ReadFields(mesh, objects, sstFlds);
502 
504  ReadFields(mesh, objects, ssymtFlds);
505 
507  ReadFields(mesh, objects, stFlds);
508 
509 
510  // Topo change container
511  polyTopoChange meshMod(mesh.boundaryMesh().size());
512 
513  // Mesh dualiser engine
514  meshDualiser dualMaker(mesh);
515 
516  // Insert all commands into polyTopoChange to create dual of mesh. This does
517  // all the hard work.
518  dualMaker.setRefinement
519  (
520  splitAllFaces,
521  featureFaces,
522  featureEdges,
523  singleCellFeaturePoints,
524  multiCellFeaturePoints,
525  meshMod
526  );
527 
528  // Create mesh, return map from old to new mesh.
529  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
530 
531  // Update fields
532  mesh.updateMesh(map);
533 
534  // Optionally inflate mesh
535  if (map().hasMotionPoints())
536  {
537  mesh.movePoints(map().preMotionPoints());
538  }
539 
540  if (!overwrite)
541  {
542  runTime++;
543  }
544  else
545  {
546  mesh.setInstance(oldInstance);
547  }
548 
549  Info<< "Writing dual mesh to " << runTime.timeName() << endl;
550 
551  mesh.write();
552 
553  Info<< "End\n" << endl;
554 
555  return 0;
556 }
557 
558 
559 // ************************************************************************* //
Foam::surfaceFields.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: polyMesh.C:1048
wordList ReadFields(const Mesh &mesh, const IOobjectList &objects, PtrList< GeoField > &fields, const bool syncPar=true)
Read all fields of the specified type.
#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
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:50
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Unit conversion functions.
Output to file stream.
Definition: OFstream.H:81
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static void noParallel()
Remove the parallel options.
Definition: argList.C:146
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:154
label findEdge(const edgeList &edges, const labelList &candidates, const label v0, const label v1)
Return edge among candidates that uses the two vertices.
Definition: meshTools.C:337
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Field reading functions for post-processing utilities.
const vectorField & faceCentres() const
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
scalar f1
Definition: createFields.H:28
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
bool erase(const iterator &)
Erase a hashedEntry specified by given iterator.
Definition: HashTable.C:367
const word & name() const
Return name.
Definition: zone.H:150
A list of faces which address into the list of points.
A List obtained as a section of another List.
Definition: SubList.H:53
bool checkPointManifold(const bool report=false, labelHashSet *setPtr=NULL) const
Checks primitivePatch for faces sharing point but not edge.
dimensionedScalar cos(const dimensionedScalar &ds)
void set(const PackedList< 1 > &)
Set specified bits.
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
virtual void updateMesh(const mapPolyMesh &mpm)
Update the mesh corresponding to given map.
A class for handling words, derived from string.
Definition: word.H:59
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:198
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
Foam::polyBoundaryMesh.
const labelList & meshEdges() const
Return global edge index for local edges.
Definition: polyPatch.C:329
static const char nl
Definition: Ostream.H:262
Creates dual of polyMesh. Every point becomes a cell (or multiple cells for feature points)...
Definition: meshDualiser.H:66
label nEdges() const
T argRead(const label index) const
Read a value from the argument at index.
Definition: argListI.H:177
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
label nEdges() const
Return number of edges in patch.
label nFaces() const
A bit-packed bool list.
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:62
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:772
label nInternalEdges() const
Number of internal edges.
Direct mesh changes based on v1.3 polyTopoChange syntax.
virtual bool write() const
Write using setting from DB.
messageStream Info
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
Definition: argList.C:83
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
Foam::argList args(argc, argv)
const labelListList & faceEdges() const
label nInternalFaces() const
Namespace for OpenFOAM.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29