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