coupleSlidingInterface.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "slidingInterface.H"
#include <dynamicMesh/polyTopoChange.H>
#include <OpenFOAM/polyMesh.H>
#include <OpenFOAM/primitiveMesh.H>
#include <dynamicMesh/enrichedPatch.H>
#include <OpenFOAM/DynamicList.H>
#include <OpenFOAM/pointHit.H>
#include <OpenFOAM/triPointRef.H>
#include <OpenFOAM/plane.H>
#include <dynamicMesh/polyTopoChanger.H>
#include <dynamicMesh/polyAddPoint.H>
#include <dynamicMesh/polyRemovePoint.H>
#include <dynamicMesh/polyAddFace.H>
#include <dynamicMesh/polyModifyPoint.H>
#include <dynamicMesh/polyModifyFace.H>
#include <dynamicMesh/polyRemoveFace.H>

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

const Foam::scalar Foam::slidingInterface::edgeCoPlanarTolDefault_ = 0.8;

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

// Index of debug signs:
// Index of debug signs:
// . - loop of the edge-to-face interaction detection
// x - reversal of direction in edge-to-face interaction detection
// + - complete edge-to-face interaction detection
// z - incomplete edge-to-face interaction detection.  This may be OK for edges
//     crossing from one to the other side of multiply connected master patch

// e - adding a point on edge
// f - adding a point on face
// . - collecting edges off another face for edge-to-edge cut
// + - finished collection of edges
// * - cut both master and slave
// n - cutting new edge
// t - intersection exists but it is outside of tolerance
// x - missed slave edge in cut
// - - missed master edge in cut
// u - edge already used in cutting

