distributedTriSurfaceMesh.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "distributedTriSurfaceMesh.H"
#include "mapDistribute.H"
#include "Random.H"
#include "addToRunTimeSelectionTable.H"
#include "triangleFuncs.H"
#include "matchPoints.H"
#include "globalIndex.H"
#include "Time.H"

#include "IFstream.H"
#include "decompositionMethod.H"
#include "geomDecomp.H"
#include "vectorList.H"
#include "PackedBoolList.H"
#include "PatchTools.H"
#include "OBJstream.H"

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

namespace Foam
{
    defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
    addToRunTimeSelectionTable
    (
        searchableSurface,
        distributedTriSurfaceMesh,
        dict
    );

    template<>
    const char* Foam::NamedEnum
    <
        Foam::distributedTriSurfaceMesh::distributionType,
        3
    >::names[] =
    {
        "follow",
        "independent",
        "frozen"
    };
}


const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
    Foam::distributedTriSurfaceMesh::distributionTypeNames_;


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

// Read my additional data from the dictionary
bool Foam::distributedTriSurfaceMesh::read()
{
    // Get bb of all domains.
    procBb_.setSize(Pstream::nProcs());

    procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
    Pstream::gatherList(procBb_);
    Pstream::scatterList(procBb_);

    // Distribution type
    distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));

    // Merge distance
    mergeDist_ = readScalar(dict_.lookup("mergeDistance"));

    return true;
}


// Is segment fully local?
bool Foam::distributedTriSurfaceMesh::isLocal
(
    const List<treeBoundBox>& myBbs,
    const point& start,
    const point& end
)
{
    forAll(myBbs, bbI)
    {
        if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
        {
            return true;
        }
    }
    return false;
}


//void Foam::distributedTriSurfaceMesh::splitSegment
//(
//    const label segmentI,
//    const point& start,
//    const point& end,
//    const treeBoundBox& bb,
//
//    DynamicList<segment>& allSegments,
//    DynamicList<label>& allSegmentMap,
//    DynamicList<label> sendMap
//) const
//{
//    // Work points
//    point clipPt0, clipPt1;
//
//    if (bb.contains(start))
//    {
//        // start within, trim end to bb
//        bool clipped = bb.intersects(end, start, clipPt0);
//
//        if (clipped)
//        {
//            // segment from start to clippedStart passes
//            // through proc.
//            sendMap[procI].append(allSegments.size());
//            allSegmentMap.append(segmentI);
//            allSegments.append(segment(start, clipPt0));
//        }
//    }
//    else if (bb.contains(end))
//    {
//        // end within, trim start to bb
//        bool clipped = bb.intersects(start, end, clipPt0);
//
//        if (clipped)
//        {
//            sendMap[procI].append(allSegments.size());
//            allSegmentMap.append(segmentI);
//            allSegments.append(segment(clipPt0, end));
//        }
//    }
//    else
//    {
//        // trim both
//        bool clippedStart = bb.intersects(start, end, clipPt0);
//
//        if (clippedStart)
//        {
//            bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
//
//            if (clippedEnd)
//            {
//                // middle part of segment passes through proc.
//                sendMap[procI].append(allSegments.size());
//                allSegmentMap.append(segmentI);
//                allSegments.append(segment(clipPt0, clipPt1));
//            }
//        }
//    }
//}


