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