void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
{
    if (debug)
    {
        Pout<< "void slidingInterface::coupleInterface"
            << "(polyTopoChange& ref) : "
            << "Coupling sliding interface " << name() << endl;
    }

    const polyMesh& mesh = topoChanger().mesh();

    const pointField& points = mesh.points();
    const faceList& faces = mesh.faces();

    const labelList& own = mesh.faceOwner();
    const labelList& nei = mesh.faceNeighbour();
    const faceZoneMesh& faceZones = mesh.faceZones();

    const primitiveFacePatch& masterPatch =
        faceZones[masterFaceZoneID_.index()]();

    const labelList& masterPatchAddr = faceZones[masterFaceZoneID_.index()];

    const boolList& masterPatchFlip =
        faceZones[masterFaceZoneID_.index()].flipMap();

    const primitiveFacePatch& slavePatch =
        faceZones[slaveFaceZoneID_.index()]();

    const labelList& slavePatchAddr = faceZones[slaveFaceZoneID_.index()];

    const boolList& slavePatchFlip =
        faceZones[slaveFaceZoneID_.index()].flipMap();

    const edgeList& masterEdges = masterPatch.edges();
    const labelListList& masterPointEdges = masterPatch.pointEdges();
    const labelList& masterMeshPoints = masterPatch.meshPoints();
    const pointField& masterLocalPoints = masterPatch.localPoints();
    const labelListList& masterFaceFaces = masterPatch.faceFaces();
    const labelListList& masterFaceEdges = masterPatch.faceEdges();
    const Map<label>& masterMeshPointMap = masterPatch.meshPointMap();

    const edgeList& slaveEdges = slavePatch.edges();
    const labelListList& slavePointEdges = slavePatch.pointEdges();
    const labelList& slaveMeshPoints = slavePatch.meshPoints();
    const pointField& slaveLocalPoints = slavePatch.localPoints();
    const Map<label>& slaveMeshPointMap = slavePatch.meshPointMap();
    const vectorField& slavePointNormals = slavePatch.pointNormals();

    // Collect projection addressing
    if
    (
        !(
            slavePointPointHitsPtr_
         && slavePointEdgeHitsPtr_
         && slavePointFaceHitsPtr_
         && masterPointEdgeHitsPtr_
        )
    )
    {
        FatalErrorIn
        (
            "void slidingInterface::coupleInterface("
            "polyTopoChange& ref) const"
        )   << "Point projection addressing not available."
            << abort(FatalError);
    }

    const labelList& slavePointPointHits = *slavePointPointHitsPtr_;
    const labelList& slavePointEdgeHits = *slavePointEdgeHitsPtr_;
    const List<objectHit>& slavePointFaceHits = *slavePointFaceHitsPtr_;
    const labelList& masterPointEdgeHits = *masterPointEdgeHitsPtr_;
    const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;

    // Create enriched patch
    enrichedPatch cutPatch
    (
        masterPatch,
        slavePatch,
        slavePointPointHits,
        slavePointEdgeHits,
        slavePointFaceHits
    );

    // Get reference to list of added point for the enriched patch.
    // This will be used for point addition
    Map<point>& pointMap = cutPatch.pointMap();

    // Get reference to the list of merged points
    Map<label>& pointMergeMap = cutPatch.pointMergeMap();

    // Create mapping for every merged point of the slave patch
    forAll (slavePointPointHits, pointI)
    {
        if (slavePointPointHits[pointI] >= 0)
        {
// Pout << "Inserting point merge pair: " << slaveMeshPoints[pointI] << " : " << masterMeshPoints[slavePointPointHits[pointI]] << endl;
            pointMergeMap.insert
            (
                slaveMeshPoints[pointI],
                masterMeshPoints[slavePointPointHits[pointI]]
            );
        }
    }

    // Collect the list of used edges for every slave edge

    List<labelHashSet> usedMasterEdges(slaveEdges.size());

    // Collect of slave point hits
    forAll (slavePointPointHits, pointI)
    {
        // For point hits, add all point-edges from master side to all point
        const labelList& curSlaveEdges = slavePointEdges[pointI];

        if (slavePointPointHits[pointI] > -1)
        {
            const labelList& curMasterEdges =
                masterPointEdges[slavePointPointHits[pointI]];

            // Mark all current master edges as used for all the current slave
            // edges
            forAll (curSlaveEdges, slaveEdgeI)
            {
                labelHashSet& sm = usedMasterEdges[curSlaveEdges[slaveEdgeI]];

                forAll (curMasterEdges, masterEdgeI)
                {
                    sm.insert(curMasterEdges[masterEdgeI]);
                }
            }
        }
        else if (slavePointEdgeHits[pointI] > -1)
        {
            // For edge hits, add the master edge
            forAll (curSlaveEdges, slaveEdgeI)
            {
                usedMasterEdges[curSlaveEdges[slaveEdgeI]].insert
                (
                    slavePointEdgeHits[pointI]
                );
            }
        }
    }

    // Collect off master point hits
    // For every master point that has hit an edge, add all edges coming from
    // that point to the slave edge that has been hit
    forAll (masterPointEdgeHits, masterPointI)
    {
        if (masterPointEdgeHits[masterPointI] > -1)
        {
            const labelList& curMasterEdges = masterPointEdges[masterPointI];

            labelHashSet& sm =
                usedMasterEdges[masterPointEdgeHits[masterPointI]];

            forAll (curMasterEdges, masterEdgeI)
            {
                sm.insert(curMasterEdges[masterEdgeI]);
            }
        }
    }

//     Pout << "used edges: " << endl;
//     forAll (usedMasterEdges, edgeI)
//     {
//         Pout << "edge: " << edgeI << " used: " << usedMasterEdges[edgeI].toc() << endl;
//     }

    // For every master and slave edge make a list of points to be added into
    // that edge.
    List<DynamicList<label> > pointsIntoMasterEdges(masterEdges.size());
    List<DynamicList<label> > pointsIntoSlaveEdges(slaveEdges.size());

    // Add all points from the slave patch that have hit the edge
    forAll (slavePointEdgeHits, pointI)
    {
        if (slavePointEdgeHits[pointI] > -1)
        {
            // Create a new point on the master edge

            point edgeCutPoint =
                masterEdges[slavePointEdgeHits[pointI]].line
                (
                    masterLocalPoints
                ).nearestDist(projectedSlavePoints[pointI]).hitPoint();

            label newPoint =
                ref.setAction
                (
                    polyAddPoint
                    (
                        edgeCutPoint,             // point
                        slaveMeshPoints[pointI],  // master point
                        cutPointZoneID_.index(),  // zone for point
                        true                      // supports a cell
                    )
                );
//             Pout << "Inserting merge pair off edge: " << slaveMeshPoints[pointI] << " " << newPoint << " cut point: " << edgeCutPoint << " orig: " << slaveLocalPoints[pointI] << " proj: " << projectedSlavePoints[pointI] << endl;
            // Add the new edge point into the merge map
            pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);

            pointsIntoMasterEdges[slavePointEdgeHits[pointI]].append
            (
                newPoint
            );

            // Add the point into the enriched patch map
            pointMap.insert
            (
                newPoint,
                edgeCutPoint
            );

            if (debug)
            {
                Pout<< "e";
//                 Pout<< newPoint << " = " << edgeCutPoint << endl;
            }
        }
    }

    if (debug)
    {
        Pout<< endl;
    }

    // Add all points from the slave patch that have hit a face
    forAll (slavePointFaceHits, pointI)
    {
        if
        (
            slavePointPointHits[pointI] < 0
         && slavePointEdgeHits[pointI] < 0
         && slavePointFaceHits[pointI].hit()
        )
        {
            label newPoint =
                ref.setAction
                (
                    polyAddPoint
                    (
                        projectedSlavePoints[pointI],   // point
                        slaveMeshPoints[pointI],        // master point
                        cutPointZoneID_.index(),        // zone for point
                        true                            // supports a cell
                    )
                );
//             Pout << "Inserting merge pair off face: " << slaveMeshPoints[pointI] << " " << newPoint << endl;
            // Add the new edge point into the merge map
            pointMergeMap.insert(slaveMeshPoints[pointI], newPoint);

            // Add the point into the enriched patch map
            pointMap.insert
            (
                newPoint,
                projectedSlavePoints[pointI]
            );

            if (debug)
            {
                Pout<< "f: " << newPoint << " = "
                    << projectedSlavePoints[pointI] << endl;
            }
        }
    }

    forAll (masterPointEdgeHits, pointI)
    {
        if (masterPointEdgeHits[pointI] > -1)
        {
            pointsIntoSlaveEdges[masterPointEdgeHits[pointI]].append
            (
                masterMeshPoints[pointI]
            );
        }
    }

    // Cut all slave edges.
    // Collect all master edges the slave edge interacts with.  Skip
    // all the edges already marked as used.  For every unused edge,
    // calculate the cut and insert the new point into the master and
    // slave edge.
    // For the edge selection algorithm, see, comment in
    // slidingInterfaceProjectPoints.C.
    // Edge cutting algoritm:
    // As the master patch defines the cutting surface, the newly
    // inserted point needs to be on the master edge.  Also, in 3-D
    // the pair of edges generally misses each other rather than
    // intersect.  Therefore, the intersection is calculated using the
    // plane the slave edge defines during projection.  The plane is
    // defined by the centrepoint of the edge in the original
    // configuration and the projected end points.  In case of flat
    // geometries (when the three points are colinear), the plane is
    // defined by the two projected end-points and the slave edge
    // normal used as the in-plane vector.  When the intersection
    // point is created, it is added as a new point for both the
    // master and the slave edge.
    // In order to be able to re-create the points on edges in mesh
    // motion without the topology change, the edge pair used to
    // create the cut point needs to be recorded in
    // cutPointEdgePairMap

    // Note.  "Processing slave edges" code is repeated twice in the
    // sliding intergace coupling in order to allow the point
    // projection to be done separately from the actual cutting.
    // Please change consistently with slidingInterfaceProjectPoints.C
    //
    if (debug)
    {
        Pout << "Processing slave edges " << endl;
    }

    if (!cutPointEdgePairMapPtr_)
    {
        FatalErrorIn
        (
            "void slidingInterface::coupleInterface("
            "polyTopoChange& ref) const"
        )   << "Cut point edge pair map pointer not set."
            << abort(FatalError);
    }

    Map<Pair<edge> >& addToCpepm = *cutPointEdgePairMapPtr_;

    // Clear the old map
    addToCpepm.clear();

    // Create a map of faces the edge can interact with
    labelHashSet curFaceMap
    (
        nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
    );

    labelHashSet addedFaces(2*primitiveMesh::edgesPerFace_);

    forAll (slaveEdges, edgeI)
    {
        const edge& curEdge = slaveEdges[edgeI];

        if
        (
            slavePointFaceHits[curEdge.start()].hit()
         || slavePointFaceHits[curEdge.end()].hit()
        )
        {
            labelHashSet& curUme = usedMasterEdges[edgeI];
//             Pout<< "Doing edge " << edgeI << " curEdge: " << curEdge << " curUme: " << curUme << endl;
            // Clear the maps
            curFaceMap.clear();
            addedFaces.clear();

            // Grab the faces for start and end points.
            const label startFace =
                slavePointFaceHits[curEdge.start()].hitObject();
            const label endFace = slavePointFaceHits[curEdge.end()].hitObject();
//             Pout << "startFace: " << slavePointFaceHits[curEdge.start()] << " endFace: " << slavePointFaceHits[curEdge.end()] << endl;
            // Insert the start face into the list
            curFaceMap.insert(startFace);
            addedFaces.insert(startFace);
//             Pout << "curFaceMap: " << curFaceMap.toc() << endl;
            label nSweeps = 0;
            bool completed = false;

            while (nSweeps < edgeFaceEscapeLimit_)
            {
                nSweeps++;

                if (addedFaces.found(endFace))
                {
                    completed = true;
                }

                // Add all face neighbours of face in the map
                const labelList cf = addedFaces.toc();
                addedFaces.clear();

                forAll (cf, cfI)
                {
                    const labelList& curNbrs = masterFaceFaces[cf[cfI]];

                    forAll (curNbrs, nbrI)
                    {
                        if (!curFaceMap.found(curNbrs[nbrI]))
                        {
                            curFaceMap.insert(curNbrs[nbrI]);
                            addedFaces.insert(curNbrs[nbrI]);
                        }
                    }
                }

                if (completed) break;

                if (debug)
                {
                    Pout << ".";
                }
            }

            if (!completed)
            {
                if (debug)
                {
                    Pout << "x";
                }

                // It is impossible to reach the end from the start, probably
                // due to disconnected domain.  Do search in opposite direction

                label nReverseSweeps = 0;

                addedFaces.clear();
                addedFaces.insert(endFace);

                while (nReverseSweeps < edgeFaceEscapeLimit_)
                {
                    nReverseSweeps++;

                    if (addedFaces.found(startFace))
                    {
                        completed = true;
                    }

                    // Add all face neighbours of face in the map
                    const labelList cf = addedFaces.toc();
                    addedFaces.clear();

                    forAll (cf, cfI)
                    {
                        const labelList& curNbrs = masterFaceFaces[cf[cfI]];

                        forAll (curNbrs, nbrI)
                        {
                            if (!curFaceMap.found(curNbrs[nbrI]))
                            {
                                curFaceMap.insert(curNbrs[nbrI]);
                                addedFaces.insert(curNbrs[nbrI]);
                            }
                        }
                    }

                    if (completed) break;

                    if (debug)
                    {
                        Pout << ".";
                    }
                }
            }

            if (completed)
            {
                if (debug)
                {
                    Pout << "+ ";
                }
            }
            else
            {
                if (debug)
                {
                    Pout << "z ";
                }
            }

            // Collect the edges

            // Create a map of edges the edge can interact with
            labelHashSet curMasterEdgesMap
            (
                nFacesPerSlaveEdge_*primitiveMesh::edgesPerFace_
            );

            const labelList curFaces = curFaceMap.toc();
//             Pout << "curFaces: " << curFaces << endl;
            forAll (curFaces, faceI)
            {
//                 Pout<< "face: " << curFaces[faceI] << " "
//                     << masterPatch[curFaces[faceI]]
//                     << " local: "
//                     << masterPatch.localFaces()[curFaces[faceI]]
//                     << endl;
                const labelList& me = masterFaceEdges[curFaces[faceI]];

                forAll (me, meI)
                {
                    curMasterEdgesMap.insert(me[meI]);
                }
            }

            const labelList curMasterEdges = curMasterEdgesMap.toc();

            // For all master edges to intersect, skip the ones
            // already used and cut the rest with a cutting plane.  If
            // the intersection point, falls inside of both edges, it
            // is valid.

            // Note.
            // The edge cutting code is repeated in
            // slidingInterface::modifyMotionPoints.  This is done for
            // efficiency reasons and avoids multiple creation of cutting
            // planes.  Please update both simultaneously.

            const point& a = projectedSlavePoints[curEdge.start()];
            const point& b = projectedSlavePoints[curEdge.end()];

            point c =
                0.5*
                (
                    slaveLocalPoints[curEdge.start()]
                  + slavePointNormals[curEdge.start()] // Add start normal
                  + slaveLocalPoints[curEdge.end()]
                  + slavePointNormals[curEdge.end()] // Add end normal
                );

            // Create the plane
            plane cutPlane(a, b, c);
//             Pout << "a: " << a << " b: " << b << " c: " << c << " plane: " << cutPlane << endl;

            linePointRef curSlaveLine = curEdge.line(projectedSlavePoints);
            const scalar curSlaveLineMag = curSlaveLine.mag();
//             Pout << "curSlaveLine: " << curSlaveLine << endl;
            forAll (curMasterEdges, masterEdgeI)
            {
                if (!curUme.found(curMasterEdges[masterEdgeI]))
                {
                    // New edge
                    if (debug)
                    {
                        Pout << "n";
                    }

                    const label cmeIndex = curMasterEdges[masterEdgeI];
                    const edge& cme = masterEdges[cmeIndex];
//                     Pout<< "Edge " << cmeIndex << " cme: " << cme << " line: " << cme.line(masterLocalPoints) << endl;
                    scalar cutOnMaster =
                        cutPlane.lineIntersect
                        (
                            cme.line(masterLocalPoints)
                        );

                    if
                    (
                        cutOnMaster > edgeEndCutoffTol_
                     && cutOnMaster < 1.0 - edgeEndCutoffTol_
                    )
                    {
                        // Master is cut, check the slave
                        point masterCutPoint =
                            masterLocalPoints[cme.start()]
                          + cutOnMaster*cme.vec(masterLocalPoints);

                        pointHit slaveCut =
                            curSlaveLine.nearestDist(masterCutPoint);

                        if (slaveCut.hit())
                        {
                            // Strict checking of slave cut to avoid capturing
                            // end points.
                            scalar cutOnSlave =
                                (
                                    (
                                        slaveCut.hitPoint()
                                      - curSlaveLine.start()
                                    ) & curSlaveLine.vec()
                                )/sqr(curSlaveLineMag);

                            // Calculate merge tolerance from the
                            // target edge length
                            scalar mergeTol =
                                edgeCoPlanarTol_*mag(b - a);
//                             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;
                            if
                            (
                                cutOnSlave > edgeEndCutoffTol_
                             && cutOnSlave < 1.0 - edgeEndCutoffTol_
                             && slaveCut.distance() < mergeTol
                            )
                            {
                                // Cut both master and slave.  Add point
                                // to edge points The point is nominally
                                // added from the start of the master edge
                                // and added to the cut point zone
                                label newPoint =
                                    ref.setAction
                                    (
                                        polyAddPoint
                                        (
                                            masterCutPoint,           // point
                                            masterMeshPoints[cme.start()],// m p
                                            cutPointZoneID_.index(),  // zone
                                            true                      // active
                                        )
                                    );
//                                 Pout << "Inserting point: " << newPoint << " as edge to edge intersection.  Slave edge: " << edgeI << " " << curEdge << " master edge: " << cmeIndex << " " << cme << endl;
                                pointsIntoSlaveEdges[edgeI].append(newPoint);
                                pointsIntoMasterEdges[cmeIndex].append(newPoint);

                                // Add the point into the enriched patch map
                                pointMap.insert
                                (
                                    newPoint,
                                    masterCutPoint
                                );

                                // Record which two edges intersect to
                                // create cut point
                                addToCpepm.insert
                                (
                                    newPoint,    // Cut point index
                                    Pair<edge>
                                    (
                                        edge
                                        (
                                            masterMeshPoints[cme.start()],
                                            masterMeshPoints[cme.end()]
                                        ),    // Master edge
                                        edge
                                        (
                                            slaveMeshPoints[curEdge.start()],
                                            slaveMeshPoints[curEdge.end()]
                                        )// Slave edge
                                    )
                                );

                                if (debug)
                                {
                                    Pout<< " " << newPoint << " = "
                                        << masterCutPoint << " ";
                                }
                            }
                            else
                            {
                                if (debug)
                                {
                                    // Intersection exists but it is too far
                                    Pout << "t";
                                }
                            }
                        }
                        else
                        {
                            if (debug)
                            {
                                // Missed slave edge
                                Pout << "x";
                            }
                        }
                    }
                    else
                    {
                        if (debug)
                        {
                            // Missed master edge
                            Pout << "-";
                        }
                    }
                }
                else
                {
                    if (debug)
                    {
                        Pout << "u";
                    }
                }
            }

            if (debug)
            {
                Pout << endl;
            }
        } // End if both ends missing
    } // End for all slave edges