void Foam::distributedTriSurfaceMesh::distributeSegment
(
    const label segmentI,
    const point& start,
    const point& end,

    DynamicList<segment>& allSegments,
    DynamicList<label>& allSegmentMap,
    List<DynamicList<label>>& sendMap
) const
{
    // 1. Fully local already handled outside. Note: retest is cheap.
    if (isLocal(procBb_[Pstream::myProcNo()], start, end))
    {
        return;
    }


    // 2. If fully inside one other processor, then only need to send
    // to that one processor even if it intersects another. Rare occurrence
    // but cheap to test.
    forAll(procBb_, procI)
    {
        if (procI != Pstream::myProcNo())
        {
            const List<treeBoundBox>& bbs = procBb_[procI];

            if (isLocal(bbs, start, end))
            {
                sendMap[procI].append(allSegments.size());
                allSegmentMap.append(segmentI);
                allSegments.append(segment(start, end));
                return;
            }
        }
    }


    // 3. If not contained in single processor send to all intersecting
    // processors.
    forAll(procBb_, procI)
    {
        const List<treeBoundBox>& bbs = procBb_[procI];

        forAll(bbs, bbI)
        {
            const treeBoundBox& bb = bbs[bbI];

            // Scheme a: any processor that intersects the segment gets
            // the segment.

            // Intersection point
            point clipPt;

            if (bb.intersects(start, end, clipPt))
            {
                sendMap[procI].append(allSegments.size());
                allSegmentMap.append(segmentI);
                allSegments.append(segment(start, end));
            }

            // Alternative: any processor only gets clipped bit of
            // segment. This gives small problems with additional
            // truncation errors.
            //splitSegment
            //(
            //    segmentI,
            //    start,
            //    end,
            //    bb,
            //
            //    allSegments,
            //    allSegmentMap,
            //   sendMap[procI]
            //);
        }
    }
}


Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::distributeSegments
(
    const pointField& start,
    const pointField& end,

    List<segment>& allSegments,
    labelList& allSegmentMap
) const
{
    // Determine send map
    // ~~~~~~~~~~~~~~~~~~

    labelListList sendMap(Pstream::nProcs());

    {
        // Since intersection test is quite expensive compared to memory
        // allocation we use DynamicList to immediately store the segment
        // in the correct bin.

        // Segments to test
        DynamicList<segment> dynAllSegments(start.size());
        // Original index of segment
        DynamicList<label> dynAllSegmentMap(start.size());
        // Per processor indices into allSegments to send
        List<DynamicList<label>> dynSendMap(Pstream::nProcs());

        forAll(start, segmentI)
        {
            distributeSegment
            (
                segmentI,
                start[segmentI],
                end[segmentI],

                dynAllSegments,
                dynAllSegmentMap,
                dynSendMap
            );
        }

        // Convert dynamicList to labelList
        sendMap.setSize(Pstream::nProcs());
        forAll(sendMap, procI)
        {
            dynSendMap[procI].shrink();
            sendMap[procI].transfer(dynSendMap[procI]);
        }

        allSegments.transfer(dynAllSegments.shrink());
        allSegmentMap.transfer(dynAllSegmentMap.shrink());
    }


    // Send over how many I need to receive.
    labelListList sendSizes(Pstream::nProcs());
    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
    forAll(sendMap, procI)
    {
        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
    }
    Pstream::gatherList(sendSizes);
    Pstream::scatterList(sendSizes);


    // Determine order of receiving
    labelListList constructMap(Pstream::nProcs());

    // My local segments first
    constructMap[Pstream::myProcNo()] = identity
    (
        sendMap[Pstream::myProcNo()].size()
    );

    label segmentI = constructMap[Pstream::myProcNo()].size();
    forAll(constructMap, procI)
    {
        if (procI != Pstream::myProcNo())
        {
            // What I need to receive is what other processor is sending to me.
            label nRecv = sendSizes[procI][Pstream::myProcNo()];
            constructMap[procI].setSize(nRecv);

            for (label i = 0; i < nRecv; i++)
            {
                constructMap[procI][i] = segmentI++;
            }
        }
    }

    return autoPtr<mapDistribute>
    (
        new mapDistribute
        (
            segmentI,       // size after construction
            sendMap.xfer(),
            constructMap.xfer()
        )
    );
}


