oldCyclicPolyPatch.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "oldCyclicPolyPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "polyBoundaryMesh.H"
#include "polyMesh.H"
#include "demandDrivenData.H"
#include "OFstream.H"
#include "patchZones.H"
#include "matchPoints.H"
#include "Time.H"
#include "transformList.H"
#include "cyclicPolyPatch.H"

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

namespace Foam
{
    defineTypeNameAndDebug(oldCyclicPolyPatch, 0);

    addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, word);
    addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, dictionary);

    template<>
    const char* Foam::NamedEnum
    <
        Foam::oldCyclicPolyPatch::transformType,
        3
    >::names[] =
    {
        "unknown",
        "rotational",
        "translational"
    };
}

const Foam::NamedEnum<Foam::oldCyclicPolyPatch::transformType, 3>
    Foam::oldCyclicPolyPatch::transformTypeNames;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Foam::pointField Foam::oldCyclicPolyPatch::calcFaceCentres
(
    const UList<face>& faces,
    const pointField& points
)
{
    pointField ctrs(faces.size());

    forAll(faces, faceI)
    {
        ctrs[faceI] = faces[faceI].centre(points);
    }

    return ctrs;
}


Foam::pointField Foam::oldCyclicPolyPatch::getAnchorPoints
(
    const UList<face>& faces,
    const pointField& points
)
{
    pointField anchors(faces.size());

    forAll(faces, faceI)
    {
        anchors[faceI] = points[faces[faceI][0]];
    }

    return anchors;
}


Foam::label Foam::oldCyclicPolyPatch::findMaxArea
(
    const pointField& points,
    const faceList& faces
)
{
    label maxI = -1;
    scalar maxAreaSqr = -GREAT;

    forAll(faces, faceI)
    {
        scalar areaSqr = magSqr(faces[faceI].normal(points));

        if (areaSqr > maxAreaSqr)
        {
            maxAreaSqr = areaSqr;
            maxI = faceI;
        }
    }
    return maxI;
}


