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

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