combinePatchFaces.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-2015 OpenFOAM Foundation
9 Copyright (C) 2016-2020 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 combinePatchFaces
29
30Group
31 grpMeshAdvancedUtilities
32
33Description
34 Checks for multiple patch faces on the same cell and combines them.
35 Multiple patch faces can result from e.g. removal of refined
36 neighbouring cells, leaving 4 exposed faces with same owner.
37
38 Rules for merging:
39 - only boundary faces (since multiple internal faces between two cells
40 not allowed anyway)
41 - faces have to have same owner
42 - faces have to be connected via edge which are not features (so angle
43 between them < feature angle)
44 - outside of faces has to be single loop
45 - outside of face should not be (or just slightly) concave (so angle
46 between consecutive edges < concaveangle
47
48 E.g. to allow all faces on same patch to be merged:
49
50 combinePatchFaces 180 -concaveAngle 90
51
52\*---------------------------------------------------------------------------*/
53
54#include "argList.H"
55#include "Time.H"
56#include "polyTopoChange.H"
57#include "polyModifyFace.H"
58#include "polyAddFace.H"
59#include "combineFaces.H"
60#include "removePoints.H"
61#include "polyMesh.H"
62#include "mapPolyMesh.H"
63#include "unitConversion.H"
64#include "motionSmoother.H"
65#include "topoSet.H"
66#include "processorMeshes.H"
67#include "PstreamReduceOps.H"
68
69using namespace Foam;
70
71// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
72
73// Merge faces on the same patch (usually from exposing refinement)
74// Can undo merges if these cause problems.
75label mergePatchFaces
76(
77 const scalar minCos,
78 const scalar concaveSin,
79 const autoPtr<IOdictionary>& qualDictPtr,
80 const Time& runTime,
82)
83{
84 // Patch face merging engine
85 combineFaces faceCombiner(mesh);
86
87 // Get all sets of faces that can be merged
88 labelListList allFaceSets(faceCombiner.getMergeSets(minCos, concaveSin));
89
90 label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
91
92 Info<< "Merging " << nFaceSets << " sets of faces." << endl;
93
94 if (nFaceSets > 0)
95 {
96 // Store the faces of the face sets
97 List<faceList> allFaceSetsFaces(allFaceSets.size());
98 forAll(allFaceSets, seti)
99 {
100 allFaceSetsFaces[seti] = UIndirectList<face>
101 (
102 mesh.faces(),
103 allFaceSets[seti]
104 );
105 }
106
108 {
109 // Topology changes container
110 polyTopoChange meshMod(mesh);
111
112 // Merge all faces of a set into the first face of the set.
113 faceCombiner.setRefinement(allFaceSets, meshMod);
114
115 // Change the mesh (no inflation)
116 map = meshMod.changeMesh(mesh, false, true);
117
118 // Update fields
119 mesh.updateMesh(map());
120
121 // Move mesh (since morphing does not do this)
122 if (map().hasMotionPoints())
123 {
124 mesh.movePoints(map().preMotionPoints());
125 }
126 else
127 {
128 // Delete mesh volumes. No other way to do this?
129 mesh.clearOut();
130 }
131 }
132
133
134 // Check for errors and undo
135 // ~~~~~~~~~~~~~~~~~~~~~~~~~
136
137 // Faces in error.
138 labelHashSet errorFaces;
139
140 if (qualDictPtr)
141 {
142 motionSmoother::checkMesh(false, mesh, *qualDictPtr, errorFaces);
143 }
144 else
145 {
146 mesh.checkFacePyramids(false, -SMALL, &errorFaces);
147 }
148
149 // Sets where the master is in error
150 labelHashSet errorSets;
151
152 forAll(allFaceSets, seti)
153 {
154 label newMasterI = map().reverseFaceMap()[allFaceSets[seti][0]];
155
156 if (errorFaces.found(newMasterI))
157 {
158 errorSets.insert(seti);
159 }
160 }
161 label nErrorSets = returnReduce(errorSets.size(), sumOp<label>());
162
163 Info<< "Detected " << nErrorSets
164 << " error faces on boundaries that have been merged."
165 << " These will be restored to their original faces."
166 << endl;
167
168 if (nErrorSets)
169 {
170 // Renumber stored faces to new vertex numbering.
171 for (const label seti : errorSets)
172 {
173 faceList& setFaceVerts = allFaceSetsFaces[seti];
174
175 forAll(setFaceVerts, i)
176 {
177 inplaceRenumber(map().reversePointMap(), setFaceVerts[i]);
178
179 // Debug: check that all points are still there.
180 forAll(setFaceVerts[i], j)
181 {
182 label newVertI = setFaceVerts[i][j];
183
184 if (newVertI < 0)
185 {
187 << "In set:" << seti << " old face labels:"
188 << allFaceSets[seti] << " new face vertices:"
189 << setFaceVerts[i] << " are unmapped vertices!"
190 << abort(FatalError);
191 }
192 }
193 }
194 }
195
196
197 // Topology changes container
198 polyTopoChange meshMod(mesh);
199
200
201 // Restore faces
202 for (const label seti : errorSets)
203 {
204 const labelList& setFaces = allFaceSets[seti];
205 const faceList& setFaceVerts = allFaceSetsFaces[seti];
206
207 label newMasterI = map().reverseFaceMap()[setFaces[0]];
208
209 // Restore. Get face properties.
210
211 label own = mesh.faceOwner()[newMasterI];
212 label zoneID = mesh.faceZones().whichZone(newMasterI);
213 bool zoneFlip = false;
214 if (zoneID >= 0)
215 {
216 const faceZone& fZone = mesh.faceZones()[zoneID];
217 zoneFlip = fZone.flipMap()[fZone.whichFace(newMasterI)];
218 }
219 label patchID = mesh.boundaryMesh().whichPatch(newMasterI);
220
221 Pout<< "Restoring new master face " << newMasterI
222 << " to vertices " << setFaceVerts[0] << endl;
223
224 // Modify the master face.
225 meshMod.setAction
226 (
228 (
229 setFaceVerts[0], // original face
230 newMasterI, // label of face
231 own, // owner
232 -1, // neighbour
233 false, // face flip
234 patchID, // patch for face
235 false, // remove from zone
236 zoneID, // zone for face
237 zoneFlip // face flip in zone
238 )
239 );
240
241
242 // Add the previously removed faces
243 for (label i = 1; i < setFaces.size(); ++i)
244 {
245 Pout<< "Restoring removed face " << setFaces[i]
246 << " with vertices " << setFaceVerts[i] << endl;
247
248 meshMod.setAction
249 (
251 (
252 setFaceVerts[i], // vertices
253 own, // owner,
254 -1, // neighbour,
255 -1, // masterPointID,
256 -1, // masterEdgeID,
257 newMasterI, // masterFaceID,
258 false, // flipFaceFlux,
259 patchID, // patchID,
260 zoneID, // zoneID,
261 zoneFlip // zoneFlip
262 )
263 );
264 }
265 }
266
267 // Change the mesh (no inflation)
268 map = meshMod.changeMesh(mesh, false, true);
269
270 // Update fields
271 mesh.updateMesh(map());
272
273 // Move mesh (since morphing does not do this)
274 if (map().hasMotionPoints())
275 {
276 mesh.movePoints(map().preMotionPoints());
277 }
278 else
279 {
280 // Delete mesh volumes. No other way to do this?
281 mesh.clearOut();
282 }
283 }
284 }
285 else
286 {
287 Info<< "No faces merged ..." << endl;
288 }
289
290 return nFaceSets;
291}
292
293
294// Remove points not used by any face or points used by only two faces where
295// the edges are in line
296label mergeEdges(const scalar minCos, polyMesh& mesh)
297{
298 Info<< "Merging all points on surface that" << nl
299 << "- are used by only two boundary faces and" << nl
300 << "- make an angle with a cosine of more than " << minCos
301 << "." << nl << endl;
302
303 // Point removal analysis engine
304 removePoints pointRemover(mesh);
305
306 // Count usage of points
307 boolList pointCanBeDeleted;
308 label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
309
310 if (nRemove > 0)
311 {
312 Info<< "Removing " << nRemove
313 << " straight edge points ..." << endl;
314
315 // Topology changes container
316 polyTopoChange meshMod(mesh);
317
318 pointRemover.setRefinement(pointCanBeDeleted, meshMod);
319
320 // Change the mesh (no inflation)
321 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
322
323 // Update fields
324 mesh.updateMesh(map());
325
326 // Move mesh (since morphing does not do this)
327 if (map().hasMotionPoints())
328 {
329 mesh.movePoints(map().preMotionPoints());
330 }
331 else
332 {
333 // Delete mesh volumes. No other way to do this?
334 mesh.clearOut();
335 }
336 }
337 else
338 {
339 Info<< "No straight edges simplified and no points removed ..." << endl;
340 }
341
342 return nRemove;
343}
344
345
346
347int main(int argc, char *argv[])
348{
349 argList::addNote
350 (
351 "Checks for multiple patch faces on the same cell and combines them."
352 );
353
354 #include "addOverwriteOption.H"
355
356 argList::addArgument
357 (
358 "featureAngle",
359 "in degrees [0-180]"
360 );
361 argList::addOption
362 (
363 "concaveAngle",
364 "degrees",
365 "Specify concave angle [0..180] (default: 30 degrees)"
366 );
367 argList::addBoolOption
368 (
369 "meshQuality",
370 "Read user-defined mesh quality criteria from system/meshQualityDict"
371 );
372
373 argList::noFunctionObjects(); // Never use function objects
374
375 #include "setRootCase.H"
376 #include "createTime.H"
377 #include "createPolyMesh.H"
378
379 const word oldInstance = mesh.pointsInstance();
380
381 const scalar featureAngle = args.get<scalar>(1);
382 const scalar minCos = Foam::cos(degToRad(featureAngle));
383
384 // Sin of angle between two consecutive edges on a face.
385 // If sin(angle) larger than this the face will be considered concave.
386 const scalar concaveAngle = args.getOrDefault<scalar>("concaveAngle", 30);
387
388 const scalar concaveSin = Foam::sin(degToRad(concaveAngle));
389
390 const bool overwrite = args.found("overwrite");
391 const bool meshQuality = args.found("meshQuality");
392
393 Info<< "Merging all faces of a cell" << nl
394 << " - which are on the same patch" << nl
395 << " - which make an angle < " << featureAngle << " degrees"
396 << nl
397 << " (cos:" << minCos << ')' << nl
398 << " - even when resulting face becomes concave by more than "
399 << concaveAngle << " degrees" << nl
400 << " (sin:" << concaveSin << ')' << nl
401 << endl;
402
403 autoPtr<IOdictionary> qualDict;
404 if (meshQuality)
405 {
406 Info<< "Enabling user-defined geometry checks." << nl << endl;
407
408 qualDict.reset
409 (
410 new IOdictionary
411 (
413 (
414 "meshQualityDict",
415 mesh.time().system(),
416 mesh,
417 IOobject::MUST_READ,
418 IOobject::NO_WRITE
419 )
420 )
421 );
422 }
423
424
425 if (!overwrite)
426 {
427 ++runTime;
428 }
429
430
431
432 // Merge faces on same patch
433 label nChanged = mergePatchFaces
434 (
435 minCos,
436 concaveSin,
437 qualDict,
438 runTime,
439 mesh
440 );
441
442 // Merge points on straight edges and remove unused points
443 if (qualDict)
444 {
445 Info<< "Merging all 'loose' points on surface edges, "
446 << "regardless of the angle they make." << endl;
447
448 // Surface bnound to be used to extrude. Merge all loose points.
449 nChanged += mergeEdges(-1, mesh);
450 }
451 else
452 {
453 nChanged += mergeEdges(minCos, mesh);
454 }
455
456 if (nChanged > 0)
457 {
458 if (overwrite)
459 {
460 mesh.setInstance(oldInstance);
461 }
462
463 Info<< "Writing morphed mesh to time " << runTime.timeName() << endl;
464
465 mesh.write();
466 topoSet::removeFiles(mesh);
467 processorMeshes::removeFiles(mesh);
468 }
469 else
470 {
471 Info<< "Mesh unchanged." << endl;
472 }
473
474 Info<< "\nEnd\n" << endl;
475
476 return 0;
477}
478
479
480// ************************************************************************* //
Inter-processor communication reduction functions.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
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
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
Combines boundary faces into single face. The faces get the patch of the first face ('the master')
Definition: combineFaces.H:59
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:372
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
A face addition data class. A face can be inflated either from a point or from another face and can e...
Definition: polyAddFace.H:56
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Class describing modification of a face.
Direct mesh changes based on v1.3 polyTopoChange syntax.
Removes selected points from mesh and updates faces using these points.
Definition: removePoints.H:61
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelIOList & zoneID
void inplaceRenumber(const labelUList &oldToNew, IntListType &lists)
Inplace renumber the values (not the indices) of a list of lists.
Definition: ensightOutput.H:90
Namespace for OpenFOAM.
dimensionedScalar sin(const dimensionedScalar &ds)
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
errorManip< error > abort(error &err)
Definition: errorManip.H:144
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.