// Get geometric zones of patch by looking at normals.
// Method 1: any edge with sharpish angle is edge between two halves.
//           (this will handle e.g. wedge geometries).
//           Also two fully disconnected regions will be handled this way.
// Method 2: sort faces into two halves based on face normal.
bool Foam::oldCyclicPolyPatch::getGeometricHalves
(
    const primitivePatch& pp,
    labelList& half0ToPatch,
    labelList& half1ToPatch
) const
{
    // Calculate normals
    const vectorField& faceNormals = pp.faceNormals();

    // Find edges with sharp angles.
    boolList regionEdge(pp.nEdges(), false);

    const labelListList& edgeFaces = pp.edgeFaces();

    label nRegionEdges = 0;

    forAll(edgeFaces, edgeI)
    {
        const labelList& eFaces = edgeFaces[edgeI];

        // Check manifold edges for sharp angle.
        // (Non-manifold already handled by patchZones)
        if (eFaces.size() == 2)
        {
            if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
            {
                regionEdge[edgeI] = true;

                nRegionEdges++;
            }
        }
    }


    // For every face determine zone it is connected to (without crossing
    // any regionEdge)
    patchZones ppZones(pp, regionEdge);

    if (debug)
    {
        Pout<< "oldCyclicPolyPatch::getGeometricHalves : "
            << "Found " << nRegionEdges << " edges on patch " << name()
            << " where the cos of the angle between two connected faces"
            << " was less than " << featureCos_ << nl
            << "Patch divided by these and by single sides edges into "
            << ppZones.nZones() << " parts." << endl;


        // Dumping zones to obj files.

        labelList nZoneFaces(ppZones.nZones());

        for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
        {
            OFstream stream
            (
                boundaryMesh().mesh().time().path()
               /name()+"_zone_"+Foam::name(zoneI)+".obj"
            );
            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing zone "
                << zoneI << " face centres to OBJ file " << stream.name()
                << endl;

            labelList zoneFaces(findIndices(ppZones, zoneI));

            forAll(zoneFaces, i)
            {
                writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
            }

            nZoneFaces[zoneI] = zoneFaces.size();
        }
    }


    if (ppZones.nZones() == 2)
    {
        half0ToPatch = findIndices(ppZones, 0);
        half1ToPatch = findIndices(ppZones, 1);
    }
    else
    {
        if (debug)
        {
            Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
                << " falling back to face-normal comparison" << endl;
        }
        label n0Faces = 0;
        half0ToPatch.setSize(pp.size());

        label n1Faces = 0;
        half1ToPatch.setSize(pp.size());

        // Compare to face 0 normal.
        forAll(faceNormals, faceI)
        {
            if ((faceNormals[faceI] & faceNormals[0]) > 0)
            {
                half0ToPatch[n0Faces++] = faceI;
            }
            else
            {
                half1ToPatch[n1Faces++] = faceI;
            }
        }
        half0ToPatch.setSize(n0Faces);
        half1ToPatch.setSize(n1Faces);

        if (debug)
        {
            Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
                << " Number of faces per zone:("
                << n0Faces << ' ' << n1Faces << ')' << endl;
        }
    }

    if (half0ToPatch.size() != half1ToPatch.size())
    {
        fileName casePath(boundaryMesh().mesh().time().path());

        // Dump halves
        {
            fileName nm0(casePath/name()+"_half0_faces.obj");
            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
                << " faces to OBJ file " << nm0 << endl;
            writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());

            fileName nm1(casePath/name()+"_half1_faces.obj");
            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
                << " faces to OBJ file " << nm1 << endl;
            writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
        }

        // Dump face centres
        {
            OFstream str0(casePath/name()+"_half0.obj");
            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
                << " face centres to OBJ file " << str0.name() << endl;

            forAll(half0ToPatch, i)
            {
                writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
            }

            OFstream str1(casePath/name()+"_half1.obj");
            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
                << " face centres to OBJ file " << str1.name() << endl;
            forAll(half1ToPatch, i)
            {
                writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
            }
        }

        SeriousErrorIn
        (
            "oldCyclicPolyPatch::getGeometricHalves"
            "(const primitivePatch&, labelList&, labelList&) const"
        )   << "Patch " << name() << " gets decomposed in two zones of"
            << "inequal size: " << half0ToPatch.size()
            << " and " << half1ToPatch.size() << endl
            << "This means that the patch is either not two separate regions"
            << " or one region where the angle between the different regions"
            << " is not sufficiently sharp." << endl
            << "Please adapt the featureCos setting." << endl
            << "Continuing with incorrect face ordering from now on!" << endl;

        return false;
    }
    else
    {
        return true;
    }
}


