src/dynamicMesh/dynamicTopoFvMesh/edgeBisect.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "Stack.H"
  27 #include "triFace.H"
  28 #include "objectMap.H"
  29 #include "changeMap.H"
  30 #include "coupledInfo.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 //    -3: Bisection failed since edge was on a noRefinement patch.
  47 const changeMap dynamicTopoFvMesh::bisectQuadFace
  48 (
  49     const label fIndex,
  50     const changeMap& masterMap,
  51     bool checkOnly,
  52     bool forceOp
  53 )
  54 {
  55     // Quad-face bisection performs the following operations:
  56     //      [1] Add two points at the middle of the face
  57     //      [2] Create a new internal face for each bisected cell
  58     //      [3] Modify existing face and create a new half-face
  59     //      [4] Modify triangular boundary faces, and create new ones as well
  60     //      [5] Create edges for new faces
  61     //      Update faceEdges and edgeFaces information
  62
  63     // Figure out which thread this is...
  64     label tIndex = self();
  65
  66     // Prepare the changeMaps
  67     changeMap map;
  68     List<changeMap> slaveMaps;
  69     bool bisectingSlave = false;
  70
  71     if
  72     (
  73         (status(TOTAL_MODIFICATIONS) > maxModifications_) &&
  74         (maxModifications_ > -1)
  75     )
  76     {
  77         // Reached the max allowable topo-changes.
  78         stack(tIndex).clear();
  79
  80         return map;
  81     }
  82
  83     // Check if edgeRefinements are to be avoided on patch.
  84     if (baseMesh_.lengthEstimator().checkRefinementPatch(whichPatch(fIndex)))
  85     {
  86         map.type() = -3;
  87
  88         return map;
  89     }
  90
  91     // Sanity check: Is the index legitimate?
  92     if (fIndex < 0)
  93     {
  94         FatalErrorIn
  95         (
  96             "\n"
  97             "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
  98             "(\n"
  99             "    const label fIndex,\n"
 100             "    const changeMap& masterMap,\n"
 101             "    bool checkOnly,\n"
 102             "    bool forceOp\n"
 103             ")\n"
 104         )
 105             << " Invalid index: " << fIndex << nl
 106             << " nFaces: " << nFaces_
 107             << abort(FatalError);
 108     }
 109
 110     bool found;
 111     label replaceFace = -1, retainFace = -1;
 112     face tmpQuadFace(4), tmpTriFace(3);
 113     labelList tmpQFEdges(4, -1), tmpTFEdges(3, -1);
 114     FixedList<label,7> newFaceIndex(-1), newEdgeIndex(-1);
 115     FixedList<edge,4> commonEdges(edge(-1, -1));
 116     FixedList<label,4> cornerEdgeIndex(-1), commonEdgeIndex(-1);
 117     FixedList<label,4> commonFaceIndex(-1);
 118     FixedList<label,2> newPointIndex(-1), newCellIndex(-1);
 119     FixedList<label,4> otherEdgeIndex(-1), otherEdgePoint(-1);
 120     FixedList<label,4> otherPointIndex(-1), nextToOtherPoint(-1);
 121     FixedList<label,2> c0BdyIndex(-1), c0IntIndex(-1);
 122     FixedList<label,2> c1BdyIndex(-1), c1IntIndex(-1);
 123     FixedList<face,2> c0BdyFace(face(3)), c0IntFace(face(4));
 124     FixedList<face,2> c1BdyFace(face(3)), c1IntFace(face(4));
 125
 126     // Get the two cells on either side...
 127     label c0 = owner_[fIndex], c1 = neighbour_[fIndex];
 128
 129     // Keep track of old / new cells
 130     FixedList<cell, 2> oldCells(cell(5));
 131     FixedList<cell, 2> newCells(cell(5));
 132
 133     // Find the prism faces for cell[0].
 134     oldCells[0] = cells_[c0];
 135
 136     meshOps::findPrismFaces
 137     (
 138         fIndex,
 139         c0,
 140         faces_,
 141         cells_,
 142         neighbour_,
 143         c0BdyFace,
 144         c0BdyIndex,
 145         c0IntFace,
 146         c0IntIndex
 147     );
 148
 149     // Check for resulting quality
 150     if (checkBisection(fIndex, c0BdyIndex[0], forceOp))
 151     {
 152         map.type() = -2;
 153         return map;
 154     }
 155
 156     if (c1 != -1)
 157     {
 158         // Find the prism faces for cell[1].
 159         meshOps::findPrismFaces
 160         (
 161             fIndex,
 162             c1,
 163             faces_,
 164             cells_,
 165             neighbour_,
 166             c1BdyFace,
 167             c1BdyIndex,
 168             c1IntFace,
 169             c1IntIndex
 170         );
 171
 172         // Check for resulting quality
 173         if (checkBisection(fIndex, c1BdyIndex[0], forceOp))
 174         {
 175             map.type() = -2;
 176             return map;
 177         }
 178     }
 179
 180     // Find the common-edge between the triangular boundary faces
 181     // and the face under consideration.
 182     meshOps::findCommonEdge
 183     (
 184         c0BdyIndex[0],
 185         fIndex,
 186         faceEdges_,
 187         commonEdgeIndex[0]
 188     );
 189
 190     meshOps::findCommonEdge
 191     (
 192         c0BdyIndex[1],
 193         fIndex,
 194         faceEdges_,
 195         commonEdgeIndex[1]
 196     );
 197
 198     commonEdges[0] = edges_[commonEdgeIndex[0]];
 199     commonEdges[1] = edges_[commonEdgeIndex[1]];
 200
 201     // If coupled modification is set, and this is a
 202     // master edge, bisect its slaves as well.
 203     bool localCouple = false, procCouple = false;
 204
 205     if (coupledModification_)
 206     {
 207         const label faceEnum = coupleMap::FACE;
 208         const label pointEnum = coupleMap::POINT;
 209
 210         // Is this a locally coupled face (either master or slave)?
 211         if (locallyCoupledEntity(fIndex, true))
 212         {
 213             localCouple = true;
 214             procCouple = false;
 215         }
 216         else
 217         if (processorCoupledEntity(fIndex))
 218         {
 219             procCouple = true;
 220             localCouple = false;
 221         }
 222
 223         if (localCouple && !procCouple)
 224         {
 225             // Determine the slave index.
 226             label sIndex = -1, pIndex = -1;
 227
 228             forAll(patchCoupling_, patchI)
 229             {
 230                 if (patchCoupling_(patchI))
 231                 {
 232                     const coupleMap& cMap = patchCoupling_[patchI].map();
 233
 234                     if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
 235                     {
 236                         pIndex = patchI;
 237
 238                         break;
 239                     }
 240
 241                     // The following bit happens only during the sliver
 242                     // exudation process, since slave faces are
 243                     // usually not added to the coupled face-stack.
 244                     if ((sIndex = cMap.findMaster(faceEnum, fIndex)) > -1)
 245                     {
 246                         pIndex = patchI;
 247
 248                         // Notice that we are collapsing a slave face first.
 249                         bisectingSlave = true;
 250
 251                         break;
 252                     }
 253                 }
 254             }
 255
 256             if (sIndex == -1)
 257             {
 258                 FatalErrorIn
 259                 (
 260                     "\n"
 261                     "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
 262                     "(\n"
 263                     "    const label fIndex,\n"
 264                     "    const changeMap& masterMap,\n"
 265                     "    bool checkOnly,\n"
 266                     "    bool forceOp"
 267                     ")\n"
 268                 )
 269                     << "Coupled maps were improperly specified." << nl
 270                     << " Slave index not found for: " << nl
 271                     << " Face: " << fIndex << ": " << faces_[fIndex] << nl
 272                     << abort(FatalError);
 273             }
 274             else
 275             {
 276                 // If we've found the slave, size up the list
 277                 meshOps::sizeUpList
 278                 (
 279                     changeMap(),
 280                     slaveMaps
 281                 );
 282
 283                 // Save index and patch for posterity
 284                 slaveMaps[0].index() = sIndex;
 285                 slaveMaps[0].patchIndex() = pIndex;
 286             }
 287
 288             if (debug > 1)
 289             {
 290                 Pout<< nl << " >> Bisecting slave face: " << sIndex
 291                     << " for master face: " << fIndex << endl;
 292             }
 293         }
 294         else
 295         if (procCouple && !localCouple)
 296         {
 297             // If this is a new entity, bail out for now.
 298             // This will be handled at the next time-step.
 299             if (fIndex >= nOldFaces_)
 300             {
 301                 return map;
 302             }
 303
 304             // Check slaves
 305             forAll(procIndices_, pI)
 306             {
 307                 // Fetch reference to subMesh
 308                 const coupledMesh& recvMesh = recvMeshes_[pI];
 309                 const coupleMap& cMap = recvMesh.map();
 310
 311                 label sIndex = -1;
 312
 313                 if ((sIndex = cMap.findSlave(faceEnum, fIndex)) > -1)
 314                 {
 315                     if (debug > 3)
 316                     {
 317                         Pout<< "Checking slave face: " << sIndex
 318                             << " on proc: " << procIndices_[pI]
 319                             << " for master face: " << fIndex
 320                             << endl;
 321                     }
 322
 323                     // Check if a lower-ranked processor is
 324                     // handling this edge
 325                     if
 326                     (
 327                         priority
 328                         (
 329                             procIndices_[pI],
 330                             lessOp<label>(),
 331                             Pstream::myProcNo()
 332                         )
 333                     )
 334                     {
 335                         if (debug > 3)
 336                         {
 337                             Pout<< "Face: " << fIndex
 338                                 << " is handled by proc: "
 339                                 << procIndices_[pI]
 340                                 << ", so bailing out."
 341                                 << endl;
 342                         }
 343
 344                         return map;
 345                     }
 346
 347                     label curIndex = slaveMaps.size();
 348
 349                     // Size up the list
 350                     meshOps::sizeUpList
 351                     (
 352                         changeMap(),
 353                         slaveMaps
 354                     );
 355
 356                     // Save index and patch for posterity
 357                     slaveMaps[curIndex].index() = sIndex;
 358                     slaveMaps[curIndex].patchIndex() = pI;
 359
 360                     // Only one slave coupling is possible, so bail out
 361                     break;
 362                 }
 363             }
 364         }
 365         else
 366         {
 367             // Something's wrong with coupling maps
 368             FatalErrorIn
 369             (
 370                 "\n"
 371                 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
 372                 "(\n"
 373                 "    const label fIndex,\n"
 374                 "    const changeMap& masterMap,\n"
 375                 "    bool checkOnly,\n"
 376                 "    bool forceOp"
 377                 ")\n"
 378             )
 379                 << "Coupled maps were improperly specified." << nl
 380                 << " localCouple: " << localCouple << nl
 381                 << " procCouple: " << procCouple << nl
 382                 << " Face: " << fIndex << ": " << faces_[fIndex] << nl
 383                 << abort(FatalError);
 384         }
 385
 386         // Alias for convenience...
 387         changeMap& slaveMap = slaveMaps[0];
 388
 389         label sIndex = slaveMap.index();
 390         label pI = slaveMap.patchIndex();
 391         const coupleMap* cMapPtr = NULL;
 392
 393         // Temporarily turn off coupledModification
 394         unsetCoupledModification();
 395
 396         if (localCouple)
 397         {
 398             cMapPtr = &(patchCoupling_[pI].map());
 399
 400             // First check the slave for bisection feasibility.
 401             slaveMap = bisectQuadFace(sIndex, changeMap(), true, forceOp);
 402         }
 403         else
 404         if (procCouple)
 405         {
 406             cMapPtr = &(recvMeshes_[pI].map());
 407
 408             coupledMesh& recvMesh = recvMeshes_[pI];
 409
 410             // First check the slave for bisection feasibility.
 411             slaveMap =
 412             (
 413                 recvMesh.subMesh().bisectQuadFace
 414                 (
 415                     sIndex,
 416                     changeMap(),
 417                     true,
 418                     forceOp
 419                 )
 420             );
 421         }
 422
 423         // Turn it back on.
 424         setCoupledModification();
 425
 426         if (slaveMap.type() <= 0)
 427         {
 428             // Slave couldn't perform bisection.
 429             map.type() = -2;
 430
 431             return map;
 432         }
 433
 434         // Save index and patch for posterity
 435         slaveMap.index() = sIndex;
 436         slaveMap.patchIndex() = pI;
 437
 438         // Alias for convenience..
 439         const coupleMap& cMap = *cMapPtr;
 440
 441         // Temporarily turn off coupledModification
 442         unsetCoupledModification();
 443
 444         // First check the master for bisection feasibility.
 445         changeMap masterMap = bisectQuadFace(fIndex, changeMap(), true);
 446
 447         // Turn it back on.
 448         setCoupledModification();
 449
 450         // Master couldn't perform bisection
 451         if (masterMap.type() <= 0)
 452         {
 453             return masterMap;
 454         }
 455
 456         // Fill the masterMap with points that
 457         // we seek maps for...
 458         FixedList<labelList, 2> slaveEdges(labelList(2, -1));
 459
 460         forAll(slaveEdges, edgeI)
 461         {
 462             labelList& slaveEdge = slaveEdges[edgeI];
 463
 464             // Renumber to slave indices
 465             forAll(slaveEdge, pointI)
 466             {
 467                 slaveEdge[pointI] =
 468                 (
 469                     cMap.findSlave
 470                     (
 471                         pointEnum,
 472                         commonEdges[edgeI][pointI]
 473                     )
 474                 );
 475             }
 476
 477             masterMap.addPoint(-1, slaveEdge);
 478         }
 479
 480         // Temporarily turn off coupledModification
 481         unsetCoupledModification();
 482
 483         if (localCouple)
 484         {
 485             // Bisect the local slave.
 486             slaveMap = bisectQuadFace(sIndex, masterMap, false, forceOp);
 487         }
 488         else
 489         {
 490             coupledMesh& recvMesh = recvMeshes_[pI];
 491
 492             // Bisect the slave face
 493             slaveMap =
 494             (
 495                 recvMesh.subMesh().bisectQuadFace
 496                 (
 497                     sIndex,
 498                     masterMap,
 499                     false,
 500                     forceOp
 501                 )
 502             );
 503         }
 504
 505         // Turn coupledModification back on.
 506         setCoupledModification();
 507
 508         // The final operation has to succeed.
 509         if (slaveMap.type() <= 0)
 510         {
 511             FatalErrorIn
 512             (
 513                 "\n"
 514                 "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
 515                 "(\n"
 516                 "    const label fIndex,\n"
 517                 "    const changeMap& masterMap,\n"
 518                 "    bool checkOnly,\n"
 519                 "    bool forceOp"
 520                 ")\n"
 521             )
 522                 << "Coupled topo-change for slave failed."
 523                 << " Master face: " << fIndex << nl
 524                 << " Slave face: " << sIndex << nl
 525                 << " Patch index: " << pI << nl
 526                 << " Type: " << slaveMap.type() << nl
 527                 << abort(FatalError);
 528         }
 529
 530         // Save index and patch for posterity
 531         slaveMap.index() = sIndex;
 532         slaveMap.patchIndex() = pI;
 533     }
 534
 535     // Are we performing only checks?
 536     if (checkOnly)
 537     {
 538         map.type() = 1;
 539         return map;
 540     }
 541
 542     if (debug > 1)
 543     {
 544         Pout<< nl << nl
 545             << "Face: " << fIndex
 546             << ": " << faces_[fIndex]
 547             << " is to be bisected. " << endl;
 548
 549         if (debug > 2)
 550         {
 551             Pout<< " On SubMesh: " << isSubMesh_ << nl;
 552             Pout<< " coupledModification: " << coupledModification_ << nl;
 553
 554             const polyBoundaryMesh& boundary = boundaryMesh();
 555
 556             label epIndex = whichPatch(fIndex);
 557
 558             Pout<< " Patch: ";
 559
 560             if (epIndex == -1)
 561             {
 562                 Pout<< "Internal" << endl;
 563             }
 564             else
 565             if (epIndex < boundary.size())
 566             {
 567                 Pout<< boundary[epIndex].name() << endl;
 568             }
 569             else
 570             {
 571                 Pout<< " New patch: " << epIndex << endl;
 572             }
 573
 574             Pout<< "Cell[0]: " << c0 << ": " << oldCells[0] << endl;
 575
 576             forAll(oldCells[0], faceI)
 577             {
 578                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
 579
 580                 Pout<< oldCells[0][faceI] << ": "
 581                     << faces_[oldCells[0][faceI]]
 582                     << " fE: " << fE
 583                     << endl;
 584
 585                 forAll(fE, edgeI)
 586                 {
 587                     const labelList& eF = edgeFaces_[fE[edgeI]];
 588
 589                     Pout<< '\t' << fE[edgeI]
 590                         << ": " << edges_[fE[edgeI]]
 591                         << " eF: " << eF
 592                         << endl;
 593                 }
 594             }
 595         }
 596
 597         // Write out VTK files prior to change
 598         if (debug > 3)
 599         {
 600             labelList cellHull(2, -1);
 601
 602             cellHull[0] = owner_[fIndex];
 603             cellHull[1] = neighbour_[fIndex];
 604
 605             writeVTK
 606             (
 607                 Foam::name(fIndex)
 608               + "_Bisect_0",
 609                 cellHull
 610             );
 611         }
 612     }
 613
 614     // Find the isolated point on both boundary faces of cell[0]
 615     meshOps::findIsolatedPoint
 616     (
 617         c0BdyFace[0],
 618         commonEdges[0],
 619         otherPointIndex[0],
 620         nextToOtherPoint[0]
 621     );
 622
 623     meshOps::findIsolatedPoint
 624     (
 625         c0BdyFace[1],
 626         commonEdges[1],
 627         otherPointIndex[1],
 628         nextToOtherPoint[1]
 629     );
 630
 631     // For convenience...
 632     otherEdgePoint[0] = commonEdges[0].otherVertex(nextToOtherPoint[0]);
 633     otherEdgePoint[1] = commonEdges[1].otherVertex(nextToOtherPoint[1]);
 634
 635     labelList mP(2, -1);
 636
 637     // Set mapping for this point
 638     mP[0] = commonEdges[0][0];
 639     mP[1] = commonEdges[0][1];
 640
 641     // Add two new points to the end of the list
 642     newPointIndex[0] =
 643     (
 644         insertPoint
 645         (
 646             0.5 * (points_[mP[0]] + points_[mP[1]]),
 647             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
 648             mP
 649         )
 650     );
 651
 652     // Set mapping for this point
 653     mP[0] = commonEdges[1][0];
 654     mP[1] = commonEdges[1][1];
 655
 656     newPointIndex[1] =
 657     (
 658         insertPoint
 659         (
 660             0.5 * (points_[mP[0]] + points_[mP[1]]),
 661             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
 662             mP
 663         )
 664     );
 665
 666     // Add the points to the map. Since this might require master mapping,
 667     // first check to see if a slave is being bisected.
 668     if (masterMap.addedPointList().size())
 669     {
 670         const List<objectMap>& pMap =
 671         (
 672             masterMap.addedPointList()
 673         );
 674
 675         edge firstEdge
 676         (
 677             pMap[0].masterObjects()[0],
 678             pMap[0].masterObjects()[1]
 679         );
 680
 681         edge secondEdge
 682         (
 683             pMap[1].masterObjects()[0],
 684             pMap[1].masterObjects()[1]
 685         );
 686
 687         if (firstEdge == commonEdges[0] && secondEdge == commonEdges[1])
 688         {
 689             // Update in conventional order
 690             map.addPoint(newPointIndex[0]);
 691             map.addPoint(newPointIndex[1]);
 692         }
 693         else
 694         if (firstEdge == commonEdges[1] && secondEdge == commonEdges[0])
 695         {
 696             // Update in reverse order
 697             map.addPoint(newPointIndex[1]);
 698             map.addPoint(newPointIndex[0]);
 699         }
 700         else
 701         {
 702             // We have a problem
 703             FatalErrorIn
 704             (
 705                 "\n"
 706                 "const changeMap "
 707                 "dynamicTopoFvMesh::bisectQuadFace\n"
 708                 "(\n"
 709                 "    const label fIndex,\n"
 710                 "    const changeMap& masterMap,\n"
 711                 "    bool checkOnly,\n"
 712                 "    bool forceOp\n"
 713                 ")\n"
 714             )
 715                 << "Coupled topo-change for slave failed."
 716                 << " firstEdge: " << firstEdge << nl
 717                 << " secondEdge: " << secondEdge << nl
 718                 << " commonEdges[0]: " << commonEdges[0] << nl
 719                 << " commonEdges[1]: " << commonEdges[1] << nl
 720                 << abort(FatalError);
 721         }
 722     }
 723     else
 724     {
 725         map.addPoint(newPointIndex[0]);
 726         map.addPoint(newPointIndex[1]);
 727     }
 728
 729     // Add a new prism cell to the end of the list.
 730     // Currently invalid, but will be updated later.
 731     newCellIndex[0] = insertCell(newCells[0], lengthScale_[c0]);
 732
 733     // Add this cell to the map.
 734     map.addCell(newCellIndex[0]);
 735
 736     // Modify the two existing triangle boundary faces
 737
 738     // Zeroth boundary face - Owner = c[0] & Neighbour [-1] (unchanged)
 739     meshOps::replaceLabel
 740     (
 741         otherEdgePoint[0],
 742         newPointIndex[0],
 743         c0BdyFace[0]
 744     );
 745
 746     // First boundary face - Owner = newCell[0], Neighbour = -1
 747     meshOps::replaceLabel
 748     (
 749         otherEdgePoint[1],
 750         newPointIndex[1],
 751         c0BdyFace[1]
 752     );
 753
 754     // Update faces.
 755     faces_[c0BdyIndex[0]] = c0BdyFace[0];
 756     faces_[c0BdyIndex[1]] = c0BdyFace[1];
 757
 758     owner_[c0BdyIndex[1]] = newCellIndex[0];
 759
 760     meshOps::replaceLabel
 761     (
 762         c0BdyIndex[1],
 763         -1,
 764         oldCells[0]
 765     );
 766
 767     // Detect edges other than commonEdges
 768     const labelList& fEdges = faceEdges_[fIndex];
 769
 770     forAll(fEdges, edgeI)
 771     {
 772         if
 773         (
 774             fEdges[edgeI] != commonEdgeIndex[0] &&
 775             fEdges[edgeI] != commonEdgeIndex[1]
 776         )
 777         {
 778             // Obtain a reference to this edge
 779             const edge& eThis = edges_[fEdges[edgeI]];
 780
 781             if
 782             (
 783                 eThis[0] == nextToOtherPoint[0]
 784              || eThis[1] == nextToOtherPoint[0]
 785             )
 786             {
 787                 otherEdgeIndex[0] = fEdges[edgeI];
 788             }
 789             else
 790             {
 791                 otherEdgeIndex[1] = fEdges[edgeI];
 792             }
 793         }
 794     }
 795
 796     // Modify point-labels on the quad face under consideration
 797     meshOps::replaceLabel
 798     (
 799         otherEdgePoint[0],
 800         newPointIndex[0],
 801         faces_[fIndex]
 802     );
 803
 804     meshOps::replaceLabel
 805     (
 806         nextToOtherPoint[1],
 807         newPointIndex[1],
 808         faces_[fIndex]
 809     );
 810
 811     // Add this face to the map.
 812     // Although this face isn't technically 'added', it's
 813     // required for coupled patch mapping.
 814     map.addFace(fIndex);
 815
 816     if (debug > 1)
 817     {
 818         Pout<< "Modified face: " << fIndex
 819             << ": " << faces_[fIndex] << endl;
 820
 821         if (debug > 2)
 822         {
 823             Pout<< "Common edges: " << nl
 824                 << commonEdgeIndex[0] << ": " << commonEdges[0] << nl
 825                 << commonEdgeIndex[1] << ": " << commonEdges[1]
 826                 << endl;
 827         }
 828     }
 829
 830     // Find the quad face that contains otherEdgeIndex[1]
 831     found = false;
 832
 833     const labelList& e1 = faceEdges_[c0IntIndex[0]];
 834
 835     forAll(e1, edgeI)
 836     {
 837         if (e1[edgeI] == otherEdgeIndex[1])
 838         {
 839             meshOps::replaceLabel
 840             (
 841                 c0IntIndex[0],
 842                 -1,
 843                 oldCells[0]
 844             );
 845
 846             replaceFace = c0IntIndex[0];
 847             retainFace = c0IntIndex[1];
 848             found = true;
 849             break;
 850         }
 851     }
 852
 853     if (!found)
 854     {
 855         // The edge was not found before
 856         meshOps::replaceLabel
 857         (
 858             c0IntIndex[1],
 859             -1,
 860             oldCells[0]
 861         );
 862
 863         replaceFace = c0IntIndex[1];
 864         retainFace = c0IntIndex[0];
 865     }
 866
 867     // Check if face reversal is necessary for the replacement
 868     if (owner_[replaceFace] == c0)
 869     {
 870         if (neighbour_[replaceFace] == -1)
 871         {
 872             // Change the owner
 873             owner_[replaceFace] = newCellIndex[0];
 874         }
 875         else
 876         {
 877             // This face has to be reversed
 878             faces_[replaceFace] = faces_[replaceFace].reverseFace();
 879             owner_[replaceFace] = neighbour_[replaceFace];
 880             neighbour_[replaceFace] = newCellIndex[0];
 881
 882             setFlip(replaceFace);
 883         }
 884     }
 885     else
 886     {
 887         // Keep owner, but change neighbour
 888         neighbour_[replaceFace] = newCellIndex[0];
 889     }
 890
 891     // Define the faces for the new cell
 892     newCells[0][0] = c0BdyIndex[1];
 893     newCells[0][1] = replaceFace;
 894
 895     // Define the set of new faces and insert...
 896
 897     // New interior face; Owner = cell[0] & Neighbour = newCell[0]
 898     tmpQuadFace[0] = otherPointIndex[0];
 899     tmpQuadFace[1] = newPointIndex[0];
 900     tmpQuadFace[2] = newPointIndex[1];
 901     tmpQuadFace[3] = otherPointIndex[1];
 902
 903     newFaceIndex[0] =
 904     (
 905         insertFace
 906         (
 907             -1,
 908             tmpQuadFace,
 909             c0,
 910             newCellIndex[0],
 911             tmpQFEdges
 912         )
 913     );
 914
 915     // Add this face to the map.
 916     map.addFace(newFaceIndex[0]);
 917
 918     // Find the common edge between quad/quad faces...
 919     meshOps::findCommonEdge
 920     (
 921         c0IntIndex[0],
 922         c0IntIndex[1],
 923         faceEdges_,
 924         otherEdgeIndex[2]
 925     );
 926
 927     // ... and size-up edgeFaces for the edge
 928     meshOps::sizeUpList
 929     (
 930         newFaceIndex[0],
 931         edgeFaces_[otherEdgeIndex[2]]
 932     );
 933
 934     meshOps::replaceLabel
 935     (
 936         -1,
 937         newFaceIndex[0],
 938         newCells[0]
 939     );
 940
 941     meshOps::replaceLabel
 942     (
 943         -1,
 944         newFaceIndex[0],
 945         oldCells[0]
 946     );
 947
 948     // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
 949     tmpTriFace[0] = otherPointIndex[0];
 950     tmpTriFace[1] = newPointIndex[0];
 951     tmpTriFace[2] = otherEdgePoint[0];
 952
 953     newFaceIndex[1] =
 954     (
 955         insertFace
 956         (
 957             whichPatch(c0BdyIndex[0]),
 958             tmpTriFace,
 959             newCellIndex[0],
 960             -1,
 961             tmpTFEdges
 962         )
 963     );
 964
 965     // Add this face to the map.
 966     map.addFace(newFaceIndex[1]);
 967
 968     meshOps::replaceLabel
 969     (
 970         -1,
 971         newFaceIndex[1],
 972         newCells[0]
 973     );
 974
 975     // Third boundary face; Owner = c[0] & Neighbour = [-1]
 976     tmpTriFace[0] = otherPointIndex[1];
 977     tmpTriFace[1] = newPointIndex[1];
 978     tmpTriFace[2] = otherEdgePoint[1];
 979
 980     newFaceIndex[2] =
 981     (
 982         insertFace
 983         (
 984             whichPatch(c0BdyIndex[1]),
 985             tmpTriFace,
 986             c0,
 987             -1,
 988             tmpTFEdges
 989         )
 990     );
 991
 992     // Add this face to the map.
 993     map.addFace(newFaceIndex[2]);
 994
 995     meshOps::replaceLabel
 996     (
 997         -1,
 998         newFaceIndex[2],
 999         oldCells[0]
1000     );
1001
1002     // Create / modify edges...
1003     labelList tmpTriEdgeFaces(3, -1);
1004
1005     // The edge bisecting the zeroth boundary triangular face
1006     tmpTriEdgeFaces[0] = c0BdyIndex[0];
1007     tmpTriEdgeFaces[1] = newFaceIndex[0];
1008     tmpTriEdgeFaces[2] = newFaceIndex[1];
1009
1010     newEdgeIndex[1] =
1011     (
1012         insertEdge
1013         (
1014             whichEdgePatch(commonEdgeIndex[0]),
1015             edge(newPointIndex[0], otherPointIndex[0]),
1016             tmpTriEdgeFaces
1017         )
1018     );
1019
1020     // Add this edge to the map.
1021     map.addEdge(newEdgeIndex[1]);
1022
1023     // Find the common edge between the quad/tri faces...
1024     meshOps::findCommonEdge
1025     (
1026         c0BdyIndex[0],
1027         replaceFace,
1028         faceEdges_,
1029         cornerEdgeIndex[0]
1030     );
1031
1032     // ...and correct faceEdges / edgeFaces
1033     meshOps::replaceLabel
1034     (
1035         cornerEdgeIndex[0],
1036         newEdgeIndex[1],
1037         faceEdges_[c0BdyIndex[0]]
1038     );
1039
1040     meshOps::replaceLabel
1041     (
1042         c0BdyIndex[0],
1043         newFaceIndex[1],
1044         edgeFaces_[cornerEdgeIndex[0]]
1045     );
1046
1047     // The edge bisecting the first boundary triangular face
1048     tmpTriEdgeFaces[0] = c0BdyIndex[1];
1049     tmpTriEdgeFaces[1] = newFaceIndex[0];
1050     tmpTriEdgeFaces[2] = newFaceIndex[2];
1051
1052     newEdgeIndex[2] =
1053     (
1054         insertEdge
1055         (
1056             whichEdgePatch(commonEdgeIndex[1]),
1057             edge(newPointIndex[1], otherPointIndex[1]),
1058             tmpTriEdgeFaces
1059         )
1060     );
1061
1062     // Add this edge to the map.
1063     map.addEdge(newEdgeIndex[2]);
1064
1065     // Find the common edge between the quad/tri faces...
1066     meshOps::findCommonEdge
1067     (
1068         c0BdyIndex[1],
1069         retainFace,
1070         faceEdges_,
1071         cornerEdgeIndex[1]
1072     );
1073
1074     // ...and correct faceEdges / edgeFaces
1075     meshOps::replaceLabel
1076     (
1077         cornerEdgeIndex[1],
1078         newEdgeIndex[2],
1079         faceEdges_[c0BdyIndex[1]]
1080     );
1081
1082     meshOps::replaceLabel
1083     (
1084         c0BdyIndex[1],
1085         newFaceIndex[2],
1086         edgeFaces_[cornerEdgeIndex[1]]
1087     );
1088
1089     if (c1 == -1)
1090     {
1091         // The quad boundary face resulting from bisection;
1092         // Owner = newCell[0] & Neighbour = [-1]
1093         tmpQuadFace[0] = newPointIndex[1];
1094         tmpQuadFace[1] = nextToOtherPoint[1];
1095         tmpQuadFace[2] = otherEdgePoint[0];
1096         tmpQuadFace[3] = newPointIndex[0];
1097
1098         newFaceIndex[3] =
1099         (
1100             insertFace
1101             (
1102                 whichPatch(fIndex),
1103                 tmpQuadFace,
1104                 newCellIndex[0],
1105                 -1,
1106                 tmpQFEdges
1107             )
1108         );
1109
1110         // Add this face to the map.
1111         map.addFace(newFaceIndex[3]);
1112
1113         // Correct edgeFaces for otherEdgeIndex[1]
1114         meshOps::replaceLabel
1115         (
1116             fIndex,
1117             newFaceIndex[3],
1118             edgeFaces_[otherEdgeIndex[1]]
1119         );
1120
1121         meshOps::replaceLabel
1122         (
1123             -1,
1124             newFaceIndex[3],
1125             newCells[0]
1126         );
1127
1128         labelList tmpBiEdgeFaces(2, -1);
1129
1130         // The edge bisecting the face
1131         tmpTriEdgeFaces[0] = newFaceIndex[3];
1132         tmpTriEdgeFaces[1] = newFaceIndex[0];
1133         tmpTriEdgeFaces[2] = fIndex;
1134
1135         newEdgeIndex[0] =
1136         (
1137             insertEdge
1138             (
1139                 whichEdgePatch(otherEdgeIndex[0]),
1140                 edge(newPointIndex[0], newPointIndex[1]),
1141                 tmpTriEdgeFaces
1142             )
1143         );
1144
1145         // Add this edge to the map.
1146         map.addEdge(newEdgeIndex[0]);
1147
1148         // Replace an edge on the bisected face
1149         meshOps::replaceLabel
1150         (
1151             otherEdgeIndex[1],
1152             newEdgeIndex[0],
1153             faceEdges_[fIndex]
1154         );
1155
1156         // Create / replace side edges created from face bisection
1157         tmpBiEdgeFaces[0] = newFaceIndex[1];
1158         tmpBiEdgeFaces[1] = newFaceIndex[3];
1159
1160         newEdgeIndex[3] =
1161         (
1162             insertEdge
1163             (
1164                 whichEdgePatch(commonEdgeIndex[0]),
1165                 edge(newPointIndex[0], otherEdgePoint[0]),
1166                 tmpBiEdgeFaces
1167             )
1168         );
1169
1170         // Add this edge to the map.
1171         map.addEdge(newEdgeIndex[3]);
1172
1173         tmpBiEdgeFaces[0] = c0BdyIndex[1];
1174         tmpBiEdgeFaces[1] = newFaceIndex[3];
1175
1176         newEdgeIndex[4] =
1177         (
1178             insertEdge
1179             (
1180                 whichEdgePatch(commonEdgeIndex[1]),
1181                 edge(newPointIndex[1], nextToOtherPoint[1]),
1182                 tmpBiEdgeFaces
1183             )
1184         );
1185
1186         // Add this edge to the map.
1187         map.addEdge(newEdgeIndex[4]);
1188
1189         // Now that edges are defined, configure faceEdges
1190         // for all new faces
1191
1192         // The quad interior face; Owner = cell[0] & Neighbour = newCell[0]
1193         tmpQFEdges[0] = newEdgeIndex[0];
1194         tmpQFEdges[1] = newEdgeIndex[1];
1195         tmpQFEdges[2] = otherEdgeIndex[2];
1196         tmpQFEdges[3] = newEdgeIndex[2];
1197         faceEdges_[newFaceIndex[0]] = tmpQFEdges;
1198
1199         // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
1200         tmpTFEdges[0] = newEdgeIndex[3];
1201         tmpTFEdges[1] = cornerEdgeIndex[0];
1202         tmpTFEdges[2] = newEdgeIndex[1];
1203         faceEdges_[newFaceIndex[1]] = tmpTFEdges;
1204
1205         // Third boundary face; Owner = c[0] & Neighbour = [-1]
1206         tmpTFEdges[0] = newEdgeIndex[2];
1207         tmpTFEdges[1] = cornerEdgeIndex[1];
1208         tmpTFEdges[2] = commonEdgeIndex[1];
1209         faceEdges_[newFaceIndex[2]] = tmpTFEdges;
1210
1211         // The quad face from bisection:
1212         tmpQFEdges[0] = otherEdgeIndex[1];
1213         tmpQFEdges[1] = newEdgeIndex[3];
1214         tmpQFEdges[2] = newEdgeIndex[0];
1215         tmpQFEdges[3] = newEdgeIndex[4];
1216         faceEdges_[newFaceIndex[3]] = tmpQFEdges;
1217
1218         meshOps::replaceLabel
1219         (
1220             commonEdgeIndex[1],
1221             newEdgeIndex[4],
1222             faceEdges_[c0BdyIndex[1]]
1223         );
1224
1225         meshOps::replaceLabel
1226         (
1227             c0BdyIndex[1],
1228             newFaceIndex[2],
1229             edgeFaces_[commonEdgeIndex[1]]
1230         );
1231
1232         if (debug > 2)
1233         {
1234             Pout<< "Modified Cell[0]: "
1235                 << c0 << ": " << oldCells[0] << endl;
1236
1237             forAll(oldCells[0], faceI)
1238             {
1239                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
1240
1241                 Pout<< oldCells[0][faceI]
1242                     << ": " << faces_[oldCells[0][faceI]]
1243                     << " fE: " << fE
1244                     << endl;
1245
1246                 forAll(fE, edgeI)
1247                 {
1248                     const labelList& eF = edgeFaces_[fE[edgeI]];
1249
1250                     Pout<< '\t' << fE[edgeI]
1251                         << ": " << edges_[fE[edgeI]]
1252                         << " eF: " << eF
1253                         << endl;
1254                 }
1255             }
1256
1257             Pout<< "New Cell[0]: " << newCellIndex[0]
1258                 << ": " << newCells[0] << endl;
1259
1260             forAll(newCells[0], faceI)
1261             {
1262                 const labelList& fE = faceEdges_[newCells[0][faceI]];
1263
1264                 Pout<< newCells[0][faceI]
1265                     << ": " << faces_[newCells[0][faceI]]
1266                     << " fE: " << fE
1267                     << endl;
1268
1269                 forAll(fE, edgeI)
1270                 {
1271                     const labelList& eF = edgeFaces_[fE[edgeI]];
1272
1273                     Pout<< '\t' << fE[edgeI]
1274                         << ": " << edges_[fE[edgeI]]
1275                         << " eF: " << eF
1276                         << endl;
1277                 }
1278             }
1279         }
1280     }
1281     else
1282     {
1283         oldCells[1] = cells_[c1];
1284
1285         // Add a new prism cell to the end of the list.
1286         // Currently invalid, but will be updated later.
1287         newCellIndex[1] = insertCell(newCells[1], lengthScale_[c1]);
1288
1289         // Add this cell to the map.
1290         map.addCell(newCellIndex[1]);
1291
1292         if (debug > 2)
1293         {
1294             Pout<< "Cell[1]: " << c1 << ": " << oldCells[1] << endl;
1295
1296             forAll(oldCells[1], faceI)
1297             {
1298                 const labelList& fE = faceEdges_[oldCells[1][faceI]];
1299
1300                 Pout<< oldCells[1][faceI] << ": "
1301                     << faces_[oldCells[1][faceI]]
1302                     << " fE: " << fE
1303                     << endl;
1304
1305                 forAll(fE, edgeI)
1306                 {
1307                     const labelList& eF = edgeFaces_[fE[edgeI]];
1308
1309                     Pout<< '\t' << fE[edgeI]
1310                         << ": " << edges_[fE[edgeI]]
1311                         << " eF: " << eF
1312                         << endl;
1313                 }
1314             }
1315         }
1316
1317         // Find the interior face that contains otherEdgeIndex[1]
1318         found = false;
1319
1320         const labelList& e2 = faceEdges_[c1IntIndex[0]];
1321
1322         forAll(e2, edgeI)
1323         {
1324             if (e2[edgeI] == otherEdgeIndex[1])
1325             {
1326                 meshOps::replaceLabel
1327                 (
1328                     c1IntIndex[0],
1329                     -1,
1330                     oldCells[1]
1331                 );
1332
1333                 replaceFace = c1IntIndex[0];
1334                 retainFace = c1IntIndex[1];
1335                 found = true;
1336                 break;
1337             }
1338         }
1339
1340         if (!found)
1341         {
1342             // The edge was not found before
1343             meshOps::replaceLabel
1344             (
1345                 c1IntIndex[1],
1346                 -1,
1347                 oldCells[1]
1348             );
1349
1350             replaceFace = c1IntIndex[1];
1351             retainFace = c1IntIndex[0];
1352         }
1353
1354         // Check if face reversal is necessary for the replacement
1355         if (owner_[replaceFace] == c1)
1356         {
1357             if (neighbour_[replaceFace] == -1)
1358             {
1359                 // Change the owner
1360                 owner_[replaceFace] = newCellIndex[1];
1361             }
1362             else
1363             {
1364                 // This face has to be reversed
1365                 faces_[replaceFace] = faces_[replaceFace].reverseFace();
1366                 owner_[replaceFace] = neighbour_[replaceFace];
1367                 neighbour_[replaceFace] = newCellIndex[1];
1368
1369                 setFlip(replaceFace);
1370             }
1371         }
1372         else
1373         {
1374             // Keep owner, but change neighbour
1375             neighbour_[replaceFace] = newCellIndex[1];
1376         }
1377
1378         // Define attributes for the new prism cell
1379         newCells[1][0] = replaceFace;
1380
1381         // The quad interior face resulting from bisection;
1382         // Owner = newCell[0] & Neighbour = newCell[1]
1383         tmpQuadFace[0] = newPointIndex[1];
1384         tmpQuadFace[1] = nextToOtherPoint[1];
1385         tmpQuadFace[2] = otherEdgePoint[0];
1386         tmpQuadFace[3] = newPointIndex[0];
1387
1388         newFaceIndex[3] =
1389         (
1390             insertFace
1391             (
1392                 -1,
1393                 tmpQuadFace,
1394                 newCellIndex[0],
1395                 newCellIndex[1],
1396                 tmpQFEdges
1397             )
1398         );
1399
1400         // Add this face to the map.
1401         map.addFace(newFaceIndex[3]);
1402
1403         // Correct edgeFaces for otherEdgeIndex[1]
1404         meshOps::replaceLabel
1405         (
1406             fIndex,
1407             newFaceIndex[3],
1408             edgeFaces_[otherEdgeIndex[1]]
1409         );
1410
1411         meshOps::replaceLabel
1412         (
1413             -1,
1414             newFaceIndex[3],
1415             newCells[0]
1416         );
1417
1418         meshOps::replaceLabel
1419         (
1420             -1,
1421             newFaceIndex[3],
1422             newCells[1]
1423         );
1424
1425         newCells[1][1] = newFaceIndex[3];
1426
1427         // Check for common edges among the two boundary faces
1428         // Find the isolated point on both boundary faces of cell[1]
1429         if
1430         (
1431             meshOps::findCommonEdge
1432             (
1433                 c1BdyIndex[0],
1434                 c0BdyIndex[0],
1435                 faceEdges_,
1436                 commonEdgeIndex[2]
1437             )
1438         )
1439         {
1440             meshOps::findCommonEdge
1441             (
1442                 c1BdyIndex[1],
1443                 c0BdyIndex[1],
1444                 faceEdges_,
1445                 commonEdgeIndex[3]
1446             );
1447
1448             commonFaceIndex[2] = c1BdyIndex[0];
1449             commonFaceIndex[3] = c1BdyIndex[1];
1450         }
1451         else
1452         {
1453             meshOps::findCommonEdge
1454             (
1455                 c1BdyIndex[0],
1456                 c0BdyIndex[1],
1457                 faceEdges_,
1458                 commonEdgeIndex[3]
1459             );
1460
1461             meshOps::findCommonEdge
1462             (
1463                 c1BdyIndex[1],
1464                 c0BdyIndex[0],
1465                 faceEdges_,
1466                 commonEdgeIndex[2]
1467             );
1468
1469             commonFaceIndex[2] = c1BdyIndex[1];
1470             commonFaceIndex[3] = c1BdyIndex[0];
1471         }
1472
1473         commonEdges[2] = edges_[commonEdgeIndex[2]];
1474         commonEdges[3] = edges_[commonEdgeIndex[3]];
1475
1476         if (debug > 2)
1477         {
1478             Pout<< "Common edges: " << nl
1479                 << commonEdgeIndex[2] << ": " << commonEdges[2] << nl
1480                 << commonEdgeIndex[3] << ": " << commonEdges[3]
1481                 << endl;
1482         }
1483
1484         meshOps::findIsolatedPoint
1485         (
1486             faces_[commonFaceIndex[2]],
1487             commonEdges[2],
1488             otherPointIndex[2],
1489             nextToOtherPoint[2]
1490         );
1491
1492         meshOps::findIsolatedPoint
1493         (
1494             faces_[commonFaceIndex[3]],
1495             commonEdges[3],
1496             otherPointIndex[3],
1497             nextToOtherPoint[3]
1498         );
1499
1500         // For convenience...
1501         otherEdgePoint[2] = commonEdges[2].otherVertex(nextToOtherPoint[2]);
1502         otherEdgePoint[3] = commonEdges[3].otherVertex(nextToOtherPoint[3]);
1503
1504         // Modify the two existing triangle boundary faces
1505
1506         // Zeroth boundary face - Owner = newCell[1], Neighbour = -1
1507         meshOps::replaceLabel
1508         (
1509             otherEdgePoint[2],
1510             newPointIndex[0],
1511             faces_[commonFaceIndex[2]]
1512         );
1513
1514         owner_[commonFaceIndex[2]] = newCellIndex[1];
1515
1516         meshOps::replaceLabel
1517         (
1518             commonFaceIndex[2],
1519             -1,
1520             oldCells[1]
1521         );
1522
1523         newCells[1][2] = commonFaceIndex[2];
1524
1525         // First boundary face - Owner = c[1] & Neighbour [-1] (unchanged)
1526         meshOps::replaceLabel
1527         (
1528             otherEdgePoint[3],
1529             newPointIndex[1],
1530             faces_[commonFaceIndex[3]]
1531         );
1532
1533         // New interior face; Owner = cell[1] & Neighbour = newCell[1]
1534         tmpQuadFace[0] = newPointIndex[0];
1535         tmpQuadFace[1] = otherPointIndex[2];
1536         tmpQuadFace[2] = otherPointIndex[3];
1537         tmpQuadFace[3] = newPointIndex[1];
1538
1539         newFaceIndex[4] =
1540         (
1541             insertFace
1542             (
1543                 -1,
1544                 tmpQuadFace,
1545                 c1,
1546                 newCellIndex[1],
1547                 tmpQFEdges
1548             )
1549         );
1550
1551         // Add this face to the map.
1552         map.addFace(newFaceIndex[4]);
1553
1554         // Find the common edge between quad/quad faces...
1555         meshOps::findCommonEdge
1556         (
1557             c1IntIndex[0],
1558             c1IntIndex[1],
1559             faceEdges_,
1560             otherEdgeIndex[3]
1561         );
1562
1563         // ... and size-up edgeFaces for the edge
1564         meshOps::sizeUpList
1565         (
1566             newFaceIndex[4],
1567             edgeFaces_[otherEdgeIndex[3]]
1568         );
1569
1570         meshOps::replaceLabel
1571         (
1572             -1,
1573             newFaceIndex[4],
1574             newCells[1]
1575         );
1576
1577         meshOps::replaceLabel
1578         (
1579             -1,
1580             newFaceIndex[4],
1581             oldCells[1]
1582         );
1583
1584         // Second boundary face; Owner = cell[1] & Neighbour [-1]
1585         tmpTriFace[0] = otherPointIndex[2];
1586         tmpTriFace[1] = newPointIndex[0];
1587         tmpTriFace[2] = otherEdgePoint[2];
1588
1589         newFaceIndex[5] =
1590         (
1591             insertFace
1592             (
1593                 whichPatch(commonFaceIndex[2]),
1594                 tmpTriFace,
1595                 c1,
1596                 -1,
1597                 tmpTFEdges
1598             )
1599         );
1600
1601         // Add this face to the map.
1602         map.addFace(newFaceIndex[5]);
1603
1604         meshOps::replaceLabel
1605         (
1606             -1,
1607             newFaceIndex[5],
1608             oldCells[1]
1609         );
1610
1611         // Third boundary face; Owner = newCell[1] & Neighbour [-1]
1612         tmpTriFace[0] = otherPointIndex[3];
1613         tmpTriFace[1] = newPointIndex[1];
1614         tmpTriFace[2] = otherEdgePoint[3];
1615
1616         newFaceIndex[6] =
1617         (
1618             insertFace
1619             (
1620                 whichPatch(commonFaceIndex[3]),
1621                 tmpTriFace,
1622                 newCellIndex[1],
1623                 -1,
1624                 tmpTFEdges
1625             )
1626         );
1627
1628         // Add this face to the map.
1629         map.addFace(newFaceIndex[6]);
1630
1631         meshOps::replaceLabel
1632         (
1633             -1,
1634             newFaceIndex[6],
1635             newCells[1]
1636         );
1637
1638         // Create / modify edges...
1639         labelList tmpQuadEdgeFaces(4, -1);
1640
1641         // The internal edge bisecting the face
1642         tmpQuadEdgeFaces[0] = fIndex;
1643         tmpQuadEdgeFaces[1] = newFaceIndex[0];
1644         tmpQuadEdgeFaces[2] = newFaceIndex[3];
1645         tmpQuadEdgeFaces[3] = newFaceIndex[4];
1646
1647         newEdgeIndex[0] =
1648         (
1649             insertEdge
1650             (
1651                 -1,
1652                 edge(newPointIndex[0], newPointIndex[1]),
1653                 tmpQuadEdgeFaces
1654             )
1655         );
1656
1657         // Add this edge to the map.
1658         map.addEdge(newEdgeIndex[0]);
1659
1660         // Replace an edge on the bisected face
1661         meshOps::replaceLabel
1662         (
1663             otherEdgeIndex[1],
1664             newEdgeIndex[0],
1665             faceEdges_[fIndex]
1666         );
1667
1668         // Create / replace side edges created from face bisection
1669         tmpTriEdgeFaces[0] = commonFaceIndex[2];
1670         tmpTriEdgeFaces[1] = newFaceIndex[3];
1671         tmpTriEdgeFaces[2] = newFaceIndex[1];
1672
1673         newEdgeIndex[3] =
1674         (
1675             insertEdge
1676             (
1677                 whichEdgePatch(commonEdgeIndex[2]),
1678                 edge(newPointIndex[0], nextToOtherPoint[2]),
1679                 tmpTriEdgeFaces
1680             )
1681         );
1682
1683         // Add this edge to the map.
1684         map.addEdge(newEdgeIndex[3]);
1685
1686         tmpTriEdgeFaces[0] = c0BdyIndex[1];
1687         tmpTriEdgeFaces[1] = newFaceIndex[3];
1688         tmpTriEdgeFaces[2] = newFaceIndex[6];
1689
1690         newEdgeIndex[4] =
1691         (
1692             insertEdge
1693             (
1694                 whichEdgePatch(commonEdgeIndex[3]),
1695                 edge(newPointIndex[1], otherEdgePoint[3]),
1696                 tmpTriEdgeFaces
1697             )
1698         );
1699
1700         // Add this edge to the map.
1701         map.addEdge(newEdgeIndex[4]);
1702
1703         // The edge bisecting the second boundary triangular face
1704         tmpTriEdgeFaces[0] = commonFaceIndex[2];
1705         tmpTriEdgeFaces[1] = newFaceIndex[4];
1706         tmpTriEdgeFaces[2] = newFaceIndex[5];
1707
1708         newEdgeIndex[5] =
1709         (
1710             insertEdge
1711             (
1712                 whichEdgePatch(commonEdgeIndex[2]),
1713                 edge(newPointIndex[0], otherPointIndex[2]),
1714                 tmpTriEdgeFaces
1715             )
1716         );
1717
1718         // Add this edge to the map.
1719         map.addEdge(newEdgeIndex[5]);
1720
1721         // Find the common edge between the quad/tri faces...
1722         meshOps::findCommonEdge
1723         (
1724             commonFaceIndex[2],
1725             retainFace,
1726             faceEdges_,
1727             cornerEdgeIndex[2]
1728         );
1729
1730         // ...and correct faceEdges / edgeFaces
1731         meshOps::replaceLabel
1732         (
1733             cornerEdgeIndex[2],
1734             newEdgeIndex[5],
1735             faceEdges_[commonFaceIndex[2]]
1736         );
1737
1738         meshOps::replaceLabel
1739         (
1740             commonFaceIndex[2],
1741             newFaceIndex[5],
1742             edgeFaces_[cornerEdgeIndex[2]]
1743         );
1744
1745         // The edge bisecting the third boundary triangular face
1746         tmpTriEdgeFaces[0] = commonFaceIndex[3];
1747         tmpTriEdgeFaces[1] = newFaceIndex[4];
1748         tmpTriEdgeFaces[2] = newFaceIndex[6];
1749
1750         newEdgeIndex[6] =
1751         (
1752             insertEdge
1753             (
1754                 whichEdgePatch(commonEdgeIndex[3]),
1755                 edge(newPointIndex[1], otherPointIndex[3]),
1756                 tmpTriEdgeFaces
1757             )
1758         );
1759
1760         // Add this edge to the map.
1761         map.addEdge(newEdgeIndex[6]);
1762
1763         // Find the common edge between the quad/tri faces...
1764         meshOps::findCommonEdge
1765         (
1766             commonFaceIndex[3],
1767             replaceFace,
1768             faceEdges_,
1769             cornerEdgeIndex[3]
1770         );
1771
1772         // ...and correct faceEdges / edgeFaces
1773         meshOps::replaceLabel
1774         (
1775             cornerEdgeIndex[3],
1776             newEdgeIndex[6],
1777             faceEdges_[commonFaceIndex[3]]
1778         );
1779
1780         meshOps::replaceLabel
1781         (
1782             commonFaceIndex[3],
1783             newFaceIndex[6],
1784             edgeFaces_[cornerEdgeIndex[3]]
1785         );
1786
1787         // Now that edges are defined, configure faceEdges
1788         // for all new faces
1789
1790         // The quad interior face; Owner = c[0] & Neighbour = newCell[0]
1791         tmpQFEdges[0] = newEdgeIndex[0];
1792         tmpQFEdges[1] = newEdgeIndex[1];
1793         tmpQFEdges[2] = otherEdgeIndex[2];
1794         tmpQFEdges[3] = newEdgeIndex[2];
1795         faceEdges_[newFaceIndex[0]] = tmpQFEdges;
1796
1797         // Second boundary face; Owner = newCell[0] & Neighbour = [-1]
1798         tmpTFEdges[0] = newEdgeIndex[3];
1799         tmpTFEdges[1] = cornerEdgeIndex[0];
1800         tmpTFEdges[2] = newEdgeIndex[1];
1801         faceEdges_[newFaceIndex[1]] = tmpTFEdges;
1802
1803         // Third boundary face; Owner = c[0] & Neighbour = [-1]
1804         tmpTFEdges[0] = newEdgeIndex[2];
1805         tmpTFEdges[1] = cornerEdgeIndex[1];
1806         tmpTFEdges[2] = commonEdgeIndex[3];
1807         faceEdges_[newFaceIndex[2]] = tmpTFEdges;
1808
1809         // The quad face from bisection:
1810         tmpQFEdges[0] = otherEdgeIndex[1];
1811         tmpQFEdges[1] = newEdgeIndex[3];
1812         tmpQFEdges[2] = newEdgeIndex[0];
1813         tmpQFEdges[3] = newEdgeIndex[4];
1814         faceEdges_[newFaceIndex[3]] = tmpQFEdges;
1815
1816         // The quad interior face; Owner = c[1] & Neighbour = newCell[1]
1817         tmpQFEdges[0] = newEdgeIndex[0];
1818         tmpQFEdges[1] = newEdgeIndex[5];
1819         tmpQFEdges[2] = otherEdgeIndex[3];
1820         tmpQFEdges[3] = newEdgeIndex[6];
1821         faceEdges_[newFaceIndex[4]] = tmpQFEdges;
1822
1823         // Second boundary face; Owner = c[1] & Neighbour = [-1]
1824         tmpTFEdges[0] = commonEdgeIndex[2];
1825         tmpTFEdges[1] = cornerEdgeIndex[2];
1826         tmpTFEdges[2] = newEdgeIndex[5];
1827         faceEdges_[newFaceIndex[5]] = tmpTFEdges;
1828
1829         // Third boundary face; Owner = newCell[1] & Neighbour = [-1]
1830         tmpTFEdges[0] = newEdgeIndex[4];
1831         tmpTFEdges[1] = cornerEdgeIndex[3];
1832         tmpTFEdges[2] = newEdgeIndex[6];
1833         faceEdges_[newFaceIndex[6]] = tmpTFEdges;
1834
1835         meshOps::replaceLabel
1836         (
1837             commonEdgeIndex[1],
1838             newEdgeIndex[4],
1839             faceEdges_[c0BdyIndex[1]]
1840         );
1841
1842         meshOps::replaceLabel
1843         (
1844             c0BdyIndex[1],
1845             newFaceIndex[2],
1846             edgeFaces_[commonEdgeIndex[1]]
1847         );
1848
1849         meshOps::replaceLabel
1850         (
1851             commonEdgeIndex[2],
1852             newEdgeIndex[3],
1853             faceEdges_[commonFaceIndex[2]]
1854         );
1855
1856         meshOps::replaceLabel
1857         (
1858             commonFaceIndex[2],
1859             newFaceIndex[5],
1860             edgeFaces_[commonEdgeIndex[2]]
1861         );
1862
1863         if (debug > 2)
1864         {
1865             Pout<< nl << "Modified Cell[0]: "
1866                 << c0 << ": " << oldCells[0] << endl;
1867
1868             forAll(oldCells[0], faceI)
1869             {
1870                 const labelList& fE = faceEdges_[oldCells[0][faceI]];
1871
1872                 Pout<< oldCells[0][faceI]
1873                     << ": " << faces_[oldCells[0][faceI]]
1874                     << " fE: " << fE
1875                     << endl;
1876
1877                 forAll(fE, edgeI)
1878                 {
1879                     const labelList& eF = edgeFaces_[fE[edgeI]];
1880
1881                     Pout<< '\t' << fE[edgeI]
1882                         << ": " << edges_[fE[edgeI]]
1883                         << " eF: " << eF
1884                         << endl;
1885                 }
1886             }
1887
1888             Pout<< "New Cell[0]: "
1889                 << newCellIndex[0] << ": " << newCells[0] << endl;
1890
1891             forAll(newCells[0], faceI)
1892             {
1893                 const labelList& fE = faceEdges_[newCells[0][faceI]];
1894
1895                 Pout<< newCells[0][faceI] << ": "
1896                     << faces_[newCells[0][faceI]]
1897                     << " fE: " << fE
1898                     << endl;
1899
1900                 forAll(fE, edgeI)
1901                 {
1902                     const labelList& eF = edgeFaces_[fE[edgeI]];
1903
1904                     Pout<< '\t' << fE[edgeI]
1905                         << ": " << edges_[fE[edgeI]]
1906                         << " eF: " << eF
1907                         << endl;
1908                 }
1909             }
1910
1911             Pout<< nl << "Modified Cell[1]: "
1912                 << c1 << ": " << oldCells[1] << endl;
1913
1914             forAll(oldCells[1], faceI)
1915             {
1916                 const labelList& fE = faceEdges_[oldCells[1][faceI]];
1917
1918                 Pout<< oldCells[1][faceI] << ": "
1919                     << faces_[oldCells[1][faceI]]
1920                     << " fE: " << fE
1921                     << endl;
1922
1923                 forAll(fE, edgeI)
1924                 {
1925                     const labelList& eF = edgeFaces_[fE[edgeI]];
1926
1927                     Pout<< '\t' << fE[edgeI]
1928                         << ": " << edges_[fE[edgeI]]
1929                         << " eF: " << eF
1930                         << endl;
1931                 }
1932             }
1933
1934             Pout<< "New Cell[1]: "
1935                 << newCellIndex[1] << ": " << newCells[1] << endl;
1936
1937             forAll(newCells[1], faceI)
1938             {
1939                 const labelList& fE = faceEdges_[newCells[1][faceI]];
1940
1941                 Pout<< newCells[1][faceI] << ": "
1942                     << faces_[newCells[1][faceI]]
1943                     << " fE: " << fE
1944                     << endl;
1945
1946                 forAll(fE, edgeI)
1947                 {
1948                     const labelList& eF = edgeFaces_[fE[edgeI]];
1949
1950                     Pout<< '\t' << fE[edgeI]
1951                         << ": " << edges_[fE[edgeI]]
1952                         << " eF: " << eF
1953                         << endl;
1954                 }
1955             }
1956         }
1957
1958         // Update the cell list.
1959         cells_[c1] = oldCells[1];
1960         cells_[newCellIndex[1]] = newCells[1];
1961     }
1962
1963     // Update the cell list.
1964     cells_[c0] = oldCells[0];
1965     cells_[newCellIndex[0]] = newCells[0];
1966
1967     // Modify point labels for common edges
1968     if (edges_[commonEdgeIndex[0]].start() == otherEdgePoint[0])
1969     {
1970         edges_[commonEdgeIndex[0]].start() = newPointIndex[0];
1971     }
1972     else
1973     {
1974         edges_[commonEdgeIndex[0]].end() = newPointIndex[0];
1975     }
1976
1977     if (edges_[commonEdgeIndex[1]].start() == nextToOtherPoint[1])
1978     {
1979         edges_[commonEdgeIndex[1]].start() = newPointIndex[1];
1980     }
1981     else
1982     {
1983         edges_[commonEdgeIndex[1]].end() = newPointIndex[1];
1984     }
1985
1986     if (coupledModification_)
1987     {
1988         // Alias for convenience...
1989         const changeMap& slaveMap = slaveMaps[0];
1990
1991         const label pI = slaveMap.patchIndex();
1992
1993         // Fetch the appropriate coupleMap
1994         const coupleMap* cMapPtr = NULL;
1995
1996         if (localCouple && !procCouple)
1997         {
1998             cMapPtr = &(patchCoupling_[pI].map());
1999         }
2000         else
2001         if (procCouple && !localCouple)
2002         {
2003             cMapPtr = &(recvMeshes_[pI].map());
2004         }
2005
2006         // Alias for convenience
2007         const coupleMap& cMap = *cMapPtr;
2008
2009         // Add the new points to the coupling map
2010         const List<objectMap>& apList = slaveMap.addedPointList();
2011
2012         if (bisectingSlave)
2013         {
2014             // Update reverse pointMap
2015             cMap.mapMaster
2016             (
2017                 coupleMap::POINT,
2018                 newPointIndex[0],
2019                 apList[0].index()
2020             );
2021
2022             cMap.mapMaster
2023             (
2024                 coupleMap::POINT,
2025                 newPointIndex[1],
2026                 apList[1].index()
2027             );
2028
2029             // Update pointMap
2030             cMap.mapSlave
2031             (
2032                 coupleMap::POINT,
2033                 apList[0].index(),
2034                 newPointIndex[0]
2035             );
2036
2037             cMap.mapSlave
2038             (
2039                 coupleMap::POINT,
2040                 apList[1].index(),
2041                 newPointIndex[1]
2042             );
2043         }
2044         else
2045         {
2046             // Update pointMap
2047             cMap.mapSlave
2048             (
2049                 coupleMap::POINT,
2050                 newPointIndex[0],
2051                 apList[0].index()
2052             );
2053
2054             cMap.mapSlave
2055             (
2056                 coupleMap::POINT,
2057                 newPointIndex[1],
2058                 apList[1].index()
2059             );
2060
2061             // Update reverse pointMap
2062             cMap.mapMaster
2063             (
2064                 coupleMap::POINT,
2065                 apList[0].index(),
2066                 newPointIndex[0]
2067             );
2068
2069             cMap.mapMaster
2070             (
2071                 coupleMap::POINT,
2072                 apList[1].index(),
2073                 newPointIndex[1]
2074             );
2075         }
2076
2077         // Create a master/slave entry for the new face on the patch.
2078         FixedList<bool, 2> foundMatch(false);
2079         FixedList<label, 2> checkFaces(-1);
2080
2081         // Fill in the faces to check for...
2082         checkFaces[0] = fIndex;
2083         checkFaces[1] = newFaceIndex[3];
2084
2085         const List<objectMap>& afList = slaveMap.addedFaceList();
2086
2087         forAll (checkFaces, indexI)
2088         {
2089             const face& mFace = faces_[checkFaces[indexI]];
2090
2091             label sFaceIndex = -1;
2092
2093             forAll(afList, sfI)
2094             {
2095                 const face* facePtr = NULL;
2096
2097                 if (localCouple && !procCouple)
2098                 {
2099                     facePtr = &(faces_[afList[sfI].index()]);
2100                 }
2101                 else
2102                 if (procCouple && !localCouple)
2103                 {
2104                     const dynamicTopoFvMesh& mesh = recvMeshes_[pI].subMesh();
2105
2106                     facePtr = &(mesh.faces_[afList[sfI].index()]);
2107                 }
2108
2109                 const face& tFace = *facePtr;
2110                 FixedList<bool, 4> cP(false);
2111
2112                 forAll(mFace, pointI)
2113                 {
2114                     const Map<label>& pointMap =
2115                     (
2116                         cMap.entityMap(coupleMap::POINT)
2117                     );
2118
2119                     if (tFace.which(pointMap[mFace[pointI]]) > -1)
2120                     {
2121                         cP[pointI] = true;
2122                     }
2123                 }
2124
2125                 if (cP[0] && cP[1] && cP[2] && cP[3])
2126                 {
2127                     sFaceIndex = afList[sfI].index();
2128                     foundMatch[indexI] = true;
2129                     break;
2130                 }
2131             }
2132
2133             if (foundMatch[indexI])
2134             {
2135                 if (bisectingSlave)
2136                 {
2137                     cMap.mapMaster
2138                     (
2139                         coupleMap::FACE,
2140                         checkFaces[indexI],
2141                         sFaceIndex
2142                     );
2143
2144                     cMap.mapSlave
2145                     (
2146                         coupleMap::FACE,
2147                         sFaceIndex,
2148                         checkFaces[indexI]
2149                     );
2150                 }
2151                 else
2152                 {
2153                     cMap.mapSlave
2154                     (
2155                         coupleMap::FACE,
2156                         checkFaces[indexI],
2157                         sFaceIndex
2158                     );
2159
2160                     cMap.mapMaster
2161                     (
2162                         coupleMap::FACE,
2163                         sFaceIndex,
2164                         checkFaces[indexI]
2165                     );
2166                 }
2167             }
2168             else
2169             {
2170                 FatalErrorIn
2171                 (
2172                     "\n"
2173                     "const changeMap "
2174                     "dynamicTopoFvMesh::bisectQuadFace\n"
2175                     "(\n"
2176                     "    const label fIndex,\n"
2177                     "    const changeMap& masterMap,\n"
2178                     "    bool checkOnly,\n"
2179                     "    bool forceOp\n"
2180                     ")\n"
2181                 )
2182                     << "Failed to build coupled maps." << nl
2183                     << " masterFace: " << mFace << nl
2184                     << " Added slave faces: " << afList << nl
2185                     << abort(FatalError);
2186             }
2187         }
2188
2189         // Push operation into coupleMap
2190         cMap.pushOperation
2191         (
2192             slaveMap.index(),
2193             coupleMap::BISECTION
2194         );
2195     }
2196
2197     // Write out VTK files after change
2198     if (debug > 3)
2199     {
2200         labelList cellHull(4, -1);
2201
2202         cellHull[0] = owner_[fIndex];
2203         cellHull[1] = neighbour_[fIndex];
2204         cellHull[2] = owner_[newFaceIndex[3]];
2205         cellHull[3] = neighbour_[newFaceIndex[3]];
2206
2207         writeVTK
2208         (
2209             Foam::name(fIndex)
2210           + "_Bisect_1",
2211             cellHull
2212         );
2213     }
2214
2215     // Fill-in mapping information
2216     FixedList<label, 4> mapCells(-1);
2217
2218     mapCells[0] = c0;
2219     mapCells[1] = newCellIndex[0];
2220
2221     if (c1 != -1)
2222     {
2223         mapCells[2] = c1;
2224         mapCells[3] = newCellIndex[1];
2225     }
2226
2227     labelList mC(1, c0);
2228
2229     forAll(mapCells, cellI)
2230     {
2231         if (mapCells[cellI] == -1)
2232         {
2233             continue;
2234         }
2235
2236         // Set the mapping for this cell
2237         setCellMapping(mapCells[cellI], mC);
2238     }
2239
2240     // Set fill-in mapping information for the modified face.
2241     if (c1 == -1)
2242     {
2243         // Set the mapping for this face
2244         setFaceMapping(fIndex, labelList(1, fIndex));
2245     }
2246     else
2247     {
2248         // Internal face. Default mapping.
2249         setFaceMapping(fIndex);
2250     }
2251
2252     forAll(newFaceIndex, faceI)
2253     {
2254         if (newFaceIndex[faceI] == -1)
2255         {
2256             continue;
2257         }
2258
2259         // Check for boundary faces
2260         if (neighbour_[newFaceIndex[faceI]] == -1)
2261         {
2262             // Boundary face. Compute mapping.
2263             labelList mC;
2264
2265             if (faces_[newFaceIndex[faceI]].size() == 4)
2266             {
2267                 // Quad-face on boundary
2268                 mC.setSize(1, fIndex);
2269             }
2270             else
2271             if (faces_[newFaceIndex[faceI]].size() == 3)
2272             {
2273                 label triFacePatch = whichPatch(newFaceIndex[faceI]);
2274
2275                 // Fetch face-normals
2276                 vector tfNorm, f0Norm, f1Norm;
2277
2278                 tfNorm = faces_[newFaceIndex[faceI]].normal(oldPoints_);
2279                 f0Norm = faces_[c0BdyIndex[0]].normal(oldPoints_);
2280                 f1Norm = faces_[c0BdyIndex[1]].normal(oldPoints_);
2281
2282                 // Tri-face on boundary. Perform normal checks
2283                 // also, because of empty patches.
2284                 if
2285                 (
2286                     (whichPatch(c0BdyIndex[0]) == triFacePatch) &&
2287                     ((tfNorm & f0Norm) > 0.0)
2288                 )
2289                 {
2290                     mC.setSize(1, c0BdyIndex[0]);
2291                 }
2292                 else
2293                 if
2294                 (
2295                     (whichPatch(c0BdyIndex[1]) == triFacePatch) &&
2296                     ((tfNorm & f1Norm) > 0.0)
2297                 )
2298                 {
2299                     mC.setSize(1, c0BdyIndex[1]);
2300                 }
2301                 else
2302                 {
2303                     FatalErrorIn
2304                     (
2305                         "\n"
2306                         "const changeMap dynamicTopoFvMesh::bisectQuadFace\n"
2307                         "(\n"
2308                         "    const label fIndex,\n"
2309                         "    const changeMap& masterMap,\n"
2310                         "    bool checkOnly\n"
2311                         ")\n"
2312                     )
2313                         << " Unable to find patch for face: "
2314                         << newFaceIndex[faceI] << ":: "
2315                         << faces_[newFaceIndex[faceI]] << nl
2316                         << " Patch: " << triFacePatch << nl
2317                         << abort(FatalError);
2318                 }
2319             }
2320
2321             // Set the mapping for this face
2322             setFaceMapping(newFaceIndex[faceI], mC);
2323         }
2324         else
2325         {
2326             // Internal quad-faces get default mapping.
2327             setFaceMapping(newFaceIndex[faceI]);
2328         }
2329     }
2330
2331     // Set the flag
2332     topoChangeFlag_ = true;
2333
2334     // Increment the counter
2335     status(TOTAL_BISECTIONS)++;
2336
2337     // Increment surface-counter
2338     if (c1 == -1)
2339     {
2340         // Do not update stats for processor patches
2341         if (!processorCoupledEntity(fIndex))
2342         {
2343             status(SURFACE_BISECTIONS)++;
2344         }
2345     }
2346
2347     // Increment the number of modifications
2348     status(TOTAL_MODIFICATIONS)++;
2349
2350     // Specify that the operation was successful
2351     map.type() = 1;
2352
2353     // Return the changeMap
2354     return map;
2355 }
2356
2357
2358 // Method for the bisection of an edge in 3D
2359 // - Returns a changeMap with a type specifying:
2360 //     1: Bisection was successful
2361 //    -1: Bisection failed since max number of topo-changes was reached.
2362 //    -2: Bisection failed since resulting quality would be unacceptable.
2363 //    -3: Bisection failed since edge was on a noRefinement patch.
2364 // - AddedPoints contain the index of the newly added point.
2365 const changeMap dynamicTopoFvMesh::bisectEdge
2366 (
2367     const label eIndex,
2368     bool checkOnly,
2369     bool forceOp
2370 )
2371 {
2372     // Edge bisection performs the following operations:
2373     //      [1] Add a point at middle of the edge
2374     //      [2] Bisect all faces surrounding this edge
2375     //      [3] Bisect all cells surrounding this edge
2376     //      [4] Create internal/external edges for each bisected face
2377     //      [5] Create internal faces for each bisected cell
2378     //      Update faceEdges and edgeFaces information
2379
2380     // For 2D meshes, perform face-bisection
2381     if (is2D())
2382     {
2383         return bisectQuadFace(eIndex, changeMap(), checkOnly);
2384     }
2385
2386     // Figure out which thread this is...
2387     label tIndex = self();
2388
2389     // Prepare the changeMaps
2390     changeMap map;
2391     List<changeMap> slaveMaps;
2392     bool bisectingSlave = false;
2393
2394     if
2395     (
2396         (status(TOTAL_MODIFICATIONS) > maxModifications_) &&
2397         (maxModifications_ > -1)
2398     )
2399     {
2400         // Reached the max allowable topo-changes.
2401         stack(tIndex).clear();
2402
2403         return map;
2404     }
2405
2406     // Check if edgeRefinements are to be avoided on patch.
2407     const labelList& eF = edgeFaces_[eIndex];
2408
2409     forAll(eF, fI)
2410     {
2411         label fPatch = whichPatch(eF[fI]);
2412
2413         if (baseMesh_.lengthEstimator().checkRefinementPatch(fPatch))
2414         {
2415             map.type() = -3;
2416
2417             return map;
2418         }
2419     }
2420
2421     // Sanity check: Is the index legitimate?
2422     if (eIndex < 0)
2423     {
2424         FatalErrorIn
2425         (
2426             "\n"
2427             "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2428             "(\n"
2429             "    const label eIndex,\n"
2430             "    bool checkOnly,\n"
2431             "    bool forceOp\n"
2432             ")\n"
2433         )
2434             << " Invalid index: " << eIndex
2435             << abort(FatalError);
2436     }
2437
2438     // If coupled modification is set, and this is a
2439     // master edge, bisect its slaves as well.
2440     bool localCouple = false, procCouple = false;
2441
2442     if (coupledModification_)
2443     {
2444         const edge& eCheck = edges_[eIndex];
2445
2446         const label edgeEnum = coupleMap::EDGE;
2447
2448         // Is this a locally coupled edge (either master or slave)?
2449         if (locallyCoupledEntity(eIndex, true))
2450         {
2451             localCouple = true;
2452             procCouple = false;
2453         }
2454         else
2455         if (processorCoupledEntity(eIndex))
2456         {
2457             procCouple = true;
2458             localCouple = false;
2459         }
2460
2461         if (localCouple && !procCouple)
2462         {
2463             // Determine the slave index.
2464             label sIndex = -1, pIndex = -1;
2465
2466             forAll(patchCoupling_, patchI)
2467             {
2468                 if (patchCoupling_(patchI))
2469                 {
2470                     const coupleMap& cMap = patchCoupling_[patchI].map();
2471
2472                     if ((sIndex = cMap.findSlave(edgeEnum, eIndex)) > -1)
2473                     {
2474                         pIndex = patchI;
2475
2476                         break;
2477                     }
2478
2479                     // The following bit happens only during the sliver
2480                     // exudation process, since slave edges are
2481                     // usually not added to the coupled edge-stack.
2482                     if ((sIndex = cMap.findMaster(edgeEnum, eIndex)) > -1)
2483                     {
2484                         pIndex = patchI;
2485
2486                         // Notice that we are bisecting a slave edge first.
2487                         bisectingSlave = true;
2488
2489                         break;
2490                     }
2491                 }
2492             }
2493
2494             if (sIndex == -1)
2495             {
2496                 FatalErrorIn
2497                 (
2498                     "\n"
2499                     "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2500                     "(\n"
2501                     "    const label eIndex,\n"
2502                     "    bool checkOnly,\n"
2503                     "    bool forceOp\n"
2504                     ")\n"
2505                 )
2506                     << "Coupled maps were improperly specified." << nl
2507                     << " Slave index not found for: " << nl
2508                     << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2509                     << abort(FatalError);
2510             }
2511             else
2512             {
2513                 // If we've found the slave, size up the list
2514                 meshOps::sizeUpList
2515                 (
2516                     changeMap(),
2517                     slaveMaps
2518                 );
2519
2520                 // Save index and patch for posterity
2521                 slaveMaps[0].index() = sIndex;
2522                 slaveMaps[0].patchIndex() = pIndex;
2523             }
2524
2525             if (debug > 1)
2526             {
2527                 Pout<< nl << " >> Bisecting slave edge: " << sIndex
2528                     << " for master edge: " << eIndex << endl;
2529             }
2530         }
2531         else
2532         if (procCouple && !localCouple)
2533         {
2534             // If this is a new entity, bail out for now.
2535             // This will be handled at the next time-step.
2536             if (eIndex >= nOldEdges_)
2537             {
2538                 return map;
2539             }
2540
2541             // Check slaves
2542             forAll(procIndices_, pI)
2543             {
2544                 // Fetch reference to subMeshes
2545                 const coupledMesh& sendMesh = sendMeshes_[pI];
2546                 const coupledMesh& recvMesh = recvMeshes_[pI];
2547
2548                 const coupleMap& scMap = sendMesh.map();
2549                 const coupleMap& rcMap = recvMesh.map();
2550
2551                 // If this edge was sent to a lower-ranked
2552                 // processor, skip it.
2553                 if
2554                 (
2555                     priority
2556                     (
2557                         procIndices_[pI],
2558                         lessOp<label>(),
2559                         Pstream::myProcNo()
2560                     )
2561                 )
2562                 {
2563                     if (scMap.reverseEntityMap(edgeEnum).found(eIndex))
2564                     {
2565                         if (debug > 3)
2566                         {
2567                             Pout<< "Edge: " << eIndex
2568                                 << "::" << eCheck
2569                                 << " was sent to proc: "
2570                                 << procIndices_[pI]
2571                                 << ", so bailing out."
2572                                 << endl;
2573                         }
2574
2575                         return map;
2576                     }
2577                 }
2578
2579                 label sIndex = -1;
2580
2581                 if ((sIndex = rcMap.findSlave(edgeEnum, eIndex)) > -1)
2582                 {
2583                     if (debug > 3)
2584                     {
2585                         Pout<< "Checking slave edge: " << sIndex
2586                             << " on proc: " << procIndices_[pI]
2587                             << " for master edge: " << eIndex
2588                             << endl;
2589                     }
2590
2591                     // Check if a lower-ranked processor is
2592                     // handling this edge
2593                     if
2594                     (
2595                         priority
2596                         (
2597                             procIndices_[pI],
2598                             lessOp<label>(),
2599                             Pstream::myProcNo()
2600                         )
2601                     )
2602                     {
2603                         if (debug > 3)
2604                         {
2605                             Pout<< "Edge: " << eIndex
2606                                 << " is handled by proc: "
2607                                 << procIndices_[pI]
2608                                 << ", so bailing out."
2609                                 << endl;
2610                         }
2611
2612                         return map;
2613                     }
2614
2615                     label curIndex = slaveMaps.size();
2616
2617                     // Size up the list
2618                     meshOps::sizeUpList
2619                     (
2620                         changeMap(),
2621                         slaveMaps
2622                     );
2623
2624                     // Save index and patch for posterity
2625                     slaveMaps[curIndex].index() = sIndex;
2626                     slaveMaps[curIndex].patchIndex() = pI;
2627                 }
2628             }
2629         }
2630         else
2631         {
2632             // Something's wrong with coupling maps
2633             FatalErrorIn
2634             (
2635                 "\n"
2636                 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2637                 "(\n"
2638                 "    const label eIndex,\n"
2639                 "    bool checkOnly,\n"
2640                 "    bool forceOp\n"
2641                 ")\n"
2642             )
2643                 << "Coupled maps were improperly specified." << nl
2644                 << " localCouple: " << localCouple << nl
2645                 << " procCouple: " << procCouple << nl
2646                 << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2647                 << abort(FatalError);
2648         }
2649
2650         // Temporarily turn off coupledModification
2651         unsetCoupledModification();
2652
2653         // First check the master for bisection feasibility.
2654         changeMap masterMap = bisectEdge(eIndex, true, forceOp);
2655
2656         // Turn it back on.
2657         setCoupledModification();
2658
2659         // Master couldn't perform bisection
2660         if (masterMap.type() <= 0)
2661         {
2662             return masterMap;
2663         }
2664
2665         // Now check each of the slaves for bisection feasibility
2666         forAll(slaveMaps, slaveI)
2667         {
2668             // Alias for convenience...
2669             changeMap& slaveMap = slaveMaps[slaveI];
2670
2671             label sIndex = slaveMap.index();
2672             label pI = slaveMap.patchIndex();
2673
2674             // If the edge is mapped onto itself, skip check
2675             // (occurs for cyclic edges)
2676             if ((sIndex == eIndex) && localCouple)
2677             {
2678                 continue;
2679             }
2680
2681             // Temporarily turn off coupledModification
2682             unsetCoupledModification();
2683
2684             // Test the slave edge
2685             if (localCouple)
2686             {
2687                 slaveMap = bisectEdge(sIndex, true, forceOp);
2688             }
2689             else
2690             if (procCouple)
2691             {
2692                 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
2693
2694                 if (debug > 3)
2695                 {
2696                     Pout<< "Checking slave edge: " << sIndex
2697                         << "::" << sMesh.edges_[sIndex]
2698                         << " on proc: " << procIndices_[pI]
2699                         << " for master edge: " << eIndex
2700                         << endl;
2701                 }
2702
2703                 slaveMap = sMesh.bisectEdge(sIndex, true, forceOp);
2704             }
2705
2706             // Turn it back on.
2707             setCoupledModification();
2708
2709             if (slaveMap.type() <= 0)
2710             {
2711                 // Slave couldn't perform bisection.
2712                 map.type() = -2;
2713
2714                 return map;
2715             }
2716
2717             // Save index and patch for posterity
2718             slaveMap.index() = sIndex;
2719             slaveMap.patchIndex() = pI;
2720         }
2721
2722         // Next bisect each slave edge..
2723         forAll(slaveMaps, slaveI)
2724         {
2725             // Alias for convenience...
2726             changeMap& slaveMap = slaveMaps[slaveI];
2727
2728             label sIndex = slaveMap.index();
2729             label pI = slaveMap.patchIndex();
2730
2731             // If the edge is mapped onto itself, skip modification
2732             // (occurs for cyclic edges)
2733             if ((sIndex == eIndex) && localCouple)
2734             {
2735                 continue;
2736             }
2737
2738             // Temporarily turn off coupledModification
2739             unsetCoupledModification();
2740
2741             // Bisect the slave edge
2742             if (localCouple)
2743             {
2744                 slaveMap = bisectEdge(sIndex, false, forceOp);
2745             }
2746             else
2747             {
2748                 dynamicTopoFvMesh& sMesh = recvMeshes_[pI].subMesh();
2749
2750                 slaveMap = sMesh.bisectEdge(sIndex, false, forceOp);
2751             }
2752
2753             // Turn it back on.
2754             setCoupledModification();
2755
2756             // The final operation has to succeed.
2757             if (slaveMap.type() <= 0)
2758             {
2759                 FatalErrorIn
2760                 (
2761                     "\n"
2762                     "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2763                     "(\n"
2764                     "    const label eIndex,\n"
2765                     "    bool checkOnly,\n"
2766                     "    bool forceOp\n"
2767                     ")\n"
2768                 )
2769                     << "Coupled topo-change for slave failed." << nl
2770                     << " Edge: " << eIndex << ": " << edges_[eIndex] << nl
2771                     << " Slave index: " << sIndex << nl
2772                     << " Patch index: " << pI << nl
2773                     << " Type: " << slaveMap.type() << nl
2774                     << abort(FatalError);
2775             }
2776
2777             // Save index and patch for posterity
2778             slaveMap.index() = sIndex;
2779             slaveMap.patchIndex() = pI;
2780         }
2781     }
2782
2783     // Before we bisect this edge, check whether the operation will
2784     // yield an acceptable cell-quality.
2785     scalar minQ = 0.0;
2786
2787     if ((minQ = computeBisectionQuality(eIndex)) < sliverThreshold_)
2788     {
2789         // Check if the quality is actually valid before forcing it.
2790         if (forceOp && (minQ < 0.0))
2791         {
2792             FatalErrorIn
2793             (
2794                 "\n"
2795                 "const changeMap dynamicTopoFvMesh::bisectEdge\n"
2796                 "(\n"
2797                 "    const label eIndex,\n"
2798                 "    bool checkOnly,\n"
2799                 "    bool forceOp\n"
2800                 ")\n"
2801             )
2802                 << " Forcing bisection on edge: " << eIndex
2803                 << " will yield an invalid cell."
2804                 << abort(FatalError);
2805         }
2806         else
2807         if (!forceOp)
2808         {
2809             map.type() = -2;
2810             return map;
2811         }
2812     }
2813
2814     // Are we performing only checks?
2815     if (checkOnly)
2816     {
2817         if (debug > 3 && isSubMesh_)
2818         {
2819             Pout<< "  Slave edge: " << eIndex
2820                 << " can be bisected."
2821                 << endl;
2822         }
2823
2824         map.type() = 1;
2825         return map;
2826     }
2827
2828     // Update number of surface bisections, if necessary.
2829     if (whichEdgePatch(eIndex) > -1)
2830     {
2831         // Do not update stats for processor patches
2832         if (!processorCoupledEntity(eIndex))
2833         {
2834             status(SURFACE_BISECTIONS)++;
2835         }
2836     }
2837
2838     // Hull variables
2839     face tmpTriFace(3);
2840     edge origEdge(edges_[eIndex]);
2841     labelList tmpEdgeFaces(3,-1);
2842     labelList tmpIntEdgeFaces(4,-1);
2843     labelList tmpFaceEdges(3,-1);
2844
2845     // Build vertexHull for this edge
2846     labelList vertexHull;
2847     buildVertexHull(eIndex, vertexHull);
2848
2849     // Size up the hull lists
2850     label m = vertexHull.size();
2851     labelList cellHull(m, -1);
2852     labelList faceHull(m, -1);
2853     labelList edgeHull(m, -1);
2854     labelListList ringEntities(4, labelList(m, -1));
2855
2856     // Construct a hull around this edge
2857     meshOps::constructHull
2858     (
2859         eIndex,
2860         faces_,
2861         edges_,
2862         cells_,
2863         owner_,
2864         neighbour_,
2865         faceEdges_,
2866         edgeFaces_,
2867         vertexHull,
2868         edgeHull,
2869         faceHull,
2870         cellHull,
2871         ringEntities
2872     );
2873
2874     if (debug > 1)
2875     {
2876         Pout<< nl << nl
2877             << "Edge: " << eIndex
2878             << ": " << origEdge
2879             << " is to be bisected. " << endl;
2880
2881         Pout<< " On SubMesh: " << isSubMesh_ << nl;
2882         Pout<< " coupledModification: " << coupledModification_ << nl;
2883
2884         label epIndex = whichEdgePatch(eIndex);
2885
2886         const polyBoundaryMesh& boundary = boundaryMesh();
2887
2888         Pout<< " Patch: ";
2889
2890         if (epIndex == -1)
2891         {
2892             Pout<< "Internal" << endl;
2893         }
2894         else
2895         if (epIndex < boundary.size())
2896         {
2897             Pout<< boundary[epIndex].name() << endl;
2898         }
2899         else
2900         {
2901             Pout<< " New patch: " << epIndex << endl;
2902         }
2903
2904         // Write out VTK files prior to change
2905         //  - Using old-points for convenience in post-processing
2906         if (debug > 3)
2907         {
2908             writeVTK
2909             (
2910                 Foam::name(eIndex)
2911               + '(' + Foam::name(origEdge[0])
2912               + ',' + Foam::name(origEdge[1]) + ')'
2913               + "_Bisect_0",
2914                 cellHull,
2915                 3, false, true
2916             );
2917         }
2918     }
2919
2920     labelList mP(2, -1);
2921
2922     // Set mapping for this point
2923     mP[0] = edges_[eIndex][0];
2924     mP[1] = edges_[eIndex][1];
2925
2926     // Add a new point to the end of the list
2927     label newPointIndex =
2928     (
2929         insertPoint
2930         (
2931             0.5 * (points_[mP[0]] + points_[mP[1]]),
2932             0.5 * (oldPoints_[mP[0]] + oldPoints_[mP[1]]),
2933             mP
2934         )
2935     );
2936
2937     // Add this point to the map.
2938     map.addPoint(newPointIndex);
2939
2940     // New edges can lie on a bounding curve between
2941     // coupled and non-coupled faces. Preferentially
2942     // add edges to coupled-patches, if they exist.
2943     // This makes it convenient for coupled patch matching.
2944     label nePatch = -1;
2945
2946     if (locallyCoupledEntity(eIndex, true))
2947     {
2948         nePatch = locallyCoupledEdgePatch(eIndex);
2949     }
2950     else
2951     {
2952         nePatch = whichEdgePatch(eIndex);
2953     }
2954
2955     // Add a new edge to the end of the list
2956     label newEdgeIndex =
2957     (
2958         insertEdge
2959         (
2960             nePatch,
2961             edge(newPointIndex,edges_[eIndex][1]),
2962             labelList(faceHull.size(),-1)
2963         )
2964     );
2965
2966     // Add this edge to the map.
2967     map.addEdge(newEdgeIndex);
2968
2969     // Remove the existing edge from the pointEdges list
2970     // of the modified point, and add it to the new point
2971     meshOps::sizeDownList(eIndex, pointEdges_[edges_[eIndex][1]]);
2972     meshOps::sizeUpList(eIndex, pointEdges_[newPointIndex]);
2973
2974     // Modify the existing edge
2975     edges_[eIndex][1] = newPointIndex;
2976
2977     // Add this edge to the map.
2978     // Although this edge isn't technically 'added', it's
2979     // required for coupled patch mapping.
2980     map.addEdge(eIndex);
2981
2982     // Keep track of added entities
2983     labelList addedCellIndices(cellHull.size(),-1);
2984     labelList addedFaceIndices(faceHull.size(),-1);
2985     labelList addedEdgeIndices(faceHull.size(),-1);
2986     labelList addedIntFaceIndices(faceHull.size(),-1);
2987
2988     // Now loop through the hull and bisect individual entities
2989     forAll(vertexHull, indexI)
2990     {
2991         // Modify the existing face
2992         meshOps::replaceLabel
2993         (
2994             edges_[newEdgeIndex][1],
2995             newPointIndex,
2996             faces_[faceHull[indexI]]
2997         );
2998
2999         // Obtain circular indices
3000         label nextI = vertexHull.fcIndex(indexI);
3001         label prevI = vertexHull.rcIndex(indexI);
3002
3003         // Check if this is an interior/boundary face
3004         if (cellHull[indexI] != -1)
3005         {
3006             // Create a new cell. Add it for now, but update later.
3007             cell newCell(4);
3008
3009             addedCellIndices[indexI] =
3010             (
3011                 insertCell
3012                 (
3013                     newCell,
3014                     lengthScale_[cellHull[indexI]]
3015                 )
3016             );
3017
3018             // Add this cell to the map.
3019             map.addCell(addedCellIndices[indexI]);
3020
3021             // Configure the interior face
3022             tmpTriFace[0] = vertexHull[nextI];
3023             tmpTriFace[1] = vertexHull[indexI];
3024             tmpTriFace[2] = newPointIndex;
3025
3026             // Insert the face
3027             addedIntFaceIndices[indexI] =
3028             (
3029                 insertFace
3030                 (
3031                     -1,
3032                     tmpTriFace,
3033                     cellHull[indexI],
3034                     addedCellIndices[indexI],
3035                     tmpFaceEdges
3036                 )
3037             );
3038
3039             // Add this face to the map.
3040             map.addFace(addedIntFaceIndices[indexI]);
3041
3042             // Add to the new cell
3043             newCell[0] = addedIntFaceIndices[indexI];
3044
3045             // Modify the existing ring face connected to newEdge[1]
3046             label replaceFace = ringEntities[3][indexI];
3047
3048             // Check if face reversal is necessary
3049             if (owner_[replaceFace] == cellHull[indexI])
3050             {
3051                 if (neighbour_[replaceFace] == -1)
3052                 {
3053                     // Change the owner
3054                     owner_[replaceFace] = addedCellIndices[indexI];
3055                 }
3056                 else
3057                 {
3058                     // This face has to be reversed
3059                     faces_[replaceFace] = faces_[replaceFace].reverseFace();
3060                     owner_[replaceFace] = neighbour_[replaceFace];
3061                     neighbour_[replaceFace] = addedCellIndices[indexI];
3062
3063                     setFlip(replaceFace);
3064                 }
3065             }
3066             else
3067             {
3068                 // Keep owner, but change neighbour
3069                 neighbour_[replaceFace] = addedCellIndices[indexI];
3070             }
3071
3072             // Modify the edge on the ring.
3073             // Add the new interior face to edgeFaces.
3074             meshOps::sizeUpList
3075             (
3076                 addedIntFaceIndices[indexI],
3077                 edgeFaces_[edgeHull[indexI]]
3078             );
3079
3080             // Add this edge to faceEdges for the new interior face
3081             faceEdges_[addedIntFaceIndices[indexI]][0] = edgeHull[indexI];
3082
3083             // Replace face labels
3084             meshOps::replaceLabel
3085             (
3086                 replaceFace,
3087                 addedIntFaceIndices[indexI],
3088                 cells_[cellHull[indexI]]
3089             );
3090
3091             // Add to the new cell
3092             newCell[1] = replaceFace;
3093
3094             // Check if this is a boundary face
3095             if (cellHull[prevI] == -1)
3096             {
3097                 // Configure the boundary face
3098                 tmpTriFace[0] = newPointIndex;
3099                 tmpTriFace[1] = edges_[newEdgeIndex][1];
3100                 tmpTriFace[2] = vertexHull[indexI];
3101
3102                 // Insert the face
3103                 addedFaceIndices[indexI] =
3104                 (
3105                     insertFace
3106                     (
3107                         whichPatch(faceHull[indexI]),
3108                         tmpTriFace,
3109                         addedCellIndices[indexI],
3110                         -1,
3111                         labelList(3, -1)
3112                     )
3113                 );
3114
3115                 // Add this face to the map.
3116                 map.addFace(addedFaceIndices[indexI]);
3117
3118                 // Configure edgeFaces
3119                 tmpEdgeFaces[0] = faceHull[indexI];
3120                 tmpEdgeFaces[1] = addedIntFaceIndices[indexI];
3121                 tmpEdgeFaces[2] = addedFaceIndices[indexI];
3122
3123                 // Add an edge
3124                 addedEdgeIndices[indexI] =
3125                 (
3126                     insertEdge
3127                     (
3128                         whichPatch(faceHull[indexI]),
3129                         edge(newPointIndex,vertexHull[indexI]),
3130                         tmpEdgeFaces
3131                     )
3132                 );
3133
3134                 // Add this edge to the map.
3135                 map.addEdge(addedEdgeIndices[indexI]);
3136
3137                 // Add this edge to the interior-face faceEdges entry
3138                 faceEdges_[addedIntFaceIndices[indexI]][1] =
3139                 (
3140                     addedEdgeIndices[indexI]
3141                 );
3142
3143                 // Configure faceEdges for this boundary face
3144                 tmpFaceEdges[0] = addedEdgeIndices[indexI];
3145                 tmpFaceEdges[1] = newEdgeIndex;
3146                 tmpFaceEdges[2] = ringEntities[2][indexI];
3147
3148                 // Modify faceEdges for the hull face
3149                 meshOps::replaceLabel
3150                 (
3151                     ringEntities[2][indexI],
3152                     addedEdgeIndices[indexI],
3153                     faceEdges_[faceHull[indexI]]
3154                 );
3155
3156                 // Modify edgeFaces for the edge connected to newEdge[1]
3157                 meshOps::replaceLabel
3158                 (
3159                     faceHull[indexI],
3160                     addedFaceIndices[indexI],
3161                     edgeFaces_[ringEntities[2][indexI]]
3162                 );
3163
3164                 // Add the faceEdges entry
3165                 faceEdges_[addedFaceIndices[indexI]] = tmpFaceEdges;
3166
3167                 // Add an entry to edgeFaces
3168                 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3169
3170                 // Add an entry for this cell
3171                 newCell[2] = addedFaceIndices[indexI];
3172             }
3173             else
3174             // Check if a cell was added before this
3175             if (addedCellIndices[prevI] != -1)
3176             {
3177                 // Configure the interior face
3178                 tmpTriFace[0] = vertexHull[indexI];
3179                 tmpTriFace[1] = edges_[newEdgeIndex][1];
3180                 tmpTriFace[2] = newPointIndex;
3181
3182                 // Insert the face
3183                 addedFaceIndices[indexI] =
3184                 (
3185                     insertFace
3186                     (
3187                         -1,
3188                         tmpTriFace,
3189                         addedCellIndices[prevI],
3190                         addedCellIndices[indexI],
3191                         labelList(3, -1)
3192                     )
3193                 );
3194
3195                 // Add this face to the map.
3196                 map.addFace(addedFaceIndices[indexI]);
3197
3198                 // Configure edgeFaces
3199                 tmpIntEdgeFaces[0] = faceHull[indexI];
3200                 tmpIntEdgeFaces[1] = addedIntFaceIndices[indexI];
3201                 tmpIntEdgeFaces[2] = addedFaceIndices[indexI];
3202                 tmpIntEdgeFaces[3] = addedIntFaceIndices[prevI];
3203
3204                 // Add an internal edge
3205                 addedEdgeIndices[indexI] =
3206                 (
3207                     insertEdge
3208                     (
3209                         -1,
3210                         edge(newPointIndex,vertexHull[indexI]),
3211                         tmpIntEdgeFaces
3212                     )
3213                 );
3214
3215                 // Add this edge to the map.
3216                 map.addEdge(addedEdgeIndices[indexI]);
3217
3218                 // Add this edge to the interior-face faceEdges entry..
3219                 faceEdges_[addedIntFaceIndices[indexI]][1] =
3220                 (
3221                     addedEdgeIndices[indexI]
3222                 );
3223
3224                 // ... and to the previous interior face as well
3225                 faceEdges_[addedIntFaceIndices[prevI]][2] =
3226                 (
3227                     addedEdgeIndices[indexI]
3228                 );
3229
3230                 // Configure faceEdges for this split interior face
3231                 tmpFaceEdges[0] = addedEdgeIndices[indexI];
3232                 tmpFaceEdges[1] = newEdgeIndex;
3233                 tmpFaceEdges[2] = ringEntities[2][indexI];
3234
3235                 // Modify faceEdges for the hull face
3236                 meshOps::replaceLabel
3237                 (
3238                     ringEntities[2][indexI],
3239                     addedEdgeIndices[indexI],
3240                     faceEdges_[faceHull[indexI]]
3241                 );
3242
3243                 // Modify edgeFaces for the edge connected to newEdge[1]
3244                 meshOps::replaceLabel
3245                 (
3246                     faceHull[indexI],
3247                     addedFaceIndices[indexI],
3248                     edgeFaces_[ringEntities[2][indexI]]
3249                 );
3250
3251                 // Add the faceEdges entry
3252                 faceEdges_[addedFaceIndices[indexI]] = tmpFaceEdges;
3253
3254                 // Add an entry to edgeFaces
3255                 edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3256
3257                 // Add an entry for this cell
3258                 newCell[2] = addedFaceIndices[indexI];
3259
3260                 // Make the final entry for the previous cell
3261                 cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
3262             }
3263
3264             // Do the first interior face at the end
3265             if (indexI == vertexHull.size() - 1)
3266             {
3267                 // Configure the interior face
3268                 tmpTriFace[0] = newPointIndex;
3269                 tmpTriFace[1] = edges_[newEdgeIndex][1];
3270                 tmpTriFace[2] = vertexHull[0];
3271
3272                 // Insert the face
3273                 addedFaceIndices[0] =
3274                 (
3275                     insertFace
3276                     (
3277                         -1,
3278                         tmpTriFace,
3279                         addedCellIndices[0],
3280                         addedCellIndices[indexI],
3281                         labelList(3, -1)
3282                     )
3283                 );
3284
3285                 // Add this face to the map.
3286                 map.addFace(addedFaceIndices[0]);
3287
3288                 // Configure edgeFaces
3289                 tmpIntEdgeFaces[0] = faceHull[0];
3290                 tmpIntEdgeFaces[1] = addedIntFaceIndices[0];
3291                 tmpIntEdgeFaces[2] = addedFaceIndices[0];
3292                 tmpIntEdgeFaces[3] = addedIntFaceIndices[indexI];
3293
3294                 // Add an internal edge
3295                 addedEdgeIndices[0] =
3296                 (
3297                     insertEdge
3298                     (
3299                         -1,
3300                         edge(newPointIndex,vertexHull[0]),
3301                         tmpIntEdgeFaces
3302                     )
3303                 );
3304
3305                 // Add this edge to the map.
3306                 map.addEdge(addedEdgeIndices[0]);
3307
3308                 // Add this edge to the interior-face faceEdges entry..
3309                 faceEdges_[addedIntFaceIndices[0]][1] =
3310                 (
3311                     addedEdgeIndices[0]
3312                 );
3313
3314                 // ... and to the previous interior face as well
3315                 faceEdges_[addedIntFaceIndices[indexI]][2] =
3316                 (
3317                     addedEdgeIndices[0]
3318                 );
3319
3320                 // Configure faceEdges for the first split face
3321                 tmpFaceEdges[0] = addedEdgeIndices[0];
3322                 tmpFaceEdges[1] = newEdgeIndex;
3323                 tmpFaceEdges[2] = ringEntities[2][0];
3324
3325                 // Modify faceEdges for the hull face
3326                 meshOps::replaceLabel
3327                 (
3328                     ringEntities[2][0],
3329                     addedEdgeIndices[0],
3330                     faceEdges_[faceHull[0]]
3331                 );
3332
3333                 // Modify edgeFaces for the edge connected to newEdge[1]
3334                 meshOps::replaceLabel
3335                 (
3336                     faceHull[0],
3337                     addedFaceIndices[0],
3338                     edgeFaces_[ringEntities[2][0]]
3339                 );
3340
3341                 // Add the faceEdges entry
3342                 faceEdges_[addedFaceIndices[0]] = tmpFaceEdges;
3343
3344                 // Add an entry to edgeFaces
3345                 edgeFaces_[newEdgeIndex][0] = addedFaceIndices[0];
3346
3347                 // Add an entry for this cell
3348                 newCell[3] = addedFaceIndices[0];
3349
3350                 // Make the final entry for the first cell
3351                 cells_[addedCellIndices[0]][2] = addedFaceIndices[0];
3352             }
3353
3354             // Update the cell list with the new cell.
3355             cells_[addedCellIndices[indexI]] = newCell;
3356         }
3357         else
3358         {
3359             // Configure the final boundary face
3360             tmpTriFace[0] = vertexHull[indexI];
3361             tmpTriFace[1] = edges_[newEdgeIndex][1];
3362             tmpTriFace[2] = newPointIndex;
3363
3364             // Insert the face
3365             addedFaceIndices[indexI] =
3366             (
3367                 insertFace
3368                 (
3369                     whichPatch(faceHull[indexI]),
3370                     tmpTriFace,
3371                     addedCellIndices[prevI],
3372                     -1,
3373                     labelList(3, -1)
3374                 )
3375             );
3376
3377             // Add this face to the map.
3378             map.addFace(addedFaceIndices[indexI]);
3379
3380             // Configure edgeFaces
3381             tmpEdgeFaces[0] = addedFaceIndices[indexI];
3382             tmpEdgeFaces[1] = addedIntFaceIndices[prevI];
3383             tmpEdgeFaces[2] = faceHull[indexI];
3384
3385             // Add an edge
3386             addedEdgeIndices[indexI] =
3387             (
3388                 insertEdge
3389                 (
3390                     whichPatch(faceHull[indexI]),
3391                     edge(newPointIndex,vertexHull[indexI]),
3392                     tmpEdgeFaces
3393                 )
3394             );
3395
3396             // Add this edge to the map.
3397             map.addEdge(addedEdgeIndices[indexI]);
3398
3399             // Add a faceEdges entry to the previous interior face
3400             faceEdges_[addedIntFaceIndices[prevI]][2] =
3401             (
3402                 addedEdgeIndices[indexI]
3403             );
3404
3405             // Configure faceEdges for the final boundary face
3406             tmpFaceEdges[0] = addedEdgeIndices[indexI];
3407             tmpFaceEdges[1] = newEdgeIndex;
3408             tmpFaceEdges[2] = ringEntities[2][indexI];
3409
3410             // Modify faceEdges for the hull face
3411             meshOps::replaceLabel
3412             (
3413                 ringEntities[2][indexI],
3414                 addedEdgeIndices[indexI],
3415                 faceEdges_[faceHull[indexI]]
3416             );
3417
3418             // Modify edgeFaces for the edge connected to newEdge[1]
3419             meshOps::replaceLabel
3420             (
3421                 faceHull[indexI],
3422                 addedFaceIndices[indexI],
3423                 edgeFaces_[ringEntities[2][indexI]]
3424             );
3425
3426             // Add the faceEdges entry
3427             faceEdges_[addedFaceIndices[indexI]] = tmpFaceEdges;
3428
3429             // Add an entry to edgeFaces
3430             edgeFaces_[newEdgeIndex][indexI] = addedFaceIndices[indexI];
3431
3432             // Make the final entry for the previous cell
3433             cells_[addedCellIndices[prevI]][3] = addedFaceIndices[indexI];
3434         }
3435     }
3436
3437     // Now that all old / new cells possess correct connectivity,
3438     // compute mapping information.
3439     forAll(cellHull, indexI)
3440     {
3441         if (cellHull[indexI] == -1)
3442         {
3443             continue;
3444         }
3445
3446         // Set mapping for both new and modified cells.
3447         FixedList<label, 2> cmIndex;
3448
3449         cmIndex[0] = cellHull[indexI];
3450         cmIndex[1] = addedCellIndices[indexI];
3451
3452         // Fill-in candidate mapping information
3453         labelList mC(1, cellHull[indexI]);
3454
3455         forAll(cmIndex, cmI)
3456         {
3457             // Set the mapping for this cell
3458             setCellMapping(cmIndex[cmI], mC);
3459         }
3460     }
3461
3462     // Set mapping information for old / new faces
3463     forAll(faceHull, indexI)
3464     {
3465         // Interior faces get default mapping
3466         if (addedIntFaceIndices[indexI] > -1)
3467         {
3468             setFaceMapping(addedIntFaceIndices[indexI]);
3469         }
3470
3471         // Decide between default / weighted mapping
3472         // based on boundary information
3473         if (whichPatch(faceHull[indexI]) == -1)
3474         {
3475             // Interior faces get default mapping
3476             setFaceMapping(faceHull[indexI]);
3477             setFaceMapping(addedFaceIndices[indexI]);
3478         }
3479         else
3480         {
3481             // Compute mapping weights for boundary faces
3482             FixedList<label, 2> fmIndex;
3483
3484             fmIndex[0] = faceHull[indexI];
3485             fmIndex[1] = addedFaceIndices[indexI];
3486
3487             // Fill-in candidate mapping information
3488             labelList mF(1, faceHull[indexI]);
3489
3490             forAll(fmIndex, fmI)
3491             {
3492                 // Set the mapping for this face
3493                 setFaceMapping(fmIndex[fmI], mF);
3494             }
3495         }
3496     }
3497
3498     // If modification is coupled, update mapping info.
3499     if (coupledModification_)
3500     {
3501         // Build a list of boundary edges / faces for mapping
3502         dynamicLabelList checkEdges(8), checkFaces(4);
3503
3504         const labelList& oeFaces = edgeFaces_[eIndex];
3505         const labelList& neFaces = edgeFaces_[newEdgeIndex];
3506
3507         forAll(oeFaces, faceI)
3508         {
3509             FixedList<label, 2> check(-1);
3510
3511             check[0] = oeFaces[faceI];
3512             check[1] = neFaces[faceI];
3513
3514             forAll(check, indexI)
3515             {
3516                 label cPatch = whichPatch(check[indexI]);
3517
3518                 if (localCouple && !procCouple)
3519                 {
3520                     if (!locallyCoupledEntity(check[indexI], true, false, true))
3521                     {
3522                         continue;
3523                     }
3524
3525                     // Check for cyclics
3526                     if (boundaryMesh()[cPatch].type() == "cyclic")
3527                     {
3528                         // Check if this is a master face
3529                         const coupleMap& cM = patchCoupling_[cPatch].map();
3530                         const Map<label>& fM = cM.entityMap(coupleMap::FACE);
3531
3532                         // Only add master faces
3533                         // (belonging to the first half)
3534                         //  - Only check[0] will be present,
3535                         //    so check for that alone
3536                         if (!fM.found(check[0]))
3537                         {
3538                             continue;
3539                         }
3540                     }
3541                 }
3542                 else
3543                 if (procCouple && !localCouple)
3544                 {
3545                     if (getNeighbourProcessor(cPatch) == -1)
3546                     {
3547                         continue;
3548                     }
3549                 }
3550
3551                 // Add face and its edges for checking
3552                 if (findIndex(checkFaces, check[indexI]) == -1)
3553                 {
3554                     // Add this face
3555                     checkFaces.append(check[indexI]);
3556
3557                     const labelList& fEdges = faceEdges_[check[indexI]];
3558
3559                     forAll(fEdges, edgeI)
3560                     {
3561                         if (findIndex(checkEdges, fEdges[edgeI]) == -1)
3562                         {
3563                             checkEdges.append(fEdges[edgeI]);
3564                         }
3565                     }
3566                 }
3567             }
3568         }
3569
3570         // Prepare a checklist
3571         boolList matchedFaces(checkFaces.size(), false);
3572         boolList matchedEdges(checkEdges.size(), false);
3573
3574         // Output check entities
3575         if (debug > 4)
3576         {
3577             writeVTK
3578             (
3579                 "checkEdges_" + Foam::name(eIndex),
3580                 checkEdges, 1, false, true
3581             );
3582
3583             writeVTK
3584             (
3585                 "checkFaces_" + Foam::name(eIndex),
3586                 checkFaces, 2, false, true
3587             );
3588         }
3589
3590         forAll(slaveMaps, slaveI)
3591         {
3592             // Alias for convenience...
3593             const changeMap& slaveMap = slaveMaps[slaveI];
3594
3595             const label pI = slaveMap.patchIndex();
3596
3597             // Fetch the appropriate coupleMap / mesh
3598             const coupleMap* cMapPtr = NULL;
3599             const dynamicTopoFvMesh* sMeshPtr = NULL;
3600
3601             if (localCouple && !procCouple)
3602             {
3603                 cMapPtr = &(patchCoupling_[pI].map());
3604                 sMeshPtr = this;
3605             }
3606             else
3607             if (procCouple && !localCouple)
3608             {
3609                 cMapPtr = &(recvMeshes_[pI].map());
3610                 sMeshPtr = &(recvMeshes_[pI].subMesh());
3611             }
3612
3613             // Alias for convenience
3614             const coupleMap& cMap = *cMapPtr;
3615             const dynamicTopoFvMesh& sMesh = *sMeshPtr;
3616
3617             // Add the new points to the coupling map
3618             const List<objectMap>& apList = slaveMap.addedPointList();
3619
3620             // Fetch the slave point
3621             label slavePoint = -1;
3622
3623             if ((slaveMap.index() == eIndex) && localCouple)
3624             {
3625                 // Cyclic edge at axis
3626                 slavePoint = newPointIndex;
3627             }
3628             else
3629             {
3630                 slavePoint = apList[0].index();
3631             }
3632
3633             // Map points
3634             if (bisectingSlave)
3635             {
3636                 // Update reverse pointMap
3637                 cMap.mapMaster
3638                 (
3639                     coupleMap::POINT,
3640                     newPointIndex,
3641                     slavePoint
3642                 );
3643
3644                 // Update pointMap
3645                 cMap.mapSlave
3646                 (
3647                     coupleMap::POINT,
3648                     slavePoint,
3649                     newPointIndex
3650                 );
3651             }
3652             else
3653             {
3654                 // Update pointMap
3655                 cMap.mapSlave
3656                 (
3657                     coupleMap::POINT,
3658                     newPointIndex,
3659                     slavePoint
3660                 );
3661
3662                 // Update reverse pointMap
3663                 cMap.mapMaster
3664                 (
3665                     coupleMap::POINT,
3666                     slavePoint,
3667                     newPointIndex
3668                 );
3669             }
3670
3671             if (debug > 2)
3672             {
3673                 Pout<< " Adding point: " << slavePoint
3674                     << " for point: " << newPointIndex
3675                     << endl;
3676             }
3677
3678             // Obtain point maps
3679             const Map<label>& pointMap = cMap.entityMap(coupleMap::POINT);
3680
3681             // Update face mapping
3682             const label faceEnum = coupleMap::FACE;
3683
3684             // Obtain references
3685             Map<label>& faceMap = cMap.entityMap(faceEnum);
3686             Map<label>& rFaceMap = cMap.reverseEntityMap(faceEnum);
3687
3688             forAll(checkFaces, faceI)
3689             {
3690                 label mfIndex = checkFaces[faceI];
3691
3692                 const face& mF = faces_[mfIndex];
3693
3694                 label mfPatch = whichPatch(mfIndex);
3695
3696                 // Check for processor match
3697                 label neiProc = getNeighbourProcessor(mfPatch);
3698
3699                 if (procCouple && !localCouple)
3700                 {
3701                     if (neiProc != procIndices_[pI])
3702                     {
3703                         continue;
3704                     }
3705                 }
3706
3707                 triFace cF
3708                 (
3709                     pointMap.found(mF[0]) ? pointMap[mF[0]] : -1,
3710                     pointMap.found(mF[1]) ? pointMap[mF[1]] : -1,
3711                     pointMap.found(mF[2]) ? pointMap[mF[2]] : -1
3712                 );
3713
3714                 // Skip mapping if all points were not found
3715                 if (cF[0] == -1 || cF[1] == -1 || cF[2] == -1)
3716                 {
3717                     continue;
3718                 }
3719
3720                 // Fetch edges connected to the slave point
3721                 const labelList& spEdges = sMesh.pointEdges_[slavePoint];
3722
3723                 forAll(spEdges, edgeI)
3724                 {
3725                     label seIndex = spEdges[edgeI];
3726
3727                     if (sMesh.whichEdgePatch(seIndex) == -1)
3728                     {
3729                         continue;
3730                     }
3731
3732                     const labelList& seFaces = sMesh.edgeFaces_[seIndex];
3733
3734                     forAll(seFaces, faceJ)
3735                     {
3736                         label sfIndex = seFaces[faceJ];
3737
3738                         if (sMesh.whichPatch(sfIndex) == -1)
3739                         {
3740                             continue;
3741                         }
3742
3743                         const face& sF = sMesh.faces_[sfIndex];
3744
3745                         if
3746                         (
3747                             triFace::compare
3748                             (
3749                                 triFace(sF[0], sF[1], sF[2]), cF
3750                             )
3751                         )
3752                         {
3753                             if (debug > 2)
3754                             {
3755                                 word pN(boundaryMesh()[mfPatch].name());
3756
3757                                 Pout<< " Found face: " << sfIndex
3758                                     << " :: " << sF
3759                                     << " with mfIndex: " << mfIndex
3760                                     << " :: " << mF
3761                                     << " Patch: " << pN
3762                                     << endl;
3763                             }
3764
3765                             if (rFaceMap.found(sfIndex))
3766                             {
3767                                 rFaceMap[sfIndex] = mfIndex;
3768                             }
3769                             else
3770                             {
3771                                 rFaceMap.insert(sfIndex, mfIndex);
3772                             }
3773
3774                             if (faceMap.found(mfIndex))
3775                             {
3776                                 faceMap[mfIndex] = sfIndex;
3777                             }
3778                             else
3779                             {
3780                                 faceMap.insert(mfIndex, sfIndex);
3781                             }
3782
3783                             matchedFaces[faceI] = true;
3784
3785                             break;
3786                         }
3787                     }
3788
3789                     if (matchedFaces[faceI])
3790                     {
3791                         break;
3792                     }
3793                 }
3794
3795                 if ((debug > 4) && !matchedFaces[faceI])
3796                 {
3797                     sMesh.writeVTK
3798                     (
3799                         "failedFacePoints_"
3800                       + Foam::name(mfIndex),
3801                         labelList(cF), 0, false, true
3802                     );
3803
3804                     writeVTK
3805                     (
3806                         "checkFaces_" + Foam::name(eIndex),
3807                         checkFaces, 2, false, true
3808                     );
3809
3810                     Pout<< " Failed to match face: "
3811                         << mfIndex << " :: " << mF
3812                         << " masterPatch: " << mfPatch
3813                         << " using comparison face: " << cF
3814                         << " on proc: " << procIndices_[pI]
3815                         << endl;
3816                 }
3817             }
3818
3819             // Update edge mapping
3820             const label edgeEnum = coupleMap::EDGE;
3821
3822             // Obtain references
3823             Map<label>& edgeMap = cMap.entityMap(edgeEnum);
3824             Map<label>& rEdgeMap = cMap.reverseEntityMap(edgeEnum);
3825
3826             forAll(checkEdges, edgeI)
3827             {
3828                 label meIndex = checkEdges[edgeI];
3829
3830                 const edge& mE = edges_[meIndex];
3831
3832                 label mePatch = whichEdgePatch(meIndex);
3833                 label neiProc = getNeighbourProcessor(mePatch);
3834
3835                 edge cE
3836                 (
3837                     pointMap.found(mE[0]) ? pointMap[mE[0]] : -1,
3838                     pointMap.found(mE[1]) ? pointMap[mE[1]] : -1
3839                 );
3840
3841                 // Skip mapping if all points were not found
3842                 if (cE[0] == -1 || cE[1] == -1)
3843                 {
3844                     continue;
3845                 }
3846
3847                 // Fetch edges connected to the slave point
3848                 const labelList& spEdges = sMesh.pointEdges_[cE[0]];
3849
3850                 forAll(spEdges, edgeJ)
3851                 {
3852                     label seIndex = spEdges[edgeJ];
3853
3854                     const edge& sE = sMesh.edges_[seIndex];
3855
3856                     if (sE == cE)
3857                     {
3858                         if (debug > 2)
3859                         {
3860                             Pout<< " Found edge: " << seIndex
3861                                 << " :: " << sE
3862                                 << " with meIndex: " << meIndex
3863                                 << " :: " << mE
3864                                 << " on proc: " << procIndices_[pI]
3865                                 << endl;
3866                         }
3867
3868                         // Update reverse map
3869                         if (rEdgeMap.found(seIndex))
3870                         {
3871                             rEdgeMap[seIndex] = meIndex;
3872                         }
3873                         else
3874                         {
3875                             rEdgeMap.insert(seIndex, meIndex);
3876                         }
3877
3878                         // Update map
3879                         if (edgeMap.found(meIndex))
3880                         {
3881                             edgeMap[meIndex] = seIndex;
3882                         }
3883                         else
3884                         {
3885                             edgeMap.insert(meIndex, seIndex);
3886                         }
3887
3888                         matchedEdges[edgeI] = true;
3889
3890                         break;
3891                     }
3892                 }
3893
3894                 if (!matchedEdges[edgeI])
3895                 {
3896                     if (procCouple && !localCouple)
3897                     {
3898                         // Rare occassion where both points
3899                         // of the edge lie on processor, but
3900                         // not the edge itself.
3901                         if (neiProc != procIndices_[pI])
3902                         {
3903                             if (debug > 2)
3904                             {
3905                                 Pout<< " Edge: " << meIndex
3906                                     << " :: " << mE
3907                                     << " with comparison: " << cE
3908                                     << " has points on processor: " << neiProc
3909                                     << " but no edge. Marking as matched."
3910                                     << endl;
3911                             }
3912
3913                             matchedEdges[edgeI] = true;
3914                         }
3915                     }
3916                 }
3917
3918                 if ((debug > 4) && !matchedEdges[edgeI])
3919                 {
3920                     sMesh.writeVTK
3921                     (
3922                         "failedEdge_"
3923                       + Foam::name(meIndex),
3924                         labelList(cE), 0, false, true
3925                     );
3926
3927                     writeVTK
3928                     (
3929                         "checkEdges_" + Foam::name(eIndex),
3930                         checkEdges, 1, false, true
3931                     );
3932
3933                     Pout<< " Failed to match edge: "
3934                         << meIndex << " :: " << mE
3935                         << " using comparison edge: " << cE
3936                         << " on proc: " << procIndices_[pI]
3937                         << endl;
3938                 }
3939             }
3940
3941             // Push operation into coupleMap
3942             if (procCouple && !localCouple)
3943             {
3944                 cMap.pushOperation
3945                 (
3946                     slaveMap.index(),
3947                     coupleMap::BISECTION
3948                 );
3949             }
3950         }
3951
3952         // Ensure that all entities were matched
3953         label nFailFace = 0, nFailEdge = 0;
3954
3955         forAll(matchedFaces, faceI)
3956         {
3957             if (!matchedFaces[faceI])
3958             {
3959                 ++nFailFace;
3960             }
3961         }
3962
3963         forAll(matchedEdges, edgeI)
3964         {
3965             if (!matchedEdges[edgeI])
3966             {
3967                 ++nFailEdge;
3968             }
3969         }
3970
3971         if (nFailFace || nFailEdge)
3972         {
3973             Pout<< " Failed to match all entities. " << nl
3974                 << "  Faces: " << nFailFace << nl
3975                 << "  Edges: " << nFailEdge << nl
3976                 << abort(FatalError);
3977         }
3978     }
3979
3980     if (debug > 2)
3981     {
3982         label bPatch = whichEdgePatch(eIndex);
3983
3984         const polyBoundaryMesh& boundary = boundaryMesh();
3985
3986         if (bPatch == -1)
3987         {
3988             Pout<< "Patch: Internal" << nl;
3989         }
3990         else
3991         if (bPatch < boundary.size())
3992         {
3993             Pout<< "Patch: " << boundary[bPatch].name() << nl;
3994         }
3995         else
3996         {
3997             Pout<< " New patch: " << bPatch << endl;
3998         }
3999
4000         Pout<< "vertexHull: " << vertexHull << nl
4001             << "Edges: " << edgeHull << nl
4002             << "Faces: " << faceHull << nl
4003             << "Cells: " << cellHull << nl;
4004
4005         Pout<< "Modified cells: " << nl;
4006
4007         forAll(cellHull, cellI)
4008         {
4009             if (cellHull[cellI] == -1)
4010             {
4011                 continue;
4012             }
4013
4014             Pout<< cellHull[cellI] << ":: "
4015                 << cells_[cellHull[cellI]]
4016                 << nl;
4017         }
4018
4019         Pout<< nl << "Added cells: " << nl;
4020
4021         forAll(addedCellIndices, cellI)
4022         {
4023             if (addedCellIndices[cellI] == -1)
4024             {
4025                 continue;
4026             }
4027
4028             Pout<< addedCellIndices[cellI] << ":: "
4029                 << cells_[addedCellIndices[cellI]] << nl
4030                 << "lengthScale: " << lengthScale_[addedCellIndices[cellI]]
4031                 << nl;
4032         }
4033
4034         Pout<< nl << "Modified faces: " << nl;
4035
4036         forAll(faceHull, faceI)
4037         {
4038             Pout<< faceHull[faceI] << ":: "
4039                 << faces_[faceHull[faceI]] << ": "
4040                 << owner_[faceHull[faceI]] << ": "
4041                 << neighbour_[faceHull[faceI]] << " "
4042                 << "faceEdges:: " << faceEdges_[faceHull[faceI]]
4043                 << nl;
4044         }
4045
4046         Pout<< nl << "Added faces: " << nl;
4047
4048         forAll(addedFaceIndices, faceI)
4049         {
4050             Pout<< addedFaceIndices[faceI] << ":: "
4051                 << faces_[addedFaceIndices[faceI]] << ": "
4052                 << owner_[addedFaceIndices[faceI]] << ": "
4053                 << neighbour_[addedFaceIndices[faceI]] << " "
4054                 << "faceEdges:: " << faceEdges_[addedFaceIndices[faceI]]
4055                 << nl;
4056         }
4057
4058         forAll(addedIntFaceIndices, faceI)
4059         {
4060             if (addedIntFaceIndices[faceI] == -1)
4061             {
4062                 continue;
4063             }
4064
4065             Pout<< addedIntFaceIndices[faceI] << ":: "
4066                 << faces_[addedIntFaceIndices[faceI]] << ": "
4067                 << owner_[addedIntFaceIndices[faceI]] << ": "
4068                 << neighbour_[addedIntFaceIndices[faceI]] << " "
4069                 << "faceEdges:: " << faceEdges_[addedIntFaceIndices[faceI]]
4070                 << nl;
4071         }
4072
4073         Pout<< nl << "New edge:: " << newEdgeIndex
4074             << ": " << edges_[newEdgeIndex] << nl
4075             << " edgeFaces:: " << edgeFaces_[newEdgeIndex]
4076             << nl;
4077
4078         Pout<< nl << "Added edges: " << nl;
4079
4080         forAll(addedEdgeIndices, edgeI)
4081         {
4082             Pout<< addedEdgeIndices[edgeI]
4083                 << ":: " << edges_[addedEdgeIndices[edgeI]] << nl
4084                 << " edgeFaces:: " << edgeFaces_[addedEdgeIndices[edgeI]]
4085                 << nl;
4086         }
4087
4088         Pout<< "New Point:: " << newPointIndex << nl
4089             << "pointEdges:: " << pointEdges_[newPointIndex] << nl;
4090
4091         // Flush buffer
4092         Pout<< endl;
4093
4094         // Write out VTK files after change
4095         if (debug > 3)
4096         {
4097             labelList newHull(cellHull.size() + addedCellIndices.size(), 0);
4098
4099             // Combine both lists into one.
4100             forAll(cellHull, i)
4101             {
4102                 newHull[i] = cellHull[i];
4103             }
4104
4105             label start = cellHull.size();
4106
4107             for(label i = start; i < newHull.size(); i++)
4108             {
4109                 newHull[i] = addedCellIndices[i - start];
4110             }
4111
4112             // Write out VTK files after change
4113             //  - Using old-points for convenience in post-processing
4114             writeVTK
4115             (
4116                 Foam::name(eIndex)
4117               + '(' + Foam::name(origEdge[0])
4118               + ',' + Foam::name(origEdge[1]) + ')'
4119               + "_Bisect_1",
4120                 newHull,
4121                 3, false, true
4122             );
4123         }
4124     }
4125
4126     // Set the flag
4127     topoChangeFlag_ = true;
4128
4129     // Increment the counter
4130     status(TOTAL_BISECTIONS)++;
4131
4132     // Increment the number of modifications
4133     status(TOTAL_MODIFICATIONS)++;
4134
4135     // Specify that the operation was successful
4136     map.type() = 1;
4137
4138     // Return the changeMap
4139     return map;
4140 }
4141
4142
4143 // Utility method to compute the quality of a
4144 // vertex hull around an edge after bisection.
4145 scalar dynamicTopoFvMesh::computeBisectionQuality
4146 (
4147     const label eIndex
4148 ) const
4149 {
4150     scalar minQuality = GREAT, minVolume = GREAT;
4151     scalar cQuality = 0.0, oldVolume = 0.0;
4152
4153     // Obtain a reference to this edge
4154     const edge& edgeToCheck = edges_[eIndex];
4155
4156     // Obtain point references
4157     const point& a = points_[edgeToCheck[0]];
4158     const point& c = points_[edgeToCheck[1]];
4159
4160     const point& aOld = oldPoints_[edgeToCheck[0]];
4161     const point& cOld = oldPoints_[edgeToCheck[1]];
4162
4163     // Compute the mid-point of the edge
4164     point midPoint = 0.5 * (a + c);
4165     point oldPoint = 0.5 * (aOld + cOld);
4166
4167     dynamicLabelList eCells(10);
4168
4169     const labelList& eFaces = edgeFaces_[eIndex];
4170
4171     // Accumulate cells connected to this edge
4172     forAll(eFaces, faceI)
4173     {
4174         label own = owner_[eFaces[faceI]];
4175         label nei = neighbour_[eFaces[faceI]];
4176
4177         if (findIndex(eCells, own) == -1)
4178         {
4179             eCells.append(own);
4180         }
4181
4182         if (nei == -1)
4183         {
4184             continue;
4185         }
4186
4187         if (findIndex(eCells, nei) == -1)
4188         {
4189             eCells.append(nei);
4190         }
4191     }
4192
4193     // Loop through all cells and compute quality
4194     forAll(eCells, cellI)
4195     {
4196         label cellIndex = eCells[cellI];
4197
4198         const cell& checkCell = cells_[cellIndex];
4199
4200         // Find two faces that don't contain the edge
4201         forAll(checkCell, faceI)
4202         {
4203             const face& checkFace = faces_[checkCell[faceI]];
4204
4205             bool check0 = (findIndex(checkFace, edgeToCheck[0]) == -1);
4206             bool check1 = (findIndex(checkFace, edgeToCheck[1]) == -1);
4207
4208             if ((check0 && !check1) || (!check0 && check1))
4209             {
4210                 // Check orientation
4211                 if (owner_[checkCell[faceI]] == cellIndex)
4212                 {
4213                     cQuality =
4214                     (
4215                         tetMetric_
4216                         (
4217                             points_[checkFace[2]],
4218                             points_[checkFace[1]],
4219                             points_[checkFace[0]],
4220                             midPoint
4221                         )
4222                     );
4223
4224                     oldVolume =
4225                     (
4226                         tetPointRef
4227                         (
4228                             oldPoints_[checkFace[2]],
4229                             oldPoints_[checkFace[1]],
4230                             oldPoints_[checkFace[0]],
4231                             oldPoint
4232                         ).mag()
4233                     );
4234                 }
4235                 else
4236                 {
4237                     cQuality =
4238                     (
4239                         tetMetric_
4240                         (
4241                             points_[checkFace[0]],
4242                             points_[checkFace[1]],
4243                             points_[checkFace[2]],
4244                             midPoint
4245                         )
4246                     );
4247
4248                     oldVolume =
4249                     (
4250                         tetPointRef
4251                         (
4252                             oldPoints_[checkFace[0]],
4253                             oldPoints_[checkFace[1]],
4254                             oldPoints_[checkFace[2]],
4255                             oldPoint
4256                         ).mag()
4257                     );
4258                 }
4259
4260                 // Check if the quality is worse
4261                 minQuality = Foam::min(cQuality, minQuality);
4262                 minVolume = Foam::min(oldVolume, minVolume);
4263             }
4264         }
4265     }
4266
4267     // Ensure that the mesh is valid
4268     if (minQuality < sliverThreshold_)
4269     {
4270         if (debug > 3 && minQuality < 0.0)
4271         {
4272             writeVTK(Foam::name(eIndex) + "_iCells", eCells);
4273         }
4274
4275         if (debug > 2)
4276         {
4277             InfoIn
4278             (
4279                 "scalar dynamicTopoFvMesh::computeBisectionQuality"
4280                 "(const label eIndex) const"
4281             )
4282                 << " Bisecting edge will fall below the"
4283                 << " sliver threshold of: " << sliverThreshold_ << nl
4284                 << " Edge: " << eIndex << ": " << edgeToCheck << nl
4285                 << " Minimum Quality: " << minQuality << nl
4286                 << " Minimum Volume: " << minVolume << nl
4287                 << " Mid point: " << midPoint << nl
4288                 << " Old point: " << oldPoint << nl
4289                 << endl;
4290         }
4291     }
4292
4293     // If a negative old-volume was encountered,
4294     // return an invalid quality.
4295     if (minVolume < 0.0)
4296     {
4297         return minVolume;
4298     }
4299
4300     return minQuality;
4301 }
4302
4303
4304 // Slice the mesh at a particular location
4305 void dynamicTopoFvMesh::sliceMesh
4306 (
4307     const labelPair& pointPair
4308 )
4309 {
4310     if (debug > 1)
4311     {
4312         Pout<< nl << nl
4313             << "Pair: " << pointPair
4314             << " is to be used for mesh slicing. " << endl;
4315     }
4316
4317     label patchIndex = -1;
4318     scalar dx = 0.0;
4319     vector gCentre = vector::zero;
4320     FixedList<vector, 2> fC(vector::zero);
4321
4322     if (is2D())
4323     {
4324         patchIndex = whichPatch(pointPair.first());
4325
4326         // Fetch face centres
4327         fC[0] = faces_[pointPair.first()].centre(points_);
4328         fC[1] = faces_[pointPair.second()].centre(points_);
4329     }
4330     else
4331     {
4332         // Find the patch that the edge-vertex is connected to.
4333         const labelList& pEdges = pointEdges_[pointPair.first()];
4334
4335         forAll(pEdges, edgeI)
4336         {
4337             if ((patchIndex = whichEdgePatch(pEdges[edgeI])) > -1)
4338             {
4339                 break;
4340             }
4341         }
4342
4343         fC[0] = points_[pointPair.first()];
4344         fC[1] = points_[pointPair.second()];
4345     }
4346
4347     linePointRef lpr(fC[0], fC[1]);
4348
4349     // Specify the centre.
4350     gCentre = lpr.centre();
4351
4352     // Specify a search distance
4353     dx = lpr.mag();
4354
4355     // Is this edge in the vicinity of a previous slice-point?
4356     if (lengthEstimator().checkOldSlices(gCentre))
4357     {
4358         if (debug > 1)
4359         {
4360             Pout<< nl << nl
4361                 << "Pair: " << pointPair
4362                 << " is too close to another slice point. "
4363                 << endl;
4364         }
4365
4366         // Too close to another slice-point. Bail out.
4367         return;
4368     }
4369
4370     // Choose a box around the centre and scan all
4371     // surface entities that fall into this region.
4372     boundBox bBox
4373     (
4374         gCentre - vector(dx, dx, dx),
4375         gCentre + vector(dx, dx, dx)
4376     );
4377
4378     vector p = vector::zero, N = vector::zero;
4379     Map<vector> checkPoints, surfFaces;
4380     Map<edge> checkEdges;
4381
4382     if (is2D())
4383     {
4384         // Assign plane point / normal
4385         p = gCentre;
4386
4387         vector gNorm = faces_[pointPair.first()].normal(points_);
4388
4389         gNorm /= (mag(gNorm) + VSMALL);
4390
4391         // Since this is 2D, assume XY-plane here.
4392         N = (gNorm ^ vector(0.0, 0.0, 1.0));
4393     }
4394     else
4395     {
4396         // Prepare surface points / edges for Dijkstra's algorithm
4397         for (label edgeI = nOldInternalEdges_; edgeI < nEdges_; edgeI++)
4398         {
4399             if (edgeFaces_[edgeI].empty())
4400             {
4401                 continue;
4402             }
4403
4404             if (whichEdgePatch(edgeI) == patchIndex)
4405             {
4406                 const edge& surfaceEdge = edges_[edgeI];
4407
4408                 if
4409                 (
4410                     (bBox.contains(points_[surfaceEdge[0]])) &&
4411                     (bBox.contains(points_[surfaceEdge[1]]))
4412                 )
4413                 {
4414                     checkEdges.insert(edgeI, surfaceEdge);
4415
4416                     if (!checkPoints.found(surfaceEdge[0]))
4417                     {
4418                         checkPoints.insert
4419                         (
4420                             surfaceEdge[0],
4421                             points_[surfaceEdge[0]]
4422                         );
4423                     }
4424
4425                     if (!checkPoints.found(surfaceEdge[1]))
4426                     {
4427                         checkPoints.insert
4428                         (
4429                             surfaceEdge[1],
4430                             points_[surfaceEdge[1]]
4431                         );
4432                     }
4433
4434                     // Add surface faces as well.
4435                     const labelList& eFaces = edgeFaces_[edgeI];
4436
4437                     forAll(eFaces, faceI)
4438                     {
4439                         if
4440                         (
4441                             (neighbour_[eFaces[faceI]] == -1) &&
4442                             (!surfFaces.found(eFaces[faceI]))
4443                         )
4444                         {
4445                             vector surfNorm =
4446                             (
4447                                 faces_[eFaces[faceI]].normal(points_)
4448                             );
4449
4450                             surfFaces.insert(eFaces[faceI], surfNorm);
4451                         }
4452                     }
4453                 }
4454             }
4455         }
4456
4457         if (debug > 1)
4458         {
4459             Pout<< nl << nl
4460                 << " Point [0]: " << points_[pointPair.first()] << nl
4461                 << " Point [1]: " << points_[pointPair.second()] << endl;
4462
4463             if (debug > 3)
4464             {
4465                 writeVTK("slicePoints", checkPoints.toc(), 0);
4466                 writeVTK("sliceEdges", checkEdges.toc(), 1);
4467             }
4468         }
4469
4470         // Find the shortest path using Dijkstra's algorithm.
4471         Map<label> shortestPath;
4472
4473         bool foundPath =
4474         (
4475             meshOps::Dijkstra
4476             (
4477                 checkPoints,
4478                 checkEdges,
4479                 pointPair.first(),
4480                 pointPair.second(),
4481                 shortestPath
4482             )
4483         );
4484
4485         // First normalize all face-normals
4486         forAllIter(Map<vector>, surfFaces, sIter)
4487         {
4488             sIter() /= (mag(sIter()) + VSMALL);
4489         }
4490
4491         if (foundPath)
4492         {
4493             // Next, take cross-products with every other
4494             // vector in the list, and accumulate.
4495             forAllIter(Map<vector>, surfFaces, sIterI)
4496             {
4497                 forAllIter(Map<vector>, surfFaces, sIterJ)
4498                 {
4499                     if (sIterI.key() != sIterJ.key())
4500                     {
4501                         vector n = (sIterI() ^ sIterJ());
4502
4503                         n /= (mag(n) + VSMALL);
4504
4505                         // Reverse the vector if necessary
4506                         if ((N & n) < 0.0)
4507                         {
4508                             n = -n;
4509                         }
4510
4511                         N += n;
4512                     }
4513                 }
4514             }
4515
4516             N /= (mag(N) + VSMALL);
4517
4518             // Obtain point position
4519             p = gCentre;
4520         }
4521         else
4522         {
4523             // Probably a membrane-type configuration.
4524             labelHashSet checkCells;
4525
4526             // Prepare a bounding cylinder with radius dx.
4527             forAllIter(Map<vector>, surfFaces, sIter)
4528             {
4529                 const face& thisFace = faces_[sIter.key()];
4530
4531                 if (thisFace.which(pointPair.first()) > -1)
4532                 {
4533                     N += sIter();
4534                 }
4535             }
4536
4537             // Normalize and reverse.
4538             N /= -(mag(N) + VSMALL);
4539
4540             vector a0 = points_[pointPair.first()];
4541             vector a1 = points_[pointPair.second()];
4542             scalar dist = mag(a1 - a0);
4543
4544             forAll(cells_, cellI)
4545             {
4546                 if (cells_[cellI].empty())
4547                 {
4548                     continue;
4549                 }
4550
4551                 scalar v = 0.0;
4552                 vector x = vector::zero;
4553
4554                 meshOps::cellCentreAndVolume
4555                 (
4556                     cellI,
4557                     points_,
4558                     faces_,
4559                     cells_,
4560                     owner_,
4561                     x,
4562                     v
4563                 );
4564
4565                 vector rx = (x - a0);
4566                 vector ra = (rx & N)*N;
4567
4568                 // Check if point falls off cylinder ends.
4569                 if (mag(ra) > dist || mag(ra) < 0.0)
4570                 {
4571                     continue;
4572                 }
4573
4574                 vector r = (rx - ra);
4575
4576                 // Check if the magnitude of 'r' is within radius.
4577                 if (mag(r) < dx)
4578                 {
4579                     checkCells.insert(cellI);
4580                 }
4581             }
4582
4583             labelList cList = checkCells.toc();
4584
4585             if (debug > 1)
4586             {
4587                 Pout<< "Dijkstra's algorithm could not find a path." << endl;
4588
4589                 if (debug > 3)
4590                 {
4591                     writeVTK("checkCells", cList, 3);
4592                 }
4593             }
4594
4595             changeMap sliceMap =
4596             (
4597                 removeCells
4598                 (
4599                     cList,
4600                     patchIndex,
4601                     "Slice_"
4602                   + Foam::name(pointPair.first()) + '_'
4603                   + Foam::name(pointPair.second())
4604                 )
4605             );
4606
4607             if (debug)
4608             {
4609                 checkConnectivity(10);
4610             }
4611
4612             // Accumulate a list of projection points
4613             labelHashSet pPoints;
4614
4615             const List<objectMap>& afList = sliceMap.addedFaceList();
4616
4617             forAll(afList, faceI)
4618             {
4619                 const face& thisFace = faces_[afList[faceI].index()];
4620
4621                 forAll(thisFace, pointI)
4622                 {
4623                     if (!pPoints.found(thisFace[pointI]))
4624                     {
4625                         pPoints.insert(thisFace[pointI]);
4626                     }
4627                 }
4628             }
4629
4630             // Now project points in a series of intermediate steps.
4631
4632             // Add an entry to sliceBoxes.
4633             lengthEstimator().appendBox(bBox);
4634
4635             return;
4636         }
4637     }
4638
4639     if (debug > 1)
4640     {
4641         Pout<< nl << nl
4642             << " Plane point: " << p << nl
4643             << " Plane normal: " << N << endl;
4644     }
4645
4646     // Mark cells and interior faces that fall
4647     // within the bounding box.
4648     labelHashSet checkCells, checkFaces, splitFaces;
4649     Map<bool> cellColors;
4650
4651     forAll(faces_, faceI)
4652     {
4653         if (faces_[faceI].empty())
4654         {
4655             continue;
4656         }
4657
4658         if (is2D() && faces_[faceI].size() == 3)
4659         {
4660             continue;
4661         }
4662
4663         vector fCentre = faces_[faceI].centre(points_);
4664
4665         FixedList<label, 2> cellsToCheck(-1);
4666         cellsToCheck[0] = owner_[faceI];
4667         cellsToCheck[1] = neighbour_[faceI];
4668
4669         if (bBox.contains(fCentre) && cellsToCheck[1] != -1)
4670         {
4671             // Add this internal face to the list.
4672             checkFaces.insert(faceI);
4673
4674             scalar volume = 0.0;
4675             vector centre = vector::zero;
4676
4677             forAll(cellsToCheck, cellI)
4678             {
4679                 if (!checkCells.found(cellsToCheck[cellI]))
4680                 {
4681                     meshOps::cellCentreAndVolume
4682                     (
4683                         cellsToCheck[cellI],
4684                         points_,
4685                         faces_,
4686                         cells_,
4687                         owner_,
4688                         centre,
4689                         volume
4690                     );
4691
4692                     checkCells.insert(cellsToCheck[cellI]);
4693
4694                     if (((centre - p) & N) > 0.0)
4695                     {
4696                         cellColors.insert(cellsToCheck[cellI], true);
4697                     }
4698                     else
4699                     {
4700                         cellColors.insert(cellsToCheck[cellI], false);
4701                     }
4702                 }
4703             }
4704         }
4705     }
4706
4707     // Prepare a list of internal faces for mesh splitting.
4708     forAllIter(labelHashSet, checkFaces, fIter)
4709     {
4710         if
4711         (
4712             cellColors[owner_[fIter.key()]]
4713          != cellColors[neighbour_[fIter.key()]]
4714         )
4715         {
4716             splitFaces.insert(fIter.key());
4717         }
4718
4719         // Loop through all points (and associated pointEdges)
4720         // for this face, and check if connected cells are also
4721         // present in the checkCells/cellColors list
4722         if (is2D())
4723         {
4724             const labelList& fEdges = faceEdges_[fIter.key()];
4725
4726             forAll(fEdges, edgeI)
4727             {
4728                 const labelList& eFaces = edgeFaces_[fEdges[edgeI]];
4729
4730                 forAll(eFaces, faceI)
4731                 {
4732                     label own = owner_[eFaces[faceI]];
4733                     label nei = neighbour_[eFaces[faceI]];
4734
4735                     if (!checkCells.found(own))
4736                     {
4737                         scalar volume = 0.0;
4738                         vector centre = vector::zero;
4739
4740                         meshOps::cellCentreAndVolume
4741                         (
4742                             own,
4743                             points_,
4744                             faces_,
4745                             cells_,
4746                             owner_,
4747                             centre,
4748                             volume
4749                         );
4750
4751                         checkCells.insert(own);
4752
4753                         if (((centre - p) & N) > 0.0)
4754                         {
4755                             cellColors.insert(own, true);
4756                         }
4757                         else
4758                         {
4759                             cellColors.insert(own, false);
4760                         }
4761                     }
4762
4763                     if (!checkCells.found(nei) && nei != -1)
4764                     {
4765                         scalar volume = 0.0;
4766                         vector centre = vector::zero;
4767
4768                         meshOps::cellCentreAndVolume
4769                         (
4770                             nei,
4771                             points_,
4772                             faces_,
4773                             cells_,
4774                             owner_,
4775                             centre,
4776                             volume
4777                         );
4778
4779                         checkCells.insert(nei);
4780
4781                         if (((centre - p) & N) > 0.0)
4782                         {
4783                             cellColors.insert(nei, true);
4784                         }
4785                         else
4786                         {
4787                             cellColors.insert(nei, false);
4788                         }
4789                     }
4790                 }
4791             }
4792         }
4793         else
4794         {
4795             const face& faceToCheck = faces_[fIter.key()];
4796
4797             forAll(faceToCheck, pointI)
4798             {
4799                 const labelList& pEdges = pointEdges_[faceToCheck[pointI]];
4800
4801                 forAll(pEdges, edgeI)
4802                 {
4803                     const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
4804
4805                     forAll(eFaces, faceI)
4806                     {
4807                         label own = owner_[eFaces[faceI]];
4808                         label nei = neighbour_[eFaces[faceI]];
4809
4810                         if (!checkCells.found(own))
4811                         {
4812                             scalar volume = 0.0;
4813                             vector centre = vector::zero;
4814
4815                             meshOps::cellCentreAndVolume
4816                             (
4817                                 own,
4818                                 points_,
4819                                 faces_,
4820                                 cells_,
4821                                 owner_,
4822                                 centre,
4823                                 volume
4824                             );
4825
4826                             checkCells.insert(own);
4827
4828                             if (((centre - p) & N) > 0.0)
4829                             {
4830                                 cellColors.insert(own, true);
4831                             }
4832                             else
4833                             {
4834                                 cellColors.insert(own, false);
4835                             }
4836                         }
4837
4838                         if (!checkCells.found(nei) && nei != -1)
4839                         {
4840                             scalar volume = 0.0;
4841                             vector centre = vector::zero;
4842
4843                             meshOps::cellCentreAndVolume
4844                             (
4845                                 nei,
4846                                 points_,
4847                                 faces_,
4848                                 cells_,
4849                                 owner_,
4850                                 centre,
4851                                 volume
4852                             );
4853
4854                             checkCells.insert(nei);
4855
4856                             if (((centre - p) & N) > 0.0)
4857                             {
4858                                 cellColors.insert(nei, true);
4859                             }
4860                             else
4861                             {
4862                                 cellColors.insert(nei, false);
4863                             }
4864                         }
4865                     }
4866                 }
4867             }
4868         }
4869     }
4870
4871     if (debug > 3)
4872     {
4873         writeVTK("splitFaces", splitFaces.toc(), 2);
4874         writeVTK("checkCells", checkCells.toc(), 3);
4875     }
4876
4877     // Pass this info into the splitInternalFaces routine.
4878     splitInternalFaces
4879     (
4880         patchIndex,
4881         splitFaces.toc(),
4882         cellColors
4883     );
4884
4885     // Add an entry to sliceBoxes.
4886     lengthEstimator().appendBox(bBox);
4887 }
4888
4889
4890 // Add cell layer above specified patch
4891 const changeMap dynamicTopoFvMesh::addCellLayer
4892 (
4893     const label patchID
4894 )
4895 {
4896     changeMap map;
4897
4898     // Maps for added entities
4899     Map<label> addedPoints;
4900     Map<label> addedHEdges, addedVEdges, currentVEdges;
4901     Map<label> addedHFaces, addedVFaces, currentVFaces;
4902     Map<labelPair> addedCells;
4903
4904     // Allocate a list of patch faces
4905     dynamicLabelList patchFaces(patchSizes_[patchID]);
4906
4907     // Loop through all patch faces and create a cell for each
4908     for (label faceI = nOldInternalFaces_; faceI < faces_.size(); faceI++)
4909     {
4910         label pIndex = whichPatch(faceI);
4911
4912         if (pIndex != patchID)
4913         {
4914             continue;
4915         }
4916
4917         // Add face to the list
4918         patchFaces.append(faceI);
4919
4920         // Add a new cell
4921         label cIndex = owner_[faceI];
4922         scalar newLengthScale = -1.0;
4923         const cell& ownCell = cells_[cIndex];
4924
4925         if (edgeRefinement_)
4926         {
4927             newLengthScale = lengthScale_[cIndex];
4928         }
4929
4930         label newCellIndex =
4931         (
4932             insertCell
4933             (
4934                 cell(ownCell.size()),
4935                 newLengthScale
4936             )
4937         );
4938
4939         // Update maps
4940         map.addCell(newCellIndex, labelList(1, cIndex));
4941         addedCells.insert(cIndex, labelPair(newCellIndex, 0));
4942     }
4943
4944     labelList mP(2, -1);
4945
4946     forAll(patchFaces, indexI)
4947     {
4948         label faceI = patchFaces[indexI];
4949         label cIndex = owner_[faceI];
4950
4951         // Fetch appropriate face / cell
4952         //  - Make copies, since holding references
4953         //    to data within this loop is unsafe.
4954         const face bFace = faces_[faceI];
4955         const cell ownCell = cells_[cIndex];
4956
4957         // Configure a new face for insertion
4958         face newHFace(bFace);
4959         labelList newHFaceEdges(bFace.size(), -1);
4960
4961         // Get the opposing face from the cell
4962         oppositeFace oFace = ownCell.opposingFace(faceI, faces_);
4963
4964         if (!oFace.found())
4965         {
4966             // Something's wrong here.
4967             FatalErrorIn
4968             (
4969                 "const changeMap dynamicTopoFvMesh::addCellLayer"
4970                 "(const label patchID)"
4971             )
4972                 << " Face: " << faceI << " :: " << bFace << nl
4973                 << " has no opposing face in cell: "
4974                 << cIndex << " :: " << ownCell << nl
4975                 << abort(FatalError);
4976         }
4977
4978         // Create points
4979         forAll(bFace, pointI)
4980         {
4981             label pIndex = bFace[pointI];
4982
4983             // Skip if we've added this already
4984             if (addedPoints.found(pIndex))
4985             {
4986                 continue;
4987             }
4988
4989             // Set master points
4990             mP[0] = pIndex;
4991             mP[1] = oFace[pointI];
4992
4993             label newPointIndex =
4994             (
4995                 insertPoint
4996                 (
4997                     0.5 * (points_[mP[0]] + points_[mP[1]]),
4998                     oldPoints_[mP[0]],
4999                     mP
5000                 )
5001             );
5002
5003             // Update maps
5004             map.addPoint(newPointIndex, mP);
5005             addedPoints.insert(pIndex, newPointIndex);
5006         }
5007
5008         // Fetch faceEdges from opposite faces.
5009         //  - Make copies, since holding references is unsafe
5010         const labelList bfEdges = faceEdges_[faceI];
5011         const labelList ofEdges = faceEdges_[oFace.oppositeIndex()];
5012
5013         // Create edges for each edge of the new horizontal face
5014         forAll(bfEdges, edgeI)
5015         {
5016             label beIndex = bfEdges[edgeI];
5017
5018             // Skip if we've added this already
5019             if (addedHEdges.found(beIndex))
5020             {
5021                 // Update face edges for the new horizontal face
5022                 newHFaceEdges[edgeI] = addedHEdges[beIndex];
5023
5024                 continue;
5025             }
5026
5027             // Configure the new edge
5028             label oeIndex = -1;
5029             const edge bEdge = edges_[beIndex];
5030
5031             // Build an edge for comparison
5032             edge cEdge
5033             (
5034                 oFace[bFace.which(bEdge[0])],
5035                 oFace[bFace.which(bEdge[1])]
5036             );
5037
5038             forAll(ofEdges, edgeJ)
5039             {
5040                 const edge& ofEdge = edges_[ofEdges[edgeJ]];
5041
5042                 if (cEdge == ofEdge)
5043                 {
5044                     oeIndex = ofEdges[edgeJ];
5045                     break;
5046                 }
5047             }
5048
5049             if (oeIndex < 0)
5050             {
5051                 FatalErrorIn
5052                 (
5053                     "const changeMap dynamicTopoFvMesh::addCellLayer"
5054                     "(const label patchID)"
5055                 )
5056                     << " Could not find comparison edge: " << cEdge
5057                     << " for edge: " << bEdge
5058                     << abort(FatalError);
5059             }
5060
5061             // Fetch patch information
5062             label hEdgePatch = whichEdgePatch(oeIndex);
5063
5064             // Set indices
5065             edge newHEdge
5066             (
5067                 addedPoints[bEdge[0]],
5068                 addedPoints[bEdge[1]]
5069             );
5070
5071             // Insert a new edge with empty edgeFaces
5072             label newHEdgeIndex =
5073             (
5074                 insertEdge
5075                 (
5076                     hEdgePatch,
5077                     newHEdge,
5078                     labelList(0)
5079                 )
5080             );
5081
5082             // Update maps
5083             map.addEdge(newHEdgeIndex);
5084             addedHEdges.insert(beIndex, newHEdgeIndex);
5085
5086             // Update face edges for the new horizontal face
5087             newHFaceEdges[edgeI] = newHEdgeIndex;
5088
5089             // Add a new vertical face for this edge
5090             label vFaceIndex = -1;
5091
5092             // Find a vertical face that contains both edges
5093             const labelList& beFaces = edgeFaces_[beIndex];
5094
5095             forAll(beFaces, faceJ)
5096             {
5097                 const labelList& testEdges = faceEdges_[beFaces[faceJ]];
5098
5099                 if
5100                 (
5101                     (findIndex(testEdges, beIndex) > -1) &&
5102                     (findIndex(testEdges, oeIndex) > -1)
5103                 )
5104                 {
5105                     vFaceIndex = beFaces[faceJ];
5106                     break;
5107                 }
5108             }
5109
5110             if (vFaceIndex < 0)
5111             {
5112                 FatalErrorIn
5113                 (
5114                     "const changeMap dynamicTopoFvMesh::addCellLayer"
5115                     "(const label patchID)"
5116                 )
5117                     << " Could not find an appropriate vertical face"
5118                     << " containing edges: " << oeIndex
5119                     << " and " << beIndex
5120                     << abort(FatalError);
5121             }
5122
5123             // Find two vertical edges on this face
5124             const labelList& vfEdges = faceEdges_[vFaceIndex];
5125
5126             forAll(vfEdges, edgeJ)
5127             {
5128                 const edge& vfEdge = edges_[vfEdges[edgeJ]];
5129
5130                 forAll(bEdge, i)
5131                 {
5132                     if (vfEdge == edge(bEdge[i], cEdge[i]))
5133                     {
5134                         // Skip if we've added this already
5135                         if (addedVEdges.found(bEdge[i]))
5136                         {
5137                             continue;
5138                         }
5139
5140                         label veIndex = vfEdges[edgeJ];
5141
5142                         // Fetch edge patch information
5143                         label vEdgePatch = whichEdgePatch(veIndex);
5144
5145                         // Set indices
5146                         edge newVEdge
5147                         (
5148                             bEdge[i],
5149                             addedPoints[bEdge[i]]
5150                         );
5151
5152                         // Insert a new edge with empty edgeFaces
5153                         label newVEdgeIndex =
5154                         (
5155                             insertEdge
5156                             (
5157                                 vEdgePatch,
5158                                 newVEdge,
5159                                 labelList(0)
5160                             )
5161                         );
5162
5163                         // Update maps
5164                         map.addEdge(newVEdgeIndex);
5165                         addedVEdges.insert(bEdge[i], newVEdgeIndex);
5166
5167                         // Note edge indices for later renumbering
5168                         currentVEdges.insert(bEdge[i], veIndex);
5169                     }
5170                 }
5171             }
5172
5173             // Configure the new vertical face
5174             face newVFace(faces_[vFaceIndex]);
5175             label newOwner = -1, newNeighbour = -1;
5176
5177             label oldOwner = owner_[vFaceIndex];
5178             label oldNeighbour = neighbour_[vFaceIndex];
5179
5180             // Fetch owner / neighbour
5181             newOwner = addedCells[oldOwner].first();
5182
5183             if (oldNeighbour > -1)
5184             {
5185                 newNeighbour = addedCells[oldNeighbour].first();
5186             }
5187
5188             // Replace point indices on the new face
5189             forAll(bEdge, i)
5190             {
5191                 meshOps::replaceLabel
5192                 (
5193                     cEdge[i],
5194                     addedPoints[bEdge[i]],
5195                     newVFace
5196                 );
5197             }
5198
5199             // Note face indices for later renumbering
5200             currentVFaces.insert(beIndex, vFaceIndex);
5201
5202             // Check if reversal is necessary
5203             if ((newNeighbour < newOwner) && (newNeighbour > -1))
5204             {
5205                 // Flip face
5206                 newVFace = newVFace.reverseFace();
5207
5208                 // Swap addressing
5209                 Foam::Swap(newOwner, newNeighbour);
5210                 Foam::Swap(oldOwner, oldNeighbour);
5211             }
5212
5213             // Configure faceEdges for the new vertical face
5214             labelList newVFaceEdges(4, -1);
5215
5216             newVFaceEdges[0] = beIndex;
5217             newVFaceEdges[1] = newHEdgeIndex;
5218             newVFaceEdges[2] = addedVEdges[bEdge[0]];
5219             newVFaceEdges[3] = addedVEdges[bEdge[1]];
5220
5221             // Add the new vertical face
5222             label newVFaceIndex =
5223             (
5224                 insertFace
5225                 (
5226                     whichPatch(vFaceIndex),
5227                     newVFace,
5228                     newOwner,
5229                     newNeighbour,
5230                     newVFaceEdges
5231                 )
5232             );
5233
5234             // Update maps
5235             map.addFace(newVFaceIndex, labelList(1, vFaceIndex));
5236             addedVFaces.insert(beIndex, newVFaceIndex);
5237
5238             // Update face count on the new cells
5239             cells_[newOwner][addedCells[oldOwner].second()++] =
5240             (
5241                 newVFaceIndex
5242             );
5243
5244             if (newNeighbour > -1)
5245             {
5246                 cells_[newNeighbour][addedCells[oldNeighbour].second()++] =
5247                 (
5248                     newVFaceIndex
5249                 );
5250             }
5251
5252             // Size up edgeFaces for each edge
5253             forAll(newVFaceEdges, edgeJ)
5254             {
5255                 label vfeIndex = newVFaceEdges[edgeJ];
5256
5257                 meshOps::sizeUpList
5258                 (
5259                     newVFaceIndex,
5260                     edgeFaces_[vfeIndex]
5261                 );
5262             }
5263         }
5264
5265         // Add a new interior face, with identical orientation
5266         forAll(newHFace, pointI)
5267         {
5268             newHFace[pointI] = addedPoints[newHFace[pointI]];
5269         }
5270
5271         // Add the new horizontal face
5272         label newHFaceIndex =
5273         (
5274             insertFace
5275             (
5276                 -1,
5277                 newHFace,
5278                 cIndex,
5279                 addedCells[cIndex].first(),
5280                 newHFaceEdges
5281             )
5282         );
5283
5284         // Update maps
5285         map.addFace(newHFaceIndex, labelList(1, faceI));
5286         addedHFaces.insert(faceI, newHFaceIndex);
5287
5288         // Replace index on the old cell
5289         meshOps::replaceLabel
5290         (
5291             faceI,
5292             newHFaceIndex,
5293             cells_[cIndex]
5294         );
5295
5296         // Update face count on the new cell
5297         label newCellIndex = addedCells[cIndex].first();
5298
5299         // Modify owner for the existing boundary face
5300         owner_[faceI] = newCellIndex;
5301
5302         cells_[newCellIndex][addedCells[cIndex].second()++] = faceI;
5303         cells_[newCellIndex][addedCells[cIndex].second()++] = newHFaceIndex;
5304
5305         // Size up edgeFaces for each edge
5306         forAll(newHFaceEdges, edgeI)
5307         {
5308             label hfeIndex = newHFaceEdges[edgeI];
5309
5310             meshOps::sizeUpList
5311             (
5312                 newHFaceIndex,
5313                 edgeFaces_[hfeIndex]
5314             );
5315         }
5316     }
5317
5318     // Renumber vertical edges
5319     forAllConstIter(Map<label>, currentVEdges, eIter)
5320     {
5321         // Fetch reference to edge
5322         edge& curEdge = edges_[eIter()];
5323
5324         if (curEdge[0] == eIter.key())
5325         {
5326             curEdge[0] = addedPoints[eIter.key()];
5327         }
5328
5329         if (curEdge[1] == eIter.key())
5330         {
5331             curEdge[1] = addedPoints[eIter.key()];
5332         }
5333
5334         // Size down pointEdges
5335         if (is3D())
5336         {
5337             meshOps::sizeDownList
5338             (
5339                 eIter(),
5340                 pointEdges_[eIter.key()]
5341             );
5342
5343             meshOps::sizeUpList
5344             (
5345                 eIter(),
5346                 pointEdges_[addedPoints[eIter.key()]]
5347             );
5348         }
5349     }
5350
5351     // Renumber vertical faces
5352     forAllConstIter(Map<label>, currentVFaces, fIter)
5353     {
5354         // Fetch reference to existing edge
5355         const edge& bEdge = edges_[fIter.key()];
5356
5357         // Replace point indices on vertical face
5358         forAll(bEdge, i)
5359         {
5360             meshOps::replaceLabel
5361             (
5362                 bEdge[i],
5363                 addedPoints[bEdge[i]],
5364                 faces_[fIter()]
5365             );
5366         }
5367
5368         // Replace edge on the existing vertical face
5369         meshOps::replaceLabel
5370         (
5371             fIter.key(),
5372             addedHEdges[fIter.key()],
5373             faceEdges_[fIter()]
5374         );
5375
5376         // Remove old face on existing boundary edge
5377         meshOps::sizeDownList
5378         (
5379             fIter(),
5380             edgeFaces_[fIter.key()]
5381         );
5382
5383         // Add old face to new horizontal edge
5384         meshOps::sizeUpList
5385         (
5386             fIter(),
5387             edgeFaces_[addedHEdges[fIter.key()]]
5388         );
5389     }
5390
5391     // Now that all old / new cells possess correct connectivity,
5392     // compute mapping information.
5393     const List<objectMap>& afList = map.addedFaceList();
5394
5395     forAll(afList, indexI)
5396     {
5397         label parent = afList[indexI].masterObjects()[0];
5398
5399         if (whichPatch(afList[indexI].index()) == -1)
5400         {
5401             // Interior faces get default mapping
5402             if (whichPatch(parent) == -1)
5403             {
5404                 setFaceMapping(parent);
5405             }
5406
5407             setFaceMapping(afList[indexI].index());
5408         }
5409     }
5410
5411     // Set the flag
5412     topoChangeFlag_ = true;
5413
5414     // Specify that the operation was successful
5415     map.type() = 1;
5416
5417     // Return the changeMap
5418     return map;
5419 }
5420
5421
5422 // Split a set of internal faces into boundary faces
5423 //   - Add boundary faces and edges to the patch specified by 'patchIndex'
5424 //   - Cell color should specify a binary value dictating either side
5425 //     of the split face.
5426 void dynamicTopoFvMesh::splitInternalFaces
5427 (
5428     const label patchIndex,
5429     const labelList& internalFaces,
5430     const Map<bool>& cellColors
5431 )
5432 {
5433     Map<label> mirrorPointLabels;
5434     FixedList<Map<label>, 2> mirrorEdgeLabels, mirrorFaceLabels;
5435
5436     // First loop through the list and accumulate a list of
5437     // points and edges that need to be duplicated.
5438     forAll(internalFaces, faceI)
5439     {
5440         const face& faceToCheck = faces_[internalFaces[faceI]];
5441
5442         forAll(faceToCheck, pointI)
5443         {
5444             if (!mirrorPointLabels.found(faceToCheck[pointI]))
5445             {
5446                 mirrorPointLabels.insert(faceToCheck[pointI], -1);
5447             }
5448         }
5449
5450         const labelList& fEdges = faceEdges_[internalFaces[faceI]];
5451
5452         forAll(fEdges, edgeI)
5453         {
5454             if (!mirrorEdgeLabels[0].found(fEdges[edgeI]))
5455             {
5456                 mirrorEdgeLabels[0].insert(fEdges[edgeI], -1);
5457             }
5458         }
5459     }
5460
5461     // Now for every point in the list, add a new one.
5462     // Add a mapping entry as well.
5463     forAllIter(Map<label>, mirrorPointLabels, pIter)
5464     {
5465         // Obtain a copy of the point before adding it,
5466         // since the reference might become invalid during list resizing.
5467         point newPoint = points_[pIter.key()];
5468         point oldPoint = oldPoints_[pIter.key()];
5469
5470         pIter() = insertPoint(newPoint, oldPoint, labelList(1, pIter.key()));
5471
5472         if (is3D())
5473         {
5474             const labelList& pEdges = pointEdges_[pIter.key()];
5475
5476             labelHashSet edgesToRemove;
5477
5478             forAll(pEdges, edgeI)
5479             {
5480                 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
5481
5482                 bool allTrue = true;
5483
5484                 forAll(eFaces, faceI)
5485                 {
5486                     label own = owner_[eFaces[faceI]];
5487                     label nei = neighbour_[eFaces[faceI]];
5488
5489                     // Check if an owner/neighbour cell is false
5490                     if (!cellColors[own])
5491                     {
5492                         allTrue = false;
5493                         break;
5494                     }
5495
5496                     if (nei != -1)
5497                     {
5498                         if (!cellColors[nei])
5499                         {
5500                             allTrue = false;
5501                             break;
5502                         }
5503                     }
5504                 }
5505
5506                 if (allTrue)
5507                 {
5508                     // Mark this edge label to be discarded later
5509                     edgesToRemove.insert(pEdges[edgeI]);
5510                 }
5511             }
5512
5513             // It is dangerous to use the pointEdges references,
5514             // so call it using array-lookup instead.
5515             forAllIter(labelHashSet, edgesToRemove, hsIter)
5516             {
5517                 // Add the edge to the mirror point list
5518                 meshOps::sizeUpList
5519                 (
5520                     hsIter.key(),
5521                     pointEdges_[pIter()]
5522                 );
5523
5524                 // Remove the edge from the original point list
5525                 meshOps::sizeDownList
5526                 (
5527                     hsIter.key(),
5528                     pointEdges_[pIter.key()]
5529                 );
5530             }
5531         }
5532     }
5533
5534     if (debug > 3)
5535     {
5536         label i = 0;
5537         labelList mPoints(mirrorPointLabels.size());
5538
5539         if (is3D())
5540         {
5541             forAllIter(Map<label>, mirrorPointLabels, pIter)
5542             {
5543                 writeVTK
5544                 (
5545                     "pEdges_o_" + Foam::name(pIter.key()) + '_',
5546                     pointEdges_[pIter.key()],
5547                     1
5548                 );
5549
5550                 writeVTK
5551                 (
5552                     "pEdges_m_" + Foam::name(pIter()) + '_',
5553                     pointEdges_[pIter()],
5554                     1
5555                 );
5556
5557                 mPoints[i++] = pIter();
5558             }
5559
5560             writeVTK
5561             (
5562                 "points_o_",
5563                 mirrorPointLabels.toc(),
5564                 0
5565             );
5566
5567             writeVTK
5568             (
5569                 "points_m_",
5570                 mPoints,
5571                 0
5572             );
5573         }
5574     }
5575
5576     // For every internal face, add a new one.
5577     //  - Stick to the rule:
5578     //    [1] Every cell marked false keeps the existing entities.
5579     //    [2] Every cell marked true gets new points/edges/faces.
5580     //  - If faces are improperly oriented, reverse them.
5581     forAll(internalFaces, faceI)
5582     {
5583         FixedList<face, 2> newFace;
5584         FixedList<label, 2> newFaceIndex(-1);
5585         FixedList<label, 2> newOwner(-1);
5586
5587         label oldOwn = owner_[internalFaces[faceI]];
5588         label oldNei = neighbour_[internalFaces[faceI]];
5589
5590         if (cellColors[oldOwn] && !cellColors[oldNei])
5591         {
5592             // The owner gets a new boundary face.
5593             // Note that orientation is already correct.
5594             newFace[0] = faces_[internalFaces[faceI]];
5595
5596             // The neighbour needs to have its face reversed
5597             // and moved to the boundary patch, thereby getting
5598             // deleted in the process.
5599             newFace[1] = newFace[0].reverseFace();
5600
5601             newOwner[0] = oldOwn;
5602             newOwner[1] = oldNei;
5603         }
5604         else
5605         if (!cellColors[oldOwn] && cellColors[oldNei])
5606         {
5607             // The neighbour gets a new boundary face.
5608             // The face is oriented in the opposite sense, however.
5609             newFace[0] = faces_[internalFaces[faceI]].reverseFace();
5610
5611             // The owner keeps the existing face and orientation.
5612             // But it also needs to be moved to the boundary.
5613             newFace[1] = faces_[internalFaces[faceI]];
5614
5615             newOwner[0] = oldNei;
5616             newOwner[1] = oldOwn;
5617         }
5618         else
5619         {
5620             // Something's wrong here.
5621             FatalErrorIn
5622             (
5623                 "dynamicTopoFvMesh::splitInternalFaces\n"
5624                 "(\n"
5625                 "    const label patchIndex,\n"
5626                 "    const labelList& internalFaces,\n"
5627                 "    const Map<bool>& cellColors\n"
5628                 ")\n"
5629             )
5630                 << nl << " Face: "
5631                 << internalFaces[faceI]
5632                 << " has cells which are improperly marked: " << nl
5633                 << oldOwn << ":: " << cellColors[oldOwn] << nl
5634                 << oldNei << ":: " << cellColors[oldNei]
5635                 << abort(FatalError);
5636         }
5637
5638         // Renumber point labels for the first new face.
5639         forAll(newFace[0], pointI)
5640         {
5641             newFace[0][pointI] = mirrorPointLabels[newFace[0][pointI]];
5642         }
5643
5644         // Insert the new boundary faces.
5645         forAll(newFace, indexI)
5646         {
5647             // Make an identical faceEdges entry.
5648             // This will be renumbered once new edges are added.
5649             labelList newFaceEdges(faceEdges_[internalFaces[faceI]]);
5650
5651             newFaceIndex[indexI] =
5652             (
5653                 insertFace
5654                 (
5655                     patchIndex,
5656                     newFace[indexI],
5657                     newOwner[indexI],
5658                     -1,
5659                     newFaceEdges
5660                 )
5661             );
5662
5663             // Replace face labels on cells
5664             meshOps::replaceLabel
5665             (
5666                 internalFaces[faceI],
5667                 newFaceIndex[indexI],
5668                 cells_[newOwner[indexI]]
5669             );
5670
5671             // Make face mapping entries for posterity.
5672             mirrorFaceLabels[indexI].insert
5673             (
5674                 internalFaces[faceI],
5675                 newFaceIndex[indexI]
5676             );
5677         }
5678     }
5679
5680     // For every edge in the list, add a new one.
5681     // We'll deal with correcting edgeFaces later.
5682     forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
5683     {
5684         // Obtain copies for the append method
5685         edge origEdge = edges_[eIter.key()];
5686         labelList newEdgeFaces(edgeFaces_[eIter.key()]);
5687
5688         eIter() =
5689         (
5690             insertEdge
5691             (
5692                 patchIndex,
5693                 edge
5694                 (
5695                     mirrorPointLabels[origEdge[0]],
5696                     mirrorPointLabels[origEdge[1]]
5697                 ),
5698                 newEdgeFaces
5699             )
5700         );
5701
5702         // Is the original edge an internal one?
5703         // If it is, we need to move it to the boundary.
5704         if (whichEdgePatch(eIter.key()) == -1)
5705         {
5706             label rplEdgeIndex =
5707             (
5708                 insertEdge
5709                 (
5710                     patchIndex,
5711                     origEdge,
5712                     newEdgeFaces
5713                 )
5714             );
5715
5716             // Map the new entry.
5717             mirrorEdgeLabels[1].insert(eIter.key(), rplEdgeIndex);
5718         }
5719         else
5720         {
5721             // This is already a boundary edge.
5722             // Make an identical map.
5723             mirrorEdgeLabels[1].insert(eIter.key(), eIter.key());
5724         }
5725     }
5726
5727     // Correct edgeFaces for all new edges
5728     forAll(mirrorEdgeLabels, indexI)
5729     {
5730         forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
5731         {
5732             labelList& eFaces = edgeFaces_[eIter()];
5733
5734             labelHashSet facesToRemove;
5735
5736             forAll(eFaces, faceI)
5737             {
5738                 bool remove = false;
5739
5740                 label own = owner_[eFaces[faceI]];
5741                 label nei = neighbour_[eFaces[faceI]];
5742
5743                 if
5744                 (
5745                     (!cellColors[own] && !indexI) ||
5746                     ( cellColors[own] &&  indexI)
5747                 )
5748                 {
5749                     remove = true;
5750                 }
5751
5752                 if (nei != -1)
5753                 {
5754                     if
5755                     (
5756                         (!cellColors[nei] && !indexI) ||
5757                         ( cellColors[nei] &&  indexI)
5758                     )
5759                     {
5760                         remove = true;
5761                     }
5762                 }
5763
5764                 if (mirrorFaceLabels[indexI].found(eFaces[faceI]))
5765                 {
5766                     // Perform a replacement instead of a removal.
5767                     eFaces[faceI] = mirrorFaceLabels[indexI][eFaces[faceI]];
5768
5769                     remove = false;
5770                 }
5771
5772                 if (remove)
5773                 {
5774                     facesToRemove.insert(eFaces[faceI]);
5775                 }
5776             }
5777
5778             // Now successively size down edgeFaces.
5779             // It is dangerous to use the eFaces reference,
5780             // so call it using array-lookup.
5781             forAllIter(labelHashSet, facesToRemove, hsIter)
5782             {
5783                 meshOps::sizeDownList(hsIter.key(), edgeFaces_[eIter()]);
5784             }
5785         }
5786     }
5787
5788     // Renumber faceEdges for all faces connected to new edges
5789     forAll(mirrorEdgeLabels, indexI)
5790     {
5791         forAllIter(Map<label>, mirrorEdgeLabels[indexI], eIter)
5792         {
5793             const labelList& eFaces = edgeFaces_[eIter()];
5794
5795             forAll(eFaces, faceI)
5796             {
5797                 labelList& fEdges = faceEdges_[eFaces[faceI]];
5798
5799                 forAll(fEdges, edgeI)
5800                 {
5801                     if (mirrorEdgeLabels[indexI].found(fEdges[edgeI]))
5802                     {
5803                         fEdges[edgeI] =
5804                         (
5805                             mirrorEdgeLabels[indexI][fEdges[edgeI]]
5806                         );
5807                     }
5808                 }
5809             }
5810         }
5811     }
5812
5813     if (is2D())
5814     {
5815         // Renumber edges and faces
5816         forAllIter(Map<label>, mirrorEdgeLabels[0], eIter)
5817         {
5818             const labelList& eFaces = edgeFaces_[eIter()];
5819
5820             // Two levels of indirection to ensure
5821             // that all entities we renumbered.
5822             // A flip-side for the lack of a pointEdges list in 2D.
5823             forAll(eFaces, faceI)
5824             {
5825                 const labelList& fEdges = faceEdges_[eFaces[faceI]];
5826
5827                 forAll(fEdges, edgeI)
5828                 {
5829                     // Renumber this edge.
5830                     edge& edgeToCheck = edges_[fEdges[edgeI]];
5831
5832                     forAll(edgeToCheck, pointI)
5833                     {
5834                         if (mirrorPointLabels.found(edgeToCheck[pointI]))
5835                         {
5836                             edgeToCheck[pointI] =
5837                             (
5838                                 mirrorPointLabels[edgeToCheck[pointI]]
5839                             );
5840                         }
5841                     }
5842
5843                     // Also renumber faces connected to this edge.
5844                     const labelList& efFaces = edgeFaces_[fEdges[edgeI]];
5845
5846                     forAll(efFaces, faceJ)
5847                     {
5848                         face& faceToCheck = faces_[efFaces[faceJ]];
5849
5850                         forAll(faceToCheck, pointI)
5851                         {
5852                             if
5853                             (
5854                                 mirrorPointLabels.found(faceToCheck[pointI])
5855                             )
5856                             {
5857                                 faceToCheck[pointI] =
5858                                 (
5859                                     mirrorPointLabels[faceToCheck[pointI]]
5860                                 );
5861                             }
5862                         }
5863                     }
5864                 }
5865             }
5866         }
5867     }
5868     else
5869     {
5870         // Point renumbering of entities connected to mirror points
5871         forAllIter(Map<label>, mirrorPointLabels, pIter)
5872         {
5873             const labelList& pEdges = pointEdges_[pIter()];
5874
5875             forAll(pEdges, edgeI)
5876             {
5877                 // Renumber this edge.
5878                 edge& edgeToCheck = edges_[pEdges[edgeI]];
5879
5880                 forAll(edgeToCheck, pointI)
5881                 {
5882                     if (edgeToCheck[pointI] == pIter.key())
5883                     {
5884                         edgeToCheck[pointI] = pIter();
5885                     }
5886                 }
5887
5888                 // Also renumber faces connected to this edge.
5889                 const labelList& eFaces = edgeFaces_[pEdges[edgeI]];
5890
5891                 forAll(eFaces, faceI)
5892                 {
5893                     face& faceToCheck = faces_[eFaces[faceI]];
5894
5895                     forAll(faceToCheck, pointI)
5896                     {
5897                         if (faceToCheck[pointI] == pIter.key())
5898                         {
5899                             faceToCheck[pointI] = pIter();
5900                         }
5901                     }
5902                 }
5903             }
5904         }
5905     }
5906
5907     // Now that we're done with the internal faces, remove them.
5908     forAll(internalFaces, faceI)
5909     {
5910         removeFace(internalFaces[faceI]);
5911     }
5912
5913     // Remove old internal edges as well.
5914     forAllIter(Map<label>, mirrorEdgeLabels[1], eIter)
5915     {
5916         if (eIter.key() != eIter())
5917         {
5918             removeEdge(eIter.key());
5919         }
5920     }
5921
5922     // Set the flag
5923     topoChangeFlag_ = true;
5924 }
5925
5926 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
5927
5928 } // End namespace Foam
5929
5930 // ************************************************************************* //