src/dynamicMesh/polyMeshModifiers/slidingInterface/coupleSlidingInterface.C

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