src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 \*---------------------------------------------------------------------------*/
  26
  27 #include "dynamicTopoFvMesh.H"
  28
  29 #include "triFace.H"
  30 #include "objectMap.H"
  31 #include "changeMap.H"
  32 #include "triPointRef.H"
  33 #include "linePointRef.H"
  34 #include "multiThreader.H"
  35 #include "coupledInfo.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  43
  44 // Perform a Delaunay test on an internal face
  45 // - Returns 'true' if the test failed
  46 bool dynamicTopoFvMesh::testDelaunay
  47 (
  48     const label fIndex
  49 ) const
  50 {
  51     bool failed = false, procCoupled = false;
  52     label eIndex = -1, pointIndex = -1, fLabel = -1;
  53     label sIndex = -1, pIndex = -1;
  54     FixedList<triFace, 2> triFaces(triFace(-1, -1, -1));
  55     FixedList<bool, 2> foundTriFace(false);
  56
  57     // If this entity was deleted, skip it.
  58     if (faces_[fIndex].empty())
  59     {
  60         return failed;
  61     }
  62
  63     // If face is not shared by prism cells, skip it.
  64     label own = owner_[fIndex], nei = neighbour_[fIndex];
  65
  66     if (cells_[own].size() != 5)
  67     {
  68         return failed;
  69     }
  70
  71     if (nei > -1)
  72     {
  73         if (cells_[nei].size() != 5)
  74         {
  75             return failed;
  76         }
  77     }
  78
  79     // Boundary faces are discarded.
  80     if (whichPatch(fIndex) > -1)
  81     {
  82         procCoupled = processorCoupledEntity(fIndex);
  83
  84         if (!procCoupled)
  85         {
  86             return failed;
  87         }
  88     }
  89
  90     const labelList& fEdges = faceEdges_[fIndex];
  91
  92     forAll(fEdges, edgeI)
  93     {
  94         // Break out if both triangular faces are found
  95         if (foundTriFace[0] && foundTriFace[1])
  96         {
  97             break;
  98         }
  99
 100         // Obtain edgeFaces for this edge
 101         const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
 102
 103         forAll(eFaces, faceI)
 104         {
 105             const face& thisFace = faces_[eFaces[faceI]];
 106
 107             if (thisFace.size() == 3)
 108             {
 109                 if (foundTriFace[0])
 110                 {
 111                     // Update the second face.
 112                     triFaces[1][0] = thisFace[0];
 113                     triFaces[1][1] = thisFace[1];
 114                     triFaces[1][2] = thisFace[2];
 115
 116                     foundTriFace[1] = true;
 117                 }
 118                 else
 119                 {
 120                     // Update the first face.
 121                     triFaces[0][0] = thisFace[0];
 122                     triFaces[0][1] = thisFace[1];
 123                     triFaces[0][2] = thisFace[2];
 124
 125                     foundTriFace[0] = true;
 126                     fLabel = eFaces[faceI];
 127
 128                     // Take this edge
 129                     eIndex = fEdges[edgeI];
 130
 131                     if (procCoupled)
 132                     {
 133                         // Stop searching for processor boundary cases
 134                         foundTriFace[1] = true;
 135                     }
 136                 }
 137             }
 138         }
 139     }
 140
 141     const edge& checkEdge = edges_[eIndex];
 142
 143     // Configure the comparison edge
 144     edge cEdge(-1, -1);
 145
 146     if (procCoupled)
 147     {
 148         // If this is a new entity, bail out for now.
 149         // This will be handled at the next time-step.
 150         if (fIndex >= nOldFaces_)
 151         {
 152             return failed;
 153         }
 154
 155         const label faceEnum = coupleMap::FACE;
 156         const label pointEnum = coupleMap::POINT;
 157
 158         // Check slaves
 159         bool foundEdge = false;
 160
 161         forAll(procIndices_, pI)
 162         {
 163             // Fetch non-const reference to subMeshes
 164             const coupledInfo& recvMesh = recvMeshes_[pI];
 165
 166             const coupleMap& cMap = recvMesh.map();
 167             const dynamicTopoFvMesh& sMesh = recvMesh.subMesh();
 168
 169             if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
 170             {
 171                 // Find equivalent points on the slave
 172                 cEdge[0] = cMap.findSlave(pointEnum, checkEdge[0]);
 173                 cEdge[1] = cMap.findSlave(pointEnum, checkEdge[1]);
 174
 175                 // Find a triangular face containing cEdge
 176                 const labelList& sfE = sMesh.faceEdges_[sIndex];
 177
 178                 forAll(sfE, edgeI)
 179                 {
 180                     // Obtain edgeFaces for this edge
 181                     const labelList& seF = sMesh.edgeFaces_[sfE[edgeI]];
 182
 183                     forAll(seF, faceI)
 184                     {
 185                         const face& tF = sMesh.faces_[seF[faceI]];
 186
 187                         if (tF.size() == 3)
 188                         {
 189                             if
 190                             (
 191                                 (findIndex(tF, cEdge[0]) > -1) &&
 192                                 (findIndex(tF, cEdge[1]) > -1)
 193                             )
 194                             {
 195                                 triFaces[1][0] = tF[0];
 196                                 triFaces[1][1] = tF[1];
 197                                 triFaces[1][2] = tF[2];
 198
 199                                 foundEdge = true;
 200
 201                                 break;
 202                             }
 203                         }
 204                     }
 205
 206                     if (foundEdge)
 207                     {
 208                         break;
 209                     }
 210                 }
 211
 212                 // Store patch index for later
 213                 pIndex = pI;
 214
 215                 break;
 216             }
 217         }
 218
 219         if (sIndex == -1 || !foundEdge)
 220         {
 221             // Write out for post-processing
 222             writeVTK("coupledDelaunayFace_" + Foam::name(fIndex), fIndex, 2);
 223
 224             FatalErrorIn
 225             (
 226                 "bool dynamicTopoFvMesh::testDelaunay"
 227                 "(const label fIndex) const"
 228             )
 229                 << "Coupled maps were improperly specified." << nl
 230                 << " Slave index not found for: " << nl
 231                 << " Face: " << fIndex << nl
 232                 << " checkEdge: " << checkEdge << nl
 233                 << " cEdge: " << cEdge << nl
 234                 << abort(FatalError);
 235         }
 236     }
 237     else
 238     {
 239         // Non-coupled case
 240         cEdge = checkEdge;
 241     }
 242
 243     // Obtain point references for the first face
 244     point a = points_[triFaces[0][0]];
 245
 246     const face& faceToCheck = faces_[fLabel];
 247
 248     vector cCentre =
 249     (
 250         triPointRef
 251         (
 252             points_[faceToCheck[0]],
 253             points_[faceToCheck[1]],
 254             points_[faceToCheck[2]]
 255         ).circumCentre()
 256     );
 257
 258     scalar rSquared = (a - cCentre) & (a - cCentre);
 259
 260     // Find the isolated point on the second face
 261     point otherPoint = vector::zero;
 262
 263     // Check the first point
 264     if (triFaces[1][0] != cEdge[0] && triFaces[1][0] != cEdge[1])
 265     {
 266         pointIndex = triFaces[1][0];
 267     }
 268
 269     // Check the second point
 270     if (triFaces[1][1] != cEdge[0] && triFaces[1][1] != cEdge[1])
 271     {
 272         pointIndex = triFaces[1][1];
 273     }
 274
 275     // Check the third point
 276     if (triFaces[1][2] != cEdge[0] && triFaces[1][2] != cEdge[1])
 277     {
 278         pointIndex = triFaces[1][2];
 279     }
 280
 281     if (procCoupled)
 282     {
 283         otherPoint = recvMeshes_[pIndex].subMesh().points_[pointIndex];
 284     }
 285     else
 286     {
 287         otherPoint = points_[pointIndex];
 288     }
 289
 290     // ...and determine whether it lies in this circle
 291     if (((otherPoint - cCentre) & (otherPoint - cCentre)) < rSquared)
 292     {
 293         // Failed the test.
 294         failed = true;
 295     }
 296
 297     return failed;
 298 }
 299
 300
 301 // Method for the swapping of a quad-face in 2D
 302 // - Returns a changeMap with a type specifying:
 303 //     1: Swap sequence was successful
 304 //    -1: Swap sequence failed
 305 const changeMap dynamicTopoFvMesh::swapQuadFace
 306 (
 307     const label fIndex
 308 )
 309 {
 310     changeMap map;
 311
 312     face f;
 313     bool found = false;
 314     label commonIndex = -1;
 315     FixedList<label,4> otherPointIndex(-1), nextToOtherPoint(-1);
 316     FixedList<label,2> c0BdyIndex(-1), c0IntIndex(-1);
 317     FixedList<label,2> c1BdyIndex(-1), c1IntIndex(-1);
 318     FixedList<face,2> c0BdyFace(face(3)), c0IntFace(face(4));
 319     FixedList<face,2> c1BdyFace(face(3)), c1IntFace(face(4));
 320     FixedList<label,2> commonEdgeIndex(-1);
 321     FixedList<edge,2> commonEdges(edge(-1, -1));
 322     FixedList<label,4> otherEdgeIndex(-1);
 323     FixedList<label,4> commonFaceIndex(-1), cornerEdgeIndex(-1);
 324     FixedList<face,4> commonFaces(face(3)), commonIntFaces(face(4));
 325     FixedList<label,4> commonIntFaceIndex(-1);
 326     FixedList<bool,2> foundTriFace0(false), foundTriFace1(false);
 327     FixedList<face,2> triFaces0(face(3)), triFaces1(face(3));
 328
 329     if (coupledModification_)
 330     {
 331         if (processorCoupledEntity(fIndex))
 332         {
 333             // When swapping across processors, we need to add the
 334             // prism-cell from (as well as delete on) the patchSubMesh
 335             const changeMap faceMap = insertCells(fIndex);
 336
 337             // Bail out if entity is handled elsewhere
 338             if (faceMap.type() == -2)
 339             {
 340                 return faceMap;
 341             }
 342
 343             if (faceMap.type() != 1)
 344             {
 345                 FatalErrorIn
 346                 (
 347                     "const changeMap dynamicTopoFvMesh::swapQuadFace"
 348                     "(const label fIndex)"
 349                 )
 350                     << " Could not insert cells around face: " << fIndex
 351                     << " :: " << faces_[fIndex] << nl
 352                     << abort(FatalError);
 353             }
 354
 355             // Figure out the new internal face index.
 356             // This should not be a coupled face anymore.
 357             label newIndex = faceMap.index();
 358
 359             // Temporarily turn off coupledModification
 360             unsetCoupledModification();
 361
 362             // Recursively call this function for the new face
 363             map = swapQuadFace(newIndex);
 364
 365             // Turn it back on.
 366             setCoupledModification();
 367
 368             return map;
 369         }
 370     }
 371
 372     // Get the two cells on either side...
 373     label c0 = owner_[fIndex];
 374     label c1 = neighbour_[fIndex];
 375
 376     // Prepare two new cells
 377     FixedList<cell, 2> newCell(cell(5));
 378
 379     // Need to find common faces and edges...
 380     // At the end of this loop, commonFaces [0] & [1] share commonEdge[0]
 381     // and commonFaces [2] & [3] share commonEdge[1]
 382     // Also, commonFaces[0] & [2] are connected to cell[0],
 383     // and commonFaces[1] & [3] are connected to cell[1]
 384     const labelList& fEdges = faceEdges_[fIndex];
 385
 386     forAll(fEdges, edgeI)
 387     {
 388         // Break out if all triangular faces are found
 389         if
 390         (
 391             foundTriFace0[0] && foundTriFace0[1] &&
 392             foundTriFace1[0] && foundTriFace1[1]
 393         )
 394         {
 395             break;
 396         }
 397
 398         // Obtain edgeFaces for this edge
 399         const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
 400
 401         forAll(eFaces, faceI)
 402         {
 403             const face& eFace = faces_[eFaces[faceI]];
 404
 405             if (eFace.size() == 3)
 406             {
 407                 // Found a triangular face. Determine which cell it belongs to.
 408                 if (owner_[eFaces[faceI]] == c0)
 409                 {
 410                     if (foundTriFace0[0])
 411                     {
 412                         // Update the second face on cell[0].
 413                         commonIndex = 2;
 414                         foundTriFace0[1] = true;
 415
 416                         if (foundTriFace1[1])
 417                         {
 418                             commonEdgeIndex[1] = fEdges[edgeI];
 419                             commonEdges[1] = edges_[fEdges[edgeI]];
 420                         }
 421                     }
 422                     else
 423                     {
 424                         // Update the first face on cell[0].
 425                         commonIndex = 0;
 426                         foundTriFace0[0] = true;
 427
 428                         if (foundTriFace1[0])
 429                         {
 430                             commonEdgeIndex[0] = fEdges[edgeI];
 431                             commonEdges[0] = edges_[fEdges[edgeI]];
 432                         }
 433                     }
 434                 }
 435                 else
 436                 {
 437                     if (foundTriFace1[0])
 438                     {
 439                         // Update the second face on cell[1].
 440                         commonIndex = 3;
 441                         foundTriFace1[1] = true;
 442
 443                         if (foundTriFace0[1])
 444                         {
 445                             commonEdgeIndex[1] = fEdges[edgeI];
 446                             commonEdges[1] = edges_[fEdges[edgeI]];
 447                         }
 448                     }
 449                     else
 450                     {
 451                         // Update the first face on cell[1].
 452                         commonIndex = 1;
 453                         foundTriFace1[0] = true;
 454
 455                         if (foundTriFace0[0])
 456                         {
 457                             commonEdgeIndex[0] = fEdges[edgeI];
 458                             commonEdges[0] = edges_[fEdges[edgeI]];
 459                         }
 460                     }
 461                 }
 462
 463                 // Store the face and index
 464                 commonFaces[commonIndex][0] = eFace[0];
 465                 commonFaces[commonIndex][1] = eFace[1];
 466                 commonFaces[commonIndex][2] = eFace[2];
 467
 468                 commonFaceIndex[commonIndex] = eFaces[faceI];
 469             }
 470         }
 471     }
 472
 473     // Find the interior/boundary faces.
 474     meshOps::findPrismFaces
 475     (
 476         fIndex,
 477         c0,
 478         faces_,
 479         cells_,
 480         neighbour_,
 481         c0BdyFace,
 482         c0BdyIndex,
 483         c0IntFace,
 484         c0IntIndex
 485     );
 486
 487     meshOps::findPrismFaces
 488     (
 489         fIndex,
 490         c1,
 491         faces_,
 492         cells_,
 493         neighbour_,
 494         c1BdyFace,
 495         c1BdyIndex,
 496         c1IntFace,
 497         c1IntIndex
 498     );
 499
 500     // Find the points that don't lie on shared edges
 501     // and the points next to them (for orientation)
 502     meshOps::findIsolatedPoint
 503     (
 504         commonFaces[1],
 505         commonEdges[0],
 506         otherPointIndex[1],
 507         nextToOtherPoint[1]
 508     );
 509
 510     meshOps::findIsolatedPoint
 511     (
 512         commonFaces[0],
 513         commonEdges[0],
 514         otherPointIndex[0],
 515         nextToOtherPoint[0]
 516     );
 517
 518     meshOps::findIsolatedPoint
 519     (
 520         commonFaces[2],
 521         commonEdges[1],
 522         otherPointIndex[2],
 523         nextToOtherPoint[2]
 524     );
 525
 526     meshOps::findIsolatedPoint
 527     (
 528         commonFaces[3],
 529         commonEdges[1],
 530         otherPointIndex[3],
 531         nextToOtherPoint[3]
 532     );
 533
 534     // Ensure that the configuration will be valid
 535     // using old point-positions. A simple area-based
 536     // calculation should suffice.
 537     FixedList<face, 2> triFaceOldPoints(face(3));
 538
 539     triFaceOldPoints[0][0] = otherPointIndex[0];
 540     triFaceOldPoints[0][1] = nextToOtherPoint[0];
 541     triFaceOldPoints[0][2] = otherPointIndex[1];
 542
 543     triFaceOldPoints[1][0] = otherPointIndex[1];
 544     triFaceOldPoints[1][1] = nextToOtherPoint[1];
 545     triFaceOldPoints[1][2] = otherPointIndex[0];
 546
 547     // Assume XY plane here
 548     vector n = vector(0,0,1);
 549
 550     forAll(triFaceOldPoints, faceI)
 551     {
 552         // Assume that centre-plane passes through the origin.
 553         vector xf, nf;
 554
 555         xf = triFaceOldPoints[faceI].centre(oldPoints_);
 556         nf = triFaceOldPoints[faceI].normal(oldPoints_);
 557
 558         if ((((xf & n) * n) & nf) < 0.0)
 559         {
 560             // This will yield an inverted cell. Bail out.
 561             return map;
 562         }
 563     }
 564
 565     // Find the other two edges on the face being flipped
 566     forAll(fEdges, edgeI)
 567     {
 568         if
 569         (
 570             fEdges[edgeI] != commonEdgeIndex[0] &&
 571             fEdges[edgeI] != commonEdgeIndex[1]
 572         )
 573         {
 574             // Obtain a reference to this edge
 575             const edge& eThis = edges_[fEdges[edgeI]];
 576
 577             if
 578             (
 579                 eThis[0] == nextToOtherPoint[0]
 580              || eThis[1] == nextToOtherPoint[0]
 581             )
 582             {
 583                 otherEdgeIndex[0] = fEdges[edgeI];
 584             }
 585             else
 586             {
 587                 otherEdgeIndex[1] = fEdges[edgeI];
 588             }
 589         }
 590     }
 591
 592     // At the end of this loop, commonIntFaces [0] & [1] share otherEdges[0]
 593     // and commonIntFaces [2] & [3] share the otherEdges[1],
 594     // where [0],[2] lie on cell[0] and [1],[3] lie on cell[1]
 595     found = false;
 596
 597     const labelList& e1 = faceEdges_[c0IntIndex[0]];
 598
 599     forAll(e1,edgeI)
 600     {
 601         if (e1[edgeI] == otherEdgeIndex[0])
 602         {
 603             commonIntFaces[0] = c0IntFace[0];
 604             commonIntFaces[2] = c0IntFace[1];
 605             commonIntFaceIndex[0] = c0IntIndex[0];
 606             commonIntFaceIndex[2] = c0IntIndex[1];
 607             found = true; break;
 608         }
 609     }
 610
 611     if (!found)
 612     {
 613         // The edge was not found before
 614         commonIntFaces[0] = c0IntFace[1];
 615         commonIntFaces[2] = c0IntFace[0];
 616         commonIntFaceIndex[0] = c0IntIndex[1];
 617         commonIntFaceIndex[2] = c0IntIndex[0];
 618     }
 619
 620     found = false;
 621
 622     const labelList& e3 = faceEdges_[c1IntIndex[0]];
 623
 624     forAll(e3,edgeI)
 625     {
 626         if (e3[edgeI] == otherEdgeIndex[0])
 627         {
 628             commonIntFaces[1] = c1IntFace[0];
 629             commonIntFaces[3] = c1IntFace[1];
 630             commonIntFaceIndex[1] = c1IntIndex[0];
 631             commonIntFaceIndex[3] = c1IntIndex[1];
 632             found = true; break;
 633         }
 634     }
 635
 636     if (!found)
 637     {
 638         // The edge was not found before
 639         commonIntFaces[1] = c1IntFace[1];
 640         commonIntFaces[3] = c1IntFace[0];
 641         commonIntFaceIndex[1] = c1IntIndex[1];
 642         commonIntFaceIndex[3] = c1IntIndex[0];
 643     }
 644
 645     // Find two common edges between quad/quad faces
 646     meshOps::findCommonEdge
 647     (
 648         c0IntIndex[0],
 649         c0IntIndex[1],
 650         faceEdges_,
 651         otherEdgeIndex[2]
 652     );
 653
 654     meshOps::findCommonEdge
 655     (
 656         c1IntIndex[0],
 657         c1IntIndex[1],
 658         faceEdges_,
 659         otherEdgeIndex[3]
 660     );
 661
 662     // Find four common edges between quad/tri faces
 663     meshOps::findCommonEdge
 664     (
 665         commonFaceIndex[1],
 666         commonIntFaceIndex[1],
 667         faceEdges_,
 668         cornerEdgeIndex[0]
 669     );
 670
 671     meshOps::findCommonEdge
 672     (
 673         commonFaceIndex[3],
 674         commonIntFaceIndex[1],
 675         faceEdges_,
 676         cornerEdgeIndex[1]
 677     );
 678
 679     meshOps::findCommonEdge
 680     (
 681         commonFaceIndex[0],
 682         commonIntFaceIndex[2],
 683         faceEdges_,
 684         cornerEdgeIndex[2]
 685     );
 686
 687     meshOps::findCommonEdge
 688     (
 689         commonFaceIndex[2],
 690         commonIntFaceIndex[2],
 691         faceEdges_,
 692         cornerEdgeIndex[3]
 693     );
 694
 695     // Perform a few debug calls before modifications
 696     if (debug > 1)
 697     {
 698         Pout<< nl << nl << "Face: " << fIndex
 699             << " needs to be flipped. " << endl;
 700
 701         if (debug > 2)
 702         {
 703             Pout<< " Cell[0]: " << c0 << ": " << cells_[c0] << nl
 704                 << " Cell[1]: " << c1 << ": " << cells_[c1] << nl;
 705
 706             Pout<< " Common Faces: Set 1: "
 707                 << commonFaceIndex[0] << ": " << commonFaces[0] << ", "
 708                 << commonFaceIndex[1] << ": " << commonFaces[1] << nl;
 709
 710             Pout<< " Common Faces: Set 2: "
 711                 << commonFaceIndex[2] << ": " << commonFaces[2] << ", "
 712                 << commonFaceIndex[3] << ": " << commonFaces[3] << nl;
 713
 714             Pout<< " Old face: " << faces_[fIndex] << nl
 715                 << " Old faceEdges: " << faceEdges_[fIndex] << endl;
 716         }
 717
 718         // Write out VTK files before change
 719         if (debug > 3)
 720         {
 721             labelList cellHull(2, -1);
 722
 723             cellHull[0] = c0;
 724             cellHull[1] = c1;
 725
 726             writeVTK(Foam::name(fIndex) + "_Swap_0", cellHull);
 727         }
 728     }
 729
 730     // Modify the five faces belonging to this hull
 731     face newFace = faces_[fIndex];
 732     labelList newFEdges = faceEdges_[fIndex];
 733     FixedList<face, 4> newBdyFace(face(3));
 734     FixedList<edge, 2> newEdges;
 735
 736     // Assign current faces
 737     forAll(newBdyFace, faceI)
 738     {
 739         newBdyFace[faceI] = faces_[commonFaceIndex[faceI]];
 740     }
 741
 742     label c0count=0, c1count=0;
 743
 744     // Size down edgeFaces for the original face
 745     meshOps::sizeDownList
 746     (
 747         fIndex,
 748         edgeFaces_[otherEdgeIndex[0]]
 749     );
 750
 751     meshOps::sizeDownList
 752     (
 753         fIndex,
 754         edgeFaces_[otherEdgeIndex[1]]
 755     );
 756
 757     // Size up edgeFaces for the face after flipping
 758     meshOps::sizeUpList
 759     (
 760         fIndex,
 761         edgeFaces_[otherEdgeIndex[2]]
 762     );
 763
 764     meshOps::sizeUpList
 765     (
 766         fIndex,
 767         edgeFaces_[otherEdgeIndex[3]]
 768     );
 769
 770     // Replace edgeFaces and faceEdges for the
 771     // four (out of 8 total) corner edges of this hull.
 772     meshOps::replaceLabel
 773     (
 774         cornerEdgeIndex[0],
 775         cornerEdgeIndex[2],
 776         faceEdges_[commonFaceIndex[1]]
 777     );
 778
 779     meshOps::replaceLabel
 780     (
 781         commonFaceIndex[1],
 782         commonFaceIndex[0],
 783         edgeFaces_[cornerEdgeIndex[0]]
 784     );
 785
 786     meshOps::replaceLabel
 787     (
 788         cornerEdgeIndex[1],
 789         cornerEdgeIndex[3],
 790         faceEdges_[commonFaceIndex[3]]
 791     );
 792
 793     meshOps::replaceLabel
 794     (
 795         commonFaceIndex[3],
 796         commonFaceIndex[2],
 797         edgeFaces_[cornerEdgeIndex[1]]
 798     );
 799
 800     meshOps::replaceLabel
 801     (
 802         cornerEdgeIndex[2],
 803         cornerEdgeIndex[0],
 804         faceEdges_[commonFaceIndex[0]]
 805     );
 806
 807     meshOps::replaceLabel
 808     (
 809         commonFaceIndex[0],
 810         commonFaceIndex[1],
 811         edgeFaces_[cornerEdgeIndex[2]]
 812     );
 813
 814     meshOps::replaceLabel
 815     (
 816         cornerEdgeIndex[3],
 817         cornerEdgeIndex[1],
 818         faceEdges_[commonFaceIndex[2]]
 819     );
 820
 821     meshOps::replaceLabel
 822     (
 823         commonFaceIndex[2],
 824         commonFaceIndex[3],
 825         edgeFaces_[cornerEdgeIndex[3]]
 826     );
 827
 828     // Define parameters for the new flipped face
 829     newFace[0] = otherPointIndex[0];
 830     newFace[1] = otherPointIndex[1];
 831     newFace[2] = otherPointIndex[3];
 832     newFace[3] = otherPointIndex[2];
 833
 834     newFEdges[0] = otherEdgeIndex[2];
 835     newFEdges[1] = commonEdgeIndex[0];
 836     newFEdges[2] = otherEdgeIndex[3];
 837     newFEdges[3] = commonEdgeIndex[1];
 838
 839     newCell[0][c0count++] = fIndex;
 840     newCell[1][c1count++] = fIndex;
 841
 842     owner_[fIndex] = c0;
 843     neighbour_[fIndex] = c1;
 844
 845     // Both commonEdges need to be renumbered.
 846     newEdges[0][0] = otherPointIndex[0];
 847     newEdges[0][1] = otherPointIndex[1];
 848
 849     newEdges[1][0] = otherPointIndex[3];
 850     newEdges[1][1] = otherPointIndex[2];
 851
 852     // Four modified boundary faces need to be constructed,
 853     // but right-handedness is also important.
 854     // Take a cue from the existing boundary-face orientation
 855
 856     // Zeroth boundary face - Owner c[0], Neighbour -1
 857     newBdyFace[0][0] = otherPointIndex[0];
 858     newBdyFace[0][1] = nextToOtherPoint[0];
 859     newBdyFace[0][2] = otherPointIndex[1];
 860     newCell[0][c0count++] = commonFaceIndex[0];
 861     owner_[commonFaceIndex[0]] = c0;
 862     neighbour_[commonFaceIndex[0]] = -1;
 863
 864     // First boundary face - Owner c[1], Neighbour -1
 865     newBdyFace[1][0] = otherPointIndex[1];
 866     newBdyFace[1][1] = nextToOtherPoint[1];
 867     newBdyFace[1][2] = otherPointIndex[0];
 868     newCell[1][c1count++] = commonFaceIndex[1];
 869     owner_[commonFaceIndex[1]] = c1;
 870     neighbour_[commonFaceIndex[1]] = -1;
 871
 872     // Second boundary face - Owner c[0], Neighbour -1
 873     newBdyFace[2][0] = otherPointIndex[3];
 874     newBdyFace[2][1] = nextToOtherPoint[3];
 875     newBdyFace[2][2] = otherPointIndex[2];
 876     newCell[0][c0count++] = commonFaceIndex[2];
 877     owner_[commonFaceIndex[2]] = c0;
 878     neighbour_[commonFaceIndex[2]] = -1;
 879
 880     // Third boundary face - Owner c[1], Neighbour -1
 881     newBdyFace[3][0] = otherPointIndex[2];
 882     newBdyFace[3][1] = nextToOtherPoint[2];
 883     newBdyFace[3][2] = otherPointIndex[3];
 884     newCell[1][c1count++] = commonFaceIndex[3];
 885     owner_[commonFaceIndex[3]] = c1;
 886     neighbour_[commonFaceIndex[3]] = -1;
 887
 888     if (debug > 1)
 889     {
 890         Pout<< "New flipped face: " << newFace << nl;
 891
 892         if (debug > 2)
 893         {
 894             forAll(newBdyFace, faceI)
 895             {
 896                 Pout<< "New boundary face[" << faceI << "]: "
 897                     << commonFaceIndex[faceI]
 898                     << ": " << newBdyFace[faceI] << nl;
 899             }
 900         }
 901
 902         Pout<< endl;
 903     }
 904
 905     // Check the orientation of the two quad faces, and modify as necessary
 906     label newOwn=0, newNei=0;
 907
 908     // The quad face belonging to cell[1] now becomes a part of cell[0]
 909     if (neighbour_[commonIntFaceIndex[1]] == -1)
 910     {
 911         // Boundary face
 912         // Face doesn't need to be flipped, just update the owner
 913         f = commonIntFaces[1];
 914         newOwn = c0;
 915         newNei = -1;
 916     }
 917     else
 918     if (owner_[commonIntFaceIndex[1]] == c1)
 919     {
 920         // This face is on the interior, check for previous owner
 921         // Upper-triangular ordering has to be maintained, however...
 922         if (c0 > neighbour_[commonIntFaceIndex[1]])
 923         {
 924             // Flip is necessary
 925             f = commonIntFaces[1].reverseFace();
 926             newOwn = neighbour_[commonIntFaceIndex[1]];
 927             newNei = c0;
 928
 929             setFlip(commonIntFaceIndex[1]);
 930         }
 931         else
 932         {
 933             // Flip isn't necessary, just change the owner
 934             f = commonIntFaces[1];
 935             newOwn = c0;
 936             newNei = neighbour_[commonIntFaceIndex[1]];
 937         }
 938     }
 939     else
 940     if (neighbour_[commonIntFaceIndex[1]] == c1)
 941     {
 942         // This face is on the interior, check for previous neighbour
 943         // Upper-triangular ordering has to be maintained, however...
 944         if (c0 < owner_[commonIntFaceIndex[1]])
 945         {
 946             // Flip is necessary
 947             f = commonIntFaces[1].reverseFace();
 948             newOwn = c0;
 949             newNei = owner_[commonIntFaceIndex[1]];
 950
 951             setFlip(commonIntFaceIndex[1]);
 952         }
 953         else
 954         {
 955             // Flip isn't necessary, just change the neighbour
 956             f = commonIntFaces[1];
 957             newOwn = owner_[commonIntFaceIndex[1]];
 958             newNei = c0;
 959         }
 960     }
 961
 962     faces_[commonIntFaceIndex[1]] = f;
 963     newCell[0][c0count++] = commonIntFaceIndex[0];
 964     newCell[0][c0count++] = commonIntFaceIndex[1];
 965     owner_[commonIntFaceIndex[1]] = newOwn;
 966     neighbour_[commonIntFaceIndex[1]] = newNei;
 967
 968     // The quad face belonging to cell[0] now becomes a part of cell[1]
 969     if (neighbour_[commonIntFaceIndex[2]] == -1)
 970     {
 971         // Boundary face
 972         // Face doesn't need to be flipped, just update the owner
 973         f = commonIntFaces[2];
 974         newOwn = c1;
 975         newNei = -1;
 976     }
 977     else
 978     if (owner_[commonIntFaceIndex[2]] == c0)
 979     {
 980         // This face is on the interior, check for previous owner
 981         // Upper-triangular ordering has to be maintained, however...
 982         if (c1 > neighbour_[commonIntFaceIndex[2]])
 983         {
 984             // Flip is necessary
 985             f = commonIntFaces[2].reverseFace();
 986             newOwn = neighbour_[commonIntFaceIndex[2]];
 987             newNei = c1;
 988
 989             setFlip(commonIntFaceIndex[2]);
 990         }
 991         else
 992         {
 993             // Flip isn't necessary, just change the owner
 994             f = commonIntFaces[2];
 995             newOwn = c1;
 996             newNei = neighbour_[commonIntFaceIndex[2]];
 997         }
 998     }
 999     else
1000     if (neighbour_[commonIntFaceIndex[2]] == c0)
1001     {
1002         // This face is on the interior, check for previous neighbour
1003         // Upper-triangular ordering has to be maintained, however...
1004         if (c1 < owner_[commonIntFaceIndex[2]])
1005         {
1006             // Flip is necessary
1007             f = commonIntFaces[2].reverseFace();
1008             newOwn = c1;
1009             newNei = owner_[commonIntFaceIndex[2]];
1010
1011             setFlip(commonIntFaceIndex[2]);
1012         }
1013         else
1014         {
1015             // Flip isn't necessary, just change the neighbour
1016             f = commonIntFaces[2];
1017             newOwn = owner_[commonIntFaceIndex[2]];
1018             newNei = c1;
1019         }
1020     }
1021
1022     faces_[commonIntFaceIndex[2]] = f;
1023     newCell[1][c1count++] = commonIntFaceIndex[2];
1024     newCell[1][c1count++] = commonIntFaceIndex[3];
1025     owner_[commonIntFaceIndex[2]] = newOwn;
1026     neighbour_[commonIntFaceIndex[2]] = newNei;
1027
1028     // Now that all entities are properly configured,
1029     // overwrite the entries in connectivity lists.
1030     cells_[c0] = newCell[0];
1031     cells_[c1] = newCell[1];
1032
1033     faces_[fIndex] = newFace;
1034     faceEdges_[fIndex] = newFEdges;
1035
1036     forAll(newBdyFace, faceI)
1037     {
1038         faces_[commonFaceIndex[faceI]] = newBdyFace[faceI];
1039     }
1040
1041     edges_[commonEdgeIndex[0]] = newEdges[0];
1042     edges_[commonEdgeIndex[1]] = newEdges[1];
1043
1044     // Fill-in mapping information
1045     labelList mC(2, -1);
1046     mC[0] = c0; mC[1] = c1;
1047
1048     forAll(mC, cellI)
1049     {
1050         // Set the mapping for this cell
1051         setCellMapping(mC[cellI], mC);
1052     }
1053
1054     // Interpolate new fluxes for the flipped face.
1055     setFaceMapping(fIndex);
1056
1057     // Write out VTK files after change
1058     if (debug > 3)
1059     {
1060         labelList cellHull(2, -1);
1061
1062         cellHull[0] = c0;
1063         cellHull[1] = c1;
1064
1065         writeVTK
1066         (
1067             Foam::name(fIndex)
1068           + "_Swap_1",
1069             cellHull
1070         );
1071     }
1072
1073     // Set the flag
1074     topoChangeFlag_ = true;
1075
1076     // Increment the counter
1077     statistics_[1]++;
1078
1079     // Return a successful operation.
1080     map.type() = 1;
1081
1082     return map;
1083 }
1084
1085
1086 // Allocate dynamic programming tables
1087 void dynamicTopoFvMesh::initTables
1088 (
1089     labelList& m,
1090     PtrList<scalarListList>& Q,
1091     PtrList<labelListList>& K,
1092     PtrList<labelListList>& triangulations,
1093     const label checkIndex
1094 ) const
1095 {
1096     label mMax = maxTetsPerEdge_;
1097
1098     // Check if resizing is necessary only for a particular index.
1099     if (checkIndex != -1)
1100     {
1101         m[checkIndex] = -1;
1102         Q[checkIndex].setSize((mMax-2),scalarList(mMax,-1.0));
1103         K[checkIndex].setSize((mMax-2),labelList(mMax,-1));
1104         triangulations[checkIndex].setSize(3,labelList((mMax-2),-1));
1105
1106         return;
1107     }
1108
1109     // Size all elements by default.
1110     label numIndices = coupledModification_ ? 2 : 1;
1111
1112     m.setSize(numIndices, -1);
1113     Q.setSize(numIndices);
1114     K.setSize(numIndices);
1115     triangulations.setSize(numIndices);
1116
1117     forAll(Q, indexI)
1118     {
1119         Q.set
1120         (
1121             indexI,
1122             new scalarListList((mMax-2),scalarList(mMax,-1.0))
1123         );
1124
1125         K.set
1126         (
1127             indexI,
1128             new labelListList((mMax-2),labelList(mMax,-1))
1129         );
1130
1131         triangulations.set
1132         (
1133             indexI,
1134             new labelListList(3,labelList((mMax-2),-1))
1135         );
1136     }
1137 }
1138
1139
1140 // Utility method to fill the dynamic programming tables
1141 //  - Returns true if the operation completed successfully.
1142 //  - Returns false if tables could not be resized.
1143 bool dynamicTopoFvMesh::fillTables
1144 (
1145     const label eIndex,
1146     scalar& minQuality,
1147     labelList& m,
1148     labelList& hullVertices,
1149     PtrList<scalarListList>& Q,
1150     PtrList<labelListList>& K,
1151     PtrList<labelListList>& triangulations,
1152     const label checkIndex
1153 ) const
1154 {
1155     // Obtain a reference to this edge
1156     const labelList& edgeFaces = edgeFaces_[eIndex];
1157
1158     // If this entity was deleted, skip it.
1159     if (edgeFaces.empty())
1160     {
1161         return false;
1162     }
1163
1164     // Ensure that edge is surrounded by triangles
1165     forAll(edgeFaces, faceI)
1166     {
1167         if (faces_[edgeFaces[faceI]].size() != 3)
1168         {
1169             return false;
1170         }
1171     }
1172
1173     const edge& edgeToCheck = edges_[eIndex];
1174
1175     if (coupledModification_)
1176     {
1177         return coupledFillTables(eIndex, minQuality, m, Q, K, triangulations);
1178     }
1179
1180     // Fill in the size
1181     m[checkIndex] = hullVertices.size();
1182
1183     // Check if a table-resize is necessary
1184     if (m[checkIndex] > maxTetsPerEdge_)
1185     {
1186         if (allowTableResize_)
1187         {
1188             // Resize the tables to account for
1189             // more tets per edge
1190             label& mtpe = const_cast<label&>(maxTetsPerEdge_);
1191
1192             mtpe = m[checkIndex];
1193
1194             // Clear tables for this index.
1195             Q[checkIndex].clear();
1196             K[checkIndex].clear();
1197             triangulations[checkIndex].clear();
1198
1199             // Resize for this index.
1200             initTables(m, Q, K, triangulations, checkIndex);
1201         }
1202         else
1203         {
1204             // Can't resize. Bail out.
1205             return false;
1206         }
1207     }
1208
1209     // Fill dynamic programming tables
1210     fillTables
1211     (
1212         edgeToCheck,
1213         minQuality,
1214         m[checkIndex],
1215         hullVertices,
1216         points_,
1217         Q[checkIndex],
1218         K[checkIndex],
1219         triangulations[checkIndex]
1220     );
1221
1222     // Print out tables for debugging
1223     if (debug > 5)
1224     {
1225         printTables(m, Q, K, checkIndex);
1226     }
1227
1228     return true;
1229 }
1230
1231
1232 // Fill tables given addressing
1233 void dynamicTopoFvMesh::fillTables
1234 (
1235     const edge& edgeToCheck,
1236     const scalar minQuality,
1237     const label m,
1238     const labelList& hullVertices,
1239     const UList<point>& points,
1240     scalarListList& Q,
1241     labelListList& K,
1242     labelListList& triangulations
1243 ) const
1244 {
1245     for (label i = (m - 3); i >= 0; i--)
1246     {
1247         for (label j = i + 2; j < m; j++)
1248         {
1249             for (label k = i + 1; k < j; k++)
1250             {
1251                 scalar q = (*tetMetric_)
1252                 (
1253                     points[hullVertices[i]],
1254                     points[hullVertices[k]],
1255                     points[hullVertices[j]],
1256                     points[edgeToCheck[0]]
1257                 );
1258
1259                 // For efficiency, check the bottom triangulation
1260                 // only when the top one if less than the hull quality.
1261                 if (q > minQuality)
1262                 {
1263                     q =
1264                     (
1265                         Foam::min
1266                         (
1267                             q,
1268                             (*tetMetric_)
1269                             (
1270                                 points[hullVertices[j]],
1271                                 points[hullVertices[k]],
1272                                 points[hullVertices[i]],
1273                                 points[edgeToCheck[1]]
1274                             )
1275                         )
1276                     );
1277                 }
1278
1279                 if (k < j - 1)
1280                 {
1281                     q = Foam::min(q, Q[k][j]);
1282                 }
1283
1284                 if (k > i + 1)
1285                 {
1286                     q = Foam::min(q, Q[i][k]);
1287                 }
1288
1289                 if ((k == i + 1) || (q > Q[i][j]))
1290                 {
1291                     Q[i][j] = q;
1292                     K[i][j] = k;
1293                 }
1294             }
1295         }
1296     }
1297 }
1298
1299
1300 // Remove the edge according to the swap sequence.
1301 // - Returns a changeMap with a type specifying:
1302 //     1: Swap sequence was successful
1303 //    -1: Swap sequence failed
1304 const changeMap dynamicTopoFvMesh::removeEdgeFlips
1305 (
1306     const label eIndex,
1307     const scalar minQuality,
1308     const labelList& vertexHull,
1309     PtrList<scalarListList>& Q,
1310     PtrList<labelListList>& K,
1311     PtrList<labelListList>& triangulations,
1312     const label checkIndex
1313 )
1314 {
1315     changeMap map, slaveMap;
1316
1317     if (debug > 2)
1318     {
1319         Pout<< nl << " Removing edge : " << eIndex << " by flipping."
1320             << " Edge: " << edges_[eIndex]
1321             << " minQuality: " << minQuality
1322             << endl;
1323     }
1324
1325     if (coupledModification_)
1326     {
1327         if (processorCoupledEntity(eIndex))
1328         {
1329             const changeMap edgeMap = insertCells(eIndex);
1330
1331             // Bail out if entity is handled elsewhere
1332             if (edgeMap.type() == -2)
1333             {
1334                 return edgeMap;
1335             }
1336
1337             if (edgeMap.type() != 1)
1338             {
1339                 FatalErrorIn
1340                 (
1341                     "\n"
1342                     "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1343                     "(\n"
1344                     "    const label eIndex,\n"
1345                     "    const scalar minQuality,\n"
1346                     "    const labelList& vertexHull,\n"
1347                     "    PtrList<scalarListList>& Q,\n"
1348                     "    PtrList<labelListList>& K,\n"
1349                     "    PtrList<labelListList>& triangulations,\n"
1350                     "    const label checkIndex\n"
1351                     ")\n"
1352                 )
1353                     << " Could not insert cells around edge: " << eIndex
1354                     << " :: " << edges_[eIndex] << nl
1355                     << abort(FatalError);
1356             }
1357
1358             // Temporarily turn off coupledModification
1359             unsetCoupledModification();
1360
1361             // Re-fill tables for the reconfigured edge.
1362             // This should not be a coupled edge anymore.
1363             label newIndex = edgeMap.index();
1364
1365             if (processorCoupledEntity(newIndex))
1366             {
1367                 // Write out edge connectivity
1368                 writeEdgeConnectivity(newIndex);
1369
1370                 FatalErrorIn
1371                 (
1372                     "\n"
1373                     "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1374                     "(\n"
1375                     "    const label eIndex,\n"
1376                     "    const scalar minQuality,\n"
1377                     "    const labelList& vertexHull,\n"
1378                     "    PtrList<scalarListList>& Q,\n"
1379                     "    PtrList<labelListList>& K,\n"
1380                     "    PtrList<labelListList>& triangulations,\n"
1381                     "    const label checkIndex\n"
1382                     ")\n"
1383                 )
1384                     << " Edge: " << newIndex
1385                     << " :: " << edges_[newIndex] << nl
1386                     << " is still processor-coupled. "
1387                     << abort(FatalError);
1388             }
1389
1390             const edge& newEdge = edges_[newIndex];
1391
1392             // Build vertexHull for this edge
1393             labelList newVertexHull;
1394             buildVertexHull(newIndex, newVertexHull);
1395
1396             fillTables
1397             (
1398                 newEdge,
1399                 minQuality,
1400                 newVertexHull.size(),
1401                 newVertexHull,
1402                 points_,
1403                 Q[0],
1404                 K[0],
1405                 triangulations[0]
1406             );
1407
1408             // Recursively call this function for the new edge
1409             map =
1410             (
1411                 removeEdgeFlips
1412                 (
1413                     newIndex,
1414                     minQuality,
1415                     newVertexHull,
1416                     Q,
1417                     K,
1418                     triangulations
1419                 )
1420             );
1421
1422             // Turn it back on.
1423             setCoupledModification();
1424
1425             return map;
1426         }
1427     }
1428
1429     label m = vertexHull.size();
1430     labelList hullFaces(m, -1);
1431     labelList hullCells(m, -1);
1432     labelList hullEdges(m, -1);
1433     labelListList ringEntities(4, labelList(m, -1));
1434
1435     // Construct the hull
1436     meshOps::constructHull
1437     (
1438         eIndex,
1439         faces_,
1440         edges_,
1441         cells_,
1442         owner_,
1443         neighbour_,
1444         faceEdges_,
1445         edgeFaces_,
1446         vertexHull,
1447         hullEdges,
1448         hullFaces,
1449         hullCells,
1450         ringEntities
1451     );
1452
1453     label numTriangulations = 0, isolatedVertex = -1;
1454
1455     // Extract the appropriate triangulations
1456     extractTriangulation
1457     (
1458         0,
1459         m-1,
1460         K[checkIndex],
1461         numTriangulations,
1462         triangulations[checkIndex]
1463     );
1464
1465     // Check old-volumes for the configuration.
1466     if
1467     (
1468         checkTriangulationVolumes
1469         (
1470             eIndex,
1471             vertexHull,
1472             triangulations[checkIndex]
1473         )
1474     )
1475     {
1476         // Reset all triangulations and bail out
1477         triangulations[checkIndex][0] = -1;
1478         triangulations[checkIndex][1] = -1;
1479         triangulations[checkIndex][2] = -1;
1480
1481         return map;
1482     }
1483
1484     // Determine the final swap triangulation
1485     label tF =
1486     (
1487         identify32Swap
1488         (
1489             eIndex,
1490             vertexHull,
1491             triangulations[checkIndex]
1492         )
1493     );
1494
1495     // Check that the triangulation is valid
1496     label pIndex = -1;
1497
1498     if (tF == -1)
1499     {
1500         const edge& edgeToCheck = edges_[eIndex];
1501
1502         Pout<< " All triangulations: " << nl
1503             << ' ' << triangulations[checkIndex][0] << nl
1504             << ' ' << triangulations[checkIndex][1] << nl
1505             << ' ' << triangulations[checkIndex][2] << nl
1506             << endl;
1507
1508         FatalErrorIn
1509         (
1510             "\n"
1511             "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1512             "(\n"
1513             "    const label eIndex,\n"
1514             "    const scalar minQuality,\n"
1515             "    const labelList& vertexHull,\n"
1516             "    PtrList<scalarListList>& Q,\n"
1517             "    PtrList<labelListList>& K,\n"
1518             "    PtrList<labelListList>& triangulations,\n"
1519             "    const label checkIndex\n"
1520             ")\n"
1521         )
1522             << "Could not determine 3-2 swap triangulation." << nl
1523             << "Edge: " << edgeToCheck << nl
1524             << "Edge Points: "
1525             << points_[edgeToCheck[0]] << ","
1526             << points_[edgeToCheck[1]] << nl
1527             << abort(FatalError);
1528
1529         // Reset all triangulations and bail out
1530         triangulations[checkIndex][0] = -1;
1531         triangulations[checkIndex][1] = -1;
1532         triangulations[checkIndex][2] = -1;
1533
1534         return map;
1535     }
1536
1537     if (debug > 2)
1538     {
1539         Pout<< " Identified tF as: " << tF << nl;
1540
1541         Pout<< " Triangulation: "
1542             << triangulations[checkIndex][0][tF] << " "
1543             << triangulations[checkIndex][1][tF] << " "
1544             << triangulations[checkIndex][2][tF] << " "
1545             << nl;
1546
1547         Pout<< " All triangulations: " << nl
1548             << ' ' << triangulations[checkIndex][0] << nl
1549             << ' ' << triangulations[checkIndex][1] << nl
1550             << ' ' << triangulations[checkIndex][2] << nl
1551             << endl;
1552     }
1553
1554     if (coupledModification_)
1555     {
1556         if (locallyCoupledEntity(eIndex))
1557         {
1558             // Flip the slave edge as well.
1559             label sIndex = -1;
1560
1561             // Determine the slave index.
1562             forAll(patchCoupling_, patchI)
1563             {
1564                 if (patchCoupling_(patchI))
1565                 {
1566                     const label edgeEnum  = coupleMap::EDGE;
1567                     const coupleMap& cMap = patchCoupling_[patchI].map();
1568
1569                     if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
1570                     {
1571                         break;
1572                     }
1573                 }
1574             }
1575
1576             if (sIndex == -1)
1577             {
1578                 FatalErrorIn
1579                 (
1580                     "\n"
1581                     "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1582                     "(\n"
1583                     "    const label eIndex,\n"
1584                     "    const scalar minQuality,\n"
1585                     "    const labelList& vertexHull,\n"
1586                     "    PtrList<scalarListList>& Q,\n"
1587                     "    PtrList<labelListList>& K,\n"
1588                     "    PtrList<labelListList>& triangulations,\n"
1589                     "    const label checkIndex\n"
1590                     ")\n"
1591                 )
1592                     << "Coupled maps were improperly specified." << nl
1593                     << " Slave index not found for: " << nl
1594                     << " Edge: " << eIndex << nl
1595                     << abort(FatalError);
1596             }
1597
1598             if (debug > 2)
1599             {
1600                 Pout<< nl << "Removing slave edge: " << sIndex
1601                     << " for master edge: " << eIndex << endl;
1602             }
1603
1604             // Build vertexHull for this edge
1605             labelList slaveVertexHull;
1606             buildVertexHull(sIndex, slaveVertexHull);
1607
1608             // Turn off switch temporarily.
1609             unsetCoupledModification();
1610
1611             // Recursively call for the slave edge.
1612             slaveMap =
1613             (
1614                 removeEdgeFlips
1615                 (
1616                     sIndex,
1617                     minQuality,
1618                     slaveVertexHull,
1619                     Q,
1620                     K,
1621                     triangulations,
1622                     1
1623                 )
1624             );
1625
1626             // Turn it back on.
1627             setCoupledModification();
1628
1629             // Bail out if the slave failed.
1630             if (slaveMap.type() == -1)
1631             {
1632                 // Reset all triangulations and bail out
1633                 triangulations[checkIndex][0] = -1;
1634                 triangulations[checkIndex][1] = -1;
1635                 triangulations[checkIndex][2] = -1;
1636
1637                 return slaveMap;
1638             }
1639         }
1640     }
1641
1642     // Perform a series of 2-3 swaps
1643     label numSwaps = 0;
1644
1645     while (numSwaps < (m-3))
1646     {
1647         for (label i = 0; i < (m-2); i++)
1648         {
1649             if ( (i != tF) && (triangulations[checkIndex][0][i] != -1) )
1650             {
1651                 // Check if triangulation is on the boundary
1652                 if
1653                 (
1654                     boundaryTriangulation
1655                     (
1656                         i,
1657                         isolatedVertex,
1658                         triangulations[checkIndex]
1659                     )
1660                 )
1661                 {
1662                     // Perform 2-3 swap
1663                     map =
1664                     (
1665                         swap23
1666                         (
1667                             isolatedVertex,
1668                             eIndex,
1669                             i,
1670                             numTriangulations,
1671                             triangulations[checkIndex],
1672                             vertexHull,
1673                             hullFaces,
1674                             hullCells
1675                         )
1676                     );
1677
1678                     // Done with this face, so reset it
1679                     triangulations[checkIndex][0][i] = -1;
1680                     triangulations[checkIndex][1][i] = -1;
1681                     triangulations[checkIndex][2][i] = -1;
1682
1683                     numSwaps++;
1684                 }
1685             }
1686         }
1687
1688         if (numSwaps == 0)
1689         {
1690             Pout<< "Triangulations: " << nl;
1691
1692             forAll(triangulations[checkIndex], row)
1693             {
1694                 Pout<< triangulations[checkIndex][row] << nl;
1695             }
1696
1697             Pout<<endl;
1698
1699             // Should have performed at least one swap
1700             FatalErrorIn
1701             (
1702                 "\n"
1703                 "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1704                 "(\n"
1705                 "    const label eIndex,\n"
1706                 "    const scalar minQuality,\n"
1707                 "    const labelList& vertexHull,\n"
1708                 "    PtrList<scalarListList>& Q,\n"
1709                 "    PtrList<labelListList>& K,\n"
1710                 "    PtrList<labelListList>& triangulations,\n"
1711                 "    const label checkIndex\n"
1712                 ")\n"
1713             )
1714                 << "Did not perform any 2-3 swaps" << nl
1715                 << abort(FatalError);
1716         }
1717     }
1718
1719     // Perform the final 3-2 / 2-2 swap
1720     map =
1721     (
1722         swap32
1723         (
1724             eIndex,
1725             tF,
1726             numTriangulations,
1727             triangulations[checkIndex],
1728             vertexHull,
1729             hullFaces,
1730             hullCells
1731         )
1732     );
1733
1734     // Done with this face, so reset it
1735     triangulations[checkIndex][0][tF] = -1;
1736     triangulations[checkIndex][1][tF] = -1;
1737     triangulations[checkIndex][2][tF] = -1;
1738
1739     // Update the coupled map
1740     if (coupledModification_)
1741     {
1742         // Create a mapping entry for the new edge.
1743         const coupleMap& cMap = patchCoupling_[pIndex].map();
1744
1745         if (locallyCoupledEntity(map.addedEdgeList()[0].index()))
1746         {
1747             cMap.mapSlave
1748             (
1749                 coupleMap::EDGE,
1750                 map.addedEdgeList()[0].index(),
1751                 slaveMap.addedEdgeList()[0].index()
1752             );
1753
1754             cMap.mapMaster
1755             (
1756                 coupleMap::EDGE,
1757                 slaveMap.addedEdgeList()[0].index(),
1758                 map.addedEdgeList()[0].index()
1759             );
1760         }
1761
1762         // Add a mapping entry for two new faces as well.
1763         face cF(3);
1764
1765         const List<objectMap>& amfList = map.addedFaceList();
1766         const List<objectMap>& asfList = slaveMap.addedFaceList();
1767
1768         forAll(amfList, mfI)
1769         {
1770             // Configure a face for comparison.
1771             const face& mF = faces_[amfList[mfI].index()];
1772
1773             forAll(mF, pointI)
1774             {
1775                 cF[pointI] = cMap.entityMap(coupleMap::POINT)[mF[pointI]];
1776             }
1777
1778             bool matched = false;
1779
1780             forAll(asfList, sfI)
1781             {
1782                 const face& sF = faces_[asfList[sfI].index()];
1783
1784                 if (triFace::compare(triFace(cF), triFace(sF)))
1785                 {
1786                     cMap.mapSlave
1787                     (
1788                         coupleMap::FACE,
1789                         amfList[mfI].index(),
1790                         asfList[sfI].index()
1791                     );
1792
1793                     cMap.mapMaster
1794                     (
1795                         coupleMap::FACE,
1796                         asfList[sfI].index(),
1797                         amfList[mfI].index()
1798                     );
1799
1800                     matched = true;
1801
1802                     break;
1803                 }
1804             }
1805
1806             if (!matched)
1807             {
1808                 Pout<< "masterFaces: " << nl
1809                     << amfList << endl;
1810
1811                 Pout<< "slaveFaces: " << nl
1812                     << asfList << endl;
1813
1814                 forAll(amfList, mfI)
1815                 {
1816                     Pout<< amfList[mfI].index() << ": "
1817                         << faces_[amfList[mfI].index()]
1818                         << nl;
1819                 }
1820
1821                 forAll(asfList, sfI)
1822                 {
1823                     Pout<< asfList[sfI].index() << ": "
1824                         << faces_[asfList[sfI].index()]
1825                         << nl;
1826                 }
1827
1828                 Pout<< endl;
1829
1830                 FatalErrorIn
1831                 (
1832                     "\n"
1833                     "const changeMap dynamicTopoFvMesh::removeEdgeFlips\n"
1834                     "(\n"
1835                     "    const label eIndex,\n"
1836                     "    const scalar minQuality,\n"
1837                     "    const labelList& vertexHull,\n"
1838                     "    PtrList<scalarListList>& Q,\n"
1839                     "    PtrList<labelListList>& K,\n"
1840                     "    PtrList<labelListList>& triangulations,\n"
1841                     "    const label checkIndex\n"
1842                     ")\n"
1843                 )
1844                     << "Failed to build coupled face maps."
1845                     << abort(FatalError);
1846             }
1847         }
1848     }
1849
1850     // Finally remove the edge
1851     removeEdge(eIndex);
1852
1853     // Update map
1854     map.removeEdge(eIndex);
1855
1856     // Increment the counter
1857     statistics_[1]++;
1858
1859     // Set the flag
1860     topoChangeFlag_ = true;
1861
1862     // Return a successful operation.
1863     return map;
1864 }
1865
1866
1867 // Extract triangulations from the programming table
1868 void dynamicTopoFvMesh::extractTriangulation
1869 (
1870     const label i,
1871     const label j,
1872     const labelListList& K,
1873     label& numTriangulations,
1874     labelListList& triangulations
1875 ) const
1876 {
1877     if ( j >= (i+2) )
1878     {
1879         label k = K[i][j];
1880
1881         // Fill in the triangulation list
1882         triangulations[0][numTriangulations] = i;
1883         triangulations[1][numTriangulations] = k;
1884         triangulations[2][numTriangulations] = j;
1885
1886         // Increment triangulation count
1887         numTriangulations++;
1888
1889         // Recursively call the function for the two sub-triangulations
1890         extractTriangulation(i,k,K,numTriangulations,triangulations);
1891         extractTriangulation(k,j,K,numTriangulations,triangulations);
1892     }
1893 }
1894
1895
1896 // Identify the 3-2 swap from the triangulation sequence
1897 //  - Use an edge-plane intersection formula
1898 label dynamicTopoFvMesh::identify32Swap
1899 (
1900     const label eIndex,
1901     const labelList& hullVertices,
1902     const labelListList& triangulations,
1903     bool output
1904 ) const
1905 {
1906     label m = hullVertices.size();
1907     const edge& edgeToCheck = edges_[eIndex];
1908
1909     // Obtain intersection point.
1910     linePointRef segment
1911     (
1912         points_[edgeToCheck.start()],
1913         points_[edgeToCheck.end()]
1914     );
1915
1916     vector intPt = vector::zero;
1917
1918     // Configure a face with triangulation
1919     for (label i = 0; i < (m-2); i++)
1920     {
1921         bool intersects =
1922         (
1923             meshOps::segmentTriFaceIntersection
1924             (
1925                 triPointRef
1926                 (
1927                     points_[hullVertices[triangulations[0][i]]],
1928                     points_[hullVertices[triangulations[1][i]]],
1929                     points_[hullVertices[triangulations[2][i]]]
1930                 ),
1931                 segment,
1932                 intPt
1933             )
1934         );
1935
1936         if (intersects)
1937         {
1938             return i;
1939         }
1940     }
1941
1942     if (debug > 1 || output)
1943     {
1944         Pout<< nl << nl << "Hull Vertices: " << nl;
1945
1946         forAll(hullVertices, vertexI)
1947         {
1948             Pout<< hullVertices[vertexI] << ": "
1949                 << points_[hullVertices[vertexI]]
1950                 << nl;
1951         }
1952
1953         InfoIn
1954         (
1955             "\n"
1956             "label dynamicTopoFvMesh::identify32Swap\n"
1957             "(\n"
1958             "    const label eIndex,\n"
1959             "    const labelList& hullVertices,\n"
1960             "    const labelListList& triangulations,\n"
1961             "    bool output\n"
1962             ") const\n"
1963         )   << nl
1964             << "Could not determine 3-2 swap triangulation." << nl
1965             << "Edge: " << edgeToCheck << nl
1966             << "Edge Points: "
1967             << segment.start() << "," << segment.end() << nl
1968             << endl;
1969     }
1970
1971     // Could not find an intersecting triangulation.
1972     //  - If this is a boundary edge, a curved surface-mesh
1973     //    was probably the reason why this failed.
1974     //  - If so, declare the nearest triangulation instead.
1975     vector eCentre = segment.centre();
1976
1977     scalarField dist(m-2, 0.0);
1978
1979     label mT = -1, ePatch = whichEdgePatch(eIndex);
1980     bool foundTriangulation = false;
1981
1982     for (label i = 0; i < (m-2); i++)
1983     {
1984         // Compute edge to face-centre distance.
1985         dist[i] =
1986         (
1987             mag
1988             (
1989                 eCentre
1990               - triPointRef
1991                 (
1992                     points_[hullVertices[triangulations[0][i]]],
1993                     points_[hullVertices[triangulations[1][i]]],
1994                     points_[hullVertices[triangulations[2][i]]]
1995                 ).centre()
1996             )
1997         );
1998     }
1999
2000     while (!foundTriangulation)
2001     {
2002         mT = findMin(dist);
2003
2004         // Check validity for boundary edges
2005         if
2006         (
2007             (ePatch > -1) &&
2008             ((triangulations[0][mT] != 0) || (triangulations[2][mT] != m-1))
2009         )
2010         {
2011             // This is a 2-3 triangulation. Try again.
2012             dist[mT] = GREAT;
2013         }
2014         else
2015         {
2016             foundTriangulation = true;
2017         }
2018     }
2019
2020     if (debug > 1 || output)
2021     {
2022         Pout<< " All distances :" << dist << nl
2023             << " Triangulation index: " << mT
2024             << endl;
2025     }
2026
2027     // Return the index
2028     return mT;
2029 }
2030
2031
2032 // Routine to check whether the triangulation at the
2033 // index lies on the boundary of the vertex ring.
2034 bool dynamicTopoFvMesh::boundaryTriangulation
2035 (
2036     const label index,
2037     label& isolatedVertex,
2038     labelListList& triangulations
2039 ) const
2040 {
2041     label first = 0, second = 0, third = 0;
2042
2043     // Count for occurrences
2044     forAll(triangulations, row)
2045     {
2046         forAll(triangulations[row], col)
2047         {
2048             if (triangulations[row][col] == triangulations[0][index])
2049             {
2050                 first++;
2051             }
2052
2053             if (triangulations[row][col] == triangulations[1][index])
2054             {
2055                 second++;
2056             }
2057
2058             if (triangulations[row][col] == triangulations[2][index])
2059             {
2060                 third++;
2061             }
2062         }
2063     }
2064
2065     if (first == 1)
2066     {
2067         isolatedVertex = triangulations[0][index];
2068         return true;
2069     }
2070
2071     if (second == 1)
2072     {
2073         isolatedVertex = triangulations[1][index];
2074         return true;
2075     }
2076
2077     if (third == 1)
2078     {
2079         isolatedVertex = triangulations[2][index];
2080         return true;
2081     }
2082
2083     // This isn't a boundary triangulation
2084     return false;
2085 }
2086
2087
2088 // Utility method to compute the minimum quality of a vertex hull
2089 scalar dynamicTopoFvMesh::computeMinQuality
2090 (
2091     const label eIndex,
2092     labelList& hullVertices
2093 ) const
2094 {
2095     scalar minQuality = GREAT;
2096
2097     // Obtain a reference to this edge
2098     const edge& edgeToCheck = edges_[eIndex];
2099     const labelList& edgeFaces = edgeFaces_[eIndex];
2100
2101     // If this entity was deleted, skip it.
2102     if (edgeFaces.empty())
2103     {
2104         return minQuality;
2105     }
2106
2107     // Ensure that edge is surrounded by triangles
2108     forAll(edgeFaces, faceI)
2109     {
2110         if (faces_[edgeFaces[faceI]].size() != 3)
2111         {
2112             return minQuality;
2113         }
2114     }
2115
2116     // Build vertexHull for this edge
2117     buildVertexHull(eIndex, hullVertices);
2118
2119     if (coupledModification_)
2120     {
2121         if (locallyCoupledEntity(eIndex))
2122         {
2123             // Compute the minimum quality of the slave edge as well.
2124             label sIndex = -1;
2125
2126             // Determine the slave index.
2127             forAll(patchCoupling_, patchI)
2128             {
2129                 if (patchCoupling_(patchI))
2130                 {
2131                     const label edgeEnum  = coupleMap::EDGE;
2132                     const coupleMap& cMap = patchCoupling_[patchI].map();
2133
2134                     if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
2135                     {
2136                         break;
2137                     }
2138                 }
2139             }
2140
2141             if (sIndex == -1)
2142             {
2143                 FatalErrorIn
2144                 (
2145                     "scalar dynamicTopoFvMesh::computeMinQuality"
2146                     "(const label eIndex, labelList& hullVertices) const"
2147                 )
2148                     << nl << "Coupled maps were improperly specified." << nl
2149                     << " Slave index not found for: " << nl
2150                     << " Edge: " << eIndex << nl
2151                     << abort(FatalError);
2152             }
2153
2154             // Build vertexHull for this edge
2155             labelList slaveVertexHull;
2156             buildVertexHull(eIndex, slaveVertexHull);
2157
2158             // Temporarily turn off coupledModification
2159             unsetCoupledModification();
2160
2161             scalar slaveQuality = computeMinQuality(sIndex, slaveVertexHull);
2162
2163             minQuality = Foam::min(slaveQuality, minQuality);
2164
2165             // Turn it back on.
2166             setCoupledModification();
2167         }
2168         else
2169         if (processorCoupledEntity(eIndex))
2170         {
2171             // Don't compute minQuality here, but do it in
2172             // fillTables instead, while calculating hullPoints
2173             return minQuality;
2174         }
2175     }
2176
2177     // Compute minQuality
2178     minQuality =
2179     (
2180         Foam::min
2181         (
2182             computeMinQuality
2183             (
2184                 edgeToCheck,
2185                 hullVertices,
2186                 points_,
2187                 (whichEdgePatch(eIndex) < 0)
2188             ),
2189             minQuality
2190         )
2191     );
2192
2193     // Ensure that the mesh is valid
2194     if (minQuality < 0.0)
2195     {
2196         // Write out faces and cells for post processing.
2197         labelHashSet iFaces, iCells, bFaces;
2198
2199         const labelList& eFaces = edgeFaces_[eIndex];
2200
2201         forAll(eFaces, faceI)
2202         {
2203             iFaces.insert(eFaces[faceI]);
2204
2205             if (!iCells.found(owner_[eFaces[faceI]]))
2206             {
2207                 iCells.insert(owner_[eFaces[faceI]]);
2208             }
2209
2210             if (!iCells.found(neighbour_[eFaces[faceI]]))
2211             {
2212                 iCells.insert(neighbour_[eFaces[faceI]]);
2213             }
2214         }
2215
2216         writeVTK(Foam::name(eIndex) + "_iCells", iCells.toc());
2217         writeVTK(Foam::name(eIndex) + "_iFaces", iFaces.toc(), 2);
2218
2219         // Write out the boundary patches (for post-processing reference)
2220         for
2221         (
2222             label faceI = nOldInternalFaces_;
2223             faceI < faces_.size();
2224             faceI++
2225         )
2226         {
2227             if (faces_[faceI].empty())
2228             {
2229                 continue;
2230             }
2231
2232             label pIndex = whichPatch(faceI);
2233
2234             if (pIndex != -1)
2235             {
2236                 bFaces.insert(faceI);
2237             }
2238         }
2239
2240         writeVTK(Foam::name(eIndex) + "_bFaces", bFaces.toc(), 2);
2241
2242         FatalErrorIn
2243         (
2244             "scalar dynamicTopoFvMesh::computeMinQuality"
2245             "(const label eIndex, labelList& hullVertices) const"
2246         )
2247             << "Encountered negative cell-quality!" << nl
2248             << "Edge: " << eIndex << ": " << edgeToCheck << nl
2249             << "vertexHull: " << hullVertices << nl
2250             << "Minimum Quality: " << minQuality
2251             << abort(FatalError);
2252     }
2253
2254     return minQuality;
2255 }
2256
2257
2258 // Compute minQuality given addressing
2259 scalar dynamicTopoFvMesh::computeMinQuality
2260 (
2261     const edge& edgeToCheck,
2262     const labelList& hullVertices,
2263     const UList<point>& points,
2264     bool closedRing
2265 ) const
2266 {
2267     scalar cQuality = 0.0;
2268     scalar minQuality = GREAT;
2269
2270     // Obtain point references
2271     const point& a = points[edgeToCheck[0]];
2272     const point& c = points[edgeToCheck[1]];
2273
2274     label start = (closedRing ? 0 : 1);
2275
2276     for (label indexJ = start; indexJ < hullVertices.size(); indexJ++)
2277     {
2278         label indexI = hullVertices.rcIndex(indexJ);
2279
2280         // Pick vertices off the list
2281         const point& b = points[hullVertices[indexI]];
2282         const point& d = points[hullVertices[indexJ]];
2283
2284         // Compute the quality
2285         cQuality = tetMetric_(a, b, c, d);
2286
2287         // Check if the quality is worse
2288         minQuality = Foam::min(cQuality, minQuality);
2289     }
2290
2291     return minQuality;
2292 }
2293
2294
2295 // Method used to perform a 2-3 swap in 3D
2296 // - Returns a changeMap with a type specifying:
2297 //     1: Swap was successful
2298 // - The index of the triangulated face in map.index()
2299 const changeMap dynamicTopoFvMesh::swap23
2300 (
2301     const label isolatedVertex,
2302     const label eIndex,
2303     const label triangulationIndex,
2304     const label numTriangulations,
2305     const labelListList& triangulations,
2306     const labelList& hullVertices,
2307     const labelList& hullFaces,
2308     const labelList& hullCells
2309 )
2310 {
2311     // A 2-3 swap performs the following operations:
2312     //      [1] Remove face: [ edgeToCheck[0] edgeToCheck[1] isolatedVertex ]
2313     //      [2] Remove two cells on either side of removed face
2314     //      [3] Add one edge
2315     //      [4] Add three new faces
2316     //      [5] Add three new cells
2317     //      Update faceEdges and edgeFaces information
2318
2319     changeMap map;
2320
2321     // Obtain a copy of the edge
2322     edge edgeToCheck = edges_[eIndex];
2323
2324     label faceForRemoval = hullFaces[isolatedVertex];
2325     label vertexForRemoval = hullVertices[isolatedVertex];
2326
2327     // Determine the two cells to be removed
2328     FixedList<label,2> cellsForRemoval;
2329     cellsForRemoval[0] = owner_[faceForRemoval];
2330     cellsForRemoval[1] = neighbour_[faceForRemoval];
2331
2332     if (debug > 1)
2333     {
2334         // Print out arguments
2335         Pout<< nl
2336             << "== Swapping 2-3 ==" << nl
2337             << "Edge: " << eIndex << ": " << edgeToCheck << endl;
2338
2339         if (debug > 2)
2340         {
2341             Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
2342             Pout<< " coupledModification: " << coupledModification_ << nl;
2343
2344             label bPatch = whichEdgePatch(eIndex);
2345
2346             const polyBoundaryMesh& boundary = boundaryMesh();
2347
2348             if (bPatch == -1)
2349             {
2350                 Pout<< " Patch: Internal" << nl;
2351             }
2352             else
2353             if (bPatch < boundary.size())
2354             {
2355                 Pout<< " Patch: " << boundary[bPatch].name() << nl;
2356             }
2357             else
2358             {
2359                 Pout<< " New patch: " << bPatch << endl;
2360             }
2361
2362             Pout<< " Ring: " << hullVertices << nl
2363                 << " Faces: " << hullFaces << nl
2364                 << " Cells: " << hullCells << nl
2365                 << " Triangulation: "
2366                 << triangulations[0][triangulationIndex] << " "
2367                 << triangulations[1][triangulationIndex] << " "
2368                 << triangulations[2][triangulationIndex] << " "
2369                 << nl
2370                 << " Isolated vertex: " << isolatedVertex << endl;
2371         }
2372
2373         if (debug > 3)
2374         {
2375             writeVTK
2376             (
2377                 Foam::name(eIndex)
2378               + '(' + Foam::name(edgeToCheck[0])
2379               + ',' + Foam::name(edgeToCheck[1]) + ')'
2380               + "_beforeSwap_"
2381               + Foam::name(numTriangulations) + "_"
2382               + Foam::name(triangulationIndex),
2383                 cellsForRemoval,
2384                 3, false, true
2385             );
2386         }
2387     }
2388
2389     // Check if this is an internal face
2390     if (cellsForRemoval[1] == -1)
2391     {
2392         // Write out for post-processing
2393         Pout<< " isolatedVertex: " << isolatedVertex << nl
2394             << " triangulations: " << triangulations << nl
2395             << " numTriangulations: " << numTriangulations << nl
2396             << " triangulationIndex: " << triangulationIndex << endl;
2397
2398         writeVTK("Edge23_" + Foam::name(eIndex), eIndex, 1);
2399         writeVTK("Cells23_" + Foam::name(eIndex), hullCells, 3);
2400
2401         // Write out identify32Swap output for diagnostics
2402         identify32Swap(eIndex, hullVertices, triangulations, true);
2403
2404         FatalErrorIn
2405         (
2406             "\n"
2407             "const changeMap dynamicTopoFvMesh::swap23\n"
2408             "(\n"
2409             "    const label isolatedVertex,\n"
2410             "    const label eIndex,\n"
2411             "    const label triangulationIndex,\n"
2412             "    const label numTriangulations,\n"
2413             "    const labelListList& triangulations,\n"
2414             "    const labelList& hullVertices,\n"
2415             "    const labelList& hullFaces,\n"
2416             "    const labelList& hullCells\n"
2417             ")\n"
2418         )
2419             << " Expected an internal face,"
2420             << " but found a boundary one instead." << nl
2421             << " Looks like identify32Swap couldn't correctly identify"
2422             << " the 2-2 swap triangulation." << nl
2423             << abort(FatalError);
2424     }
2425
2426     // Add three new cells to the end of the cell list
2427     FixedList<label,3> newCellIndex(-1);
2428     FixedList<cell, 3> newTetCell(cell(4));
2429
2430     forAll(newCellIndex, cellI)
2431     {
2432         scalar avgScale = -1.0;
2433
2434         if (edgeRefinement_)
2435         {
2436             avgScale =
2437             (
2438                 0.5 *
2439                 (
2440                     lengthScale_[cellsForRemoval[0]]
2441                   + lengthScale_[cellsForRemoval[1]]
2442                 )
2443             );
2444         }
2445
2446         // Insert the cell
2447         newCellIndex[cellI] = insertCell(newTetCell[cellI], avgScale);
2448
2449         // Add this cell to the map.
2450         map.addCell(newCellIndex[cellI]);
2451     }
2452
2453     // Obtain point-ordering for the other vertices
2454     // otherVertices[0] is the point before isolatedVertex
2455     // otherVertices[1] is the point after isolatedVertex
2456     FixedList<label,2> otherVertices;
2457
2458     if (triangulations[0][triangulationIndex] == isolatedVertex)
2459     {
2460         otherVertices[0] = hullVertices[triangulations[2][triangulationIndex]];
2461         otherVertices[1] = hullVertices[triangulations[1][triangulationIndex]];
2462     }
2463     else
2464     if (triangulations[1][triangulationIndex] == isolatedVertex)
2465     {
2466         otherVertices[0] = hullVertices[triangulations[0][triangulationIndex]];
2467         otherVertices[1] = hullVertices[triangulations[2][triangulationIndex]];
2468     }
2469     else
2470     {
2471         otherVertices[0] = hullVertices[triangulations[1][triangulationIndex]];
2472         otherVertices[1] = hullVertices[triangulations[0][triangulationIndex]];
2473     }
2474
2475     // Insert three new internal faces
2476     FixedList<label,3> newFaceIndex;
2477     face tmpTriFace(3);
2478
2479     // First face: The actual triangulation
2480     tmpTriFace[0] = otherVertices[0];
2481     tmpTriFace[1] = vertexForRemoval;
2482     tmpTriFace[2] = otherVertices[1];
2483
2484     newFaceIndex[0] =
2485     (
2486         insertFace
2487         (
2488             -1,
2489             tmpTriFace,
2490             newCellIndex[0],
2491             newCellIndex[1]
2492         )
2493     );
2494
2495     // Add this face to the map.
2496     map.addFace(newFaceIndex[0]);
2497
2498     // Note the triangulation face in index()
2499     map.index() = newFaceIndex[0];
2500
2501     // Second face: Triangle involving edgeToCheck[0]
2502     tmpTriFace[0] = otherVertices[0];
2503     tmpTriFace[1] = edgeToCheck[0];
2504     tmpTriFace[2] = otherVertices[1];
2505
2506     newFaceIndex[1] =
2507     (
2508         insertFace
2509         (
2510             -1,
2511             tmpTriFace,
2512             newCellIndex[1],
2513             newCellIndex[2]
2514         )
2515     );
2516
2517     // Add this face to the map.
2518     map.addFace(newFaceIndex[1]);
2519
2520     // Third face: Triangle involving edgeToCheck[1]
2521     tmpTriFace[0] = otherVertices[1];
2522     tmpTriFace[1] = edgeToCheck[1];
2523     tmpTriFace[2] = otherVertices[0];
2524
2525     newFaceIndex[2] =
2526     (
2527         insertFace
2528         (
2529             -1,
2530             tmpTriFace,
2531             newCellIndex[0],
2532             newCellIndex[2]
2533         )
2534     );
2535
2536     // Add this face to the map.
2537     map.addFace(newFaceIndex[2]);
2538
2539     // Append three dummy faceEdges entries.
2540     for (label i = 0; i < 3; i++)
2541     {
2542         faceEdges_.append(labelList(0));
2543     }
2544
2545     // Add an entry to edgeFaces
2546     labelList newEdgeFaces(3);
2547     newEdgeFaces[0] = newFaceIndex[0];
2548     newEdgeFaces[1] = newFaceIndex[1];
2549     newEdgeFaces[2] = newFaceIndex[2];
2550
2551     // Add a new internal edge to the mesh
2552     label newEdgeIndex =
2553     (
2554         insertEdge
2555         (
2556             -1,
2557             edge
2558             (
2559                 otherVertices[0],
2560                 otherVertices[1]
2561             ),
2562             newEdgeFaces
2563         )
2564     );
2565
2566     // Add this edge to the map.
2567     map.addEdge(newEdgeIndex);
2568
2569     // Define the six edges to check while building faceEdges:
2570     FixedList<edge,6> check;
2571
2572     check[0][0] = vertexForRemoval; check[0][1] = otherVertices[0];
2573     check[1][0] = vertexForRemoval; check[1][1] = otherVertices[1];
2574     check[2][0] = edgeToCheck[0];   check[2][1] = otherVertices[0];
2575     check[3][0] = edgeToCheck[1];   check[3][1] = otherVertices[0];
2576     check[4][0] = edgeToCheck[0];   check[4][1] = otherVertices[1];
2577     check[5][0] = edgeToCheck[1];   check[5][1] = otherVertices[1];
2578
2579     // Add three new entries to faceEdges
2580     label nE0 = 0, nE1 = 0, nE2 = 0;
2581     FixedList<labelList,3> newFaceEdges(labelList(3));
2582
2583     newFaceEdges[0][nE0++] = newEdgeIndex;
2584     newFaceEdges[1][nE1++] = newEdgeIndex;
2585     newFaceEdges[2][nE2++] = newEdgeIndex;
2586
2587     // Fill-in information for the three new cells,
2588     // and correct info on existing neighbouring cells
2589     label nF0 = 0, nF1 = 0, nF2 = 0;
2590     FixedList<bool,2> foundEdge;
2591
2592     // Add the newly created faces to cells
2593     newTetCell[0][nF0++] = newFaceIndex[0];
2594     newTetCell[0][nF0++] = newFaceIndex[2];
2595     newTetCell[1][nF1++] = newFaceIndex[0];
2596     newTetCell[1][nF1++] = newFaceIndex[1];
2597     newTetCell[2][nF2++] = newFaceIndex[1];
2598     newTetCell[2][nF2++] = newFaceIndex[2];
2599
2600     forAll(cellsForRemoval, cellI)
2601     {
2602         label cellIndex = cellsForRemoval[cellI];
2603
2604         forAll(cells_[cellIndex], faceI)
2605         {
2606             label faceIndex = cells_[cellIndex][faceI];
2607
2608             foundEdge[0] = false; foundEdge[1] = false;
2609
2610             // Check if face contains edgeToCheck[0]
2611             if
2612             (
2613                 (faces_[faceIndex][0] == edgeToCheck[0])
2614              || (faces_[faceIndex][1] == edgeToCheck[0])
2615              || (faces_[faceIndex][2] == edgeToCheck[0])
2616             )
2617             {
2618                 foundEdge[0] = true;
2619             }
2620
2621             // Check if face contains edgeToCheck[1]
2622             if
2623             (
2624                 (faces_[faceIndex][0] == edgeToCheck[1])
2625              || (faces_[faceIndex][1] == edgeToCheck[1])
2626              || (faces_[faceIndex][2] == edgeToCheck[1])
2627             )
2628             {
2629                 foundEdge[1] = true;
2630             }
2631
2632             // Face is connected to edgeToCheck[0]
2633             if (foundEdge[0] && !foundEdge[1])
2634             {
2635                 // Check if a face-flip is necessary
2636                 if (owner_[faceIndex] == cellIndex)
2637                 {
2638                     if (neighbour_[faceIndex] == -1)
2639                     {
2640                         // Change the owner
2641                         owner_[faceIndex] = newCellIndex[1];
2642                     }
2643                     else
2644                     {
2645                         // Flip this face
2646                         faces_[faceIndex] = faces_[faceIndex].reverseFace();
2647                         owner_[faceIndex] = neighbour_[faceIndex];
2648                         neighbour_[faceIndex] = newCellIndex[1];
2649
2650                         setFlip(faceIndex);
2651                     }
2652                 }
2653                 else
2654                 {
2655                     // Flip is unnecessary. Just update neighbour
2656                     neighbour_[faceIndex] = newCellIndex[1];
2657                 }
2658
2659                 // Add this face to the cell
2660                 newTetCell[1][nF1++] = faceIndex;
2661
2662                 // Update faceEdges and edgeFaces
2663                 forAll(faceEdges_[faceIndex], edgeI)
2664                 {
2665                     if (edges_[faceEdges_[faceIndex][edgeI]] == check[0])
2666                     {
2667                         newFaceEdges[0][nE0++] = faceEdges_[faceIndex][edgeI];
2668
2669                         meshOps::sizeUpList
2670                         (
2671                             newFaceIndex[0],
2672                             edgeFaces_[faceEdges_[faceIndex][edgeI]]
2673                         );
2674                     }
2675
2676                     if (edges_[faceEdges_[faceIndex][edgeI]] == check[1])
2677                     {
2678                         newFaceEdges[0][nE0++] = faceEdges_[faceIndex][edgeI];
2679
2680                         meshOps::sizeUpList
2681                         (
2682                             newFaceIndex[0],
2683                             edgeFaces_[faceEdges_[faceIndex][edgeI]]
2684                         );
2685                     }
2686
2687                     if (edges_[faceEdges_[faceIndex][edgeI]] == check[2])
2688                     {
2689                         newFaceEdges[1][nE1++] = faceEdges_[faceIndex][edgeI];
2690
2691                         meshOps::sizeUpList
2692                         (
2693                             newFaceIndex[1],
2694                             edgeFaces_[faceEdges_[faceIndex][edgeI]]
2695                         );
2696                     }
2697
2698                     if (edges_[faceEdges_[faceIndex][edgeI]] == check[4])
2699                     {
2700                         newFaceEdges[1][nE1++] = faceEdges_[faceIndex][edgeI];
2701
2702                         meshOps::sizeUpList
2703                         (
2704                             newFaceIndex[1],
2705                             edgeFaces_[faceEdges_[faceIndex][edgeI]]
2706                         );
2707                     }
2708                 }
2709             }
2710
2711             // Face is connected to edgeToCheck[1]
2712             if (!foundEdge[0] && foundEdge[1])
2713             {
2714                 // Check if a face-flip is necessary
2715                 if (owner_[faceIndex] == cellIndex)
2716                 {
2717                     if (neighbour_[faceIndex] == -1)
2718                     {
2719                         // Change the owner
2720                         owner_[faceIndex] = newCellIndex[0];
2721                     }
2722                     else
2723                     {
2724                         // Flip this face
2725                         faces_[faceIndex] = faces_[faceIndex].reverseFace();
2726                         owner_[faceIndex] = neighbour_[faceIndex];
2727                         neighbour_[faceIndex] = newCellIndex[0];
2728
2729                         setFlip(faceIndex);
2730                     }
2731                 }
2732                 else
2733                 {
2734                     // Flip is unnecessary. Just update neighbour
2735                     neighbour_[faceIndex] = newCellIndex[0];
2736                 }
2737
2738                 // Add this face to the cell
2739                 newTetCell[0][nF0++] = faceIndex;
2740
2741                 // Update faceEdges and edgeFaces
2742                 const labelList& fEdges = faceEdges_[faceIndex];
2743
2744                 forAll(fEdges, edgeI)
2745                 {
2746                     if (edges_[fEdges[edgeI]] == check[3])
2747                     {
2748                         newFaceEdges[2][nE2++] = fEdges[edgeI];
2749
2750                         meshOps::sizeUpList
2751                         (
2752                             newFaceIndex[2],
2753                             edgeFaces_[fEdges[edgeI]]
2754                         );
2755                     }
2756
2757                     if (edges_[fEdges[edgeI]] == check[5])
2758                     {
2759                         newFaceEdges[2][nE2++] = fEdges[edgeI];
2760
2761                         meshOps::sizeUpList
2762                         (
2763                             newFaceIndex[2],
2764                             edgeFaces_[fEdges[edgeI]]
2765                         );
2766                     }
2767                 }
2768             }
2769
2770             // Face is connected to both edgeToCheck [0] and [1]
2771             if
2772             (
2773                 (foundEdge[0] && foundEdge[1]) &&
2774                 (faceIndex != faceForRemoval)
2775             )
2776             {
2777                 // Check if a face-flip is necessary
2778                 if (owner_[faceIndex] == cellIndex)
2779                 {
2780                     if (neighbour_[faceIndex] == -1)
2781                     {
2782                         // Change the owner
2783                         owner_[faceIndex] = newCellIndex[2];
2784                     }
2785                     else
2786                     {
2787                         // Flip this face
2788                         faces_[faceIndex] = faces_[faceIndex].reverseFace();
2789                         owner_[faceIndex] = neighbour_[faceIndex];
2790                         neighbour_[faceIndex] = newCellIndex[2];
2791
2792                         setFlip(faceIndex);
2793                     }
2794                 }
2795                 else
2796                 {
2797                     // Flip is unnecessary. Just update neighbour
2798                     neighbour_[faceIndex] = newCellIndex[2];
2799                 }
2800
2801                 // Add this face to the cell
2802                 newTetCell[2][nF2++] = faceIndex;
2803             }
2804         }
2805     }
2806
2807     // Now update faceEdges for the three new faces
2808     forAll(newFaceEdges, faceI)
2809     {
2810         faceEdges_[newFaceIndex[faceI]] = newFaceEdges[faceI];
2811     }
2812
2813     // Update edgeFaces for edges of the removed face
2814     forAll(faceEdges_[faceForRemoval], edgeI)
2815     {
2816         label edgeIndex = faceEdges_[faceForRemoval][edgeI];
2817
2818         meshOps::sizeDownList
2819         (
2820             faceForRemoval,
2821             edgeFaces_[edgeIndex]
2822         );
2823     }
2824
2825     // Remove the face
2826     removeFace(faceForRemoval);
2827
2828     // Update map
2829     map.removeFace(faceForRemoval);
2830
2831     forAll(cellsForRemoval, cellI)
2832     {
2833         removeCell(cellsForRemoval[cellI]);
2834
2835         // Update map
2836         map.removeCell(cellsForRemoval[cellI]);
2837     }
2838
2839     // Update the cell list with newly configured cells.
2840     forAll(newCellIndex, cellI)
2841     {
2842         cells_[newCellIndex[cellI]] = newTetCell[cellI];
2843
2844         if (cellI == 2)
2845         {
2846             // Skip mapping for the intermediate cell.
2847             setCellMapping(newCellIndex[cellI], hullCells, false);
2848         }
2849         else
2850         {
2851             // Set the mapping for this cell
2852             setCellMapping(newCellIndex[cellI], hullCells);
2853         }
2854     }
2855
2856     // Fill in mapping information for three new faces.
2857     // Since they're all internal, interpolate fluxes by default.
2858     forAll(newFaceIndex, faceI)
2859     {
2860         setFaceMapping(newFaceIndex[faceI]);
2861     }
2862
2863     if (debug > 2)
2864     {
2865         Pout<< "Added edge: " << nl;
2866
2867         Pout<< newEdgeIndex << ":: "
2868             << edges_[newEdgeIndex]
2869             << " edgeFaces: "
2870             << edgeFaces_[newEdgeIndex]
2871             << nl;
2872
2873         Pout<< "Added faces: " << nl;
2874
2875         forAll(newFaceIndex, faceI)
2876         {
2877             Pout<< newFaceIndex[faceI] << ":: "
2878                 << faces_[newFaceIndex[faceI]]
2879                 << " faceEdges: "
2880                 << faceEdges_[newFaceIndex[faceI]]
2881                 << nl;
2882         }
2883
2884         Pout<< "Added cells: " << nl;
2885
2886         forAll(newCellIndex, cellI)
2887         {
2888             Pout<< newCellIndex[cellI] << ":: "
2889                 << cells_[newCellIndex[cellI]]
2890                 << nl;
2891         }
2892
2893         Pout<< endl;
2894
2895         if (debug > 3)
2896         {
2897             writeVTK
2898             (
2899                 Foam::name(eIndex)
2900               + '(' + Foam::name(edgeToCheck[0])
2901               + ',' + Foam::name(edgeToCheck[1]) + ')'
2902               + "_afterSwap_"
2903               + Foam::name(numTriangulations) + "_"
2904               + Foam::name(triangulationIndex),
2905                 newCellIndex,
2906                 3, false, true
2907             );
2908         }
2909     }
2910
2911     // Specify that the swap was successful
2912     map.type() = 1;
2913
2914     // Return the changeMap
2915     return map;
2916 }
2917
2918
2919 // Method used to perform a 2-2 / 3-2 swap in 3D
2920 // - Returns a changeMap with a type specifying:
2921 //     1: Swap was successful
2922 // - The index of the triangulated face in map.index()
2923 const changeMap dynamicTopoFvMesh::swap32
2924 (
2925     const label eIndex,
2926     const label triangulationIndex,
2927     const label numTriangulations,
2928     const labelListList& triangulations,
2929     const labelList& hullVertices,
2930     const labelList& hullFaces,
2931     const labelList& hullCells
2932 )
2933 {
2934     // A 2-2 / 3-2 swap performs the following operations:
2935     //      [1] Remove three faces surrounding edgeToCheck
2936     //      [2] Remove two (2-2 swap) or three(3-2 swap)
2937     //          cells surrounding edgeToCheck
2938     //      [3] Add one internal face
2939     //      [4] Add two new cells
2940     //      [5] If edgeToCheck is on a boundary,
2941     //          add two boundary faces and a boundary edge (2-2 swap)
2942     //      eIndex is removed later by removeEdgeFlips
2943     //      Update faceEdges and edgeFaces information
2944
2945     changeMap map;
2946
2947     // Obtain a copy of the edge
2948     edge edgeToCheck = edges_[eIndex];
2949
2950     // Determine the patch this edge belongs to
2951     label edgePatch = whichEdgePatch(eIndex);
2952
2953     // Determine the three faces to be removed
2954     FixedList<label,3> facesForRemoval;
2955     DynamicList<label> cellRemovalList(3);
2956
2957     forAll(facesForRemoval, faceI)
2958     {
2959         facesForRemoval[faceI] =
2960         (
2961             hullFaces[triangulations[faceI][triangulationIndex]]
2962         );
2963
2964         label own = owner_[facesForRemoval[faceI]];
2965         label nei = neighbour_[facesForRemoval[faceI]];
2966
2967         // Check and add cells as well
2968         if (findIndex(cellRemovalList, own) == -1)
2969         {
2970             cellRemovalList.append(own);
2971         }
2972
2973         if (nei != -1)
2974         {
2975             if (findIndex(cellRemovalList, nei) == -1)
2976             {
2977                 cellRemovalList.append(nei);
2978             }
2979         }
2980     }
2981
2982     if (debug > 1)
2983     {
2984         // Print out arguments
2985         Pout<< nl;
2986
2987         if (edgePatch < 0)
2988         {
2989             Pout<< "== Swapping 3-2 ==" << nl;
2990         }
2991         else
2992         {
2993             Pout<< "== Swapping 2-2 ==" << nl;
2994         }
2995
2996         Pout<< " Edge: " << eIndex << ": " << edgeToCheck << endl;
2997
2998         if (debug > 2)
2999         {
3000             Pout<< " On SubMesh: " << Switch::asText(isSubMesh_) << nl;
3001             Pout<< " coupledModification: " << coupledModification_ << nl;
3002
3003             label bPatch = whichEdgePatch(eIndex);
3004
3005             const polyBoundaryMesh& boundary = boundaryMesh();
3006
3007             if (bPatch == -1)
3008             {
3009                 Pout<< " Patch: Internal" << nl;
3010             }
3011             else
3012             if (bPatch < boundary.size())
3013             {
3014                 Pout<< " Patch: " << boundary[bPatch].name() << nl;
3015             }
3016             else
3017             {
3018                 Pout<< " New patch: " << bPatch << endl;
3019             }
3020
3021             Pout<< " Ring: " << hullVertices << nl
3022                 << " Faces: " << hullFaces << nl
3023                 << " Cells: " << hullCells << nl
3024                 << " Triangulation: "
3025                 << triangulations[0][triangulationIndex] << " "
3026                 << triangulations[1][triangulationIndex] << " "
3027                 << triangulations[2][triangulationIndex] << " "
3028                 << endl;
3029         }
3030
3031         if (debug > 3)
3032         {
3033             writeVTK
3034             (
3035                 Foam::name(eIndex)
3036               + '(' + Foam::name(edgeToCheck[0])
3037               + ',' + Foam::name(edgeToCheck[1]) + ')'
3038               + "_beforeSwap_"
3039               + Foam::name(numTriangulations) + "_"
3040               + Foam::name(triangulationIndex),
3041                 cellRemovalList,
3042                 3, false, true
3043             );
3044         }
3045     }
3046
3047     // Add two new cells to the end of the cell list
3048     FixedList<label,2> newCellIndex(-1);
3049     FixedList<cell, 2> newTetCell(cell(4));
3050
3051     forAll(newCellIndex, cellI)
3052     {
3053         scalar avgScale = 0.0;
3054
3055         if (edgeRefinement_)
3056         {
3057             forAll(cellRemovalList, indexI)
3058             {
3059                 avgScale += lengthScale_[cellRemovalList[indexI]];
3060             }
3061
3062             avgScale /= cellRemovalList.size();
3063         }
3064
3065         // Insert the cell
3066         newCellIndex[cellI] = insertCell(newTetCell[cellI], avgScale);
3067
3068         // Add this cell to the map.
3069         map.addCell(newCellIndex[cellI]);
3070     }
3071
3072     // Insert a new internal face
3073     face newTriFace(3);
3074
3075     newTriFace[0] = hullVertices[triangulations[0][triangulationIndex]];
3076     newTriFace[1] = hullVertices[triangulations[1][triangulationIndex]];
3077     newTriFace[2] = hullVertices[triangulations[2][triangulationIndex]];
3078
3079     label newFaceIndex =
3080     (
3081         insertFace
3082         (
3083             -1,
3084             newTriFace,
3085             newCellIndex[0],
3086             newCellIndex[1]
3087         )
3088     );
3089
3090     // Add this face to the map.
3091     map.addFace(newFaceIndex);
3092
3093     // Note the triangulation face in index()
3094     map.index() = newFaceIndex;
3095
3096     // Add faceEdges for the new face as well.
3097     faceEdges_.append(labelList(3));
3098
3099     // Define the three edges to check while building faceEdges:
3100     FixedList<edge,3> check;
3101
3102     check[0][0] = newTriFace[0]; check[0][1] = newTriFace[1];
3103     check[1][0] = newTriFace[1]; check[1][1] = newTriFace[2];
3104     check[2][0] = newTriFace[2]; check[2][1] = newTriFace[0];
3105
3106     // For 2-2 swaps, two faces are introduced
3107     label nE = 0, nBf = 0;
3108     FixedList<label,2> nBE(0);
3109     FixedList<labelList,2> bdyFaceEdges(labelList(3, -1));
3110
3111     // Fill-in information for the two new cells,
3112     // and correct info on existing neighbouring cells
3113     label nF0 = 0, nF1 = 0;
3114     label otherPoint = -1, nextPoint = -1;
3115     FixedList<bool,2> foundEdge;
3116
3117     // For a 2-2 swap on a boundary edge,
3118     // add two boundary faces and an edge
3119     label newEdgeIndex = -1;
3120     labelList oldBdyFaceIndex(2, -1), newBdyFaceIndex(2, -1);
3121
3122     if (edgePatch > -1)
3123     {
3124         // Temporary local variables
3125         label facePatch = -1;
3126         edge newEdge(-1, -1);
3127         FixedList<label,2> nBEdge(0);
3128         FixedList<FixedList<label,2>,2> bdyEdges;
3129         FixedList<face,2> newBdyTriFace(face(3));
3130
3131         // Get a cue for face orientation from existing faces
3132         forAll(facesForRemoval, faceI)
3133         {
3134             if (neighbour_[facesForRemoval[faceI]] == -1)
3135             {
3136                 facePatch = whichPatch(facesForRemoval[faceI]);
3137
3138                 // Record this face-index for mapping.
3139                 oldBdyFaceIndex[nBf++] = facesForRemoval[faceI];
3140
3141                 meshOps::findIsolatedPoint
3142                 (
3143                     faces_[facesForRemoval[faceI]],
3144                     edgeToCheck,
3145                     otherPoint,
3146                     nextPoint
3147                 );
3148
3149                 if (nextPoint == edgeToCheck[0])
3150                 {
3151                     newEdge[1] = otherPoint;
3152                     newBdyTriFace[0][0] = otherPoint;
3153                     newBdyTriFace[0][1] = edgeToCheck[0];
3154                     newBdyTriFace[1][2] = otherPoint;
3155                 }
3156                 else
3157                 {
3158                     newEdge[0] = otherPoint;
3159                     newBdyTriFace[1][0] = otherPoint;
3160                     newBdyTriFace[1][1] = edgeToCheck[1];
3161                     newBdyTriFace[0][2] = otherPoint;
3162                 }
3163
3164                 // Also update faceEdges for the new boundary faces
3165                 forAll(faceEdges_[facesForRemoval[faceI]], edgeI)
3166                 {
3167                     if
3168                     (
3169                         edges_[faceEdges_[facesForRemoval[faceI]][edgeI]]
3170                      == edge(edgeToCheck[0], otherPoint)
3171                     )
3172                     {
3173                         bdyFaceEdges[0][nBE[0]++] =
3174                         (
3175                             faceEdges_[facesForRemoval[faceI]][edgeI]
3176                         );
3177
3178                         bdyEdges[0][nBEdge[0]++] =
3179                         (
3180                             faceEdges_[facesForRemoval[faceI]][edgeI]
3181                         );
3182                     }
3183
3184                     if
3185                     (
3186                         edges_[faceEdges_[facesForRemoval[faceI]][edgeI]]
3187                      == edge(edgeToCheck[1], otherPoint)
3188                     )
3189                     {
3190                         bdyFaceEdges[1][nBE[1]++] =
3191                         (
3192                             faceEdges_[facesForRemoval[faceI]][edgeI]
3193                         );
3194
3195                         bdyEdges[1][nBEdge[1]++] =
3196                         (
3197                             faceEdges_[facesForRemoval[faceI]][edgeI]
3198                         );
3199                     }
3200                 }
3201             }
3202         }
3203
3204         // Insert the first of two new faces
3205         newBdyFaceIndex[0] =
3206         (
3207             insertFace
3208             (
3209                 facePatch,
3210                 newBdyTriFace[0],
3211                 newCellIndex[1],
3212                 -1
3213             )
3214         );
3215
3216         // Add this face to the map.
3217         map.addFace(newBdyFaceIndex[0]);
3218
3219         // Insert the second of two new faces
3220         newBdyFaceIndex[1] =
3221         (
3222             insertFace
3223             (
3224                 facePatch,
3225                 newBdyTriFace[1],
3226                 newCellIndex[0],
3227                 -1
3228             )
3229         );
3230
3231         // Add this face to the map.
3232         map.addFace(newBdyFaceIndex[1]);
3233
3234         // Update the new cells
3235         newTetCell[0][nF0++] = newBdyFaceIndex[1];
3236         newTetCell[1][nF1++] = newBdyFaceIndex[0];
3237
3238         // Add an edgeFaces entry
3239         labelList newBdyEdgeFaces(3, -1);
3240         newBdyEdgeFaces[0] = newBdyFaceIndex[0];
3241         newBdyEdgeFaces[1] = newFaceIndex;
3242         newBdyEdgeFaces[2] = newBdyFaceIndex[1];
3243
3244         // Find the point other than the new edge
3245         // on the new triangular face
3246         meshOps::findIsolatedPoint
3247         (
3248             newTriFace,
3249             newEdge,
3250             otherPoint,
3251             nextPoint
3252         );
3253
3254         // Insert the edge
3255         newEdgeIndex =
3256         (
3257             insertEdge
3258             (
3259                 edgePatch,
3260                 newEdge,
3261                 newBdyEdgeFaces
3262             )
3263         );
3264
3265         // Add this edge to the map.
3266         map.addEdge(newEdgeIndex);
3267
3268         // Update faceEdges with the new edge
3269         faceEdges_[newFaceIndex][nE++] = newEdgeIndex;
3270         bdyFaceEdges[0][nBE[0]++] = newEdgeIndex;
3271         bdyFaceEdges[1][nBE[1]++] = newEdgeIndex;
3272
3273         // Update edgeFaces with the two new faces
3274         forAll(bdyEdges[0], edgeI)
3275         {
3276             meshOps::sizeUpList
3277             (
3278                 newBdyFaceIndex[0],
3279                 edgeFaces_[bdyEdges[0][edgeI]]
3280             );
3281
3282             meshOps::sizeUpList
3283             (
3284                 newBdyFaceIndex[1],
3285                 edgeFaces_[bdyEdges[1][edgeI]]
3286             );
3287         }
3288
3289         // Add faceEdges for the two new boundary faces
3290         faceEdges_.append(bdyFaceEdges[0]);
3291         faceEdges_.append(bdyFaceEdges[1]);
3292
3293         // Update the number of surface swaps.
3294         statistics_[2]++;
3295     }
3296
3297     newTetCell[0][nF0++] = newFaceIndex;
3298     newTetCell[1][nF1++] = newFaceIndex;
3299
3300     forAll(cellRemovalList, cellI)
3301     {
3302         label cellIndex = cellRemovalList[cellI];
3303
3304         forAll(cells_[cellIndex], faceI)
3305         {
3306             label faceIndex = cells_[cellIndex][faceI];
3307
3308             foundEdge[0] = false; foundEdge[1] = false;
3309
3310             // Find the face that contains only
3311             // edgeToCheck[0] or edgeToCheck[1]
3312             forAll(faces_[faceIndex], pointI)
3313             {
3314                 if (faces_[faceIndex][pointI] == edgeToCheck[0])
3315                 {
3316                     foundEdge[0] = true;
3317                 }
3318
3319                 if (faces_[faceIndex][pointI] == edgeToCheck[1])
3320                 {
3321                     foundEdge[1] = true;
3322                 }
3323             }
3324
3325             // Face is connected to edgeToCheck[0]
3326             if (foundEdge[0] && !foundEdge[1])
3327             {
3328                 // Check if a face-flip is necessary
3329                 if (owner_[faceIndex] == cellIndex)
3330                 {
3331                     if (neighbour_[faceIndex] == -1)
3332                     {
3333                         // Change the owner
3334                         owner_[faceIndex] = newCellIndex[1];
3335                     }
3336                     else
3337                     {
3338                         // Flip this face
3339                         faces_[faceIndex] = faces_[faceIndex].reverseFace();
3340                         owner_[faceIndex] = neighbour_[faceIndex];
3341                         neighbour_[faceIndex] = newCellIndex[1];
3342
3343                         setFlip(faceIndex);
3344                     }
3345                 }
3346                 else
3347                 {
3348                     // Flip is unnecessary. Just update neighbour
3349                     neighbour_[faceIndex] = newCellIndex[1];
3350                 }
3351
3352                 // Add this face to the cell
3353                 newTetCell[1][nF1++] = faceIndex;
3354
3355                 // Update faceEdges and edgeFaces
3356                 forAll(faceEdges_[faceIndex], edgeI)
3357                 {
3358                     if
3359                     (
3360                         (edges_[faceEdges_[faceIndex][edgeI]] == check[0])
3361                      || (edges_[faceEdges_[faceIndex][edgeI]] == check[1])
3362                      || (edges_[faceEdges_[faceIndex][edgeI]] == check[2])
3363                     )
3364                     {
3365                         faceEdges_[newFaceIndex][nE++] =
3366                         (
3367                             faceEdges_[faceIndex][edgeI]
3368                         );
3369
3370                         meshOps::sizeUpList
3371                         (
3372                             newFaceIndex,
3373                             edgeFaces_[faceEdges_[faceIndex][edgeI]]
3374                         );
3375
3376                         break;
3377                     }
3378                 }
3379             }
3380
3381             // Face is connected to edgeToCheck[1]
3382             if (!foundEdge[0] && foundEdge[1])
3383             {
3384                 // Check if a face-flip is necessary
3385                 if (owner_[faceIndex] == cellIndex)
3386                 {
3387                     if (neighbour_[faceIndex] == -1)
3388                     {
3389                         // Change the owner
3390                         owner_[faceIndex] = newCellIndex[0];
3391                     }
3392                     else
3393                     {
3394                         // Flip this face
3395                         faces_[faceIndex] = faces_[faceIndex].reverseFace();
3396                         owner_[faceIndex] = neighbour_[faceIndex];
3397                         neighbour_[faceIndex] = newCellIndex[0];
3398
3399                         setFlip(faceIndex);
3400                     }
3401                 }
3402                 else
3403                 {
3404                     // Flip is unnecessary. Just update neighbour
3405                     neighbour_[faceIndex] = newCellIndex[0];
3406                 }
3407
3408                 // Add this face to the cell
3409                 newTetCell[0][nF0++] = faceIndex;
3410             }
3411         }
3412     }
3413
3414     // Remove the faces and update associated edges
3415     forAll(facesForRemoval, faceI)
3416     {
3417         // Update edgeFaces
3418         forAll(faceEdges_[facesForRemoval[faceI]], edgeI)
3419         {
3420             label edgeIndex = faceEdges_[facesForRemoval[faceI]][edgeI];
3421
3422             if (edgeIndex != eIndex)
3423             {
3424                 meshOps::sizeDownList
3425                 (
3426                     facesForRemoval[faceI],
3427                     edgeFaces_[edgeIndex]
3428                 );
3429             }
3430         }
3431
3432         // Now remove the face...
3433         removeFace(facesForRemoval[faceI]);
3434
3435         // Update map
3436         map.removeFace(facesForRemoval[faceI]);
3437     }
3438
3439     forAll(cellRemovalList, cellI)
3440     {
3441         removeCell(cellRemovalList[cellI]);
3442
3443         // Update map
3444         map.removeCell(cellRemovalList[cellI]);
3445     }
3446
3447     // Update the cell list with newly configured cells.
3448     forAll(newCellIndex, cellI)
3449     {
3450         cells_[newCellIndex[cellI]] = newTetCell[cellI];
3451
3452         // Set the mapping for this cell
3453         setCellMapping(newCellIndex[cellI], hullCells);
3454     }
3455
3456     // Set fill-in mapping for two new boundary faces
3457     if (edgePatch > -1)
3458     {
3459         forAll(newBdyFaceIndex, i)
3460         {
3461             // Set the mapping for this face
3462             setFaceMapping(newBdyFaceIndex[i], oldBdyFaceIndex);
3463         }
3464     }
3465
3466     // Fill in mapping information for the new face.
3467     // Since it is internal, interpolate fluxes by default.
3468     setFaceMapping(newFaceIndex);
3469
3470     if (debug > 2)
3471     {
3472         if (edgePatch > -1)
3473         {
3474             Pout<< "Added edge: "
3475                 << newEdgeIndex << ":: "
3476                 << edges_[newEdgeIndex]
3477                 << " edgeFaces: "
3478                 << edgeFaces_[newEdgeIndex]
3479                 << endl;
3480         }
3481
3482         Pout<< "Added face(s): " << nl;
3483
3484         Pout<< newFaceIndex << ":: "
3485             << faces_[newFaceIndex];
3486
3487         Pout<< " faceEdges: "
3488             << faceEdges_[newFaceIndex]
3489             << endl;
3490
3491         if (edgePatch > -1)
3492         {
3493             forAll(newBdyFaceIndex, faceI)
3494             {
3495                 Pout<< newBdyFaceIndex[faceI] << ":: "
3496                     << faces_[newBdyFaceIndex[faceI]]
3497                     << " faceEdges: "
3498                     << faceEdges_[newBdyFaceIndex[faceI]]
3499                     << nl;
3500             }
3501         }
3502
3503         Pout<< "Added cells: " << nl;
3504
3505         forAll(newCellIndex, cellI)
3506         {
3507             Pout<< newCellIndex[cellI] << ":: "
3508                 << cells_[newCellIndex[cellI]]
3509                 << nl;
3510         }
3511
3512         Pout<< endl;
3513
3514         if (debug > 3)
3515         {
3516             writeVTK
3517             (
3518                 Foam::name(eIndex)
3519               + '(' + Foam::name(edgeToCheck[0])
3520               + ',' + Foam::name(edgeToCheck[1]) + ')'
3521               + "_afterSwap_"
3522               + Foam::name(numTriangulations) + "_"
3523               + Foam::name(triangulationIndex),
3524                 newCellIndex,
3525                 3, false, true
3526             );
3527         }
3528     }
3529
3530     // Specify that the swap was successful
3531     map.type() = 1;
3532
3533     // Return the changeMap
3534     return map;
3535 }
3536
3537 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3538
3539 } // End namespace Foam
3540
3541 // ************************************************************************* //