applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C

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