// Given a split of faces into left and right half calculate the centres
// and anchor points. Transform the left points so they align with the
// right ones.
void Foam::oldCyclicPolyPatch::getCentresAndAnchors
(
    const primitivePatch& pp,
    const faceList& half0Faces,
    const faceList& half1Faces,

    pointField& ppPoints,
    pointField& half0Ctrs,
    pointField& half1Ctrs,
    pointField& anchors0,
    scalarField& tols
) const
{
    // Get geometric data on both halves.
    half0Ctrs = calcFaceCentres(half0Faces, pp.points());
    anchors0 = getAnchorPoints(half0Faces, pp.points());
    half1Ctrs = calcFaceCentres(half1Faces, pp.points());

    switch (transform_)
    {
        case ROTATIONAL:
        {
            label face0 = getConsistentRotationFace(half0Ctrs);
            label face1 = getConsistentRotationFace(half1Ctrs);

            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
            n0 /= mag(n0) + VSMALL;
            n1 /= mag(n1) + VSMALL;

            if (debug)
            {
                Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
                    << " Specified rotation :"
                    << " n0:" << n0 << " n1:" << n1 << endl;
            }

            // Rotation (around origin)
            const tensor reverseT(rotationTensor(n0, -n1));

            // Rotation
            forAll(half0Ctrs, faceI)
            {
                half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
                anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
            }

            ppPoints = Foam::transform(reverseT, pp.points());

            break;
        }
        //- Problem: usually specified translation is not accurate enough
        //- to get proper match so keep automatic determination over here.
        //case TRANSLATIONAL:
        //{
        //    // Transform 0 points.
        //
        //    if (debug)
        //    {
        //        Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
        //            << "Specified translation : " << separationVector_
        //            << endl;
        //    }
        //
        //    half0Ctrs += separationVector_;
        //    anchors0 += separationVector_;
        //    break;
        //}
        default:
        {
            // Assumes that cyclic is planar. This is also the initial
            // condition for patches without faces.

            // Determine the face with max area on both halves. These
            // two faces are used to determine the transformation tensors
            label max0I = findMaxArea(pp.points(), half0Faces);
            vector n0 = half0Faces[max0I].normal(pp.points());
            n0 /= mag(n0) + VSMALL;

            label max1I = findMaxArea(pp.points(), half1Faces);
            vector n1 = half1Faces[max1I].normal(pp.points());
            n1 /= mag(n1) + VSMALL;

            if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
            {
                if (debug)
                {
                    Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
                        << " Detected rotation :"
                        << " n0:" << n0 << " n1:" << n1 << endl;
                }

                // Rotation (around origin)
                const tensor reverseT(rotationTensor(n0, -n1));

                // Rotation
                forAll(half0Ctrs, faceI)
                {
                    half0Ctrs[faceI] = Foam::transform
                    (
                        reverseT,
                        half0Ctrs[faceI]
                    );
                    anchors0[faceI] = Foam::transform
                    (
                        reverseT,
                        anchors0[faceI]
                    );
                }
                ppPoints = Foam::transform(reverseT, pp.points());
            }
            else
            {
                // Parallel translation. Get average of all used points.

                primitiveFacePatch half0(half0Faces, pp.points());
                const pointField& half0Pts = half0.localPoints();
                const point ctr0(sum(half0Pts)/half0Pts.size());

                primitiveFacePatch half1(half1Faces, pp.points());
                const pointField& half1Pts = half1.localPoints();
                const point ctr1(sum(half1Pts)/half1Pts.size());

                if (debug)
                {
                    Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
                        << " Detected translation :"
                        << " n0:" << n0 << " n1:" << n1
                        << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
                }

                half0Ctrs += ctr1 - ctr0;
                anchors0 += ctr1 - ctr0;
                ppPoints = pp.points() + ctr1 - ctr0;
            }
            break;
        }
    }


    // Calculate typical distance per face
    tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
}