//     Pout << "pointsIntoMasterEdges: " << pointsIntoMasterEdges << endl;
//     Pout << "pointsIntoSlaveEdges: " << pointsIntoSlaveEdges << endl;

    // Re-pack the points into edges lists
    labelListList pime(pointsIntoMasterEdges.size());

    forAll (pointsIntoMasterEdges, i)
    {
        pime[i].transfer(pointsIntoMasterEdges[i]);
    }

    labelListList pise(pointsIntoSlaveEdges.size());

    forAll (pointsIntoSlaveEdges, i)
    {
        pise[i].transfer(pointsIntoSlaveEdges[i]);
    }

    // Prepare the enriched faces
    cutPatch.calcEnrichedFaces
    (
        pime,
        pise,
        projectedSlavePoints
    );

    // Demand driven calculate the cut faces. Apart from the
    // cutFaces/cutFaceMaster/cutFaceSlave no information from the cutPatch
    // is used anymore!
    const faceList& cutFaces = cutPatch.cutFaces();
    const labelList& cutFaceMaster = cutPatch.cutFaceMaster();
    const labelList& cutFaceSlave = cutPatch.cutFaceSlave();

    const labelList& masterFc = masterFaceCells();
    const labelList& slaveFc = slaveFaceCells();

    // Couple the interface

    // Algorithm:
    // 1) Go through the cut faces and check if the cut face is the same as the
    //    defining master or slave.  If so, modify the appropriate
    //    face and mark the other for relegation into the face zone.
    // 2) If different, mark both sides for relegation and insert the new face


    boolList orphanedMaster(masterPatch.size(), false);
    boolList orphanedSlave(slavePatch.size(), false);

    forAll (cutFaces, faceI)
    {
        const face& curCutFace = cutFaces[faceI];
        const label curMaster = cutFaceMaster[faceI];
        const label curSlave = cutFaceSlave[faceI];

//         Pout << "Doing insertion of face " << faceI << ": ";

        // Check if the face has changed topologically
        bool insertedFace = false;

        if (curMaster >= 0)
        {
            // Face has got a master
            if (curCutFace == masterPatch[curMaster])
            {
                // Face is equal to master.  Modify master face.
//                 Pout << "Face is equal to master and is ";

                // If the face has got both master and slave, it is an
                // internal face; otherwise it is a patch face in the
                // master patch.  Keep it in the master face zone.

                if (curSlave >= 0)
                {
//                     Pout << "internal" << endl;
                    if (masterFc[curMaster] < slaveFc[curSlave])
                    {
                        // Cut face should point into slave.
                        // Be careful about flips in zone!
                        ref.setAction
                        (
                            polyModifyFace
                            (
                                curCutFace,                  // new face
                                masterPatchAddr[curMaster],  // master face id
                                masterFc[curMaster],         // owner
                                slaveFc[curSlave],           // neighbour
                                false,                       // flux flip
                                -1,                          // patch ID
                                false,                       // remove from zone
                                masterFaceZoneID_.index(),   // zone ID
                                masterPatchFlip[curMaster]   // zone flip
                            )
                        );
//                         Pout << "modifying master face. Old master: " << masterPatch[curMaster] << " new face: " << curCutFace.reverseFace() << " own: " << masterFc[curMaster] << " nei: " << slaveFc[curSlave] << endl;
                    }
                    else
                    {
                        // Cut face should point into master.  Flip required.
                        // Be careful about flips in zone!
                        ref.setAction
                        (
                            polyModifyFace
                            (
                                curCutFace.reverseFace(),    // new face
                                masterPatchAddr[curMaster],  // master face id
                                slaveFc[curSlave],           // owner
                                masterFc[curMaster],         // neighbour
                                true,                        // flux flip
                                -1,                          // patch ID
                                false,                       // remove from zone
                                masterFaceZoneID_.index(),   // zone ID
                                !masterPatchFlip[curMaster]  // zone flip
                            )
                        );
                    }

                    // Orphan the slave
                    orphanedSlave[curSlave] = true;
                }
                else
                {
//                     Pout << "master boundary" << endl;
                    ref.setAction
                    (
                        polyModifyFace
                        (
                            curCutFace,                  // new face
                            masterPatchAddr[curMaster],  // master face index
                            masterFc[curMaster],         // owner
                            -1,                          // neighbour
                            false,                       // flux flip
                            masterPatchID_.index(),      // patch ID
                            false,                       // remove from zone
                            masterFaceZoneID_.index(),   // zone ID
                            masterPatchFlip[curMaster]   // zone flip
                        )
                    );
                }

                insertedFace = true;
            }
        }
        else if (curSlave >= 0)
        {
            // Face has got a slave

            // Renumber the slave face into merged labels
            face rsf(slavePatch[curSlave]);

            forAll (rsf, i)
            {
                Map<label>::const_iterator mpIter = pointMergeMap.find(rsf[i]);

                if (mpIter != pointMergeMap.end())
                {
                    rsf[i] = mpIter();
                }
            }

            if (curCutFace == rsf)
            {
                // Face is equal to slave.  Modify slave face.
//                 Pout << "Face is equal to slave and is ";

                if (curMaster >= 0)
                {
//                     Pout << "regular internal" << endl;
                    if (masterFc[curMaster] < slaveFc[curSlave])
                    {
                        ref.setAction
                        (
                            polyModifyFace
                            (
                                curCutFace,                  // new face
                                slavePatchAddr[curSlave],    // master face id
                                masterFc[curMaster],         // owner
                                slaveFc[curSlave],           // neighbour
                                true,                        // flux flip
                                -1,                          // patch ID
                                false,                       // remove from zone
                                slaveFaceZoneID_.index(),    // zone ID
                                !slavePatchFlip[curMaster]   // zone flip
                            )
                        );
                    }
                    else
                    {
                        // Cut face should point into master.
                        // Be careful about flips in zone!
//                         Pout << "flipped internal" << endl;
                        ref.setAction
                        (
                            polyModifyFace
                            (
                                curCutFace.reverseFace(),    // new face
                                slavePatchAddr[curSlave],    // master face id
                                slaveFc[curSlave],           // owner
                                masterFc[curMaster],         // neighbour
                                false,                       // flux flip
                                -1,                          // patch ID
                                false,                       // remove from zone
                                slaveFaceZoneID_.index(),    // zone ID
                                slavePatchFlip[curSlave]     // zone flip
                            )
                        );
                    }

                    // Orphan the master
                    orphanedMaster[curMaster] = true;
                }
                else
                {
//                     Pout << "slave boundary" << endl;
                    ref.setAction
                    (
                        polyModifyFace
                        (
                            curCutFace,                  // new face
                            slavePatchAddr[curSlave],    // master face index
                            slaveFc[curSlave],           // owner
                            -1,                          // neighbour
                            false,                       // flux flip
                            slavePatchID_.index(),       // patch ID
                            false,                       // remove from zone
                            slaveFaceZoneID_.index(),    // zone ID
                            slavePatchFlip[curSlave]     // zone flip
                        )
                    );
                }

                insertedFace = true;
            }
        }
        else
        {
            FatalErrorIn
            (
                "void slidingInterface::coupleInterface("
                "polyTopoChange& ref) const"
            )   << "Face " << faceI << " in cut faces has neither a master "
                << "nor a slave.  Error in the cutting algorithm on modify."
                << abort(FatalError);
        }

        if (!insertedFace)
        {
            // Face is different from both master and slave
//             Pout << "Face different from both master and slave" << endl;

            if (curMaster >= 0)
            {
                if (curSlave >= 0)
                {
                    // Add internal face
                    if (masterFc[curMaster] < slaveFc[curSlave])
                    {
//                         Pout << "Adding internal face " << curCutFace << " owner: " << masterFc[curMaster] << " slave: " << slaveFc[curSlave] << " master face: " << masterPatchAddr[curMaster] << endl;
                        // Cut face should point into slave.
                        ref.setAction
                        (
                            polyAddFace
                            (
                                curCutFace,                  // new face
                                masterFc[curMaster],         // owner
                                slaveFc[curSlave],           // neighbour
                                -1,                          // master point
                                -1,                          // master edge
                                masterPatchAddr[curMaster],  // master face id
                                false,                       // flux flip
                                -1,                          // patch ID
                                cutFaceZoneID_.index(),      // zone ID
                                false                        // zone flip
                            )
                        );
                    }
                    else
                    {
                        // Cut face should point into master.  Flip required.
                        ref.setAction
                        (
                            polyAddFace
                            (
                                curCutFace.reverseFace(),    // new face
                                slaveFc[curSlave],           // owner
                                masterFc[curMaster],         // neighbour
                                -1,                          // master point
                                -1,                          // master edge
                                masterPatchAddr[curMaster],  // master face id
                                true,                        // flux flip
                                -1,                          // patch ID
                                cutFaceZoneID_.index(),      // zone ID
                                true                         // zone flip
                            )
                        );
                    }

                    // Orphan slave
                    orphanedSlave[curSlave] = true;
                }
                else
                {
//                 Pout << "Adding solo master face " << curCutFace << " owner: " << masterFc[curMaster] << " master face: " << masterPatchAddr[curMaster] << endl;
                    // Add master patch face
                    ref.setAction
                    (
                        polyAddFace
                        (
                            curCutFace,                  // new face
                            masterFc[curMaster],         // owner
                            -1,                          // neighbour
                            -1,                          // master point
                            -1,                          // master edge
                            masterPatchAddr[curMaster],  // master face index
                            false,                       // flux flip
                            masterPatchID_.index(),      // patch ID
                            cutFaceZoneID_.index(),      // zone ID
                            false                        // zone flip
                        )
                    );
                }

                // Orphan master
                orphanedMaster[curMaster] = true;
            }
            else if (curSlave >= 0)
            {
//                 Pout << "Adding solo slave face " << curCutFace << " owner: " << slaveFc[curSlave] << " master face: " << slavePatchAddr[curSlave] << endl;
                // Add slave patch face
                ref.setAction
                (
                    polyAddFace
                    (
                        curCutFace,                  // new face
                        slaveFc[curSlave],           // owner
                        -1,                          // neighbour
                        -1,                          // master point
                        -1,                          // master edge
                        slavePatchAddr[curSlave],    // master face index
                        false,                       // flux flip
                        slavePatchID_.index(),       // patch ID
                        cutFaceZoneID_.index(),      // zone ID
                        false                        // zone flip
                    )
                );

                // Orphan slave
                orphanedSlave[curSlave] = true;
            }
            else
            {
                FatalErrorIn
                (
                    "void slidingInterface::coupleInterface("
                    "polyTopoChange& ref) const"
                )   << "Face " << faceI << " in cut faces has neither a master "
                    << "nor a slave.  Error in the cutting algorithm on add."
                    << abort(FatalError);
            }
        }
    }

    // Move the orphaned faces into the face zone
