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