FreeFOAM The Cross-Platform CFD Toolkit
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 | 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  combinePatchFaces
26 
27 Description
28  Checks for multiple patch faces on same cell and combines them.
29 
30  These result from e.g. refined neighbouring cells getting removed, leaving
31  4 exposed faces with same owner.
32 
33  Rules for merging:
34  - only boundary faces (since multiple internal faces between two cells
35  not allowed anyway)
36  - faces have to have same owner
37  - faces have to be connected via edge which are not features (so angle
38  between them < feature angle)
39  - outside of faces has to be single loop
40  - outside of face should not be (or just slightly) concave (so angle
41  between consecutive edges < concaveangle
42 
43  E.g. to allow all faces on same patch to be merged:
44 
45  combinePatchFaces .. cavity 180 -concaveAngle 90
46 
47 Usage
48 
49  - combinePatchFaces [OPTIONS] <feature angle [0..180]>
50 
51  @param <feature angle [0..180]> \n
52  @todo Detailed description of argument.
53 
54  @param -snapMesh \n
55  Remove loose points on edges.
56 
57  @param -concaveAngle <angle [0..180]>\n
58  Maximum allowed concave angle.
59 
60  @param -overwrite \n
61  Overwrite existing data.
62 
63  @param -case <dir>\n
64  Case directory.
65 
66  @param -parallel \n
67  Run in parallel.
68 
69  @param -help \n
70  Display help message.
71 
72  @param -doc \n
73  Display Doxygen API documentation page for this application.
74 
75  @param -srcDoc \n
76  Display Doxygen source documentation page for this application.
77 
78 \*---------------------------------------------------------------------------*/
79 
81 #include <OpenFOAM/argList.H>
82 #include <OpenFOAM/Time.H>
88 #include <OpenFOAM/polyMesh.H>
89 #include <OpenFOAM/mapPolyMesh.H>
91 
92 using namespace Foam;
93 
94 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
95 
96 // Sin of angle between two consecutive edges on a face. If sin(angle) larger
97 // than this the face will be considered concave.
98 const scalar defaultConcaveAngle = 30;
99 
100 
101 // Same check as snapMesh
102 void checkSnapMesh
103 (
104  const Time& runTime,
105  const polyMesh& mesh,
106  labelHashSet& wrongFaces
107 )
108 {
109  IOdictionary snapDict
110  (
111  IOobject
112  (
113  "snapMeshDict",
114  runTime.system(),
115  mesh,
118  )
119  );
120 
121  // Max nonorthogonality allowed
122  scalar maxNonOrtho(readScalar(snapDict.lookup("maxNonOrtho")));
123  // Max concaveness allowed.
124  scalar maxConcave(readScalar(snapDict.lookup("maxConcave")));
125  // Min volume allowed (factor of minimum cellVolume)
126  scalar relMinVol(readScalar(snapDict.lookup("minVol")));
127  const scalar minCellVol = min(mesh.cellVolumes());
128  const scalar minPyrVol = relMinVol*minCellVol;
129  // Min area
130  scalar minArea(readScalar(snapDict.lookup("minArea")));
131 
132  if (maxNonOrtho < 180.0-SMALL)
133  {
134  Pout<< "Checking non orthogonality" << endl;
135 
136  label nOldSize = wrongFaces.size();
137  mesh.setNonOrthThreshold(maxNonOrtho);
138  mesh.checkFaceOrthogonality(false, &wrongFaces);
139 
140  Pout<< "Detected " << wrongFaces.size() - nOldSize
141  << " faces with non-orthogonality > " << maxNonOrtho << " degrees"
142  << endl;
143  }
144 
145  if (minPyrVol > -GREAT)
146  {
147  Pout<< "Checking face pyramids" << endl;
148 
149  label nOldSize = wrongFaces.size();
150  mesh.checkFacePyramids(false, minPyrVol, &wrongFaces);
151  Pout<< "Detected additional " << wrongFaces.size() - nOldSize
152  << " faces with illegal face pyramids" << endl;
153  }
154 
155  if (maxConcave < 180.0-SMALL)
156  {
157  Pout<< "Checking face angles" << endl;
158 
159  label nOldSize = wrongFaces.size();
160  mesh.checkFaceAngles(false, maxConcave, &wrongFaces);
161  Pout<< "Detected additional " << wrongFaces.size() - nOldSize
162  << " faces with concavity > " << maxConcave << " degrees"
163  << endl;
164  }
165 
166  if (minArea > -SMALL)
167  {
168  Pout<< "Checking face areas" << endl;
169 
170  label nOldSize = wrongFaces.size();
171 
172  const scalarField magFaceAreas = mag(mesh.faceAreas());
173 
174  forAll(magFaceAreas, faceI)
175  {
176  if (magFaceAreas[faceI] < minArea)
177  {
178  wrongFaces.insert(faceI);
179  }
180  }
181  Pout<< "Detected additional " << wrongFaces.size() - nOldSize
182  << " faces with area < " << minArea << " m^2" << endl;
183  }
184 }
185 
186 
187 // Merge faces on the same patch (usually from exposing refinement)
188 // Can undo merges if these cause problems.
189 label mergePatchFaces
190 (
191  const scalar minCos,
192  const scalar concaveSin,
193  const bool snapMeshDict,
194  const Time& runTime,
195  polyMesh& mesh
196 )
197 {
198  // Patch face merging engine
199  combineFaces faceCombiner(mesh);
200 
201  // Get all sets of faces that can be merged
202  labelListList allFaceSets(faceCombiner.getMergeSets(minCos, concaveSin));
203 
204  label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
205 
206  Info<< "Merging " << nFaceSets << " sets of faces." << endl;
207 
208  if (nFaceSets > 0)
209  {
210  // Store the faces of the face sets
211  List<faceList> allFaceSetsFaces(allFaceSets.size());
212  forAll(allFaceSets, setI)
213  {
214  allFaceSetsFaces[setI] = UIndirectList<face>
215  (
216  mesh.faces(),
217  allFaceSets[setI]
218  );
219  }
220 
222  {
223  // Topology changes container
224  polyTopoChange meshMod(mesh);
225 
226  // Merge all faces of a set into the first face of the set.
227  faceCombiner.setRefinement(allFaceSets, meshMod);
228 
229  // Change the mesh (no inflation)
230  map = meshMod.changeMesh(mesh, false, true);
231 
232  // Update fields
233  mesh.updateMesh(map);
234 
235  // Move mesh (since morphing does not do this)
236  if (map().hasMotionPoints())
237  {
238  mesh.movePoints(map().preMotionPoints());
239  }
240  else
241  {
242  // Delete mesh volumes. No other way to do this?
243  mesh.clearOut();
244  }
245  }
246 
247 
248  // Check for errors and undo
249  // ~~~~~~~~~~~~~~~~~~~~~~~~~
250 
251  // Faces in error.
252  labelHashSet errorFaces;
253 
254  if (snapMeshDict)
255  {
256  checkSnapMesh(runTime, mesh, errorFaces);
257  }
258  else
259  {
260  mesh.checkFacePyramids(false, -SMALL, &errorFaces);
261  }
262 
263  // Sets where the master is in error
264  labelHashSet errorSets;
265 
266  forAll(allFaceSets, setI)
267  {
268  label newMasterI = map().reverseFaceMap()[allFaceSets[setI][0]];
269 
270  if (errorFaces.found(newMasterI))
271  {
272  errorSets.insert(setI);
273  }
274  }
275  label nErrorSets = returnReduce(errorSets.size(), sumOp<label>());
276 
277  Info<< "Detected " << nErrorSets
278  << " error faces on boundaries that have been merged."
279  << " These will be restored to their original faces."
280  << endl;
281 
282  if (nErrorSets > 0)
283  {
284  // Renumber stored faces to new vertex numbering.
285  forAllConstIter(labelHashSet, errorSets, iter)
286  {
287  label setI = iter.key();
288 
289  faceList& setFaceVerts = allFaceSetsFaces[setI];
290 
291  forAll(setFaceVerts, i)
292  {
293  inplaceRenumber(map().reversePointMap(), setFaceVerts[i]);
294 
295  // Debug: check that all points are still there.
296  forAll(setFaceVerts[i], j)
297  {
298  label newVertI = setFaceVerts[i][j];
299 
300  if (newVertI < 0)
301  {
302  FatalErrorIn("mergePatchFaces")
303  << "In set:" << setI << " old face labels:"
304  << allFaceSets[setI] << " new face vertices:"
305  << setFaceVerts[i] << " are unmapped vertices!"
306  << abort(FatalError);
307  }
308  }
309  }
310  }
311 
312 
313  // Topology changes container
314  polyTopoChange meshMod(mesh);
315 
316 
317  // Restore faces
318  forAllConstIter(labelHashSet, errorSets, iter)
319  {
320  label setI = iter.key();
321 
322  const labelList& setFaces = allFaceSets[setI];
323  const faceList& setFaceVerts = allFaceSetsFaces[setI];
324 
325  label newMasterI = map().reverseFaceMap()[setFaces[0]];
326 
327  // Restore. Get face properties.
328 
329  label own = mesh.faceOwner()[newMasterI];
330  label zoneID = mesh.faceZones().whichZone(newMasterI);
331  bool zoneFlip = false;
332  if (zoneID >= 0)
333  {
334  const faceZone& fZone = mesh.faceZones()[zoneID];
335  zoneFlip = fZone.flipMap()[fZone.whichFace(newMasterI)];
336  }
337  label patchI = mesh.boundaryMesh().whichPatch(newMasterI);
338 
339 
340  Pout<< "Restoring new master face " << newMasterI
341  << " to vertices " << setFaceVerts[0] << endl;
342 
343  // Modify the master face.
344  meshMod.setAction
345  (
347  (
348  setFaceVerts[0], // original face
349  newMasterI, // label of face
350  own, // owner
351  -1, // neighbour
352  false, // face flip
353  patchI, // patch for face
354  false, // remove from zone
355  zoneID, // zone for face
356  zoneFlip // face flip in zone
357  )
358  );
359 
360 
361  // Add the previously removed faces
362  for (label i = 1; i < setFaces.size(); i++)
363  {
364  Pout<< "Restoring removed face " << setFaces[i]
365  << " with vertices " << setFaceVerts[i] << endl;
366 
367  meshMod.setAction
368  (
370  (
371  setFaceVerts[i], // vertices
372  own, // owner,
373  -1, // neighbour,
374  -1, // masterPointID,
375  -1, // masterEdgeID,
376  newMasterI, // masterFaceID,
377  false, // flipFaceFlux,
378  patchI, // patchID,
379  zoneID, // zoneID,
380  zoneFlip // zoneFlip
381  )
382  );
383  }
384  }
385 
386  // Change the mesh (no inflation)
387  map = meshMod.changeMesh(mesh, false, true);
388 
389  // Update fields
390  mesh.updateMesh(map);
391 
392  // Move mesh (since morphing does not do this)
393  if (map().hasMotionPoints())
394  {
395  mesh.movePoints(map().preMotionPoints());
396  }
397  else
398  {
399  // Delete mesh volumes. No other way to do this?
400  mesh.clearOut();
401  }
402  }
403  }
404  else
405  {
406  Info<< "No faces merged ..." << endl;
407  }
408 
409  return nFaceSets;
410 }
411 
412 
413 // Remove points not used by any face or points used by only two faces where
414 // the edges are in line
415 label mergeEdges(const scalar minCos, polyMesh& mesh)
416 {
417  Info<< "Merging all points on surface that" << nl
418  << "- are used by only two boundary faces and" << nl
419  << "- make an angle with a cosine of more than " << minCos
420  << "." << nl << endl;
421 
422  // Point removal analysis engine
423  removePoints pointRemover(mesh);
424 
425  // Count usage of points
426  boolList pointCanBeDeleted;
427  label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
428 
429  if (nRemove > 0)
430  {
431  Info<< "Removing " << nRemove
432  << " straight edge points ..." << endl;
433 
434  // Topology changes container
435  polyTopoChange meshMod(mesh);
436 
437  pointRemover.setRefinement(pointCanBeDeleted, meshMod);
438 
439  // Change the mesh (no inflation)
440  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
441 
442  // Update fields
443  mesh.updateMesh(map);
444 
445  // Move mesh (since morphing does not do this)
446  if (map().hasMotionPoints())
447  {
448  mesh.movePoints(map().preMotionPoints());
449  }
450  else
451  {
452  // Delete mesh volumes. No other way to do this?
453  mesh.clearOut();
454  }
455  }
456  else
457  {
458  Info<< "No straight edges simplified and no points removed ..." << endl;
459  }
460 
461  return nRemove;
462 }
463 
464 
465 // Main program:
466 
467 int main(int argc, char *argv[])
468 {
469  argList::validArgs.append("feature angle [0..180]");
470  argList::validOptions.insert("concaveAngle", "[0..180]");
471  argList::validOptions.insert("snapMesh", "");
472  argList::validOptions.insert("overwrite", "");
473 
474 # include <OpenFOAM/setRootCase.H>
475 # include <OpenFOAM/createTime.H>
476  runTime.functionObjects().off();
477 # include <OpenFOAM/createPolyMesh.H>
478  const word oldInstance = mesh.pointsInstance();
479 
480  scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
481 
482  scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.0);
483 
484  scalar concaveAngle = defaultConcaveAngle;
485  args.optionReadIfPresent("concaveAngle", concaveAngle);
486 
487  scalar concaveSin = Foam::sin(concaveAngle*mathematicalConstant::pi/180.0);
488 
489  bool snapMeshDict = args.optionFound("snapMesh");
490  bool overwrite = args.optionFound("overwrite");
491 
492  Info<< "Merging all faces of a cell" << nl
493  << " - which are on the same patch" << nl
494  << " - which make an angle < " << featureAngle << " degrees"
495  << nl
496  << " (cos:" << minCos << ')' << nl
497  << " - even when resulting face becomes concave by more than "
498  << concaveAngle << " degrees" << nl
499  << " (sin:" << concaveSin << ')' << nl
500  << endl;
501 
502  if (!overwrite)
503  {
504  runTime++;
505  }
506 
507  // Merge faces on same patch
508  label nChanged = mergePatchFaces
509  (
510  minCos,
511  concaveSin,
512  snapMeshDict,
513  runTime,
514  mesh
515  );
516 
517  // Merge points on straight edges and remove unused points
518  if (snapMeshDict)
519  {
520  Info<< "Merging all 'loose' points on surface edges"
521  << ", regardless of the angle they make." << endl;
522 
523  // Surface bnound to be used to extrude. Merge all loose points.
524  nChanged += mergeEdges(-1, mesh);
525  }
526  else
527  {
528  nChanged += mergeEdges(minCos, mesh);
529  }
530 
531  if (nChanged > 0)
532  {
533  if (overwrite)
534  {
535  mesh.setInstance(oldInstance);
536  }
537 
538  Info<< "Writing morphed mesh to time " << runTime.timeName() << endl;
539 
540  mesh.write();
541  }
542  else
543  {
544  Info<< "Mesh unchanged." << endl;
545  }
546 
547  Info << "End\n" << endl;
548 
549  return 0;
550 }
551 
552 
553 // ************************ vim: set sw=4 sts=4 et: ************************ //