src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "Stack.H"
  28 #include "triFace.H"
  29 #include "objectMap.H"
  30 #include "changeMap.H"
  31 #include "multiThreader.H"
  32 #include "dynamicTopoFvMesh.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  40
  41 // Method for the bisection of a quad-face in 2D
  42 // - Returns a changeMap with a type specifying:
  43 //     1: Bisection was successful
  44 //    -1: Bisection failed since max number of topo-changes was reached.
  45 //    -2: Bisection failed since resulting quality would be unacceptable.
  46 const changeMap dynamicTopoFvMesh::bisectQuadFace
  47 (
  48     const label fIndex,
  49     const changeMap& masterMap,
  50     bool checkOnly,
  51     bool forceOp
  52 )
  53 {
  54     // Quad-face bisection performs the following operations:
  55     //      [1] Add two points at the middle of the face
  56     //      [2] Create a new internal face for each bisected cell
  57     //      [3] Modify existing face and create a new half-face
  58     //      [4] Modify triangular boundary faces, and create new ones as well
  59     //      [5] Create edges for new faces
  60     //      Update faceEdges and edgeFaces information
  61
  62     // Figure out which thread this is...
  63     label tIndex = self();
  64
  65     // Prepare the changeMaps
  66     changeMap map, slaveMap;
  67
  68     if
  69     (
  70         (statistics_[0] > maxModifications_) &&
  71         (maxModifications_ > -1)
  72     )
  73     {
  74         // Reached the max allowable topo-changes.
  75         stack(tIndex).clear();
  76
  77         return map;
  78     }
  79
  80     // Check if edgeRefinements are to be avoided on patch.
  81     if (lengthEstimator().checkRefinementPatch(whichPatch(fIndex)))
  82     {
  83         return map;
  84     }
  85
  86     // Sanity check: Is the index legitimate?
  87     if (fIndex < 0)
  88     {
  89         FatalErrorIn
  90         (
  91             "\n"
  92             "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
  93             "(\n"
  94             "    const label fIndex,\n"
  95             "    const changeMap& masterMap,\n"
  96             "    bool checkOnly\n"
  97             ")\n"
  98         )
  99             << " Invalid index: " << fIndex << nl
 100             << " nFaces: " << nFaces_
 101             << abort(FatalError);
 102     }
 103
 104     bool found;
 105     label replaceFace = -1, retainFace = -1;
 106     face tmpQuadFace(4), tmpTriFace(3);
 107     labelList tmpQFEdges(4, -1), tmpTFEdges(3, -1);
 108     FixedList<label,7> newFaceIndex(-1), newEdgeIndex(-1);
 109     FixedList<edge,4> commonEdges;
 110     FixedList<label,4> cornerEdgeIndex(-1), commonEdgeIndex(-1);
 111     FixedList<label,4> commonFaceIndex(-1);
 112     FixedList<label,2> newPointIndex(-1), newCellIndex(-1);
 113     FixedList<label,4> otherEdgeIndex(-1), otherEdgePoint(-1);
 114     FixedList<label,4> otherPointIndex(-1), nextToOtherPoint(-1);
 115     FixedList<label,2> c0BdyIndex, c0IntIndex, c1BdyIndex, c1IntIndex;
 116     FixedList<face,2>  c0BdyFace,  c0IntFace,  c1BdyFace,  c1IntFace;
 117
 118     // Get the two cells on either side...
 119     label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
 120
 121     // Keep track of old / new cells
 122     FixedList<cell, 2> oldCells(cell(5));
 123     FixedList<cell, 2> newCells(cell(5));
 124
 125     // Find the prism faces for cell[0].
 126     oldCells[0] = cells_[c0];
 127
 128     meshOps::findPrismFaces
 129     (
 130         fIndex,
 131         c0,
 132         faces_,
 133         cells_,
 134         neighbour_,
 135         c0BdyFace,
 136         c0BdyIndex,
 137         c0IntFace,
 138         c0IntIndex
 139     );
 140
 141     // Check for resulting quality
 142     if (checkBisection(fIndex, c0BdyIndex[0], forceOp))
 143     {
 144         map.type() = -2;
 145         return map;
 146     }
 147
 148     if (c1 != -1)
 149     {
 150         // Find the prism faces for cell[1].
 151         meshOps::findPrismFaces
 152         (
 153             fIndex,
 154             c1,
 155             faces_,
 156             cells_,
 157             neighbour_,
 158             c1BdyFace,
 159             c1BdyIndex,
 160             c1IntFace,
 161             c1IntIndex
 162         );
 163
 164         // Check for resulting quality
 165         if (checkBisection(fIndex, c1BdyIndex[0], forceOp))
 166         {
 167             map.type() = -2;
 168             return map;
 169         }
 170     }
 171
 172     if (debug > 1)
 173     {
 174         Info << nl << nl << "Face: " << fIndex
 175              << ": " << faces_[fIndex] << " is to be bisected. " << endl;
 176
 177         label epIndex = whichPatch(fIndex);
 178
 179         Info << "Patch: ";
 180
 181         if (epIndex == -1)
 182         {
 183             Info << "Internal" << endl;
 184         }
 185         else
 186         {
 187             Info << boundaryMesh()[epIndex].name() << endl;
 188         }
 189
 190         if (debug > 2)
 191         {
 192             Info << "Cell[0]: " << c0 << ": " << oldCells[0] << endl;
 193
 194             forAll(oldCells[0], faceI)
 195             {
 196                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
 197
 198                 Info << oldCells[0][faceI] << ": "
 199                      << faces_[oldCells[0][faceI]]
 200                      << " fE: " << fE
 201                      << endl;
 202
 203                 forAll(fE, edgeI)
 204                 {
 205                     const labelList& eF = edgeFaces_[fE[edgeI]];
 206
 207                     Info << '\t' << fE[edgeI]
 208                          << ": " << edges_[fE[edgeI]]
 209                          << " eF: " << eF
 210                          << endl;
 211                 }
 212             }
 213         }
 214
 215         // Write out VTK files prior to change
 216         if (debug > 3)
 217         {
 218             labelList cellHull(2, -1);
 219
 220             cellHull[0] = owner_[fIndex];
 221             cellHull[1] = neighbour_[fIndex];
 222
 223             writeVTK
 224             (
 225                 Foam::name(fIndex)
 226               + "_Bisect_0",
 227                 cellHull
 228             );
 229         }
 230     }
 231
 232     // Find the common-edge between the triangular boundary faces
 233     // and the face under consideration.
 234     meshOps::findCommonEdge
 235     (
 236         c0BdyIndex[0],
 237         fIndex,
 238         faceEdges_,
 239         commonEdgeIndex[0]
 240     );
 241
 242     meshOps::findCommonEdge
 243     (
 244         c0BdyIndex[1],
 245         fIndex,
 246         faceEdges_,
 247         commonEdgeIndex[1]
 248     );
 249
 250     commonEdges[0] = edges_[commonEdgeIndex[0]];
 251     commonEdges[1] = edges_[commonEdgeIndex[1]];
 252
 253     // Are we performing only checks?
 254     if (checkOnly)
 255     {
 256         map.type() = 1;
 257         return map;
 258     }
 259
 260     // Find the isolated point on both boundary faces of cell[0]
 261     meshOps::findIsolatedPoint
 262     (
 263         c0BdyFace[0],
 264         commonEdges[0],
 265         otherPointIndex[0],
 266         nextToOtherPoint[0]
 267     );
 268
 269     meshOps::findIsolatedPoint
 270     (
 271         c0BdyFace[1],
 272         commonEdges[1],
 273         otherPointIndex[1],
 274         nextToOtherPoint[1]
 275     );
 276
 277     // For convenience...
 278     otherEdgePoint[0] = commonEdges[0].otherVertex(nextToOtherPoint[0]);
 279     otherEdgePoint[1] = commonEdges[1].otherVertex(nextToOtherPoint[1]);
 280
 281     labelList mP(2, -1);
 282
 283     // Set mapping for this point
 284     mP[0] = commonEdges[0][0];
 285     mP[1] = commonEdges[0][1];
 286
 287     // Add two new points to the end of the list
 288     newPointIndex[0] =
 289     (
 290         insertPoint
 291         (
 292             0.5 * (points_[mP[0]] + points_[mP[1]]),
 293             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
 294             mP
 295         )
 296     );
 297
 298     // Set mapping for this point
 299     mP[0] = commonEdges[1][0];
 300     mP[1] = commonEdges[1][1];
 301
 302     newPointIndex[1] =
 303     (
 304         insertPoint
 305         (
 306             0.5 * (points_[mP[0]] + points_[mP[1]]),
 307             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
 308             mP
 309         )
 310     );
 311
 312     {
 313         map.addPoint(newPointIndex[0]);
 314         map.addPoint(newPointIndex[1]);
 315     }
 316
 317     // Add a new prism cell to the end of the list.
 318     // Currently invalid, but will be updated later.
 319     newCellIndex[0] = insertCell(newCells[0], lengthScale_[c0]);
 320
 321     // Modify the two existing triangle boundary faces
 322
 323     // Zeroth boundary face - Owner = c[0] & Neighbour [-1] (unchanged)
 324     meshOps::replaceLabel
 325     (
 326         otherEdgePoint[0],
 327         newPointIndex[0],
 328         c0BdyFace[0]
 329     );
 330
 331     // First boundary face - Owner = newCell[0], Neighbour = -1
 332     meshOps::replaceLabel
 333     (
 334         otherEdgePoint[1],
 335         newPointIndex[1],
 336         c0BdyFace[1]
 337     );
 338
 339     // Update faces.
 340     faces_[c0BdyIndex[0]] = c0BdyFace[0];
 341     faces_[c0BdyIndex[1]] = c0BdyFace[1];
 342
 343     owner_[c0BdyIndex[1]] = newCellIndex[0];
 344
 345     meshOps::replaceLabel
 346     (
 347         c0BdyIndex[1],
 348         -1,
 349         oldCells[0]
 350     );
 351
 352     // Detect edges other than commonEdges
 353     const labelList& fEdges = faceEdges_[fIndex];
 354
 355     forAll(fEdges, edgeI)
 356     {
 357         if
 358         (
 359             fEdges[edgeI] != commonEdgeIndex[0] &&
 360             fEdges[edgeI] != commonEdgeIndex[1]
 361         )
 362         {
 363             // Obtain a reference to this edge
 364             const edge& eThis = edges_[fEdges[edgeI]];
 365
 366             if
 367             (
 368                 eThis[0] == nextToOtherPoint[0]
 369              || eThis[1] == nextToOtherPoint[0]
 370             )
 371             {
 372                 otherEdgeIndex[0] = fEdges[edgeI];
 373             }
 374             else
 375             {
 376                 otherEdgeIndex[1] = fEdges[edgeI];
 377             }
 378         }
 379     }
 380
 381     // Modify point-labels on the quad face under consideration
 382     meshOps::replaceLabel
 383     (
 384         otherEdgePoint[0],
 385         newPointIndex[0],
 386         faces_[fIndex]
 387     );
 388
 389     meshOps::replaceLabel
 390     (
 391         nextToOtherPoint[1],
 392         newPointIndex[1],
 393         faces_[fIndex]
 394     );
 395
 396     // Add this face to the map.
 397     // Although this face isn't technically 'added', it's
 398     // required for coupled patch mapping.
 399     map.addFace(fIndex);
 400
 401     if (debug > 1)
 402     {
 403         Info << "Modified face: " << fIndex
 404              << ": " << faces_[fIndex] << endl;
 405
 406         if (debug > 2)
 407         {
 408             Info << "Common edges: " << nl
 409                  << commonEdgeIndex[0] << ": " << commonEdges[0] << nl
 410                  << commonEdgeIndex[1] << ": " << commonEdges[1]
 411                  << endl;
 412         }
 413     }
 414
 415     // Find the quad face that contains otherEdgeIndex[1]
 416     found = false;
 417
 418     const labelList& e1 = faceEdges_[c0IntIndex[0]];
 419
 420     forAll(e1, edgeI)
 421     {
 422         if (e1[edgeI] == otherEdgeIndex[1])
 423         {
 424             meshOps::replaceLabel
 425             (
 426                 c0IntIndex[0],
 427                 -1,
 428                 oldCells[0]
 429             );
 430
 431             replaceFace = c0IntIndex[0];
 432             retainFace = c0IntIndex[1];
 433             found = true;
 434             break;
 435         }
 436     }
 437
 438     if (!found)
 439     {
 440         // The edge was not found before
 441         meshOps::replaceLabel
 442         (
 443             c0IntIndex[1],
 444             -1,
 445             oldCells[0]
 446         );
 447
 448         replaceFace = c0IntIndex[1];
 449         retainFace = c0IntIndex[0];
 450     }
 451
 452     // Check if face reversal is necessary for the replacement
 453     if (owner_[replaceFace] == c0)
 454     {
 455         if (neighbour_[replaceFace] == -1)
 456         {
 457             // Change the owner
 458             owner_[replaceFace] = newCellIndex[0];
 459         }
 460         else
 461         {
 462             // This face has to be reversed
 463             faces_[replaceFace] = faces_[replaceFace].reverseFace();
 464             owner_[replaceFace] = neighbour_[replaceFace];
 465             neighbour_[replaceFace] = newCellIndex[0];
 466
 467             setFlip(replaceFace);
 468         }
 469     }
 470     else
 471     {
 472         // Keep owner, but change neighbour
 473         neighbour_[replaceFace] = newCellIndex[0];
 474     }
 475
 476     // Define the faces for the new cell
 477     newCells[0][0] = c0BdyIndex[1];
 478     newCells[0][1] = replaceFace;
 479
 480     // Define the set of new faces and insert...
 481
 482     // New interior face; Owner = cell[0] & Neighbour = newCell[0]
 483     tmpQuadFace[0] = otherPointIndex[0];
 484     tmpQuadFace[1] = newPointIndex[0];
 485     tmpQuadFace[2] = newPointIndex[1];
 486     tmpQuadFace[3] = otherPointIndex[1];
 487
 488     newFaceIndex[0] =
 489     (
 490         insertFace
 491         (
 492             -1,
 493             tmpQuadFace,
 494             c0,
 495             newCellIndex[0]
 496         )
 497     );
 498
 499     // Add a faceEdges entry as well
 500     faceEdges_.append(tmpQFEdges);
 501
 502     // Find the common edge between quad/quad faces...
 503     meshOps::findCommonEdge
 504     (
 505         c0IntIndex[0],
 506         c0IntIndex[1],
 507         faceEdges_,
 508         otherEdgeIndex[2]
 509     );
 510
 511     // ... and size-up edgeFaces for the edge
 512     meshOps::sizeUpList
 513     (
 514         newFaceIndex[0],
 515         edgeFaces_[otherEdgeIndex[2]]
 516     );
 517
 518     meshOps::replaceLabel
 519     (
 520         -1,
 521         newFaceIndex[0],
 522         newCells[0]
 523     );
 524
 525     meshOps::replaceLabel
 526     (
 527         -1,
 528         newFaceIndex[0],
 529         oldCells[0]
 530     );
 531
 532     // remove2DSliver requires this face index for removal
 533     map.addFace(newFaceIndex[0]);
 534
 535     // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
 536     tmpTriFace[0] = otherPointIndex[0];
 537     tmpTriFace[1] = newPointIndex[0];
 538     tmpTriFace[2] = otherEdgePoint[0];
 539
 540     newFaceIndex[1] =
 541     (
 542         insertFace
 543         (
 544             whichPatch(c0BdyIndex[0]),
 545             tmpTriFace,
 546             newCellIndex[0],
 547             -1
 548         )
 549     );
 550
 551     // Add a faceEdges entry as well
 552     faceEdges_.append(tmpTFEdges);
 553
 554     meshOps::replaceLabel
 555     (
 556         -1,
 557         newFaceIndex[1],
 558         newCells[0]
 559     );
 560
 561     // Third boundary face; Owner = c[0] & Neighbour = [-1]
 562     tmpTriFace[0] = otherPointIndex[1];
 563     tmpTriFace[1] = newPointIndex[1];
 564     tmpTriFace[2] = otherEdgePoint[1];
 565
 566     newFaceIndex[2] =
 567     (
 568         insertFace
 569         (
 570             whichPatch(c0BdyIndex[1]),
 571             tmpTriFace,
 572             c0,
 573             -1
 574         )
 575     );
 576
 577     // Add a faceEdges entry as well
 578     faceEdges_.append(tmpTFEdges);
 579
 580     meshOps::replaceLabel
 581     (
 582         -1,
 583         newFaceIndex[2],
 584         oldCells[0]
 585     );
 586
 587     // Create / modify edges...
 588     labelList tmpTriEdgeFaces(3, -1);
 589
 590     // The edge bisecting the zeroth boundary triangular face
 591     tmpTriEdgeFaces[0] = c0BdyIndex[0];
 592     tmpTriEdgeFaces[1] = newFaceIndex[0];
 593     tmpTriEdgeFaces[2] = newFaceIndex[1];
 594
 595     newEdgeIndex[1] =
 596     (
 597         insertEdge
 598         (
 599             whichEdgePatch(commonEdgeIndex[0]),
 600             edge(newPointIndex[0], otherPointIndex[0]),
 601             tmpTriEdgeFaces
 602         )
 603     );
 604
 605     // Find the common edge between the quad/tri faces...
 606     meshOps::findCommonEdge
 607     (
 608         c0BdyIndex[0],
 609         replaceFace,
 610         faceEdges_,
 611         cornerEdgeIndex[0]
 612     );
 613
 614     // ...and correct faceEdges / edgeFaces
 615     meshOps::replaceLabel
 616     (
 617         cornerEdgeIndex[0],
 618         newEdgeIndex[1],
 619         faceEdges_[c0BdyIndex[0]]
 620     );
 621
 622     meshOps::replaceLabel
 623     (
 624         c0BdyIndex[0],
 625         newFaceIndex[1],
 626         edgeFaces_[cornerEdgeIndex[0]]
 627     );
 628
 629     // The edge bisecting the first boundary triangular face
 630     tmpTriEdgeFaces[0] = c0BdyIndex[1];
 631     tmpTriEdgeFaces[1] = newFaceIndex[0];
 632     tmpTriEdgeFaces[2] = newFaceIndex[2];
 633
 634     newEdgeIndex[2] =
 635     (
 636         insertEdge
 637         (
 638             whichEdgePatch(commonEdgeIndex[1]),
 639             edge(newPointIndex[1], otherPointIndex[1]),
 640             tmpTriEdgeFaces
 641         )
 642     );
 643
 644     // Find the common edge between the quad/tri faces...
 645     meshOps::findCommonEdge
 646     (
 647         c0BdyIndex[1],
 648         retainFace,
 649         faceEdges_,
 650         cornerEdgeIndex[1]
 651     );
 652
 653     // ...and correct faceEdges / edgeFaces
 654     meshOps::replaceLabel
 655     (
 656         cornerEdgeIndex[1],
 657         newEdgeIndex[2],
 658         faceEdges_[c0BdyIndex[1]]
 659     );
 660
 661     meshOps::replaceLabel
 662     (
 663         c0BdyIndex[1],
 664         newFaceIndex[2],
 665         edgeFaces_[cornerEdgeIndex[1]]
 666     );
 667
 668     if (c1 == -1)
 669     {
 670         // The quad boundary face resulting from bisection;
 671         // Owner = newCell[0] & Neighbour = [-1]
 672         tmpQuadFace[0] = newPointIndex[1];
 673         tmpQuadFace[1] = nextToOtherPoint[1];
 674         tmpQuadFace[2] = otherEdgePoint[0];
 675         tmpQuadFace[3] = newPointIndex[0];
 676
 677         newFaceIndex[3] =
 678         (
 679             insertFace
 680             (
 681                 whichPatch(fIndex),
 682                 tmpQuadFace,
 683                 newCellIndex[0],
 684                 -1
 685             )
 686         );
 687
 688         // Add this face to the map.
 689         map.addFace(newFaceIndex[3]);
 690
 691         // Add a faceEdges entry as well
 692         faceEdges_.append(tmpQFEdges);
 693
 694         // Correct edgeFaces for otherEdgeIndex[1]
 695         meshOps::replaceLabel
 696         (
 697             fIndex,
 698             newFaceIndex[3],
 699             edgeFaces_[otherEdgeIndex[1]]
 700         );
 701
 702         meshOps::replaceLabel
 703         (
 704             -1,
 705             newFaceIndex[3],
 706             newCells[0]
 707         );
 708
 709         labelList tmpBiEdgeFaces(2, -1);
 710
 711         // The edge bisecting the face
 712         tmpTriEdgeFaces[0] = newFaceIndex[3];
 713         tmpTriEdgeFaces[1] = newFaceIndex[0];
 714         tmpTriEdgeFaces[2] = fIndex;
 715
 716         newEdgeIndex[0] =
 717         (
 718             insertEdge
 719             (
 720                 whichEdgePatch(otherEdgeIndex[0]),
 721                 edge(newPointIndex[0], newPointIndex[1]),
 722                 tmpTriEdgeFaces
 723             )
 724         );
 725
 726         // Replace an edge on the bisected face
 727         meshOps::replaceLabel
 728         (
 729             otherEdgeIndex[1],
 730             newEdgeIndex[0],
 731             faceEdges_[fIndex]
 732         );
 733
 734         // Create / replace side edges created from face bisection
 735         tmpBiEdgeFaces[0] = newFaceIndex[1];
 736         tmpBiEdgeFaces[1] = newFaceIndex[3];
 737
 738         newEdgeIndex[3] =
 739         (
 740             insertEdge
 741             (
 742                 whichEdgePatch(commonEdgeIndex[0]),
 743                 edge(newPointIndex[0], otherEdgePoint[0]),
 744                 tmpBiEdgeFaces
 745             )
 746         );
 747
 748         tmpBiEdgeFaces[0] = c0BdyIndex[1];
 749         tmpBiEdgeFaces[1] = newFaceIndex[3];
 750
 751         newEdgeIndex[4] =
 752         (
 753             insertEdge
 754             (
 755                 whichEdgePatch(commonEdgeIndex[1]),
 756                 edge(newPointIndex[1], nextToOtherPoint[1]),
 757                 tmpBiEdgeFaces
 758             )
 759         );
 760
 761         // Now that edges are defined, configure faceEdges
 762         // for all new faces
 763
 764         // The quad interior face; Owner = cell[0] & Neighbour = newCell[0]
 765         tmpQFEdges[0] = newEdgeIndex[0];
 766         tmpQFEdges[1] = newEdgeIndex[1];
 767         tmpQFEdges[2] = otherEdgeIndex[2];
 768         tmpQFEdges[3] = newEdgeIndex[2];
 769         faceEdges_[newFaceIndex[0]] = tmpQFEdges;
 770
 771         // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
 772         tmpTFEdges[0] = newEdgeIndex[3];
 773         tmpTFEdges[1] = cornerEdgeIndex[0];
 774         tmpTFEdges[2] = newEdgeIndex[1];
 775         faceEdges_[newFaceIndex[1]] = tmpTFEdges;
 776
 777         // Third boundary face; Owner = c[0] & Neighbour = [-1]
 778         tmpTFEdges[0] = newEdgeIndex[2];
 779         tmpTFEdges[1] = cornerEdgeIndex[1];
 780         tmpTFEdges[2] = commonEdgeIndex[1];
 781         faceEdges_[newFaceIndex[2]] = tmpTFEdges;
 782
 783         // The quad face from bisection:
 784         tmpQFEdges[0] = otherEdgeIndex[1];
 785         tmpQFEdges[1] = newEdgeIndex[3];
 786         tmpQFEdges[2] = newEdgeIndex[0];
 787         tmpQFEdges[3] = newEdgeIndex[4];
 788         faceEdges_[newFaceIndex[3]] = tmpQFEdges;
 789
 790         meshOps::replaceLabel
 791         (
 792             commonEdgeIndex[1],
 793             newEdgeIndex[4],
 794             faceEdges_[c0BdyIndex[1]]
 795         );
 796
 797         meshOps::replaceLabel
 798         (
 799             c0BdyIndex[1],
 800             newFaceIndex[2],
 801             edgeFaces_[commonEdgeIndex[1]]
 802         );
 803
 804         if (debug > 2)
 805         {
 806             Info << "Modified Cell[0]: "
 807                  << c0 << ": " << oldCells[0] << endl;
 808
 809             forAll(oldCells[0], faceI)
 810             {
 811                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
 812
 813                 Info << oldCells[0][faceI]
 814                      << ": " << faces_[oldCells[0][faceI]]
 815                      << " fE: " << fE
 816                      << endl;
 817
 818                 forAll(fE, edgeI)
 819                 {
 820                     const labelList& eF = edgeFaces_[fE[edgeI]];
 821
 822                     Info << '\t' << fE[edgeI]
 823                          << ": " << edges_[fE[edgeI]]
 824                          << " eF: " << eF
 825                          << endl;
 826                 }
 827             }
 828
 829             Info << "New Cell[0]: " << newCellIndex[0]
 830                  << ": " << newCells[0] << endl;
 831
 832             forAll(newCells[0], faceI)
 833             {
 834                 const labelList& fE = faceEdges_[newCells[0][faceI]];
 835
 836                 Info << newCells[0][faceI]
 837                      << ": " << faces_[newCells[0][faceI]]
 838                      << " fE: " << fE
 839                      << endl;
 840
 841                 forAll(fE, edgeI)
 842                 {
 843                     const labelList& eF = edgeFaces_[fE[edgeI]];
 844
 845                     Info << '\t' << fE[edgeI]
 846                          << ": " << edges_[fE[edgeI]]
 847                          << " eF: " << eF
 848                          << endl;
 849                 }
 850             }
 851         }
 852     }
 853     else
 854     {
 855         oldCells[1] = cells_[c1];
 856
 857         newCellIndex[1] = insertCell(newCells[1], lengthScale_[c1]);
 858
 859         if (debug > 2)
 860         {
 861             Info << "Cell[1]: " << c1 << ": " << oldCells[1] << endl;
 862
 863             forAll(oldCells[1], faceI)
 864             {
 865                 const labelList& fE = faceEdges_[oldCells[1][faceI]];
 866
 867                 Info << oldCells[1][faceI] << ": "
 868                      << faces_[oldCells[1][faceI]]
 869                      << " fE: " << fE
 870                      << endl;
 871
 872                 forAll(fE, edgeI)
 873                 {
 874                     const labelList& eF = edgeFaces_[fE[edgeI]];
 875
 876                     Info << '\t' << fE[edgeI]
 877                          << ": " << edges_[fE[edgeI]]
 878                          << " eF: " << eF
 879                          << endl;
 880                 }
 881             }
 882         }
 883
 884         // Find the interior face that contains otherEdgeIndex[1]
 885         found = false;
 886
 887         const labelList& e2 = faceEdges_[c1IntIndex[0]];
 888
 889         forAll(e2, edgeI)
 890         {
 891             if (e2[edgeI] == otherEdgeIndex[1])
 892             {
 893                 meshOps::replaceLabel
 894                 (
 895                     c1IntIndex[0],
 896                     -1,
 897                     oldCells[1]
 898                 );
 899
 900                 replaceFace = c1IntIndex[0];
 901                 retainFace = c1IntIndex[1];
 902                 found = true;
 903                 break;
 904             }
 905         }
 906
 907         if (!found)
 908         {
 909             // The edge was not found before
 910             meshOps::replaceLabel
 911             (
 912                 c1IntIndex[1],
 913                 -1,
 914                 oldCells[1]
 915             );
 916
 917             replaceFace = c1IntIndex[1];
 918             retainFace = c1IntIndex[0];
 919         }
 920
 921         // Check if face reversal is necessary for the replacement
 922         if (owner_[replaceFace] == c1)
 923         {
 924             if (neighbour_[replaceFace] == -1)
 925             {
 926                 // Change the owner
 927                 owner_[replaceFace] = newCellIndex[1];
 928             }
 929             else
 930             {
 931                 // This face has to be reversed
 932                 faces_[replaceFace] = faces_[replaceFace].reverseFace();
 933                 owner_[replaceFace] = neighbour_[replaceFace];
 934                 neighbour_[replaceFace] = newCellIndex[1];
 935
 936                 setFlip(replaceFace);
 937             }
 938         }
 939         else
 940         {
 941             // Keep owner, but change neighbour
 942             neighbour_[replaceFace] = newCellIndex[1];
 943         }
 944
 945         // Define attributes for the new prism cell
 946         newCells[1][0] = replaceFace;
 947
 948         // The quad interior face resulting from bisection;
 949         // Owner = newCell[0] & Neighbour = newCell[1]
 950         tmpQuadFace[0] = newPointIndex[1];
 951         tmpQuadFace[1] = nextToOtherPoint[1];
 952         tmpQuadFace[2] = otherEdgePoint[0];
 953         tmpQuadFace[3] = newPointIndex[0];
 954
 955         newFaceIndex[3] =
 956         (
 957             insertFace
 958             (
 959                 -1,
 960                 tmpQuadFace,
 961                 newCellIndex[0],
 962                 newCellIndex[1]
 963             )
 964         );
 965
 966         // Add a faceEdges entry as well
 967         faceEdges_.append(tmpQFEdges);
 968
 969         // Correct edgeFaces for otherEdgeIndex[1]
 970         meshOps::replaceLabel
 971         (
 972             fIndex,
 973             newFaceIndex[3],
 974             edgeFaces_[otherEdgeIndex[1]]
 975         );
 976
 977         meshOps::replaceLabel
 978         (
 979             -1,
 980             newFaceIndex[3],
 981             newCells[0]
 982         );
 983
 984         meshOps::replaceLabel
 985         (
 986             -1,
 987             newFaceIndex[3],
 988             newCells[1]
 989         );
 990
 991         newCells[1][1] = newFaceIndex[3];
 992
 993         // Check for common edges among the two boundary faces
 994         // Find the isolated point on both boundary faces of cell[1]
 995         if
 996         (
 997             meshOps::findCommonEdge
 998             (
 999                 c1BdyIndex[0],
1000                 c0BdyIndex[0],
1001                 faceEdges_,
1002                 commonEdgeIndex[2]
1003             )
1004         )
1005         {
1006             meshOps::findCommonEdge
1007             (
1008                 c1BdyIndex[1],
1009                 c0BdyIndex[1],
1010                 faceEdges_,
1011                 commonEdgeIndex[3]
1012             );
1013
1014             commonFaceIndex[2] = c1BdyIndex[0];
1015             commonFaceIndex[3] = c1BdyIndex[1];
1016         }
1017         else
1018         {
1019             meshOps::findCommonEdge
1020             (
1021                 c1BdyIndex[0],
1022                 c0BdyIndex[1],
1023                 faceEdges_,
1024                 commonEdgeIndex[3]
1025             );
1026
1027             meshOps::findCommonEdge
1028             (
1029                 c1BdyIndex[1],
1030                 c0BdyIndex[0],
1031                 faceEdges_,
1032                 commonEdgeIndex[2]
1033             );
1034
1035             commonFaceIndex[2] = c1BdyIndex[1];
1036             commonFaceIndex[3] = c1BdyIndex[0];
1037         }
1038
1039         commonEdges[2] = edges_[commonEdgeIndex[2]];
1040         commonEdges[3] = edges_[commonEdgeIndex[3]];
1041
1042         if (debug > 2)
1043         {
1044             Info << "Common edges: " << nl
1045                  << commonEdgeIndex[2] << ": " << commonEdges[2] << nl
1046                  << commonEdgeIndex[3] << ": " << commonEdges[3]
1047                  << endl;
1048         }
1049
1050         meshOps::findIsolatedPoint
1051         (
1052             faces_[commonFaceIndex[2]],
1053             commonEdges[2],
1054             otherPointIndex[2],
1055             nextToOtherPoint[2]
1056         );
1057
1058         meshOps::findIsolatedPoint
1059         (
1060             faces_[commonFaceIndex[3]],
1061             commonEdges[3],
1062             otherPointIndex[3],
1063             nextToOtherPoint[3]
1064         );
1065
1066         // For convenience...
1067         otherEdgePoint[2] = commonEdges[2].otherVertex(nextToOtherPoint[2]);
1068         otherEdgePoint[3] = commonEdges[3].otherVertex(nextToOtherPoint[3]);
1069
1070         // Modify the two existing triangle boundary faces
1071
1072         // Zeroth boundary face - Owner = newCell[1], Neighbour = -1
1073         meshOps::replaceLabel
1074         (
1075             otherEdgePoint[2],
1076             newPointIndex[0],
1077             faces_[commonFaceIndex[2]]
1078         );
1079
1080         owner_[commonFaceIndex[2]] = newCellIndex[1];
1081
1082         meshOps::replaceLabel
1083         (
1084             commonFaceIndex[2],
1085             -1,
1086             oldCells[1]
1087         );
1088
1089         newCells[1][2] = commonFaceIndex[2];
1090
1091         // First boundary face - Owner = c[1] & Neighbour [-1] (unchanged)
1092         meshOps::replaceLabel
1093         (
1094             otherEdgePoint[3],
1095             newPointIndex[1],
1096             faces_[commonFaceIndex[3]]
1097         );
1098
1099         // New interior face; Owner = cell[1] & Neighbour = newCell[1]
1100         tmpQuadFace[0] = newPointIndex[0];
1101         tmpQuadFace[1] = otherPointIndex[2];
1102         tmpQuadFace[2] = otherPointIndex[3];
1103         tmpQuadFace[3] = newPointIndex[1];
1104
1105         newFaceIndex[4] =
1106         (
1107             insertFace
1108             (
1109                 -1,
1110                 tmpQuadFace,
1111                 c1,
1112                 newCellIndex[1]
1113             )
1114         );
1115
1116         // Add a faceEdges entry as well
1117         faceEdges_.append(tmpQFEdges);
1118
1119         // remove2DSliver requires this face index for removal
1120         map.addFace(newFaceIndex[4]);
1121
1122         // Find the common edge between quad/quad faces...
1123         meshOps::findCommonEdge
1124         (
1125             c1IntIndex[0],
1126             c1IntIndex[1],
1127             faceEdges_,
1128             otherEdgeIndex[3]
1129         );
1130
1131         // ... and size-up edgeFaces for the edge
1132         meshOps::sizeUpList
1133         (
1134             newFaceIndex[4],
1135             edgeFaces_[otherEdgeIndex[3]]
1136         );
1137
1138         meshOps::replaceLabel
1139         (
1140             -1,
1141             newFaceIndex[4],
1142             newCells[1]
1143         );
1144
1145         meshOps::replaceLabel
1146         (
1147             -1,
1148             newFaceIndex[4],
1149             oldCells[1]
1150         );
1151
1152         // Second boundary face; Owner = cell[1] & Neighbour [-1]
1153         tmpTriFace[0] = otherPointIndex[2];
1154         tmpTriFace[1] = newPointIndex[0];
1155         tmpTriFace[2] = otherEdgePoint[2];
1156
1157         newFaceIndex[5] =
1158         (
1159             insertFace
1160             (
1161                 whichPatch(commonFaceIndex[2]),
1162                 tmpTriFace,
1163                 c1,
1164                 -1
1165             )
1166         );
1167
1168         // Add a faceEdges entry as well
1169         faceEdges_.append(tmpTFEdges);
1170
1171         meshOps::replaceLabel
1172         (
1173             -1,
1174             newFaceIndex[5],
1175             oldCells[1]
1176         );
1177
1178         // Third boundary face; Owner = newCell[1] & Neighbour [-1]
1179         tmpTriFace[0] = otherPointIndex[3];
1180         tmpTriFace[1] = newPointIndex[1];
1181         tmpTriFace[2] = otherEdgePoint[3];
1182
1183         newFaceIndex[6] =
1184         (
1185             insertFace
1186             (
1187                 whichPatch(commonFaceIndex[3]),
1188                 tmpTriFace,
1189                 newCellIndex[1],
1190                 -1
1191             )
1192         );
1193
1194         // Add a faceEdges entry as well
1195         faceEdges_.append(tmpTFEdges);
1196
1197         meshOps::replaceLabel
1198         (
1199             -1,
1200             newFaceIndex[6],
1201             newCells[1]
1202         );
1203
1204         // Create / modify edges...
1205         labelList tmpQuadEdgeFaces(4, -1);
1206
1207         // The internal edge bisecting the face
1208         tmpQuadEdgeFaces[0] = fIndex;
1209         tmpQuadEdgeFaces[1] = newFaceIndex[0];
1210         tmpQuadEdgeFaces[2] = newFaceIndex[3];
1211         tmpQuadEdgeFaces[3] = newFaceIndex[4];
1212
1213         newEdgeIndex[0] =
1214         (
1215             insertEdge
1216             (
1217                 -1,
1218                 edge(newPointIndex[0], newPointIndex[1]),
1219                 tmpQuadEdgeFaces
1220             )
1221         );
1222
1223         // Replace an edge on the bisected face
1224         meshOps::replaceLabel
1225         (
1226             otherEdgeIndex[1],
1227             newEdgeIndex[0],
1228             faceEdges_[fIndex]
1229         );
1230
1231         // Create / replace side edges created from face bisection
1232         tmpTriEdgeFaces[0] = commonFaceIndex[2];
1233         tmpTriEdgeFaces[1] = newFaceIndex[3];
1234         tmpTriEdgeFaces[2] = newFaceIndex[1];
1235
1236         newEdgeIndex[3] =
1237         (
1238             insertEdge
1239             (
1240                 whichEdgePatch(commonEdgeIndex[2]),
1241                 edge(newPointIndex[0], nextToOtherPoint[2]),
1242                 tmpTriEdgeFaces
1243             )
1244         );
1245
1246         tmpTriEdgeFaces[0] = c0BdyIndex[1];
1247         tmpTriEdgeFaces[1] = newFaceIndex[3];
1248         tmpTriEdgeFaces[2] = newFaceIndex[6];
1249
1250         newEdgeIndex[4] =
1251         (
1252             insertEdge
1253             (
1254                 whichEdgePatch(commonEdgeIndex[3]),
1255                 edge(newPointIndex[1], otherEdgePoint[3]),
1256                 tmpTriEdgeFaces
1257             )
1258         );
1259
1260         // The edge bisecting the second boundary triangular face
1261         tmpTriEdgeFaces[0] = commonFaceIndex[2];
1262         tmpTriEdgeFaces[1] = newFaceIndex[4];
1263         tmpTriEdgeFaces[2] = newFaceIndex[5];
1264
1265         newEdgeIndex[5] =
1266         (
1267             insertEdge
1268             (
1269                 whichEdgePatch(commonEdgeIndex[2]),
1270                 edge(newPointIndex[0], otherPointIndex[2]),
1271                 tmpTriEdgeFaces
1272             )
1273         );
1274
1275         // Find the common edge between the quad/tri faces...
1276         meshOps::findCommonEdge
1277         (
1278             commonFaceIndex[2],
1279             retainFace,
1280             faceEdges_,
1281             cornerEdgeIndex[2]
1282         );
1283
1284         // ...and correct faceEdges / edgeFaces
1285         meshOps::replaceLabel
1286         (
1287             cornerEdgeIndex[2],
1288             newEdgeIndex[5],
1289             faceEdges_[commonFaceIndex[2]]
1290         );
1291
1292         meshOps::replaceLabel
1293         (
1294             commonFaceIndex[2],
1295             newFaceIndex[5],
1296             edgeFaces_[cornerEdgeIndex[2]]
1297         );
1298
1299         // The edge bisecting the third boundary triangular face
1300         tmpTriEdgeFaces[0] = commonFaceIndex[3];
1301         tmpTriEdgeFaces[1] = newFaceIndex[4];
1302         tmpTriEdgeFaces[2] = newFaceIndex[6];
1303
1304         newEdgeIndex[6] =
1305         (
1306             insertEdge
1307             (
1308                 whichEdgePatch(commonEdgeIndex[3]),
1309                 edge(newPointIndex[1], otherPointIndex[3]),
1310                 tmpTriEdgeFaces
1311             )
1312         );
1313
1314         // Find the common edge between the quad/tri faces...
1315         meshOps::findCommonEdge
1316         (
1317             commonFaceIndex[3],
1318             replaceFace,
1319             faceEdges_,
1320             cornerEdgeIndex[3]
1321         );
1322
1323         // ...and correct faceEdges / edgeFaces
1324         meshOps::replaceLabel
1325         (
1326             cornerEdgeIndex[3],
1327             newEdgeIndex[6],
1328             faceEdges_[commonFaceIndex[3]]
1329         );
1330
1331         meshOps::replaceLabel
1332         (
1333             commonFaceIndex[3],
1334             newFaceIndex[6],
1335             edgeFaces_[cornerEdgeIndex[3]]
1336         );
1337
1338         // Now that edges are defined, configure faceEdges
1339         // for all new faces
1340
1341         // The quad interior face; Owner = c[0] & Neighbour = newCell[0]
1342         tmpQFEdges[0] = newEdgeIndex[0];
1343         tmpQFEdges[1] = newEdgeIndex[1];
1344         tmpQFEdges[2] = otherEdgeIndex[2];
1345         tmpQFEdges[3] = newEdgeIndex[2];
1346         faceEdges_[newFaceIndex[0]] = tmpQFEdges;
1347
1348         // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
1349         tmpTFEdges[0] = newEdgeIndex[3];
1350         tmpTFEdges[1] = cornerEdgeIndex[0];
1351         tmpTFEdges[2] = newEdgeIndex[1];
1352         faceEdges_[newFaceIndex[1]] = tmpTFEdges;
1353
1354         // Third boundary face; Owner = c[0] & Neighbour = [-1]
1355         tmpTFEdges[0] = newEdgeIndex[2];
1356         tmpTFEdges[1] = cornerEdgeIndex[1];
1357         tmpTFEdges[2] = commonEdgeIndex[3];
1358         faceEdges_[newFaceIndex[2]] = tmpTFEdges;
1359
1360         // The quad face from bisection:
1361         tmpQFEdges[0] = otherEdgeIndex[1];
1362         tmpQFEdges[1] = newEdgeIndex[3];
1363         tmpQFEdges[2] = newEdgeIndex[0];
1364         tmpQFEdges[3] = newEdgeIndex[4];
1365         faceEdges_[newFaceIndex[3]] = tmpQFEdges;
1366
1367         // The quad interior face; Owner = c[1] & Neighbour = newCell[1]
1368         tmpQFEdges[0] = newEdgeIndex[0];
1369         tmpQFEdges[1] = newEdgeIndex[5];
1370         tmpQFEdges[2] = otherEdgeIndex[3];
1371         tmpQFEdges[3] = newEdgeIndex[6];
1372         faceEdges_[newFaceIndex[4]] = tmpQFEdges;
1373
1374         // Second boundary face; Owner = c[1] & Neighbour = [-1]
1375         tmpTFEdges[0] = commonEdgeIndex[2];
1376         tmpTFEdges[1] = cornerEdgeIndex[2];
1377         tmpTFEdges[2] = newEdgeIndex[5];
1378         faceEdges_[newFaceIndex[5]] = tmpTFEdges;
1379
1380         // Third boundary face; Owner = newCell[1] & Neighbour = [-1]
1381         tmpTFEdges[0] = newEdgeIndex[4];
1382         tmpTFEdges[1] = cornerEdgeIndex[3];
1383         tmpTFEdges[2] = newEdgeIndex[6];
1384         faceEdges_[newFaceIndex[6]] = tmpTFEdges;
1385
1386         meshOps::replaceLabel
1387         (
1388             commonEdgeIndex[1],
1389             newEdgeIndex[4],
1390             faceEdges_[c0BdyIndex[1]]
1391         );
1392
1393         meshOps::replaceLabel
1394         (
1395             c0BdyIndex[1],
1396             newFaceIndex[2],
1397             edgeFaces_[commonEdgeIndex[1]]
1398         );
1399
1400         meshOps::replaceLabel
1401         (
1402             commonEdgeIndex[2],
1403             newEdgeIndex[3],
1404             faceEdges_[commonFaceIndex[2]]
1405         );
1406
1407         meshOps::replaceLabel
1408         (
1409             commonFaceIndex[2],
1410             newFaceIndex[5],
1411             edgeFaces_[commonEdgeIndex[2]]
1412         );
1413
1414         if (debug > 2)
1415         {
1416             Info << nl << "Modified Cell[0]: "
1417                  << c0 << ": " << oldCells[0] << endl;
1418
1419             forAll(oldCells[0], faceI)
1420             {
1421                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
1422
1423                 Info << oldCells[0][faceI]
1424                      << ": " << faces_[oldCells[0][faceI]]
1425                      << " fE: " << fE
1426                      << endl;
1427
1428                 forAll(fE, edgeI)
1429                 {
1430                     const labelList& eF = edgeFaces_[fE[edgeI]];
1431
1432                     Info << '\t' << fE[edgeI]
1433                          << ": " << edges_[fE[edgeI]]
1434                          << " eF: " << eF
1435                          << endl;
1436                 }
1437             }
1438
1439             Info << "New Cell[0]: "
1440                  << newCellIndex[0] << ": " << newCells[0] << endl;
1441
1442             forAll(newCells[0], faceI)
1443             {
1444                 const labelList& fE = faceEdges_[newCells[0][faceI]];
1445
1446                 Info << newCells[0][faceI] << ": "
1447                      << faces_[newCells[0][faceI]]
1448                      << " fE: " << fE
1449                      << endl;
1450
1451                 forAll(fE, edgeI)
1452                 {
1453                     const labelList& eF = edgeFaces_[fE[edgeI]];
1454
1455                     Info << '\t' << fE[edgeI]
1456                          << ": " << edges_[fE[edgeI]]
1457                          << " eF: " << eF
1458                          << endl;
1459                 }
1460             }
1461
1462             Info << nl << "Modified Cell[1]: "
1463                  << c1 << ": " << oldCells[1] << endl;
1464
1465             forAll(oldCells[1], faceI)
1466             {
1467                 const labelList& fE = faceEdges_[oldCells[1][faceI]];
1468
1469                 Info << oldCells[1][faceI] << ": "
1470                      << faces_[oldCells[1][faceI]]
1471                      << " fE: " << fE
1472                      << endl;
1473
1474                 forAll(fE, edgeI)
1475                 {
1476                     const labelList& eF = edgeFaces_[fE[edgeI]];
1477
1478                     Info << '\t' << fE[edgeI]
1479                          << ": " << edges_[fE[edgeI]]
1480                          << " eF: " << eF
1481                          << endl;
1482                 }
1483             }
1484
1485             Info << "New Cell[1]: "
1486                  << newCellIndex[1] << ": " << newCells[1] << endl;
1487
1488             forAll(newCells[1], faceI)
1489             {
1490                 const labelList& fE = faceEdges_[newCells[1][faceI]];
1491
1492                 Info << newCells[1][faceI] << ": "
1493                      << faces_[newCells[1][faceI]]
1494                      << " fE: " << fE
1495                      << endl;
1496
1497                 forAll(fE, edgeI)
1498                 {
1499                     const labelList& eF = edgeFaces_[fE[edgeI]];
1500
1501                     Info << '\t' << fE[edgeI]
1502                          << ": " << edges_[fE[edgeI]]
1503                          << " eF: " << eF
1504                          << endl;
1505                 }
1506             }
1507         }
1508
1509         // Update the cell list.
1510         cells_[c1] = oldCells[1];
1511         cells_[newCellIndex[1]] = newCells[1];
1512     }
1513
1514     // Update the cell list.
1515     cells_[c0] = oldCells[0];
1516     cells_[newCellIndex[0]] = newCells[0];
1517
1518     // Modify point labels for common edges
1519     if (edges_[commonEdgeIndex[0]].start() == otherEdgePoint[0])
1520     {
1521         edges_[commonEdgeIndex[0]].start() = newPointIndex[0];
1522     }
1523     else
1524     {
1525         edges_[commonEdgeIndex[0]].end() = newPointIndex[0];
1526     }
1527
1528     if (edges_[commonEdgeIndex[1]].start() == nextToOtherPoint[1])
1529     {
1530         edges_[commonEdgeIndex[1]].start() = newPointIndex[1];
1531     }
1532     else
1533     {
1534         edges_[commonEdgeIndex[1]].end() = newPointIndex[1];
1535     }
1536
1537     // Write out VTK files after change
1538     if (debug > 3)
1539     {
1540         labelList cellHull(4, -1);
1541
1542         cellHull[0] = owner_[fIndex];
1543         cellHull[1] = neighbour_[fIndex];
1544         cellHull[2] = owner_[newFaceIndex[3]];
1545         cellHull[3] = neighbour_[newFaceIndex[3]];
1546
1547         writeVTK
1548         (
1549             Foam::name(fIndex)
1550           + "_Bisect_1",
1551             cellHull
1552         );
1553     }
1554
1555     // Fill-in mapping information
1556     FixedList<label, 4> mapCells(-1);
1557
1558     mapCells[0] = c0;
1559     mapCells[1] = newCellIndex[0];
1560
1561     if (c1 != -1)
1562     {
1563         mapCells[2] = c1;
1564         mapCells[3] = newCellIndex[1];
1565     }
1566
1567     labelList mC(1, c0);
1568
1569     forAll(mapCells, cellI)
1570     {
1571         if (mapCells[cellI] == -1)
1572         {
1573             continue;
1574         }
1575
1576         // Set the mapping for this cell
1577         setCellMapping(mapCells[cellI], mC);
1578     }
1579
1580     // Set fill-in mapping information for the modified face.
1581     if (c1 == -1)
1582     {
1583         // Set the mapping for this face
1584         setFaceMapping(fIndex, labelList(1, fIndex));
1585     }
1586     else
1587     {
1588         // Internal face. Default mapping.
1589         setFaceMapping(fIndex);
1590     }
1591
1592     forAll(newFaceIndex, faceI)
1593     {
1594         if (newFaceIndex[faceI] == -1)
1595         {
1596             continue;
1597         }
1598
1599         // Check for boundary faces
1600         if (neighbour_[newFaceIndex[faceI]] == -1)
1601         {
1602             // Boundary face. Compute mapping.
1603             labelList mC;
1604
1605             if (faces_[newFaceIndex[faceI]].size() == 4)
1606             {
1607                 // Quad-face on boundary
1608                 mC.setSize(1, fIndex);
1609             }
1610             else
1611             if (faces_[newFaceIndex[faceI]].size() == 3)
1612             {
1613                 label triFacePatch = whichPatch(newFaceIndex[faceI]);
1614
1615                 // Fetch face-normals
1616                 vector tfNorm, f0Norm, f1Norm;
1617
1618                 tfNorm = faces_[newFaceIndex[faceI]].normal(oldPoints_);
1619                 f0Norm = faces_[c0BdyIndex[0]].normal(oldPoints_);
1620                 f1Norm = faces_[c0BdyIndex[1]].normal(oldPoints_);
1621
1622                 // Tri-face on boundary. Perform normal checks
1623                 // also, because of empty patches.
1624                 if
1625                 (
1626                     (whichPatch(c0BdyIndex[0]) == triFacePatch) &&
1627                     ((tfNorm & f0Norm) > 0.0)
1628                 )
1629                 {
1630                     mC.setSize(1, c0BdyIndex[0]);
1631                 }
1632                 else
1633                 if
1634                 (
1635                     (whichPatch(c0BdyIndex[1]) == triFacePatch) &&
1636                     ((tfNorm & f1Norm) > 0.0)
1637                 )
1638                 {
1639                     mC.setSize(1, c0BdyIndex[1]);
1640                 }
1641                 else
1642                 {
1643                     FatalErrorIn
1644                     (
1645                         "\n"
1646                         "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
1647                         "(\n"
1648                         "    const label fIndex,\n"
1649                         "    const changeMap& masterMap,\n"
1650                         "    bool checkOnly\n"
1651                         ")\n"
1652                     )
1653                         << " Unable to find patch for face: "
1654                         << newFaceIndex[faceI] << ":: "
1655                         << faces_[newFaceIndex[faceI]] << nl
1656                         << " Patch: " << triFacePatch << nl
1657                         << abort(FatalError);
1658                 }
1659             }
1660
1661             // Set the mapping for this face
1662             setFaceMapping(newFaceIndex[faceI], mC);
1663         }
1664         else
1665         {
1666             // Internal quad-faces get default mapping.
1667             setFaceMapping(newFaceIndex[faceI]);
1668         }
1669     }
1670
1671     // Set the flag
1672     topoChangeFlag_ = true;
1673
1674     // Increment the counter
1675     statistics_[3]++;
1676
1677     // Increment surface-counter
1678     if (c1 == -1)
1679     {
1680         statistics_[5]++;
1681     }
1682
1683     // Increment the number of modifications
1684     statistics_[0]++;
1685
1686     // Specify that the operation was successful
1687     map.type() = 1;
1688
1689     // Return the changeMap
1690     return map;
1691 }
1692
1693 // Method for the bisection of an edge in 3D
1694 // - Returns a changeMap with a type specifying:
1695 //     1: Bisection was successful
1696 //    -1: Bisection failed since max number of topo-changes was reached.
1697 //    -2: Bisection failed since resulting quality would be unacceptable.
1698 // - AddedPoints contain the index of the newly added point.
1699 const changeMap dynamicTopoFvMesh::bisectEdge
1700 (
1701     const label eIndex,
1702     bool checkOnly,
1703     bool forceOp
1704 )
1705 {
1706     // Edge bisection performs the following operations:
1707     //      [1] Add a point at middle of the edge
1708     //      [2] Bisect all faces surrounding this edge
1709     //      [3] Bisect all cells surrounding this edge
1710     //      [4] Create internal/external edges for each bisected face
1711     //      [5] Create internal faces for each bisected cell
1712     //      Update faceEdges, edgeFaces and edgePoints information
1713
1714     // For 2D meshes, perform face-bisection
1715     if (twoDMesh_)
1716     {
1717         return bisectQuadFace(eIndex, changeMap(), checkOnly);
1718     }
1719
1720     // Figure out which thread this is...
1721     label tIndex = self();
1722
1723     // Prepare the changeMaps
1724     changeMap map, slaveMap;
1725
1726     if
1727     (
1728         (statistics_[0] > maxModifications_) &&
1729         (maxModifications_ > -1)
1730     )
1731     {
1732         // Reached the max allowable topo-changes.
1733         stack(tIndex).clear();
1734
1735         return map;
1736     }
1737
1738     // Check if edgeRefinements are to be avoided on patch.
1739     if (lengthEstimator().checkRefinementPatch(whichEdgePatch(eIndex)))
1740     {
1741         return map;
1742     }
1743
1744     // Sanity check: Is the index legitimate?
1745     if (eIndex < 0)
1746     {
1747         FatalErrorIn
1748         (
1749             "\n"
1750             "const changeMap dynamicTopoFvMesh::bisectEdge\n"
1751             "(\n"
1752             "    const label eIndex,\n"
1753             "    bool checkOnly,\n"
1754             "    bool forceOp\n"
1755             ")\n"
1756         )
1757             << " Invalid index: " << eIndex
1758             << " nEdges: " << nEdges_
1759             << abort(FatalError);
1760     }
1761
1762     // Before we bisect this edge, check whether the operation will
1763     // yield an acceptable cell-quality.
1764     scalar minQ = 0.0;
1765
1766     if ((minQ = computeBisectionQuality(eIndex)) < sliverThreshold_)
1767     {
1768         // Check if the quality is actually valid before forcing it.
1769         if (forceOp && (minQ < 0.0))
1770         {
1771             FatalErrorIn
1772             (
1773                 "\n"
1774                 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
1775                 "(\n"
1776                 "    const label eIndex,\n"
1777                 "    bool checkOnly,\n"
1778                 "    bool forceOp\n"
1779                 ")\n"
1780             )
1781                 << " Forcing bisection on edge: " << eIndex
1782                 << " will yield an invalid cell."
1783                 << abort(FatalError);
1784         }
1785         else
1786         if (!forceOp)
1787         {
1788             map.type() = -2;
1789             return map;
1790         }
1791     }
1792
1793     // Are we performing only checks?
1794     if (checkOnly)
1795     {
1796         map.type() = 1;
1797         return map;
1798     }
1799
1800     // Update number of surface bisections, if necessary.
1801     if (whichEdgePatch(eIndex) > -1)
1802     {
1803         statistics_[5]++;
1804     }
1805
1806     // Hull variables
1807     face tmpTriFace(3);
1808     labelList tmpEdgeFaces(3,-1);
1809     labelList tmpIntEdgeFaces(4,-1);
1810     labelList tmpEdgePoints(3,-1);
1811     labelList tmpIntEdgePoints(4,-1);
1812     labelList tmpFaceEdges(3,-1);
1813
1814     // Make a copy of existing entities
1815     const labelList vertexHull = edgePoints_[eIndex];
1816     label m = vertexHull.size();
1817
1818     // Size up the hull lists
1819     labelList cellHull(m, -1);
1820     labelList faceHull(m, -1);
1821     labelList edgeHull(m, -1);
1822     labelListList ringEntities(4, labelList(m, -1));
1823
1824     // Construct a hull around this edge
1825     meshOps::constructHull
1826     (
1827         eIndex,
1828         faces_,
1829         edges_,
1830         cells_,
1831         owner_,
1832         neighbour_,
1833         faceEdges_,
1834         edgeFaces_,
1835         edgePoints_,
1836         edgeHull,
1837         faceHull,
1838         cellHull,
1839         ringEntities
1840     );
1841
1842     if (debug > 1)
1843     {
1844         Info << nl << nl
1845              << "Edge: " << eIndex
1846              << ": " << edges_[eIndex]
1847              << " is to be bisected. " << endl;
1848
1849         label epIndex = whichEdgePatch(eIndex);
1850
1851         Info << "Patch: ";
1852
1853         if (epIndex == -1)
1854         {
1855             Info << "Internal" << endl;
1856         }
1857         else
1858         {
1859             Info << boundaryMesh()[epIndex].name() << endl;
1860         }
1861
1862         // Write out VTK files prior to change
1863         if (debug > 3)
1864         {
1865             writeVTK
1866             (
1867                 Foam::name(eIndex)
1868               + "_Bisect_0",
1869                 cellHull
1870             );
1871         }
1872     }
1873
1874     labelList mP(2, -1);
1875
1876     // Set mapping for this point
1877     mP[0] = edges_[eIndex][0];
1878     mP[1] = edges_[eIndex][1];
1879
1880     // Add a new point to the end of the list
1881     label newPointIndex =
1882     (
1883         insertPoint
1884         (
1885             0.5 * (points_[mP[0]] + points_[mP[1]]),
1886             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
1887             mP
1888         )
1889     );
1890
1891     // Add this point to the map.
1892     map.addPoint(newPointIndex);
1893
1894     // New edges can lie on a bounding curve between
1895     // coupled and non-coupled faces. Preferentially
1896     // add edges to coupled-patches, if they exist.
1897     // This makes it convenient for coupled patch matching.
1898     label nePatch = -1;
1899
1900     {
1901         nePatch = whichEdgePatch(eIndex);
1902     }
1903
1904     // Add a new edge to the end of the list
1905     label newEdgeIndex =
1906     (
1907         insertEdge
1908         (
1909             nePatch,
1910             edge(newPointIndex,edges_[eIndex][1]),
1911             labelList(faceHull.size(),-1),
1912             vertexHull
1913         )
1914     );
1915
1916     // Add this edge to the map.
1917     map.addEdge(newEdgeIndex);
1918
1919     // Remove the existing edge from the pointEdges list
1920     // of the modified point, and add it to the new point
1921     meshOps::sizeDownList(eIndex, pointEdges_[edges_[eIndex][1]]);
1922     meshOps::sizeUpList(eIndex, pointEdges_[newPointIndex]);
1923
1924     // Modify the existing edge
1925     edges_[eIndex][1] = newPointIndex;
1926
1927     // Add this edge to the map.
1928     // Although this edge isn't technically 'added', it's
1929     // required for coupled patch mapping.
1930     map.addEdge(eIndex);
1931
1932     // Keep track of added entities
1933     labelList addedCellIndices(cellHull.size(),-1);
1934     labelList addedFaceIndices(faceHull.size(),-1);
1935     labelList addedEdgeIndices(faceHull.size(),-1);
1936     labelList addedIntFaceIndices(faceHull.size(),-1);
1937
1938     // Now loop through the hull and bisect individual entities
1939     forAll(vertexHull, indexI)
1940     {
1941         // Modify the existing face
1942         meshOps::replaceLabel
1943         (
1944             edges_[newEdgeIndex][1],
1945             newPointIndex,
1946             faces_[faceHull[indexI]]
1947         );
1948
1949         // Modify edgePoints for the edge
1950         meshOps::replaceLabel
1951         (
1952             edges_[newEdgeIndex][1],
1953             newPointIndex,
1954             edgePoints_[ringEntities[0][indexI]]
1955         );
1956
1957         // Obtain circular indices
1958         label nextI = vertexHull.fcIndex(indexI);
1959         label prevI = vertexHull.rcIndex(indexI);
1960
1961         // Check if this is an interior/boundary face
1962         if (cellHull[indexI] != -1)
1963         {
1964             // Create a new cell. Add it for now, but update later.
1965             cell newCell(4);
1966
1967             addedCellIndices[indexI] =
1968             (
1969                 insertCell(newCell, lengthScale_[cellHull[indexI]])
1970             );
1971
1972             // Add this cell to the map.
1973             map.addCell(addedCellIndices[indexI]);
1974
1975             // Configure the interior face
1976             tmpTriFace[0] = vertexHull[nextI];
1977             tmpTriFace[1] = vertexHull[indexI];
1978             tmpTriFace[2] = newPointIndex;
1979
1980             // Insert the face
1981             addedIntFaceIndices[indexI] =
1982             (
1983                 insertFace
1984                 (
1985                     -1,
1986                     tmpTriFace,
1987                     cellHull[indexI],
1988                     addedCellIndices[indexI]
1989                 )
1990             );
1991
1992             // Add a faceEdges entry as well
1993             faceEdges_.append(tmpFaceEdges);
1994
1995             // Add this face to the map.
1996             map.addFace(addedIntFaceIndices[indexI]);
1997
1998             // Add to the new cell
1999             newCell[0] = addedIntFaceIndices[indexI];
2000
2001             // Modify the existing ring face connected to newEdge[1]
2002             label replaceFace = ringEntities[3][indexI];
2003
2004             // Check if face reversal is necessary
2005             if (owner_[replaceFace] == cellHull[indexI])
2006             {
2007                 if (neighbour_[replaceFace] == -1)
2008                 {
2009                     // Change the owner
2010                     owner_[replaceFace] = addedCellIndices[indexI];
2011                 }
2012                 else
2013                 {
2014                     // This face has to be reversed
2015                     faces_[replaceFace] = faces_[replaceFace].reverseFace();
2016                     owner_[replaceFace] = neighbour_[replaceFace];
2017                     neighbour_[replaceFace] = addedCellIndices[indexI];
2018
2019                     setFlip(replaceFace);
2020                 }
2021             }
2022             else
2023             {
2024                 // Keep owner, but change neighbour
2025                 neighbour_[replaceFace] = addedCellIndices[indexI];
2026             }
2027
2028             // Modify the edge on the ring.
2029             // Add the new interior face to edgeFaces.
2030             meshOps::sizeUpList
2031             (
2032                 addedIntFaceIndices[indexI],
2033                 edgeFaces_[edgeHull[indexI]]
2034             );
2035
2036             // Insert the new point to edgePoints for the ring edge
2037             meshOps::insertLabel
2038             (
2039                 newPointIndex,
2040                 edges_[eIndex][0],
2041                 edges_[newEdgeIndex][1],
2042                 edgePoints_[edgeHull[indexI]]
2043             );
2044
2045             // Add this edge to faceEdges for the new interior face
2046             faceEdges_[addedIntFaceIndices[indexI]][0] = edgeHull[indexI];
2047
2048             // Replace face labels
2049             meshOps::replaceLabel
2050             (
2051                 replaceFace,
2052                 addedIntFaceIndices[indexI],
2053                 cells_[cellHull[indexI]]
2054             );
2055
2056             // Add to the new cell
2057             newCell[1] = replaceFace;
2058
2059             // Check if this is a boundary face
2060             if (cellHull[prevI] == -1)
2061             {
2062                 // Configure the boundary face
2063                 tmpTriFace[0] = newPointIndex;
2064                 tmpTriFace[1] = edges_[newEdgeIndex][1];
2065                 tmpTriFace[2] = vertexHull[indexI];
2066
2067                 // Insert the face
2068                 addedFaceIndices[indexI] =
2069                 (
2070                     insertFace
2071                     (
2072                         whichPatch(faceHull[indexI]),
2073                         tmpTriFace,
2074                         addedCellIndices[indexI],
2075                         -1
2076                     )
2077                 );
2078
2079                 // Add this face to the map.
2080                 map.addFace(addedFaceIndices[indexI]);
2081
2082                 // Configure edgeFaces
2083                 tmpEdgeFaces[0] = faceHull[indexI];
2084                 tmpEdgeFaces[1] = addedIntFaceIndices[indexI];
2085                 tmpEdgeFaces[2] = addedFaceIndices[indexI];
2086
2087                 // Configure edgePoints
2088                 tmpEdgePoints[0] = edges_[eIndex][0];
2089                 tmpEdgePoints[1] = vertexHull[nextI];
2090                 tmpEdgePoints[2] = edges_[newEdgeIndex][1];
2091
2092                 // Add an edge
2093                 addedEdgeIndices[indexI] =
2094                 (
2095                     insertEdge
2096                     (
2097                         whichPatch(faceHull[indexI]),
2098                         edge(newPointIndex,vertexHull[indexI]),
2099                         tmpEdgeFaces,
2100                         tmpEdgePoints
2101                     )
2102                 );
2103
2104                 // Add this edge to the map.
2105                 map.addEdge(addedEdgeIndices[indexI]);
2106
2107                 // Add this edge to the interior-face faceEdges entry
2108                 faceEdges_[addedIntFaceIndices[indexI]][1] =
2109                 (
2110                     addedEdgeIndices[indexI]
2111                 );
2112
2113                 // Configure faceEdges for this boundary face
2114                 tmpFaceEdges[0] = addedEdgeIndices[indexI];
2115                 tmpFaceEdges[1] = newEdgeIndex;
2116                 tmpFaceEdges[2] = ringEntities[2][indexI];
2117
2118                 // Modify faceEdges for the hull face
2119                 meshOps::replaceLabel
2120                 (
2121                     ringEntities[2][indexI],
2122                     addedEdgeIndices[indexI],
2123                     faceEdges_[faceHull[indexI]]
2124                 );
2125
2126                 // Modify edgeFaces for the edge connected to newEdge[1]
2127                 meshOps::replaceLabel
2128                 (
2129                     faceHull[indexI],
2130                     addedFaceIndices[indexI],
2131                     edgeFaces_[ringEntities[2][indexI]]
2132                 );
2133
2134                 // Modify edgePoints for the edge connected to newEdge[1]
2135                 meshOps::replaceLabel
2136                 (
2137                     edges_[eIndex][0],
2138                     newPointIndex,
2139                     edgePoints_[ringEntities[2][indexI]]
2140                 );
2141
2142                 // Add the faceEdges entry
2143                 faceEdges_.append(tmpFaceEdges);
2144
2145                 // Add an entry to edgeFaces
2146                 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
2147
2148                 // Add an entry for this cell
2149                 newCell[2] = addedFaceIndices[indexI];
2150             }
2151             else
2152             // Check if a cell was added before this
2153             if (addedCellIndices[prevI] != -1)
2154             {
2155                 // Configure the interior face
2156                 tmpTriFace[0] = vertexHull[indexI];
2157                 tmpTriFace[1] = edges_[newEdgeIndex][1];
2158                 tmpTriFace[2] = newPointIndex;
2159
2160                 // Insert the face
2161                 addedFaceIndices[indexI] =
2162                 (
2163                     insertFace
2164                     (
2165                         -1,
2166                         tmpTriFace,
2167                         addedCellIndices[prevI],
2168                         addedCellIndices[indexI]
2169                     )
2170                 );
2171
2172                 // Add this face to the map.
2173                 map.addFace(addedFaceIndices[indexI]);
2174
2175                 // Configure edgeFaces
2176                 tmpIntEdgeFaces[0] = faceHull[indexI];
2177                 tmpIntEdgeFaces[1] = addedIntFaceIndices[indexI];
2178                 tmpIntEdgeFaces[2] = addedFaceIndices[indexI];
2179                 tmpIntEdgeFaces[3] = addedIntFaceIndices[prevI];
2180
2181                 // Configure edgePoints
2182                 tmpIntEdgePoints[0] = edges_[eIndex][0];
2183                 tmpIntEdgePoints[1] = vertexHull[nextI];
2184                 tmpIntEdgePoints[2] = edges_[newEdgeIndex][1];
2185                 tmpIntEdgePoints[3] = vertexHull[prevI];
2186
2187                 // Add an internal edge
2188                 addedEdgeIndices[indexI] =
2189                 (
2190                     insertEdge
2191                     (
2192                         -1,
2193                         edge(newPointIndex,vertexHull[indexI]),
2194                         tmpIntEdgeFaces,
2195                         tmpIntEdgePoints
2196                     )
2197                 );
2198
2199                 // Add this edge to the map.
2200                 map.addEdge(addedEdgeIndices[indexI]);
2201
2202                 // Add this edge to the interior-face faceEdges entry..
2203                 faceEdges_[addedIntFaceIndices[indexI]][1] =
2204                 (
2205                     addedEdgeIndices[indexI]
2206                 );
2207
2208                 // ... and to the previous interior face as well
2209                 faceEdges_[addedIntFaceIndices[prevI]][2] =
2210                 (
2211                     addedEdgeIndices[indexI]
2212                 );
2213
2214                 // Configure faceEdges for this split interior face
2215                 tmpFaceEdges[0] = addedEdgeIndices[indexI];
2216                 tmpFaceEdges[1] = newEdgeIndex;
2217                 tmpFaceEdges[2] = ringEntities[2][indexI];
2218
2219                 // Modify faceEdges for the hull face
2220                 meshOps::replaceLabel
2221                 (
2222                     ringEntities[2][indexI],
2223                     addedEdgeIndices[indexI],
2224                     faceEdges_[faceHull[indexI]]
2225                 );
2226
2227                 // Modify edgeFaces for the edge connected to newEdge[1]
2228                 meshOps::replaceLabel
2229                 (
2230                     faceHull[indexI],
2231                     addedFaceIndices[indexI],
2232                     edgeFaces_[ringEntities[2][indexI]]
2233                 );
2234
2235                 // Modify edgePoints for the edge connected to newEdge[1]
2236                 meshOps::replaceLabel
2237                 (
2238                     edges_[eIndex][0],
2239                     newPointIndex,
2240                     edgePoints_[ringEntities[2][indexI]]
2241                 );
2242
2243                 // Add the faceEdges entry
2244                 faceEdges_.append(tmpFaceEdges);
2245
2246                 // Add an entry to edgeFaces
2247                 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
2248
2249                 // Add an entry for this cell
2250                 newCell[2] = addedFaceIndices[indexI];
2251
2252                 // Make the final entry for the previous cell
2253                 cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
2254             }
2255
2256             // Do the first interior face at the end
2257             if (indexI == vertexHull.size() - 1)
2258             {
2259                 // Configure the interior face
2260                 tmpTriFace[0] = newPointIndex;
2261                 tmpTriFace[1] = edges_[newEdgeIndex][1];
2262                 tmpTriFace[2] = vertexHull[0];
2263
2264                 // Insert the face
2265                 addedFaceIndices[0] =
2266                 (
2267                     insertFace
2268                     (
2269                         -1,
2270                         tmpTriFace,
2271                         addedCellIndices[0],
2272                         addedCellIndices[indexI]
2273                     )
2274                 );
2275
2276                 // Add this face to the map.
2277                 map.addFace(addedFaceIndices[0]);
2278
2279                 // Configure edgeFaces
2280                 tmpIntEdgeFaces[0] = faceHull[0];
2281                 tmpIntEdgeFaces[1] = addedIntFaceIndices[0];
2282                 tmpIntEdgeFaces[2] = addedFaceIndices[0];
2283                 tmpIntEdgeFaces[3] = addedIntFaceIndices[indexI];
2284
2285                 // Configure edgePoints
2286                 tmpIntEdgePoints[0] = edges_[eIndex][0];
2287                 tmpIntEdgePoints[1] = vertexHull[1];
2288                 tmpIntEdgePoints[2] = edges_[newEdgeIndex][1];
2289                 tmpIntEdgePoints[3] = vertexHull[indexI];
2290
2291                 // Add an internal edge
2292                 addedEdgeIndices[0] =
2293                 (
2294                     insertEdge
2295                     (
2296                         -1,
2297                         edge(newPointIndex,vertexHull[0]),
2298                         tmpIntEdgeFaces,
2299                         tmpIntEdgePoints
2300                     )
2301                 );
2302
2303                 // Add this edge to the map.
2304                 map.addEdge(addedEdgeIndices[0]);
2305
2306                 // Add this edge to the interior-face faceEdges entry..
2307                 faceEdges_[addedIntFaceIndices[0]][1] =
2308                 (
2309                     addedEdgeIndices[0]
2310                 );
2311
2312                 // ... and to the previous interior face as well
2313                 faceEdges_[addedIntFaceIndices[indexI]][2] =
2314                 (
2315                     addedEdgeIndices[0]
2316                 );
2317
2318                 // Configure faceEdges for the first split face
2319                 tmpFaceEdges[0] = addedEdgeIndices[0];
2320                 tmpFaceEdges[1] = newEdgeIndex;
2321                 tmpFaceEdges[2] = ringEntities[2][0];
2322
2323                 // Modify faceEdges for the hull face
2324                 meshOps::replaceLabel
2325                 (
2326                     ringEntities[2][0],
2327                     addedEdgeIndices[0],
2328                     faceEdges_[faceHull[0]]
2329                 );
2330
2331                 // Modify edgeFaces for the edge connected to newEdge[1]
2332                 meshOps::replaceLabel
2333                 (
2334                     faceHull[0],
2335                     addedFaceIndices[0],
2336                     edgeFaces_[ringEntities[2][0]]
2337                 );
2338
2339                 // Modify edgePoints for the edge connected to newEdge[1]
2340                 meshOps::replaceLabel
2341                 (
2342                     edges_[eIndex][0],
2343                     newPointIndex,
2344                     edgePoints_[ringEntities[2][0]]
2345                 );
2346
2347                 // Add the faceEdges entry
2348                 faceEdges_.append(tmpFaceEdges);
2349
2350                 // Add an entry to edgeFaces
2351                 edgeFaces_[newEdgeIndex][0] = addedFaceIndices[0];
2352
2353                 // Add an entry for this cell
2354                 newCell[3] = addedFaceIndices[0];
2355
2356                 // Make the final entry for the first cell
2357                 cells_[addedCellIndices[0]][2] = addedFaceIndices[0];
2358             }
2359
2360             // Update the cell list with the new cell.
2361             cells_[addedCellIndices[indexI]] = newCell;
2362         }
2363         else
2364         {
2365             // Configure the final boundary face
2366             tmpTriFace[0] = vertexHull[indexI];
2367             tmpTriFace[1] = edges_[newEdgeIndex][1];
2368             tmpTriFace[2] = newPointIndex;
2369
2370             // Insert the face
2371             addedFaceIndices[indexI] =
2372             (
2373                 insertFace
2374                 (
2375                     whichPatch(faceHull[indexI]),
2376                     tmpTriFace,
2377                     addedCellIndices[prevI],
2378                     -1
2379                 )
2380             );
2381
2382             // Add this face to the map.
2383             map.addFace(addedFaceIndices[indexI]);
2384
2385             // Configure edgeFaces
2386             tmpEdgeFaces[0] = addedFaceIndices[indexI];
2387             tmpEdgeFaces[1] = addedIntFaceIndices[prevI];
2388             tmpEdgeFaces[2] = faceHull[indexI];
2389
2390             // Configure edgePoints
2391             tmpEdgePoints[0] = edges_[newEdgeIndex][1];
2392             tmpEdgePoints[1] = vertexHull[prevI];
2393             tmpEdgePoints[2] = edges_[eIndex][0];
2394
2395             // Add an edge
2396             addedEdgeIndices[indexI] =
2397             (
2398                 insertEdge
2399                 (
2400                     whichPatch(faceHull[indexI]),
2401                     edge(newPointIndex,vertexHull[indexI]),
2402                     tmpEdgeFaces,
2403                     tmpEdgePoints
2404                 )
2405             );
2406
2407             // Add this edge to the map.
2408             map.addEdge(addedEdgeIndices[indexI]);
2409
2410             // Add a faceEdges entry to the previous interior face
2411             faceEdges_[addedIntFaceIndices[prevI]][2] =
2412             (
2413                 addedEdgeIndices[indexI]
2414             );
2415
2416             // Configure faceEdges for the final boundary face
2417             tmpFaceEdges[0] = addedEdgeIndices[indexI];
2418             tmpFaceEdges[1] = newEdgeIndex;
2419             tmpFaceEdges[2] = ringEntities[2][indexI];
2420
2421             // Modify faceEdges for the hull face
2422             meshOps::replaceLabel
2423             (
2424                 ringEntities[2][indexI],
2425                 addedEdgeIndices[indexI],
2426                 faceEdges_[faceHull[indexI]]
2427             );
2428
2429             // Modify edgeFaces for the edge connected to newEdge[1]
2430             meshOps::replaceLabel
2431             (
2432                 faceHull[indexI],
2433                 addedFaceIndices[indexI],
2434                 edgeFaces_[ringEntities[2][indexI]]
2435             );
2436
2437             // Modify edgePoints for the edge connected to newEdge[1]
2438             meshOps::replaceLabel
2439             (
2440                 edges_[eIndex][0],
2441                 newPointIndex,
2442                 edgePoints_[ringEntities[2][indexI]]
2443             );
2444
2445             // Add the faceEdges entry
2446             faceEdges_.append(tmpFaceEdges);
2447
2448             // Add an entry to edgeFaces
2449             edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
2450
2451             // Make the final entry for the previous cell
2452             cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
2453         }
2454     }
2455
2456     // Now that all old / new cells possess correct connectivity,
2457     // compute mapping information.
2458     forAll(cellHull, indexI)
2459     {
2460         if (cellHull[indexI] == -1)
2461         {
2462             continue;
2463         }
2464
2465         // Set mapping for both new and modified cells.
2466         FixedList<label, 2> cmIndex;
2467
2468         cmIndex[0] = cellHull[indexI];
2469         cmIndex[1] = addedCellIndices[indexI];
2470
2471         // Fill-in candidate mapping information
2472         labelList mC(1, cellHull[indexI]);
2473
2474         forAll(cmIndex, cmI)
2475         {
2476             // Set the mapping for this cell
2477             setCellMapping(cmIndex[cmI], mC);
2478         }
2479     }
2480
2481     // Set mapping information for old / new faces
2482     forAll(faceHull, indexI)
2483     {
2484         // Interior faces get default mapping
2485         if (addedIntFaceIndices[indexI] > -1)
2486         {
2487             setFaceMapping(addedIntFaceIndices[indexI]);
2488         }
2489
2490         // Decide between default / weighted mapping
2491         // based on boundary information
2492         if (whichPatch(faceHull[indexI]) == -1)
2493         {
2494             // Interior faces get default mapping
2495             setFaceMapping(faceHull[indexI]);
2496             setFaceMapping(addedFaceIndices[indexI]);
2497         }
2498         else
2499         {
2500             // Compute mapping weights for boundary faces
2501             FixedList<label, 2> fmIndex;
2502
2503             fmIndex[0] = faceHull[indexI];
2504             fmIndex[1] = addedFaceIndices[indexI];
2505
2506             // Fill-in candidate mapping information
2507             labelList mF(1, faceHull[indexI]);
2508
2509             forAll(fmIndex, fmI)
2510             {
2511                 // Set the mapping for this face
2512                 setFaceMapping(fmIndex[fmI], mF);
2513             }
2514         }
2515     }
2516
2517     if (debug > 2)
2518     {
2519         label bPatch = whichEdgePatch(eIndex);
2520
2521         if (bPatch == -1)
2522         {
2523             Info << "Patch: Internal" << endl;
2524         }
2525         else
2526         {
2527             Info << "Patch: " << boundaryMesh()[bPatch].name() << endl;
2528         }
2529
2530         Info << "EdgePoints: " << vertexHull << endl;
2531         Info << "Edges: " << edgeHull << endl;
2532         Info << "Faces: " << faceHull << endl;
2533         Info << "Cells: " << cellHull << endl;
2534
2535         Info << "Modified cells: " << endl;
2536
2537         forAll(cellHull, cellI)
2538         {
2539             if (cellHull[cellI] == -1)
2540             {
2541                 continue;
2542             }
2543
2544             Info << cellHull[cellI] << ":: "
2545                  << cells_[cellHull[cellI]]
2546                  << endl;
2547         }
2548
2549         Info << "Added cells: " << endl;
2550
2551         forAll(addedCellIndices, cellI)
2552         {
2553             if (addedCellIndices[cellI] == -1)
2554             {
2555                 continue;
2556             }
2557
2558             Info << addedCellIndices[cellI] << ":: "
2559                  << cells_[addedCellIndices[cellI]] << nl
2560                  << "lengthScale: " << lengthScale_[addedCellIndices[cellI]]
2561                  << endl;
2562         }
2563
2564         Info << "Modified faces: " << endl;
2565
2566         forAll(faceHull, faceI)
2567         {
2568             Info << faceHull[faceI] << ":: "
2569                  << faces_[faceHull[faceI]] << ": "
2570                  << owner_[faceHull[faceI]] << ": "
2571                  << neighbour_[faceHull[faceI]] << " "
2572                  << "faceEdges:: " << faceEdges_[faceHull[faceI]]
2573                  << endl;
2574         }
2575
2576         Info << "Added faces: " << endl;
2577
2578         forAll(addedFaceIndices, faceI)
2579         {
2580             Info << addedFaceIndices[faceI] << ":: "
2581                  << faces_[addedFaceIndices[faceI]] << ": "
2582                  << owner_[addedFaceIndices[faceI]] << ": "
2583                  << neighbour_[addedFaceIndices[faceI]] << " "
2584                  << "faceEdges:: " << faceEdges_[addedFaceIndices[faceI]]
2585                  << endl;
2586         }
2587
2588         forAll(addedIntFaceIndices, faceI)
2589         {
2590             if (addedIntFaceIndices[faceI] == -1)
2591             {
2592                 continue;
2593             }
2594
2595             Info << addedIntFaceIndices[faceI] << ":: "
2596                  << faces_[addedIntFaceIndices[faceI]] << ": "
2597                  << owner_[addedIntFaceIndices[faceI]] << ": "
2598                  << neighbour_[addedIntFaceIndices[faceI]] << " "
2599                  << "faceEdges:: "
2600                  << faceEdges_[addedIntFaceIndices[faceI]]
2601                  << endl;
2602         }
2603
2604         Info << "New edge:: " << newEdgeIndex
2605              << ": " << edges_[newEdgeIndex] << nl
2606              << " edgeFaces:: " << edgeFaces_[newEdgeIndex] << nl
2607              << " edgePoints:: " << edgePoints_[newEdgeIndex]
2608              << endl;
2609
2610         Info << "Added edges: " << endl;
2611
2612         forAll(addedEdgeIndices, edgeI)
2613         {
2614             Info << addedEdgeIndices[edgeI]
2615                  << ":: " << edges_[addedEdgeIndices[edgeI]] << nl
2616                  << " edgeFaces:: " << edgeFaces_[addedEdgeIndices[edgeI]] << nl
2617                  << " edgePoints:: " << edgePoints_[addedEdgeIndices[edgeI]]
2618                  << endl;
2619         }
2620
2621         Info << "New Point:: " << newPointIndex << endl;
2622         Info << "pointEdges:: " << pointEdges_[newPointIndex] << endl;
2623
2624         // Write out VTK files after change
2625         if (debug > 3)
2626         {
2627             labelList newHull(cellHull.size() + addedCellIndices.size(), 0);
2628
2629             // Combine both lists into one.
2630             forAll(cellHull, i)
2631             {
2632                 newHull[i] = cellHull[i];
2633             }
2634
2635             label start = cellHull.size();
2636
2637             for(label i = start; i < newHull.size(); i++)
2638             {
2639                 newHull[i] = addedCellIndices[i - start];
2640             }
2641
2642             writeVTK
2643             (
2644                 Foam::name(eIndex)
2645               + "_Bisect_1",
2646                 newHull
2647             );
2648         }
2649     }
2650
2651     // Set the flag
2652     topoChangeFlag_ = true;
2653
2654     // Increment the counter
2655     statistics_[3]++;
2656
2657     // Increment the number of modifications
2658     statistics_[0]++;
2659
2660     // Specify that the operation was successful
2661     map.type() = 1;
2662
2663     // Return the changeMap
2664     return map;
2665 }
2666
2667 // Method for the trisection of a face in 3D
2668 // - Returns a changeMap with a type specifying:
2669 //     1: Trisection was successful
2670 //    -1: Trisection failed since max number of topo-changes was reached.
2671 //    -2: Trisection failed since resulting quality would be really bad.
2672 // - AddedPoint is the index of the newly added point.
2673 const changeMap dynamicTopoFvMesh::trisectFace
2674 (
2675     const label fIndex,
2676     bool checkOnly,
2677     bool forceOp
2678 )
2679 {
2680     // Face trisection performs the following operations:
2681     //      [1] Add a point at middle of the face
2682     //      [2] Remove the face and add three new faces in place.
2683     //      [3] Add three cells for each trisected cell (remove the originals).
2684     //      [4] Create one internal edge for each trisected cell.
2685     //      [5] Create three edges for the trisected face.
2686     //      [6] Create three internal faces for each trisected cell.
2687     //      Update faceEdges, edgeFaces and edgePoints information.
2688
2689     // Figure out which thread this is...
2690     label tIndex = self();
2691
2692     // Prepare the changeMaps
2693     changeMap map, slaveMap;
2694
2695     if
2696     (
2697         (statistics_[0] > maxModifications_) &&
2698         (maxModifications_ > -1)
2699     )
2700     {
2701         // Reached the max allowable topo-changes.
2702         stack(tIndex).clear();
2703
2704         return map;
2705     }
2706
2707     // Sanity check: Is the index legitimate?
2708     if (fIndex < 0)
2709     {
2710         FatalErrorIn
2711         (
2712             "const changeMap dynamicTopoFvMesh::trisectFace\n"
2713             "(\n"
2714             "    const label fIndex,\n"
2715             "    bool checkOnly,\n"
2716             "    bool forceOp\n"
2717             ")\n"
2718         )
2719             << " Invalid index: " << fIndex
2720             << abort(FatalError);
2721     }
2722
2723     // Before we trisect this face, check whether the operation will
2724     // yield an acceptable cell-quality.
2725     scalar minQ = 0.0;
2726
2727     if ((minQ = computeTrisectionQuality(fIndex)) < sliverThreshold_)
2728     {
2729         // Check if the quality is actually valid before forcing it.
2730         if (forceOp && (minQ < 0.0))
2731         {
2732             FatalErrorIn
2733             (
2734                 "const changeMap dynamicTopoFvMesh::trisectFace\n"
2735                 "(\n"
2736                 "    const label fIndex,\n"
2737                 "    bool checkOnly,\n"
2738                 "    bool forceOp\n"
2739                 ")\n"
2740             )
2741                 << " Forcing trisection on face: " << fIndex
2742                 << " will yield an invalid cell."
2743                 << abort(FatalError);
2744         }
2745         else
2746         if (!forceOp)
2747         {
2748             map.type() = -2;
2749             return map;
2750         }
2751     }
2752
2753     // Are we performing only checks?
2754     if (checkOnly)
2755     {
2756         map.type() = 1;
2757         return map;
2758     }
2759
2760     // Update number of surface bisections, if necessary.
2761     if (whichPatch(fIndex) > -1)
2762     {
2763         statistics_[5]++;
2764     }
2765
2766     // Hull variables
2767     face tmpTriFace(3);
2768     labelList newTriEdgeFaces(3), newTriEdgePoints(3);
2769     labelList newQuadEdgeFaces(4), newQuadEdgePoints(4);
2770
2771     FixedList<label,2> apexPoint(-1);
2772     FixedList<face, 3> checkFace(face(3));
2773     FixedList<label,5> newEdgeIndex(-1);
2774     FixedList<label,9> newFaceIndex(-1);
2775     FixedList<label,6> newCellIndex(-1);
2776     FixedList<cell, 6> newTetCell(cell(4));
2777     FixedList<labelList, 9> newFaceEdges(labelList(3));
2778
2779     // Counters for entities
2780     FixedList<label, 9> nE(0);
2781     FixedList<label, 6> nF(0);
2782
2783     // Determine the two cells to be removed
2784     FixedList<label,2> cellsForRemoval;
2785     cellsForRemoval[0] = owner_[fIndex];
2786     cellsForRemoval[1] = neighbour_[fIndex];
2787
2788     if (debug > 1)
2789     {
2790         Info << nl << nl
2791              << "Face: " << fIndex
2792              << ": " << faces_[fIndex]
2793              << " is to be trisected. " << endl;
2794
2795         // Write out VTK files prior to change
2796         if (debug > 3)
2797         {
2798             labelList vtkCells;
2799
2800             if (neighbour_[fIndex] == -1)
2801             {
2802                 vtkCells.setSize(1);
2803                 vtkCells[0] = owner_[fIndex];
2804             }
2805             else
2806             {
2807                 vtkCells.setSize(2);
2808                 vtkCells[0] = owner_[fIndex];
2809                 vtkCells[1] = neighbour_[fIndex];
2810             }
2811
2812             writeVTK
2813             (
2814                 Foam::name(fIndex)
2815               + "Trisect_0",
2816                 vtkCells
2817             );
2818         }
2819     }
2820
2821     labelList mP(3, -1);
2822
2823     // Fill in mapping information
2824     mP[0] = faces_[fIndex][0];
2825     mP[1] = faces_[fIndex][1];
2826     mP[2] = faces_[fIndex][2];
2827
2828     // Add a new point to the end of the list
2829     scalar oT = (1.0/3.0);
2830
2831     label newPointIndex =
2832     (
2833         insertPoint
2834         (
2835             oT * (points_[mP[0]] + points_[mP[1]] + points_[mP[2]]),
2836             oT * (oldPoints_[mP[0]] + oldPoints_[mP[1]] + oldPoints_[mP[2]]),
2837             mP
2838         )
2839     );
2840
2841     // Add this point to the map.
2842     map.addPoint(newPointIndex);
2843
2844     // Add three new cells to the end of the cell list
2845     for (label i = 0; i < 3; i++)
2846     {
2847         scalar parentScale = -1.0;
2848
2849         if (edgeRefinement_)
2850         {
2851             parentScale = lengthScale_[cellsForRemoval[0]];
2852         }
2853
2854         newCellIndex[i] = insertCell(newTetCell[i], parentScale);
2855
2856         // Add cells to the map
2857         map.addCell(newCellIndex[i]);
2858     }
2859
2860     // Find the apex point for this cell
2861     apexPoint[0] =
2862     (
2863         meshOps::tetApexPoint
2864         (
2865             owner_[fIndex],
2866             fIndex,
2867             faces_,
2868             cells_
2869         )
2870     );
2871
2872     // Insert three new internal faces
2873
2874     // First face: Owner: newCellIndex[0], Neighbour: newCellIndex[1]
2875     tmpTriFace[0] = newPointIndex;
2876     tmpTriFace[1] = faces_[fIndex][0];
2877     tmpTriFace[2] = apexPoint[0];
2878
2879     newFaceIndex[0] =
2880     (
2881         insertFace
2882         (
2883             -1,
2884             tmpTriFace,
2885             newCellIndex[0],
2886             newCellIndex[1]
2887         )
2888     );
2889
2890     // Second face: Owner: newCellIndex[1], Neighbour: newCellIndex[2]
2891     tmpTriFace[0] = newPointIndex;
2892     tmpTriFace[1] = faces_[fIndex][1];
2893     tmpTriFace[2] = apexPoint[0];
2894
2895     newFaceIndex[1] =
2896     (
2897         insertFace
2898         (
2899             -1,
2900             tmpTriFace,
2901             newCellIndex[1],
2902             newCellIndex[2]
2903         )
2904     );
2905
2906     // Third face: Owner: newCellIndex[0], Neighbour: newCellIndex[2]
2907     tmpTriFace[0] = newPointIndex;
2908     tmpTriFace[1] = apexPoint[0];
2909     tmpTriFace[2] = faces_[fIndex][2];
2910
2911     newFaceIndex[2] =
2912     (
2913         insertFace
2914         (
2915             -1,
2916             tmpTriFace,
2917             newCellIndex[0],
2918             newCellIndex[2]
2919         )
2920     );
2921
2922     // Add an entry to edgeFaces
2923     newTriEdgeFaces[0] = newFaceIndex[0];
2924     newTriEdgeFaces[1] = newFaceIndex[1];
2925     newTriEdgeFaces[2] = newFaceIndex[2];
2926
2927     // Add an entry for edgePoints as well
2928     newTriEdgePoints[0] = faces_[fIndex][0];
2929     newTriEdgePoints[1] = faces_[fIndex][1];
2930     newTriEdgePoints[2] = faces_[fIndex][2];
2931
2932     // Add a new internal edge to the mesh
2933     newEdgeIndex[0] =
2934     (
2935         insertEdge
2936         (
2937             -1,
2938             edge
2939             (
2940                newPointIndex,
2941                apexPoint[0]
2942             ),
2943             newTriEdgeFaces,
2944             newTriEdgePoints
2945         )
2946     );
2947
2948     // Configure faceEdges with the new internal edge
2949     newFaceEdges[0][nE[0]++] = newEdgeIndex[0];
2950     newFaceEdges[1][nE[1]++] = newEdgeIndex[0];
2951     newFaceEdges[2][nE[2]++] = newEdgeIndex[0];
2952
2953     // Add the newly created faces to cells
2954     newTetCell[0][nF[0]++] = newFaceIndex[0];
2955     newTetCell[0][nF[0]++] = newFaceIndex[2];
2956     newTetCell[1][nF[1]++] = newFaceIndex[0];
2957     newTetCell[1][nF[1]++] = newFaceIndex[1];
2958     newTetCell[2][nF[2]++] = newFaceIndex[1];
2959     newTetCell[2][nF[2]++] = newFaceIndex[2];
2960
2961     // Define the three faces to check for orientation:
2962     checkFace[0][0] = faces_[fIndex][2];
2963     checkFace[0][1] = apexPoint[0];
2964     checkFace[0][2] = faces_[fIndex][0];
2965
2966     checkFace[1][0] = faces_[fIndex][0];
2967     checkFace[1][1] = apexPoint[0];
2968     checkFace[1][2] = faces_[fIndex][1];
2969
2970     checkFace[2][0] = faces_[fIndex][1];
2971     checkFace[2][1] = apexPoint[0];
2972     checkFace[2][2] = faces_[fIndex][2];
2973
2974     // Check the orientation of faces on the first cell.
2975     forAll(cells_[owner_[fIndex]], faceI)
2976     {
2977         label faceIndex = cells_[owner_[fIndex]][faceI];
2978
2979         if (faceIndex == fIndex)
2980         {
2981             continue;
2982         }
2983
2984         const face& faceToCheck = faces_[faceIndex];
2985         label cellIndex = cellsForRemoval[0];
2986         label newIndex = -1;
2987
2988         // Check against faces.
2989         if (triFace::compare(triFace(faceToCheck), triFace(checkFace[0])))
2990         {
2991             newIndex = newCellIndex[0];
2992             newTetCell[0][nF[0]++] = faceIndex;
2993         }
2994         else
2995         if (triFace::compare(triFace(faceToCheck), triFace(checkFace[1])))
2996         {
2997             newIndex = newCellIndex[1];
2998             newTetCell[1][nF[1]++] = faceIndex;
2999         }
3000         else
3001         if (triFace::compare(triFace(faceToCheck), triFace(checkFace[2])))
3002         {
3003             newIndex = newCellIndex[2];
3004             newTetCell[2][nF[2]++] = faceIndex;
3005         }
3006         else
3007         {
3008             // Something's terribly wrong.
3009             FatalErrorIn
3010             (
3011                 "const changeMap dynamicTopoFvMesh::trisectFace\n"
3012                 "(\n"
3013                 "    const label fIndex,\n"
3014                 "    bool checkOnly,\n"
3015                 "    bool forceOp\n"
3016                 ")\n"
3017             )
3018                 << "Failed to determine a face match."
3019                 << abort(FatalError);
3020         }
3021
3022         // Check if a face-flip is necessary
3023         if (owner_[faceIndex] == cellIndex)
3024         {
3025             if (neighbour_[faceIndex] == -1)
3026             {
3027                 // Change the owner
3028                 owner_[faceIndex] = newIndex;
3029             }
3030             else
3031             {
3032                 // Flip this face
3033                 faces_[faceIndex] = faceToCheck.reverseFace();
3034                 owner_[faceIndex] = neighbour_[faceIndex];
3035                 neighbour_[faceIndex] = newIndex;
3036
3037                 setFlip(faceIndex);
3038             }
3039         }
3040         else
3041         {
3042             // Flip is unnecessary. Just update neighbour
3043             neighbour_[faceIndex] = newIndex;
3044         }
3045     }
3046
3047     if (cellsForRemoval[1] == -1)
3048     {
3049         // Boundary face. Determine its patch.
3050         label facePatch = whichPatch(fIndex);
3051
3052         // Add three new boundary faces.
3053
3054         // Fourth face: Owner: newCellIndex[0], Neighbour: -1
3055         tmpTriFace[0] = newPointIndex;
3056         tmpTriFace[1] = faces_[fIndex][2];
3057         tmpTriFace[2] = faces_[fIndex][0];
3058
3059         newFaceIndex[3] =
3060         (
3061             insertFace
3062             (
3063                 facePatch,
3064                 tmpTriFace,
3065                 newCellIndex[0],
3066                 -1
3067             )
3068         );
3069
3070         // Fifth face: Owner: newCellIndex[1], Neighbour: -1
3071         tmpTriFace[0] = newPointIndex;
3072         tmpTriFace[1] = faces_[fIndex][0];
3073         tmpTriFace[2] = faces_[fIndex][1];
3074
3075         newFaceIndex[4] =
3076         (
3077             insertFace
3078             (
3079                 facePatch,
3080                 tmpTriFace,
3081                 newCellIndex[1],
3082                 -1
3083             )
3084         );
3085
3086         // Sixth face: Owner: newCellIndex[2], Neighbour: -1
3087         tmpTriFace[0] = newPointIndex;
3088         tmpTriFace[1] = faces_[fIndex][1];
3089         tmpTriFace[2] = faces_[fIndex][2];
3090
3091         newFaceIndex[5] =
3092         (
3093             insertFace
3094             (
3095                 facePatch,
3096                 tmpTriFace,
3097                 newCellIndex[2],
3098                 -1
3099             )
3100         );
3101
3102         // Add the newly created faces to cells
3103         newTetCell[0][nF[0]++] = newFaceIndex[3];
3104         newTetCell[1][nF[1]++] = newFaceIndex[4];
3105         newTetCell[2][nF[2]++] = newFaceIndex[5];
3106
3107         // Configure edgeFaces and edgePoints for three new boundary edges.
3108         newTriEdgeFaces[0] = newFaceIndex[4];
3109         newTriEdgeFaces[1] = newFaceIndex[0];
3110         newTriEdgeFaces[2] = newFaceIndex[3];
3111
3112         newTriEdgePoints[0] = faces_[fIndex][1];
3113         newTriEdgePoints[1] = apexPoint[0];
3114         newTriEdgePoints[2] = faces_[fIndex][2];
3115
3116         newEdgeIndex[1] =
3117         (
3118             insertEdge
3119             (
3120                 facePatch,
3121                 edge
3122                 (
3123                    newPointIndex,
3124                    faces_[fIndex][0]
3125                 ),
3126                 newTriEdgeFaces,
3127                 newTriEdgePoints
3128             )
3129         );
3130
3131         newTriEdgeFaces[0] = newFaceIndex[5];
3132         newTriEdgeFaces[1] = newFaceIndex[1];
3133         newTriEdgeFaces[2] = newFaceIndex[4];
3134
3135         newTriEdgePoints[0] = faces_[fIndex][2];
3136         newTriEdgePoints[1] = apexPoint[0];
3137         newTriEdgePoints[2] = faces_[fIndex][0];
3138
3139         newEdgeIndex[2] =
3140         (
3141             insertEdge
3142             (
3143                 facePatch,
3144                 edge
3145                 (
3146                    newPointIndex,
3147                    faces_[fIndex][1]
3148                 ),
3149                 newTriEdgeFaces,
3150                 newTriEdgePoints
3151             )
3152         );
3153
3154         newTriEdgeFaces[0] = newFaceIndex[3];
3155         newTriEdgeFaces[1] = newFaceIndex[2];
3156         newTriEdgeFaces[2] = newFaceIndex[5];
3157
3158         newTriEdgePoints[0] = faces_[fIndex][0];
3159         newTriEdgePoints[1] = apexPoint[0];
3160         newTriEdgePoints[2] = faces_[fIndex][1];
3161
3162         newEdgeIndex[3] =
3163         (
3164             insertEdge
3165             (
3166                 facePatch,
3167                 edge
3168                 (
3169                    newPointIndex,
3170                    faces_[fIndex][2]
3171                 ),
3172                 newTriEdgeFaces,
3173                 newTriEdgePoints
3174             )
3175         );
3176
3177         // Configure faceEdges with the three new edges.
3178         newFaceEdges[0][nE[0]++] = newEdgeIndex[1];
3179         newFaceEdges[1][nE[1]++] = newEdgeIndex[2];
3180         newFaceEdges[2][nE[2]++] = newEdgeIndex[3];
3181
3182         newFaceEdges[3][nE[3]++] = newEdgeIndex[1];
3183         newFaceEdges[3][nE[3]++] = newEdgeIndex[3];
3184         newFaceEdges[4][nE[4]++] = newEdgeIndex[1];
3185         newFaceEdges[4][nE[4]++] = newEdgeIndex[2];
3186         newFaceEdges[5][nE[5]++] = newEdgeIndex[2];
3187         newFaceEdges[5][nE[5]++] = newEdgeIndex[3];
3188
3189         // Define the six edges to check while building faceEdges:
3190         FixedList<edge,6> check;
3191
3192         check[0][0] = apexPoint[0]; check[0][1] = faces_[fIndex][0];
3193         check[1][0] = apexPoint[0]; check[1][1] = faces_[fIndex][1];
3194         check[2][0] = apexPoint[0]; check[2][1] = faces_[fIndex][2];
3195
3196         check[3][0] = faces_[fIndex][2]; check[3][1] = faces_[fIndex][0];
3197         check[4][0] = faces_[fIndex][0]; check[4][1] = faces_[fIndex][1];
3198         check[5][0] = faces_[fIndex][1]; check[5][1] = faces_[fIndex][2];
3199
3200         // Build a list of cellEdges
3201         labelHashSet cellEdges;
3202
3203         forAll(cells_[owner_[fIndex]], faceI)
3204         {
3205             const labelList& fEdges =
3206             (
3207                 faceEdges_[cells_[owner_[fIndex]][faceI]]
3208             );
3209
3210             forAll(fEdges, edgeI)
3211             {
3212                 if (!cellEdges.found(fEdges[edgeI]))
3213                 {
3214                     cellEdges.insert(fEdges[edgeI]);
3215                 }
3216             }
3217         }
3218
3219         // Loop through cellEdges, and perform appropriate actions.
3220         forAllIter(labelHashSet, cellEdges, eIter)
3221         {
3222             const edge& edgeToCheck = edges_[eIter.key()];
3223
3224             // Check against the specified edges.
3225             if (edgeToCheck == check[0])
3226             {
3227                 meshOps::insertLabel
3228                 (
3229                     newPointIndex,
3230                     faces_[fIndex][1],
3231                     faces_[fIndex][2],
3232                     edgePoints_[eIter.key()]
3233                 );
3234
3235                 meshOps::sizeUpList(newFaceIndex[0], edgeFaces_[eIter.key()]);
3236                 newFaceEdges[0][nE[0]++] = eIter.key();
3237             }
3238
3239             if (edgeToCheck == check[1])
3240             {
3241                 meshOps::insertLabel
3242                 (
3243                     newPointIndex,
3244                     faces_[fIndex][0],
3245                     faces_[fIndex][2],
3246                     edgePoints_[eIter.key()]
3247                 );
3248
3249                 meshOps::sizeUpList(newFaceIndex[1], edgeFaces_[eIter.key()]);
3250                 newFaceEdges[1][nE[1]++] = eIter.key();
3251             }
3252
3253             if (edgeToCheck == check[2])
3254             {
3255                 meshOps::insertLabel
3256                 (
3257                     newPointIndex,
3258                     faces_[fIndex][0],
3259                     faces_[fIndex][1],
3260                     edgePoints_[eIter.key()]
3261                 );
3262
3263                 meshOps::sizeUpList(newFaceIndex[2], edgeFaces_[eIter.key()]);
3264                 newFaceEdges[2][nE[2]++] = eIter.key();
3265             }
3266
3267             if (edgeToCheck == check[3])
3268             {
3269                 meshOps::replaceLabel
3270                 (
3271                     faces_[fIndex][1],
3272                     newPointIndex,
3273                     edgePoints_[eIter.key()]
3274                 );
3275
3276                 meshOps::replaceLabel
3277                 (
3278                     fIndex,
3279                     newFaceIndex[3],
3280                     edgeFaces_[eIter.key()]
3281                 );
3282
3283                 newFaceEdges[3][nE[3]++] = eIter.key();
3284             }
3285
3286             if (edgeToCheck == check[4])
3287             {
3288                 meshOps::replaceLabel
3289                 (
3290                     faces_[fIndex][2],
3291                     newPointIndex,
3292                     edgePoints_[eIter.key()]
3293                 );
3294
3295                 meshOps::replaceLabel
3296                 (
3297                     fIndex,
3298                     newFaceIndex[4],
3299                     edgeFaces_[eIter.key()]
3300                 );
3301
3302                 newFaceEdges[4][nE[4]++] = eIter.key();
3303             }
3304
3305             if (edgeToCheck == check[5])
3306             {
3307                 meshOps::replaceLabel
3308                 (
3309                     faces_[fIndex][0],
3310                     newPointIndex,
3311                     edgePoints_[eIter.key()]
3312                 );
3313
3314                 meshOps::replaceLabel
3315                 (
3316                     fIndex,
3317                     newFaceIndex[5],
3318                     edgeFaces_[eIter.key()]
3319                 );
3320
3321                 newFaceEdges[5][nE[5]++] = eIter.key();
3322             }
3323         }
3324
3325         // Now that faceEdges has been configured, append them to the list.
3326         for (label i = 0; i < 6; i++)
3327         {
3328             faceEdges_.append(newFaceEdges[i]);
3329
3330             // Add faces to the map.
3331             map.addFace(newFaceIndex[i]);
3332         }
3333     }
3334     else
3335     {
3336         // Add three new cells to the end of the cell list
3337         for (label i = 3; i < 6; i++)
3338         {
3339             scalar parentScale = -1.0;
3340
3341             if (edgeRefinement_)
3342             {
3343                 parentScale = lengthScale_[cellsForRemoval[1]];
3344             }
3345
3346             newCellIndex[i] = insertCell(newTetCell[i], parentScale);
3347
3348             // Add to the map.
3349             map.addCell(newCellIndex[i]);
3350         }
3351
3352         // Find the apex point for this cell
3353         apexPoint[1] =
3354         (
3355             meshOps::tetApexPoint
3356             (
3357                 neighbour_[fIndex],
3358                 fIndex,
3359                 faces_,
3360                 cells_
3361             )
3362         );
3363
3364         // Add six new interior faces.
3365
3366         // Fourth face: Owner: newCellIndex[0], Neighbour: newCellIndex[3]
3367         tmpTriFace[0] = newPointIndex;
3368         tmpTriFace[1] = faces_[fIndex][2];
3369         tmpTriFace[2] = faces_[fIndex][0];
3370
3371         newFaceIndex[3] =
3372         (
3373             insertFace
3374             (
3375                 -1,
3376                 tmpTriFace,
3377                 newCellIndex[0],
3378                 newCellIndex[3]
3379             )
3380         );
3381
3382         // Fifth face: Owner: newCellIndex[1], Neighbour: newCellIndex[4]
3383         tmpTriFace[0] = newPointIndex;
3384         tmpTriFace[1] = faces_[fIndex][0];
3385         tmpTriFace[2] = faces_[fIndex][1];
3386
3387         newFaceIndex[4] =
3388         (
3389             insertFace
3390             (
3391                 -1,
3392                 tmpTriFace,
3393                 newCellIndex[1],
3394                 newCellIndex[4]
3395             )
3396         );
3397
3398         // Sixth face: Owner: newCellIndex[2], Neighbour: newCellIndex[5]
3399         tmpTriFace[0] = newPointIndex;
3400         tmpTriFace[1] = faces_[fIndex][1];
3401         tmpTriFace[2] = faces_[fIndex][2];
3402
3403         newFaceIndex[5] =
3404         (
3405             insertFace
3406             (
3407                 -1,
3408                 tmpTriFace,
3409                 newCellIndex[2],
3410                 newCellIndex[5]
3411             )
3412         );
3413
3414         // Seventh face: Owner: newCellIndex[3], Neighbour: newCellIndex[4]
3415         tmpTriFace[0] = newPointIndex;
3416         tmpTriFace[1] = apexPoint[1];
3417         tmpTriFace[2] = faces_[fIndex][0];
3418
3419         newFaceIndex[6] =
3420         (
3421             insertFace
3422             (
3423                 -1,
3424                 tmpTriFace,
3425                 newCellIndex[3],
3426                 newCellIndex[4]
3427             )
3428         );
3429
3430         // Eighth face: Owner: newCellIndex[4], Neighbour: newCellIndex[5]
3431         tmpTriFace[0] = newPointIndex;
3432         tmpTriFace[1] = apexPoint[1];
3433         tmpTriFace[2] = faces_[fIndex][1];
3434
3435         newFaceIndex[7] =
3436         (
3437             insertFace
3438             (
3439                 -1,
3440                 tmpTriFace,
3441                 newCellIndex[4],
3442                 newCellIndex[5]
3443             )
3444         );
3445
3446         // Ninth face: Owner: newCellIndex[3], Neighbour: newCellIndex[5]
3447         tmpTriFace[0] = newPointIndex;
3448         tmpTriFace[1] = faces_[fIndex][2];
3449         tmpTriFace[2] = apexPoint[1];
3450
3451         newFaceIndex[8] =
3452         (
3453             insertFace
3454             (
3455                 -1,
3456                 tmpTriFace,
3457                 newCellIndex[3],
3458                 newCellIndex[5]
3459             )
3460         );
3461
3462         // Add the newly created faces to cells
3463         newTetCell[3][nF[3]++] = newFaceIndex[6];
3464         newTetCell[3][nF[3]++] = newFaceIndex[8];
3465         newTetCell[4][nF[4]++] = newFaceIndex[6];
3466         newTetCell[4][nF[4]++] = newFaceIndex[7];
3467         newTetCell[5][nF[5]++] = newFaceIndex[7];
3468         newTetCell[5][nF[5]++] = newFaceIndex[8];
3469
3470         newTetCell[0][nF[0]++] = newFaceIndex[3];
3471         newTetCell[1][nF[1]++] = newFaceIndex[4];
3472         newTetCell[2][nF[2]++] = newFaceIndex[5];
3473
3474         newTetCell[3][nF[3]++] = newFaceIndex[3];
3475         newTetCell[4][nF[4]++] = newFaceIndex[4];
3476         newTetCell[5][nF[5]++] = newFaceIndex[5];
3477
3478         // Define the three faces to check for orientation:
3479         checkFace[0][0] = faces_[fIndex][2];
3480         checkFace[0][1] = apexPoint[1];
3481         checkFace[0][2] = faces_[fIndex][0];
3482
3483         checkFace[1][0] = faces_[fIndex][0];
3484         checkFace[1][1] = apexPoint[1];
3485         checkFace[1][2] = faces_[fIndex][1];
3486
3487         checkFace[2][0] = faces_[fIndex][1];
3488         checkFace[2][1] = apexPoint[1];
3489         checkFace[2][2] = faces_[fIndex][2];
3490
3491         // Check the orientation of faces on the second cell.
3492         forAll(cells_[neighbour_[fIndex]], faceI)
3493         {
3494             label faceIndex = cells_[neighbour_[fIndex]][faceI];
3495
3496             if (faceIndex == fIndex)
3497             {
3498                 continue;
3499             }
3500
3501             const face& faceToCheck = faces_[faceIndex];
3502             label cellIndex = cellsForRemoval[1];
3503             label newIndex = -1;
3504
3505             // Check against faces.
3506             if (triFace::compare(triFace(faceToCheck), triFace(checkFace[0])))
3507             {
3508                 newIndex = newCellIndex[3];
3509                 newTetCell[3][nF[3]++] = faceIndex;
3510             }
3511             else
3512             if (triFace::compare(triFace(faceToCheck), triFace(checkFace[1])))
3513             {
3514                 newIndex = newCellIndex[4];
3515                 newTetCell[4][nF[4]++] = faceIndex;
3516             }
3517             else
3518             if (triFace::compare(triFace(faceToCheck), triFace(checkFace[2])))
3519             {
3520                 newIndex = newCellIndex[5];
3521                 newTetCell[5][nF[5]++] = faceIndex;
3522             }
3523             else
3524             {
3525                 // Something's terribly wrong.
3526                 FatalErrorIn
3527                 (
3528                     "const changeMap dynamicTopoFvMesh::trisectFace\n"
3529                     "(\n"
3530                     "    const label fIndex,\n"
3531                     "    bool checkOnly,\n"
3532                     "    bool forceOp\n"
3533                     ")\n"
3534                 )
3535                     << "Failed to determine a face match."
3536                     << abort(FatalError);
3537             }
3538
3539             // Check if a face-flip is necessary
3540             if (owner_[faceIndex] == cellIndex)
3541             {
3542                 if (neighbour_[faceIndex] == -1)
3543                 {
3544                     // Change the owner
3545                     owner_[faceIndex] = newIndex;
3546                 }
3547                 else
3548                 {
3549                     // Flip this face
3550                     faces_[faceIndex] = faceToCheck.reverseFace();
3551                     owner_[faceIndex] = neighbour_[faceIndex];
3552                     neighbour_[faceIndex] = newIndex;
3553
3554                     setFlip(faceIndex);
3555                 }
3556             }
3557             else
3558             {
3559                 // Flip is unnecessary. Just update neighbour
3560                 neighbour_[faceIndex] = newIndex;
3561             }
3562         }
3563
3564         // Configure edgeFaces and edgePoints for four new interior edges.
3565         newQuadEdgeFaces[0] = newFaceIndex[4];
3566         newQuadEdgeFaces[1] = newFaceIndex[0];
3567         newQuadEdgeFaces[2] = newFaceIndex[3];
3568         newQuadEdgeFaces[3] = newFaceIndex[6];
3569
3570         newQuadEdgePoints[0] = faces_[fIndex][1];
3571         newQuadEdgePoints[1] = apexPoint[0];
3572         newQuadEdgePoints[2] = faces_[fIndex][2];
3573         newQuadEdgePoints[3] = apexPoint[1];
3574
3575         newEdgeIndex[1] =
3576         (
3577             insertEdge
3578             (
3579                 -1,
3580                 edge
3581                 (
3582                    newPointIndex,
3583                    faces_[fIndex][0]
3584                 ),
3585                 newQuadEdgeFaces,
3586                 newQuadEdgePoints
3587             )
3588         );
3589
3590         newQuadEdgeFaces[0] = newFaceIndex[5];
3591         newQuadEdgeFaces[1] = newFaceIndex[1];
3592         newQuadEdgeFaces[2] = newFaceIndex[4];
3593         newQuadEdgeFaces[3] = newFaceIndex[7];
3594
3595         newQuadEdgePoints[0] = faces_[fIndex][2];
3596         newQuadEdgePoints[1] = apexPoint[0];
3597         newQuadEdgePoints[2] = faces_[fIndex][0];
3598         newQuadEdgePoints[3] = apexPoint[1];
3599
3600         newEdgeIndex[2] =
3601         (
3602             insertEdge
3603             (
3604                 -1,
3605                 edge
3606                 (
3607                    newPointIndex,
3608                    faces_[fIndex][1]
3609                 ),
3610                 newQuadEdgeFaces,
3611                 newQuadEdgePoints
3612             )
3613         );
3614
3615         newQuadEdgeFaces[0] = newFaceIndex[3];
3616         newQuadEdgeFaces[1] = newFaceIndex[2];
3617         newQuadEdgeFaces[2] = newFaceIndex[5];
3618         newQuadEdgeFaces[3] = newFaceIndex[8];
3619
3620         newQuadEdgePoints[0] = faces_[fIndex][0];
3621         newQuadEdgePoints[1] = apexPoint[0];
3622         newQuadEdgePoints[2] = faces_[fIndex][1];
3623         newQuadEdgePoints[3] = apexPoint[1];
3624
3625         newEdgeIndex[3] =
3626         (
3627             insertEdge
3628             (
3629                 -1,
3630                 edge
3631                 (
3632                    newPointIndex,
3633                    faces_[fIndex][2]
3634                 ),
3635                 newQuadEdgeFaces,
3636                 newQuadEdgePoints
3637             )
3638         );
3639
3640         newTriEdgeFaces[0] = newFaceIndex[6];
3641         newTriEdgeFaces[1] = newFaceIndex[7];
3642         newTriEdgeFaces[2] = newFaceIndex[8];
3643
3644         newTriEdgePoints[0] = faces_[fIndex][0];
3645         newTriEdgePoints[1] = faces_[fIndex][1];
3646         newTriEdgePoints[2] = faces_[fIndex][2];
3647
3648         newEdgeIndex[4] =
3649         (
3650             insertEdge
3651             (
3652                 -1,
3653                 edge
3654                 (
3655                    apexPoint[1],
3656                    newPointIndex
3657                 ),
3658                 newTriEdgeFaces,
3659                 newTriEdgePoints
3660             )
3661         );
3662
3663         // Configure faceEdges with the new internal edges
3664         newFaceEdges[0][nE[0]++] = newEdgeIndex[1];
3665         newFaceEdges[1][nE[1]++] = newEdgeIndex[2];
3666         newFaceEdges[2][nE[2]++] = newEdgeIndex[3];
3667
3668         newFaceEdges[3][nE[3]++] = newEdgeIndex[1];
3669         newFaceEdges[3][nE[3]++] = newEdgeIndex[3];
3670         newFaceEdges[4][nE[4]++] = newEdgeIndex[1];
3671         newFaceEdges[4][nE[4]++] = newEdgeIndex[2];
3672         newFaceEdges[5][nE[5]++] = newEdgeIndex[2];
3673         newFaceEdges[5][nE[5]++] = newEdgeIndex[3];
3674
3675         newFaceEdges[6][nE[6]++] = newEdgeIndex[1];
3676         newFaceEdges[7][nE[7]++] = newEdgeIndex[2];
3677         newFaceEdges[8][nE[8]++] = newEdgeIndex[3];
3678
3679         newFaceEdges[6][nE[6]++] = newEdgeIndex[4];
3680         newFaceEdges[7][nE[7]++] = newEdgeIndex[4];
3681         newFaceEdges[8][nE[8]++] = newEdgeIndex[4];
3682
3683         // Define the nine edges to check while building faceEdges:
3684         FixedList<edge,9> check;
3685
3686         check[0][0] = apexPoint[0]; check[0][1] = faces_[fIndex][0];
3687         check[1][0] = apexPoint[0]; check[1][1] = faces_[fIndex][1];
3688         check[2][0] = apexPoint[0]; check[2][1] = faces_[fIndex][2];
3689
3690         check[3][0] = faces_[fIndex][2]; check[3][1] = faces_[fIndex][0];
3691         check[4][0] = faces_[fIndex][0]; check[4][1] = faces_[fIndex][1];
3692         check[5][0] = faces_[fIndex][1]; check[5][1] = faces_[fIndex][2];
3693
3694         check[6][0] = apexPoint[1]; check[6][1] = faces_[fIndex][0];
3695         check[7][0] = apexPoint[1]; check[7][1] = faces_[fIndex][1];
3696         check[8][0] = apexPoint[1]; check[8][1] = faces_[fIndex][2];
3697
3698         // Build a list of cellEdges
3699         labelHashSet cellEdges;
3700
3701         forAll(cells_[owner_[fIndex]], faceI)
3702         {
3703             const labelList& fEdges =
3704             (
3705                 faceEdges_[cells_[owner_[fIndex]][faceI]]
3706             );
3707
3708             forAll(fEdges, edgeI)
3709             {
3710                 if (!cellEdges.found(fEdges[edgeI]))
3711                 {
3712                     cellEdges.insert(fEdges[edgeI]);
3713                 }
3714             }
3715         }
3716
3717         forAll(cells_[neighbour_[fIndex]], faceI)
3718         {
3719             const labelList& fEdges =
3720             (
3721                 faceEdges_[cells_[neighbour_[fIndex]][faceI]]
3722             );
3723
3724             forAll(fEdges, edgeI)
3725             {
3726                 if (!cellEdges.found(fEdges[edgeI]))
3727                 {
3728                     cellEdges.insert(fEdges[edgeI]);
3729                 }
3730             }
3731         }
3732
3733         // Loop through cellEdges, and perform appropriate actions.
3734         forAllIter(labelHashSet, cellEdges, eIter)
3735         {
3736             const edge& edgeToCheck = edges_[eIter.key()];
3737
3738             // Check against the specified edges.
3739             if (edgeToCheck == check[0])
3740             {
3741                 meshOps::insertLabel
3742                 (
3743                     newPointIndex,
3744                     faces_[fIndex][1],
3745                     faces_[fIndex][2],
3746                     edgePoints_[eIter.key()]
3747                 );
3748
3749                 meshOps::sizeUpList(newFaceIndex[0], edgeFaces_[eIter.key()]);
3750                 newFaceEdges[0][nE[0]++] = eIter.key();
3751             }
3752
3753             if (edgeToCheck == check[1])
3754             {
3755                 meshOps::insertLabel
3756                 (
3757                     newPointIndex,
3758                     faces_[fIndex][0],
3759                     faces_[fIndex][2],
3760                     edgePoints_[eIter.key()]
3761                 );
3762
3763                 meshOps::sizeUpList(newFaceIndex[1], edgeFaces_[eIter.key()]);
3764                 newFaceEdges[1][nE[1]++] = eIter.key();
3765             }
3766
3767             if (edgeToCheck == check[2])
3768             {
3769                 meshOps::insertLabel
3770                 (
3771                     newPointIndex,
3772                     faces_[fIndex][0],
3773                     faces_[fIndex][1],
3774                     edgePoints_[eIter.key()]
3775                 );
3776
3777                 meshOps::sizeUpList(newFaceIndex[2], edgeFaces_[eIter.key()]);
3778                 newFaceEdges[2][nE[2]++] = eIter.key();
3779             }
3780
3781             if (edgeToCheck == check[3])
3782             {
3783                 meshOps::replaceLabel
3784                 (
3785                     faces_[fIndex][1],
3786                     newPointIndex,
3787                     edgePoints_[eIter.key()]
3788                 );
3789
3790                 meshOps::replaceLabel
3791                 (
3792                     fIndex,
3793                     newFaceIndex[3],
3794                     edgeFaces_[eIter.key()]
3795                 );
3796
3797                 newFaceEdges[3][nE[3]++] = eIter.key();
3798             }
3799
3800             if (edgeToCheck == check[4])
3801             {
3802                 meshOps::replaceLabel
3803                 (
3804                     faces_[fIndex][2],
3805                     newPointIndex,
3806                     edgePoints_[eIter.key()]
3807                 );
3808
3809                 meshOps::replaceLabel
3810                 (
3811                     fIndex,
3812                     newFaceIndex[4],
3813                     edgeFaces_[eIter.key()]
3814                 );
3815
3816                 newFaceEdges[4][nE[4]++] = eIter.key();
3817             }
3818
3819             if (edgeToCheck == check[5])
3820             {
3821                 meshOps::replaceLabel
3822                 (
3823                     faces_[fIndex][0],
3824                     newPointIndex,
3825                     edgePoints_[eIter.key()]
3826                 );
3827
3828                 meshOps::replaceLabel
3829                 (
3830                     fIndex,
3831                     newFaceIndex[5],
3832                     edgeFaces_[eIter.key()]
3833                 );
3834
3835                 newFaceEdges[5][nE[5]++] = eIter.key();
3836             }
3837
3838             if (edgeToCheck == check[6])
3839             {
3840                 meshOps::insertLabel
3841                 (
3842                     newPointIndex,
3843                     faces_[fIndex][1],
3844                     faces_[fIndex][2],
3845                     edgePoints_[eIter.key()]
3846                 );
3847
3848                 meshOps::sizeUpList(newFaceIndex[6], edgeFaces_[eIter.key()]);
3849                 newFaceEdges[6][nE[6]++] = eIter.key();
3850             }
3851
3852             if (edgeToCheck == check[7])
3853             {
3854                 meshOps::insertLabel
3855                 (
3856                     newPointIndex,
3857                     faces_[fIndex][0],
3858                     faces_[fIndex][2],
3859                     edgePoints_[eIter.key()]
3860                 );
3861
3862                 meshOps::sizeUpList(newFaceIndex[7], edgeFaces_[eIter.key()]);
3863                 newFaceEdges[7][nE[7]++] = eIter.key();
3864             }
3865
3866             if (edgeToCheck == check[8])
3867             {
3868                 meshOps::insertLabel
3869                 (
3870                     newPointIndex,
3871                     faces_[fIndex][0],
3872                     faces_[fIndex][1],
3873                     edgePoints_[eIter.key()]
3874                 );
3875
3876                 meshOps::sizeUpList(newFaceIndex[8], edgeFaces_[eIter.key()]);
3877                 newFaceEdges[8][nE[8]++] = eIter.key();
3878             }
3879         }
3880
3881         // Now that faceEdges has been configured, append them to the list.
3882         for (label i = 0; i < 9; i++)
3883         {
3884             faceEdges_.append(newFaceEdges[i]);
3885
3886             // Add faces to the map.
3887             map.addFace(newFaceIndex[i]);
3888         }
3889     }
3890
3891     // Added edges are those connected to the new point
3892     const labelList& pointEdges = pointEdges_[newPointIndex];
3893
3894     forAll(pointEdges, edgeI)
3895     {
3896         map.addEdge(pointEdges[edgeI]);
3897     }
3898
3899     // Now generate mapping info and remove entities.
3900     forAll(cellsForRemoval, cellI)
3901     {
3902         label cIndex = cellsForRemoval[cellI];
3903
3904         if (cIndex == -1)
3905         {
3906             continue;
3907         }
3908
3909         // Fill-in mapping information
3910         labelList mC(1, cellsForRemoval[cellI]);
3911
3912         if (cellI == 0)
3913         {
3914             for (label i = 0; i < 3; i++)
3915             {
3916                 // Update the cell list with newly configured cells.
3917                 cells_[newCellIndex[i]] = newTetCell[i];
3918
3919                 setCellMapping(newCellIndex[i], mC);
3920             }
3921         }
3922         else
3923         {
3924             for (label i = 3; i < 6; i++)
3925             {
3926                 // Update the cell list with newly configured cells.
3927                 cells_[newCellIndex[i]] = newTetCell[i];
3928
3929                 setCellMapping(newCellIndex[i], mC);
3930             }
3931         }
3932
3933         removeCell(cIndex);
3934     }
3935
3936     // Set default mapping for interior faces.
3937     for (label i = 0; i < 3; i++)
3938     {
3939         setFaceMapping(newFaceIndex[i]);
3940     }
3941
3942     if (cellsForRemoval[1] == -1)
3943     {
3944         // Set mapping for boundary faces.
3945         for (label i = 3; i < 6; i++)
3946         {
3947             setFaceMapping(newFaceIndex[i], labelList(1, fIndex));
3948         }
3949     }
3950     else
3951     {
3952         // Set default mapping for interior faces.
3953         for (label i = 3; i < 9; i++)
3954         {
3955             setFaceMapping(newFaceIndex[i]);
3956         }
3957     }
3958
3959     // Now finally remove the face...
3960     removeFace(fIndex);
3961
3962     if (debug > 2)
3963     {
3964         Info << "New Point:: " << newPointIndex << endl;
3965
3966         const labelList& pEdges = pointEdges_[newPointIndex];
3967
3968         Info << "pointEdges:: " << pEdges << endl;
3969
3970         Info << "Added edges: " << endl;
3971         forAll(pEdges, edgeI)
3972         {
3973             Info << pEdges[edgeI]
3974                  << ":: " << edges_[pEdges[edgeI]] << nl
3975                  << " edgeFaces:: " << edgeFaces_[pEdges[edgeI]] << nl
3976                  << " edgePoints:: " << edgePoints_[pEdges[edgeI]]
3977                  << endl;
3978         }
3979
3980         Info << "Added faces: " << endl;
3981         forAll(newFaceIndex, faceI)
3982         {
3983             if (newFaceIndex[faceI] == -1)
3984             {
3985                 continue;
3986             }
3987
3988             Info << newFaceIndex[faceI] << ":: "
3989                  << faces_[newFaceIndex[faceI]]
3990                  << endl;
3991         }
3992
3993         Info << "Added cells: " << endl;
3994         forAll(newCellIndex, cellI)
3995         {
3996             if (newCellIndex[cellI] == -1)
3997             {
3998                 continue;
3999             }
4000
4001             Info << newCellIndex[cellI] << ":: "
4002                  << cells_[newCellIndex[cellI]]
4003                  << endl;
4004         }
4005
4006         // Write out VTK files after change
4007         if (debug > 3)
4008         {
4009             labelList vtkCells;
4010
4011             if (cellsForRemoval[1] == -1)
4012             {
4013                 vtkCells.setSize(3);
4014
4015                 // Fill in cell indices
4016                 vtkCells[0] = newCellIndex[0];
4017                 vtkCells[1] = newCellIndex[1];
4018                 vtkCells[2] = newCellIndex[2];
4019             }
4020             else
4021             {
4022                 vtkCells.setSize(6);
4023
4024                 // Fill in cell indices
4025                 forAll(newCellIndex, indexI)
4026                 {
4027                     vtkCells[indexI] = newCellIndex[indexI];
4028                 }
4029             }
4030
4031             writeVTK
4032             (
4033                 Foam::name(fIndex)
4034               + "Trisect_1",
4035                 vtkCells
4036             );
4037         }
4038     }
4039
4040     // Set the flag
4041     topoChangeFlag_ = true;
4042
4043     // Increment the counter
4044     statistics_[3]++;
4045
4046     // Increment the number of modifications
4047     statistics_[0]++;
4048
4049     // Specify that the operation was successful
4050     map.type() = 1;
4051
4052     // Return the changeMap
4053     return map;
4054 }
4055
4056
4057 // Utility method to compute the quality of a
4058 // vertex hull around an edge after bisection.
4059 scalar dynamicTopoFvMesh::computeBisectionQuality
4060 (
4061     const label eIndex
4062 ) const
4063 {
4064     scalar minQuality = GREAT, minVolume = GREAT;
4065     scalar cQuality = 0.0, oldVolume = 0.0;
4066
4067     // Obtain a reference to this edge and corresponding edgePoints
4068     const edge& edgeToCheck = edges_[eIndex];
4069     const labelList& hullVertices = edgePoints_[eIndex];
4070
4071     // Obtain point references
4072     const point& a = points_[edgeToCheck[0]];
4073     const point& c = points_[edgeToCheck[1]];
4074
4075     const point& aOld = oldPoints_[edgeToCheck[0]];
4076     const point& cOld = oldPoints_[edgeToCheck[1]];
4077
4078     // Compute the mid-point of the edge
4079     point midPoint = 0.5*(a + c);
4080     point oldPoint = 0.5*(aOld + cOld);
4081
4082     if (whichEdgePatch(eIndex) < 0)
4083     {
4084         // Internal edge.
4085         forAll(hullVertices, indexI)
4086         {
4087             label prevIndex = hullVertices.rcIndex(indexI);
4088
4089             // Pick vertices off the list
4090             const point& b = points_[hullVertices[prevIndex]];
4091             const point& d = points_[hullVertices[indexI]];
4092
4093             const point& bOld = oldPoints_[hullVertices[prevIndex]];
4094             const point& dOld = oldPoints_[hullVertices[indexI]];
4095
4096             // Compute the quality of the upper half.
4097             cQuality = tetMetric_(a, b, midPoint, d);
4098
4099             // Compute old volume of the upper half.
4100             oldVolume = tetPointRef(aOld, bOld, oldPoint, dOld).mag();
4101
4102             // Check if the volume / quality is worse
4103             minQuality = Foam::min(cQuality, minQuality);
4104             minVolume = Foam::min(oldVolume, minVolume);
4105
4106             // Compute the quality of the lower half.
4107             cQuality = tetMetric_(midPoint, b, c, d);
4108
4109             // Compute old volume of the lower half.
4110             oldVolume = tetPointRef(oldPoint, bOld, cOld, dOld).mag();
4111
4112             // Check if the quality is worse
4113             minQuality = Foam::min(cQuality, minQuality);
4114             minVolume = Foam::min(oldVolume, minVolume);
4115         }
4116     }
4117     else
4118     {
4119         // Boundary edge.
4120         for(label indexI = 1; indexI < hullVertices.size(); indexI++)
4121         {
4122             // Pick vertices off the list
4123             const point& b = points_[hullVertices[indexI-1]];
4124             const point& d = points_[hullVertices[indexI]];
4125
4126             const point& bOld = oldPoints_[hullVertices[indexI-1]];
4127             const point& dOld = oldPoints_[hullVertices[indexI]];
4128
4129             // Compute the quality of the upper half.
4130             cQuality = tetMetric_(a, b, midPoint, d);
4131
4132             // Compute old volume of the upper half.
4133             oldVolume = tetPointRef(aOld, bOld, oldPoint, dOld).mag();
4134
4135             // Check if the quality is worse
4136             minQuality = Foam::min(cQuality, minQuality);
4137             minVolume = Foam::min(oldVolume, minVolume);
4138
4139             // Compute the quality of the lower half.
4140             cQuality = tetMetric_(midPoint, b, c, d);
4141
4142             // Compute old volume of the lower half.
4143             oldVolume = tetPointRef(oldPoint, bOld, cOld, dOld).mag();
4144
4145             // Check if the quality is worse
4146             minQuality = Foam::min(cQuality, minQuality);
4147             minVolume = Foam::min(oldVolume, minVolume);
4148         }
4149     }
4150
4151     // Ensure that the mesh is valid
4152     if (minQuality < sliverThreshold_)
4153     {
4154         if (debug > 3 && minQuality < 0.0)
4155         {
4156             // Write out cells for post processing.
4157             labelHashSet iCells;
4158
4159             const labelList& eFaces = edgeFaces_[eIndex];
4160
4161             forAll(eFaces, faceI)
4162             {
4163                 if (!iCells.found(owner_[eFaces[faceI]]))
4164                 {
4165                     iCells.insert(owner_[eFaces[faceI]]);
4166                 }
4167
4168                 if (!iCells.found(neighbour_[eFaces[faceI]]))
4169                 {
4170                     iCells.insert(neighbour_[eFaces[faceI]]);
4171                 }
4172             }
4173
4174             writeVTK(Foam::name(eIndex) + "_iCells", iCells.toc());
4175         }
4176
4177         if (debug > 2)
4178         {
4179             InfoIn
4180             (
4181                 "scalar dynamicTopoFvMesh::computeBisectionQuality"
4182                 "(const label eIndex) const"
4183             )
4184                 << "Bisecting edge will fall below the "
4185                 << "sliver threshold of: " << sliverThreshold_ << nl
4186                 << "Edge: " << eIndex << ": " << edgeToCheck << nl
4187                 << "EdgePoints: " << hullVertices << nl
4188                 << "Minimum Quality: " << minQuality << nl
4189                 << "Mid point: " << midPoint
4190                 << endl;
4191         }
4192     }
4193
4194     // If a negative old-volume was encountered,
4195     // return an invalid quality.
4196     if (minVolume < 0.0)
4197     {
4198         return minVolume;
4199     }
4200
4201     return minQuality;
4202 }
4203
4204
4205 // Utility method to compute the quality of cells
4206 // around a face after trisection.
4207 scalar dynamicTopoFvMesh::computeTrisectionQuality
4208 (
4209     const label fIndex
4210 ) const
4211 {
4212     scalar minQuality = GREAT;
4213     scalar cQuality = 0.0;
4214
4215     point midPoint;
4216
4217     // Fetch the midPoint
4218     midPoint = faces_[fIndex].centre(points_);
4219
4220     FixedList<label,2> apexPoint(-1);
4221
4222     // Find the apex point
4223     apexPoint[0] =
4224     (
4225         meshOps::tetApexPoint
4226         (
4227             owner_[fIndex],
4228             fIndex,
4229             faces_,
4230             cells_
4231         )
4232     );
4233
4234     const face& faceToCheck = faces_[fIndex];
4235
4236     forAll(faceToCheck, pointI)
4237     {
4238         // Pick vertices off the list
4239         const point& b = points_[faceToCheck[pointI]];
4240         const point& c = points_[apexPoint[0]];
4241         const point& d = points_[faceToCheck[faceToCheck.fcIndex(pointI)]];
4242
4243         // Compute the quality of the upper half.
4244         cQuality = tetMetric_(midPoint, b, c, d);
4245
4246         // Check if the quality is worse
4247         minQuality = Foam::min(cQuality, minQuality);
4248     }
4249
4250     if (whichPatch(fIndex) == -1)
4251     {
4252         apexPoint[1] =
4253         (
4254             meshOps::tetApexPoint
4255             (
4256                 neighbour_[fIndex],
4257                 fIndex,
4258                 faces_,
4259                 cells_
4260             )
4261         );
4262
4263         forAll(faceToCheck, pointI)
4264         {
4265             // Pick vertices off the list
4266             const point& b = points_[faceToCheck[pointI]];
4267             const point& c = points_[apexPoint[1]];
4268             const point& d = points_[faceToCheck[faceToCheck.rcIndex(pointI)]];
4269
4270             // Compute the quality of the upper half.
4271             cQuality = tetMetric_(midPoint, b, c, d);
4272
4273             // Check if the quality is worse
4274             minQuality = Foam::min(cQuality, minQuality);
4275         }
4276     }
4277
4278     return minQuality;
4279 }
4280
4281
4282 // Split a set of internal faces into boundary faces
4283 //   - Add boundary faces and edges to the patch specified by 'patchIndex'
4284 //   - Cell color should specify a binary value dictating either side
4285 //     of the split face.
4286 void dynamicTopoFvMesh::splitInternalFaces
4287 (
4288     const label patchIndex,
4289     const labelList& internalFaces,
4290     const Map<bool>& cellColors
4291 )
4292 {
4293     Map<label> mirrorPointLabels;
4294     FixedList<Map<label>, 2> mirrorEdgeLabels, mirrorFaceLabels;
4295
4296     // First loop through the list and accumulate a list of
4297     // points and edges that need to be duplicated.
4298     forAll(internalFaces, faceI)
4299     {
4300         const face& faceToCheck = faces_[internalFaces[faceI]];
4301
4302         forAll(faceToCheck, pointI)
4303         {
4304             if (!mirrorPointLabels.found(faceToCheck[pointI]))
4305             {
4306                 mirrorPointLabels.insert(faceToCheck[pointI], -1);
4307             }
4308         }
4309
4310         const labelList& fEdges = faceEdges_[internalFaces[faceI]];
4311
4312         forAll(fEdges, edgeI)
4313         {
4314             if (!mirrorEdgeLabels[0].found(fEdges[edgeI]))
4315             {
4316                 mirrorEdgeLabels[0].insert(fEdges[edgeI], -1);
4317             }
4318         }
4319     }
4320
4321     // Now for every point in the list, add a new one.
4322     // Add a mapping entry as well.
4323     forAllIter(Map<label>, mirrorPointLabels, pIter)
4324     {
4325         // Obtain a copy of the point before adding it,
4326         // since the reference might become invalid during list resizing.
4327         point newPoint = points_[pIter.key()];
4328         point oldPoint = oldPoints_[pIter.key()];
4329
4330         pIter() = insertPoint(newPoint, oldPoint, labelList(1, pIter.key()));
4331
4332         if (!twoDMesh_)
4333         {
4334             const labelList& pEdges = pointEdges_[pIter.key()];
4335
4336             labelHashSet edgesToRemove;
4337
4338             forAll(pEdges, edgeI)
4339             {
4340                 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
4341
4342                 bool allTrue = true;
4343
4344                 forAll(eFaces, faceI)
4345                 {
4346                     label own = owner_[eFaces[faceI]];
4347                     label nei = neighbour_[eFaces[faceI]];
4348
4349                     // Check if an owner/neighbour cell is false
4350                     if (!cellColors[own])
4351                     {
4352                         allTrue = false;
4353                         break;
4354                     }
4355
4356                     if (nei != -1)
4357                     {
4358                         if (!cellColors[nei])
4359                         {
4360                             allTrue = false;
4361                             break;
4362                         }
4363                     }
4364                 }
4365
4366                 if (allTrue)
4367                 {
4368                     // Mark this edge label to be discarded later
4369                     edgesToRemove.insert(pEdges[edgeI]);
4370                 }
4371             }
4372
4373             // It is dangerous to use the pointEdges references,
4374             // so call it using array-lookup instead.
4375             forAllIter(labelHashSet, edgesToRemove, hsIter)
4376             {
4377                 // Add the edge to the mirror point list
4378                 meshOps::sizeUpList
4379                 (
4380                     hsIter.key(),
4381                     pointEdges_[pIter()]
4382                 );
4383
4384                 // Remove the edge from the original point list
4385                 meshOps::sizeDownList
4386                 (
4387                     hsIter.key(),
4388                     pointEdges_[pIter.key()]
4389                 );
4390             }
4391         }
4392     }
4393
4394     if (debug > 3)
4395     {
4396         label i = 0;
4397         labelList mPoints(mirrorPointLabels.size());
4398
4399         if (!twoDMesh_)
4400         {
4401             forAllIter(Map<label>, mirrorPointLabels, pIter)
4402             {
4403                 writeVTK
4404                 (
4405                     "pEdges_o_" + Foam::name(pIter.key()) + '_',
4406                     pointEdges_[pIter.key()],
4407                     1
4408                 );
4409
4410                 writeVTK
4411                 (
4412                     "pEdges_m_" + Foam::name(pIter()) + '_',
4413                     pointEdges_[pIter()],
4414                     1
4415                 );
4416
4417                 mPoints[i++] = pIter();
4418             }
4419
4420             writeVTK
4421             (
4422                 "points_o_",
4423                 mirrorPointLabels.toc(),
4424                 0
4425             );
4426
4427             writeVTK
4428             (
4429                 "points_m_",
4430                 mPoints,
4431                 0
4432             );
4433         }
4434     }
4435
4436     // For every internal face, add a new one.
4437     //  - Stick to the rule:
4438     //    [1] Every cell marked false keeps the existing entities.
4439     //    [2] Every cell marked true gets new points/edges/faces.
4440     //  - If faces are improperly oriented, reverse them.
4441     forAll(internalFaces, faceI)
4442     {
4443         FixedList<face, 2> newFace;
4444         FixedList<label, 2> newFaceIndex(-1);
4445         FixedList<label, 2> newOwner(-1);
4446
4447         label oldOwn = owner_[internalFaces[faceI]];
4448         label oldNei = neighbour_[internalFaces[faceI]];
4449
4450         if (cellColors[oldOwn] && !cellColors[oldNei])
4451         {
4452             // The owner gets a new boundary face.
4453             // Note that orientation is already correct.
4454             newFace[0] = faces_[internalFaces[faceI]];
4455
4456             // The neighbour needs to have its face reversed
4457             // and moved to the boundary patch, thereby getting
4458             // deleted in the process.
4459             newFace[1] = newFace[0].reverseFace();
4460
4461             newOwner[0] = oldOwn;
4462             newOwner[1] = oldNei;
4463         }
4464         else
4465         if (!cellColors[oldOwn] && cellColors[oldNei])
4466         {
4467             // The neighbour gets a new boundary face.
4468             // The face is oriented in the opposite sense, however.
4469             newFace[0] = faces_[internalFaces[faceI]].reverseFace();
4470
4471             // The owner keeps the existing face and orientation.
4472             // But it also needs to be moved to the boundary.
4473             newFace[1] = faces_[internalFaces[faceI]];
4474
4475             newOwner[0] = oldNei;
4476             newOwner[1] = oldOwn;
4477         }
4478         else
4479         {
4480             // Something's wrong here.
4481             FatalErrorIn
4482             (
4483                 "dynamicTopoFvMesh::splitInternalFaces\n"
4484                 "(\n"
4485                 "    const label patchIndex,\n"
4486                 "    const labelList& internalFaces,\n"
4487                 "    const Map<bool>& cellColors\n"
4488                 ")\n"
4489             )
4490                 << nl << " Face: "
4491                 << internalFaces[faceI]
4492                 << " has cells which are improperly marked: " << nl
4493                 << oldOwn << ":: " << cellColors[oldOwn] << nl
4494                 << oldNei << ":: " << cellColors[oldNei]
4495                 << abort(FatalError);
4496         }
4497
4498         // Renumber point labels for the first new face.
4499         forAll(newFace[0], pointI)
4500         {
4501             newFace[0][pointI] = mirrorPointLabels[newFace[0][pointI]];
4502         }
4503
4504         // Insert the new boundary faces.
4505         forAll(newFace, indexI)
4506         {
4507             newFaceIndex[indexI] =
4508             (
4509                 insertFace
4510                 (
4511                     patchIndex,
4512                     newFace[indexI],
4513                     newOwner[indexI],
4514                     -1
4515                 )
4516             );
4517
4518             // Make an identical faceEdges entry.
4519             // This will be renumbered once new edges are added.
4520             labelList newFaceEdges(faceEdges_[internalFaces[faceI]]);
4521
4522             faceEdges_.append(newFaceEdges);
4523
4524             // Replace face labels on cells
4525             meshOps::replaceLabel
4526             (
4527                 internalFaces[faceI],
4528                 newFaceIndex[indexI],
4529                 cells_[newOwner[indexI]]
4530             );
4531
4532             // Make face mapping entries for posterity.
4533             mirrorFaceLabels[indexI].insert
4534             (
4535                 internalFaces[faceI],
4536                 newFaceIndex[indexI]
4537             );
4538         }
4539     }
4540
4541     // For every edge in the list, add a new one.
4542     // We'll deal with correcting edgeFaces and edgePoints later.
4543     forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
4544     {
4545         // Obtain copies for the append method
4546         edge origEdge = edges_[eIter.key()];
4547         labelList newEdgeFaces(edgeFaces_[eIter.key()]);
4548
4549         eIter() =
4550         (
4551             insertEdge
4552             (
4553                 patchIndex,
4554                 edge
4555                 (
4556                     mirrorPointLabels[origEdge[0]],
4557                     mirrorPointLabels[origEdge[1]]
4558                 ),
4559                 newEdgeFaces,
4560                 labelList(0)
4561             )
4562         );
4563
4564         // Is the original edge an internal one?
4565         // If it is, we need to move it to the boundary.
4566         if (whichEdgePatch(eIter.key()) == -1)
4567         {
4568             label rplEdgeIndex =
4569             (
4570                 insertEdge
4571                 (
4572                     patchIndex,
4573                     origEdge,
4574                     newEdgeFaces,
4575                     labelList(0)
4576                 )
4577             );
4578
4579             // Map the new entry.
4580             mirrorEdgeLabels[1].insert(eIter.key(), rplEdgeIndex);
4581         }
4582         else
4583         {
4584             // This is already a boundary edge.
4585             // Make an identical map.
4586             mirrorEdgeLabels[1].insert(eIter.key(), eIter.key());
4587         }
4588     }
4589
4590     // Correct edgeFaces for all new edges
4591     forAll(mirrorEdgeLabels, indexI)
4592     {
4593         forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
4594         {
4595             labelList& eFaces = edgeFaces_[eIter()];
4596
4597             labelHashSet facesToRemove;
4598
4599             forAll(eFaces, faceI)
4600             {
4601                 bool remove = false;
4602
4603                 label own = owner_[eFaces[faceI]];
4604                 label nei = neighbour_[eFaces[faceI]];
4605
4606                 if
4607                 (
4608                     (!cellColors[own] && !indexI) ||
4609                     ( cellColors[own] &&  indexI)
4610                 )
4611                 {
4612                     remove = true;
4613                 }
4614
4615                 if (nei != -1)
4616                 {
4617                     if
4618                     (
4619                         (!cellColors[nei] && !indexI) ||
4620                         ( cellColors[nei] &&  indexI)
4621                     )
4622                     {
4623                         remove = true;
4624                     }
4625                 }
4626
4627                 if (mirrorFaceLabels[indexI].found(eFaces[faceI]))
4628                 {
4629                     // Perform a replacement instead of a removal.
4630                     eFaces[faceI] = mirrorFaceLabels[indexI][eFaces[faceI]];
4631
4632                     remove = false;
4633                 }
4634
4635                 if (remove)
4636                 {
4637                     facesToRemove.insert(eFaces[faceI]);
4638                 }
4639             }
4640
4641             // Now successively size down edgeFaces.
4642             // It is dangerous to use the eFaces reference,
4643             // so call it using array-lookup.
4644             forAllIter(labelHashSet, facesToRemove, hsIter)
4645             {
4646                 meshOps::sizeDownList(hsIter.key(), edgeFaces_[eIter()]);
4647             }
4648         }
4649     }
4650
4651     // Renumber faceEdges for all faces connected to new edges
4652     forAll(mirrorEdgeLabels, indexI)
4653     {
4654         forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
4655         {
4656             const labelList& eFaces = edgeFaces_[eIter()];
4657
4658             forAll(eFaces, faceI)
4659             {
4660                 labelList& fEdges = faceEdges_[eFaces[faceI]];
4661
4662                 forAll(fEdges, edgeI)
4663                 {
4664                     if (mirrorEdgeLabels[indexI].found(fEdges[edgeI]))
4665                     {
4666                         fEdges[edgeI] =
4667                         (
4668                             mirrorEdgeLabels[indexI][fEdges[edgeI]]
4669                         );
4670                     }
4671                 }
4672             }
4673         }
4674     }
4675
4676     if (twoDMesh_)
4677     {
4678         // Renumber edges and faces
4679         forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
4680         {
4681             const labelList& eFaces = edgeFaces_[eIter()];
4682
4683             // Two levels of indirection to ensure
4684             // that all entities we renumbered.
4685             // A flip-side for the lack of a pointEdges list in 2D.
4686             forAll(eFaces, faceI)
4687             {
4688                 const labelList& fEdges = faceEdges_[eFaces[faceI]];
4689
4690                 forAll(fEdges, edgeI)
4691                 {
4692                     // Renumber this edge.
4693                     edge& edgeToCheck = edges_[fEdges[edgeI]];
4694
4695                     forAll(edgeToCheck, pointI)
4696                     {
4697                         if (mirrorPointLabels.found(edgeToCheck[pointI]))
4698                         {
4699                             edgeToCheck[pointI] =
4700                             (
4701                                 mirrorPointLabels[edgeToCheck[pointI]]
4702                             );
4703                         }
4704                     }
4705
4706                     // Also renumber faces connected to this edge.
4707                     const labelList& efFaces = edgeFaces_[fEdges[edgeI]];
4708
4709                     forAll(efFaces, faceJ)
4710                     {
4711                         face& faceToCheck = faces_[efFaces[faceJ]];
4712
4713                         forAll(faceToCheck, pointI)
4714                         {
4715                             if
4716                             (
4717                                 mirrorPointLabels.found(faceToCheck[pointI])
4718                             )
4719                             {
4720                                 faceToCheck[pointI] =
4721                                 (
4722                                     mirrorPointLabels[faceToCheck[pointI]]
4723                                 );
4724                             }
4725                         }
4726                     }
4727                 }
4728             }
4729         }
4730     }
4731     else
4732     {
4733         // Point renumbering of entities connected to mirror points
4734         forAllIter(Map<label>, mirrorPointLabels, pIter)
4735         {
4736             const labelList& pEdges = pointEdges_[pIter()];
4737
4738             forAll(pEdges, edgeI)
4739             {
4740                 // Renumber this edge.
4741                 edge& edgeToCheck = edges_[pEdges[edgeI]];
4742
4743                 forAll(edgeToCheck, pointI)
4744                 {
4745                     if (edgeToCheck[pointI] == pIter.key())
4746                     {
4747                         edgeToCheck[pointI] = pIter();
4748                     }
4749                 }
4750
4751                 // Also renumber faces connected to this edge.
4752                 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
4753
4754                 forAll(eFaces, faceI)
4755                 {
4756                     face& faceToCheck = faces_[eFaces[faceI]];
4757
4758                     forAll(faceToCheck, pointI)
4759                     {
4760                         if (faceToCheck[pointI] == pIter.key())
4761                         {
4762                             faceToCheck[pointI] = pIter();
4763                         }
4764                     }
4765                 }
4766             }
4767         }
4768
4769         // Scan edges connected to mirror points,
4770         // and correct any edgePoints entries for
4771         // edges connected to the other vertex.
4772         forAllIter(Map<label>, mirrorPointLabels, pIter)
4773         {
4774             const labelList& pEdges = pointEdges_[pIter()];
4775
4776             forAll(pEdges, edgeI)
4777             {
4778                 label otherVertex = edges_[pEdges[edgeI]].otherVertex(pIter());
4779
4780                 // Scan edgePoints for edges connected to this point
4781                 const labelList& opEdges = pointEdges_[otherVertex];
4782
4783                 forAll(opEdges, edgeJ)
4784                 {
4785                     labelList& ePoints = edgePoints_[opEdges[edgeJ]];
4786
4787                     forAll(ePoints, pointI)
4788                     {
4789                         if (mirrorPointLabels.found(ePoints[pointI]))
4790                         {
4791                             // Replace this point with the mirror point
4792                             ePoints[pointI] =
4793                             (
4794                                 mirrorPointLabels[ePoints[pointI]]
4795                             );
4796                         }
4797                     }
4798                 }
4799             }
4800         }
4801
4802         // Build edgePoints for new edges
4803         forAll(mirrorEdgeLabels, indexI)
4804         {
4805             forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
4806             {
4807                 buildEdgePoints(eIter());
4808             }
4809         }
4810     }
4811
4812     // Now that we're done with the internal faces, remove them.
4813     forAll(internalFaces, faceI)
4814     {
4815         removeFace(internalFaces[faceI]);
4816     }
4817
4818     // Remove old internal edges as well.
4819     forAllIter(Map<label>, mirrorEdgeLabels[1], eIter)
4820     {
4821         if (eIter.key() != eIter())
4822         {
4823             removeEdge(eIter.key());
4824         }
4825     }
4826
4827     // Set the flag
4828     topoChangeFlag_ = true;
4829 }
4830
4831 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
4832
4833 } // End namespace Foam
4834
4835 // ************************************************************************* //