applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.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     Calculate the dual of a polyMesh. Adheres to all the feature and patch
  26     edges.
  27
  28 Usage
  29
  30     - polyDualMesh featureAngle
  31
  32     Detects any boundary edge > angle and creates multiple boundary faces
  33     for it. Normal behaviour is to have each point become a cell
  34     (1.5 behaviour)
  35
  36     \param -concaveMultiCells
  37     Creates multiple cells for each point on a concave edge. Might limit
  38     the amount of distortion on some meshes.
  39
  40     \param -splitAllFaces
  41     Normally only constructs a single face between two cells. This single face
  42     might be too distorted. splitAllFaces will create a single face for every
  43     original cell the face passes through. The mesh will thus have
  44     multiple faces inbetween two cells! (so is not strictly upper-triangular
  45     anymore - checkMesh will complain)
  46
  47     \param -doNotPreserveFaceZones:
  48     By default all faceZones are preserved by marking all faces, edges and
  49     points on them as features. The -doNotPreserveFaceZones disables this
  50     behaviour.
  51
  52     Note: is just a driver for meshDualiser. Substitute your own
  53     simpleMarkFeatures to have different behaviour.
  54
  55 \*---------------------------------------------------------------------------*/
  56
  57 #include "argList.H"
  58 #include "Time.H"
  59 #include "timeSelector.H"
  60 #include "fvMesh.H"
  61 #include "unitConversion.H"
  62 #include "polyTopoChange.H"
  63 #include "mapPolyMesh.H"
  64 #include "PackedBoolList.H"
  65 #include "meshTools.H"
  66 #include "OFstream.H"
  67 #include "meshDualiser.H"
  68 #include "ReadFields.H"
  69 #include "volFields.H"
  70 #include "surfaceFields.H"
  71
  72 using namespace Foam;
  73
  74 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  75
  76 // Naive feature detection. All boundary edges with angle > featureAngle become
  77 // feature edges. All points on feature edges become feature points. All
  78 // boundary faces become feature faces.
  79 void simpleMarkFeatures
  80 (
  81     const polyMesh& mesh,
  82     const PackedBoolList& isBoundaryEdge,
  83     const scalar featureAngle,
  84     const bool concaveMultiCells,
  85     const bool doNotPreserveFaceZones,
  86
  87     labelList& featureFaces,
  88     labelList& featureEdges,
  89     labelList& singleCellFeaturePoints,
  90     labelList& multiCellFeaturePoints
  91 )
  92 {
  93     scalar minCos = Foam::cos(degToRad(featureAngle));
  94
  95     const polyBoundaryMesh& patches = mesh.boundaryMesh();
  96
  97     // Working sets
  98     labelHashSet featureEdgeSet;
  99     labelHashSet singleCellFeaturePointSet;
 100     labelHashSet multiCellFeaturePointSet;
 101
 102
 103     // 1. Mark all edges between patches
 104     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 105
 106     forAll(patches, patchI)
 107     {
 108         const polyPatch& pp = patches[patchI];
 109         const labelList& meshEdges = pp.meshEdges();
 110
 111         // All patch corner edges. These need to be feature points & edges!
 112         for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
 113         {
 114             label meshEdgeI = meshEdges[edgeI];
 115             featureEdgeSet.insert(meshEdgeI);
 116             singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][0]);
 117             singleCellFeaturePointSet.insert(mesh.edges()[meshEdgeI][1]);
 118         }
 119     }
 120
 121
 122
 123     // 2. Mark all geometric feature edges
 124     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 125     // Make distinction between convex features where the boundary point becomes
 126     // a single cell and concave features where the boundary point becomes
 127     // multiple 'half' cells.
 128
 129     // Addressing for all outside faces
 130     primitivePatch allBoundary
 131     (
 132         SubList<face>
 133         (
 134             mesh.faces(),
 135             mesh.nFaces()-mesh.nInternalFaces(),
 136             mesh.nInternalFaces()
 137         ),
 138         mesh.points()
 139     );
 140
 141     // Check for non-manifold points (surface pinched at point)
 142     allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);
 143
 144     // Check for non-manifold edges (surface pinched at edge)
 145     const labelListList& edgeFaces = allBoundary.edgeFaces();
 146     const labelList& meshPoints = allBoundary.meshPoints();
 147
 148     forAll(edgeFaces, edgeI)
 149     {
 150         const labelList& eFaces = edgeFaces[edgeI];
 151
 152         if (eFaces.size() > 2)
 153         {
 154             const edge& e = allBoundary.edges()[edgeI];
 155
 156             //Info<< "Detected non-manifold boundary edge:" << edgeI
 157             //    << " coords:"
 158             //    << allBoundary.points()[meshPoints[e[0]]]
 159             //    << allBoundary.points()[meshPoints[e[1]]] << endl;
 160
 161             singleCellFeaturePointSet.insert(meshPoints[e[0]]);
 162             singleCellFeaturePointSet.insert(meshPoints[e[1]]);
 163         }
 164     }
 165
 166     // Check for features.
 167     forAll(edgeFaces, edgeI)
 168     {
 169         const labelList& eFaces = edgeFaces[edgeI];
 170
 171         if (eFaces.size() == 2)
 172         {
 173             label f0 = eFaces[0];
 174             label f1 = eFaces[1];
 175
 176             // check angle
 177             const vector& n0 = allBoundary.faceNormals()[f0];
 178             const vector& n1 = allBoundary.faceNormals()[f1];
 179
 180             if ((n0 & n1) < minCos)
 181             {
 182                 const edge& e = allBoundary.edges()[edgeI];
 183                 label v0 = meshPoints[e[0]];
 184                 label v1 = meshPoints[e[1]];
 185
 186                 label meshEdgeI = meshTools::findEdge(mesh, v0, v1);
 187                 featureEdgeSet.insert(meshEdgeI);
 188
 189                 // Check if convex or concave by looking at angle
 190                 // between face centres and normal
 191                 vector c1c0
 192                 (
 193                     allBoundary[f1].centre(allBoundary.points())
 194                   - allBoundary[f0].centre(allBoundary.points())
 195                 );
 196
 197                 if (concaveMultiCells && (c1c0 & n0) > SMALL)
 198                 {
 199                     // Found concave edge. Make into multiCell features
 200                     Info<< "Detected concave feature edge:" << edgeI
 201                         << " cos:" << (c1c0 & n0)
 202                         << " coords:"
 203                         << allBoundary.points()[v0]
 204                         << allBoundary.points()[v1]
 205                         << endl;
 206
 207                     singleCellFeaturePointSet.erase(v0);
 208                     multiCellFeaturePointSet.insert(v0);
 209                     singleCellFeaturePointSet.erase(v1);
 210                     multiCellFeaturePointSet.insert(v1);
 211                 }
 212                 else
 213                 {
 214                     // Convex. singleCell feature.
 215                     if (!multiCellFeaturePointSet.found(v0))
 216                     {
 217                         singleCellFeaturePointSet.insert(v0);
 218                     }
 219                     if (!multiCellFeaturePointSet.found(v1))
 220                     {
 221                         singleCellFeaturePointSet.insert(v1);
 222                     }
 223                 }
 224             }
 225         }
 226     }
 227
 228
 229     // 3. Mark all feature faces
 230     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 231
 232     // Face centres that need inclusion in the dual mesh
 233     labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces());
 234     // A. boundary faces.
 235     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 236     {
 237         featureFaceSet.insert(faceI);
 238     }
 239
 240     // B. face zones.
 241     const faceZoneMesh& faceZones = mesh.faceZones();
 242
 243     if (doNotPreserveFaceZones)
 244     {
 245         if (faceZones.size() > 0)
 246         {
 247             WarningIn("simpleMarkFeatures(..)")
 248                 << "Detected " << faceZones.size()
 249                 << " faceZones. These will not be preserved."
 250                 << endl;
 251         }
 252     }
 253     else
 254     {
 255         if (faceZones.size() > 0)
 256         {
 257             Info<< "Detected " << faceZones.size()
 258                 << " faceZones. Preserving these by marking their"
 259                 << " points, edges and faces as features." << endl;
 260         }
 261
 262         forAll(faceZones, zoneI)
 263         {
 264             const faceZone& fz = faceZones[zoneI];
 265
 266             Info<< "Inserting all faces in faceZone " << fz.name()
 267                 << " as features." << endl;
 268
 269             forAll(fz, i)
 270             {
 271                 label faceI = fz[i];
 272                 const face& f = mesh.faces()[faceI];
 273                 const labelList& fEdges = mesh.faceEdges()[faceI];
 274
 275                 featureFaceSet.insert(faceI);
 276                 forAll(f, fp)
 277                 {
 278                     // Mark point as multi cell point (since both sides of
 279                     // face should have different cells)
 280                     singleCellFeaturePointSet.erase(f[fp]);
 281                     multiCellFeaturePointSet.insert(f[fp]);
 282
 283                     // Make sure there are points on the edges.
 284                     featureEdgeSet.insert(fEdges[fp]);
 285                 }
 286             }
 287         }
 288     }
 289
 290     // Transfer to arguments
 291     featureFaces = featureFaceSet.toc();
 292     featureEdges = featureEdgeSet.toc();
 293     singleCellFeaturePoints = singleCellFeaturePointSet.toc();
 294     multiCellFeaturePoints = multiCellFeaturePointSet.toc();
 295 }
 296
 297
 298 // Dump features to .obj files
 299 void dumpFeatures
 300 (
 301     const polyMesh& mesh,
 302     const labelList& featureFaces,
 303     const labelList& featureEdges,
 304     const labelList& singleCellFeaturePoints,
 305     const labelList& multiCellFeaturePoints
 306 )
 307 {
 308     {
 309         OFstream str("featureFaces.obj");
 310         Info<< "Dumping centres of featureFaces to obj file " << str.name()
 311             << endl;
 312         forAll(featureFaces, i)
 313         {
 314             meshTools::writeOBJ(str, mesh.faceCentres()[featureFaces[i]]);
 315         }
 316     }
 317     {
 318         OFstream str("featureEdges.obj");
 319         Info<< "Dumping featureEdges to obj file " << str.name() << endl;
 320         label vertI = 0;
 321
 322         forAll(featureEdges, i)
 323         {
 324             const edge& e = mesh.edges()[featureEdges[i]];
 325             meshTools::writeOBJ(str, mesh.points()[e[0]]);
 326             vertI++;
 327             meshTools::writeOBJ(str, mesh.points()[e[1]]);
 328             vertI++;
 329             str<< "l " << vertI-1 << ' ' << vertI << nl;
 330         }
 331     }
 332     {
 333         OFstream str("singleCellFeaturePoints.obj");
 334         Info<< "Dumping featurePoints that become a single cell to obj file "
 335             << str.name() << endl;
 336         forAll(singleCellFeaturePoints, i)
 337         {
 338             meshTools::writeOBJ(str, mesh.points()[singleCellFeaturePoints[i]]);
 339         }
 340     }
 341     {
 342         OFstream str("multiCellFeaturePoints.obj");
 343         Info<< "Dumping featurePoints that become multiple cells to obj file "
 344             << str.name() << endl;
 345         forAll(multiCellFeaturePoints, i)
 346         {
 347             meshTools::writeOBJ(str, mesh.points()[multiCellFeaturePoints[i]]);
 348         }
 349     }
 350 }
 351
 352
 353 int main(int argc, char *argv[])
 354 {
 355 #   include "addOverwriteOption.H"
 356     argList::noParallel();
 357     timeSelector::addOptions(true, false);
 358
 359     argList::validArgs.append("featureAngle [0-180]");
 360     argList::addBoolOption
 361     (
 362         "splitAllFaces",
 363         "have multiple faces inbetween cells"
 364     );
 365     argList::addBoolOption
 366     (
 367         "concaveMultiCells",
 368         "split cells on concave boundary edges into multiple cells"
 369     );
 370     argList::addBoolOption
 371     (
 372         "doNotPreserveFaceZones",
 373         "disable the default behaviour of preserving faceZones by having"
 374         " multiple faces inbetween cells"
 375     );
 376
 377 #   include "setRootCase.H"
 378 #   include "createTime.H"
 379
 380     instantList timeDirs = timeSelector::select0(runTime, args);
 381
 382 #   include "createMesh.H"
 383
 384     const word oldInstance = mesh.pointsInstance();
 385
 386     // Mark boundary edges and points.
 387     // (Note: in 1.4.2 we can use the built-in mesh point ordering
 388     //  facility instead)
 389     PackedBoolList isBoundaryEdge(mesh.nEdges());
 390     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 391     {
 392         const labelList& fEdges = mesh.faceEdges()[faceI];
 393
 394         forAll(fEdges, i)
 395         {
 396             isBoundaryEdge.set(fEdges[i], 1);
 397         }
 398     }
 399
 400     const scalar featureAngle = args.argRead<scalar>(1);
 401     const scalar minCos = Foam::cos(degToRad(featureAngle));
 402
 403     Info<< "Feature:" << featureAngle << endl
 404         << "minCos :" << minCos << endl
 405         << endl;
 406
 407
 408     const bool splitAllFaces = args.optionFound("splitAllFaces");
 409     if (splitAllFaces)
 410     {
 411         Info<< "Splitting all internal faces to create multiple faces"
 412             << " between two cells." << nl
 413             << endl;
 414     }
 415
 416     const bool overwrite = args.optionFound("overwrite");
 417     const bool doNotPreserveFaceZones = args.optionFound
 418     (
 419         "doNotPreserveFaceZones"
 420     );
 421     const bool concaveMultiCells = args.optionFound("concaveMultiCells");
 422     if (concaveMultiCells)
 423     {
 424         Info<< "Generating multiple cells for points on concave feature edges."
 425             << nl << endl;
 426     }
 427
 428
 429     // Face(centre)s that need inclusion in the dual mesh
 430     labelList featureFaces;
 431     // Edge(centre)s  ,,
 432     labelList featureEdges;
 433     // Points (that become a single cell) that need inclusion in the dual mesh
 434     labelList singleCellFeaturePoints;
 435     // Points (that become a multiple cells)        ,,
 436     labelList multiCellFeaturePoints;
 437
 438     // Sample implementation of feature detection.
 439     simpleMarkFeatures
 440     (
 441         mesh,
 442         isBoundaryEdge,
 443         featureAngle,
 444         concaveMultiCells,
 445         doNotPreserveFaceZones,
 446
 447         featureFaces,
 448         featureEdges,
 449         singleCellFeaturePoints,
 450         multiCellFeaturePoints
 451     );
 452
 453     // If we want to split all polyMesh faces into one dualface per cell
 454     // we are passing through we also need a point
 455     // at the polyMesh facecentre and edgemid of the faces we want to
 456     // split.
 457     if (splitAllFaces)
 458     {
 459         featureEdges = identity(mesh.nEdges());
 460         featureFaces = identity(mesh.nFaces());
 461     }
 462
 463     // Write obj files for debugging
 464     dumpFeatures
 465     (
 466         mesh,
 467         featureFaces,
 468         featureEdges,
 469         singleCellFeaturePoints,
 470         multiCellFeaturePoints
 471     );
 472
 473
 474
 475     // Read objects in time directory
 476     IOobjectList objects(mesh, runTime.timeName());
 477
 478     // Read vol fields.
 479     PtrList<volScalarField> vsFlds;
 480     ReadFields(mesh, objects, vsFlds);
 481
 482     PtrList<volVectorField> vvFlds;
 483     ReadFields(mesh, objects, vvFlds);
 484
 485     PtrList<volSphericalTensorField> vstFlds;
 486     ReadFields(mesh, objects, vstFlds);
 487
 488     PtrList<volSymmTensorField> vsymtFlds;
 489     ReadFields(mesh, objects, vsymtFlds);
 490
 491     PtrList<volTensorField> vtFlds;
 492     ReadFields(mesh, objects, vtFlds);
 493
 494     // Read surface fields.
 495     PtrList<surfaceScalarField> ssFlds;
 496     ReadFields(mesh, objects, ssFlds);
 497
 498     PtrList<surfaceVectorField> svFlds;
 499     ReadFields(mesh, objects, svFlds);
 500
 501     PtrList<surfaceSphericalTensorField> sstFlds;
 502     ReadFields(mesh, objects, sstFlds);
 503
 504     PtrList<surfaceSymmTensorField> ssymtFlds;
 505     ReadFields(mesh, objects, ssymtFlds);
 506
 507     PtrList<surfaceTensorField> stFlds;
 508     ReadFields(mesh, objects, stFlds);
 509
 510
 511     // Topo change container
 512     polyTopoChange meshMod(mesh.boundaryMesh().size());
 513
 514     // Mesh dualiser engine
 515     meshDualiser dualMaker(mesh);
 516
 517     // Insert all commands into polyTopoChange to create dual of mesh. This does
 518     // all the hard work.
 519     dualMaker.setRefinement
 520     (
 521         splitAllFaces,
 522         featureFaces,
 523         featureEdges,
 524         singleCellFeaturePoints,
 525         multiCellFeaturePoints,
 526         meshMod
 527     );
 528
 529     // Create mesh, return map from old to new mesh.
 530     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
 531
 532     // Update fields
 533     mesh.updateMesh(map);
 534
 535     // Optionally inflate mesh
 536     if (map().hasMotionPoints())
 537     {
 538         mesh.movePoints(map().preMotionPoints());
 539     }
 540
 541     if (!overwrite)
 542     {
 543         runTime++;
 544     }
 545     else
 546     {
 547         mesh.setInstance(oldInstance);
 548     }
 549
 550     Info<< "Writing dual mesh to " << runTime.timeName() << endl;
 551
 552     mesh.write();
 553
 554     Info<< "End\n" << endl;
 555
 556     return 0;
 557 }
 558
 559
 560 // ************************************************************************* //