src/dynamicMesh/slidingInterface/coupleSlidingInterface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "slidingInterface.H"
  28 #include "polyTopoChange.H"
  29 #include "polyMesh.H"
  30 #include "primitiveMesh.H"
  31 #include "enrichedPatch.H"
  32 #include "DynamicList.H"
  33 #include "pointHit.H"
  34 #include "triPointRef.H"
  35 #include "plane.H"
  36 #include "polyTopoChanger.H"
  37 #include "polyAddPoint.H"
  38 #include "polyRemovePoint.H"
  39 #include "polyAddFace.H"
  40 #include "polyModifyPoint.H"
  41 #include "polyModifyFace.H"
  42 #include "polyRemoveFace.H"
  43
  44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  45
  46 const Foam::scalar Foam::slidingInterface::edgeCoPlanarTol_ = 0.8;
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 // Index of debug signs:
  51 // Index of debug signs:
  52 // . - loop of the edge-to-face interaction detection
  53 // x - reversal of direction in edge-to-face interaction detection
  54 // + - complete edge-to-face interaction detection
  55 // z - incomplete edge-to-face interaction detection.  This may be OK for edges
  56 //     crossing from one to the other side of multiply connected master patch
  57
  58 // e - adding a point on edge
  59 // f - adding a point on face
  60 // . - collecting edges off another face for edge-to-edge cut
  61 // + - finished collection of edges
  62 // * - cut both master and slave
  63 // n - cutting new edge
  64 // t - intersection exists but it is outside of tolerance
  65 // x - missed slave edge in cut
  66 // - - missed master edge in cut
  67 // u - edge already used in cutting
  68
  69 void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
  70 {
  71     if (debug)
  72     {
  73         Pout<< "void slidingInterface::coupleInterface"
  74             << "(polyTopoChange& ref) : "
  75             << "Coupling sliding interface " << name() << endl;
  76     }
  77
  78     const polyMesh& mesh = topoChanger().mesh();
  79
  80     const pointField& points = mesh.points();
  81     const faceList& faces = mesh.faces();
  82
  83     const labelList& own = mesh.faceOwner();
  84     const labelList& nei = mesh.faceNeighbour();
  85     const faceZoneMesh& faceZones = mesh.faceZones();
  86
  87     const primitiveFacePatch& masterPatch =
  88         faceZones[masterFaceZoneID_.index()]();
  89
  90     const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];
  91
  92     const boolList& masterPatchFlip =
  93         faceZones[masterFaceZoneID_.index()].flipMap();
  94
  95     const primitiveFacePatch& slavePatch =
  96         faceZones[slaveFaceZoneID_.index()]();
  97
  98     const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];
  99
 100     const boolList& slavePatchFlip =
 101         faceZones[slaveFaceZoneID_.index()].flipMap();
 102
 103     const edgeList& masterEdges = masterPatch.edges();
 104     const labelListList& masterPointEdges = masterPatch.pointEdges();
 105     const labelList& masterMeshPoints = masterPatch.meshPoints();
 106     const pointField& masterLocalPoints = masterPatch.localPoints();
 107     const labelListList& masterFaceFaces = masterPatch.faceFaces();
 108     const labelListList& masterFaceEdges = masterPatch.faceEdges();
 109     const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();
 110
 111     const edgeList& slaveEdges = slavePatch.edges();
 112     const labelListList& slavePointEdges = slavePatch.pointEdges();
 113     const labelList& slaveMeshPoints = slavePatch.meshPoints();
 114     const pointField& slaveLocalPoints = slavePatch.localPoints();
 115     const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
 116     const vectorField& slavePointNormals = slavePatch.pointNormals();
 117
 118     // Collect projection addressing
 119     if
 120     (
 121         !(
 122             slavePointPointHitsPtr_
 123          && slavePointEdgeHitsPtr_
 124          && slavePointFaceHitsPtr_
 125          && masterPointEdgeHitsPtr_
 126         )
 127     )
 128     {
 129         FatalErrorIn
 130         (
 131             "void slidingInterface::coupleInterface("
 132             "polyTopoChange& ref) const"
 133         )   << "Point projection addressing not available."
 134             << abort(FatalError);
 135     }
 136
 137     const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
 138     const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
 139     const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
 140     const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
 141     const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
 142
 143     // Create enriched patch
 144     enrichedPatch cutPatch
 145     (
 146         masterPatch,
 147         slavePatch,
 148         slavePointPointHits,
 149         slavePointEdgeHits,
 150         slavePointFaceHits
 151     );
 152
 153     // Get reference to list of added point for the enriched patch.
 154     // This will be used for point addition
 155     Map<point>& pointMap = cutPatch.pointMap();
 156
 157     // Get reference to the list of merged points
 158     Map<label>& pointMergeMap = cutPatch.pointMergeMap();
 159
 160     // Create mapping for every merged point of the slave patch
 161     forAll (slavePointPointHits, pointI)
 162     {
 163         if (slavePointPointHits[pointI] >= 0)
 164         {
 165 // Pout << "Inserting point merge pair: " << slaveMeshPoints[pointI] << " : " << masterMeshPoints[slavePointPointHits[pointI]] << endl;
 166             pointMergeMap.insert
 167             (
 168                 slaveMeshPoints[pointI],
 169                 masterMeshPoints[slavePointPointHits[pointI]]
 170             );
 171         }
 172     }
 173
 174     // Collect the list of used edges for every slave edge
 175
 176     List<labelHashSet> usedMasterEdges(slaveEdges.size());
 177
 178     // Collect of slave point hits
 179     forAll (slavePointPointHits, pointI)
 180     {
 181         // For point hits, add all point-edges from master side to all point
 182         const labelList& curSlaveEdges = slavePointEdges[pointI];
 183
 184         if (slavePointPointHits[pointI] > -1)
 185         {
 186             const labelList& curMasterEdges =
 187                 masterPointEdges[slavePointPointHits[pointI]];
 188
 189             // Mark all current master edges as used for all the current slave
 190             // edges
 191             forAll (curSlaveEdges, slaveEdgeI)
 192             {
 193                 labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];
 194
 195                 forAll (curMasterEdges, masterEdgeI)
 196                 {
 197                     sm.insert(curMasterEdges[masterEdgeI]);
 198                 }
 199             }
 200         }
 201         else if (slavePointEdgeHits[pointI] > -1)
 202         {
 203             // For edge hits, add the master edge
 204             forAll (curSlaveEdges, slaveEdgeI)
 205             {
 206                 usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
 207                 (
 208                     slavePointEdgeHits[pointI]
 209                 );
 210             }
 211         }
 212     }
 213
 214     // Collect off master point hits
 215     // For every master point that has hit an edge, add all edges coming from
 216     // that point to the slave edge that has been hit
 217     forAll (masterPointEdgeHits, masterPointI)
 218     {
 219         if (masterPointEdgeHits[masterPointI] > -1)
 220         {
 221             const labelList& curMasterEdges = masterPointEdges[masterPointI];
 222
 223             labelHashSet& sm =
 224                 usedMasterEdges[masterPointEdgeHits[masterPointI]];
 225
 226             forAll (curMasterEdges, masterEdgeI)
 227             {
 228                 sm.insert(curMasterEdges[masterEdgeI]);
 229             }
 230         }
 231     }
 232
 233 //     Pout << "used edges: " << endl;
 234 //     forAll (usedMasterEdges, edgeI)
 235 //     {
 236 //         Pout << "edge: " << edgeI << " used: " << usedMasterEdges[edgeI].toc() << endl;
 237 //     }
 238
 239     // For every master and slave edge make a list of points to be added into
 240     // that edge.
 241     List<DynamicList<label> > pointsIntoMasterEdges(masterEdges.size());
 242     List<DynamicList<label> > pointsIntoSlaveEdges(slaveEdges.size());
 243
 244     // Add all points from the slave patch that have hit the edge
 245     forAll (slavePointEdgeHits, pointI)
 246     {
 247         if (slavePointEdgeHits[pointI] > -1)
 248         {
 249             // Create a new point on the master edge
 250
 251             point edgeCutPoint =
 252                 masterEdges[slavePointEdgeHits[pointI]].line
 253                 (
 254                     masterLocalPoints
 255                 ).nearestDist(projectedSlavePoints[pointI]).hitPoint();
 256
 257             label newPoint =
 258                 ref.setAction
 259                 (
 260                     polyAddPoint
 261                     (
 262                         edgeCutPoint,             // point
 263                         slaveMeshPoints[pointI],  // master point
 264                         cutPointZoneID_.index(),  // zone for point
 265                         true                      // supports a cell
 266                     )
 267                 );
 268 //             Pout << "Inserting merge pair off edge: " << slaveMeshPoints[pointI] << " " << newPoint << " cut point: " << edgeCutPoint << " orig: " << slaveLocalPoints[pointI] << " proj: " << projectedSlavePoints[pointI] << endl;
 269             // Add the new edge point into the merge map
 270             pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
 271
 272             pointsIntoMasterEdges[slavePointEdgeHits[pointI]].append
 273             (
 274                 newPoint
 275             );
 276
 277             // Add the point into the enriched patch map
 278             pointMap.insert
 279             (
 280                 newPoint,
 281                 edgeCutPoint
 282             );
 283
 284             if (debug)
 285             {
 286                 Pout<< "e";
 287 //                 Pout<< newPoint << " = " << edgeCutPoint << endl;
 288             }
 289         }
 290     }
 291
 292     if (debug)
 293     {
 294         Pout<< endl;
 295     }
 296
 297     // Add all points from the slave patch that have hit a face
 298     forAll (slavePointFaceHits, pointI)
 299     {
 300         if
 301         (
 302             slavePointPointHits[pointI] < 0
 303          && slavePointEdgeHits[pointI] < 0
 304          && slavePointFaceHits[pointI].hit()
 305         )
 306         {
 307             label newPoint =
 308                 ref.setAction
 309                 (
 310                     polyAddPoint
 311                     (
 312                         projectedSlavePoints[pointI],   // point
 313                         slaveMeshPoints[pointI],        // master point
 314                         cutPointZoneID_.index(),        // zone for point
 315                         true                            // supports a cell
 316                     )
 317                 );
 318 //             Pout << "Inserting merge pair off face: " << slaveMeshPoints[pointI] << " " << newPoint << endl;
 319             // Add the new edge point into the merge map
 320             pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);
 321
 322             // Add the point into the enriched patch map
 323             pointMap.insert
 324             (
 325                 newPoint,
 326                 projectedSlavePoints[pointI]
 327             );
 328
 329             if (debug)
 330             {
 331                 Pout<< "f: " << newPoint << " = "
 332                     << projectedSlavePoints[pointI] << endl;
 333             }
 334         }
 335     }
 336
 337     forAll (masterPointEdgeHits, pointI)
 338     {
 339         if (masterPointEdgeHits[pointI] > -1)
 340         {
 341             pointsIntoSlaveEdges[masterPointEdgeHits[pointI]].append
 342             (
 343                 masterMeshPoints[pointI]
 344             );
 345         }
 346     }
 347
 348     // Cut all slave edges.
 349     // Collect all master edges the slave edge interacts with.  Skip
 350     // all the edges already marked as used.  For every unused edge,
 351     // calculate the cut and insert the new point into the master and
 352     // slave edge.
 353     // For the edge selection algorithm, see, comment in
 354     // slidingInterfaceProjectPoints.C.
 355     // Edge cutting algoritm:
 356     // As the master patch defines the cutting surface, the newly
 357     // inserted point needs to be on the master edge.  Also, in 3-D
 358     // the pair of edges generally misses each other rather than
 359     // intersect.  Therefore, the intersection is calculated using the
 360     // plane the slave edge defines during projection.  The plane is
 361     // defined by the centrepoint of the edge in the original
 362     // configuration and the projected end points.  In case of flat
 363     // geometries (when the three points are colinear), the plane is
 364     // defined by the two projected end-points and the slave edge
 365     // normal used as the in-plane vector.  When the intersection
 366     // point is created, it is added as a new point for both the
 367     // master and the slave edge.
 368     // In order to be able to re-create the points on edges in mesh
 369     // motion without the topology change, the edge pair used to
 370     // create the cut point needs to be recorded in
 371     // cutPointEdgePairMap
 372
 373     // Note.  "Processing slave edges" code is repeated twice in the
 374     // sliding intergace coupling in order to allow the point
 375     // projection to be done separately from the actual cutting.
 376     // Please change consistently with slidingInterfaceProjectPoints.C
 377     //
 378     if (debug)
 379     {
 380         Pout << "Processing slave edges " << endl;
 381     }
 382
 383     if (!cutPointEdgePairMapPtr_)
 384     {
 385         FatalErrorIn
 386         (
 387             "void slidingInterface::coupleInterface("
 388             "polyTopoChange& ref) const"
 389         )   << "Cut point edge pair map pointer not set."
 390             << abort(FatalError);
 391     }
 392
 393     Map<Pair<edge> >& addToCpepm = *cutPointEdgePairMapPtr_;
 394
 395     // Clear the old map
 396     addToCpepm.clear();
 397
 398     // Create a map of faces the edge can interact with
 399     labelHashSet curFaceMap
 400     (
 401         nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
 402     );
 403
 404     labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);
 405
 406     forAll (slaveEdges, edgeI)
 407     {
 408         const edge& curEdge = slaveEdges[edgeI];
 409
 410         if
 411         (
 412             slavePointFaceHits[curEdge.start()].hit()
 413          || slavePointFaceHits[curEdge.end()].hit()
 414         )
 415         {
 416             labelHashSet& curUme = usedMasterEdges[edgeI];
 417 //             Pout<< "Doing edge " << edgeI << " curEdge: " << curEdge << " curUme: " << curUme << endl;
 418             // Clear the maps
 419             curFaceMap.clear();
 420             addedFaces.clear();
 421
 422             // Grab the faces for start and end points.
 423             const label startFace =
 424                 slavePointFaceHits[curEdge.start()].hitObject();
 425             const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
 426 //             Pout << "startFace: " << slavePointFaceHits[curEdge.start()] << " endFace: " << slavePointFaceHits[curEdge.end()] << endl;
 427             // Insert the start face into the list
 428             curFaceMap.insert(startFace);
 429             addedFaces.insert(startFace);
 430 //             Pout << "curFaceMap: " << curFaceMap.toc() << endl;
 431             label nSweeps = 0;
 432             bool completed = false;
 433
 434             while (nSweeps < edgeFaceEscapeLimit_)
 435             {
 436                 nSweeps++;
 437
 438                 if (addedFaces.found(endFace))
 439                 {
 440                     completed = true;
 441                 }
 442
 443                 // Add all face neighbours of face in the map
 444                 const labelList cf = addedFaces.toc();
 445                 addedFaces.clear();
 446
 447                 forAll (cf, cfI)
 448                 {
 449                     const labelList& curNbrs = masterFaceFaces[cf[cfI]];
 450
 451                     forAll (curNbrs, nbrI)
 452                     {
 453                         if (!curFaceMap.found(curNbrs[nbrI]))
 454                         {
 455                             curFaceMap.insert(curNbrs[nbrI]);
 456                             addedFaces.insert(curNbrs[nbrI]);
 457                         }
 458                     }
 459                 }
 460
 461                 if (completed) break;
 462
 463                 if (debug)
 464                 {
 465                     Pout << ".";
 466                 }
 467             }
 468
 469             if (!completed)
 470             {
 471                 if (debug)
 472                 {
 473                     Pout << "x";
 474                 }
 475
 476                 // It is impossible to reach the end from the start, probably
 477                 // due to disconnected domain.  Do search in opposite direction
 478
 479                 label nReverseSweeps = 0;
 480
 481                 addedFaces.clear();
 482                 addedFaces.insert(endFace);
 483
 484                 while (nReverseSweeps < edgeFaceEscapeLimit_)
 485                 {
 486                     nReverseSweeps++;
 487
 488                     if (addedFaces.found(startFace))
 489                     {
 490                         completed = true;
 491                     }
 492
 493                     // Add all face neighbours of face in the map
 494                     const labelList cf = addedFaces.toc();
 495                     addedFaces.clear();
 496
 497                     forAll (cf, cfI)
 498                     {
 499                         const labelList& curNbrs = masterFaceFaces[cf[cfI]];
 500
 501                         forAll (curNbrs, nbrI)
 502                         {
 503                             if (!curFaceMap.found(curNbrs[nbrI]))
 504                             {
 505                                 curFaceMap.insert(curNbrs[nbrI]);
 506                                 addedFaces.insert(curNbrs[nbrI]);
 507                             }
 508                         }
 509                     }
 510
 511                     if (completed) break;
 512
 513                     if (debug)
 514                     {
 515                         Pout << ".";
 516                     }
 517                 }
 518             }
 519
 520             if (completed)
 521             {
 522                 if (debug)
 523                 {
 524                     Pout << "+ ";
 525                 }
 526             }
 527             else
 528             {
 529                 if (debug)
 530                 {
 531                     Pout << "z ";
 532                 }
 533             }
 534
 535             // Collect the edges
 536
 537             // Create a map of edges the edge can interact with
 538             labelHashSet curMasterEdgesMap
 539             (
 540                 nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
 541             );
 542
 543             const labelList curFaces = curFaceMap.toc();
 544 //             Pout << "curFaces: " << curFaces << endl;
 545             forAll (curFaces, faceI)
 546             {
 547 //                 Pout<< "face: " << curFaces[faceI] << " "
 548 //                     << masterPatch[curFaces[faceI]]
 549 //                     << " local: "
 550 //                     << masterPatch.localFaces()[curFaces[faceI]]
 551 //                     << endl;
 552                 const labelList& me = masterFaceEdges[curFaces[faceI]];
 553
 554                 forAll (me, meI)
 555                 {
 556                     curMasterEdgesMap.insert(me[meI]);
 557                 }
 558             }
 559
 560             const labelList curMasterEdges = curMasterEdgesMap.toc();
 561
 562             // For all master edges to intersect, skip the ones
 563             // already used and cut the rest with a cutting plane.  If
 564             // the intersection point, falls inside of both edges, it
 565             // is valid.
 566
 567             // Note.
 568             // The edge cutting code is repeated in
 569             // slidingInterface::modifyMotionPoints.  This is done for
 570             // efficiency reasons and avoids multiple creation of cutting
 571             // planes.  Please update both simultaneously.
 572
 573             const point& a = projectedSlavePoints[curEdge.start()];
 574             const point& b = projectedSlavePoints[curEdge.end()];
 575
 576             point c =
 577                 0.5*
 578                 (
 579                     slaveLocalPoints[curEdge.start()]
 580                   + slavePointNormals[curEdge.start()] // Add start normal
 581                   + slaveLocalPoints[curEdge.end()]
 582                   + slavePointNormals[curEdge.end()] // Add end normal
 583                 );
 584
 585             // Create the plane
 586             plane cutPlane(a, b, c);
 587 //             Pout << "a: " << a << " b: " << b << " c: " << c << " plane: " << cutPlane << endl;
 588
 589             linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
 590             const scalar curSlaveLineMag = curSlaveLine.mag();
 591 //             Pout << "curSlaveLine: " << curSlaveLine << endl;
 592             forAll (curMasterEdges, masterEdgeI)
 593             {
 594                 if (!curUme.found(curMasterEdges[masterEdgeI]))
 595                 {
 596                     // New edge
 597                     if (debug)
 598                     {
 599                         Pout << "n";
 600                     }
 601
 602                     const label cmeIndex = curMasterEdges[masterEdgeI];
 603                     const edge& cme = masterEdges[cmeIndex];
 604 //                     Pout<< "Edge " << cmeIndex << " cme: " << cme << " line: " << cme.line(masterLocalPoints) << endl;
 605                     scalar cutOnMaster =
 606                         cutPlane.lineIntersect
 607                         (
 608                             cme.line(masterLocalPoints)
 609                         );
 610
 611                     if
 612                     (
 613                         cutOnMaster > edgeEndCutoffTol_
 614                      && cutOnMaster < 1.0 - edgeEndCutoffTol_
 615                     )
 616                     {
 617                         // Master is cut, check the slave
 618                         point masterCutPoint =
 619                             masterLocalPoints[cme.start()]
 620                           + cutOnMaster*cme.vec(masterLocalPoints);
 621
 622                         pointHit slaveCut =
 623                             curSlaveLine.nearestDist(masterCutPoint);
 624
 625                         if (slaveCut.hit())
 626                         {
 627                             // Strict checking of slave cut to avoid capturing
 628                             // end points.
 629                             scalar cutOnSlave =
 630                                 (
 631                                     (
 632                                         slaveCut.hitPoint()
 633                                       - curSlaveLine.start()
 634                                     ) & curSlaveLine.vec()
 635                                 )/sqr(curSlaveLineMag);
 636
 637                             // Calculate merge tolerance from the
 638                             // target edge length
 639                             scalar mergeTol =
 640                                 edgeCoPlanarTol_*mag(b - a);
 641 //                             Pout << "cutOnMaster: " << cutOnMaster << " masterCutPoint: " << masterCutPoint << " slaveCutPoint: " << slaveCut.hitPoint() << " slaveCut.distance(): " << slaveCut.distance() << " slave length: " << mag(b - a) << " mergeTol: " << mergeTol << " 1: " << mag(b - a) << " 2: " << cme.line(masterLocalPoints).mag() << endl;
 642                             if
 643                             (
 644                                 cutOnSlave > edgeEndCutoffTol_
 645                              && cutOnSlave < 1.0 - edgeEndCutoffTol_
 646                              && slaveCut.distance() < mergeTol
 647                             )
 648                             {
 649                                 // Cut both master and slave.  Add point
 650                                 // to edge points The point is nominally
 651                                 // added from the start of the master edge
 652                                 // and added to the cut point zone
 653                                 label newPoint =
 654                                     ref.setAction
 655                                     (
 656                                         polyAddPoint
 657                                         (
 658                                             masterCutPoint,           // point
 659                                             masterMeshPoints[cme.start()],// m p
 660                                             cutPointZoneID_.index(),  // zone
 661                                             true                      // active
 662                                         )
 663                                     );
 664 //                                 Pout << "Inserting point: " << newPoint << " as edge to edge intersection.  Slave edge: " << edgeI << " " << curEdge << " master edge: " << cmeIndex << " " << cme << endl;
 665                                 pointsIntoSlaveEdges[edgeI].append(newPoint);
 666                                 pointsIntoMasterEdges[cmeIndex].append(newPoint);
 667
 668                                 // Add the point into the enriched patch map
 669                                 pointMap.insert
 670                                 (
 671                                     newPoint,
 672                                     masterCutPoint
 673                                 );
 674
 675                                 // Record which two edges intersect to
 676                                 // create cut point
 677                                 addToCpepm.insert
 678                                 (
 679                                     newPoint,    // Cut point index
 680                                     Pair<edge>
 681                                     (
 682                                         edge
 683                                         (
 684                                             masterMeshPoints[cme.start()],
 685                                             masterMeshPoints[cme.end()]
 686                                         ),    // Master edge
 687                                         edge
 688                                         (
 689                                             slaveMeshPoints[curEdge.start()],
 690                                             slaveMeshPoints[curEdge.end()]
 691                                         )// Slave edge
 692                                     )
 693                                 );
 694
 695                                 if (debug)
 696                                 {
 697                                     Pout<< " " << newPoint << " = "
 698                                         << masterCutPoint << " ";
 699                                 }
 700                             }
 701                             else
 702                             {
 703                                 if (debug)
 704                                 {
 705                                     // Intersection exists but it is too far
 706                                     Pout << "t";
 707                                 }
 708                             }
 709                         }
 710                         else
 711                         {
 712                             if (debug)
 713                             {
 714                                 // Missed slave edge
 715                                 Pout << "x";
 716                             }
 717                         }
 718                     }
 719                     else
 720                     {
 721                         if (debug)
 722                         {
 723                             // Missed master edge
 724                             Pout << "-";
 725                         }
 726                     }
 727                 }
 728                 else
 729                 {
 730                     if (debug)
 731                     {
 732                         Pout << "u";
 733                     }
 734                 }
 735             }
 736
 737             if (debug)
 738             {
 739                 Pout << endl;
 740             }
 741         } // End if both ends missing
 742     } // End for all slave edges
 743
 744     // Re-pack the points into edges lists
 745     labelListList pime(pointsIntoMasterEdges.size());
 746
 747     forAll (pointsIntoMasterEdges, i)
 748     {
 749         pime[i].transfer(pointsIntoMasterEdges[i].shrink());
 750     }
 751
 752     labelListList pise(pointsIntoSlaveEdges.size());
 753
 754     forAll (pointsIntoSlaveEdges, i)
 755     {
 756         pise[i].transfer(pointsIntoSlaveEdges[i].shrink());
 757     }
 758
 759     // Prepare the enriched faces
 760     cutPatch.calcEnrichedFaces
 761     (
 762         pime,
 763         pise,
 764         projectedSlavePoints
 765     );
 766
 767     // Demand driven calculate the cut faces. Apart from the
 768     // cutFaces/cutFaceMaster/cutFaceSlave no information from the cutPatch
 769     // is used anymore!
 770     const faceList& cutFaces = cutPatch.cutFaces();
 771     const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
 772     const labelList& cutFaceSlave = cutPatch.cutFaceSlave();
 773
 774     const labelList& masterFc = masterFaceCells();
 775     const labelList& slaveFc = slaveFaceCells();
 776
 777     // Couple the interface
 778
 779     // Algorithm:
 780     // 1) Go through the cut faces and check if the cut face is the same as the
 781     //    defining master or slave.  If so, modify the appropriate
 782     //    face and mark the other for relegation into the face zone.
 783     // 2) If different, mark both sides for relegation and insert the new face
 784
 785
 786     boolList orphanedMaster(masterPatch.size(), false);
 787     boolList orphanedSlave(slavePatch.size(), false);
 788
 789     forAll (cutFaces, faceI)
 790     {
 791         const face& curCutFace = cutFaces[faceI];
 792         const label curMaster = cutFaceMaster[faceI];
 793         const label curSlave = cutFaceSlave[faceI];
 794
 795 //         Pout << "Doing insertion of face " << faceI << ": ";
 796
 797         // Check if the face has changed topologically
 798         bool insertedFace = false;
 799
 800         if (curMaster >= 0)
 801         {
 802             // Face has got a master
 803             if (curCutFace == masterPatch[curMaster])
 804             {
 805                 // Face is equal to master.  Modify master face.
 806 //                 Pout << "Face is equal to master and is ";
 807
 808                 // If the face has got both master and slave, it is an
 809                 // internal face; otherwise it is a patch face in the
 810                 // master patch.  Keep it in the master face zone.
 811
 812                 if (curSlave >= 0)
 813                 {
 814 //                     Pout << "internal" << endl;
 815                     if (masterFc[curMaster] < slaveFc[curSlave])
 816                     {
 817                         // Cut face should point into slave.
 818                         // Be careful about flips in zone!
 819                         ref.setAction
 820                         (
 821                             polyModifyFace
 822                             (
 823                                 curCutFace,                  // new face
 824                                 masterPatchAddr[curMaster],  // master face id
 825                                 masterFc[curMaster],         // owner
 826                                 slaveFc[curSlave],           // neighbour
 827                                 false,                       // flux flip
 828                                 -1,                          // patch ID
 829                                 false,                       // remove from zone
 830                                 masterFaceZoneID_.index(),   // zone ID
 831                                 masterPatchFlip[curMaster]   // zone flip
 832                             )
 833                         );
 834 //                         Pout << "modifying master face. Old master: " << masterPatch[curMaster] << " new face: " << curCutFace.reverseFace() << " own: " << masterFc[curMaster] << " nei: " << slaveFc[curSlave] << endl;
 835                     }
 836                     else
 837                     {
 838                         // Cut face should point into master.  Flip required.
 839                         // Be careful about flips in zone!
 840                         ref.setAction
 841                         (
 842                             polyModifyFace
 843                             (
 844                                 curCutFace.reverseFace(),    // new face
 845                                 masterPatchAddr[curMaster],  // master face id
 846                                 slaveFc[curSlave],           // owner
 847                                 masterFc[curMaster],         // neighbour
 848                                 true,                        // flux flip
 849                                 -1,                          // patch ID
 850                                 false,                       // remove from zone
 851                                 masterFaceZoneID_.index(),   // zone ID
 852                                 !masterPatchFlip[curMaster]  // zone flip
 853                             )
 854                         );
 855                     }
 856
 857                     // Orphan the slave
 858                     orphanedSlave[curSlave] = true;
 859                 }
 860                 else
 861                 {
 862 //                     Pout << "master boundary" << endl;
 863                     ref.setAction
 864                     (
 865                         polyModifyFace
 866                         (
 867                             curCutFace,                  // new face
 868                             masterPatchAddr[curMaster],  // master face index
 869                             masterFc[curMaster],         // owner
 870                             -1,                          // neighbour
 871                             false,                       // flux flip
 872                             masterPatchID_.index(),      // patch ID
 873                             false,                       // remove from zone
 874                             masterFaceZoneID_.index(),   // zone ID
 875                             masterPatchFlip[curMaster]   // zone flip
 876                         )
 877                     );
 878                 }
 879
 880                 insertedFace = true;
 881             }
 882         }
 883         else if (curSlave >= 0)
 884         {
 885             // Face has got a slave
 886
 887             // Renumber the slave face into merged labels
 888             face rsf(slavePatch[curSlave]);
 889
 890             forAll (rsf, i)
 891             {
 892                 Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);
 893
 894                 if (mpIter != pointMergeMap.end())
 895                 {
 896                     rsf[i] = mpIter();
 897                 }
 898             }
 899
 900             if (curCutFace == rsf)
 901             {
 902                 // Face is equal to slave.  Modify slave face.
 903 //                 Pout << "Face is equal to slave and is ";
 904
 905                 if (curMaster >= 0)
 906                 {
 907 //                     Pout << "regular internal" << endl;
 908                     if (masterFc[curMaster] < slaveFc[curSlave])
 909                     {
 910                         ref.setAction
 911                         (
 912                             polyModifyFace
 913                             (
 914                                 curCutFace,                  // new face
 915                                 slavePatchAddr[curSlave],    // master face id
 916                                 masterFc[curMaster],         // owner
 917                                 slaveFc[curSlave],           // neighbour
 918                                 true,                        // flux flip
 919                                 -1,                          // patch ID
 920                                 false,                       // remove from zone
 921                                 slaveFaceZoneID_.index(),    // zone ID
 922                                 !slavePatchFlip[curMaster]   // zone flip
 923                             )
 924                         );
 925                     }
 926                     else
 927                     {
 928                         // Cut face should point into master.
 929                         // Be careful about flips in zone!
 930 //                         Pout << "flipped internal" << endl;
 931                         ref.setAction
 932                         (
 933                             polyModifyFace
 934                             (
 935                                 curCutFace.reverseFace(),    // new face
 936                                 slavePatchAddr[curSlave],    // master face id
 937                                 slaveFc[curSlave],           // owner
 938                                 masterFc[curMaster],         // neighbour
 939                                 false,                       // flux flip
 940                                 -1,                          // patch ID
 941                                 false,                       // remove from zone
 942                                 slaveFaceZoneID_.index(),    // zone ID
 943                                 slavePatchFlip[curSlave]     // zone flip
 944                             )
 945                         );
 946                     }
 947
 948                     // Orphan the master
 949                     orphanedMaster[curMaster] = true;
 950                 }
 951                 else
 952                 {
 953 //                     Pout << "slave boundary" << endl;
 954                     ref.setAction
 955                     (
 956                         polyModifyFace
 957                         (
 958                             curCutFace,                  // new face
 959                             slavePatchAddr[curSlave],    // master face index
 960                             slaveFc[curSlave],           // owner
 961                             -1,                          // neighbour
 962                             false,                       // flux flip
 963                             slavePatchID_.index(),       // patch ID
 964                             false,                       // remove from zone
 965                             slaveFaceZoneID_.index(),    // zone ID
 966                             slavePatchFlip[curSlave]     // zone flip
 967                         )
 968                     );
 969                 }
 970
 971                 insertedFace = true;
 972             }
 973         }
 974         else
 975         {
 976             FatalErrorIn
 977             (
 978                 "void slidingInterface::coupleInterface("
 979                 "polyTopoChange& ref) const"
 980             )   << "Face " << faceI << " in cut faces has neither a master "
 981                 << "nor a slave.  Error in the cutting algorithm on modify."
 982                 << abort(FatalError);
 983         }
 984
 985         if (!insertedFace)
 986         {
 987             // Face is different from both master and slave
 988 //             Pout << "Face different from both master and slave" << endl;
 989
 990             if (curMaster >= 0)
 991             {
 992                 if (curSlave >= 0)
 993                 {
 994                     // Add internal face
 995                     if (masterFc[curMaster] < slaveFc[curSlave])
 996                     {
 997 //                         Pout << "Adding internal face " << curCutFace << " owner: " << masterFc[curMaster] << " slave: " << slaveFc[curSlave] << " master face: " << masterPatchAddr[curMaster] << endl;
 998                         // Cut face should point into slave.
 999                         ref.setAction
1000                         (
1001                             polyAddFace
1002                             (
1003                                 curCutFace,                  // new face
1004                                 masterFc[curMaster],         // owner
1005                                 slaveFc[curSlave],           // neighbour
1006                                 -1,                          // master point
1007                                 -1,                          // master edge
1008                                 masterPatchAddr[curMaster],  // master face id
1009                                 false,                       // flux flip
1010                                 -1,                          // patch ID
1011                                 cutFaceZoneID_.index(),      // zone ID
1012                                 false                        // zone flip
1013                             )
1014                         );
1015                     }
1016                     else
1017                     {
1018                         // Cut face should point into master.  Flip required.
1019                         ref.setAction
1020                         (
1021                             polyAddFace
1022                             (
1023                                 curCutFace.reverseFace(),    // new face
1024                                 slaveFc[curSlave],           // owner
1025                                 masterFc[curMaster],         // neighbour
1026                                 -1,                          // master point
1027                                 -1,                          // master edge
1028                                 masterPatchAddr[curMaster],  // master face id
1029                                 true,                        // flux flip
1030                                 -1,                          // patch ID
1031                                 cutFaceZoneID_.index(),      // zone ID
1032                                 true                         // zone flip
1033                             )
1034                         );
1035                     }
1036
1037                     // Orphan slave
1038                     orphanedSlave[curSlave] = true;
1039                 }
1040                 else
1041                 {
1042 //                 Pout << "Adding solo master face " << curCutFace << " owner: " << masterFc[curMaster] << " master face: " << masterPatchAddr[curMaster] << endl;
1043                     // Add master patch face
1044                     ref.setAction
1045                     (
1046                         polyAddFace
1047                         (
1048                             curCutFace,                  // new face
1049                             masterFc[curMaster],         // owner
1050                             -1,                          // neighbour
1051                             -1,                          // master point
1052                             -1,                          // master edge
1053                             masterPatchAddr[curMaster],  // master face index
1054                             false,                       // flux flip
1055                             masterPatchID_.index(),      // patch ID
1056                             cutFaceZoneID_.index(),      // zone ID
1057                             false                        // zone flip
1058                         )
1059                     );
1060                 }
1061
1062                 // Orphan master
1063                 orphanedMaster[curMaster] = true;
1064             }
1065             else if (curSlave >= 0)
1066             {
1067 //                 Pout << "Adding solo slave face " << curCutFace << " owner: " << slaveFc[curSlave] << " master face: " << slavePatchAddr[curSlave] << endl;
1068                 // Add slave patch face
1069                 ref.setAction
1070                 (
1071                     polyAddFace
1072                     (
1073                         curCutFace,                  // new face
1074                         slaveFc[curSlave],           // owner
1075                         -1,                          // neighbour
1076                         -1,                          // master point
1077                         -1,                          // master edge
1078                         slavePatchAddr[curSlave],    // master face index
1079                         false,                       // flux flip
1080                         slavePatchID_.index(),       // patch ID
1081                         cutFaceZoneID_.index(),      // zone ID
1082                         false                        // zone flip
1083                     )
1084                 );
1085
1086                 // Orphan slave
1087                 orphanedSlave[curSlave] = true;
1088             }
1089             else
1090             {
1091                 FatalErrorIn
1092                 (
1093                     "void slidingInterface::coupleInterface("
1094                     "polyTopoChange& ref) const"
1095                 )   << "Face " << faceI << " in cut faces has neither a master "
1096                     << "nor a slave.  Error in the cutting algorithm on add."
1097                     << abort(FatalError);
1098             }
1099         }
1100     }
1101
1102     // Move the orphaned faces into the face zone
1103 //     Pout << "Orphaned master faces: " << orphanedMaster << endl;
1104 //     Pout << "Orphaned slave faces: " << orphanedSlave << endl;
1105
1106     label nOrphanedMasters = 0;
1107
1108     forAll (orphanedMaster, faceI)
1109     {
1110         if (orphanedMaster[faceI])
1111         {
1112             nOrphanedMasters++;
1113
1114             //// Recover original orientation
1115             //ref.setAction
1116             //(
1117             //    polyModifyFace
1118             //    (
1119             //        masterPatch[faceI],                 // new face
1120             //        masterPatchAddr[faceI],             // master face index
1121             //        -1,                                 // owner
1122             //        -1,                                 // neighbour
1123             //        false,                              // flux flip
1124             //        -1,                                 // patch ID
1125             //        false,                              // remove from zone
1126             //        masterFaceZoneID_.index(),          // zone ID
1127             //        false                               // zone flip
1128             //    )
1129             //);
1130
1131             //Pout<< "**MJ:deleting master face " << masterPatchAddr[faceI]
1132             //    << " old verts:" << masterPatch[faceI] << endl;
1133             ref.setAction(polyRemoveFace(masterPatchAddr[faceI]));
1134         }
1135     }
1136
1137     label nOrphanedSlaves = 0;
1138
1139     forAll (orphanedSlave, faceI)
1140     {
1141         if (orphanedSlave[faceI])
1142         {
1143             nOrphanedSlaves++;
1144
1145             //// Recover original orientation
1146             //ref.setAction
1147             //(
1148             //    polyModifyFace
1149             //    (
1150             //        slavePatch[faceI],                // new face
1151             //        slavePatchAddr[faceI],            // slave face index
1152             //        -1,                               // owner
1153             //        -1,                               // neighbour
1154             //        false,                            // flux flip
1155             //        -1,                               // patch ID
1156             //        false,                            // remove from zone
1157             //        slaveFaceZoneID_.index(),         // zone ID
1158             //        false                             // zone flip
1159             //    )
1160             //);
1161
1162             //Pout<< "**MJ:deleting slave face " << slavePatchAddr[faceI]
1163             //    << " old verts:" << slavePatch[faceI] << endl;
1164             ref.setAction(polyRemoveFace(slavePatchAddr[faceI]));
1165         }
1166     }
1167
1168     if (debug)
1169     {
1170         Pout<< "Number of orphaned faces: "
1171             << "master = " << nOrphanedMasters << " out of "
1172             << orphanedMaster.size()
1173             << " slave = " << nOrphanedSlaves << " out of "
1174             << orphanedSlave.size() << endl;
1175     }
1176
1177     // Finished coupling the plane of sliding interface.
1178
1179     // Modify faces influenced by the sliding interface
1180     // These are the faces belonging to the master and slave cell
1181     // layer which have not already been modified.
1182     // Modification comes in three types:
1183     // 1) eliminate the points that have been removed when the old sliding
1184     //    interface has been removed
1185     // 2) Merge the slave points that have hit points on the master patch
1186     // 3) Introduce new points resulting from the new sliding interface cut
1187
1188     // Collect all affected faces
1189
1190     // Master side
1191
1192     // Grab the list of faces in the layer
1193     const labelList& masterStickOuts = masterStickOutFaces();
1194
1195 //     Pout << "masterStickOuts: " << masterStickOuts << endl;
1196
1197     // Re-create the master stick-out faces
1198     forAll (masterStickOuts, faceI)
1199     {
1200         // Renumber the face and remove additional points
1201
1202         const label curFaceID = masterStickOuts[faceI];
1203
1204         const face& oldRichFace = faces[curFaceID];
1205
1206         bool changed = false;
1207
1208         // Remove removed points from the face
1209         face oldFace(oldRichFace.size());
1210         label nOldFace = 0;
1211
1212         forAll (oldRichFace, pointI)
1213         {
1214             if (ref.pointRemoved(oldRichFace[pointI]))
1215             {
1216                 changed = true;
1217             }
1218             else
1219             {
1220                 // Point off patch
1221                 oldFace[nOldFace] = oldRichFace[pointI];
1222                 nOldFace++;
1223             }
1224         }
1225
1226         oldFace.setSize(nOldFace);
1227
1228 //         Pout << "old rich master face: " << oldRichFace << " old face: " << oldFace << endl;
1229         DynamicList<label> newFaceLabels(2*oldFace.size());
1230
1231         forAll (oldFace, pointI)
1232         {
1233             if (masterMeshPointMap.found(oldFace[pointI]))
1234             {
1235                 // Point is in master patch. Add it
1236
1237                 // If the point is a direct hit, grab its label; otherwise
1238                 // keep the original
1239                 if (pointMergeMap.found(oldFace[pointI]))
1240                 {
1241                     changed = true;
1242                     newFaceLabels.append
1243                     (
1244                         pointMergeMap.find(oldFace[pointI])()
1245                     );
1246                 }
1247                 else
1248                 {
1249                     newFaceLabels.append(oldFace[pointI]);
1250                 }
1251
1252                 // Find if there are additional points inserted onto the edge
1253                 // between current point and the next point
1254                 // Algorithm:
1255                 // 1) Find all the edges in the master patch coming
1256                 //    out of the current point.
1257                 // 2) If the next point in the face to pick the right edge
1258                 const label localFirstLabel =
1259                     masterMeshPointMap.find(oldFace[pointI])();
1260
1261                 const labelList& curEdges = masterPointEdges[localFirstLabel];
1262
1263                 const label  nextLabel = oldFace.nextLabel(pointI);
1264
1265                 Map<label>::const_iterator mmpmIter =
1266                     masterMeshPointMap.find(nextLabel);
1267
1268                 if (mmpmIter != masterMeshPointMap.end())
1269                 {
1270 //                     Pout << "found label pair " << oldFace[pointI] << " and " << nextLabel;
1271                     // Find the points on the edge between them
1272                     const label localNextLabel = mmpmIter();
1273
1274                     forAll (curEdges, curEdgeI)
1275                     {
1276                         if
1277                         (
1278                             masterEdges[curEdges[curEdgeI]].otherVertex
1279                             (
1280                                 localFirstLabel
1281                             )
1282                          == localNextLabel
1283                         )
1284                         {
1285 //                             Pout << " found edge: " << curEdges[curEdgeI] << endl;
1286
1287                             // Get points on current edge
1288                             const labelList& curPime = pime[curEdges[curEdgeI]];
1289
1290                             if (curPime.size() > 0)
1291                             {
1292                                 changed = true;
1293                                 // Pout << "curPime: " << curPime << endl;
1294                                 // Insert the edge points into the face
1295                                 // in the correct order
1296                                 const point& startPoint =
1297                                     masterLocalPoints[localFirstLabel];
1298
1299                                 vector e =
1300                                     masterLocalPoints[localNextLabel]
1301                                   - startPoint;
1302
1303                                 e /= magSqr(e);
1304
1305                                 scalarField edgePointWeights(curPime.size());
1306
1307                                 forAll (curPime, curPimeI)
1308                                 {
1309                                     edgePointWeights[curPimeI] =
1310                                         (
1311                                             e
1312                                           & (
1313                                               pointMap.find(curPime[curPimeI])()
1314                                             - startPoint
1315                                             )
1316                                         );
1317                                 }
1318
1319                                 if (debug)
1320                                 {
1321                                     if
1322                                     (
1323                                         min(edgePointWeights) < 0
1324                                      || max(edgePointWeights) > 1
1325                                     )
1326                                     {
1327                                         FatalErrorIn
1328                                         (
1329                                             "void slidingInterface::"
1330                                             "coupleInterface("
1331                                             "polyTopoChange& ref) const"
1332                                         )   << "Error in master stick-out edge "
1333                                             << "point collection."
1334                                             << abort(FatalError);
1335                                     }
1336                                 }
1337
1338                                 // Go through the points and collect
1339                                 // them based on weights from lower to
1340                                 // higher.  This gives the correct
1341                                 // order of points along the edge.
1342                                 for
1343                                 (
1344                                     label passI = 0;
1345                                     passI < edgePointWeights.size();
1346                                     passI++
1347                                 )
1348                                 {
1349                                     // Max weight can only be one, so
1350                                     // the sorting is done by
1351                                     // elimination.
1352                                     label nextPoint = -1;
1353                                     scalar dist = 2;
1354
1355                                     forAll (edgePointWeights, wI)
1356                                     {
1357                                         if (edgePointWeights[wI] < dist)
1358                                         {
1359                                             dist = edgePointWeights[wI];
1360                                             nextPoint = wI;
1361                                         }
1362                                     }
1363
1364                                     // Insert the next point and reset
1365                                     // its weight to exclude it from
1366                                     // future picks
1367                                     newFaceLabels.append(curPime[nextPoint]);
1368                                     edgePointWeights[nextPoint] = GREAT;
1369                                 }
1370                             }
1371
1372                             break;
1373                         } // End of found edge
1374                     } // End of looking through current edges
1375                 }
1376             }
1377             else
1378             {
1379                 newFaceLabels.append(oldFace[pointI]);
1380             }
1381         }
1382
1383         if (changed)
1384         {
1385             if (newFaceLabels.size() < 3)
1386             {
1387                 FatalErrorIn
1388                 (
1389                     "void slidingInterface::coupleInterface("
1390                     "polyTopoChange& ref) const"
1391                 )   << "Face " << curFaceID << " reduced to less than "
1392                     << "3 points.  Topological/cutting error A." << nl
1393                     << "Old face: " << oldFace << " new face: " << newFaceLabels
1394                     << abort(FatalError);
1395             }
1396
1397             // Get face zone and its flip
1398             label modifiedFaceZone = faceZones.whichZone(curFaceID);
1399             bool modifiedFaceZoneFlip = false;
1400
1401             if (modifiedFaceZone >= 0)
1402             {
1403                 modifiedFaceZoneFlip =
1404                     faceZones[modifiedFaceZone].flipMap()
1405                     [
1406                         faceZones[modifiedFaceZone].whichFace(curFaceID)
1407                     ];
1408             }
1409
1410             face newFace;
1411             newFace.transfer(newFaceLabels.shrink());
1412
1413 //             Pout << "Modifying master stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1414
1415             // Modify the face
1416             if (mesh.isInternalFace(curFaceID))
1417             {
1418                 ref.setAction
1419                 (
1420                     polyModifyFace
1421                     (
1422                         newFace,                // modified face
1423                         curFaceID,              // label of face being modified
1424                         own[curFaceID],         // owner
1425                         nei[curFaceID],         // neighbour
1426                         false,                  // face flip
1427                         -1,                     // patch for face
1428                         false,                  // remove from zone
1429                         modifiedFaceZone,       // zone for face
1430                         modifiedFaceZoneFlip    // face flip in zone
1431                     )
1432                 );
1433             }
1434             else
1435             {
1436                 ref.setAction
1437                 (
1438                     polyModifyFace
1439                     (
1440                         newFace,                // modified face
1441                         curFaceID,              // label of face being modified
1442                         own[curFaceID],         // owner
1443                         -1,                     // neighbour
1444                         false,                  // face flip
1445                         mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1446                         false,                  // remove from zone
1447                         modifiedFaceZone,       // zone for face
1448                         modifiedFaceZoneFlip    // face flip in zone
1449                     )
1450                 );
1451             }
1452         }
1453     }
1454
1455 //     Pout << "Finished master side" << endl;
1456
1457     // Slave side
1458
1459     // Grab the list of faces in the layer
1460     const labelList& slaveStickOuts = slaveStickOutFaces();
1461
1462 //     Pout << "slaveStickOuts: " << slaveStickOuts << endl;
1463
1464     const Map<label>& rpm = retiredPointMap();
1465
1466     // Re-create the slave stick-out faces
1467
1468     forAll (slaveStickOuts, faceI)
1469     {
1470         // Renumber the face and remove additional points
1471         const label curFaceID = slaveStickOuts[faceI];
1472
1473         const face& oldRichFace = faces[curFaceID];
1474
1475         bool changed = false;
1476
1477         // Remove removed points from the face
1478         face oldFace(oldRichFace.size());
1479         label nOldFace = 0;
1480
1481         forAll (oldRichFace, pointI)
1482         {
1483             if
1484             (
1485                 rpm.found(oldRichFace[pointI])
1486              || slaveMeshPointMap.found(oldRichFace[pointI])
1487             )
1488             {
1489                 // Point definitely live. Add it
1490                 oldFace[nOldFace] = oldRichFace[pointI];
1491                 nOldFace++;
1492             }
1493             else if
1494             (
1495                 ref.pointRemoved(oldRichFace[pointI])
1496              || masterMeshPointMap.found(oldRichFace[pointI])
1497             )
1498             {
1499                 // Point removed and not on slave patch
1500                 // (first if takes care of that!) or
1501                 // point belonging to master patch
1502                 changed = true;
1503             }
1504             else
1505             {
1506                 // Point off patch
1507                 oldFace[nOldFace] = oldRichFace[pointI];
1508                 nOldFace++;
1509             }
1510         }
1511
1512         oldFace.setSize(nOldFace);
1513
1514         DynamicList<label> newFaceLabels(2*oldFace.size());
1515
1516 //         Pout << "old rich slave face: " << oldRichFace << " old face: " << oldFace << endl;
1517         forAll (oldFace, pointI)
1518         {
1519             // Try to find the point in retired points
1520             label curP = oldFace[pointI];
1521
1522             Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointI]);
1523
1524             if (rpmIter != rpm.end())
1525             {
1526                 changed = true;
1527                 curP = rpmIter();
1528             }
1529
1530             if (slaveMeshPointMap.found(curP))
1531             {
1532                 // Point is in slave patch. Add it
1533
1534                 // If the point is a direct hit, grab its label; otherwise
1535                 // keep the original
1536                 if (pointMergeMap.found(curP))
1537                 {
1538                     changed = true;
1539                     newFaceLabels.append
1540                     (
1541                         pointMergeMap.find(curP)()
1542                     );
1543                 }
1544                 else
1545                 {
1546                     newFaceLabels.append(curP);
1547                 }
1548
1549                 // Find if there are additional points inserted onto the edge
1550                 // between current point and the next point
1551                 // Algorithm:
1552                 // 1) Find all the edges in the slave patch coming
1553                 //    out of the current point.
1554                 // 2) Use the next point in the face to pick the right edge
1555
1556                 const label localFirstLabel =
1557                     slaveMeshPointMap.find(curP)();
1558
1559                 const labelList& curEdges = slavePointEdges[localFirstLabel];
1560
1561                 label nextLabel = oldFace.nextLabel(pointI);
1562
1563                 Map<label>::const_iterator rpmNextIter =
1564                     rpm.find(nextLabel);
1565
1566                 if (rpmNextIter != rpm.end())
1567                 {
1568                     nextLabel = rpmNextIter();
1569                 }
1570
1571                 Map<label>::const_iterator mmpmIter =
1572                     slaveMeshPointMap.find(nextLabel);
1573
1574                 if (mmpmIter != slaveMeshPointMap.end())
1575                 {
1576                     // Both points on the slave patch.
1577                     // Find the points on the edge between them
1578                     const label localNextLabel = mmpmIter();
1579
1580                     forAll (curEdges, curEdgeI)
1581                     {
1582                         if
1583                         (
1584                             slaveEdges[curEdges[curEdgeI]].otherVertex
1585                             (
1586                                 localFirstLabel
1587                             )
1588                          == localNextLabel
1589                         )
1590                         {
1591 //                             Pout << " found edge: " << curEdges[curEdgeI] << endl;
1592
1593                             // Get points on current edge
1594                             const labelList& curPise = pise[curEdges[curEdgeI]];
1595
1596                             if (curPise.size() > 0)
1597                             {
1598                                 changed = true;
1599 //                                 Pout << "curPise: " << curPise << endl;
1600                                 // Insert the edge points into the face
1601                                 // in the correct order
1602                                 const point& startPoint =
1603                                     projectedSlavePoints[localFirstLabel];
1604
1605                                 vector e =
1606                                     projectedSlavePoints[localNextLabel]
1607                                   - startPoint;
1608
1609                                 e /= magSqr(e);
1610
1611                                 scalarField edgePointWeights(curPise.size());
1612
1613                                 forAll (curPise, curPiseI)
1614                                 {
1615                                     edgePointWeights[curPiseI] =
1616                                     (
1617                                         e
1618                                       & (
1619                                             pointMap.find(curPise[curPiseI])()
1620                                           - startPoint
1621                                         )
1622                                     );
1623                                 }
1624
1625                                 if (debug)
1626                                 {
1627                                     if
1628                                     (
1629                                         min(edgePointWeights) < 0
1630                                      || max(edgePointWeights) > 1
1631                                     )
1632                                     {
1633                                         FatalErrorIn
1634                                         (
1635                                             "void slidingInterface::"
1636                                             "coupleInterface("
1637                                             "polyTopoChange& ref) const"
1638                                         )   << "Error in slave stick-out edge "
1639                                             << "point collection."
1640                                             << abort(FatalError);
1641                                         }
1642                                     }
1643
1644                                 // Go through the points and collect
1645                                 // them based on weights from lower to
1646                                 // higher.  This gives the correct
1647                                 // order of points along the edge.
1648                                 for
1649                                 (
1650                                     label passI = 0;
1651                                     passI < edgePointWeights.size();
1652                                     passI++
1653                                 )
1654                                 {
1655                                     // Max weight can only be one, so
1656                                     // the sorting is done by
1657                                     // elimination.
1658                                     label nextPoint = -1;
1659                                     scalar dist = 2;
1660
1661                                     forAll (edgePointWeights, wI)
1662                                     {
1663                                         if (edgePointWeights[wI] < dist)
1664                                         {
1665                                             dist = edgePointWeights[wI];
1666                                             nextPoint = wI;
1667                                         }
1668                                     }
1669
1670                                     // Insert the next point and reset
1671                                     // its weight to exclude it from
1672                                     // future picks
1673                                     newFaceLabels.append(curPise[nextPoint]);
1674                                     edgePointWeights[nextPoint] = GREAT;
1675                                 }
1676                             }
1677
1678                             break;
1679                         }
1680                     } // End of found edge
1681                 } // End of looking through current edges
1682             }
1683             else
1684             {
1685                 newFaceLabels.append(oldFace[pointI]);
1686             }
1687         }
1688
1689         if (changed)
1690         {
1691             if (newFaceLabels.size() < 3)
1692             {
1693                 FatalErrorIn
1694                 (
1695                     "void slidingInterface::coupleInterface("
1696                     "polyTopoChange& ref) const"
1697                 )   << "Face " << curFaceID << " reduced to less than "
1698                     << "3 points.  Topological/cutting error B." << nl
1699                     << "Old face: " << oldFace << " new face: " << newFaceLabels
1700                     << abort(FatalError);
1701             }
1702
1703             // Get face zone and its flip
1704             label modifiedFaceZone = faceZones.whichZone(curFaceID);
1705             bool modifiedFaceZoneFlip = false;
1706
1707             if (modifiedFaceZone >= 0)
1708             {
1709                 modifiedFaceZoneFlip =
1710                     faceZones[modifiedFaceZone].flipMap()
1711                     [
1712                         faceZones[modifiedFaceZone].whichFace(curFaceID)
1713                     ];
1714             }
1715
1716             face newFace;
1717             newFace.transfer(newFaceLabels.shrink());
1718
1719 //             Pout << "Modifying slave stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;
1720
1721             // Modify the face
1722             if (mesh.isInternalFace(curFaceID))
1723             {
1724                 ref.setAction
1725                 (
1726                     polyModifyFace
1727                     (
1728                         newFace,                // modified face
1729                         curFaceID,              // label of face being modified
1730                         own[curFaceID],         // owner
1731                         nei[curFaceID],         // neighbour
1732                         false,                  // face flip
1733                         -1,                     // patch for face
1734                         false,                  // remove from zone
1735                         modifiedFaceZone,       // zone for face
1736                         modifiedFaceZoneFlip    // face flip in zone
1737                     )
1738                 );
1739             }
1740             else
1741             {
1742                 ref.setAction
1743                 (
1744                     polyModifyFace
1745                     (
1746                         newFace,                // modified face
1747                         curFaceID,              // label of face being modified
1748                         own[curFaceID],         // owner
1749                         -1,                     // neighbour
1750                         false,                  // face flip
1751                         mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
1752                         false,                  // remove from zone
1753                         modifiedFaceZone,       // zone for face
1754                         modifiedFaceZoneFlip    // face flip in zone
1755                     )
1756                 );
1757             }
1758         }
1759     }
1760
1761     // Activate and retire slave patch points
1762     // This needs to be done last, so that the map of removed points
1763     // does not get damaged by point modifications
1764
1765     if (!retiredPointMapPtr_)
1766     {
1767         FatalErrorIn
1768         (
1769             "void slidingInterface::coupleInterface("
1770             "polyTopoChange& ref) const"
1771         )   << "Retired point map pointer not set."
1772             << abort(FatalError);
1773     }
1774
1775     Map<label>& addToRpm = *retiredPointMapPtr_;
1776
1777     // Clear the old map
1778     addToRpm.clear();
1779
1780     label nRetiredPoints = 0;
1781
1782     forAll (slaveMeshPoints, pointI)
1783     {
1784         if (pointMergeMap.found(slaveMeshPoints[pointI]))
1785         {
1786             // Retire the point - only used for supporting the detached
1787             // slave patch
1788             nRetiredPoints++;
1789
1790             //ref.setAction
1791             //(
1792             //    polyModifyPoint
1793             //    (
1794             //        slaveMeshPoints[pointI],             // point ID
1795             //        points[slaveMeshPoints[pointI]],     // point
1796             //        false,                               // remove from zone
1797             //        mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1798             //        false                                // in a cell
1799             //    )
1800             //);
1801             //Pout<< "MJ retire slave point " << slaveMeshPoints[pointI]
1802             //    << " coord " << points[slaveMeshPoints[pointI]]
1803             //    << endl;
1804             ref.setAction
1805             (
1806                 polyRemovePoint
1807                 (
1808                     slaveMeshPoints[pointI]
1809                 )
1810             );
1811
1812             addToRpm.insert
1813             (
1814                 pointMergeMap.find(slaveMeshPoints[pointI])(),
1815                 slaveMeshPoints[pointI]
1816             );
1817         }
1818         else
1819         {
1820             ref.setAction
1821             (
1822                 polyModifyPoint
1823                 (
1824                     slaveMeshPoints[pointI],             // point ID
1825                     points[slaveMeshPoints[pointI]],     // point
1826                     false,                               // remove from zone
1827                     mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
1828                     true                                 // in a cell
1829                 )
1830             );
1831         }
1832     }
1833
1834     if (debug)
1835     {
1836         Pout << "Retired " << nRetiredPoints << " out of "
1837             << slaveMeshPoints.size() << " points." << endl;
1838     }
1839
1840     // Grab cut face master and slave addressing
1841     if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
1842     cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());
1843
1844     if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
1845     cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());
1846
1847     // Finished coupling
1848     attached_ = true;
1849
1850     if (debug)
1851     {
1852         Pout<< "void slidingInterface::coupleInterface("
1853             << "polyTopoChange& ref) : "
1854             << "Finished coupling sliding interface " << name() << endl;
1855     }
1856 }
1857
1858
1859 // ************************************************************************* //