//     Pout << "Orphaned master faces: " << orphanedMaster << endl;
//     Pout << "Orphaned slave faces: " << orphanedSlave << endl;

    label nOrphanedMasters = 0;

    forAll (orphanedMaster, faceI)
    {
        if (orphanedMaster[faceI])
        {
            nOrphanedMasters++;

            //ref.setAction
            //(
            //    polyModifyFace
            //    (
            //        masterPatch[faceI],                 // new face
            //        masterPatchAddr[faceI],             // master face index
            //        -1,                                 // owner
            //        -1,                                 // neighbour
            //        false,                              // flux flip
            //        -1,                                 // patch ID
            //        false,                              // remove from zone
            //        masterFaceZoneID_.index(),          // zone ID
            //        false                               // zone flip
            //    )
            //);

            //Pout<< "**MJ:deleting master face " << masterPatchAddr[faceI]
            //    << " old verts:" << masterPatch[faceI] << endl;
            ref.setAction(polyRemoveFace(masterPatchAddr[faceI]));
        }
    }

    label nOrphanedSlaves = 0;

    forAll (orphanedSlave, faceI)
    {
        if (orphanedSlave[faceI])
        {
            nOrphanedSlaves++;

            //ref.setAction
            //(
            //    polyModifyFace
            //    (
            //        slavePatch[faceI],                // new face
            //        slavePatchAddr[faceI],            // slave face index
            //        -1,                               // owner
            //        -1,                               // neighbour
            //        false,                            // flux flip
            //        -1,                               // patch ID
            //        false,                            // remove from zone
            //        slaveFaceZoneID_.index(),         // zone ID
            //        false                             // zone flip
            //    )
            //);

            //Pout<< "**MJ:deleting slave face " << slavePatchAddr[faceI]
            //    << " old verts:" << slavePatch[faceI] << endl;
            ref.setAction(polyRemoveFace(slavePatchAddr[faceI]));
        }
    }

    if (debug)
    {
        Pout<< "Number of orphaned faces: "
            << "master = " << nOrphanedMasters << " out of "
            << orphanedMaster.size()
            << " slave = " << nOrphanedSlaves << " out of "
            << orphanedSlave.size() << endl;
    }

    // Finished coupling the plane of sliding interface.

    // Modify faces influenced by the sliding interface
    // These are the faces belonging to the master and slave cell
    // layer which have not already been modified.
    // Modification comes in three types:
    // 1) eliminate the points that have been removed when the old sliding
    //    interface has been removed
    // 2) Merge the slave points that have hit points on the master patch
    // 3) Introduce new points resulting from the new sliding interface cut

    // Collect all affected faces

    // Master side

    // Grab the list of faces in the layer
    const labelList& masterStickOuts = masterStickOutFaces();