void Foam::distributedTriSurfaceMesh::findLine
(
    const bool nearestIntersection,
    const pointField& start,
    const pointField& end,
    List<pointIndexHit>& info
) const
{
    const indexedOctree<treeDataTriSurface>& octree = tree();

    // Initialise
    info.setSize(start.size());
    forAll(info, i)
    {
        info[i].setMiss();
    }

    if (!Pstream::parRun())
    {
        forAll(start, i)
        {
            if (nearestIntersection)
            {
                info[i] = octree.findLine(start[i], end[i]);
            }
            else
            {
                info[i] = octree.findLineAny(start[i], end[i]);
            }
        }
    }
    else
    {
        // Important:force synchronised construction of indexing
        const globalIndex& triIndexer = globalTris();


        // Do any local queries
        // ~~~~~~~~~~~~~~~~~~~~

        label nLocal = 0;

        forAll(start, i)
        {
            if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
            {
                if (nearestIntersection)
                {
                    info[i] = octree.findLine(start[i], end[i]);
                }
                else
                {
                    info[i] = octree.findLineAny(start[i], end[i]);
                }

                if (info[i].hit())
                {
                    info[i].setIndex(triIndexer.toGlobal(info[i].index()));
                }
                nLocal++;
            }
        }


        if
        (
            returnReduce(nLocal, sumOp<label>())
          < returnReduce(start.size(), sumOp<label>())
        )
        {
            // Not all can be resolved locally. Build segments and map,
            // send over segments, do intersections, send back and merge.


            // Construct queries (segments)
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

            // Segments to test
            List<segment> allSegments(start.size());
            // Original index of segment
            labelList allSegmentMap(start.size());

            const autoPtr<mapDistribute> mapPtr
            (
                distributeSegments
                (
                    start,
                    end,
                    allSegments,
                    allSegmentMap
                )
            );
            const mapDistribute& map = mapPtr();

            label nOldAllSegments = allSegments.size();


            // Exchange the segments
            // ~~~~~~~~~~~~~~~~~~~~~

            map.distribute(allSegments);


            // Do tests I need to do
            // ~~~~~~~~~~~~~~~~~~~~~

            // Intersections
            List<pointIndexHit> intersections(allSegments.size());

            forAll(allSegments, i)
            {
                if (nearestIntersection)
                {
                    intersections[i] = octree.findLine
                    (
                        allSegments[i].first(),
                        allSegments[i].second()
                    );
                }
                else
                {
                    intersections[i] = octree.findLineAny
                    (
                        allSegments[i].first(),
                        allSegments[i].second()
                    );
                }

                // Convert triangle index to global numbering
                if (intersections[i].hit())
                {
                    intersections[i].setIndex
                    (
                        triIndexer.toGlobal(intersections[i].index())
                    );
                }
            }


            // Exchange the intersections (opposite to segments)
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

            map.reverseDistribute(nOldAllSegments, intersections);


            // Extract the hits
            // ~~~~~~~~~~~~~~~~

            forAll(intersections, i)
            {
                const pointIndexHit& allInfo = intersections[i];
                label segmentI = allSegmentMap[i];
                pointIndexHit& hitInfo = info[segmentI];

                if (allInfo.hit())
                {
                    if (!hitInfo.hit())
                    {
                        // No intersection yet so take this one
                        hitInfo = allInfo;
                    }
                    else if (nearestIntersection)
                    {
                        // Nearest intersection
                        if
                        (
                            magSqr(allInfo.hitPoint()-start[segmentI])
                          < magSqr(hitInfo.hitPoint()-start[segmentI])
                        )
                        {
                            hitInfo = allInfo;
                        }
                    }
                }
            }
        }
    }
}


// Exchanges indices to the processor they come from.
// - calculates exchange map
// - uses map to calculate local triangle index
Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
    const List<pointIndexHit>& info,
    labelList& triangleIndex
) const
{
    triangleIndex.setSize(info.size());

    const globalIndex& triIndexer = globalTris();


    // Determine send map
    // ~~~~~~~~~~~~~~~~~~

    // Since determining which processor the query should go to is
    // cheap we do a multi-pass algorithm to save some memory temporarily.

    // 1. Count
    labelList nSend(Pstream::nProcs(), 0);

    forAll(info, i)
    {
        if (info[i].hit())
        {
            label procI = triIndexer.whichProcID(info[i].index());
            nSend[procI]++;
        }
    }

    // 2. Size sendMap
    labelListList sendMap(Pstream::nProcs());
    forAll(nSend, procI)
    {
        sendMap[procI].setSize(nSend[procI]);
        nSend[procI] = 0;
    }

    // 3. Fill sendMap
    forAll(info, i)
    {
        if (info[i].hit())
        {
            label procI = triIndexer.whichProcID(info[i].index());
            triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
            sendMap[procI][nSend[procI]++] = i;
        }
        else
        {
            triangleIndex[i] = -1;
        }
    }


    // Send over how many I need to receive
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    labelListList sendSizes(Pstream::nProcs());
    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
    forAll(sendMap, procI)
    {
        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
    }
    Pstream::gatherList(sendSizes);
    Pstream::scatterList(sendSizes);


    // Determine receive map
    // ~~~~~~~~~~~~~~~~~~~~~

    labelListList constructMap(Pstream::nProcs());

    // My local segments first
    constructMap[Pstream::myProcNo()] = identity
    (
        sendMap[Pstream::myProcNo()].size()
    );

    label segmentI = constructMap[Pstream::myProcNo()].size();
    forAll(constructMap, procI)
    {
        if (procI != Pstream::myProcNo())
        {
            // What I need to receive is what other processor is sending to me.
            label nRecv = sendSizes[procI][Pstream::myProcNo()];
            constructMap[procI].setSize(nRecv);

            for (label i = 0; i < nRecv; i++)
            {
                constructMap[procI][i] = segmentI++;
            }
        }
    }


    // Pack into distribution map
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~

    autoPtr<mapDistribute> mapPtr
    (
        new mapDistribute
        (
            segmentI,       // size after construction
            sendMap.xfer(),
            constructMap.xfer()
        )
    );
    const mapDistribute& map = mapPtr();


    // Send over queries
    // ~~~~~~~~~~~~~~~~~

    map.distribute(triangleIndex);


    return mapPtr;
}


Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
(
    const point& centre,
    const scalar radiusSqr,
    boolList& overlaps
) const
{
    overlaps = false;
    label nOverlaps = 0;

    forAll(procBb_, procI)
    {
        const List<treeBoundBox>& bbs = procBb_[procI];

        forAll(bbs, bbI)
        {
            if (bbs[bbI].overlaps(centre, radiusSqr))
            {
                overlaps[procI] = true;
                nOverlaps++;
                break;
            }
        }
    }
    return nOverlaps;
}


