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