applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.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 Class
  25     meshDualiser
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "meshDualiser.H"
  30 #include "meshTools.H"
  31 #include "polyMesh.H"
  32 #include "directTopoChange.H"
  33 #include "mapPolyMesh.H"
  34 #include "edgeFaceCirculator.H"
  35 #include "mergePoints.H"
  36 #include "OFstream.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::meshDualiser, 0);
  41
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 void Foam::meshDualiser::checkPolyTopoChange(const directTopoChange& meshMod)
  46 {
  47     // Assume no removed points
  48     pointField points(meshMod.points().size());
  49     forAll(meshMod.points(), i)
  50     {
  51         points[i] = meshMod.points()[i];
  52     }
  53
  54     labelList oldToNew;
  55     pointField newPoints;
  56     bool hasMerged = mergePoints
  57     (
  58         points,
  59         1E-6,
  60         false,
  61         oldToNew,
  62         newPoints
  63     );
  64
  65     if (hasMerged)
  66     {
  67         labelListList newToOld(invertOneToMany(newPoints.size(), oldToNew));
  68
  69         forAll(newToOld, newI)
  70         {
  71             if (newToOld[newI].size() != 1)
  72             {
  73                 FatalErrorIn
  74                 (
  75                     "meshDualiser::checkPolyTopoChange(const directTopoChange&)"
  76                 )   << "duplicate verts:" << newToOld[newI]
  77                     << " coords:"
  78                     << UIndirectList<point>(points, newToOld[newI])()
  79                     << abort(FatalError);
  80             }
  81         }
  82     }
  83 }
  84
  85
  86 // Dump state so far.
  87 void Foam::meshDualiser::dumpPolyTopoChange
  88 (
  89     const directTopoChange& meshMod,
  90     const fileName& prefix
  91 )
  92 {
  93     OFstream str1(prefix + "Faces.obj");
  94     OFstream str2(prefix + "Edges.obj");
  95
  96     Info<< "Dumping current directTopoChange. Faces to " << str1.name()
  97         << " , points and edges to " << str2.name() << endl;
  98
  99     const DynamicList<point>& points = meshMod.points();
 100
 101     forAll(points, pointI)
 102     {
 103         meshTools::writeOBJ(str1, points[pointI]);
 104         meshTools::writeOBJ(str2, points[pointI]);
 105     }
 106
 107     const DynamicList<face>& faces = meshMod.faces();
 108
 109     forAll(faces, faceI)
 110     {
 111         const face& f = faces[faceI];
 112
 113         str1<< 'f';
 114         forAll(f, fp)
 115         {
 116             str1<< ' ' << f[fp]+1;
 117         }
 118         str1<< nl;
 119
 120         str2<< 'l';
 121         forAll(f, fp)
 122         {
 123             str2<< ' ' << f[fp]+1;
 124         }
 125         str2<< ' ' << f[0]+1 << nl;
 126     }
 127 }
 128
 129
 130 //- Given cell and point on mesh finds the corresponding dualCell. Most
 131 //  points become only one cell but the feature points might be split.
 132 Foam::label Foam::meshDualiser::findDualCell
 133 (
 134     const label cellI,
 135     const label pointI
 136 ) const
 137 {
 138     const labelList& dualCells = pointToDualCells_[pointI];
 139
 140     if (dualCells.size() == 1)
 141     {
 142         return dualCells[0];
 143     }
 144     else
 145     {
 146         label index = findIndex(mesh_.pointCells()[pointI], cellI);
 147
 148         return dualCells[index];
 149     }
 150 }
 151
 152
 153 // Helper function to generate dualpoints on all boundary edges emanating
 154 // from (boundary & feature) point
 155 void Foam::meshDualiser::generateDualBoundaryEdges
 156 (
 157     const PackedBoolList& isBoundaryEdge,
 158     const label pointI,
 159     directTopoChange& meshMod
 160 )
 161 {
 162     const labelList& pEdges = mesh_.pointEdges()[pointI];
 163
 164     forAll(pEdges, pEdgeI)
 165     {
 166         label edgeI = pEdges[pEdgeI];
 167
 168         if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.get(edgeI) == 1)
 169         {
 170             const edge& e = mesh_.edges()[edgeI];
 171
 172             edgeToDualPoint_[edgeI] = meshMod.addPoint
 173             (
 174                 e.centre(mesh_.points()),
 175                 pointI, // masterPoint
 176                 -1,     // zoneID
 177                 true    // inCell
 178             );
 179         }
 180     }
 181 }
 182
 183
 184 // Return true if point on face has same dual cells on both owner and neighbour
 185 // sides.
 186 bool Foam::meshDualiser::sameDualCell
 187 (
 188     const label faceI,
 189     const label pointI
 190 ) const
 191 {
 192     if (!mesh_.isInternalFace(faceI))
 193     {
 194         FatalErrorIn("sameDualCell(const label, const label)")
 195             << "face:" << faceI << " is not internal face."
 196             << abort(FatalError);
 197     }
 198
 199     label own = mesh_.faceOwner()[faceI];
 200     label nei = mesh_.faceNeighbour()[faceI];
 201
 202     return findDualCell(own, pointI) == findDualCell(nei, pointI);
 203 }
 204
 205
 206 Foam::label Foam::meshDualiser::addInternalFace
 207 (
 208     const label masterPointI,
 209     const label masterEdgeI,
 210     const label masterFaceI,
 211
 212     const bool edgeOrder,
 213     const label dualCell0,
 214     const label dualCell1,
 215     const DynamicList<label>& verts,
 216     directTopoChange& meshMod
 217 ) const
 218 {
 219     face newFace(verts);
 220
 221     if (edgeOrder != (dualCell0 < dualCell1))
 222     {
 223         reverse(newFace);
 224     }
 225
 226     if (debug)
 227     {
 228         pointField facePoints(meshMod.points(), newFace);
 229
 230         labelList oldToNew;
 231         pointField newPoints;
 232         bool hasMerged = mergePoints
 233         (
 234             facePoints,
 235             1E-6,
 236             false,
 237             oldToNew,
 238             newPoints
 239         );
 240
 241         if (hasMerged)
 242         {
 243             FatalErrorIn("addInternalFace(..)")
 244                 << "verts:" << verts << " newFace:" << newFace
 245                 << " face points:" << facePoints
 246                 << abort(FatalError);
 247         }
 248     }
 249
 250
 251     label zoneID = -1;
 252     bool zoneFlip = false;
 253     if (masterFaceI != -1)
 254     {
 255         zoneID = mesh_.faceZones().whichZone(masterFaceI);
 256
 257         if (zoneID != -1)
 258         {
 259             const faceZone& fZone = mesh_.faceZones()[zoneID];
 260
 261             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
 262         }
 263     }
 264
 265     label dualFaceI;
 266
 267     if (dualCell0 < dualCell1)
 268     {
 269         dualFaceI = meshMod.addFace
 270         (
 271             newFace,
 272             dualCell0,      // own
 273             dualCell1,      // nei
 274             masterPointI,   // masterPointID
 275             masterEdgeI,    // masterEdgeID
 276             masterFaceI,    // masterFaceID
 277             false,          // flipFaceFlux
 278             -1,             // patchID
 279             zoneID,         // zoneID
 280             zoneFlip        // zoneFlip
 281         );
 282
 283         //pointField dualPoints(meshMod.points());
 284         //vector n(newFace.normal(dualPoints));
 285         //n /= mag(n);
 286         //Pout<< "Generated internal dualFace:" << dualFaceI
 287         //    << " verts:" << newFace
 288         //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
 289         //    << " n:" << n
 290         //    << " between dualowner:" << dualCell0
 291         //    << " dualneigbour:" << dualCell1
 292         //    << endl;
 293     }
 294     else
 295     {
 296         dualFaceI = meshMod.addFace
 297         (
 298             newFace,
 299             dualCell1,      // own
 300             dualCell0,      // nei
 301             masterPointI,   // masterPointID
 302             masterEdgeI,    // masterEdgeID
 303             masterFaceI,    // masterFaceID
 304             false,          // flipFaceFlux
 305             -1,             // patchID
 306             zoneID,         // zoneID
 307             zoneFlip        // zoneFlip
 308         );
 309
 310         //pointField dualPoints(meshMod.points());
 311         //vector n(newFace.normal(dualPoints));
 312         //n /= mag(n);
 313         //Pout<< "Generated internal dualFace:" << dualFaceI
 314         //    << " verts:" << newFace
 315         //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
 316         //    << " n:" << n
 317         //    << " between dualowner:" << dualCell1
 318         //    << " dualneigbour:" << dualCell0
 319         //    << endl;
 320     }
 321     return dualFaceI;
 322 }
 323
 324
 325 Foam::label Foam::meshDualiser::addBoundaryFace
 326 (
 327     const label masterPointI,
 328     const label masterEdgeI,
 329     const label masterFaceI,
 330
 331     const label dualCellI,
 332     const label patchI,
 333     const DynamicList<label>& verts,
 334     directTopoChange& meshMod
 335 ) const
 336 {
 337     face newFace(verts);
 338
 339     label zoneID = -1;
 340     bool zoneFlip = false;
 341     if (masterFaceI != -1)
 342     {
 343         zoneID = mesh_.faceZones().whichZone(masterFaceI);
 344
 345         if (zoneID != -1)
 346         {
 347             const faceZone& fZone = mesh_.faceZones()[zoneID];
 348
 349             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
 350         }
 351     }
 352
 353     label dualFaceI = meshMod.addFace
 354     (
 355         newFace,
 356         dualCellI,      // own
 357         -1,             // nei
 358         masterPointI,   // masterPointID
 359         masterEdgeI,    // masterEdgeID
 360         masterFaceI,    // masterFaceID
 361         false,          // flipFaceFlux
 362         patchI,         // patchID
 363         zoneID,         // zoneID
 364         zoneFlip        // zoneFlip
 365     );
 366
 367     //pointField dualPoints(meshMod.points());
 368     //vector n(newFace.normal(dualPoints));
 369     //n /= mag(n);
 370     //Pout<< "Generated boundary dualFace:" << dualFaceI
 371     //    << " verts:" << newFace
 372     //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
 373     //    << " n:" << n
 374     //    << " on dualowner:" << dualCellI
 375     //    << endl;
 376     return dualFaceI;
 377 }
 378
 379
 380 // Walks around edgeI.
 381 // splitFace=true : creates multiple faces
 382 // splitFace=false: creates single face if same dual cells on both sides,
 383 //                  multiple faces otherwise.
 384 void Foam::meshDualiser::createFacesAroundEdge
 385 (
 386     const bool splitFace,
 387     const PackedBoolList& isBoundaryEdge,
 388     const label edgeI,
 389     const label startFaceI,
 390     directTopoChange& meshMod,
 391     boolList& doneEFaces
 392 ) const
 393 {
 394     const edge& e = mesh_.edges()[edgeI];
 395     const labelList& eFaces = mesh_.edgeFaces()[edgeI];
 396
 397     label fp = edgeFaceCirculator::getMinIndex
 398     (
 399         mesh_.faces()[startFaceI],
 400         e[0],
 401         e[1]
 402     );
 403
 404     edgeFaceCirculator ie
 405     (
 406         mesh_,
 407         startFaceI, // face
 408         true,       // ownerSide
 409         fp,         // fp
 410         isBoundaryEdge.get(edgeI) == 1  // isBoundaryEdge
 411     );
 412     ie.setCanonical();
 413
 414     bool edgeOrder = ie.sameOrder(e[0], e[1]);
 415     label startFaceLabel = ie.faceLabel();
 416
 417     //Pout<< "At edge:" << edgeI << " verts:" << e
 418     //    << " points:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
 419     //    << " started walking at face:" << ie.faceLabel()
 420     //    << " verts:" << mesh_.faces()[ie.faceLabel()]
 421     //    << " edgeOrder:" << edgeOrder
 422     //    << " in direction of cell:" << ie.cellLabel()
 423     //    << endl;
 424
 425     // Walk and collect face.
 426     DynamicList<label> verts(100);
 427
 428     if (edgeToDualPoint_[edgeI] != -1)
 429     {
 430         verts.append(edgeToDualPoint_[edgeI]);
 431     }
 432     if (faceToDualPoint_[ie.faceLabel()] != -1)
 433     {
 434         doneEFaces[findIndex(eFaces, ie.faceLabel())] = true;
 435         verts.append(faceToDualPoint_[ie.faceLabel()]);
 436     }
 437     if (cellToDualPoint_[ie.cellLabel()] != -1)
 438     {
 439         verts.append(cellToDualPoint_[ie.cellLabel()]);
 440     }
 441
 442     label currentDualCell0 = findDualCell(ie.cellLabel(), e[0]);
 443     label currentDualCell1 = findDualCell(ie.cellLabel(), e[1]);
 444
 445     ++ie;
 446
 447     while (true)
 448     {
 449         label faceI = ie.faceLabel();
 450
 451         // Mark face as visited.
 452         doneEFaces[findIndex(eFaces, faceI)] = true;
 453
 454         if (faceToDualPoint_[faceI] != -1)
 455         {
 456             verts.append(faceToDualPoint_[faceI]);
 457         }
 458
 459         label cellI = ie.cellLabel();
 460
 461         if (cellI == -1)
 462         {
 463             // At ending boundary face. We've stored the face point above
 464             // so this is the whole face.
 465             break;
 466         }
 467
 468
 469         label dualCell0 = findDualCell(cellI, e[0]);
 470         label dualCell1 = findDualCell(cellI, e[1]);
 471
 472         // Generate face. (always if splitFace=true; only if needed to
 473         // separate cells otherwise)
 474         if
 475         (
 476             splitFace
 477          || (
 478                 dualCell0 != currentDualCell0
 479              || dualCell1 != currentDualCell1
 480             )
 481         )
 482         {
 483             // Close current face.
 484             addInternalFace
 485             (
 486                 -1,         // masterPointI
 487                 edgeI,      // masterEdgeI
 488                 -1,         // masterFaceI
 489                 edgeOrder,
 490                 currentDualCell0,
 491                 currentDualCell1,
 492                 verts.shrink(),
 493                 meshMod
 494             );
 495
 496             // Restart
 497             currentDualCell0 = dualCell0;
 498             currentDualCell1 = dualCell1;
 499
 500             verts.clear();
 501             if (edgeToDualPoint_[edgeI] != -1)
 502             {
 503                 verts.append(edgeToDualPoint_[edgeI]);
 504             }
 505             if (faceToDualPoint_[faceI] != -1)
 506             {
 507                 verts.append(faceToDualPoint_[faceI]);
 508             }
 509         }
 510
 511         if (cellToDualPoint_[cellI] != -1)
 512         {
 513             verts.append(cellToDualPoint_[cellI]);
 514         }
 515
 516         ++ie;
 517
 518         if (ie == ie.end())
 519         {
 520             // Back at start face (for internal edge only). See if this needs
 521             // adding.
 522             if (isBoundaryEdge.get(edgeI) == 0)
 523             {
 524                 label startDual = faceToDualPoint_[startFaceLabel];
 525
 526                 if (startDual != -1 && findIndex(verts, startDual) == -1)
 527                 {
 528                     verts.append(startDual);
 529                 }
 530             }
 531             break;
 532         }
 533     }
 534
 535     verts.shrink();
 536     addInternalFace
 537     (
 538         -1,         // masterPointI
 539         edgeI,      // masterEdgeI
 540         -1,         // masterFaceI
 541         edgeOrder,
 542         currentDualCell0,
 543         currentDualCell1,
 544         verts,
 545         meshMod
 546     );
 547 }
 548
 549
 550 // Walks around circumference of faceI. Creates single face. Gets given
 551 // starting (feature) edge to start from. Returns ending edge. (all edges
 552 // in form of index in faceEdges)
 553 void Foam::meshDualiser::createFaceFromInternalFace
 554 (
 555     const label faceI,
 556     label& fp,
 557     directTopoChange& meshMod
 558 ) const
 559 {
 560     const face& f = mesh_.faces()[faceI];
 561     const labelList& fEdges = mesh_.faceEdges()[faceI];
 562     label own = mesh_.faceOwner()[faceI];
 563     label nei = mesh_.faceNeighbour()[faceI];
 564
 565     //Pout<< "createFaceFromInternalFace : At face:" << faceI
 566     //    << " verts:" << f
 567     //    << " points:" << UIndirectList<point>(mesh_.points(), f)()
 568     //    << " started walking at edge:" << fEdges[fp]
 569     //    << " verts:" << mesh_.edges()[fEdges[fp]]
 570     //    << endl;
 571
 572
 573     // Walk and collect face.
 574     DynamicList<label> verts(100);
 575
 576     verts.append(faceToDualPoint_[faceI]);
 577     verts.append(edgeToDualPoint_[fEdges[fp]]);
 578
 579     // Step to vertex after edge mid
 580     fp = f.fcIndex(fp);
 581
 582     // Get cells on either side of face at that point
 583     label currentDualCell0 = findDualCell(own, f[fp]);
 584     label currentDualCell1 = findDualCell(nei, f[fp]);
 585
 586     forAll(f, i)
 587     {
 588         // Check vertex
 589         if (pointToDualPoint_[f[fp]] != -1)
 590         {
 591             verts.append(pointToDualPoint_[f[fp]]);
 592         }
 593
 594         // Edge between fp and fp+1
 595         label edgeI = fEdges[fp];
 596
 597         if (edgeToDualPoint_[edgeI] != -1)
 598         {
 599             verts.append(edgeToDualPoint_[edgeI]);
 600         }
 601
 602         // Next vertex on edge
 603         label nextFp = f.fcIndex(fp);
 604
 605         // Get dual cells on nextFp to check whether face needs closing.
 606         label dualCell0 = findDualCell(own, f[nextFp]);
 607         label dualCell1 = findDualCell(nei, f[nextFp]);
 608
 609         if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
 610         {
 611             // Check: make sure that there is a midpoint on the edge.
 612             if (edgeToDualPoint_[edgeI] == -1)
 613             {
 614                 FatalErrorIn("createFacesFromInternalFace(..)")
 615                     << "face:" << faceI << " verts:" << f
 616                     << " points:" << UIndirectList<point>(mesh_.points(), f)()
 617                     << " no feature edge between " << f[fp]
 618                     << " and " << f[nextFp] << " although have different"
 619                     << " dual cells." << endl
 620                     << "point " << f[fp] << " has dual cells "
 621                     << currentDualCell0 << " and " << currentDualCell1
 622                     << " ; point "<< f[nextFp] << " has dual cells "
 623                     << dualCell0 << " and " << dualCell1
 624                     << abort(FatalError);
 625             }
 626
 627
 628             // Close current face.
 629             verts.shrink();
 630             addInternalFace
 631             (
 632                 -1,         // masterPointI
 633                 -1,         // masterEdgeI
 634                 faceI,      // masterFaceI
 635                 true,       // edgeOrder,
 636                 currentDualCell0,
 637                 currentDualCell1,
 638                 verts,
 639                 meshMod
 640             );
 641             break;
 642         }
 643
 644         fp = nextFp;
 645     }
 646 }
 647
 648
 649 // Given a point on a face converts the faces around the point.
 650 // (pointFaces()). Gets starting face and marks off visited faces in donePFaces.
 651 void Foam::meshDualiser::createFacesAroundBoundaryPoint
 652 (
 653     const label patchI,
 654     const label patchPointI,
 655     const label startFaceI,
 656     directTopoChange& meshMod,
 657     boolList& donePFaces            // pFaces visited
 658 ) const
 659 {
 660     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 661     const polyPatch& pp = patches[patchI];
 662     const labelList& pFaces = pp.pointFaces()[patchPointI];
 663     const labelList& own = mesh_.faceOwner();
 664
 665     label pointI = pp.meshPoints()[patchPointI];
 666
 667     if (pointToDualPoint_[pointI] == -1)
 668     {
 669         // Not a feature point. Loop over all connected
 670         // pointFaces.
 671
 672         // Starting face
 673         label faceI = startFaceI;
 674
 675         DynamicList<label> verts(4);
 676
 677         while (true)
 678         {
 679             label index = findIndex(pFaces, faceI-pp.start());
 680
 681             // Has face been visited already?
 682             if (donePFaces[index])
 683             {
 684                 break;
 685             }
 686             donePFaces[index] = true;
 687
 688             // Insert face centre
 689             verts.append(faceToDualPoint_[faceI]);
 690
 691             label dualCellI = findDualCell(own[faceI], pointI);
 692
 693             // Get the edge before the patchPointI
 694             const face& f = mesh_.faces()[faceI];
 695             label fp = findIndex(f, pointI);
 696             label prevFp = f.rcIndex(fp);
 697             label edgeI = mesh_.faceEdges()[faceI][prevFp];
 698
 699             if (edgeToDualPoint_[edgeI] != -1)
 700             {
 701                 verts.append(edgeToDualPoint_[edgeI]);
 702             }
 703
 704             // Get next boundary face (whilst staying on edge).
 705             edgeFaceCirculator circ
 706             (
 707                 mesh_,
 708                 faceI,
 709                 true,   // ownerSide
 710                 prevFp, // index of edge in face
 711                 true    // isBoundaryEdge
 712             );
 713
 714             do
 715             {
 716                 ++circ;
 717             }
 718             while (mesh_.isInternalFace(circ.faceLabel()));
 719
 720             // Step to next face
 721             faceI = circ.faceLabel();
 722
 723             if (faceI < pp.start() || faceI >= pp.start()+pp.size())
 724             {
 725                 FatalErrorIn
 726                 (
 727                     "meshDualiser::createFacesAroundBoundaryPoint(..)"
 728                 )   << "Walked from face on patch:" << patchI
 729                     << " to face:" << faceI
 730                     << " fc:" << mesh_.faceCentres()[faceI]
 731                     << " on patch:" << patches.whichPatch(faceI)
 732                     << abort(FatalError);
 733             }
 734
 735             // Check if different cell.
 736             if (dualCellI != findDualCell(own[faceI], pointI))
 737             {
 738                 FatalErrorIn
 739                 (
 740                     "meshDualiser::createFacesAroundBoundaryPoint(..)"
 741                 )   << "Different dual cells but no feature edge"
 742                     << " inbetween point:" << pointI
 743                     << " coord:" << mesh_.points()[pointI]
 744                     << abort(FatalError);
 745             }
 746         }
 747
 748         verts.shrink();
 749
 750         label dualCellI = findDualCell(own[faceI], pointI);
 751
 752         //Bit dodgy: create dualface from the last face (instead of from
 753         // the central point). This will also use the original faceZone to
 754         // put the new face (which might span multiple original faces) in.
 755
 756         addBoundaryFace
 757         (
 758             //pointI,     // masterPointI
 759             -1,         // masterPointI
 760             -1,         // masterEdgeI
 761             faceI,      // masterFaceI
 762             dualCellI,
 763             patchI,
 764             verts,
 765             meshMod
 766         );
 767     }
 768     else
 769     {
 770         label faceI = startFaceI;
 771
 772         // Storage for face
 773         DynamicList<label> verts(mesh_.faces()[faceI].size());
 774
 775         // Starting point.
 776         verts.append(pointToDualPoint_[pointI]);
 777
 778         // Find edge between pointI and next point on face.
 779         const labelList& fEdges = mesh_.faceEdges()[faceI];
 780         label nextEdgeI = fEdges[findIndex(mesh_.faces()[faceI], pointI)];
 781         if (edgeToDualPoint_[nextEdgeI] != -1)
 782         {
 783             verts.append(edgeToDualPoint_[nextEdgeI]);
 784         }
 785
 786         do
 787         {
 788             label index = findIndex(pFaces, faceI-pp.start());
 789
 790             // Has face been visited already?
 791             if (donePFaces[index])
 792             {
 793                 break;
 794             }
 795             donePFaces[index] = true;
 796
 797             // Face centre
 798             verts.append(faceToDualPoint_[faceI]);
 799
 800             // Find edge before pointI on faceI
 801             const labelList& fEdges = mesh_.faceEdges()[faceI];
 802             const face& f = mesh_.faces()[faceI];
 803             label prevFp = f.rcIndex(findIndex(f, pointI));
 804             label edgeI = fEdges[prevFp];
 805
 806             if (edgeToDualPoint_[edgeI] != -1)
 807             {
 808                 // Feature edge. Close any face so far. Note: uses face to
 809                 // create dualFace from. Could use pointI instead.
 810                 verts.append(edgeToDualPoint_[edgeI]);
 811                 addBoundaryFace
 812                 (
 813                     -1,     // masterPointI
 814                     -1,     // masterEdgeI
 815                     faceI,  // masterFaceI
 816                     findDualCell(own[faceI], pointI),
 817                     patchI,
 818                     verts.shrink(),
 819                     meshMod
 820                 );
 821                 verts.clear();
 822
 823                 verts.append(pointToDualPoint_[pointI]);
 824                 verts.append(edgeToDualPoint_[edgeI]);
 825             }
 826
 827             // Cross edgeI to next boundary face
 828             edgeFaceCirculator circ
 829             (
 830                 mesh_,
 831                 faceI,
 832                 true,   // ownerSide
 833                 prevFp, // index of edge in face
 834                 true    // isBoundaryEdge
 835             );
 836
 837             do
 838             {
 839                 ++circ;
 840             }
 841             while (mesh_.isInternalFace(circ.faceLabel()));
 842
 843             // Step to next face. Quit if not on same patch.
 844             faceI = circ.faceLabel();
 845         }
 846         while
 847         (
 848             faceI != startFaceI
 849          && faceI >= pp.start()
 850          && faceI < pp.start()+pp.size()
 851         );
 852
 853         if (verts.size() > 2)
 854         {
 855             // Note: face created from face, not from pointI
 856             addBoundaryFace
 857             (
 858                 -1,             // masterPointI
 859                 -1,             // masterEdgeI
 860                 startFaceI,     // masterFaceI
 861                 findDualCell(own[faceI], pointI),
 862                 patchI,
 863                 verts.shrink(),
 864                 meshMod
 865             );
 866         }
 867     }
 868 }
 869
 870
 871 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 872
 873 Foam::meshDualiser::meshDualiser(const polyMesh& mesh)
 874 :
 875     mesh_(mesh),
 876     pointToDualCells_(mesh_.nPoints()),
 877     pointToDualPoint_(mesh_.nPoints(), -1),
 878     cellToDualPoint_(mesh_.nCells()),
 879     faceToDualPoint_(mesh_.nFaces(), -1),
 880     edgeToDualPoint_(mesh_.nEdges(), -1)
 881 {}
 882
 883
 884 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 885
 886
 887 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 888
 889 void Foam::meshDualiser::setRefinement
 890 (
 891     const bool splitFace,
 892     const labelList& featureFaces,
 893     const labelList& featureEdges,
 894     const labelList& singleCellFeaturePoints,
 895     const labelList& multiCellFeaturePoints,
 896     directTopoChange& meshMod
 897 )
 898 {
 899     const labelList& own = mesh_.faceOwner();
 900     const labelList& nei = mesh_.faceNeighbour();
 901     const vectorField& cellCentres = mesh_.cellCentres();
 902
 903     // Mark boundary edges and points.
 904     // (Note: in 1.4.2 we can use the built-in mesh point ordering
 905     //  facility instead)
 906     PackedBoolList isBoundaryEdge(mesh_.nEdges());
 907     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 908     {
 909         const labelList& fEdges = mesh_.faceEdges()[faceI];
 910
 911         forAll(fEdges, i)
 912         {
 913             isBoundaryEdge.set(fEdges[i], 1);
 914         }
 915     }
 916
 917
 918     if (splitFace)
 919     {
 920         // This is a special mode where whenever we are walking around an edge
 921         // every area through a cell becomes a separate dualface. So two
 922         // dual cells will probably have more than one dualface between them!
 923         // This mode implies that
 924         // - all faces have to be feature faces since there has to be a
 925         //   dualpoint at the face centre.
 926         // - all edges have to be feature edges ,,
 927         boolList featureFaceSet(mesh_.nFaces(), false);
 928         forAll(featureFaces, i)
 929         {
 930             featureFaceSet[featureFaces[i]] = true;
 931         }
 932         label faceI = findIndex(featureFaceSet, false);
 933
 934         if (faceI != -1)
 935         {
 936             FatalErrorIn
 937             (
 938                 "meshDualiser::setRefinement"
 939                 "(const labelList&, const labelList&, const labelList&"
 940                 ", const labelList, directTopoChange&)"
 941             )   << "In split-face-mode (splitFace=true) but not all faces"
 942                 << " marked as feature faces." << endl
 943                 << "First conflicting face:" << faceI
 944                 << " centre:" << mesh_.faceCentres()[faceI]
 945                 << abort(FatalError);
 946         }
 947
 948         boolList featureEdgeSet(mesh_.nEdges(), false);
 949         forAll(featureEdges, i)
 950         {
 951             featureEdgeSet[featureEdges[i]] = true;
 952         }
 953         label edgeI = findIndex(featureEdgeSet, false);
 954
 955         if (edgeI != -1)
 956         {
 957             const edge& e = mesh_.edges()[edgeI];
 958             FatalErrorIn
 959             (
 960                 "meshDualiser::setRefinement"
 961                 "(const labelList&, const labelList&, const labelList&"
 962                 ", const labelList, directTopoChange&)"
 963             )   << "In split-face-mode (splitFace=true) but not all edges"
 964                 << " marked as feature edges." << endl
 965                 << "First conflicting edge:" << edgeI
 966                 << " verts:" << e
 967                 << " coords:" << mesh_.points()[e[0]] << mesh_.points()[e[1]]
 968                 << abort(FatalError);
 969         }
 970     }
 971     else
 972     {
 973         // Check that all boundary faces are feature faces.
 974
 975         boolList featureFaceSet(mesh_.nFaces(), false);
 976         forAll(featureFaces, i)
 977         {
 978             featureFaceSet[featureFaces[i]] = true;
 979         }
 980         for
 981         (
 982             label faceI = mesh_.nInternalFaces();
 983             faceI < mesh_.nFaces();
 984             faceI++
 985         )
 986         {
 987             if (!featureFaceSet[faceI])
 988             {
 989                 FatalErrorIn
 990                 (
 991                     "meshDualiser::setRefinement"
 992                     "(const labelList&, const labelList&, const labelList&"
 993                     ", const labelList, directTopoChange&)"
 994                 )   << "Not all boundary faces marked as feature faces."
 995                     << endl
 996                     << "First conflicting face:" << faceI
 997                     << " centre:" << mesh_.faceCentres()[faceI]
 998                     << abort(FatalError);
 999             }
1000         }
1001     }
1002
1003
1004
1005
1006     // Start creating cells, points, and faces (in that order)
1007
1008
1009     // 1. Mark which cells to create
1010     // Mostly every point becomes one cell but sometimes (for feature points)
1011     // all cells surrounding a feature point become cells. Also a non-manifold
1012     // point can create two cells! So a dual cell is uniquely defined by a
1013     // mesh point + cell (as in pointCells index)
1014
1015     // 2. Mark which face centres to create
1016
1017     // 3. Internal faces can now consist of
1018     //      - only cell centres of walk around edge
1019     //      - cell centres + face centres of walk around edge
1020     //      - same but now other side is not a single cell
1021
1022     // 4. Boundary faces (or internal faces between cell zones!) now consist of
1023     //      - walk around boundary point.
1024
1025
1026
1027     autoPtr<OFstream> dualCcStr;
1028     if (debug)
1029     {
1030         dualCcStr.reset(new OFstream("dualCc.obj"));
1031         Pout<< "Dumping centres of dual cells to " << dualCcStr().name()
1032             << endl;
1033     }
1034
1035
1036     // Dual cells (from points)
1037     // ~~~~~~~~~~~~~~~~~~~~~~~~
1038
1039     // pointToDualCells_[pointI]
1040     // - single entry : all cells surrounding point all become the same
1041     //                  cell.
1042     // - multiple entries: in order of pointCells.
1043
1044
1045     // feature points that become single cell
1046     forAll(singleCellFeaturePoints, i)
1047     {
1048         label pointI = singleCellFeaturePoints[i];
1049
1050         pointToDualPoint_[pointI] = meshMod.addPoint
1051         (
1052             mesh_.points()[pointI],
1053             pointI,                                 // masterPoint
1054             mesh_.pointZones().whichZone(pointI),   // zoneID
1055             true                                    // inCell
1056         );
1057
1058         // Generate single cell
1059         pointToDualCells_[pointI].setSize(1);
1060         pointToDualCells_[pointI][0] = meshMod.addCell
1061         (
1062             pointI, //masterPointID,
1063             -1,     //masterEdgeID,
1064             -1,     //masterFaceID,
1065             -1,     //masterCellID,
1066             -1      //zoneID
1067         );
1068         if (dualCcStr.valid())
1069         {
1070             meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
1071         }
1072     }
1073
1074     // feature points that become multiple cells
1075     forAll(multiCellFeaturePoints, i)
1076     {
1077         label pointI = multiCellFeaturePoints[i];
1078
1079         if (pointToDualCells_[pointI].size() > 0)
1080         {
1081             FatalErrorIn
1082             (
1083                 "meshDualiser::setRefinement"
1084                 "(const labelList&, const labelList&, const labelList&"
1085                 ", const labelList, directTopoChange&)"
1086             )   << "Point " << pointI << " at:" << mesh_.points()[pointI]
1087                 << " is both in singleCellFeaturePoints"
1088                 << " and multiCellFeaturePoints."
1089                 << abort(FatalError);
1090         }
1091
1092         pointToDualPoint_[pointI] = meshMod.addPoint
1093         (
1094             mesh_.points()[pointI],
1095             pointI,                                 // masterPoint
1096             mesh_.pointZones().whichZone(pointI),   // zoneID
1097             true                                    // inCell
1098         );
1099
1100         // Create dualcell for every cell connected to dual point
1101
1102         const labelList& pCells = mesh_.pointCells()[pointI];
1103
1104         pointToDualCells_[pointI].setSize(pCells.size());
1105
1106         forAll(pCells, pCellI)
1107         {
1108             pointToDualCells_[pointI][pCellI] = meshMod.addCell
1109             (
1110                 pointI,                                     //masterPointID
1111                 -1,                                         //masterEdgeID
1112                 -1,                                         //masterFaceID
1113                 -1,                                         //masterCellID
1114                 mesh_.cellZones().whichZone(pCells[pCellI]) //zoneID
1115             );
1116             if (dualCcStr.valid())
1117             {
1118                 meshTools::writeOBJ
1119                 (
1120                     dualCcStr(),
1121                     0.5*(mesh_.points()[pointI]+cellCentres[pCells[pCellI]])
1122                 );
1123             }
1124         }
1125     }
1126     // Normal points
1127     forAll(mesh_.points(), pointI)
1128     {
1129         if (pointToDualCells_[pointI].empty())
1130         {
1131             pointToDualCells_[pointI].setSize(1);
1132             pointToDualCells_[pointI][0] = meshMod.addCell
1133             (
1134                 pointI, //masterPointID,
1135                 -1,     //masterEdgeID,
1136                 -1,     //masterFaceID,
1137                 -1,     //masterCellID,
1138                 -1      //zoneID
1139             );
1140
1141             if (dualCcStr.valid())
1142             {
1143                 meshTools::writeOBJ(dualCcStr(), mesh_.points()[pointI]);
1144             }
1145         }
1146     }
1147
1148
1149     // Dual points (from cell centres, feature faces, feature edges)
1150     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1151
1152     forAll(cellToDualPoint_, cellI)
1153     {
1154         cellToDualPoint_[cellI] = meshMod.addPoint
1155         (
1156             cellCentres[cellI],
1157             mesh_.faces()[mesh_.cells()[cellI][0]][0],     // masterPoint
1158             -1,     // zoneID
1159             true    // inCell
1160         );
1161     }
1162
1163     // From face to dual point
1164
1165     forAll(featureFaces, i)
1166     {
1167         label faceI = featureFaces[i];
1168
1169         faceToDualPoint_[faceI] = meshMod.addPoint
1170         (
1171             mesh_.faceCentres()[faceI],
1172             mesh_.faces()[faceI][0],     // masterPoint
1173             -1,     // zoneID
1174             true    // inCell
1175         );
1176     }
1177     // Detect whether different dual cells on either side of a face. This
1178     // would neccesitate having a dual face built from the face and thus a
1179     // dual point at the face centre.
1180     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1181     {
1182         if (faceToDualPoint_[faceI] == -1)
1183         {
1184             const face& f = mesh_.faces()[faceI];
1185
1186             forAll(f, fp)
1187             {
1188                 label ownDualCell = findDualCell(own[faceI], f[fp]);
1189                 label neiDualCell = findDualCell(nei[faceI], f[fp]);
1190
1191                 if (ownDualCell != neiDualCell)
1192                 {
1193                     faceToDualPoint_[faceI] = meshMod.addPoint
1194                     (
1195                         mesh_.faceCentres()[faceI],
1196                         f[fp],  // masterPoint
1197                         -1,     // zoneID
1198                         true    // inCell
1199                     );
1200
1201                     break;
1202                 }
1203             }
1204         }
1205     }
1206
1207     // From edge to dual point
1208
1209     forAll(featureEdges, i)
1210     {
1211         label edgeI = featureEdges[i];
1212
1213         const edge& e = mesh_.edges()[edgeI];
1214
1215         edgeToDualPoint_[edgeI] = meshMod.addPoint
1216         (
1217             e.centre(mesh_.points()),
1218             e[0],   // masterPoint
1219             -1,     // zoneID
1220             true    // inCell
1221         );
1222     }
1223
1224     // Detect whether different dual cells on either side of an edge. This
1225     // would neccesitate having a dual face built perpendicular to the edge
1226     // and thus a dual point at the mid of the edge.
1227     // Note: not really true - the face can be built without the edge centre!
1228     const labelListList& edgeCells = mesh_.edgeCells();
1229
1230     forAll(edgeCells, edgeI)
1231     {
1232        if (edgeToDualPoint_[edgeI] == -1)
1233        {
1234             const edge& e = mesh_.edges()[edgeI];
1235
1236             // We need a point on the edge if not all cells on both sides
1237             // are the same.
1238
1239             const labelList& eCells = mesh_.edgeCells()[edgeI];
1240
1241             label dualE0 = findDualCell(eCells[0], e[0]);
1242             label dualE1 = findDualCell(eCells[0], e[1]);
1243
1244             for (label i = 1; i < eCells.size(); i++)
1245             {
1246                 label newDualE0 = findDualCell(eCells[i], e[0]);
1247
1248                 if (dualE0 != newDualE0)
1249                 {
1250                     edgeToDualPoint_[edgeI] = meshMod.addPoint
1251                     (
1252                         e.centre(mesh_.points()),
1253                         e[0],                               // masterPoint
1254                         mesh_.pointZones().whichZone(e[0]), // zoneID
1255                         true                                // inCell
1256                     );
1257
1258                     break;
1259                 }
1260
1261                 label newDualE1 = findDualCell(eCells[i], e[1]);
1262
1263                 if (dualE1 != newDualE1)
1264                 {
1265                     edgeToDualPoint_[edgeI] = meshMod.addPoint
1266                     (
1267                         e.centre(mesh_.points()),
1268                         e[1],   // masterPoint
1269                         mesh_.pointZones().whichZone(e[1]), // zoneID
1270                         true    // inCell
1271                     );
1272
1273                     break;
1274                 }
1275             }
1276         }
1277     }
1278
1279     // Make sure all boundary edges emanating from feature points are
1280     // feature edges as well.
1281     forAll(singleCellFeaturePoints, i)
1282     {
1283         generateDualBoundaryEdges
1284         (
1285             isBoundaryEdge,
1286             singleCellFeaturePoints[i],
1287             meshMod
1288         );
1289     }
1290     forAll(multiCellFeaturePoints, i)
1291     {
1292         generateDualBoundaryEdges
1293         (
1294             isBoundaryEdge,
1295             multiCellFeaturePoints[i],
1296             meshMod
1297         );
1298     }
1299
1300
1301     // Check for duplicate points
1302     if (debug)
1303     {
1304         dumpPolyTopoChange(meshMod, "generatedPoints_");
1305         checkPolyTopoChange(meshMod);
1306     }
1307
1308
1309     // Now we have all points and cells
1310     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1311     //  - pointToDualCells_ : per point a single dualCell or multiple dualCells
1312     //  - pointToDualPoint_ : per point -1 or the dual point at the coordinate
1313     //  - edgeToDualPoint_  : per edge -1 or the edge centre
1314     //  - faceToDualPoint_  : per face -1 or the face centre
1315     //  - cellToDualPoint_  : per cell the cell centre
1316     // Now we have to walk all edges and construct faces. Either single face
1317     // per edge or multiple (-if nonmanifold edge -if different dualcells)
1318
1319     const edgeList& edges = mesh_.edges();
1320
1321     forAll(edges, edgeI)
1322     {
1323         const labelList& eFaces = mesh_.edgeFaces()[edgeI];
1324
1325         boolList doneEFaces(eFaces.size(), false);
1326
1327         forAll(eFaces, i)
1328         {
1329             if (!doneEFaces[i])
1330             {
1331                 // We found a face that hasn't yet been visited. This might
1332                 // happen for non-manifold edges where a single edge can
1333                 // become multiple faces.
1334
1335                 label startFaceI = eFaces[i];
1336
1337                 //Pout<< "Walking edge:" << edgeI
1338                 //    << " points:" << mesh_.points()[e[0]]
1339                 //    << mesh_.points()[e[1]]
1340                 //    << " startFace:" << startFaceI
1341                 //    << " at:" << mesh_.faceCentres()[startFaceI]
1342                 //    << endl;
1343
1344                 createFacesAroundEdge
1345                 (
1346                     splitFace,
1347                     isBoundaryEdge,
1348                     edgeI,
1349                     startFaceI,
1350                     meshMod,
1351                     doneEFaces
1352                 );
1353             }
1354         }
1355     }
1356
1357     if (debug)
1358     {
1359         dumpPolyTopoChange(meshMod, "generatedFacesFromEdges_");
1360     }
1361
1362     // Create faces from feature faces. These can be internal or external faces.
1363     // - feature face : centre needs to be included.
1364     //      - single cells on either side: triangulate
1365     //      - multiple cells: create single face between unique cell pair. Only
1366     //                        create face where cells differ on either side.
1367     // - non-feature face : inbetween cell zones.
1368     forAll(faceToDualPoint_, faceI)
1369     {
1370         if (faceToDualPoint_[faceI] != -1 && mesh_.isInternalFace(faceI))
1371         {
1372             const face& f = mesh_.faces()[faceI];
1373             const labelList& fEdges = mesh_.faceEdges()[faceI];
1374
1375             // Starting edge
1376             label fp = 0;
1377
1378             do
1379             {
1380                 // Find edge that is in dual mesh and where the
1381                 // next point (fp+1) has different dual cells on either side.
1382                 bool foundStart = false;
1383
1384                 do
1385                 {
1386                     if
1387                     (
1388                         edgeToDualPoint_[fEdges[fp]] != -1
1389                     && !sameDualCell(faceI, f.nextLabel(fp))
1390                     )
1391                     {
1392                         foundStart = true;
1393                         break;
1394                     }
1395                     fp = f.fcIndex(fp);
1396                 }
1397                 while (fp != 0);
1398
1399                 if (!foundStart)
1400                 {
1401                     break;
1402                 }
1403
1404                 // Walk from edge fp and generate a face.
1405                 createFaceFromInternalFace
1406                 (
1407                     faceI,
1408                     fp,
1409                     meshMod
1410                 );
1411             }
1412             while(fp != 0);
1413         }
1414     }
1415
1416     if (debug)
1417     {
1418         dumpPolyTopoChange(meshMod, "generatedFacesFromFeatFaces_");
1419     }
1420
1421
1422     // Create boundary faces. Every boundary point has one or more dualcells.
1423     // These need to be closed.
1424     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1425
1426     forAll(patches, patchI)
1427     {
1428         const polyPatch& pp = patches[patchI];
1429
1430         const labelListList& pointFaces = pp.pointFaces();
1431
1432         forAll(pointFaces, patchPointI)
1433         {
1434             const labelList& pFaces = pointFaces[patchPointI];
1435
1436             boolList donePFaces(pFaces.size(), false);
1437
1438             forAll(pFaces, i)
1439             {
1440                 if (!donePFaces[i])
1441                 {
1442                     // Starting face
1443                     label startFaceI = pp.start()+pFaces[i];
1444
1445                     //Pout<< "Walking around point:" << pointI
1446                     //    << " coord:" << mesh_.points()[pointI]
1447                     //    << " on patch:" << patchI
1448                     //    << " startFace:" << startFaceI
1449                     //    << " at:" << mesh_.faceCentres()[startFaceI]
1450                     //    << endl;
1451
1452                     createFacesAroundBoundaryPoint
1453                     (
1454                         patchI,
1455                         patchPointI,
1456                         startFaceI,
1457                         meshMod,
1458                         donePFaces            // pFaces visited
1459                     );
1460                 }
1461             }
1462         }
1463     }
1464
1465     if (debug)
1466     {
1467         dumpPolyTopoChange(meshMod, "generatedFacesFromBndFaces_");
1468     }
1469 }
1470
1471
1472
1473 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1474
1475
1476 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
1477
1478
1479 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1480
1481
1482 // ************************************************************************* //