//     Pout << "masterStickOuts: " << masterStickOuts << endl;

    // Re-create the master stick-out faces
    forAll (masterStickOuts, faceI)
    {
        // Renumber the face and remove additional points

        const label curFaceID = masterStickOuts[faceI];

        const face& oldRichFace = faces[curFaceID];

        bool changed = false;

        // Remove removed points from the face
        face oldFace(oldRichFace.size());
        label nOldFace = 0;

        forAll (oldRichFace, pointI)
        {
            if (ref.pointRemoved(oldRichFace[pointI]))
            {
                changed = true;
            }
            else
            {
                // Point off patch
                oldFace[nOldFace] = oldRichFace[pointI];
                nOldFace++;
            }
        }

        oldFace.setSize(nOldFace);

//         Pout << "old rich master face: " << oldRichFace << " old face: " << oldFace << endl;
        DynamicList<label> newFaceLabels(2*oldFace.size());

        forAll (oldFace, pointI)
        {
            if (masterMeshPointMap.found(oldFace[pointI]))
            {
                // Point is in master patch. Add it

                // If the point is a direct hit, grab its label; otherwise
                // keep the original
                if (pointMergeMap.found(oldFace[pointI]))
                {
                    changed = true;
                    newFaceLabels.append
                    (
                        pointMergeMap.find(oldFace[pointI])()
                    );
                }
                else
                {
                    newFaceLabels.append(oldFace[pointI]);
                }

                // Find if there are additional points inserted onto the edge
                // between current point and the next point
                // Algorithm:
                // 1) Find all the edges in the master patch coming
                //    out of the current point.
                // 2) If the next point in the face to pick the right edge
                const label localFirstLabel =
                    masterMeshPointMap.find(oldFace[pointI])();

                const labelList& curEdges = masterPointEdges[localFirstLabel];

                const label  nextLabel = oldFace.nextLabel(pointI);

                Map<label>::const_iterator mmpmIter =
                    masterMeshPointMap.find(nextLabel);

                if (mmpmIter != masterMeshPointMap.end())
                {
//                     Pout << "found label pair " << oldFace[pointI] << " and " << nextLabel;
                    // Find the points on the edge between them
                    const label localNextLabel = mmpmIter();

                    forAll (curEdges, curEdgeI)
                    {
                        if
                        (
                            masterEdges[curEdges[curEdgeI]].otherVertex
                            (
                                localFirstLabel
                            )
                         == localNextLabel
                        )
                        {
//                             Pout << " found edge: " << curEdges[curEdgeI] << endl;

                            // Get points on current edge
                            const labelList& curPime = pime[curEdges[curEdgeI]];

                            if (curPime.size())
                            {
                                changed = true;
                                // Pout << "curPime: " << curPime << endl;
                                // Insert the edge points into the face
                                // in the correct order
                                const point& startPoint =
                                    masterLocalPoints[localFirstLabel];

                                vector e =
                                    masterLocalPoints[localNextLabel]
                                  - startPoint;

                                e /= magSqr(e);

                                scalarField edgePointWeights(curPime.size());

                                forAll (curPime, curPimeI)
                                {
                                    edgePointWeights[curPimeI] =
                                        (
                                            e
                                          & (
                                              pointMap.find(curPime[curPimeI])()
                                            - startPoint
                                            )
                                        );
                                }

                                if (debug)
                                {
                                    if
                                    (
                                        min(edgePointWeights) < 0
                                     || max(edgePointWeights) > 1
                                    )
                                    {
                                        FatalErrorIn
                                        (
                                            "void slidingInterface::"
                                            "coupleInterface("
                                            "polyTopoChange& ref) const"
                                        )   << "Error in master stick-out edge "
                                            << "point collection."
                                            << abort(FatalError);
                                    }
                                }

                                // Go through the points and collect
                                // them based on weights from lower to
                                // higher.  This gives the correct
                                // order of points along the edge.
                                for
                                (
                                    label passI = 0;
                                    passI < edgePointWeights.size();
                                    passI++
                                )
                                {
                                    // Max weight can only be one, so
                                    // the sorting is done by
                                    // elimination.
                                    label nextPoint = -1;
                                    scalar dist = 2;

                                    forAll (edgePointWeights, wI)
                                    {
                                        if (edgePointWeights[wI] < dist)
                                        {
                                            dist = edgePointWeights[wI];
                                            nextPoint = wI;
                                        }
                                    }

                                    // Insert the next point and reset
                                    // its weight to exclude it from
                                    // future picks
                                    newFaceLabels.append(curPime[nextPoint]);
                                    edgePointWeights[nextPoint] = GREAT;
                                }
                            }

                            break;
                        } // End of found edge
                    } // End of looking through current edges
                }
            }
            else
            {
                newFaceLabels.append(oldFace[pointI]);
            }
        }

        if (changed)
        {
            if (newFaceLabels.size() < 3)
            {
                FatalErrorIn
                (
                    "void slidingInterface::coupleInterface("
                    "polyTopoChange& ref) const"
                )   << "Face " << curFaceID << " reduced to less than "
                    << "3 points.  Topological/cutting error A." << nl
                    << "Old face: " << oldFace << " new face: " << newFaceLabels
                    << abort(FatalError);
            }

            // Get face zone and its flip
            label modifiedFaceZone = faceZones.whichZone(curFaceID);
            bool modifiedFaceZoneFlip = false;

            if (modifiedFaceZone >= 0)
            {
                modifiedFaceZoneFlip =
                    faceZones[modifiedFaceZone].flipMap()
                    [
                        faceZones[modifiedFaceZone].whichFace(curFaceID)
                    ];
            }

            face newFace;
            newFace.transfer(newFaceLabels);

            //Pout << "Modifying master stick-out face " << curFaceID
            //    << " old face: " << oldFace << " new face: " << newFace << endl;

            // Modify the face
            if (mesh.isInternalFace(curFaceID))
            {
                ref.setAction
                (
                    polyModifyFace
                    (
                        newFace,                // modified face
                        curFaceID,              // label of face being modified
                        own[curFaceID],         // owner
                        nei[curFaceID],         // neighbour
                        false,                  // face flip
                        -1,                     // patch for face
                        false,                  // remove from zone
                        modifiedFaceZone,       // zone for face
                        modifiedFaceZoneFlip    // face flip in zone
                    )
                );
            }
            else
            {
                ref.setAction
                (
                    polyModifyFace
                    (
                        newFace,                // modified face
                        curFaceID,              // label of face being modified
                        own[curFaceID],         // owner
                        -1,                     // neighbour
                        false,                  // face flip
                        mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
                        false,                  // remove from zone
                        modifiedFaceZone,       // zone for face
                        modifiedFaceZoneFlip    // face flip in zone
                    )
                );
            }
        }
    }

