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-2008 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.setSize(pFaces.size()+2+1);
 240         // Store dualVertex for feature edge
 241         dualFace.append(pointToDualPoint[meshPointI]);
 242     }
 243     else
 244     {
 245         dualFace.setSize(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.shrink());
 416     dualFace.clear();
 417
 418     featEdgeIndices2.transfer(featEdgeIndices.shrink());
 419     featEdgeIndices.clear();
 420
 421     if (reverseFace)
 422     {
 423         reverse(dualFace2);
 424
 425         // Correct featEdgeIndices for change in dualFace2
 426         forAll(featEdgeIndices2, i)
 427         {
 428             featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
 429         }
 430         // Reverse indices (might not be nessecary but do anyway)
 431         reverse(featEdgeIndices2);
 432     }
 433 }
 434
 435
 436 void Foam::polyDualMesh::splitFace
 437 (
 438     const polyPatch& patch,
 439     const labelList& pointToDualPoint,
 440
 441     const label pointI,
 442     const labelList& dualFace,
 443     const labelList& featEdgeIndices,
 444
 445     DynamicList<face>& dualFaces,
 446     DynamicList<label>& dualOwner,
 447     DynamicList<label>& dualNeighbour,
 448     DynamicList<label>& dualRegion
 449 )
 450 {
 451
 452     // Split face because of feature edges/points
 453     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 454     label meshPointI = patch.meshPoints()[pointI];
 455
 456     if (pointToDualPoint[meshPointI] >= 0)
 457     {
 458         // Feature point. Do face-centre decomposition.
 459
 460         if (featEdgeIndices.size() < 2)
 461         {
 462             // Feature point but no feature edges. Not handled at the moment
 463             dualFaces.append(face(dualFace));
 464             dualOwner.append(meshPointI);
 465             dualNeighbour.append(-1);
 466             dualRegion.append(patch.index());
 467         }
 468         else
 469         {
 470             // Do 'face-centre' decomposition. Start from first feature
 471             // edge create face up until next feature edge.
 472
 473             forAll(featEdgeIndices, i)
 474             {
 475                 label startFp = featEdgeIndices[i];
 476
 477                 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
 478
 479                 // Collect face from startFp to endFp
 480                 label sz = 0;
 481
 482                 if (endFp > startFp)
 483                 {
 484                     sz = endFp - startFp + 2;
 485                 }
 486                 else
 487                 {
 488                     sz = endFp + dualFace.size() - startFp + 2;
 489                 }
 490                 face subFace(sz);
 491
 492                 // feature point becomes face centre.
 493                 subFace[0] = pointToDualPoint[patch.meshPoints()[pointI]];
 494
 495                 // Copy from startFp up to endFp.
 496                 for (label subFp = 1; subFp < subFace.size(); subFp++)
 497                 {
 498                     subFace[subFp] = dualFace[startFp];
 499
 500                     startFp = (startFp+1) % dualFace.size();
 501                 }
 502
 503                 dualFaces.append(face(subFace));
 504                 dualOwner.append(meshPointI);
 505                 dualNeighbour.append(-1);
 506                 dualRegion.append(patch.index());
 507             }
 508         }
 509     }
 510     else
 511     {
 512         // No feature point. Check if feature edges.
 513         if (featEdgeIndices.size() < 2)
 514         {
 515             // Not enough feature edges. No split.
 516             dualFaces.append(face(dualFace));
 517             dualOwner.append(meshPointI);
 518             dualNeighbour.append(-1);
 519             dualRegion.append(patch.index());
 520         }
 521         else
 522         {
 523             // Multiple feature edges but no feature point.
 524             // Split face along feature edges. Gives weird result if
 525             // number of feature edges > 2.
 526
 527             // Storage for new face
 528             DynamicList<label> subFace(dualFace.size());
 529
 530             forAll(featEdgeIndices, featI)
 531             {
 532                 label startFp = featEdgeIndices[featI];
 533                 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
 534
 535                 label fp = startFp;
 536
 537                 while (true)
 538                 {
 539                     label vertI = dualFace[fp];
 540
 541                     subFace.append(vertI);
 542
 543                     if (fp == endFp)
 544                     {
 545                         break;
 546                     }
 547
 548                     fp = dualFace.fcIndex(fp);
 549                 }
 550
 551                 if (subFace.size() > 2)
 552                 {
 553                     // Enough vertices to create a face from.
 554                     subFace.shrink();
 555
 556                     dualFaces.append(face(subFace));
 557                     dualOwner.append(meshPointI);
 558                     dualNeighbour.append(-1);
 559                     dualRegion.append(patch.index());
 560
 561                     subFace.clear();
 562                 }
 563             }
 564             // Check final face.
 565             if (subFace.size() > 2)
 566             {
 567                 // Enough vertices to create a face from.
 568                 subFace.shrink();
 569
 570                 dualFaces.append(face(subFace));
 571                 dualOwner.append(meshPointI);
 572                 dualNeighbour.append(-1);
 573                 dualRegion.append(patch.index());
 574
 575                 subFace.clear();
 576             }
 577         }
 578     }
 579 }
 580
 581
 582 // Create boundary face for every point in patch
 583 void Foam::polyDualMesh::dualPatch
 584 (
 585     const polyPatch& patch,
 586     const label patchToDualOffset,
 587     const labelList& edgeToDualPoint,
 588     const labelList& pointToDualPoint,
 589
 590     const pointField& dualPoints,
 591
 592     DynamicList<face>& dualFaces,
 593     DynamicList<label>& dualOwner,
 594     DynamicList<label>& dualNeighbour,
 595     DynamicList<label>& dualRegion
 596 )
 597 {
 598     const labelList& meshEdges = patch.meshEdges();
 599
 600     // Whether edge has been done.
 601     // 0 : not
 602     // 1 : done e.start()
 603     // 2 : done e.end()
 604     // 3 : done both
 605     labelList doneEdgeSide(meshEdges.size(), 0);
 606
 607     boolList donePoint(patch.nPoints(), false);
 608
 609
 610     // Do points on edge of patch
 611     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 612
 613     forAll(doneEdgeSide, patchEdgeI)
 614     {
 615         const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
 616
 617         if (eFaces.size() == 1)
 618         {
 619             const edge& e = patch.edges()[patchEdgeI];
 620
 621             forAll(e, eI)
 622             {
 623                 label bitMask = 1<<eI;
 624
 625                 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
 626                 {
 627                     // Construct face by walking around e[eI] starting from
 628                     // patchEdgeI
 629
 630                     label pointI = e[eI];
 631
 632                     label edgeI = patchEdgeI;
 633                     labelList dualFace
 634                     (
 635                         collectPatchSideFace
 636                         (
 637                             patch,
 638                             patchToDualOffset,
 639                             edgeToDualPoint,
 640                             pointToDualPoint,
 641
 642                             pointI,
 643                             edgeI
 644                         )
 645                     );
 646
 647                     // Now edgeI is end edge. Mark as visited
 648                     if (patch.edges()[edgeI][0] == pointI)
 649                     {
 650                         doneEdgeSide[edgeI] |= 1;
 651                     }
 652                     else
 653                     {
 654                         doneEdgeSide[edgeI] |= 2;
 655                     }
 656
 657                     dualFaces.append(face(dualFace));
 658                     dualOwner.append(patch.meshPoints()[pointI]);
 659                     dualNeighbour.append(-1);
 660                     dualRegion.append(patch.index());
 661
 662                     doneEdgeSide[patchEdgeI] |= bitMask;
 663                     donePoint[pointI] = true;
 664                 }
 665             }
 666         }
 667     }
 668
 669
 670
 671     // Do patch-internal points
 672     // ~~~~~~~~~~~~~~~~~~~~~~~~
 673
 674     forAll(donePoint, pointI)
 675     {
 676         if (!donePoint[pointI])
 677         {
 678             labelList dualFace, featEdgeIndices;
 679
 680             collectPatchInternalFace
 681             (
 682                 patch,
 683                 patchToDualOffset,
 684                 edgeToDualPoint,
 685                 pointI,
 686                 patch.pointEdges()[pointI][0],  // Arbitrary starting edge
 687
 688                 dualFace,
 689                 featEdgeIndices
 690             );
 691
 692             //- Either keep in one piece or split face according to feature.
 693
 694             //// Keep face in one piece.
 695             //dualFaces.append(face(dualFace));
 696             //dualOwner.append(patch.meshPoints()[pointI]);
 697             //dualNeighbour.append(-1);
 698             //dualRegion.append(patch.index());
 699
 700             splitFace
 701             (
 702                 patch,
 703                 pointToDualPoint,
 704                 pointI,
 705                 dualFace,
 706                 featEdgeIndices,
 707
 708                 dualFaces,
 709                 dualOwner,
 710                 dualNeighbour,
 711                 dualRegion
 712             );
 713
 714             donePoint[pointI] = true;
 715         }
 716     }
 717 }
 718
 719
 720 void Foam::polyDualMesh::calcDual
 721 (
 722     const polyMesh& mesh,
 723     const labelList& featureEdges,
 724     const labelList& featurePoints
 725 )
 726 {
 727     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 728     const label nIntFaces = mesh.nInternalFaces();
 729
 730
 731     // Get patch for all of outside
 732     primitivePatch allBoundary
 733     (
 734         SubList<face>
 735         (
 736             mesh.faces(),
 737             mesh.nFaces() - nIntFaces,
 738             nIntFaces
 739         ),
 740         mesh.points()
 741     );
 742
 743     // Correspondence from patch edge to mesh edge.
 744     labelList meshEdges
 745     (
 746         allBoundary.meshEdges
 747         (
 748             mesh.edges(),
 749             mesh.cellEdges(),
 750             SubList<label>(mesh.faceOwner(), allBoundary.size(), nIntFaces)
 751         )
 752     );
 753
 754     {
 755         pointSet nonManifoldPoints
 756         (
 757             mesh,
 758             "nonManifoldPoints",
 759             meshEdges.size() / 100
 760         );
 761
 762         allBoundary.checkPointManifold(true, &nonManifoldPoints);
 763
 764         if (nonManifoldPoints.size() > 0)
 765         {
 766             nonManifoldPoints.write();
 767
 768             FatalErrorIn
 769             (
 770                 "polyDualMesh::calcDual(const polyMesh&, const labelList&"
 771                 ", const labelList&)"
 772             )   << "There are " << nonManifoldPoints.size() << " points where"
 773                 << " the outside of the mesh is non-manifold." << nl
 774                 << "Such a mesh cannot be converted to a dual." << nl
 775                 << "Writing points at non-manifold sites to pointSet "
 776                 << nonManifoldPoints.name()
 777                 << exit(FatalError);
 778         }
 779     }
 780
 781
 782     // Assign points
 783     // ~~~~~~~~~~~~~
 784
 785     // mesh label   dualMesh vertex
 786     // ----------   ---------------
 787     // cellI        cellI
 788     // faceI        nCells+faceI-nIntFaces
 789     // featEdgeI    nCells+nFaces-nIntFaces+featEdgeI
 790     // featPointI   nCells+nFaces-nIntFaces+nFeatEdges+featPointI
 791
 792     pointField dualPoints
 793     (
 794         mesh.nCells()                           // cell centres
 795       + mesh.nFaces() - nIntFaces               // boundary face centres
 796       + featureEdges.size()                     // additional boundary edges
 797       + featurePoints.size()                    // additional boundary points
 798     );
 799
 800     label dualPointI = 0;
 801
 802
 803     // Cell centres.
 804     const pointField& cellCentres = mesh.cellCentres();
 805
 806     cellPoint_.setSize(cellCentres.size());
 807
 808     forAll(cellCentres, cellI)
 809     {
 810         cellPoint_[cellI] = dualPointI;
 811         dualPoints[dualPointI++] = cellCentres[cellI];
 812     }
 813
 814     // Boundary faces centres
 815     const pointField& faceCentres = mesh.faceCentres();
 816
 817     boundaryFacePoint_.setSize(mesh.nFaces() - nIntFaces);
 818
 819     for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
 820     {
 821         boundaryFacePoint_[faceI - nIntFaces] = dualPointI;
 822         dualPoints[dualPointI++] = faceCentres[faceI];
 823     }
 824
 825     // Edge status:
 826     //  >0 : dual point label (edge is feature edge)
 827     //  -1 : is boundary edge.
 828     //  -2 : is internal edge.
 829     labelList edgeToDualPoint(mesh.nEdges(), -2);
 830
 831     forAll(meshEdges, patchEdgeI)
 832     {
 833         label edgeI = meshEdges[patchEdgeI];
 834
 835         edgeToDualPoint[edgeI] = -1;
 836     }
 837
 838     forAll(featureEdges, i)
 839     {
 840         label edgeI = featureEdges[i];
 841
 842         const edge& e = mesh.edges()[edgeI];
 843
 844         edgeToDualPoint[edgeI] = dualPointI;
 845         dualPoints[dualPointI++] = e.centre(mesh.points());
 846     }
 847
 848
 849
 850     // Point status:
 851     //  >0 : dual point label (point is feature point)
 852     //  -1 : is point on edge of patch
 853     //  -2 : is point on patch (but not on edge)
 854     //  -3 : is internal point.
 855     labelList pointToDualPoint(mesh.nPoints(), -3);
 856
 857     forAll(patches, patchI)
 858     {
 859         const labelList& meshPoints = patches[patchI].meshPoints();
 860
 861         forAll(meshPoints, i)
 862         {
 863             pointToDualPoint[meshPoints[i]] = -2;
 864         }
 865
 866         const labelListList& loops = patches[patchI].edgeLoops();
 867
 868         forAll(loops, i)
 869         {
 870             const labelList& loop = loops[i];
 871
 872             forAll(loop, j)
 873             {
 874                  pointToDualPoint[meshPoints[loop[j]]] = -1;
 875             }
 876         }
 877     }
 878
 879     forAll(featurePoints, i)
 880     {
 881         label pointI = featurePoints[i];
 882
 883         pointToDualPoint[pointI] = dualPointI;
 884         dualPoints[dualPointI++] = mesh.points()[pointI];
 885     }
 886
 887
 888     // Storage for new faces.
 889     // Dynamic sized since we don't know size.
 890
 891     DynamicList<face> dynDualFaces(mesh.nEdges());
 892     DynamicList<label> dynDualOwner(mesh.nEdges());
 893     DynamicList<label> dynDualNeighbour(mesh.nEdges());
 894     DynamicList<label> dynDualRegion(mesh.nEdges());
 895
 896
 897     // Generate faces from edges on the boundary
 898     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 899
 900     forAll(meshEdges, patchEdgeI)
 901     {
 902         label edgeI = meshEdges[patchEdgeI];
 903
 904         const edge& e = mesh.edges()[edgeI];
 905
 906         label owner = -1;
 907         label neighbour = -1;
 908
 909         if (e[0] < e[1])
 910         {
 911             owner = e[0];
 912             neighbour = e[1];
 913         }
 914         else
 915         {
 916             owner = e[1];
 917             neighbour = e[0];
 918         }
 919
 920         // Find the boundary faces using the edge.
 921         const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
 922
 923         if (patchFaces.size() != 2)
 924         {
 925             FatalErrorIn("polyDualMesh::calcDual")
 926                 << "Cannot handle edges with " << patchFaces.size()
 927                 << " connected boundary faces."
 928                 << abort(FatalError);
 929         }
 930
 931         label face0 = patchFaces[0] + nIntFaces;
 932         const face& f0 = mesh.faces()[face0];
 933
 934         label face1 = patchFaces[1] + nIntFaces;
 935
 936         // We want to start walking from patchFaces[0] or patchFaces[1],
 937         // depending on which one uses owner,neighbour in the right order.
 938
 939         label startFaceI = -1;
 940         label endFaceI = -1;
 941
 942         label index = findIndex(f0, neighbour);
 943
 944         if (f0.nextLabel(index) == owner)
 945         {
 946             startFaceI = face0;
 947             endFaceI = face1;
 948         }
 949         else
 950         {
 951             startFaceI = face1;
 952             endFaceI = face0;
 953         }
 954
 955         // Now walk from startFaceI to cell to face to cell etc. until endFaceI
 956
 957         DynamicList<label> dualFace;
 958
 959         if (edgeToDualPoint[edgeI] >= 0)
 960         {
 961             // Number of cells + 2 boundary faces + feature edge point
 962             dualFace.setSize(mesh.edgeCells()[edgeI].size()+2+1);
 963             // Store dualVertex for feature edge
 964             dualFace.append(edgeToDualPoint[edgeI]);
 965         }
 966         else
 967         {
 968             dualFace.setSize(mesh.edgeCells()[edgeI].size()+2);
 969         }
 970
 971         // Store dual vertex for starting face.
 972         dualFace.append(mesh.nCells() + startFaceI - nIntFaces);
 973
 974         label cellI = mesh.faceOwner()[startFaceI];
 975         label faceI = startFaceI;
 976
 977         while (true)
 978         {
 979             // Store dual vertex from cellI.
 980             dualFace.append(cellI);
 981
 982             // Cross cell to other face on edge.
 983             label f0, f1;
 984             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
 985
 986             if (f0 == faceI)
 987             {
 988                 faceI = f1;
 989             }
 990             else
 991             {
 992                 faceI = f0;
 993             }
 994
 995             // Cross face to other cell.
 996             if (faceI == endFaceI)
 997             {
 998                 break;
 999             }
1000
1001             if (mesh.faceOwner()[faceI] == cellI)
1002             {
1003                 cellI = mesh.faceNeighbour()[faceI];
1004             }
1005             else
1006             {
1007                 cellI = mesh.faceOwner()[faceI];
1008             }
1009         }
1010
1011         // Store dual vertex for endFace.
1012         dualFace.append(mesh.nCells() + endFaceI - nIntFaces);
1013
1014         dynDualFaces.append(face(dualFace.shrink()));
1015         dynDualOwner.append(owner);
1016         dynDualNeighbour.append(neighbour);
1017         dynDualRegion.append(-1);
1018
1019         {
1020             // Check orientation.
1021             const face& f = dynDualFaces[dynDualFaces.size()-1];
1022             vector n = f.normal(dualPoints);
1023             if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1024             {
1025                 WarningIn("calcDual") << "Incorrect orientation"
1026                     << " on boundary edge:" << edgeI
1027                     << mesh.points()[mesh.edges()[edgeI][0]]
1028                     << mesh.points()[mesh.edges()[edgeI][1]]
1029                     << endl;
1030             }
1031         }
1032     }
1033
1034
1035     // Generate faces from internal edges
1036     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1037
1038     forAll(edgeToDualPoint, edgeI)
1039     {
1040         if (edgeToDualPoint[edgeI] == -2)
1041         {
1042             // Internal edge.
1043
1044             // Find dual owner, neighbour
1045
1046             const edge& e = mesh.edges()[edgeI];
1047
1048             label owner = -1;
1049             label neighbour = -1;
1050
1051             if (e[0] < e[1])
1052             {
1053                 owner = e[0];
1054                 neighbour = e[1];
1055             }
1056             else
1057             {
1058                 owner = e[1];
1059                 neighbour = e[0];
1060             }
1061
1062             // Get a starting cell
1063             const labelList& eCells = mesh.edgeCells()[edgeI];
1064
1065             label cellI = eCells[0];
1066
1067             // Get the two faces on the cell and edge.
1068             label face0, face1;
1069             meshTools::getEdgeFaces(mesh, cellI, edgeI, face0, face1);
1070
1071             // Find the starting face by looking at the order in which
1072             // the face uses the owner, neighbour
1073             const face& f0 = mesh.faces()[face0];
1074
1075             label index = findIndex(f0, neighbour);
1076
1077             bool f0OrderOk = (f0.nextLabel(index) == owner);
1078
1079             label startFaceI = -1;
1080
1081             if (f0OrderOk == (mesh.faceOwner()[face0] == cellI))
1082             {
1083                 startFaceI = face0;
1084             }
1085             else
1086             {
1087                 startFaceI = face1;
1088             }
1089
1090             // Walk face-cell-face until starting face reached.
1091             DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
1092
1093             label faceI = startFaceI;
1094
1095             while (true)
1096             {
1097                 // Store dual vertex from cellI.
1098                 dualFace.append(cellI);
1099
1100                 // Cross cell to other face on edge.
1101                 label f0, f1;
1102                 meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
1103
1104                 if (f0 == faceI)
1105                 {
1106                     faceI = f1;
1107                 }
1108                 else
1109                 {
1110                     faceI = f0;
1111                 }
1112
1113                 // Cross face to other cell.
1114                 if (faceI == startFaceI)
1115                 {
1116                     break;
1117                 }
1118
1119                 if (mesh.faceOwner()[faceI] == cellI)
1120                 {
1121                     cellI = mesh.faceNeighbour()[faceI];
1122                 }
1123                 else
1124                 {
1125                     cellI = mesh.faceOwner()[faceI];
1126                 }
1127             }
1128
1129             dynDualFaces.append(face(dualFace.shrink()));
1130             dynDualOwner.append(owner);
1131             dynDualNeighbour.append(neighbour);
1132             dynDualRegion.append(-1);
1133
1134             {
1135                 // Check orientation.
1136                 const face& f = dynDualFaces[dynDualFaces.size()-1];
1137                 vector n = f.normal(dualPoints);
1138                 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1139                 {
1140                     WarningIn("calcDual") << "Incorrect orientation"
1141                         << " on internal edge:" << edgeI
1142                         << mesh.points()[mesh.edges()[edgeI][0]]
1143                         << mesh.points()[mesh.edges()[edgeI][1]]
1144                         << endl;
1145                 }
1146             }
1147         }
1148     }
1149
1150     // Dump faces.
1151     if (debug)
1152     {
1153         dynDualFaces.shrink();
1154         dynDualOwner.shrink();
1155         dynDualNeighbour.shrink();
1156         dynDualRegion.shrink();
1157
1158         OFstream str("dualInternalFaces.obj");
1159         Pout<< "polyDualMesh::calcDual : dumping internal faces to "
1160             << str.name() << endl;
1161
1162         forAll(dualPoints, dualPointI)
1163         {
1164             meshTools::writeOBJ(str, dualPoints[dualPointI]);
1165         }
1166         forAll(dynDualFaces, dualFaceI)
1167         {
1168             const face& f = dynDualFaces[dualFaceI];
1169
1170             str<< 'f';
1171             forAll(f, fp)
1172             {
1173                 str<< ' ' << f[fp]+1;
1174             }
1175             str<< nl;
1176         }
1177     }
1178
1179     const label nInternalFaces = dynDualFaces.size();
1180
1181     // Outside faces
1182     // ~~~~~~~~~~~~~
1183
1184     forAll(patches, patchI)
1185     {
1186         const polyPatch& pp = patches[patchI];
1187
1188         dualPatch
1189         (
1190             pp,
1191             mesh.nCells() + pp.start() - nIntFaces,
1192             edgeToDualPoint,
1193             pointToDualPoint,
1194
1195             dualPoints,
1196
1197             dynDualFaces,
1198             dynDualOwner,
1199             dynDualNeighbour,
1200             dynDualRegion
1201         );
1202     }
1203
1204
1205     // Transfer face info to straight lists
1206     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1207     faceList dualFaces(dynDualFaces.shrink(), true);
1208     dynDualFaces.clear();
1209
1210     labelList dualOwner(dynDualOwner.shrink(), true);
1211     dynDualOwner.clear();
1212
1213     labelList dualNeighbour(dynDualNeighbour.shrink(), true);
1214     dynDualNeighbour.clear();
1215
1216     labelList dualRegion(dynDualRegion.shrink(), true);
1217     dynDualRegion.clear();
1218
1219
1220
1221     // Dump faces.
1222     if (debug)
1223     {
1224         OFstream str("dualFaces.obj");
1225         Pout<< "polyDualMesh::calcDual : dumping all faces to "
1226             << str.name() << endl;
1227
1228         forAll(dualPoints, dualPointI)
1229         {
1230             meshTools::writeOBJ(str, dualPoints[dualPointI]);
1231         }
1232         forAll(dualFaces, dualFaceI)
1233         {
1234             const face& f = dualFaces[dualFaceI];
1235
1236             str<< 'f';
1237             forAll(f, fp)
1238             {
1239                 str<< ' ' << f[fp]+1;
1240             }
1241             str<< nl;
1242         }
1243     }
1244
1245     // Create cells.
1246     cellList dualCells(mesh.nPoints());
1247
1248     forAll(dualCells, cellI)
1249     {
1250         dualCells[cellI].setSize(0);
1251     }
1252
1253     forAll(dualOwner, faceI)
1254     {
1255         label cellI = dualOwner[faceI];
1256
1257         labelList& cFaces = dualCells[cellI];
1258
1259         label sz = cFaces.size();
1260         cFaces.setSize(sz+1);
1261         cFaces[sz] = faceI;
1262     }
1263     forAll(dualNeighbour, faceI)
1264     {
1265         label cellI = dualNeighbour[faceI];
1266
1267         if (cellI != -1)
1268         {
1269             labelList& cFaces = dualCells[cellI];
1270
1271             label sz = cFaces.size();
1272             cFaces.setSize(sz+1);
1273             cFaces[sz] = faceI;
1274         }
1275     }
1276
1277
1278     // Do upper-triangular face ordering. Determines face reordering map and
1279     // number of internal faces.
1280     label dummy;
1281
1282     labelList oldToNew
1283     (
1284         getFaceOrder
1285         (
1286             dualOwner,
1287             dualNeighbour,
1288             dualCells,
1289             dummy
1290         )
1291     );
1292
1293     // Reorder faces.
1294     inplaceReorder(oldToNew, dualFaces);
1295     inplaceReorder(oldToNew, dualOwner);
1296     inplaceReorder(oldToNew, dualNeighbour);
1297     inplaceReorder(oldToNew, dualRegion);
1298     forAll(dualCells, cellI)
1299     {
1300         inplaceRenumber(oldToNew, dualCells[cellI]);
1301     }
1302
1303
1304     // Create patches
1305     labelList patchSizes(patches.size(), 0);
1306
1307     forAll(dualRegion, faceI)
1308     {
1309         if (dualRegion[faceI] >= 0)
1310         {
1311             patchSizes[dualRegion[faceI]]++;
1312         }
1313     }
1314
1315     labelList patchStarts(patches.size(), 0);
1316
1317     label faceI = nInternalFaces;
1318
1319     forAll(patches, patchI)
1320     {
1321         patchStarts[patchI] = faceI;
1322
1323         faceI += patchSizes[patchI];
1324     }
1325
1326
1327     Pout<< "nFaces:" << dualFaces.size()
1328         << " patchSizes:" << patchSizes
1329         << " patchStarts:" << patchStarts
1330         << endl;
1331
1332
1333     // Add patches. First add zero sized (since mesh still 0 size)
1334     List<polyPatch*> dualPatches(patches.size());
1335
1336     forAll(patches, patchI)
1337     {
1338         const polyPatch& pp = patches[patchI];
1339
1340         dualPatches[patchI] = pp.clone
1341         (
1342             boundaryMesh(),
1343             patchI,
1344             0, //patchSizes[patchI],
1345             0  //patchStarts[patchI]
1346         ).ptr();
1347     }
1348     addPatches(dualPatches);
1349
1350     // Assign to mesh.
1351     resetPrimitives
1352     (
1353         dualFaces.size(),
1354         dualPoints,
1355         dualFaces,
1356         dualOwner,
1357         dualNeighbour,
1358         patchSizes,
1359         patchStarts
1360     );
1361 }
1362
1363
1364 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1365
1366 // Construct from components
1367 Foam::polyDualMesh::polyDualMesh(const IOobject& io)
1368 :
1369     polyMesh(io),
1370     cellPoint_
1371     (
1372         IOobject
1373         (
1374             "cellPoint",
1375             time().findInstance(meshDir(), "cellPoint"),
1376             meshSubDir,
1377             *this,
1378             IOobject::MUST_READ,
1379             IOobject::NO_WRITE
1380         )
1381     ),
1382     boundaryFacePoint_
1383     (
1384         IOobject
1385         (
1386             "boundaryFacePoint",
1387             time().findInstance(meshDir(), "boundaryFacePoint"),
1388             meshSubDir,
1389             *this,
1390             IOobject::MUST_READ,
1391             IOobject::NO_WRITE
1392         )
1393     )
1394 {}
1395
1396
1397 // Construct from polyMesh
1398 Foam::polyDualMesh::polyDualMesh
1399 (
1400     const polyMesh& mesh,
1401     const labelList& featureEdges,
1402     const labelList& featurePoints
1403 )
1404 :
1405     polyMesh
1406     (
1407         mesh,
1408         pointField(0),
1409         faceList(0),
1410         cellList(0)
1411     ),
1412     cellPoint_
1413     (
1414         IOobject
1415         (
1416             "cellPoint",
1417             time().findInstance(meshDir(), "faces"),
1418             meshSubDir,
1419             *this,
1420             IOobject::NO_READ,
1421             IOobject::AUTO_WRITE
1422         ),
1423         labelList(mesh.nCells(), -1)
1424     ),
1425     boundaryFacePoint_
1426     (
1427         IOobject
1428         (
1429             "boundaryFacePoint",
1430             time().findInstance(meshDir(), "faces"),
1431             meshSubDir,
1432             *this,
1433             IOobject::NO_READ,
1434             IOobject::AUTO_WRITE
1435         ),
1436         labelList(mesh.nFaces() - mesh.nInternalFaces())
1437     )
1438 {
1439     calcDual(mesh, featureEdges, featurePoints);
1440 }
1441
1442
1443 // Construct from polyMesh and feature angle
1444 Foam::polyDualMesh::polyDualMesh
1445 (
1446     const polyMesh& mesh,
1447     const scalar featureCos
1448 )
1449 :
1450     polyMesh
1451     (
1452         mesh,
1453         pointField(0),  // to prevent any warnings "points not allocated"
1454         faceList(0),    //            ,,            faces ,,
1455         cellList(0)
1456     ),
1457     cellPoint_
1458     (
1459         IOobject
1460         (
1461             "cellPoint",
1462             time().findInstance(meshDir(), "faces"),
1463             meshSubDir,
1464             *this,
1465             IOobject::NO_READ,
1466             IOobject::AUTO_WRITE
1467         ),
1468         labelList(mesh.nCells(), -1)
1469     ),
1470     boundaryFacePoint_
1471     (
1472         IOobject
1473         (
1474             "boundaryFacePoint",
1475             time().findInstance(meshDir(), "faces"),
1476             meshSubDir,
1477             *this,
1478             IOobject::NO_READ,
1479             IOobject::AUTO_WRITE
1480         ),
1481         labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
1482     )
1483 {
1484     labelList featureEdges, featurePoints;
1485
1486     calcFeatures(mesh, featureCos, featureEdges, featurePoints);
1487
1488     calcDual(mesh, featureEdges, featurePoints);
1489 }
1490
1491
1492 void Foam::polyDualMesh::calcFeatures
1493 (
1494     const polyMesh& mesh,
1495     const scalar featureCos,
1496     labelList& featureEdges,
1497     labelList& featurePoints
1498 )
1499 {
1500     // Create big primitivePatch for all outside.
1501     primitivePatch allBoundary
1502     (
1503         SubList<face>
1504         (
1505             mesh.faces(),
1506             mesh.nFaces() - mesh.nInternalFaces(),
1507             mesh.nInternalFaces()
1508         ),
1509         mesh.points()
1510     );
1511
1512     // For ease of use store patch number per face in allBoundary.
1513     labelList allRegion(allBoundary.size());
1514
1515     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1516
1517     forAll(patches, patchI)
1518     {
1519         const polyPatch& pp = patches[patchI];
1520
1521         forAll(pp, i)
1522         {
1523             allRegion[i + pp.start() - mesh.nInternalFaces()] = patchI;
1524         }
1525     }
1526
1527
1528     // Calculate patch/feature edges
1529     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1530
1531     const labelListList& edgeFaces = allBoundary.edgeFaces();
1532     const vectorField& faceNormals = allBoundary.faceNormals();
1533     const labelList& meshPoints = allBoundary.meshPoints();
1534
1535     boolList isFeatureEdge(edgeFaces.size(), false);
1536
1537     forAll(edgeFaces, edgeI)
1538     {
1539         const labelList& eFaces = edgeFaces[edgeI];
1540
1541         if (eFaces.size() != 2)
1542         {
1543             // Non-manifold. Problem?
1544             const edge& e = allBoundary.edges()[edgeI];
1545
1546             WarningIn("polyDualMesh::calcFeatures") << "Edge "
1547                 << meshPoints[e[0]] << ' ' << meshPoints[e[1]]
1548                 << "  coords:" << mesh.points()[meshPoints[e[0]]]
1549                 << mesh.points()[meshPoints[e[1]]]
1550                 << " has more than 2 faces connected to it:"
1551                 << eFaces.size() << endl;
1552
1553             isFeatureEdge[edgeI] = true;
1554         }
1555         else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1556         {
1557             isFeatureEdge[edgeI] = true;
1558         }
1559         else if
1560         (
1561             (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1562           < featureCos
1563         )
1564         {
1565             isFeatureEdge[edgeI] = true;
1566         }
1567     }
1568
1569
1570     // Calculate feature points
1571     // ~~~~~~~~~~~~~~~~~~~~~~~~
1572
1573     const labelListList& pointEdges = allBoundary.pointEdges();
1574
1575     DynamicList<label> allFeaturePoints(pointEdges.size());
1576
1577     forAll(pointEdges, pointI)
1578     {
1579         const labelList& pEdges = pointEdges[pointI];
1580
1581         label nFeatEdges = 0;
1582
1583         forAll(pEdges, i)
1584         {
1585             if (isFeatureEdge[pEdges[i]])
1586             {
1587                 nFeatEdges++;
1588             }
1589         }
1590         if (nFeatEdges > 2)
1591         {
1592             // Store in mesh vertex label
1593             allFeaturePoints.append(allBoundary.meshPoints()[pointI]);
1594         }
1595     }
1596     featurePoints.transfer(allFeaturePoints.shrink());
1597     allFeaturePoints.clear();
1598
1599     if (debug)
1600     {
1601         OFstream str("featurePoints.obj");
1602
1603         Pout<< "polyDualMesh::calcFeatures : dumping feature points to "
1604             << str.name() << endl;
1605
1606         forAll(featurePoints, i)
1607         {
1608             label pointI = featurePoints[i];
1609             meshTools::writeOBJ(str, mesh.points()[pointI]);
1610         }
1611     }
1612
1613
1614     // Get all feature edges.
1615     labelList meshEdges
1616     (
1617         allBoundary.meshEdges
1618         (
1619             mesh.edges(),
1620             mesh.cellEdges(),
1621             SubList<label>
1622             (
1623                 mesh.faceOwner(),
1624                 allBoundary.size(),
1625                 mesh.nInternalFaces()
1626             )
1627         )
1628     );
1629
1630     DynamicList<label> allFeatureEdges(isFeatureEdge.size());
1631     forAll(isFeatureEdge, edgeI)
1632     {
1633         if (isFeatureEdge[edgeI])
1634         {
1635             // Store in mesh edge label.
1636             allFeatureEdges.append(meshEdges[edgeI]);
1637         }
1638     }
1639     featureEdges.transfer(allFeatureEdges.shrink());
1640     allFeatureEdges.clear();
1641 }
1642
1643
1644 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1645
1646 Foam::polyDualMesh::~polyDualMesh()
1647 {}
1648
1649
1650 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1651
1652
1653 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1654
1655
1656 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
1657
1658
1659 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1660
1661
1662 // ************************************************************************* //