// Calculates faceMap and rotation. Returns true if all ok.
bool Foam::oldCyclicPolyPatch::matchAnchors
(
    const bool report,
    const primitivePatch& pp,
    const labelList& half0ToPatch,
    const pointField& anchors0,

    const labelList& half1ToPatch,
    const faceList& half1Faces,
    const labelList& from1To0,

    const scalarField& tols,

    labelList& faceMap,
    labelList& rotation
) const
{
    // Set faceMap such that half0 faces get first and corresponding half1
    // faces last.

    forAll(half0ToPatch, half0FaceI)
    {
        // Label in original patch
        label patchFaceI = half0ToPatch[half0FaceI];

        faceMap[patchFaceI] = half0FaceI;

        // No rotation
        rotation[patchFaceI] = 0;
    }

    bool fullMatch = true;

    forAll(from1To0, half1FaceI)
    {
        label patchFaceI = half1ToPatch[half1FaceI];

        // This face has to match the corresponding one on half0.
        label half0FaceI = from1To0[half1FaceI];

        label newFaceI = half0FaceI + pp.size()/2;

        faceMap[patchFaceI] = newFaceI;

        // Rotate patchFaceI such that its f[0] aligns with that of
        // the corresponding face
        // (which after shuffling will be at position half0FaceI)

        const point& wantedAnchor = anchors0[half0FaceI];

        rotation[newFaceI] = getRotation
        (
            pp.points(),
            half1Faces[half1FaceI],
            wantedAnchor,
            tols[half1FaceI]
        );

        if (rotation[newFaceI] == -1)
        {
            fullMatch = false;

            if (report)
            {
                const face& f = half1Faces[half1FaceI];
                SeriousErrorIn
                (
                    "oldCyclicPolyPatch::matchAnchors(..)"
                )   << "Patch:" << name() << " : "
                    << "Cannot find point on face " << f
                    << " with vertices:"
                    << UIndirectList<point>(pp.points(), f)()
                    << " that matches point " << wantedAnchor
                    << " when matching the halves of cyclic patch " << name()
                    << endl
                    << "Continuing with incorrect face ordering from now on!"
                    << endl;
            }
        }
    }
    return fullMatch;
}


Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
(
    const pointField& faceCentres
) const
{
    const scalarField magRadSqr =
        magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
    scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_;
    axisLen = axisLen - min(axisLen);
    const scalarField magLenSqr = magRadSqr + axisLen*axisLen;

    label rotFace = -1;
    scalar maxMagLenSqr = -GREAT;
    scalar maxMagRadSqr = -GREAT;
    forAll(faceCentres, i)
    {
        if (magLenSqr[i] >= maxMagLenSqr)
        {
            if (magRadSqr[i] > maxMagRadSqr)
            {
                rotFace = i;
                maxMagLenSqr = magLenSqr[i];
                maxMagRadSqr = magRadSqr[i];
            }
        }
    }

    if (debug)
    {
        Info<< "getConsistentRotationFace(const pointField&)" << nl
            << "    rotFace = " << rotFace << nl
            << "    point =  " << faceCentres[rotFace] << endl;
    }

    return rotFace;
}


// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //

Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
(
    const word& name,
    const label size,
    const label start,
    const label index,
    const polyBoundaryMesh& bm
)
:
    coupledPolyPatch(name, size, start, index, bm),
    featureCos_(0.9),
    transform_(UNKNOWN),
    rotationAxis_(vector::zero),
    rotationCentre_(point::zero),
    separationVector_(vector::zero)
{}


Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
(
    const word& name,
    const dictionary& dict,
    const label index,
    const polyBoundaryMesh& bm
)
:
    coupledPolyPatch(name, dict, index, bm),
    featureCos_(0.9),
    transform_(UNKNOWN),
    rotationAxis_(vector::zero),
    rotationCentre_(point::zero),
    separationVector_(vector::zero)
{
    if (dict.found("neighbourPatch"))
    {
        FatalIOErrorIn
        (
            "oldCyclicPolyPatch::oldCyclicPolyPatch\n"
            "(\n"
            "    const word& name,\n"
            "    const dictionary& dict,\n"
            "    const label index,\n"
            "    const polyBoundaryMesh& bm\n"
            ")",
            dict
        )   << "Found \"neighbourPatch\" entry when reading cyclic patch "
            << name << endl
            << "Is this mesh already with split cyclics?" << endl
            << "If so run a newer version that supports it"
            << ", if not comment out the \"neighbourPatch\" entry and re-run"
            << exit(FatalIOError);
    }

    dict.readIfPresent("featureCos", featureCos_);

    if (dict.found("transform"))
    {
        transform_ = transformTypeNames.read(dict.lookup("transform"));
        switch (transform_)
        {
            case ROTATIONAL:
            {
                dict.lookup("rotationAxis") >> rotationAxis_;
                dict.lookup("rotationCentre") >> rotationCentre_;
                break;
            }
            case TRANSLATIONAL:
            {
                dict.lookup("separationVector") >> separationVector_;
                break;
            }
            default:
            {
                // no additional info required
            }
        }
    }
}


Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
(
    const oldCyclicPolyPatch& pp,
    const polyBoundaryMesh& bm
)
:
    coupledPolyPatch(pp, bm),
    featureCos_(pp.featureCos_),
    transform_(pp.transform_),
    rotationAxis_(pp.rotationAxis_),
    rotationCentre_(pp.rotationCentre_),
    separationVector_(pp.separationVector_)
{}


Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
(
    const oldCyclicPolyPatch& pp,
    const polyBoundaryMesh& bm,
    const label index,
    const label newSize,
    const label newStart
)
:
    coupledPolyPatch(pp, bm, index, newSize, newStart),
    featureCos_(pp.featureCos_),
    transform_(pp.transform_),
    rotationAxis_(pp.rotationAxis_),
    rotationCentre_(pp.rotationCentre_),
    separationVector_(pp.separationVector_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::oldCyclicPolyPatch::~oldCyclicPolyPatch()
{}


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

void Foam::oldCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
{
    polyPatch::initGeometry(pBufs);
}


void Foam::oldCyclicPolyPatch::calcGeometry
(
    const primitivePatch& referPatch,
    const UList<point>& thisCtrs,
    const UList<point>& thisAreas,
    const UList<point>& thisCc,
    const UList<point>& nbrCtrs,
    const UList<point>& nbrAreas,
    const UList<point>& nbrCc
)
{}


void Foam::oldCyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs)
{}


