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