// Generate queries for parallel distance calculation
// - calculates exchange map
// - uses map to exchange points and radius
Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
    const pointField& centres,
    const scalarField& radiusSqr,

    pointField& allCentres,
    scalarField& allRadiusSqr,
    labelList& allSegmentMap
) const
{
    // Determine queries
    // ~~~~~~~~~~~~~~~~~

    labelListList sendMap(Pstream::nProcs());

    {
        // Queries
        DynamicList<point> dynAllCentres(centres.size());
        DynamicList<scalar> dynAllRadiusSqr(centres.size());
        // Original index of segment
        DynamicList<label> dynAllSegmentMap(centres.size());
        // Per processor indices into allSegments to send
        List<DynamicList<label>> dynSendMap(Pstream::nProcs());

        // Work array - whether processor bb overlaps the bounding sphere.
        boolList procBbOverlaps(Pstream::nProcs());

        forAll(centres, centreI)
        {
            // Find the processor this sample+radius overlaps.
            calcOverlappingProcs
            (
                centres[centreI],
                radiusSqr[centreI],
                procBbOverlaps
            );

            forAll(procBbOverlaps, procI)
            {
                if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
                {
                    dynSendMap[procI].append(dynAllCentres.size());
                    dynAllSegmentMap.append(centreI);
                    dynAllCentres.append(centres[centreI]);
                    dynAllRadiusSqr.append(radiusSqr[centreI]);
                }
            }
        }

        // Convert dynamicList to labelList
        sendMap.setSize(Pstream::nProcs());
        forAll(sendMap, procI)
        {
            dynSendMap[procI].shrink();
            sendMap[procI].transfer(dynSendMap[procI]);
        }

        allCentres.transfer(dynAllCentres.shrink());
        allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
        allSegmentMap.transfer(dynAllSegmentMap.shrink());
    }


    // Send over how many I need to receive.
    labelListList sendSizes(Pstream::nProcs());
    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
    forAll(sendMap, procI)
    {
        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
    }
    Pstream::gatherList(sendSizes);
    Pstream::scatterList(sendSizes);


    // Determine order of receiving
    labelListList constructMap(Pstream::nProcs());

    // My local segments first
    constructMap[Pstream::myProcNo()] = identity
    (
        sendMap[Pstream::myProcNo()].size()
    );

    label segmentI = constructMap[Pstream::myProcNo()].size();
    forAll(constructMap, procI)
    {
        if (procI != Pstream::myProcNo())
        {
            // What I need to receive is what other processor is sending to me.
            label nRecv = sendSizes[procI][Pstream::myProcNo()];
            constructMap[procI].setSize(nRecv);

            for (label i = 0; i < nRecv; i++)
            {
                constructMap[procI][i] = segmentI++;
            }
        }
    }

    autoPtr<mapDistribute> mapPtr
    (
        new mapDistribute
        (
            segmentI,       // size after construction
            sendMap.xfer(),
            constructMap.xfer()
        )
    );
    return mapPtr;
}