void Foam::oldCyclicPolyPatch::initMovePoints
(
    PstreamBuffers& pBufs,
    const pointField& p
)
{
    polyPatch::initMovePoints(pBufs, p);
}


void Foam::oldCyclicPolyPatch::movePoints
(
    PstreamBuffers& pBufs,
    const pointField& p
)
{
    polyPatch::movePoints(pBufs, p);
}


void Foam::oldCyclicPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
{
    polyPatch::initUpdateMesh(pBufs);
}


void Foam::oldCyclicPolyPatch::updateMesh(PstreamBuffers& pBufs)
{
    polyPatch::updateMesh(pBufs);
}


void Foam::oldCyclicPolyPatch::initOrder
(
    PstreamBuffers&,
    const primitivePatch& pp
) const
{}


//  Return new ordering. Ordering is -faceMap: for every face index
//  the new face -rotation:for every new face the clockwise shift
//  of the original face. Return false if nothing changes (faceMap
//  is identity, rotation is 0)
bool Foam::oldCyclicPolyPatch::order
(
    PstreamBuffers&,
    const primitivePatch& pp,
    labelList& faceMap,
    labelList& rotation
) const
{
    faceMap.setSize(pp.size());
    faceMap = -1;

    rotation.setSize(pp.size());
    rotation = 0;

    if (pp.empty())
    {
        // No faces, nothing to change.
        return false;
    }

    if (pp.size()&1)
    {
        FatalErrorIn("oldCyclicPolyPatch::order(..)")
            << "Size of cyclic " << name() << " should be a multiple of 2"
            << ". It is " << pp.size() << abort(FatalError);
    }

    label halfSize = pp.size()/2;

    // Supplied primitivePatch already with new points.
    // Cyclics are limited to one transformation tensor
    // currently anyway (i.e. straight plane) so should not be too big a
    // problem.


    // Indices of faces on half0
    labelList half0ToPatch;
    // Indices of faces on half1
    labelList half1ToPatch;


    // 1. Test if already correctly oriented by starting from trivial ordering.
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    half0ToPatch = identity(halfSize);
    half1ToPatch = half0ToPatch + halfSize;

    // Get faces
    faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
    faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));

    // Get geometric quantities
    pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
    scalarField tols;
    getCentresAndAnchors
    (
        pp,
        half0Faces,
        half1Faces,

        ppPoints,
        half0Ctrs,
        half1Ctrs,
        anchors0,
        tols
    );

    // Geometric match of face centre vectors
    labelList from1To0;
    bool matchedAll = matchPoints
    (
        half1Ctrs,
        half0Ctrs,
        tols,
        false,
        from1To0
    );

    if (debug)
    {
        Pout<< "oldCyclicPolyPatch::order : test if already ordered:"
            << matchedAll << endl;

        // Dump halves
        fileName nm0("match1_"+name()+"_half0_faces.obj");
        Pout<< "oldCyclicPolyPatch::order : Writing half0"
            << " faces to OBJ file " << nm0 << endl;
        writeOBJ(nm0, half0Faces, ppPoints);

        fileName nm1("match1_"+name()+"_half1_faces.obj");
        Pout<< "oldCyclicPolyPatch::order : Writing half1"
            << " faces to OBJ file " << nm1 << endl;
        writeOBJ(nm1, half1Faces, pp.points());

        OFstream ccStr
        (
            boundaryMesh().mesh().time().path()
           /"match1_"+ name() + "_faceCentres.obj"
        );
        Pout<< "oldCyclicPolyPatch::order : "
            << "Dumping currently found cyclic match as lines between"
            << " corresponding face centres to file " << ccStr.name()
            << endl;

        // Recalculate untransformed face centres
        //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
        label vertI = 0;

        forAll(half1Ctrs, i)
        {
            //if (from1To0[i] != -1)
            {
                // Write edge between c1 and c0
                //const point& c0 = rawHalf0Ctrs[from1To0[i]];
                //const point& c0 = half0Ctrs[from1To0[i]];
                const point& c0 = half0Ctrs[i];
                const point& c1 = half1Ctrs[i];
                writeOBJ(ccStr, c0, c1, vertI);
            }
        }
    }


    // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if (!matchedAll)
    {
        label faceI = 0;
        for (label i = 0; i < halfSize; i++)
        {
            half0ToPatch[i] = faceI++;
            half1ToPatch[i] = faceI++;
        }

        // And redo all matching
        half0Faces = UIndirectList<face>(pp, half0ToPatch);
        half1Faces = UIndirectList<face>(pp, half1ToPatch);

        getCentresAndAnchors
        (
            pp,
            half0Faces,
            half1Faces,

            ppPoints,
            half0Ctrs,
            half1Ctrs,
            anchors0,
            tols
        );

        // Geometric match of face centre vectors
        matchedAll = matchPoints
        (
            half1Ctrs,
            half0Ctrs,
            tols,
            false,
            from1To0
        );

        if (debug)
        {
            Pout<< "oldCyclicPolyPatch::order : test if pairwise ordered:"
                << matchedAll << endl;
            // Dump halves
            fileName nm0("match2_"+name()+"_half0_faces.obj");
            Pout<< "oldCyclicPolyPatch::order : Writing half0"
                << " faces to OBJ file " << nm0 << endl;
            writeOBJ(nm0, half0Faces, ppPoints);

            fileName nm1("match2_"+name()+"_half1_faces.obj");
            Pout<< "oldCyclicPolyPatch::order : Writing half1"
                << " faces to OBJ file " << nm1 << endl;
            writeOBJ(nm1, half1Faces, pp.points());

            OFstream ccStr
            (
                boundaryMesh().mesh().time().path()
               /"match2_"+name()+"_faceCentres.obj"
            );
            Pout<< "oldCyclicPolyPatch::order : "
                << "Dumping currently found cyclic match as lines between"
                << " corresponding face centres to file " << ccStr.name()
                << endl;

            // Recalculate untransformed face centres
            label vertI = 0;

            forAll(half1Ctrs, i)
            {
                if (from1To0[i] != -1)
                {
                    // Write edge between c1 and c0
                    const point& c0 = half0Ctrs[from1To0[i]];
                    const point& c1 = half1Ctrs[i];
                    writeOBJ(ccStr, c0, c1, vertI);
                }
            }
        }
    }


    // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if (!matchedAll)
    {
        label baffleI = 0;

        forAll(pp, faceI)
        {
            const face& f = pp.localFaces()[faceI];
            const labelList& pFaces = pp.pointFaces()[f[0]];

            label matchedFaceI = -1;

            forAll(pFaces, i)
            {
                label otherFaceI = pFaces[i];

                if (otherFaceI > faceI)
                {
                    const face& otherF = pp.localFaces()[otherFaceI];

                    // Note: might pick up two similar oriented faces
                    //       (but that is illegal anyway)
                    if (f == otherF)
                    {
                        matchedFaceI = otherFaceI;
                        break;
                    }
                }
            }

            if (matchedFaceI != -1)