src/dynamicMesh/dynamicTopoFvMesh/edgeSwap.C

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