src/conversion/polyDualMesh/polyDualMesh.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 InClass
  25     polyDualMesh
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "polyDualMesh.H"
  30 #include "meshTools.H"
  31 #include "OFstream.H"
  32 #include "foamTime.H"
  33 #include "SortableList.H"
  34 #include "pointSet.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(polyDualMesh, 0);
  41 }
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 // Determine order for faces:
  46 // - upper-triangular order for internal faces
  47 // - external faces after internal faces (were already so)
  48 Foam::labelList Foam::polyDualMesh::getFaceOrder
  49 (
  50     const labelList& faceOwner,
  51     const labelList& faceNeighbour,
  52     const cellList& cells,
  53     label& nInternalFaces
  54 )
  55 {
  56     labelList oldToNew(faceOwner.size(), -1);
  57
  58     // First unassigned face
  59     label newFaceI = 0;
  60
  61     forAll(cells, cellI)
  62     {
  63         const labelList& cFaces = cells[cellI];
  64
  65         SortableList<label> nbr(cFaces.size());
  66
  67         forAll(cFaces, i)
  68         {
  69             label faceI = cFaces[i];
  70
  71             label nbrCellI = faceNeighbour[faceI];
  72
  73             if (nbrCellI != -1)
  74             {
  75                 // Internal face. Get cell on other side.
  76                 if (nbrCellI == cellI)
  77                 {
  78                     nbrCellI = faceOwner[faceI];
  79                 }
  80
  81                 if (cellI < nbrCellI)
  82                 {
  83                     // CellI is master
  84                     nbr[i] = nbrCellI;
  85                 }
  86                 else
  87                 {
  88                     // nbrCell is master. Let it handle this face.
  89                     nbr[i] = -1;
  90                 }
  91             }
  92             else
  93             {
  94                 // External face. Do later.
  95                 nbr[i] = -1;
  96             }
  97         }
  98
  99         nbr.sort();
 100
 101         forAll(nbr, i)
 102         {
 103             if (nbr[i] != -1)
 104             {
 105                 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
 106             }
 107         }
 108     }
 109
 110     nInternalFaces = newFaceI;
 111
 112     Pout<< "nInternalFaces:" << nInternalFaces << endl;
 113     Pout<< "nFaces:" << faceOwner.size() << endl;
 114     Pout<< "nCells:" << cells.size() << endl;
 115
 116     // Leave patch faces intact.
 117     for (label faceI = newFaceI; faceI < faceOwner.size(); faceI++)
 118     {
 119         oldToNew[faceI] = faceI;
 120     }
 121
 122
 123     // Check done all faces.
 124     forAll(oldToNew, faceI)
 125     {
 126         if (oldToNew[faceI] == -1)
 127         {
 128             FatalErrorIn
 129             (
 130                 "polyDualMesh::getFaceOrder"
 131                 "(const labelList&, const labelList&, const label) const"
 132             )   << "Did not determine new position"
 133                 << " for face " << faceI
 134                 << abort(FatalError);
 135         }
 136     }
 137
 138     return oldToNew;
 139 }
 140
 141
 142 // Get the two edges on faceI using pointI. Returns them such that the order
 143 // - otherVertex of e0
 144 // - pointI
 145 // - otherVertex(pointI) of e1
 146 // is in face order
 147 void Foam::polyDualMesh::getPointEdges
 148 (
 149     const primitivePatch& patch,
 150     const label faceI,
 151     const label pointI,
 152     label& e0,
 153     label& e1
 154 )
 155 {
 156     const labelList& fEdges = patch.faceEdges()[faceI];
 157     const face& f = patch.localFaces()[faceI];
 158
 159     e0 = -1;
 160     e1 = -1;
 161
 162     forAll(fEdges, i)
 163     {
 164         label edgeI = fEdges[i];
 165
 166         const edge& e = patch.edges()[edgeI];
 167
 168         if (e[0] == pointI)
 169         {
 170             // One of the edges using pointI. Check which one.
 171             label index = findIndex(f, pointI);
 172
 173             if (f.nextLabel(index) == e[1])
 174             {
 175                 e1 = edgeI;
 176             }
 177             else
 178             {
 179                 e0 = edgeI;
 180             }
 181
 182             if (e0 != -1 && e1 != -1)
 183             {
 184                 return;
 185             }
 186         }
 187         else if (e[1] == pointI)
 188         {
 189             // One of the edges using pointI. Check which one.
 190             label index = findIndex(f, pointI);
 191
 192             if (f.nextLabel(index) == e[0])
 193             {
 194                 e1 = edgeI;
 195             }
 196             else
 197             {
 198                 e0 = edgeI;
 199             }
 200
 201             if (e0 != -1 && e1 != -1)
 202             {
 203                 return;
 204             }
 205         }
 206     }
 207
 208     FatalErrorIn("getPointEdges") << "Cannot find two edges on face:" << faceI
 209         << " vertices:" << patch.localFaces()[faceI]
 210         << " that uses point:" << pointI
 211         << abort(FatalError);
 212 }
 213
 214
 215 // Collect the face on an external point of the patch.
 216 Foam::labelList Foam::polyDualMesh::collectPatchSideFace
 217 (
 218     const polyPatch& patch,
 219     const label patchToDualOffset,
 220     const labelList& edgeToDualPoint,
 221     const labelList& pointToDualPoint,
 222     const label pointI,
 223
 224     label& edgeI
 225 )
 226 {
 227     // Construct face by walking around e[eI] starting from
 228     // patchEdgeI
 229
 230     label meshPointI = patch.meshPoints()[pointI];
 231     const labelList& pFaces = patch.pointFaces()[pointI];
 232
 233     dynamicLabelList dualFace;
 234
 235     if (pointToDualPoint[meshPointI] >= 0)
 236     {
 237         // Number of pFaces + 2 boundary edge + feature point
 238         dualFace.setCapacity(pFaces.size()+2+1);
 239         // Store dualVertex for feature edge
 240         dualFace.append(pointToDualPoint[meshPointI]);
 241     }
 242     else
 243     {
 244         dualFace.setCapacity(pFaces.size()+2);
 245     }
 246
 247     // Store dual vertex for starting edge.
 248     if (edgeToDualPoint[patch.meshEdges()[edgeI]] < 0)
 249     {
 250         FatalErrorIn("polyDualMesh::collectPatchSideFace") << edgeI
 251             << abort(FatalError);
 252     }
 253
 254     dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
 255
 256     label faceI = patch.edgeFaces()[edgeI][0];
 257
 258     // Check order of vertices of edgeI in face to see if we need to reverse.
 259     bool reverseFace;
 260
 261     label e0, e1;
 262     getPointEdges(patch, faceI, pointI, e0, e1);
 263
 264     if (e0 == edgeI)
 265     {
 266         reverseFace = true;
 267     }
 268     else
 269     {
 270         reverseFace = false;
 271     }
 272
 273     while (true)
 274     {
 275         // Store dual vertex for faceI.
 276         dualFace.append(faceI + patchToDualOffset);
 277
 278         // Cross face to other edge on pointI
 279         label e0, e1;
 280         getPointEdges(patch, faceI, pointI, e0, e1);
 281
 282         if (e0 == edgeI)
 283         {
 284             edgeI = e1;
 285         }
 286         else
 287         {
 288             edgeI = e0;
 289         }
 290
 291         if (edgeToDualPoint[patch.meshEdges()[edgeI]] >= 0)
 292         {
 293             // Feature edge. Insert dual vertex for edge.
 294             dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
 295         }
 296
 297         const labelList& eFaces = patch.edgeFaces()[edgeI];
 298
 299         if (eFaces.size() == 1)
 300         {
 301             // Reached other edge of patch
 302             break;
 303         }
 304
 305         // Cross edge to other face.
 306         if (eFaces[0] == faceI)
 307         {
 308             faceI = eFaces[1];
 309         }
 310         else
 311         {
 312             faceI = eFaces[0];
 313         }
 314     }
 315
 316     dualFace.shrink();
 317
 318     if (reverseFace)
 319     {
 320         reverse(dualFace);
 321     }
 322
 323     return dualFace;
 324 }
 325
 326
 327 // Collect face around pointI which is not on the outside of the patch.
 328 // Returns the vertices of the face and the indices in these vertices of
 329 // any points which are on feature edges.
 330 void Foam::polyDualMesh::collectPatchInternalFace
 331 (
 332     const polyPatch& patch,
 333     const label patchToDualOffset,
 334     const labelList& edgeToDualPoint,
 335
 336     const label pointI,
 337     const label startEdgeI,
 338
 339     labelList& dualFace2,
 340     labelList& featEdgeIndices2
 341 )
 342 {
 343     // Construct face by walking around pointI starting from startEdgeI
 344     const labelList& meshEdges = patch.meshEdges();
 345     const labelList& pFaces = patch.pointFaces()[pointI];
 346
 347     // Vertices of dualFace
 348     dynamicLabelList dualFace(pFaces.size());
 349     // Indices in dualFace of vertices added because of feature edge
 350     dynamicLabelList featEdgeIndices(dualFace.size());
 351
 352
 353     label edgeI = startEdgeI;
 354     label faceI = patch.edgeFaces()[edgeI][0];
 355
 356     // Check order of vertices of edgeI in face to see if we need to reverse.
 357     bool reverseFace;
 358
 359     label e0, e1;
 360     getPointEdges(patch, faceI, pointI, e0, e1);
 361
 362     if (e0 == edgeI)
 363     {
 364         reverseFace = true;
 365     }
 366     else
 367     {
 368         reverseFace = false;
 369     }
 370
 371     while (true)
 372     {
 373         // Insert dual vertex for face
 374         dualFace.append(faceI + patchToDualOffset);
 375
 376         // Cross face to other edge on pointI
 377         label e0, e1;
 378         getPointEdges(patch, faceI, pointI, e0, e1);
 379
 380         if (e0 == edgeI)
 381         {
 382             edgeI = e1;
 383         }
 384         else
 385         {
 386             edgeI = e0;
 387         }
 388
 389         if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
 390         {
 391             // Feature edge. Insert dual vertex for edge.
 392             dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
 393             featEdgeIndices.append(dualFace.size()-1);
 394         }
 395
 396         if (edgeI == startEdgeI)
 397         {
 398             break;
 399         }
 400
 401         // Cross edge to other face.
 402         const labelList& eFaces = patch.edgeFaces()[edgeI];
 403
 404         if (eFaces[0] == faceI)
 405         {
 406             faceI = eFaces[1];
 407         }
 408         else
 409         {
 410             faceI = eFaces[0];
 411         }
 412     }
 413
 414     dualFace2.transfer(dualFace);
 415
 416     featEdgeIndices2.transfer(featEdgeIndices);
 417
 418     if (reverseFace)
 419     {
 420         reverse(dualFace2);
 421
 422         // Correct featEdgeIndices for change in dualFace2
 423         forAll(featEdgeIndices2, i)
 424         {
 425             featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
 426         }
 427         // Reverse indices (might not be nessecary but do anyway)
 428         reverse(featEdgeIndices2);
 429     }
 430 }
 431
 432
 433 void Foam::polyDualMesh::splitFace
 434 (
 435     const polyPatch& patch,
 436     const labelList& pointToDualPoint,
 437
 438     const label pointI,
 439     const labelList& dualFace,
 440     const labelList& featEdgeIndices,
 441
 442     DynamicList<face>& dualFaces,
 443     dynamicLabelList& dualOwner,
 444     dynamicLabelList& dualNeighbour,
 445     dynamicLabelList& dualRegion
 446 )
 447 {
 448
 449     // Split face because of feature edges/points
 450     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 451     label meshPointI = patch.meshPoints()[pointI];
 452
 453     if (pointToDualPoint[meshPointI] >= 0)
 454     {
 455         // Feature point. Do face-centre decomposition.
 456
 457         if (featEdgeIndices.size() < 2)
 458         {
 459             // Feature point but no feature edges. Not handled at the moment
 460             dualFaces.append(face(dualFace));
 461             dualOwner.append(meshPointI);
 462             dualNeighbour.append(-1);
 463             dualRegion.append(patch.index());
 464         }
 465         else
 466         {
 467             // Do 'face-centre' decomposition. Start from first feature
 468             // edge create face up until next feature edge.
 469
 470             forAll(featEdgeIndices, i)
 471             {
 472                 label startFp = featEdgeIndices[i];
 473
 474                 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
 475
 476                 // Collect face from startFp to endFp
 477                 label sz = 0;
 478
 479                 if (endFp > startFp)
 480                 {
 481                     sz = endFp - startFp + 2;
 482                 }
 483                 else
 484                 {
 485                     sz = endFp + dualFace.size() - startFp + 2;
 486                 }
 487                 face subFace(sz);
 488
 489                 // feature point becomes face centre.
 490                 subFace[0] = pointToDualPoint[patch.meshPoints()[pointI]];
 491
 492                 // Copy from startFp up to endFp.
 493                 for (label subFp = 1; subFp < subFace.size(); subFp++)
 494                 {
 495                     subFace[subFp] = dualFace[startFp];
 496
 497                     startFp = (startFp+1) % dualFace.size();
 498                 }
 499
 500                 dualFaces.append(face(subFace));
 501                 dualOwner.append(meshPointI);
 502                 dualNeighbour.append(-1);
 503                 dualRegion.append(patch.index());
 504             }
 505         }
 506     }
 507     else
 508     {
 509         // No feature point. Check if feature edges.
 510         if (featEdgeIndices.size() < 2)
 511         {
 512             // Not enough feature edges. No split.
 513             dualFaces.append(face(dualFace));
 514             dualOwner.append(meshPointI);
 515             dualNeighbour.append(-1);
 516             dualRegion.append(patch.index());
 517         }
 518         else
 519         {
 520             // Multiple feature edges but no feature point.
 521             // Split face along feature edges. Gives weird result if
 522             // number of feature edges > 2.
 523
 524             // Storage for new face
 525             dynamicLabelList subFace(dualFace.size());
 526
 527             forAll(featEdgeIndices, featI)
 528             {
 529                 label startFp = featEdgeIndices[featI];
 530                 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
 531
 532                 label fp = startFp;
 533
 534                 while (true)
 535                 {
 536                     label vertI = dualFace[fp];
 537
 538                     subFace.append(vertI);
 539
 540                     if (fp == endFp)
 541                     {
 542                         break;
 543                     }
 544
 545                     fp = dualFace.fcIndex(fp);
 546                 }
 547
 548                 if (subFace.size() > 2)
 549                 {
 550                     // Enough vertices to create a face from.
 551                     subFace.shrink();
 552
 553                     dualFaces.append(face(subFace));
 554                     dualOwner.append(meshPointI);
 555                     dualNeighbour.append(-1);
 556                     dualRegion.append(patch.index());
 557
 558                     subFace.clear();
 559                 }
 560             }
 561             // Check final face.
 562             if (subFace.size() > 2)
 563             {
 564                 // Enough vertices to create a face from.
 565                 subFace.shrink();
 566
 567                 dualFaces.append(face(subFace));
 568                 dualOwner.append(meshPointI);
 569                 dualNeighbour.append(-1);
 570                 dualRegion.append(patch.index());
 571
 572                 subFace.clear();
 573             }
 574         }
 575     }
 576 }
 577
 578
 579 // Create boundary face for every point in patch
 580 void Foam::polyDualMesh::dualPatch
 581 (
 582     const polyPatch& patch,
 583     const label patchToDualOffset,
 584     const labelList& edgeToDualPoint,
 585     const labelList& pointToDualPoint,
 586
 587     const pointField& dualPoints,
 588
 589     DynamicList<face>& dualFaces,
 590     dynamicLabelList& dualOwner,
 591     dynamicLabelList& dualNeighbour,
 592     dynamicLabelList& dualRegion
 593 )
 594 {
 595     const labelList& meshEdges = patch.meshEdges();
 596
 597     // Whether edge has been done.
 598     // 0 : not
 599     // 1 : done e.start()
 600     // 2 : done e.end()
 601     // 3 : done both
 602     labelList doneEdgeSide(meshEdges.size(), 0);
 603
 604     boolList donePoint(patch.nPoints(), false);
 605
 606
 607     // Do points on edge of patch
 608     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 609
 610     forAll(doneEdgeSide, patchEdgeI)
 611     {
 612         const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
 613
 614         if (eFaces.size() == 1)
 615         {
 616             const edge& e = patch.edges()[patchEdgeI];
 617
 618             forAll(e, eI)
 619             {
 620                 label bitMask = 1<<eI;
 621
 622                 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
 623                 {
 624                     // Construct face by walking around e[eI] starting from
 625                     // patchEdgeI
 626
 627                     label pointI = e[eI];
 628
 629                     label edgeI = patchEdgeI;
 630                     labelList dualFace
 631                     (
 632                         collectPatchSideFace
 633                         (
 634                             patch,
 635                             patchToDualOffset,
 636                             edgeToDualPoint,
 637                             pointToDualPoint,
 638
 639                             pointI,
 640                             edgeI
 641                         )
 642                     );
 643
 644                     // Now edgeI is end edge. Mark as visited
 645                     if (patch.edges()[edgeI][0] == pointI)
 646                     {
 647                         doneEdgeSide[edgeI] |= 1;
 648                     }
 649                     else
 650                     {
 651                         doneEdgeSide[edgeI] |= 2;
 652                     }
 653
 654                     dualFaces.append(face(dualFace));
 655                     dualOwner.append(patch.meshPoints()[pointI]);
 656                     dualNeighbour.append(-1);
 657                     dualRegion.append(patch.index());
 658
 659                     doneEdgeSide[patchEdgeI] |= bitMask;
 660                     donePoint[pointI] = true;
 661                 }
 662             }
 663         }
 664     }
 665
 666
 667
 668     // Do patch-internal points
 669     // ~~~~~~~~~~~~~~~~~~~~~~~~
 670
 671     forAll(donePoint, pointI)
 672     {
 673         if (!donePoint[pointI])
 674         {
 675             labelList dualFace, featEdgeIndices;
 676
 677             collectPatchInternalFace
 678             (
 679                 patch,
 680                 patchToDualOffset,
 681                 edgeToDualPoint,
 682                 pointI,
 683                 patch.pointEdges()[pointI][0],  // Arbitrary starting edge
 684
 685                 dualFace,
 686                 featEdgeIndices
 687             );
 688
 689             //- Either keep in one piece or split face according to feature.
 690
 691             //// Keep face in one piece.
 692             //dualFaces.append(face(dualFace));
 693             //dualOwner.append(patch.meshPoints()[pointI]);
 694             //dualNeighbour.append(-1);
 695             //dualRegion.append(patch.index());
 696
 697             splitFace
 698             (
 699                 patch,
 700                 pointToDualPoint,
 701                 pointI,
 702                 dualFace,
 703                 featEdgeIndices,
 704
 705                 dualFaces,
 706                 dualOwner,
 707                 dualNeighbour,
 708                 dualRegion
 709             );
 710
 711             donePoint[pointI] = true;
 712         }
 713     }
 714 }
 715
 716
 717 void Foam::polyDualMesh::calcDual
 718 (
 719     const polyMesh& mesh,
 720     const labelList& featureEdges,
 721     const labelList& featurePoints
 722 )
 723 {
 724     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 725     const label nIntFaces = mesh.nInternalFaces();
 726
 727
 728     // Get patch for all of outside
 729     primitivePatch allBoundary
 730     (
 731         SubList<face>
 732         (
 733             mesh.faces(),
 734             mesh.nFaces() - nIntFaces,
 735             nIntFaces
 736         ),
 737         mesh.points()
 738     );
 739
 740     // Correspondence from patch edge to mesh edge.
 741     labelList meshEdges
 742     (
 743         allBoundary.meshEdges
 744         (
 745             mesh.edges(),
 746             mesh.pointEdges()
 747         )
 748     );
 749
 750     {
 751         pointSet nonManifoldPoints
 752         (
 753             mesh,
 754             "nonManifoldPoints",
 755             meshEdges.size() / 100
 756         );
 757
 758         allBoundary.checkPointManifold(true, &nonManifoldPoints);
 759
 760         if (nonManifoldPoints.size())
 761         {
 762             nonManifoldPoints.write();
 763
 764             FatalErrorIn
 765             (
 766                 "polyDualMesh::calcDual(const polyMesh&, const labelList&"
 767                 ", const labelList&)"
 768             )   << "There are " << nonManifoldPoints.size() << " points where"
 769                 << " the outside of the mesh is non-manifold." << nl
 770                 << "Such a mesh cannot be converted to a dual." << nl
 771                 << "Writing points at non-manifold sites to pointSet "
 772                 << nonManifoldPoints.name()
 773                 << exit(FatalError);
 774         }
 775     }
 776
 777
 778     // Assign points
 779     // ~~~~~~~~~~~~~
 780
 781     // mesh label   dualMesh vertex
 782     // ----------   ---------------
 783     // cellI        cellI
 784     // faceI        nCells+faceI-nIntFaces
 785     // featEdgeI    nCells+nFaces-nIntFaces+featEdgeI
 786     // featPointI   nCells+nFaces-nIntFaces+nFeatEdges+featPointI
 787
 788     pointField dualPoints
 789     (
 790         mesh.nCells()                           // cell centres
 791       + mesh.nFaces() - nIntFaces               // boundary face centres
 792       + featureEdges.size()                     // additional boundary edges
 793       + featurePoints.size()                    // additional boundary points
 794     );
 795
 796     label dualPointI = 0;
 797
 798
 799     // Cell centres.
 800     const pointField& cellCentres = mesh.cellCentres();
 801
 802     cellPoint_.setSize(cellCentres.size());
 803
 804     forAll(cellCentres, cellI)
 805     {
 806         cellPoint_[cellI] = dualPointI;
 807         dualPoints[dualPointI++] = cellCentres[cellI];
 808     }
 809
 810     // Boundary faces centres
 811     const pointField& faceCentres = mesh.faceCentres();
 812
 813     boundaryFacePoint_.setSize(mesh.nFaces() - nIntFaces);
 814
 815     for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
 816     {
 817         boundaryFacePoint_[faceI - nIntFaces] = dualPointI;
 818         dualPoints[dualPointI++] = faceCentres[faceI];
 819     }
 820
 821     // Edge status:
 822     //  >0 : dual point label (edge is feature edge)
 823     //  -1 : is boundary edge.
 824     //  -2 : is internal edge.
 825     labelList edgeToDualPoint(mesh.nEdges(), -2);
 826
 827     forAll(meshEdges, patchEdgeI)
 828     {
 829         label edgeI = meshEdges[patchEdgeI];
 830
 831         edgeToDualPoint[edgeI] = -1;
 832     }
 833
 834     forAll(featureEdges, i)
 835     {
 836         label edgeI = featureEdges[i];
 837
 838         const edge& e = mesh.edges()[edgeI];
 839
 840         edgeToDualPoint[edgeI] = dualPointI;
 841         dualPoints[dualPointI++] = e.centre(mesh.points());
 842     }
 843
 844
 845
 846     // Point status:
 847     //  >0 : dual point label (point is feature point)
 848     //  -1 : is point on edge of patch
 849     //  -2 : is point on patch (but not on edge)
 850     //  -3 : is internal point.
 851     labelList pointToDualPoint(mesh.nPoints(), -3);
 852
 853     forAll(patches, patchI)
 854     {
 855         const labelList& meshPoints = patches[patchI].meshPoints();
 856
 857         forAll(meshPoints, i)
 858         {
 859             pointToDualPoint[meshPoints[i]] = -2;
 860         }
 861
 862         const labelListList& loops = patches[patchI].edgeLoops();
 863
 864         forAll(loops, i)
 865         {
 866             const labelList& loop = loops[i];
 867
 868             forAll(loop, j)
 869             {
 870                  pointToDualPoint[meshPoints[loop[j]]] = -1;
 871             }
 872         }
 873     }
 874
 875     forAll(featurePoints, i)
 876     {
 877         label pointI = featurePoints[i];
 878
 879         pointToDualPoint[pointI] = dualPointI;
 880         dualPoints[dualPointI++] = mesh.points()[pointI];
 881     }
 882
 883
 884     // Storage for new faces.
 885     // Dynamic sized since we don't know size.
 886
 887     DynamicList<face> dynDualFaces(mesh.nEdges());
 888     dynamicLabelList dynDualOwner(mesh.nEdges());
 889     dynamicLabelList dynDualNeighbour(mesh.nEdges());
 890     dynamicLabelList dynDualRegion(mesh.nEdges());
 891
 892
 893     // Generate faces from edges on the boundary
 894     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 895
 896     forAll(meshEdges, patchEdgeI)
 897     {
 898         label edgeI = meshEdges[patchEdgeI];
 899
 900         const edge& e = mesh.edges()[edgeI];
 901
 902         label owner = -1;
 903         label neighbour = -1;
 904
 905         if (e[0] < e[1])
 906         {
 907             owner = e[0];
 908             neighbour = e[1];
 909         }
 910         else
 911         {
 912             owner = e[1];
 913             neighbour = e[0];
 914         }
 915
 916         // Find the boundary faces using the edge.
 917         const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
 918
 919         if (patchFaces.size() != 2)
 920         {
 921             FatalErrorIn("polyDualMesh::calcDual")
 922                 << "Cannot handle edges with " << patchFaces.size()
 923                 << " connected boundary faces."
 924                 << abort(FatalError);
 925         }
 926
 927         label face0 = patchFaces[0] + nIntFaces;
 928         const face& f0 = mesh.faces()[face0];
 929
 930         label face1 = patchFaces[1] + nIntFaces;
 931
 932         // We want to start walking from patchFaces[0] or patchFaces[1],
 933         // depending on which one uses owner,neighbour in the right order.
 934
 935         label startFaceI = -1;
 936         label endFaceI = -1;
 937
 938         label index = findIndex(f0, neighbour);
 939
 940         if (f0.nextLabel(index) == owner)
 941         {
 942             startFaceI = face0;
 943             endFaceI = face1;
 944         }
 945         else
 946         {
 947             startFaceI = face1;
 948             endFaceI = face0;
 949         }
 950
 951         // Now walk from startFaceI to cell to face to cell etc. until endFaceI
 952
 953         dynamicLabelList dualFace;
 954
 955         if (edgeToDualPoint[edgeI] >= 0)
 956         {
 957             // Number of cells + 2 boundary faces + feature edge point
 958             dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
 959             // Store dualVertex for feature edge
 960             dualFace.append(edgeToDualPoint[edgeI]);
 961         }
 962         else
 963         {
 964             dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
 965         }
 966
 967         // Store dual vertex for starting face.
 968         dualFace.append(mesh.nCells() + startFaceI - nIntFaces);
 969
 970         label cellI = mesh.faceOwner()[startFaceI];
 971         label faceI = startFaceI;
 972
 973         while (true)
 974         {
 975             // Store dual vertex from cellI.
 976             dualFace.append(cellI);
 977
 978             // Cross cell to other face on edge.
 979             label f0, f1;
 980             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
 981
 982             if (f0 == faceI)
 983             {
 984                 faceI = f1;
 985             }
 986             else
 987             {
 988                 faceI = f0;
 989             }
 990
 991             // Cross face to other cell.
 992             if (faceI == endFaceI)
 993             {
 994                 break;
 995             }
 996
 997             if (mesh.faceOwner()[faceI] == cellI)
 998             {
 999                 cellI = mesh.faceNeighbour()[faceI];
1000             }
1001             else
1002             {
1003                 cellI = mesh.faceOwner()[faceI];
1004             }
1005         }
1006
1007         // Store dual vertex for endFace.
1008         dualFace.append(mesh.nCells() + endFaceI - nIntFaces);
1009
1010         dynDualFaces.append(face(dualFace.shrink()));
1011         dynDualOwner.append(owner);
1012         dynDualNeighbour.append(neighbour);
1013         dynDualRegion.append(-1);
1014
1015         {
1016             // Check orientation.
1017             const face& f = dynDualFaces[dynDualFaces.size()-1];
1018             vector n = f.normal(dualPoints);
1019             if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1020             {
1021                 WarningIn("calcDual") << "Incorrect orientation"
1022                     << " on boundary edge:" << edgeI
1023                     << mesh.points()[mesh.edges()[edgeI][0]]
1024                     << mesh.points()[mesh.edges()[edgeI][1]]
1025                     << endl;
1026             }
1027         }
1028     }
1029
1030
1031     // Generate faces from internal edges
1032     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1033
1034     forAll(edgeToDualPoint, edgeI)
1035     {
1036         if (edgeToDualPoint[edgeI] == -2)
1037         {
1038             // Internal edge.
1039
1040             // Find dual owner, neighbour
1041
1042             const edge& e = mesh.edges()[edgeI];
1043
1044             label owner = -1;
1045             label neighbour = -1;
1046
1047             if (e[0] < e[1])
1048             {
1049                 owner = e[0];
1050                 neighbour = e[1];
1051             }
1052             else
1053             {
1054                 owner = e[1];
1055                 neighbour = e[0];
1056             }
1057
1058             // Get a starting cell
1059             const labelList& eCells = mesh.edgeCells()[edgeI];
1060
1061             label cellI = eCells[0];
1062
1063             // Get the two faces on the cell and edge.
1064             label face0, face1;
1065             meshTools::getEdgeFaces(mesh, cellI, edgeI, face0, face1);
1066
1067             // Find the starting face by looking at the order in which
1068             // the face uses the owner, neighbour
1069             const face& f0 = mesh.faces()[face0];
1070
1071             label index = findIndex(f0, neighbour);
1072
1073             bool f0OrderOk = (f0.nextLabel(index) == owner);
1074
1075             label startFaceI = -1;
1076
1077             if (f0OrderOk == (mesh.faceOwner()[face0] == cellI))
1078             {
1079                 startFaceI = face0;
1080             }
1081             else
1082             {
1083                 startFaceI = face1;
1084             }
1085
1086             // Walk face-cell-face until starting face reached.
1087             dynamicLabelList dualFace(mesh.edgeCells()[edgeI].size());
1088
1089             label faceI = startFaceI;
1090
1091             while (true)
1092             {
1093                 // Store dual vertex from cellI.
1094                 dualFace.append(cellI);
1095
1096                 // Cross cell to other face on edge.
1097                 label f0, f1;
1098                 meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
1099
1100                 if (f0 == faceI)
1101                 {
1102                     faceI = f1;
1103                 }
1104                 else
1105                 {
1106                     faceI = f0;
1107                 }
1108
1109                 // Cross face to other cell.
1110                 if (faceI == startFaceI)
1111                 {
1112                     break;
1113                 }
1114
1115                 if (mesh.faceOwner()[faceI] == cellI)
1116                 {
1117                     cellI = mesh.faceNeighbour()[faceI];
1118                 }
1119                 else
1120                 {
1121                     cellI = mesh.faceOwner()[faceI];
1122                 }
1123             }
1124
1125             dynDualFaces.append(face(dualFace.shrink()));
1126             dynDualOwner.append(owner);
1127             dynDualNeighbour.append(neighbour);
1128             dynDualRegion.append(-1);
1129
1130             {
1131                 // Check orientation.
1132                 const face& f = dynDualFaces[dynDualFaces.size()-1];
1133                 vector n = f.normal(dualPoints);
1134                 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1135                 {
1136                     WarningIn("calcDual") << "Incorrect orientation"
1137                         << " on internal edge:" << edgeI
1138                         << mesh.points()[mesh.edges()[edgeI][0]]
1139                         << mesh.points()[mesh.edges()[edgeI][1]]
1140                         << endl;
1141                 }
1142             }
1143         }
1144     }
1145
1146     // Dump faces.
1147     if (debug)
1148     {
1149         dynDualFaces.shrink();
1150         dynDualOwner.shrink();
1151         dynDualNeighbour.shrink();
1152         dynDualRegion.shrink();
1153
1154         OFstream str("dualInternalFaces.obj");
1155         Pout<< "polyDualMesh::calcDual : dumping internal faces to "
1156             << str.name() << endl;
1157
1158         forAll(dualPoints, dualPointI)
1159         {
1160             meshTools::writeOBJ(str, dualPoints[dualPointI]);
1161         }
1162         forAll(dynDualFaces, dualFaceI)
1163         {
1164             const face& f = dynDualFaces[dualFaceI];
1165
1166             str<< 'f';
1167             forAll(f, fp)
1168             {
1169                 str<< ' ' << f[fp]+1;
1170             }
1171             str<< nl;
1172         }
1173     }
1174
1175     const label nInternalFaces = dynDualFaces.size();
1176
1177     // Outside faces
1178     // ~~~~~~~~~~~~~
1179
1180     forAll(patches, patchI)
1181     {
1182         const polyPatch& pp = patches[patchI];
1183
1184         dualPatch
1185         (
1186             pp,
1187             mesh.nCells() + pp.start() - nIntFaces,
1188             edgeToDualPoint,
1189             pointToDualPoint,
1190
1191             dualPoints,
1192
1193             dynDualFaces,
1194             dynDualOwner,
1195             dynDualNeighbour,
1196             dynDualRegion
1197         );
1198     }
1199
1200
1201     // Transfer face info to straight lists
1202     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1203     faceList dualFaces(dynDualFaces.shrink(), true);
1204     dynDualFaces.clear();
1205
1206     labelList dualOwner(dynDualOwner.shrink(), true);
1207     dynDualOwner.clear();
1208
1209     labelList dualNeighbour(dynDualNeighbour.shrink(), true);
1210     dynDualNeighbour.clear();
1211
1212     labelList dualRegion(dynDualRegion.shrink(), true);
1213     dynDualRegion.clear();
1214
1215
1216
1217     // Dump faces.
1218     if (debug)
1219     {
1220         OFstream str("dualFaces.obj");
1221         Pout<< "polyDualMesh::calcDual : dumping all faces to "
1222             << str.name() << endl;
1223
1224         forAll(dualPoints, dualPointI)
1225         {
1226             meshTools::writeOBJ(str, dualPoints[dualPointI]);
1227         }
1228         forAll(dualFaces, dualFaceI)
1229         {
1230             const face& f = dualFaces[dualFaceI];
1231
1232             str<< 'f';
1233             forAll(f, fp)
1234             {
1235                 str<< ' ' << f[fp]+1;
1236             }
1237             str<< nl;
1238         }
1239     }
1240
1241     // Create cells.
1242     cellList dualCells(mesh.nPoints());
1243
1244     forAll(dualCells, cellI)
1245     {
1246         dualCells[cellI].setSize(0);
1247     }
1248
1249     forAll(dualOwner, faceI)
1250     {
1251         label cellI = dualOwner[faceI];
1252
1253         labelList& cFaces = dualCells[cellI];
1254
1255         label sz = cFaces.size();
1256         cFaces.setSize(sz+1);
1257         cFaces[sz] = faceI;
1258     }
1259     forAll(dualNeighbour, faceI)
1260     {
1261         label cellI = dualNeighbour[faceI];
1262
1263         if (cellI != -1)
1264         {
1265             labelList& cFaces = dualCells[cellI];
1266
1267             label sz = cFaces.size();
1268             cFaces.setSize(sz+1);
1269             cFaces[sz] = faceI;
1270         }
1271     }
1272
1273
1274     // Do upper-triangular face ordering. Determines face reordering map and
1275     // number of internal faces.
1276     label dummy;
1277
1278     labelList oldToNew
1279     (
1280         getFaceOrder
1281         (
1282             dualOwner,
1283             dualNeighbour,
1284             dualCells,
1285             dummy
1286         )
1287     );
1288
1289     // Reorder faces.
1290     inplaceReorder(oldToNew, dualFaces);
1291     inplaceReorder(oldToNew, dualOwner);
1292     inplaceReorder(oldToNew, dualNeighbour);
1293     inplaceReorder(oldToNew, dualRegion);
1294     forAll(dualCells, cellI)
1295     {
1296         inplaceRenumber(oldToNew, dualCells[cellI]);
1297     }
1298
1299
1300     // Create patches
1301     labelList patchSizes(patches.size(), 0);
1302
1303     forAll(dualRegion, faceI)
1304     {
1305         if (dualRegion[faceI] >= 0)
1306         {
1307             patchSizes[dualRegion[faceI]]++;
1308         }
1309     }
1310
1311     labelList patchStarts(patches.size(), 0);
1312
1313     label faceI = nInternalFaces;
1314
1315     forAll(patches, patchI)
1316     {
1317         patchStarts[patchI] = faceI;
1318
1319         faceI += patchSizes[patchI];
1320     }
1321
1322
1323     Pout<< "nFaces:" << dualFaces.size()
1324         << " patchSizes:" << patchSizes
1325         << " patchStarts:" << patchStarts
1326         << endl;
1327
1328
1329     // Add patches. First add zero sized (since mesh still 0 size)
1330     List<polyPatch*> dualPatches(patches.size());
1331
1332     forAll(patches, patchI)
1333     {
1334         const polyPatch& pp = patches[patchI];
1335
1336         dualPatches[patchI] = pp.clone
1337         (
1338             boundaryMesh(),
1339             patchI,
1340             0, //patchSizes[patchI],
1341             0  //patchStarts[patchI]
1342         ).ptr();
1343     }
1344     addPatches(dualPatches);
1345
1346     // Assign to mesh.
1347     resetPrimitives
1348     (
1349         xferMove(dualPoints),
1350         xferMove(dualFaces),
1351         xferMove(dualOwner),
1352         xferMove(dualNeighbour),
1353         patchSizes,
1354         patchStarts
1355     );
1356 }
1357
1358
1359 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1360
1361 // Construct from components
1362 Foam::polyDualMesh::polyDualMesh(const IOobject& io)
1363 :
1364     polyMesh(io),
1365     cellPoint_
1366     (
1367         IOobject
1368         (
1369             "cellPoint",
1370             time().findInstance(meshDir(), "cellPoint"),
1371             meshSubDir,
1372             *this,
1373             IOobject::MUST_READ,
1374             IOobject::NO_WRITE
1375         )
1376     ),
1377     boundaryFacePoint_
1378     (
1379         IOobject
1380         (
1381             "boundaryFacePoint",
1382             time().findInstance(meshDir(), "boundaryFacePoint"),
1383             meshSubDir,
1384             *this,
1385             IOobject::MUST_READ,
1386             IOobject::NO_WRITE
1387         )
1388     )
1389 {}
1390
1391
1392 // Construct from polyMesh
1393 Foam::polyDualMesh::polyDualMesh
1394 (
1395     const polyMesh& mesh,
1396     const labelList& featureEdges,
1397     const labelList& featurePoints
1398 )
1399 :
1400     polyMesh
1401     (
1402         mesh,
1403         xferCopy(pointField()),   // to prevent any warnings "points not allocated"
1404         xferCopy(faceList()),     // to prevent any warnings "faces  not allocated"
1405         xferCopy(cellList())
1406     ),
1407     cellPoint_
1408     (
1409         IOobject
1410         (
1411             "cellPoint",
1412             time().findInstance(meshDir(), "faces"),
1413             meshSubDir,
1414             *this,
1415             IOobject::NO_READ,
1416             IOobject::AUTO_WRITE
1417         ),
1418         labelList(mesh.nCells(), -1)
1419     ),
1420     boundaryFacePoint_
1421     (
1422         IOobject
1423         (
1424             "boundaryFacePoint",
1425             time().findInstance(meshDir(), "faces"),
1426             meshSubDir,
1427             *this,
1428             IOobject::NO_READ,
1429             IOobject::AUTO_WRITE
1430         ),
1431         labelList(mesh.nFaces() - mesh.nInternalFaces())
1432     )
1433 {
1434     calcDual(mesh, featureEdges, featurePoints);
1435 }
1436
1437
1438 // Construct from polyMesh and feature angle
1439 Foam::polyDualMesh::polyDualMesh
1440 (
1441     const polyMesh& mesh,
1442     const scalar featureCos
1443 )
1444 :
1445     polyMesh
1446     (
1447         mesh,
1448         xferCopy(pointField()),   // to prevent any warnings "points not allocated"
1449         xferCopy(faceList()),     // to prevent any warnings "faces  not allocated"
1450         xferCopy(cellList())
1451     ),
1452     cellPoint_
1453     (
1454         IOobject
1455         (
1456             "cellPoint",
1457             time().findInstance(meshDir(), "faces"),
1458             meshSubDir,
1459             *this,
1460             IOobject::NO_READ,
1461             IOobject::AUTO_WRITE
1462         ),
1463         labelList(mesh.nCells(), -1)
1464     ),
1465     boundaryFacePoint_
1466     (
1467         IOobject
1468         (
1469             "boundaryFacePoint",
1470             time().findInstance(meshDir(), "faces"),
1471             meshSubDir,
1472             *this,
1473             IOobject::NO_READ,
1474             IOobject::AUTO_WRITE
1475         ),
1476         labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
1477     )
1478 {
1479     labelList featureEdges, featurePoints;
1480
1481     calcFeatures(mesh, featureCos, featureEdges, featurePoints);
1482     calcDual(mesh, featureEdges, featurePoints);
1483 }
1484
1485
1486 void Foam::polyDualMesh::calcFeatures
1487 (
1488     const polyMesh& mesh,
1489     const scalar featureCos,
1490     labelList& featureEdges,
1491     labelList& featurePoints
1492 )
1493 {
1494     // Create big primitivePatch for all outside.
1495     primitivePatch allBoundary
1496     (
1497         SubList<face>
1498         (
1499             mesh.faces(),
1500             mesh.nFaces() - mesh.nInternalFaces(),
1501             mesh.nInternalFaces()
1502         ),
1503         mesh.points()
1504     );
1505
1506     // For ease of use store patch number per face in allBoundary.
1507     labelList allRegion(allBoundary.size());
1508
1509     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1510
1511     forAll(patches, patchI)
1512     {
1513         const polyPatch& pp = patches[patchI];
1514
1515         forAll(pp, i)
1516         {
1517             allRegion[i + pp.start() - mesh.nInternalFaces()] = patchI;
1518         }
1519     }
1520
1521
1522     // Calculate patch/feature edges
1523     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1524
1525     const labelListList& edgeFaces = allBoundary.edgeFaces();
1526     const vectorField& faceNormals = allBoundary.faceNormals();
1527     const labelList& meshPoints = allBoundary.meshPoints();
1528
1529     boolList isFeatureEdge(edgeFaces.size(), false);
1530
1531     forAll(edgeFaces, edgeI)
1532     {
1533         const labelList& eFaces = edgeFaces[edgeI];
1534
1535         if (eFaces.size() != 2)
1536         {
1537             // Non-manifold. Problem?
1538             const edge& e = allBoundary.edges()[edgeI];
1539
1540             WarningIn("polyDualMesh::calcFeatures") << "Edge "
1541                 << meshPoints[e[0]] << ' ' << meshPoints[e[1]]
1542                 << "  coords:" << mesh.points()[meshPoints[e[0]]]
1543                 << mesh.points()[meshPoints[e[1]]]
1544                 << " has more than 2 faces connected to it:"
1545                 << eFaces.size() << endl;
1546
1547             isFeatureEdge[edgeI] = true;
1548         }
1549         else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1550         {
1551             isFeatureEdge[edgeI] = true;
1552         }
1553         else if
1554         (
1555             (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1556           < featureCos
1557         )
1558         {
1559             isFeatureEdge[edgeI] = true;
1560         }
1561     }
1562
1563
1564     // Calculate feature points
1565     // ~~~~~~~~~~~~~~~~~~~~~~~~
1566
1567     const labelListList& pointEdges = allBoundary.pointEdges();
1568
1569     dynamicLabelList allFeaturePoints(pointEdges.size());
1570
1571     forAll(pointEdges, pointI)
1572     {
1573         const labelList& pEdges = pointEdges[pointI];
1574
1575         label nFeatEdges = 0;
1576
1577         forAll(pEdges, i)
1578         {
1579             if (isFeatureEdge[pEdges[i]])
1580             {
1581                 nFeatEdges++;
1582             }
1583         }
1584         if (nFeatEdges > 2)
1585         {
1586             // Store in mesh vertex label
1587             allFeaturePoints.append(allBoundary.meshPoints()[pointI]);
1588         }
1589     }
1590     featurePoints.transfer(allFeaturePoints);
1591
1592     if (debug)
1593     {
1594         OFstream str("featurePoints.obj");
1595
1596         Pout<< "polyDualMesh::calcFeatures : dumping feature points to "
1597             << str.name() << endl;
1598
1599         forAll(featurePoints, i)
1600         {
1601             label pointI = featurePoints[i];
1602             meshTools::writeOBJ(str, mesh.points()[pointI]);
1603         }
1604     }
1605
1606
1607     // Get all feature edges.
1608     labelList meshEdges
1609     (
1610         allBoundary.meshEdges
1611         (
1612             mesh.edges(),
1613             mesh.cellEdges(),
1614             SubList<label>
1615             (
1616                 mesh.faceOwner(),
1617                 allBoundary.size(),
1618                 mesh.nInternalFaces()
1619             )
1620         )
1621     );
1622
1623     dynamicLabelList allFeatureEdges(isFeatureEdge.size());
1624     forAll(isFeatureEdge, edgeI)
1625     {
1626         if (isFeatureEdge[edgeI])
1627         {
1628             // Store in mesh edge label.
1629             allFeatureEdges.append(meshEdges[edgeI]);
1630         }
1631     }
1632     featureEdges.transfer(allFeatureEdges);
1633 }
1634
1635
1636 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1637
1638 Foam::polyDualMesh::~polyDualMesh()
1639 {}
1640
1641
1642 // ************************************************************************* //