//     Pout << "Finished master side" << endl;

    // Slave side

    // Grab the list of faces in the layer
    const labelList& slaveStickOuts = slaveStickOutFaces();

//     Pout << "slaveStickOuts: " << slaveStickOuts << endl;

    const Map<label>& rpm = retiredPointMap();

    // Re-create the slave stick-out faces

    forAll (slaveStickOuts, faceI)
    {
        // Renumber the face and remove additional points
        const label curFaceID = slaveStickOuts[faceI];

        const face& oldRichFace = faces[curFaceID];

        bool changed = false;

        // Remove removed points from the face
        face oldFace(oldRichFace.size());
        label nOldFace = 0;

        forAll (oldRichFace, pointI)
        {
            if
            (
                rpm.found(oldRichFace[pointI])
             || slaveMeshPointMap.found(oldRichFace[pointI])
            )
            {
                // Point definitely live. Add it
                oldFace[nOldFace] = oldRichFace[pointI];
                nOldFace++;
            }
            else if
            (
                ref.pointRemoved(oldRichFace[pointI])
             || masterMeshPointMap.found(oldRichFace[pointI])
            )
            {
                // Point removed and not on slave patch
                // (first if takes care of that!) or
                // point belonging to master patch
                changed = true;
            }
            else
            {
                // Point off patch
                oldFace[nOldFace] = oldRichFace[pointI];
                nOldFace++;
            }
        }

        oldFace.setSize(nOldFace);

        DynamicList<label> newFaceLabels(2*oldFace.size());

//         Pout << "old rich slave face: " << oldRichFace << " old face: " << oldFace << endl;
        forAll (oldFace, pointI)
        {
            // Try to find the point in retired points
            label curP = oldFace[pointI];

            Map<label>::const_iterator rpmIter = rpm.find(oldFace[pointI]);

            if (rpmIter != rpm.end())
            {
                changed = true;
                curP = rpmIter();
            }

            if (slaveMeshPointMap.found(curP))
            {
                // Point is in slave patch. Add it

                // If the point is a direct hit, grab its label; otherwise
                // keep the original
                if (pointMergeMap.found(curP))
                {
                    changed = true;
                    newFaceLabels.append
                    (
                        pointMergeMap.find(curP)()
                    );
                }
                else
                {
                    newFaceLabels.append(curP);
                }

                // Find if there are additional points inserted onto the edge
                // between current point and the next point
                // Algorithm:
                // 1) Find all the edges in the slave patch coming
                //    out of the current point.
                // 2) Use the next point in the face to pick the right edge

                const label localFirstLabel =
                    slaveMeshPointMap.find(curP)();

                const labelList& curEdges = slavePointEdges[localFirstLabel];

                label nextLabel = oldFace.nextLabel(pointI);

                Map<label>::const_iterator rpmNextIter =
                    rpm.find(nextLabel);

                if (rpmNextIter != rpm.end())
                {
                    nextLabel = rpmNextIter();
                }

                Map<label>::const_iterator mmpmIter =
                    slaveMeshPointMap.find(nextLabel);

                if (mmpmIter != slaveMeshPointMap.end())
                {
                    // Both points on the slave patch.
                    // Find the points on the edge between them
                    const label localNextLabel = mmpmIter();

                    forAll (curEdges, curEdgeI)
                    {
                        if
                        (
                            slaveEdges[curEdges[curEdgeI]].otherVertex
                            (
                                localFirstLabel
                            )
                         == localNextLabel
                        )
                        {
//                             Pout << " found edge: " << curEdges[curEdgeI] << endl;

                            // Get points on current edge
                            const labelList& curPise = pise[curEdges[curEdgeI]];

                            if (curPise.size())
                            {
                                changed = true;
//                                 Pout << "curPise: " << curPise << endl;
                                // Insert the edge points into the face
                                // in the correct order
                                const point& startPoint =
                                    projectedSlavePoints[localFirstLabel];

                                vector e =
                                    projectedSlavePoints[localNextLabel]
                                  - startPoint;

                                e /= magSqr(e);

                                scalarField edgePointWeights(curPise.size());

                                forAll (curPise, curPiseI)
                                {
                                    edgePointWeights[curPiseI] =
                                    (
                                        e
                                      & (
                                            pointMap.find(curPise[curPiseI])()
                                          - startPoint
                                        )
                                    );
                                }

                                if (debug)
                                {
                                    if
                                    (
                                        min(edgePointWeights) < 0
                                     || max(edgePointWeights) > 1
                                    )
                                    {
                                        FatalErrorIn
                                        (
                                            "void slidingInterface::"
                                            "coupleInterface("
                                            "polyTopoChange& ref) const"
                                        )   << "Error in slave stick-out edge "
                                            << "point collection."
                                            << abort(FatalError);
                                        }
                                    }

                                // Go through the points and collect
                                // them based on weights from lower to
                                // higher.  This gives the correct
                                // order of points along the edge.
                                for
                                (
                                    label passI = 0;
                                    passI < edgePointWeights.size();
                                    passI++
                                )
                                {
                                    // Max weight can only be one, so
                                    // the sorting is done by
                                    // elimination.
                                    label nextPoint = -1;
                                    scalar dist = 2;

                                    forAll (edgePointWeights, wI)
                                    {
                                        if (edgePointWeights[wI] < dist)
                                        {
                                            dist = edgePointWeights[wI];
                                            nextPoint = wI;
                                        }
                                    }

                                    // Insert the next point and reset
                                    // its weight to exclude it from
                                    // future picks
                                    newFaceLabels.append(curPise[nextPoint]);
                                    edgePointWeights[nextPoint] = GREAT;
                                }
                            }

                            break;
                        }
                    } // End of found edge
                } // End of looking through current edges
            }
            else
            {
                newFaceLabels.append(oldFace[pointI]);
            }
        }

        if (changed)
        {
            if (newFaceLabels.size() < 3)
            {
                FatalErrorIn
                (
                    "void slidingInterface::coupleInterface("
                    "polyTopoChange& ref) const"
                )   << "Face " << curFaceID << " reduced to less than "
                    << "3 points.  Topological/cutting error B." << nl
                    << "Old face: " << oldFace << " new face: " << newFaceLabels
                    << abort(FatalError);
            }

            // Get face zone and its flip
            label modifiedFaceZone = faceZones.whichZone(curFaceID);
            bool modifiedFaceZoneFlip = false;

            if (modifiedFaceZone >= 0)
            {
                modifiedFaceZoneFlip =
                    faceZones[modifiedFaceZone].flipMap()
                    [
                        faceZones[modifiedFaceZone].whichFace(curFaceID)
                    ];
            }

            face newFace;
            newFace.transfer(newFaceLabels);

//             Pout << "Modifying slave stick-out face " << curFaceID << " old face: " << oldFace << " new face: " << newFace << endl;

            // Modify the face
            if (mesh.isInternalFace(curFaceID))
            {
                ref.setAction
                (
                    polyModifyFace
                    (
                        newFace,                // modified face
                        curFaceID,              // label of face being modified
                        own[curFaceID],         // owner
                        nei[curFaceID],         // neighbour
                        false,                  // face flip
                        -1,                     // patch for face
                        false,                  // remove from zone
                        modifiedFaceZone,       // zone for face
                        modifiedFaceZoneFlip    // face flip in zone
                    )
                );
            }
            else
            {
                ref.setAction
                (
                    polyModifyFace
                    (
                        newFace,                // modified face
                        curFaceID,              // label of face being modified
                        own[curFaceID],         // owner
                        -1,                     // neighbour
                        false,                  // face flip
                        mesh.boundaryMesh().whichPatch(curFaceID), // patch for face
                        false,                  // remove from zone
                        modifiedFaceZone,       // zone for face
                        modifiedFaceZoneFlip    // face flip in zone
                    )
                );
            }
        }
    }

    // Activate and retire slave patch points
    // This needs to be done last, so that the map of removed points
    // does not get damaged by point modifications

    if (!retiredPointMapPtr_)
    {
        FatalErrorIn
        (
            "void slidingInterface::coupleInterface("
            "polyTopoChange& ref) const"
        )   << "Retired point map pointer not set."
            << abort(FatalError);
    }

    Map<label>& addToRpm = *retiredPointMapPtr_;

    // Clear the old map
    addToRpm.clear();

    label nRetiredPoints = 0;

    forAll (slaveMeshPoints, pointI)
    {
        if (pointMergeMap.found(slaveMeshPoints[pointI]))
        {
            // Retire the point - only used for supporting the detached
            // slave patch
            nRetiredPoints++;

            //ref.setAction
            //(
            //    polyModifyPoint
            //    (
            //        slaveMeshPoints[pointI],             // point ID
            //        points[slaveMeshPoints[pointI]],     // point
            //        false,                               // remove from zone
            //        mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
            //        false                                // in a cell
            //    )
            //);
            //Pout<< "MJ retire slave point " << slaveMeshPoints[pointI]
            //    << " coord " << points[slaveMeshPoints[pointI]]
            //    << endl;
            ref.setAction
            (
                polyRemovePoint
                (
                    slaveMeshPoints[pointI]
                )
            );

            addToRpm.insert
            (
                pointMergeMap.find(slaveMeshPoints[pointI])(),
                slaveMeshPoints[pointI]
            );
        }
        else
        {
            ref.setAction
            (
                polyModifyPoint
                (
                    slaveMeshPoints[pointI],             // point ID
                    points[slaveMeshPoints[pointI]],     // point
                    false,                               // remove from zone
                    mesh.pointZones().whichZone(slaveMeshPoints[pointI]),// zone
                    true                                 // in a cell
                )
            );
        }
    }

    if (debug)
    {
        Pout << "Retired " << nRetiredPoints << " out of "
            << slaveMeshPoints.size() << " points." << endl;
    }

    // Grab cut face master and slave addressing
    if (cutFaceMasterPtr_) deleteDemandDrivenData(cutFaceMasterPtr_);
    cutFaceMasterPtr_ = new labelList(cutPatch.cutFaceMaster());

    if (cutFaceSlavePtr_) deleteDemandDrivenData(cutFaceSlavePtr_);
    cutFaceSlavePtr_ = new labelList(cutPatch.cutFaceSlave());

    // Finished coupling
    attached_ = true;

    if (debug)
    {
        Pout<< "void slidingInterface::coupleInterface("
            << "polyTopoChange& ref) : "
            << "Finished coupling sliding interface " << name() << endl;
    }
}


// ************************ vim: set sw=4 sts=4 et: ************************ //