// Find bounding boxes that guarantee a more or less uniform distribution
// of triangles. Decomposition in here is only used to get the bounding
// boxes, actual decomposition is done later on.
// Returns a per processor a list of bounding boxes that most accurately
// describe the shape. For now just a single bounding box per processor but
// optimisation might be to determine a better fitting shape.
Foam::List<Foam::List<Foam::treeBoundBox>>
Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
(
    const triSurface& s
)
{
    if (!decomposer_.valid())
    {
        // Use singleton decomposeParDict. Cannot use decompositionModel
        // here since we've only got Time and not a mesh.
        if
        (
            searchableSurface::time().foundObject<IOdictionary>
            (
                "decomposeParDict"
            )
        )
        {
            decomposer_ = decompositionMethod::New
            (
                searchableSurface::time().lookupObject<IOdictionary>
                (
                    "decomposeParDict"
                )
            );
        }
        else
        {
            if (!decomposeParDict_.valid())
            {
                decomposeParDict_.reset
                (
                    new IOdictionary
                    (
                        IOobject
                        (
                            "decomposeParDict",
                            searchableSurface::time().system(),
                            searchableSurface::time(),
                            IOobject::MUST_READ,
                            IOobject::NO_WRITE
                        )
                    )
                );
            }
            decomposer_ = decompositionMethod::New(decomposeParDict_());
        }


        if (!decomposer_().parallelAware())
        {
            FatalErrorInFunction
                << "The decomposition method " << decomposer_().typeName
                << " does not decompose in parallel."
                << " Please choose one that does." << exit(FatalError);
        }

        if (!isA<geomDecomp>(decomposer_()))
        {
            FatalErrorInFunction
                << "The decomposition method " << decomposer_().typeName
                << " is not a geometric decomposition method." << endl
                << "Only geometric decomposition methods are currently"
                << " supported."
                << exit(FatalError);
        }
    }

    // Do decomposition according to triangle centre
    pointField triCentres(s.size());
    forAll(s, triI)
    {
        triCentres[triI] = s[triI].centre(s.points());
    }


    geomDecomp& decomposer = refCast<geomDecomp>(decomposer_());

    // Do the actual decomposition
    labelList distribution(decomposer.decompose(triCentres));

    // Find bounding box for all triangles on new distribution.

    // Initialise to inverted box (VGREAT, -VGREAT)
    List<List<treeBoundBox>> bbs(Pstream::nProcs());
    forAll(bbs, procI)
    {
        bbs[procI].setSize(1);
        //bbs[procI][0] = boundBox::invertedBox;
        bbs[procI][0].min() = point( VGREAT,  VGREAT,  VGREAT);
        bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
    }

    forAll(s, triI)
    {
        point& bbMin = bbs[distribution[triI]][0].min();
        point& bbMax = bbs[distribution[triI]][0].max();

        const triSurface::FaceType& f = s[triI];
        forAll(f, fp)
        {
            const point& pt = s.points()[f[fp]];
            bbMin = ::Foam::min(bbMin, pt);
            bbMax = ::Foam::max(bbMax, pt);
        }
    }

    // Now combine for all processors and convert to correct format.
    forAll(bbs, procI)
    {
        forAll(bbs[procI], i)
        {
            reduce(bbs[procI][i].min(), minOp<point>());
            reduce(bbs[procI][i].max(), maxOp<point>());
        }
    }
    return bbs;
}


// Does any part of triangle overlap bb.
bool Foam::distributedTriSurfaceMesh::overlaps
(
    const List<treeBoundBox>& bbs,
    const point& p0,
    const point& p1,
    const point& p2
)
{
    forAll(bbs, bbI)
    {
        const treeBoundBox& bb = bbs[bbI];

        treeBoundBox triBb(p0, p0);
        triBb.min() = min(triBb.min(), p1);
        triBb.min() = min(triBb.min(), p2);

        triBb.max() = max(triBb.max(), p1);
        triBb.max() = max(triBb.max(), p2);

        // Exact test of triangle intersecting bb

        // Quick rejection. If whole bounding box of tri is outside cubeBb then
        // there will be no intersection.
        if (bb.overlaps(triBb))
        {
            // Check if one or more triangle point inside
            if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
            {
                // One or more points inside
                return true;
            }

            // Now we have the difficult case: all points are outside but
            // connecting edges might go through cube. Use fast intersection
            // of bounding box.
            bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);

            if (intersect)
            {
                return true;
            }
        }
    }
    return false;
}


void Foam::distributedTriSurfaceMesh::subsetMeshMap
(
    const triSurface& s,
    const boolList& include,
    const label nIncluded,
    labelList& newToOldPoints,
    labelList& oldToNewPoints,
    labelList& newToOldFaces
)
{
    newToOldFaces.setSize(nIncluded);
    newToOldPoints.setSize(s.points().size());
    oldToNewPoints.setSize(s.points().size());
    oldToNewPoints = -1;
    {
        label faceI = 0;
        label pointI = 0;

        forAll(include, oldFacei)
        {
            if (include[oldFacei])
            {