applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C

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