src/dynamicMesh/slidingInterface/coupleSlidingInterface.C

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