View Issue Details

IDProjectCategoryView StatusLast Update
0001479OpenFOAMBugpublic2015-11-21 21:37
ReporterGRAUPS Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformIntel64OSRHELOS Version6.4
Summary0001479: non-axis-aligned flow rate monitor in porous zone
DescriptionSometimes the flow-rate face source monitors in openfoam models are not reporting the correct flow-rates (verified through post-processing). I was able to reproduce this problem with a non-axis-aligned simple pipe.

Please find attached a folder with two pipe models. One model has a porous zone, one does not. Both have a constant flow rate at the inlet of 1 kg/s and both have the monitoring plane in the middle of the pipe. I found that the problem only seems to exist when I add a porous zone into the model. The model with the porous zone is reporting half the flow-rate that it should be at the center of the pipe, where the model without the porous zone (the second included model) reports this correctly.

Let me know if you have any questions, and Happy New Year!
Steps To ReproducePlease see the run_openfoam.sh scripts for each case
Additional InformationI have a feeling this might be related to face normals. I assumed that it would have been resolved though by making sure to put phi under the orientedFaces header. That doesn't seem to be the case though.
TagsNo tags attached.

Activities

GRAUPS

2015-01-08 16:52

reporter  

non_aligned_pipes.tar.gz (215,829 bytes)

user4

2015-03-02 12:53

  ~0003940

It seems the faceZones are not oriented correctly across the processor patches and this causes the sum to be incorrect. This could be due to any one of snappyHexMesh, createPatch, reconstructParMesh, decomposePar, createBaffles or renumberMesh. Do you see the same problems running non-parallel?

GRAUPS

2015-03-04 18:49

reporter   ~0003956

Mattijs, thanks for taking a look at this. I did some more testing and isolating, here is a summary of what I modified and uploaded in the porous setup.

1.) Converted everything to serial
2.) Removed reconstructParMesh, decomposePar, createBaffles, and renumberMesh (createBaffles was a placeholder and not really doing anything)
3.) Manually removed 0 sized patches in order to eliminate the createPatch command

After rerunning the new case using the above modifications, the issue remains. The only command left in the list from your previous note is snappyHexMesh. So the problem with faceZones not being oriented properly may be in there.

I attached my modified case, let me know if you need additional troubleshooting. Thanks!

GRAUPS

2015-03-04 18:50

reporter  

user4

2015-03-12 17:02

  ~0004096

I thought we fixed the snappyHexMesh zone orientation a while ago. I'll check.

GRAUPS

2015-03-16 16:58

reporter   ~0004147

Thanks mattijs for checking, I look forward to a possible resolution.

user4

2015-03-30 16:34

 

normals.png (68,491 bytes)   
normals.png (68,491 bytes)   

user4

2015-03-30 16:39

  ~0004552

The 'porous' faceZone consists of two parts. These are correctly oriented pointing out of the corresponding cellZone (see normals.png). The problem one is the 'int' faceZone which does not get a consistent orientation.

GRAUPS

2015-03-30 18:12

reporter   ~0004553

I would agree, the inconsistent orientation of the int facezone is likely the cause of the problem. Note, this problem of inconsistent normals only seems to present itself when I add a cell zone to the pipe and put the int faceZone inside the cellZone. If I remove the cellZone and keep everything else the same, the int faceZone gets oriented properly.

So the problem seems to appear in the presence of a cell zone, if that helps narrow it down at all.

user4

2015-04-10 15:14

 

meshRefinementBaffles.C (100,964 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     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 "meshRefinement.H"
#include "fvMesh.H"
#include "syncTools.H"
#include "Time.H"
#include "refinementSurfaces.H"
#include "pointSet.H"
#include "faceSet.H"
#include "indirectPrimitivePatch.H"
#include "polyTopoChange.H"
#include "meshTools.H"
#include "polyModifyFace.H"
#include "polyModifyCell.H"
#include "polyAddFace.H"
#include "polyRemoveFace.H"
#include "polyAddPoint.H"
#include "localPointRegion.H"
#include "duplicatePoints.H"
#include "OFstream.H"
#include "regionSplit.H"
#include "removeCells.H"
#include "unitConversion.H"
#include "OBJstream.H"
#include "patchFaceOrientation.H"
#include "PatchEdgeFaceWave.H"
#include "patchEdgeFaceRegion.H"

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

// Repatches external face or creates baffle for internal face
// with user specified patches (might be different for both sides).
// Returns label of added face.
Foam::label Foam::meshRefinement::createBaffle
(
    const label faceI,
    const label ownPatch,
    const label neiPatch,
    polyTopoChange& meshMod
) const
{
    const face& f = mesh_.faces()[faceI];
    label zoneID = mesh_.faceZones().whichZone(faceI);
    bool zoneFlip = false;

    if (zoneID >= 0)
    {
        const faceZone& fZone = mesh_.faceZones()[zoneID];
        zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
    }

    meshMod.setAction
    (
        polyModifyFace
        (
            f,                          // modified face
            faceI,                      // label of face
            mesh_.faceOwner()[faceI],   // owner
            -1,                         // neighbour
            false,                      // face flip
            ownPatch,                   // patch for face
            false,                      // remove from zone
            zoneID,                     // zone for face
            zoneFlip                    // face flip in zone
        )
    );


    label dupFaceI = -1;

    if (mesh_.isInternalFace(faceI))
    {
        if (neiPatch == -1)
        {
            FatalErrorIn
            (
                "meshRefinement::createBaffle"
                "(const label, const label, const label, polyTopoChange&)"
            )   << "No neighbour patch for internal face " << faceI
                << " fc:" << mesh_.faceCentres()[faceI]
                << " ownPatch:" << ownPatch << abort(FatalError);
        }

        bool reverseFlip = false;
        if (zoneID >= 0)
        {
            reverseFlip = !zoneFlip;
        }

        dupFaceI = meshMod.setAction
        (
            polyAddFace
            (
                f.reverseFace(),            // modified face
                mesh_.faceNeighbour()[faceI],// owner
                -1,                         // neighbour
                -1,                         // masterPointID
                -1,                         // masterEdgeID
                faceI,                      // masterFaceID,
                true,                       // face flip
                neiPatch,                   // patch for face
                zoneID,                     // zone for face
                reverseFlip                 // face flip in zone
            )
        );
    }
    return dupFaceI;
}


//// Check if we are a boundary face and normal of surface does
//// not align with test vector. In this case there'd probably be
//// a freestanding 'baffle' so we might as well not create it.
//// Note that since it is not a proper baffle we cannot detect it
//// afterwards so this code cannot be merged with the
//// filterDuplicateFaces code.
//bool Foam::meshRefinement::validBaffleTopology
//(
//    const label faceI,
//    const vector& n1,
//    const vector& n2,
//    const vector& testDir
//) const
//{
//
//    label patchI = mesh_.boundaryMesh().whichPatch(faceI);
//    if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
//    {
//        return true;
//    }
//    else if (mag(n1&n2) > cos(degToRad(30)))
//    {
//        // Both normals aligned. Check that test vector perpendicularish to
//        // surface normal
//        scalar magTestDir = mag(testDir);
//        if (magTestDir > VSMALL)
//        {
//            if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
//            {
//                //Pout<< "** disabling baffling face "
//                //    << mesh_.faceCentres()[faceI] << endl;
//                return false;
//            }
//        }
//    }
//    return true;
//}


// Determine patches for baffles on all intersected unnamed faces
void Foam::meshRefinement::getBafflePatches
(
    const labelList& globalToMasterPatch,
    const labelList& neiLevel,
    const pointField& neiCc,

    labelList& ownPatch,
    labelList& neiPatch
) const
{
    autoPtr<OFstream> str;
    label vertI = 0;
    if (debug&OBJINTERSECTIONS)
    {
        mkDir(mesh_.time().path()/timeName());
        str.reset
        (
            new OFstream
            (
                mesh_.time().path()/timeName()/"intersections.obj"
            )
        );

        Pout<< "getBafflePatches : Writing surface intersections to file "
            << str().name() << nl << endl;
    }

    const pointField& cellCentres = mesh_.cellCentres();

    // Surfaces that need to be baffled
    const labelList surfacesToBaffle
    (
        surfaceZonesInfo::getUnnamedSurfaces(surfaces_.surfZones())
    );

    ownPatch.setSize(mesh_.nFaces());
    ownPatch = -1;
    neiPatch.setSize(mesh_.nFaces());
    neiPatch = -1;


    // Collect candidate faces
    // ~~~~~~~~~~~~~~~~~~~~~~~

    labelList testFaces(intersectedFaces());

    // Collect segments
    // ~~~~~~~~~~~~~~~~

    pointField start(testFaces.size());
    pointField end(testFaces.size());

    forAll(testFaces, i)
    {
        label faceI = testFaces[i];

        label own = mesh_.faceOwner()[faceI];

        if (mesh_.isInternalFace(faceI))
        {
            start[i] = cellCentres[own];
            end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
        }
        else
        {
            start[i] = cellCentres[own];
            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
        }
    }

    // Extend segments a bit
    {
        const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
        start -= smallVec;
        end += smallVec;
    }


    // Do test for intersections
    // ~~~~~~~~~~~~~~~~~~~~~~~~~
    labelList surface1;
    List<pointIndexHit> hit1;
    labelList region1;
    labelList surface2;
    List<pointIndexHit> hit2;
    labelList region2;
    surfaces_.findNearestIntersection
    (
        surfacesToBaffle,
        start,
        end,

        surface1,
        hit1,
        region1,

        surface2,
        hit2,
        region2
    );

    forAll(testFaces, i)
    {
        label faceI = testFaces[i];

        if (hit1[i].hit() && hit2[i].hit())
        {
            if (str.valid())
            {
                meshTools::writeOBJ(str(), start[i]);
                vertI++;
                meshTools::writeOBJ(str(), hit1[i].rawPoint());
                vertI++;
                meshTools::writeOBJ(str(), hit2[i].rawPoint());
                vertI++;
                meshTools::writeOBJ(str(), end[i]);
                vertI++;
                str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
                str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
                str()<< "l " << vertI-1 << ' ' << vertI << nl;
            }

            // Pick up the patches
            ownPatch[faceI] = globalToMasterPatch
            [
                surfaces_.globalRegion(surface1[i], region1[i])
            ];
            neiPatch[faceI] = globalToMasterPatch
            [
                surfaces_.globalRegion(surface2[i], region2[i])
            ];

            if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
            {
                FatalErrorIn("getBafflePatches(..)")
                    << "problem." << abort(FatalError);
            }
        }
    }

    // No need to parallel sync since intersection data (surfaceIndex_ etc.)
    // already guaranteed to be synced...
    // However:
    // - owncc and neicc are reversed on different procs so might pick
    //   up different regions reversed? No problem. Neighbour on one processor
    //   might not be owner on the other processor but the neighbour is
    //   not used when creating baffles from proc faces.
    // - tolerances issues occasionally crop up.
    syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
    syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>());
}


Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches
(
    const bool allowBoundary,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch
) const
{
    Map<labelPair> bafflePatch(mesh_.nFaces()/1000);

    const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
    const faceZoneMesh& fZones = mesh_.faceZones();

    forAll(surfZones, surfI)
    {
        const word& faceZoneName = surfZones[surfI].faceZoneName();

        if (faceZoneName.size())
        {
            // Get zone
            label zoneI = fZones.findZoneID(faceZoneName);

            const faceZone& fZone = fZones[zoneI];

            // Get patch allocated for zone
            label globalRegionI = surfaces_.globalRegion(surfI, 0);
            labelPair zPatches
            (
                globalToMasterPatch[globalRegionI],
                globalToSlavePatch[globalRegionI]
            );

            Info<< "For zone " << fZone.name() << " found patches "
                << mesh_.boundaryMesh()[zPatches[0]].name() << " and "
                << mesh_.boundaryMesh()[zPatches[1]].name()
                << endl;

            forAll(fZone, i)
            {
                label faceI = fZone[i];

                if (allowBoundary || mesh_.isInternalFace(faceI))
                {
                    labelPair patches = zPatches;
                    if (fZone.flipMap()[i])
                    {
                       patches = reverse(patches);
                    }

                    if (!bafflePatch.insert(faceI, patches))
                    {
                        FatalErrorIn("getZoneBafflePatches(..)")
                            << "Face " << faceI
                            << " fc:" << mesh_.faceCentres()[faceI]
                            << " in zone " << fZone.name()
                            << " is in multiple zones!"
                            << abort(FatalError);
                    }
                }
            }
        }
    }
    return bafflePatch;
}


// Create baffle for every face where ownPatch != -1
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
(
    const labelList& ownPatch,
    const labelList& neiPatch
)
{
    if
    (
        ownPatch.size() != mesh_.nFaces()
     || neiPatch.size() != mesh_.nFaces()
    )
    {
        FatalErrorIn
        (
            "meshRefinement::createBaffles"
            "(const labelList&, const labelList&)"
        )   << "Illegal size :"
            << " ownPatch:" << ownPatch.size()
            << " neiPatch:" << neiPatch.size()
            << ". Should be number of faces:" << mesh_.nFaces()
            << abort(FatalError);
    }

    if (debug)
    {
        labelList syncedOwnPatch(ownPatch);
        syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>());
        labelList syncedNeiPatch(neiPatch);
        syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>());

        forAll(syncedOwnPatch, faceI)
        {
            if
            (
                (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
             || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
            )
            {
                FatalErrorIn
                (
                    "meshRefinement::createBaffles"
                    "(const labelList&, const labelList&)"
                )   << "Non synchronised at face:" << faceI
                    << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
                    << " fc:" << mesh_.faceCentres()[faceI] << endl
                    << "ownPatch:" << ownPatch[faceI]
                    << " syncedOwnPatch:" << syncedOwnPatch[faceI]
                    << " neiPatch:" << neiPatch[faceI]
                    << " syncedNeiPatch:" << syncedNeiPatch[faceI]
                    << abort(FatalError);
            }
        }
    }

    // Topochange container
    polyTopoChange meshMod(mesh_);

    label nBaffles = 0;

    forAll(ownPatch, faceI)
    {
        if (ownPatch[faceI] != -1)
        {
            // Create baffle or repatch face. Return label of inserted baffle
            // face.
            createBaffle
            (
                faceI,
                ownPatch[faceI],   // owner side patch
                neiPatch[faceI],   // neighbour side patch
                meshMod
            );
            nBaffles++;
        }
    }

    // Change the mesh (no inflation, parallel sync)
    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);

    // Update fields
    mesh_.updateMesh(map);

    // Move mesh if in inflation mode
    if (map().hasMotionPoints())
    {
        mesh_.movePoints(map().preMotionPoints());
    }
    else
    {
        // Delete mesh volumes.
        mesh_.clearOut();
    }


    // Reset the instance for if in overwrite mode
    mesh_.setInstance(timeName());

    //- Redo the intersections on the newly create baffle faces. Note that
    //  this changes also the cell centre positions.
    faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);

    const labelList& reverseFaceMap = map().reverseFaceMap();
    const labelList& faceMap = map().faceMap();

    // Pick up owner side of baffle
    forAll(ownPatch, oldFaceI)
    {
        label faceI = reverseFaceMap[oldFaceI];

        if (ownPatch[oldFaceI] != -1 && faceI >= 0)
        {
            const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];

            forAll(ownFaces, i)
            {
                baffledFacesSet.insert(ownFaces[i]);
            }
        }
    }
    // Pick up neighbour side of baffle (added faces)
    forAll(faceMap, faceI)
    {
        label oldFaceI = faceMap[faceI];

        if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
        {
            const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];

            forAll(ownFaces, i)
            {
                baffledFacesSet.insert(ownFaces[i]);
            }
        }
    }
    baffledFacesSet.sync(mesh_);

    updateMesh(map, baffledFacesSet.toc());

    return map;
}


void Foam::meshRefinement::checkZoneFaces() const
{
    const faceZoneMesh& fZones = mesh_.faceZones();

    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();

    forAll(pbm, patchI)
    {
        const polyPatch& pp = pbm[patchI];

        if (isA<processorPolyPatch>(pp))
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                label zoneI = fZones.whichZone(faceI);

                if (zoneI != -1)
                {
                    FatalErrorIn("meshRefinement::checkZoneFaces")
                        << "face:" << faceI << " on patch " << pp.name()
                        << " is in zone " << fZones[zoneI].name()
                        << exit(FatalError);
                }
            }
        }
    }
}


Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
(
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch,
    List<labelPair>& baffles
)
{
    const labelList zonedSurfaces
    (
        surfaceZonesInfo::getNamedSurfaces(surfaces_.surfZones())
    );

    autoPtr<mapPolyMesh> map;

    // No need to sync; all processors will have all same zonedSurfaces.
    if (zonedSurfaces.size())
    {
        // Split internal faces on interface surfaces
        Info<< "Converting zoned faces into baffles ..." << endl;

        // Get faces (internal only) to be baffled. Map from face to patch
        // label.
        Map<labelPair> faceToPatch
        (
            getZoneBafflePatches
            (
                false,
                globalToMasterPatch,
                globalToSlavePatch
            )
        );

        label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
        if (nZoneFaces > 0)
        {
            // Convert into labelLists
            labelList ownPatch(mesh_.nFaces(), -1);
            labelList neiPatch(mesh_.nFaces(), -1);
            forAllConstIter(Map<labelPair>, faceToPatch, iter)
            {
                ownPatch[iter.key()] = iter().first();
                neiPatch[iter.key()] = iter().second();
            }

            // Create baffles. both sides same patch.
            map = createBaffles(ownPatch, neiPatch);

            // Get pairs of faces created.
            // Just loop over faceMap and store baffle if we encounter a slave
            // face.

            baffles.setSize(faceToPatch.size());
            label baffleI = 0;

            const labelList& faceMap = map().faceMap();
            const labelList& reverseFaceMap = map().reverseFaceMap();

            forAll(faceMap, faceI)
            {
                label oldFaceI = faceMap[faceI];

                // Does face originate from face-to-patch
                Map<labelPair>::const_iterator iter = faceToPatch.find
                (
                    oldFaceI
                );

                if (iter != faceToPatch.end())
                {
                    label masterFaceI = reverseFaceMap[oldFaceI];
                    if (faceI != masterFaceI)
                    {
                        baffles[baffleI++] = labelPair(masterFaceI, faceI);
                    }
                }
            }

            if (baffleI != faceToPatch.size())
            {
                FatalErrorIn("meshRefinement::createZoneBaffles(..)")
                    << "Had " << faceToPatch.size() << " patches to create "
                    << " but encountered " << baffleI
                    << " slave faces originating from patcheable faces."
                    << abort(FatalError);
            }

            if (debug&MESH)
            {
                const_cast<Time&>(mesh_.time())++;
                Pout<< "Writing zone-baffled mesh to time " << timeName()
                    << endl;
                write
                (
                    debugType(debug),
                    writeType(writeLevel() | WRITEMESH),
                    mesh_.time().path()/"baffles"
                );
            }
        }
        Info<< "Created " << nZoneFaces << " baffles in = "
            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
    }
    return map;
}


// Extract those baffles (duplicate) faces that are on the edge of a baffle
// region. These are candidates for merging.
// Done by counting the number of baffles faces per mesh edge. If edge
// has 2 boundary faces and both are baffle faces it is the edge of a baffle
// region.
Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
(
    const List<labelPair>& couples,
    const scalar planarAngle
) const
{
    // All duplicate faces on edge of the patch are to be merged.
    // So we count for all edges of duplicate faces how many duplicate
    // faces use them.
    labelList nBafflesPerEdge(mesh_.nEdges(), 0);


    // This algorithm is quite tricky. We don't want to use edgeFaces and
    // also want it to run in parallel so it is now an algorithm over
    // all (boundary) faces instead.
    // We want to pick up any edges that are only used by the baffle
    // or internal faces but not by any other boundary faces. So
    // - increment count on an edge by 1 if it is used by any (uncoupled)
    //   boundary face.
    // - increment count on an edge by 1000000 if it is used by a baffle face
    // - sum in parallel
    //
    // So now any edge that is used by baffle faces only will have the
    // value 2*1000000+2*1.


    const label baffleValue = 1000000;


    // Count number of boundary faces per edge
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        // Count number of boundary faces. Discard coupled boundary faces.
        if (!pp.coupled())
        {
            label faceI = pp.start();

            forAll(pp, i)
            {
                const labelList& fEdges = mesh_.faceEdges(faceI);

                forAll(fEdges, fEdgeI)
                {
                    nBafflesPerEdge[fEdges[fEdgeI]]++;
                }
                faceI++;
            }
        }
    }


    DynamicList<label> fe0;
    DynamicList<label> fe1;


    // Count number of duplicate boundary faces per edge
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    forAll(couples, i)
    {
        {
            label f0 = couples[i].first();
            const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
            forAll(fEdges0, fEdgeI)
            {
                nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
            }
        }

        {
            label f1 = couples[i].second();
            const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
            forAll(fEdges1, fEdgeI)
            {
                nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
            }
        }
    }

    // Add nBaffles on shared edges
    syncTools::syncEdgeList
    (
        mesh_,
        nBafflesPerEdge,
        plusEqOp<label>(),  // in-place add
        label(0)            // initial value
    );


    // Baffles which are not next to other boundaries and baffles will have
    // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which
    // are both baffle faces)

    List<labelPair> filteredCouples(couples.size());
    label filterI = 0;

    forAll(couples, i)
    {
        const labelPair& couple = couples[i];

        if
        (
            patches.whichPatch(couple.first())
         == patches.whichPatch(couple.second())
        )
        {
            const labelList& fEdges = mesh_.faceEdges(couple.first());

            forAll(fEdges, fEdgeI)
            {
                label edgeI = fEdges[fEdgeI];

                if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
                {
                    filteredCouples[filterI++] = couple;
                    break;
                }
            }
        }
    }
    filteredCouples.setSize(filterI);


    label nFiltered = returnReduce(filteredCouples.size(), sumOp<label>());

    Info<< "freeStandingBaffles : detected "
        << nFiltered
        << " free-standing baffles out of "
        << returnReduce(couples.size(), sumOp<label>())
        << nl << endl;


    if (nFiltered > 0)
    {
        // Collect segments
        // ~~~~~~~~~~~~~~~~

        pointField start(filteredCouples.size());
        pointField end(filteredCouples.size());

        const pointField& cellCentres = mesh_.cellCentres();

        forAll(filteredCouples, i)
        {
            const labelPair& couple = filteredCouples[i];
            start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
            end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
        }

        // Extend segments a bit
        {
            const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
            start -= smallVec;
            end += smallVec;
        }


        // Do test for intersections
        // ~~~~~~~~~~~~~~~~~~~~~~~~~
        labelList surface1;
        List<pointIndexHit> hit1;
        labelList region1;
        vectorField normal1;

        labelList surface2;
        List<pointIndexHit> hit2;
        labelList region2;
        vectorField normal2;

        surfaces_.findNearestIntersection
        (
            identity(surfaces_.surfaces().size()),
            start,
            end,

            surface1,
            hit1,
            region1,
            normal1,

            surface2,
            hit2,
            region2,
            normal2
        );

        //mkDir(mesh_.time().path()/timeName());
        //OBJstream str
        //(
        //    mesh_.time().path()/timeName()/"flatBaffles.obj"
        //);

        const scalar planarAngleCos = Foam::cos(degToRad(planarAngle));

        label filterI = 0;
        forAll(filteredCouples, i)
        {
            const labelPair& couple = filteredCouples[i];

            if
            (
                hit1[i].hit()
             && hit2[i].hit()
             && (
                    surface1[i] != surface2[i]
                 || hit1[i].index() != hit2[i].index()
                )
            )
            {
                // Two different hits. Check angle.
                //str.write
                //(
                //    linePointRef(hit1[i].hitPoint(), hit2[i].hitPoint()),
                //    normal1[i],
                //    normal2[i]
                //);

                if ((normal1[i]&normal2[i]) > planarAngleCos)
                {
                    // Both normals aligned
                    vector n = end[i]-start[i];
                    scalar magN = mag(n);
                    if (magN > VSMALL)
                    {
                        filteredCouples[filterI++] = couple;
                    }
                }
            }
            else if (hit1[i].hit() || hit2[i].hit())
            {
                // Single hit. Do not include in freestanding baffles.
            }
        }

        filteredCouples.setSize(filterI);

        Info<< "freeStandingBaffles : detected "
            << returnReduce(filterI, sumOp<label>())
            << " planar (within " << planarAngle
            << " degrees) free-standing baffles out of "
            << nFiltered
            << nl << endl;
    }

    return filteredCouples;
}


// Merge baffles. Gets pairs of faces.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
(
    const List<labelPair>& couples
)
{
    // Mesh change engine
    polyTopoChange meshMod(mesh_);

    const faceList& faces = mesh_.faces();
    const labelList& faceOwner = mesh_.faceOwner();
    const faceZoneMesh& faceZones = mesh_.faceZones();

    forAll(couples, i)
    {
        label face0 = couples[i].first();
        label face1 = couples[i].second();

        // face1 < 0 signals a coupled face that has been converted to baffle.

        label own0 = faceOwner[face0];
        label own1 = faceOwner[face1];

        if (face1 < 0 || own0 < own1)
        {
            // Use face0 as the new internal face.
            label zoneID = faceZones.whichZone(face0);
            bool zoneFlip = false;

            if (zoneID >= 0)
            {
                const faceZone& fZone = faceZones[zoneID];
                zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
            }

            label nei = (face1 < 0 ? -1 : own1);

            meshMod.setAction(polyRemoveFace(face1));
            meshMod.setAction
            (
                polyModifyFace
                (
                    faces[face0],           // modified face
                    face0,                  // label of face being modified
                    own0,                   // owner
                    nei,                    // neighbour
                    false,                  // face flip
                    -1,                     // patch for face
                    false,                  // remove from zone
                    zoneID,                 // zone for face
                    zoneFlip                // face flip in zone
                )
            );
        }
        else
        {
            // Use face1 as the new internal face.
            label zoneID = faceZones.whichZone(face1);
            bool zoneFlip = false;

            if (zoneID >= 0)
            {
                const faceZone& fZone = faceZones[zoneID];
                zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
            }

            meshMod.setAction(polyRemoveFace(face0));
            meshMod.setAction
            (
                polyModifyFace
                (
                    faces[face1],           // modified face
                    face1,                  // label of face being modified
                    own1,                   // owner
                    own0,                   // neighbour
                    false,                  // face flip
                    -1,                     // patch for face
                    false,                  // remove from zone
                    zoneID,                 // zone for face
                    zoneFlip                // face flip in zone
                )
            );
        }
    }

    // Change the mesh (no inflation)
    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);

    // Update fields
    mesh_.updateMesh(map);

    // Move mesh (since morphing does not do this)
    if (map().hasMotionPoints())
    {
        mesh_.movePoints(map().preMotionPoints());
    }
    else
    {
        // Delete mesh volumes.
        mesh_.clearOut();
    }

    // Reset the instance for if in overwrite mode
    mesh_.setInstance(timeName());

    // Update intersections. Recalculate intersections on merged faces since
    // this seems to give problems? Note: should not be necessary since
    // baffles preserve intersections from when they were created.
    labelList newExposedFaces(2*couples.size());
    label newI = 0;

    forAll(couples, i)
    {
        label newFace0 = map().reverseFaceMap()[couples[i].first()];
        if (newFace0 != -1)
        {
            newExposedFaces[newI++] = newFace0;
        }

        label newFace1 = map().reverseFaceMap()[couples[i].second()];
        if (newFace1 != -1)
        {
            newExposedFaces[newI++] = newFace1;
        }
    }
    newExposedFaces.setSize(newI);
    updateMesh(map, newExposedFaces);

    return map;
}


// Finds region per cell for cells inside closed named surfaces
void Foam::meshRefinement::findCellZoneGeometric
(
    const pointField& neiCc,
    const labelList& closedNamedSurfaces,   // indices of closed surfaces
    labelList& namedSurfaceIndex,           // per face index of named surface
    const labelList& surfaceToCellZone,     // cell zone index per surface

    labelList& cellToZone
) const
{
    const pointField& cellCentres = mesh_.cellCentres();
    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();

    // Check if cell centre is inside
    labelList insideSurfaces;
    surfaces_.findInside
    (
        closedNamedSurfaces,
        cellCentres,
        insideSurfaces
    );

    forAll(insideSurfaces, cellI)
    {
        if (cellToZone[cellI] == -2)
        {
            label surfI = insideSurfaces[cellI];

            if (surfI != -1)
            {
                cellToZone[cellI] = surfaceToCellZone[surfI];
            }
        }
    }


    // Some cells with cell centres close to surface might have
    // had been put into wrong surface. Recheck with perturbed cell centre.


    // 1. Collect points

    // Count points to test.
    label nCandidates = 0;
    forAll(namedSurfaceIndex, faceI)
    {
        label surfI = namedSurfaceIndex[faceI];

        if (surfI != -1)
        {
            if (mesh_.isInternalFace(faceI))
            {
                nCandidates += 2;
            }
            else
            {
                nCandidates += 1;
            }
        }
    }

    // Collect points.
    pointField candidatePoints(nCandidates);
    nCandidates = 0;
    forAll(namedSurfaceIndex, faceI)
    {
        label surfI = namedSurfaceIndex[faceI];

        if (surfI != -1)
        {
            label own = faceOwner[faceI];
            const point& ownCc = cellCentres[own];

            if (mesh_.isInternalFace(faceI))
            {
                label nei = faceNeighbour[faceI];
                const point& neiCc = cellCentres[nei];
                // Perturbed cc
                const vector d = 1e-4*(neiCc - ownCc);
                candidatePoints[nCandidates++] = ownCc-d;
                candidatePoints[nCandidates++] = neiCc+d;
            }
            else
            {
                //const point& neiFc = mesh_.faceCentres()[faceI];
                const point& neiFc = neiCc[faceI-mesh_.nInternalFaces()];

                // Perturbed cc
                const vector d = 1e-4*(neiFc - ownCc);
                candidatePoints[nCandidates++] = ownCc-d;
            }
        }
    }


    // 2. Test points for inside

    surfaces_.findInside
    (
        closedNamedSurfaces,
        candidatePoints,
        insideSurfaces
    );


    // 3. Update zone information

    nCandidates = 0;
    forAll(namedSurfaceIndex, faceI)
    {
        label surfI = namedSurfaceIndex[faceI];

        if (surfI != -1)
        {
            label own = faceOwner[faceI];

            if (mesh_.isInternalFace(faceI))
            {
                label ownSurfI = insideSurfaces[nCandidates++];
                if (ownSurfI != -1)
                {
                    cellToZone[own] = surfaceToCellZone[ownSurfI];
                }

                label neiSurfI = insideSurfaces[nCandidates++];
                if (neiSurfI != -1)
                {
                    label nei = faceNeighbour[faceI];

                    cellToZone[nei] = surfaceToCellZone[neiSurfI];
                }
            }
            else
            {
                label ownSurfI = insideSurfaces[nCandidates++];
                if (ownSurfI != -1)
                {
                    cellToZone[own] = surfaceToCellZone[ownSurfI];
                }
            }
        }
    }


    // Adapt the namedSurfaceIndex
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // for if any cells were not completely covered.

    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
    {
        label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
        label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];

        if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
        {
            // Give face the zone of max cell zone
            namedSurfaceIndex[faceI] = findIndex
            (
                surfaceToCellZone,
                max(ownZone, neiZone)
            );
        }
    }

    labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        if (pp.coupled())
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
                neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
            }
        }
    }
    syncTools::swapBoundaryFaceList(mesh_, neiCellZone);

    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        if (pp.coupled())
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
                label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];

                if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
                {
                    // Give face the max cell zone
                    namedSurfaceIndex[faceI] = findIndex
                    (
                        surfaceToCellZone,
                        max(ownZone, neiZone)
                    );
                }
            }
        }
    }

    // Sync
    syncTools::syncFaceList(mesh_, namedSurfaceIndex, maxEqOp<label>());
}


void Foam::meshRefinement::findCellZoneInsideWalk
(
    const labelList& locationSurfaces,  // indices of surfaces with inside point
    const labelList& namedSurfaceIndex, // per face index of named surface
    const labelList& surfaceToCellZone, // cell zone index per surface

    labelList& cellToZone
) const
{
    // Analyse regions. Reuse regionsplit
    boolList blockedFace(mesh_.nFaces());
    //selectSeparatedCoupledFaces(blockedFace);

    forAll(namedSurfaceIndex, faceI)
    {
        if (namedSurfaceIndex[faceI] == -1)
        {
            blockedFace[faceI] = false;
        }
        else
        {
            blockedFace[faceI] = true;
        }
    }
    // No need to sync since namedSurfaceIndex already is synced

    // Set region per cell based on walking
    regionSplit cellRegion(mesh_, blockedFace);
    blockedFace.clear();


    // Force calculation of face decomposition (used in findCell)
    (void)mesh_.tetBasePtIs();

    const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();

    // For all locationSurface find the cell
    forAll(locationSurfaces, i)
    {
        label surfI = locationSurfaces[i];

        const point& insidePoint = surfZones[surfI].zoneInsidePoint();

        Info<< "For surface " << surfaces_.names()[surfI]
            << " finding inside point " << insidePoint
            << endl;

        // Find the region containing the insidePoint
        label keepRegionI = findRegion
        (
            mesh_,
            cellRegion,
            mergeDistance_*vector(1,1,1),
            insidePoint
        );

        Info<< "For surface " << surfaces_.names()[surfI]
            << " found point " << insidePoint
            << " in global region " << keepRegionI
            << " out of " << cellRegion.nRegions() << " regions." << endl;

        if (keepRegionI == -1)
        {
            FatalErrorIn
            (
                "meshRefinement::findCellZoneInsideWalk"
                "(const labelList&, const labelList&"
                ", const labelList&, const labelList&)"
            )   << "Point " << insidePoint
                << " is not inside the mesh." << nl
                << "Bounding box of the mesh:" << mesh_.bounds()
                << exit(FatalError);
        }

        // Set all cells with this region
        forAll(cellRegion, cellI)
        {
            if (cellRegion[cellI] == keepRegionI)
            {
                if (cellToZone[cellI] == -2)
                {
                    cellToZone[cellI] = surfaceToCellZone[surfI];
                }
                else if (cellToZone[cellI] != surfaceToCellZone[surfI])
                {
                    WarningIn
                    (
                        "meshRefinement::findCellZoneInsideWalk"
                        "(const labelList&, const labelList&"
                        ", const labelList&, const labelList&)"
                    )   << "Cell " << cellI
                        << " at " << mesh_.cellCentres()[cellI]
                        << " is inside surface " << surfaces_.names()[surfI]
                        << " but already marked as being in zone "
                        << cellToZone[cellI] << endl
                        << "This can happen if your surfaces are not"
                        << " (sufficiently) closed."
                        << endl;
                }
            }
        }
    }
}


bool Foam::meshRefinement::calcRegionToZone
(
    const label surfZoneI,
    const label ownRegion,
    const label neiRegion,

    labelList& regionToCellZone
) const
{
    bool changed = false;

    // Check whether inbetween different regions
    if (ownRegion != neiRegion)
    {
        // Jump. Change one of the sides to my type.

        // 1. Interface between my type and unset region.
        // Set region to keepRegion

        if (regionToCellZone[ownRegion] == -2)
        {
            if (regionToCellZone[neiRegion] == surfZoneI)
            {
                // Face between unset and my region. Put unset
                // region into keepRegion
                regionToCellZone[ownRegion] = -1;
                changed = true;
            }
            else if (regionToCellZone[neiRegion] != -2)
            {
                // Face between unset and other region.
                // Put unset region into my region
                regionToCellZone[ownRegion] = surfZoneI;
                changed = true;
            }
        }
        else if (regionToCellZone[neiRegion] == -2)
        {
            if (regionToCellZone[ownRegion] == surfZoneI)
            {
                // Face between unset and my region. Put unset
                // region into keepRegion
                regionToCellZone[neiRegion] = -1;
                changed = true;
            }
            else if (regionToCellZone[ownRegion] != -2)
            {
                // Face between unset and other region.
                // Put unset region into my region
                regionToCellZone[neiRegion] = surfZoneI;
                changed = true;
            }
        }
    }
    return changed;
}


// Finds region per cell. Assumes:
// - region containing keepPoint does not go into a cellZone
// - all other regions can be found by crossing faces marked in
//   namedSurfaceIndex.
void Foam::meshRefinement::findCellZoneTopo
(
    const point& keepPoint,
    const labelList& namedSurfaceIndex,
    const labelList& surfaceToCellZone,
    labelList& cellToZone
) const
{
    // Analyse regions. Reuse regionsplit
    boolList blockedFace(mesh_.nFaces());

    forAll(namedSurfaceIndex, faceI)
    {
        if (namedSurfaceIndex[faceI] == -1)
        {
            blockedFace[faceI] = false;
        }
        else
        {
            blockedFace[faceI] = true;
        }
    }
    // No need to sync since namedSurfaceIndex already is synced

    // Set region per cell based on walking
    regionSplit cellRegion(mesh_, blockedFace);
    blockedFace.clear();

    // Per mesh region the zone the cell should be put in.
    // -2   : not analysed yet
    // -1   : keepPoint region. Not put into any cellzone.
    // >= 0 : index of cellZone
    labelList regionToCellZone(cellRegion.nRegions(), -2);

    // See which cells already are set in the cellToZone (from geometric
    // searching) and use these to take over their zones.
    // Note: could be improved to count number of cells per region.
    forAll(cellToZone, cellI)
    {
        if (cellToZone[cellI] != -2)
        {
            regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
        }
    }


    // Find the region containing the keepPoint
    label keepRegionI = findRegion
    (
        mesh_,
        cellRegion,
        mergeDistance_*vector(1,1,1),
        keepPoint
    );

    Info<< "Found point " << keepPoint
        << " in global region " << keepRegionI
        << " out of " << cellRegion.nRegions() << " regions." << endl;

    if (keepRegionI == -1)
    {
        FatalErrorIn
        (
            "meshRefinement::findCellZoneTopo"
            "(const point&, const labelList&, const labelList&, labelList&)"
        )   << "Point " << keepPoint
            << " is not inside the mesh." << nl
            << "Bounding box of the mesh:" << mesh_.bounds()
            << exit(FatalError);
    }

    // Mark default region with zone -1.
    if (regionToCellZone[keepRegionI] == -2)
    {
        regionToCellZone[keepRegionI] = -1;
    }


    // Find correspondence between cell zone and surface
    // by changing cell zone every time we cross a surface.
    while (true)
    {
        // Synchronise regionToCellZone.
        // Note:
        // - region numbers are identical on all processors
        // - keepRegion is identical ,,
        // - cellZones are identical ,,
        // This done at top of loop to account for geometric matching
        // not being synchronised.
        Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
        Pstream::listCombineScatter(regionToCellZone);


        bool changed = false;

        // Internal faces

        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
        {
            label surfI = namedSurfaceIndex[faceI];

            // Connected even if no cellZone defined for surface
            if (surfI != -1)
            {
                // Calculate region to zone from cellRegions on either side
                // of internal face.
                bool changedCell = calcRegionToZone
                (
                    surfaceToCellZone[surfI],
                    cellRegion[mesh_.faceOwner()[faceI]],
                    cellRegion[mesh_.faceNeighbour()[faceI]],
                    regionToCellZone
                );

                changed = changed | changedCell;
            }
        }

        // Boundary faces

        const polyBoundaryMesh& patches = mesh_.boundaryMesh();

        // Get coupled neighbour cellRegion
        labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
        forAll(patches, patchI)
        {
            const polyPatch& pp = patches[patchI];

            if (pp.coupled())
            {
                forAll(pp, i)
                {
                    label faceI = pp.start()+i;
                    neiCellRegion[faceI-mesh_.nInternalFaces()] =
                        cellRegion[mesh_.faceOwner()[faceI]];
                }
            }
        }
        syncTools::swapBoundaryFaceList(mesh_, neiCellRegion);

        // Calculate region to zone from cellRegions on either side of coupled
        // face.
        forAll(patches, patchI)
        {
            const polyPatch& pp = patches[patchI];

            if (pp.coupled())
            {
                forAll(pp, i)
                {
                    label faceI = pp.start()+i;

                    label surfI = namedSurfaceIndex[faceI];

                    // Connected even if no cellZone defined for surface
                    if (surfI != -1)
                    {
                        bool changedCell = calcRegionToZone
                        (
                            surfaceToCellZone[surfI],
                            cellRegion[mesh_.faceOwner()[faceI]],
                            neiCellRegion[faceI-mesh_.nInternalFaces()],
                            regionToCellZone
                        );

                        changed = changed | changedCell;
                    }
                }
            }
        }

        if (!returnReduce(changed, orOp<bool>()))
        {
            break;
        }
    }


    forAll(regionToCellZone, regionI)
    {
        label zoneI = regionToCellZone[regionI];

        if (zoneI ==  -2)
        {
            FatalErrorIn
            (
                "meshRefinement::findCellZoneTopo"
                "(const point&, const labelList&, const labelList&, labelList&)"
            )   << "For region " << regionI << " haven't set cell zone."
                << exit(FatalError);
        }
    }

    if (debug)
    {
        forAll(regionToCellZone, regionI)
        {
            Pout<< "Region " << regionI
                << " becomes cellZone:" << regionToCellZone[regionI]
                << endl;
        }
    }

    // Rework into cellToZone
    forAll(cellToZone, cellI)
    {
        cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
    }
}


// Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
// faces inbetween same cell zone.
void Foam::meshRefinement::makeConsistentFaceIndex
(
    const labelList& cellToZone,
    labelList& namedSurfaceIndex
) const
{
    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();

    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
    {
        label ownZone = cellToZone[faceOwner[faceI]];
        label neiZone = cellToZone[faceNeighbour[faceI]];

        if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
        {
            namedSurfaceIndex[faceI] = -1;
        }
        else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
        {
            FatalErrorIn("meshRefinement::zonify()")
                << "Different cell zones on either side of face " << faceI
                << " at " << mesh_.faceCentres()[faceI]
                << " but face not marked with a surface."
                << abort(FatalError);
        }
    }

    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    // Get coupled neighbour cellZone
    labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        if (pp.coupled())
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                neiCellZone[faceI-mesh_.nInternalFaces()] =
                    cellToZone[mesh_.faceOwner()[faceI]];
            }
        }
    }
    syncTools::swapBoundaryFaceList(mesh_, neiCellZone);

    // Use coupled cellZone to do check
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        if (pp.coupled())
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;

                label ownZone = cellToZone[faceOwner[faceI]];
                label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];

                if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
                {
                    namedSurfaceIndex[faceI] = -1;
                }
                else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
                {
                    FatalErrorIn("meshRefinement::zonify()")
                        << "Different cell zones on either side of face "
                        << faceI << " at " << mesh_.faceCentres()[faceI]
                        << " but face not marked with a surface."
                        << abort(FatalError);
                }
            }
        }
        else
        {
            // Unzonify boundary faces
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                namedSurfaceIndex[faceI] = -1;
            }
        }
    }
}


void Foam::meshRefinement::handleSnapProblems
(
    const snapParameters& snapParams,
    const bool useTopologicalSnapDetection,
    const bool removeEdgeConnectedCells,
    const scalarField& perpendicularAngle,
    const dictionary& motionDict,
    Time& runTime,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch
)
{
    Info<< nl
        << "Introducing baffles to block off problem cells" << nl
        << "----------------------------------------------" << nl
        << endl;

    labelList facePatch;
    if (useTopologicalSnapDetection)
    {
        facePatch = markFacesOnProblemCells
        (
            motionDict,
            removeEdgeConnectedCells,
            perpendicularAngle,
            globalToMasterPatch
        );
    }
    else
    {
        facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
    }
    Info<< "Analyzed problem cells in = "
        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

    if (debug&MESH)
    {
        faceSet problemFaces(mesh_, "problemFaces", mesh_.nFaces()/100);

        forAll(facePatch, faceI)
        {
            if (facePatch[faceI] != -1)
            {
                problemFaces.insert(faceI);
            }
        }
        problemFaces.instance() = timeName();
        Pout<< "Dumping " << problemFaces.size()
            << " problem faces to " << problemFaces.objectPath() << endl;
        problemFaces.write();
    }

    Info<< "Introducing baffles to delete problem cells." << nl << endl;

    if (debug)
    {
        runTime++;
    }

    // Create baffles with same owner and neighbour for now.
    createBaffles(facePatch, facePatch);

    if (debug)
    {
        // Debug:test all is still synced across proc patches
        checkData();
    }
    Info<< "Created baffles in = "
        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

    printMeshInfo(debug, "After introducing baffles");

    if (debug&MESH)
    {
        Pout<< "Writing extra baffled mesh to time "
            << timeName() << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            runTime.path()/"extraBaffles"
        );
        Pout<< "Dumped debug data in = "
            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
    }
}


Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
(
    const labelList& surfaceToCellZone,
    const labelList& namedSurfaceIndex,
    const labelList& cellToZone,
    const labelList& neiCellZone
) const
{
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();


    DynamicList<label> faceLabels(mesh_.nFaces()/20);

    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
    {
        label surfI = namedSurfaceIndex[faceI];
        if (surfI != -1)
        {
            if
            (
                surfaceToCellZone[surfI] == -1
            ||  max
                (
                    cellToZone[faceOwner[faceI]],
                    cellToZone[faceNeighbour[faceI]]
                ) == -1
            )
            {
                faceLabels.append(faceI);
            }
        }
    }
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        forAll(pp, i)
        {
            label faceI = pp.start()+i;
            label surfI = namedSurfaceIndex[faceI];
            if (surfI != -1)
            {
                if
                (
                    surfaceToCellZone[surfI] == -1
                ||  max
                    (
                        cellToZone[faceOwner[faceI]],
                        neiCellZone[faceI-mesh_.nInternalFaces()]
                    ) == -1
                )
                {
                    faceLabels.append(faceI);
                }
            }
        }
    }
    return faceLabels.shrink();
}


void Foam::meshRefinement::calcPatchNumMasterFaces
(
    const PackedBoolList& isMasterFace,
    const indirectPrimitivePatch& patch,
    labelList& nMasterFacesPerEdge
) const
{
    // Number of (master)faces per edge
    nMasterFacesPerEdge.setSize(patch.nEdges());
    nMasterFacesPerEdge = 0;

    forAll(patch.addressing(), faceI)
    {
        const label meshFaceI = patch.addressing()[faceI];

        if (isMasterFace[meshFaceI])
        {
            const labelList& fEdges = patch.faceEdges()[faceI];
            forAll(fEdges, fEdgeI)
            {
                nMasterFacesPerEdge[fEdges[fEdgeI]]++;
            }
        }
    }

    syncTools::syncEdgeList
    (
        mesh_,
        patch.meshEdges(mesh_.edges(), mesh_.pointEdges()),
        nMasterFacesPerEdge,
        plusEqOp<label>(),
        0
    );
}


Foam::label Foam::meshRefinement::markPatchZones
(
    const indirectPrimitivePatch& patch,
    const labelList& nMasterFacesPerEdge,
    labelList& faceToZone
) const
{
    List<patchEdgeFaceRegion> allEdgeInfo(patch.nEdges());
    List<patchEdgeFaceRegion> allFaceInfo(patch.size());


    // Protect all non-manifold edges
    {
        label nProtected = 0;

        forAll(nMasterFacesPerEdge, edgeI)
        {
            if (nMasterFacesPerEdge[edgeI] > 2)
            {
                allEdgeInfo[edgeI] = -2;
                nProtected++;
            }
        }
        //Info<< "Protected from visiting "
        //    << returnReduce(nProtected, sumOp<label>())
        //    << " non-manifold edges" << nl << endl;
    }


    // Hand out zones

    DynamicList<label> changedEdges;
    DynamicList<patchEdgeFaceRegion> changedInfo;

    const scalar tol = PatchEdgeFaceWave
    <
        indirectPrimitivePatch,
        patchEdgeFaceRegion
    >::propagationTol();

    int dummyTrackData;

    const globalIndex globalFaces(patch.size());

    label faceI = 0;

    label currentZoneI = 0;

    while (true)
    {
        // Pick an unset face
        label globalSeed = labelMax;
        for (; faceI < allFaceInfo.size(); faceI++)
        {
            if (!allFaceInfo[faceI].valid(dummyTrackData))
            {
                globalSeed = globalFaces.toGlobal(faceI);
                break;
            }
        }

        reduce(globalSeed, minOp<label>());

        if (globalSeed == labelMax)
        {
            break;
        }

        label procI = globalFaces.whichProcID(globalSeed);
        label seedFaceI = globalFaces.toLocal(procI, globalSeed);

        //Info<< "Seeding zone " << currentZoneI
        //    << " from processor " << procI << " face " << seedFaceI
        //    << endl;

        if (procI == Pstream::myProcNo())
        {
            patchEdgeFaceRegion& faceInfo = allFaceInfo[seedFaceI];


            // Set face
            faceInfo = currentZoneI;

            // .. and seed its edges
            const labelList& fEdges = patch.faceEdges()[seedFaceI];
            forAll(fEdges, fEdgeI)
            {
                label edgeI = fEdges[fEdgeI];

                patchEdgeFaceRegion& edgeInfo = allEdgeInfo[edgeI];

                if
                (
                    edgeInfo.updateEdge<int>
                    (
                        mesh_,
                        patch,
                        edgeI,
                        seedFaceI,
                        faceInfo,
                        tol,
                        dummyTrackData
                    )
                )
                {
                    changedEdges.append(edgeI);
                    changedInfo.append(edgeInfo);
                }
            }
        }


        if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
        {
            break;
        }


        // Walk
        PatchEdgeFaceWave
        <
            indirectPrimitivePatch,
            patchEdgeFaceRegion
        > calc
        (
            mesh_,
            patch,
            changedEdges,
            changedInfo,
            allEdgeInfo,
            allFaceInfo,
            returnReduce(patch.nEdges(), sumOp<label>())
        );

        currentZoneI++;
    }


    faceToZone.setSize(patch.size());
    forAll(allFaceInfo, faceI)
    {
        if (!allFaceInfo[faceI].valid(dummyTrackData))
        {
            FatalErrorIn("meshRefinement::markPatchZones(..)")
                << "Problem: unvisited face " << faceI
                << " at " << patch.faceCentres()[faceI]
                << exit(FatalError);
        }
        faceToZone[faceI] = allFaceInfo[faceI].region();
    }

    return currentZoneI;
}


void Foam::meshRefinement::consistentOrientation
(
    const PackedBoolList& isMasterFace,
    const indirectPrimitivePatch& patch,
    const labelList& nMasterFacesPerEdge,
    const labelList& faceToZone,
    const Map<label>& zoneToOrientation,
    PackedBoolList& meshFlipMap
) const
{
    const polyBoundaryMesh& bm = mesh_.boundaryMesh();

    // Data on all edges and faces
    List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
    List<patchFaceOrientation> allFaceInfo(patch.size());

    // Make sure we don't walk through
    // - slaves of coupled faces
    // - non-manifold edges
    {
        label nProtected = 0;

        forAll(patch.addressing(), faceI)
        {
            const label meshFaceI = patch.addressing()[faceI];
            const label patchI = bm.whichPatch(meshFaceI);

            if
            (
                patchI != -1
             && bm[patchI].coupled()
             && !isMasterFace[meshFaceI]
            )
            {
                // Slave side. Mark so doesn't get visited.
                allFaceInfo[faceI] = orientedSurface::NOFLIP;
                nProtected++;
            }
        }
        //Info<< "Protected from visiting "
        //    << returnReduce(nProtected, sumOp<label>())
        //    << " slaves of coupled faces" << nl << endl;
    }
    {
        label nProtected = 0;

        forAll(nMasterFacesPerEdge, edgeI)
        {
            if (nMasterFacesPerEdge[edgeI] > 2)
            {
                allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
                nProtected++;
            }
        }

        Info<< "Protected from visiting "
            << returnReduce(nProtected, sumOp<label>())
            << " non-manifold edges" << nl << endl;
    }



    DynamicList<label> changedEdges;
    DynamicList<patchFaceOrientation> changedInfo;

    const scalar tol = PatchEdgeFaceWave
    <
        indirectPrimitivePatch,
        patchFaceOrientation
    >::propagationTol();

    int dummyTrackData;

    globalIndex globalFaces(patch.size());

    while (true)
    {
        // Pick an unset face
        label globalSeed = labelMax;
        forAll(allFaceInfo, faceI)
        {
            if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
            {
                globalSeed = globalFaces.toGlobal(faceI);
                break;
            }
        }

        reduce(globalSeed, minOp<label>());

        if (globalSeed == labelMax)
        {
            break;
        }

        label procI = globalFaces.whichProcID(globalSeed);
        label seedFaceI = globalFaces.toLocal(procI, globalSeed);

        //Info<< "Seeding from processor " << procI << " face " << seedFaceI
        //    << endl;

        if (procI == Pstream::myProcNo())
        {
            // Determine orientation of seedFace

            patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];

            // Start off with correct orientation
            faceInfo = orientedSurface::NOFLIP;

            if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
            {
                faceInfo.flip();
            }


            const labelList& fEdges = patch.faceEdges()[seedFaceI];
            forAll(fEdges, fEdgeI)
            {
                label edgeI = fEdges[fEdgeI];

                patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];

                if
                (
                    edgeInfo.updateEdge<int>
                    (
                        mesh_,
                        patch,
                        edgeI,
                        seedFaceI,
                        faceInfo,
                        tol,
                        dummyTrackData
                    )
                )
                {
                    changedEdges.append(edgeI);
                    changedInfo.append(edgeInfo);
                }
            }
        }


        if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
        {
            break;
        }



        // Walk
        PatchEdgeFaceWave
        <
            indirectPrimitivePatch,
            patchFaceOrientation
        > calc
        (
            mesh_,
            patch,
            changedEdges,
            changedInfo,
            allEdgeInfo,
            allFaceInfo,
            returnReduce(patch.nEdges(), sumOp<label>())
        );
    }


    // Push master zone info over to slave (since slave faces never visited)
    {
        labelList neiStatus
        (
            mesh_.nFaces()-mesh_.nInternalFaces(),
            orientedSurface::UNVISITED
        );

        forAll(patch.addressing(), i)
        {
            const label meshFaceI = patch.addressing()[i];
            if (!mesh_.isInternalFace(meshFaceI))
            {
                neiStatus[meshFaceI-mesh_.nInternalFaces()] =
                    allFaceInfo[i].flipStatus();
            }
        }
        syncTools::swapBoundaryFaceList(mesh_, neiStatus);

        forAll(patch.addressing(), i)
        {
            const label meshFaceI = patch.addressing()[i];
            const label patchI = bm.whichPatch(meshFaceI);

            if
            (
                patchI != -1
             && bm[patchI].coupled()
             && !isMasterFace[meshFaceI]
            )
            {
                // Slave side. Take flipped from neighbour
                label bFaceI = meshFaceI-mesh_.nInternalFaces();

                if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
                {
                    allFaceInfo[i] = orientedSurface::FLIP;
                }
                else if (neiStatus[bFaceI] == orientedSurface::FLIP)
                {
                    allFaceInfo[i] = orientedSurface::NOFLIP;
                }
                else
                {
                    FatalErrorIn("meshRefinement::consistentOrientation(..)")
                        << "Incorrect status for face " << meshFaceI
                        << abort(FatalError);
                }
            }
        }
    }


    // Convert to meshFlipMap and adapt faceZones

    meshFlipMap.setSize(mesh_.nFaces());
    meshFlipMap = false;

    forAll(allFaceInfo, faceI)
    {
        label meshFaceI = patch.addressing()[faceI];

        if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
        {
            meshFlipMap[meshFaceI] = false;
        }
        else if (allFaceInfo[faceI] == orientedSurface::FLIP)
        {
            meshFlipMap[meshFaceI] = true;
        }
        else
        {
            FatalErrorIn("meshRefinement::consistentOrientation(..)")
                << "Problem : unvisited face " << faceI
                << " centre:" << mesh_.faceCentres()[meshFaceI]
                << abort(FatalError);
        }
    }
}


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

// Split off unreachable areas of mesh.
void Foam::meshRefinement::baffleAndSplitMesh
(
    const bool doHandleSnapProblems,
    const snapParameters& snapParams,
    const bool useTopologicalSnapDetection,
    const bool removeEdgeConnectedCells,
    const scalarField& perpendicularAngle,
    const bool mergeFreeStanding,
    const scalar planarAngle,
    const dictionary& motionDict,
    Time& runTime,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch,
    const point& keepPoint
)
{
    // Introduce baffles
    // ~~~~~~~~~~~~~~~~~

    // Split the mesh along internal faces wherever there is a pierce between
    // two cell centres.

    Info<< "Introducing baffles for "
        << returnReduce(countHits(), sumOp<label>())
        << " faces that are intersected by the surface." << nl << endl;

    // Swap neighbouring cell centres and cell level
    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
    calcNeighbourData(neiLevel, neiCc);

    labelList ownPatch, neiPatch;
    getBafflePatches
    (
        globalToMasterPatch,
        neiLevel,
        neiCc,

        ownPatch,
        neiPatch
    );

    createBaffles(ownPatch, neiPatch);

    if (debug)
    {
        // Debug:test all is still synced across proc patches
        checkData();
    }

    Info<< "Created baffles in = "
        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

    printMeshInfo(debug, "After introducing baffles");

    if (debug&MESH)
    {
        Pout<< "Writing baffled mesh to time " << timeName()
            << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            runTime.path()/"baffles"
        );
        Pout<< "Dumped debug data in = "
            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
    }


    // Introduce baffles to delete problem cells
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Create some additional baffles where we want surface cells removed.

    if (doHandleSnapProblems)
    {
        handleSnapProblems
        (
            snapParams,
            useTopologicalSnapDetection,
            removeEdgeConnectedCells,
            perpendicularAngle,
            motionDict,
            runTime,
            globalToMasterPatch,
            globalToSlavePatch
        );
    }


    // Select part of mesh
    // ~~~~~~~~~~~~~~~~~~~

    Info<< nl
        << "Remove unreachable sections of mesh" << nl
        << "-----------------------------------" << nl
        << endl;

    if (debug)
    {
        runTime++;
    }

    splitMeshRegions(globalToMasterPatch, globalToSlavePatch, keepPoint);

    if (debug)
    {
        // Debug:test all is still synced across proc patches
        checkData();
    }
    Info<< "Split mesh in = "
        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

    printMeshInfo(debug, "After subsetting");

    if (debug&MESH)
    {
        Pout<< "Writing subsetted mesh to time " << timeName()
            << endl;
        write
        (
            debugType(debug),
            writeType(writeLevel() | WRITEMESH),
            runTime.path()/timeName()
        );
        Pout<< "Dumped debug data in = "
            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
    }


    // Merge baffles
    // ~~~~~~~~~~~~~

    if (mergeFreeStanding)
    {
        Info<< nl
            << "Merge free-standing baffles" << nl
            << "---------------------------" << nl
            << endl;


        // List of pairs of freestanding baffle faces.
        List<labelPair> couples
        (
            freeStandingBaffles    // filter out freestanding baffles
            (
                localPointRegion::findDuplicateFacePairs(mesh_),
                planarAngle
            )
        );

        label nCouples = couples.size();
        reduce(nCouples, sumOp<label>());

        Info<< "Detected free-standing baffles : " << nCouples << endl;

        if (nCouples > 0)
        {
            // Actually merge baffles. Note: not exactly parallellized. Should
            // convert baffle faces into processor faces if they resulted
            // from them.
            mergeBaffles(couples);

            // Detect any problem cells resulting from merging of baffles
            // and delete them
            handleSnapProblems
            (
                snapParams,
                useTopologicalSnapDetection,
                removeEdgeConnectedCells,
                perpendicularAngle,
                motionDict,
                runTime,
                globalToMasterPatch,
                globalToSlavePatch
            );

            if (debug)
            {
                // Debug:test all is still synced across proc patches
                checkData();
            }
        }
        Info<< "Merged free-standing baffles in = "
            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
    }
}


// Split off (with optional buffer layers) unreachable areas of mesh.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
(
    const label nBufferLayers,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch,
    const point& keepPoint
)
{
    // Determine patches to put intersections into
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Swap neighbouring cell centres and cell level
    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
    calcNeighbourData(neiLevel, neiCc);

    labelList ownPatch, neiPatch;
    getBafflePatches
    (
        globalToMasterPatch,
        neiLevel,
        neiCc,

        ownPatch,
        neiPatch
    );

    // Analyse regions. Reuse regionsplit
    boolList blockedFace(mesh_.nFaces(), false);

    forAll(ownPatch, faceI)
    {
        if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
        {
            blockedFace[faceI] = true;
        }
    }
    syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>());

    // Set region per cell based on walking
    regionSplit cellRegion(mesh_, blockedFace);
    blockedFace.clear();

    // Find the region containing the keepPoint
    const label keepRegionI = findRegion
    (
        mesh_,
        cellRegion,
        mergeDistance_*vector(1,1,1),
        keepPoint
    );

    Info<< "Found point " << keepPoint
        << " in global region " << keepRegionI
        << " out of " << cellRegion.nRegions() << " regions." << endl;

    if (keepRegionI == -1)
    {
        FatalErrorIn
        (
            "meshRefinement::splitMesh"
            "(const label, const labelList&, const point&)"
        )   << "Point " << keepPoint
            << " is not inside the mesh." << nl
            << "Bounding box of the mesh:" << mesh_.bounds()
            << exit(FatalError);
    }


    // Walk out nBufferlayers from region split
    // (modifies cellRegion, ownPatch)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Takes over face patch onto points and then back to faces and cells
    // (so cell-face-point walk)

    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();

    // Patch for exposed faces for lack of anything sensible.
    label defaultPatch = 0;
    if (globalToMasterPatch.size())
    {
        defaultPatch = globalToMasterPatch[0];
    }

    for (label i = 0; i < nBufferLayers; i++)
    {
        // 1. From cells (via faces) to points

        labelList pointBaffle(mesh_.nPoints(), -1);

        forAll(faceNeighbour, faceI)
        {
            const face& f = mesh_.faces()[faceI];

            label ownRegion = cellRegion[faceOwner[faceI]];
            label neiRegion = cellRegion[faceNeighbour[faceI]];

            if (ownRegion == keepRegionI && neiRegion != keepRegionI)
            {
                // Note max(..) since possibly regionSplit might have split
                // off extra unreachable parts of mesh. Note: or can this only
                // happen for boundary faces?
                forAll(f, fp)
                {
                    pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
                }
            }
            else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
            {
                label newPatchI = neiPatch[faceI];
                if (newPatchI == -1)
                {
                    newPatchI = max(defaultPatch, ownPatch[faceI]);
                }
                forAll(f, fp)
                {
                    pointBaffle[f[fp]] = newPatchI;
                }
            }
        }
        for
        (
            label faceI = mesh_.nInternalFaces();
            faceI < mesh_.nFaces();
            faceI++
        )
        {
            const face& f = mesh_.faces()[faceI];

            label ownRegion = cellRegion[faceOwner[faceI]];

            if (ownRegion == keepRegionI)
            {
                forAll(f, fp)
                {
                    pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
                }
            }
        }

        // Sync
        syncTools::syncPointList
        (
            mesh_,
            pointBaffle,
            maxEqOp<label>(),
            label(-1)           // null value
        );


        // 2. From points back to faces

        const labelListList& pointFaces = mesh_.pointFaces();

        forAll(pointFaces, pointI)
        {
            if (pointBaffle[pointI] != -1)
            {
                const labelList& pFaces = pointFaces[pointI];

                forAll(pFaces, pFaceI)
                {
                    label faceI = pFaces[pFaceI];

                    if (ownPatch[faceI] == -1)
                    {
                        ownPatch[faceI] = pointBaffle[pointI];
                    }
                }
            }
        }
        syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());


        // 3. From faces to cells (cellRegion) and back to faces (ownPatch)

        labelList newOwnPatch(ownPatch);

        forAll(ownPatch, faceI)
        {
            if (ownPatch[faceI] != -1)
            {
                label own = faceOwner[faceI];

                if (cellRegion[own] != keepRegionI)
                {
                    cellRegion[own] = keepRegionI;

                    const cell& ownFaces = mesh_.cells()[own];
                    forAll(ownFaces, j)
                    {
                        if (ownPatch[ownFaces[j]] == -1)
                        {
                            newOwnPatch[ownFaces[j]] = ownPatch[faceI];
                        }
                    }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = faceNeighbour[faceI];

                    if (cellRegion[nei] != keepRegionI)
                    {
                        cellRegion[nei] = keepRegionI;

                        const cell& neiFaces = mesh_.cells()[nei];
                        forAll(neiFaces, j)
                        {
                            if (ownPatch[neiFaces[j]] == -1)
                            {
                                newOwnPatch[neiFaces[j]] = ownPatch[faceI];
                            }
                        }
                    }
                }
            }
        }

        ownPatch.transfer(newOwnPatch);

        syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
    }



    // Subset
    // ~~~~~~

    // Get cells to remove
    DynamicList<label> cellsToRemove(mesh_.nCells());
    forAll(cellRegion, cellI)
    {
        if (cellRegion[cellI] != keepRegionI)
        {
            cellsToRemove.append(cellI);
        }
    }
    cellsToRemove.shrink();

    label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
    reduce(nCellsToKeep, sumOp<label>());

    Info<< "Keeping all cells in region " << keepRegionI
        << " containing point " << keepPoint << endl
        << "Selected for keeping : " << nCellsToKeep
        << " cells." << endl;


    // Remove cells
    removeCells cellRemover(mesh_);

    // Pick up patches for exposed faces
    labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
    labelList exposedPatches(exposedFaces.size());

    forAll(exposedFaces, i)
    {
        label faceI = exposedFaces[i];

        if (ownPatch[faceI] != -1)
        {
            exposedPatches[i] = ownPatch[faceI];
        }
        else
        {
            WarningIn("meshRefinement::splitMesh(..)")
                << "For exposed face " << faceI
                << " fc:" << mesh_.faceCentres()[faceI]
                << " found no patch." << endl
                << "    Taking patch " << defaultPatch
                << " instead." << endl;
            exposedPatches[i] = defaultPatch;
        }
    }

    return doRemoveCells
    (
        cellsToRemove,
        exposedFaces,
        exposedPatches,
        cellRemover
    );
}


// Find boundary points that connect to more than one cell region and
// split them.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
(
    const localPointRegion& regionSide
)
{
    // Topochange container
    polyTopoChange meshMod(mesh_);

    label nNonManifPoints = returnReduce
    (
        regionSide.meshPointMap().size(),
        sumOp<label>()
    );

    Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
        << " non-manifold points (out of "
        << mesh_.globalData().nTotalPoints()
        << ')' << endl;

    // Topo change engine
    duplicatePoints pointDuplicator(mesh_);

    // Insert changes into meshMod
    pointDuplicator.setRefinement(regionSide, meshMod);

    // Change the mesh (no inflation, parallel sync)
    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);

    // Update fields
    mesh_.updateMesh(map);

    // Move mesh if in inflation mode
    if (map().hasMotionPoints())
    {
        mesh_.movePoints(map().preMotionPoints());
    }
    else
    {
        // Delete mesh volumes.
        mesh_.clearOut();
    }

    // Reset the instance for if in overwrite mode
    mesh_.setInstance(timeName());

    // Update intersections. Is mapping only (no faces created, positions stay
    // same) so no need to recalculate intersections.
    updateMesh(map, labelList(0));

    return map;
}


// Find boundary points that connect to more than one cell region and
// split them.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
{
    // Analyse which points need to be duplicated
    localPointRegion regionSide(mesh_);

    return dupNonManifoldPoints(regionSide);
}


// Zoning
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
(
    const point& keepPoint,
    const bool allowFreeStandingZoneFaces
)
{
    const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();

    labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));

    forAll(namedSurfaces, i)
    {
        label surfI = namedSurfaces[i];

        Info<< "Surface : " << surfaces_.names()[surfI] << nl
            << "    faceZone : " << surfZones[surfI].faceZoneName() << nl
            << "    cellZone : " << surfZones[surfI].cellZoneName() << endl;
    }


    // Add zones to mesh
    labelList surfaceToFaceZone =
        surfaceZonesInfo::addFaceZonesToMesh
        (
            surfZones,
            namedSurfaces,
            mesh_
        );

    labelList surfaceToCellZone =
        surfaceZonesInfo::addCellZonesToMesh
        (
            surfZones,
            namedSurfaces,
            mesh_
        );


    const pointField& cellCentres = mesh_.cellCentres();
    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();


    // Swap neighbouring cell centres and cell level
    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
    calcNeighbourData(neiLevel, neiCc);


    // Mark faces intersecting zoned surfaces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


    // Like surfaceIndex_ but only for named surfaces.
    labelList namedSurfaceIndex(mesh_.nFaces(), -1);
    PackedBoolList posOrientation(mesh_.nFaces());

    {
        // Statistics: number of faces per faceZone
        labelList nSurfFaces(surfZones.size(), 0);

        // Note: for all internal faces? internal + coupled?
        // Since zonify is run after baffling the surfaceIndex_ on baffles is
        // not synchronised across both baffle faces. Fortunately we don't
        // do zonify baffle faces anyway (they are normal boundary faces).

        // Collect candidate faces
        // ~~~~~~~~~~~~~~~~~~~~~~~

        labelList testFaces(intersectedFaces());

        // Collect segments
        // ~~~~~~~~~~~~~~~~

        pointField start(testFaces.size());
        pointField end(testFaces.size());

        forAll(testFaces, i)
        {
            label faceI = testFaces[i];

            if (mesh_.isInternalFace(faceI))
            {
                start[i] = cellCentres[faceOwner[faceI]];
                end[i] = cellCentres[faceNeighbour[faceI]];
            }
            else
            {
                start[i] = cellCentres[faceOwner[faceI]];
                end[i] = neiCc[faceI-mesh_.nInternalFaces()];
            }
        }

        // Extend segments a bit
        {
            const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
            start -= smallVec;
            end += smallVec;
        }


        // Do test for intersections
        // ~~~~~~~~~~~~~~~~~~~~~~~~~
        // Note that we intersect all intersected faces again. Could reuse
        // the information already in surfaceIndex_.

        labelList surface1;
        List<pointIndexHit> hit1;
        vectorField normal1;
        labelList surface2;
        List<pointIndexHit> hit2;
        vectorField normal2;
        {
            labelList region1;
            labelList region2;
            surfaces_.findNearestIntersection
            (
                namedSurfaces,
                start,
                end,

                surface1,
                hit1,
                region1,
                normal1,

                surface2,
                hit2,
                region2,
                normal2
            );
        }

        forAll(testFaces, i)
        {
            label faceI = testFaces[i];
            const vector& area = mesh_.faceAreas()[faceI];

            if (surface1[i] != -1)
            {
                // If both hit should probably choose 'nearest'
                if
                (
                    surface2[i] != -1
                 && (
                        magSqr(hit2[i].hitPoint())
                      < magSqr(hit1[i].hitPoint())
                    )
                )
                {
                    namedSurfaceIndex[faceI] = surface2[i];
                    posOrientation[faceI] = ((area&normal2[i]) > 0);
                    nSurfFaces[surface2[i]]++;
                }
                else
                {
                    namedSurfaceIndex[faceI] = surface1[i];
                    posOrientation[faceI] = ((area&normal1[i]) > 0);
                    nSurfFaces[surface1[i]]++;
                }
            }
            else if (surface2[i] != -1)
            {
                namedSurfaceIndex[faceI] = surface2[i];
                posOrientation[faceI] = ((area&normal2[i]) > 0);
                nSurfFaces[surface2[i]]++;
            }
        }


        // surfaceIndex might have different surfaces on both sides if
        // there happen to be a (obviously thin) surface with different
        // regions between the cell centres. If one is on a named surface
        // and the other is not this might give problems so sync.
        syncTools::syncFaceList
        (
            mesh_,
            namedSurfaceIndex,
            maxEqOp<label>()
        );

        // Print a bit
        if (debug)
        {
            forAll(nSurfFaces, surfI)
            {
                Pout<< "Surface:"
                    << surfaces_.names()[surfI]
                    << "  nZoneFaces:" << nSurfFaces[surfI] << nl;
            }
            Pout<< endl;
        }
    }


    // Put the cells into the correct zone
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Zone per cell:
    // -2 : unset
    // -1 : not in any zone
    // >=0: zoneID
    labelList cellToZone(mesh_.nCells(), -2);


    // Set using geometric test
    // ~~~~~~~~~~~~~~~~~~~~~~~~

    // Closed surfaces with cellZone specified.
    labelList closedNamedSurfaces
    (
        surfaceZonesInfo::getClosedNamedSurfaces
        (
            surfZones,
            surfaces_.geometry(),
            surfaces_.surfaces()
        )
    );

    if (closedNamedSurfaces.size())
    {
        Info<< "Found " << closedNamedSurfaces.size()
            << " closed, named surfaces. Assigning cells in/outside"
            << " these surfaces to the corresponding cellZone."
            << nl << endl;

        findCellZoneGeometric
        (
            neiCc,
            closedNamedSurfaces,    // indices of closed surfaces
            namedSurfaceIndex,      // per face index of named surface
            surfaceToCellZone,      // cell zone index per surface

            cellToZone
        );
    }


    // Set using provided locations
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    labelList locationSurfaces
    (
        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
    );

    if (locationSurfaces.size())
    {
        Info<< "Found " << locationSurfaces.size()
            << " named surfaces with a provided inside point."
            << " Assigning cells inside these surfaces"
            << " to the corresponding cellZone."
            << nl << endl;

        findCellZoneInsideWalk
        (
            locationSurfaces,       // indices of closed surfaces
            namedSurfaceIndex,      // per face index of named surface
            surfaceToCellZone,      // cell zone index per surface

            cellToZone
        );
    }


    // Set using walking
    // ~~~~~~~~~~~~~~~~~

    {
        Info<< "Walking from location-in-mesh " << keepPoint
            << " to assign cellZones "
            << "- crossing a faceZone face changes cellZone" << nl << endl;

        // Topological walk
        findCellZoneTopo
        (
            keepPoint,
            namedSurfaceIndex,
            surfaceToCellZone,

            cellToZone
        );
    }


    // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
    if (!allowFreeStandingZoneFaces)
    {
        Info<< "Only keeping zone faces inbetween different cellZones."
            << nl << endl;

        makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
    }

    // Topochange container
    polyTopoChange meshMod(mesh_);



    // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
    labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        if (pp.coupled())
        {
            forAll(pp, i)
            {
                label faceI = pp.start()+i;
                neiCellZone[faceI-mesh_.nInternalFaces()] =
                    cellToZone[mesh_.faceOwner()[faceI]];
            }
        }
    }
    syncTools::swapBoundaryFaceList(mesh_, neiCellZone);

    // Get per face whether is it master (of a coupled set of faces)
    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));



    // faceZones
    // ~~~~~~~~~
    // Faces on faceZones come in two variants:
    // - faces on the outside of a cellZone. They will be oriented to
    //   point out of the maximum cellZone.
    // - free-standing faces. These will be oriented according to the
    //   local surface normal. We do this in a two step algorithm:
    //      - do a consistent orientation
    //      - check number of faces with consistent orientation
    //      - if <0 flip the whole patch
    PackedBoolList meshFlipMap(mesh_.nFaces(), false);
    {
        // Collect all data on zone faces without cellZones on either side.
        const indirectPrimitivePatch patch
        (
            IndirectList<face>
            (
                mesh_.faces(),
                freeStandingBaffleFaces
                (
                    surfaceToCellZone,
                    namedSurfaceIndex,
                    cellToZone,
                    neiCellZone
                )
            ),
            mesh_.points()
        );

        label nFreeStanding = returnReduce(patch.size(), sumOp<label>());
        if (nFreeStanding > 0)
        {
            Info<< "Detected " << nFreeStanding << " free-standing zone faces"
                << endl;

            // Detect non-manifold edges
            labelList nMasterFacesPerEdge;
            calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);

            // Mark zones. Even a single original surface might create multiple
            // disconnected/non-manifold-connected zones
            labelList faceToZone;
            const label nZones = markPatchZones
            (
                patch,
                nMasterFacesPerEdge,
                faceToZone
            );

            Map<label> nPosOrientation(2*nZones);
            for (label zoneI = 0; zoneI < nZones; zoneI++)
            {
                nPosOrientation.insert(zoneI, 0);
            }

            // Make orientations consistent in a topological way. This just
            // checks the first face per zone for whether nPosOrientation
            // is negative (which it never is at this point)
            consistentOrientation
            (
                isMasterFace,
                patch,
                nMasterFacesPerEdge,
                faceToZone,
                nPosOrientation,

                meshFlipMap
            );

            // Count per region the number of orientations (taking the new
            // flipMap into account)
            forAll(patch.addressing(), faceI)
            {
                label meshFaceI = patch.addressing()[faceI];

                if (isMasterFace[meshFaceI])
                {
                    label n = 1;
                    if
                    (
                        bool(posOrientation[meshFaceI])
                     == meshFlipMap[meshFaceI]
                    )
                    {
                        n = -1;
                    }

                    nPosOrientation.find(faceToZone[faceI])() += n;
                }
            }
            Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
            Pstream::mapCombineScatter(nPosOrientation);


            Info<< "Split " << nFreeStanding << " free-standing zone faces"
                << " into " << nZones << " disconnected regions with size"
                << " (negative denotes wrong orientation) :"
                << endl;

            for (label zoneI = 0; zoneI < nZones; zoneI++)
            {
                Info<< "    " << zoneI << "\t" << nPosOrientation[zoneI]
                    << endl;
            }
            Info<< endl;


            // Reapply with new counts (in nPosOrientation). This will cause
            // zones with a negative count to be flipped.
            consistentOrientation
            (
                isMasterFace,
                patch,
                nMasterFacesPerEdge,
                faceToZone,
                nPosOrientation,

                meshFlipMap
            );
        }
    }


    // Put the faces into the correct zone
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
    {
        label surfI = namedSurfaceIndex[faceI];

        if (surfI != -1)
        {
            bool flip;

            if (surfaceToCellZone[surfI] == -1)
            {
                // Surface with faceZone only (= freestanding faceZone).
                // Use geometrically derived orientation
                flip = meshFlipMap[faceI];
            }
            else
            {
                // Orient face zone to have slave cells in max cell zone.
                label ownZone = cellToZone[faceOwner[faceI]];
                label neiZone = cellToZone[faceNeighbour[faceI]];
                label maxZone = max(ownZone, neiZone);

                if (maxZone == -1)
                {
                    // free-standing face. Use geometrically derived orientation
                    flip = meshFlipMap[faceI];
                }
                else if (ownZone == maxZone)
                {
                    flip = false;
                }
                else
                {
                    flip = true;
                }
            }

            meshMod.setAction
            (
                polyModifyFace
                (
                    mesh_.faces()[faceI],           // modified face
                    faceI,                          // label of face
                    faceOwner[faceI],               // owner
                    faceNeighbour[faceI],           // neighbour
                    false,                          // face flip
                    -1,                             // patch for face
                    false,                          // remove from zone
                    surfaceToFaceZone[surfI],       // zone for face
                    flip                            // face flip in zone
                )
            );
        }
    }


    // Set owner as no-flip
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];

        label faceI = pp.start();

        forAll(pp, i)
        {
            label surfI = namedSurfaceIndex[faceI];

            if (surfI != -1)
            {
                bool flip;

                if (surfaceToCellZone[surfI] == -1)
                {
                    // Surface with faceZone only (= freestanding faceZone).
                    // Use geometrically derived orientation
                    flip = meshFlipMap[faceI];
                }
                else
                {
                    label ownZone = cellToZone[faceOwner[faceI]];
                    label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
                    label maxZone = max(ownZone, neiZone);

                    if (maxZone == -1)
                    {
                        // Surface with faceZone only (= freestanding faceZone)
                        // Use geometrically derived orientation
                        flip = meshFlipMap[faceI];
                    }
                    else if (ownZone == neiZone)
                    {
                        // Free-standing zone face or coupled boundary.
                        // Keep masterface unflipped.
                        flip = !isMasterFace[faceI];
                    }
                    else if (neiZone == maxZone)
                    {
                        flip = true;
                    }
                    else
                    {
                        flip = false;
                    }
                }

                meshMod.setAction
                (
                    polyModifyFace
                    (
                        mesh_.faces()[faceI],           // modified face
                        faceI,                          // label of face
                        faceOwner[faceI],               // owner
                        -1,                             // neighbour
                        false,                          // face flip
                        patchI,                         // patch for face
                        false,                          // remove from zone
                        surfaceToFaceZone[surfI],       // zone for face
                        flip                            // face flip in zone
                    )
                );
            }
            faceI++;
        }
    }


    // Put the cells into the correct zone
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    forAll(cellToZone, cellI)
    {
        label zoneI = cellToZone[cellI];

        if (zoneI >= 0)
        {
            meshMod.setAction
            (
                polyModifyCell
                (
                    cellI,
                    false,          // removeFromZone
                    zoneI
                )
            );
        }
    }

    // Change the mesh (no inflation, parallel sync)
    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);

    // Update fields
    mesh_.updateMesh(map);

    // Move mesh if in inflation mode
    if (map().hasMotionPoints())
    {
        mesh_.movePoints(map().preMotionPoints());
    }
    else
    {
        // Delete mesh volumes.
        mesh_.clearOut();
    }

    // Reset the instance for if in overwrite mode
    mesh_.setInstance(timeName());

    // Print some stats (note: zones are synchronised)
    if (mesh_.cellZones().size() > 0)
    {
        Info<< "CellZones:" << endl;
        forAll(mesh_.cellZones(), zoneI)
        {
            const cellZone& cz = mesh_.cellZones()[zoneI];
            Info<< "    " << cz.name()
                << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
                << endl;
        }
        Info<< endl;
    }
    if (mesh_.faceZones().size() > 0)
    {
        Info<< "FaceZones:" << endl;
        forAll(mesh_.faceZones(), zoneI)
        {
            const faceZone& fz = mesh_.faceZones()[zoneI];
            Info<< "    " << fz.name()
                << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
                << endl;
        }
        Info<< endl;
    }

    // None of the faces has changed, only the zones. Still...
    updateMesh(map, labelList());

    return map;
}


// ************************************************************************* //
meshRefinementBaffles.C (100,964 bytes)   

user4

2015-04-10 15:15

 

meshRefinement.H (37,347 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     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/>.

Class
    Foam::meshRefinement

Description
    Helper class which maintains intersections of (changing) mesh with
    (static) surfaces.

    Maintains
    - per face any intersections of the cc-cc segment with any of the surfaces

SourceFiles
    meshRefinement.C
    meshRefinementBaffles.C
    meshRefinementMerge.C
    meshRefinementProblemCells.C
    meshRefinementRefine.C

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

#ifndef meshRefinement_H
#define meshRefinement_H

#include "hexRef8.H"
#include "mapPolyMesh.H"
#include "autoPtr.H"
#include "labelPair.H"
#include "indirectPrimitivePatch.H"
#include "pointFieldsFwd.H"
#include "Tuple2.H"
#include "pointIndexHit.H"

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

namespace Foam
{

// Class forward declarations
class fvMesh;
class mapDistributePolyMesh;
class decompositionMethod;
class refinementSurfaces;
class refinementFeatures;
class shellSurfaces;
class removeCells;
class fvMeshDistribute;
class searchableSurface;
class regionSplit;
class globalIndex;
class removePoints;
class localPointRegion;

class snapParameters;

/*---------------------------------------------------------------------------*\
                           Class meshRefinement Declaration
\*---------------------------------------------------------------------------*/

class meshRefinement
{

public:

    // Public data types

        //- Enumeration for what to debug
        enum IOdebugType
        {
            IOMESH,
            //IOSCALARLEVELS,
            IOOBJINTERSECTIONS,
            IOFEATURESEEDS,
            IOATTRACTION,
            IOLAYERINFO
        };
        static const NamedEnum<IOdebugType, 5> IOdebugTypeNames;
        enum debugType
        {
            MESH = 1<<IOMESH,
            //SCALARLEVELS = 1<<IOSCALARLEVELS,
            OBJINTERSECTIONS = 1<<IOOBJINTERSECTIONS,
            FEATURESEEDS = 1<<IOFEATURESEEDS,
            ATTRACTION = 1<< IOATTRACTION,
            LAYERINFO = 1<<IOLAYERINFO
        };

        //- Enumeration for what to output
        enum IOoutputType
        {
            IOOUTPUTLAYERINFO
        };
        static const NamedEnum<IOoutputType, 1> IOoutputTypeNames;
        enum outputType
        {
            OUTPUTLAYERINFO = 1<<IOOUTPUTLAYERINFO
        };

        //- Enumeration for what to write
        enum IOwriteType
        {
            IOWRITEMESH,
            IOWRITELEVELS,
            IOWRITELAYERSETS,
            IOWRITELAYERFIELDS
        };
        static const NamedEnum<IOwriteType, 4> IOwriteTypeNames;
        enum writeType
        {
            WRITEMESH = 1<<IOWRITEMESH,
            WRITELEVELS = 1<<IOWRITELEVELS,
            WRITELAYERSETS = 1<<IOWRITELAYERSETS,
            WRITELAYERFIELDS = 1<<IOWRITELAYERFIELDS
        };

        //- Enumeration for how the userdata is to be mapped upon refinement.
        enum mapType
        {
            MASTERONLY = 1, /*!< maintain master only */
            KEEPALL = 2,    /*!< have slaves (upon refinement) from master */
            REMOVE = 4      /*!< set value to -1 any face that was refined */
        };


private:

    // Static data members

        //- Control of writing level
        static writeType writeLevel_;

        //- Control of output/log level
        static outputType outputLevel_;


    // Private data

        //- Reference to mesh
        fvMesh& mesh_;

        //- tolerance used for sorting coordinates (used in 'less' routine)
        const scalar mergeDistance_;

        //- overwrite the mesh?
        const bool overwrite_;

        //- Instance of mesh upon construction. Used when in overwrite_ mode.
        const word oldInstance_;

        //- All surface-intersection interaction
        const refinementSurfaces& surfaces_;

        //- All feature-edge interaction
        const refinementFeatures& features_;

        //- All shell-refinement interaction
        const shellSurfaces& shells_;

        //- refinement engine
        hexRef8 meshCutter_;

        //- per cc-cc vector the index of the surface hit
        labelIOList surfaceIndex_;

        //- user supplied face based data.
        List<Tuple2<mapType, labelList> > userFaceData_;

        //- Meshed patches - are treated differently. Stored as wordList since
        //  order changes.
        wordList meshedPatches_;


    // Private Member Functions

        //- Add patchfield of given type to all fields on mesh
        template<class GeoField>
        static void addPatchFields(fvMesh&, const word& patchFieldType);

        //- Reorder patchfields of all fields on mesh
        template<class GeoField>
        static void reorderPatchFields(fvMesh&, const labelList& oldToNew);

        //- Find out which faces have changed given cells (old mesh labels)
        //  that were marked for refinement.
        static labelList getChangedFaces
        (
            const mapPolyMesh&,
            const labelList& oldCellsToRefine
        );

        //- Calculate coupled boundary end vector and refinement level
        void calcNeighbourData
        (
            labelList& neiLevel,
            pointField& neiCc
        ) const;

        //- Find any intersection of surface. Store in surfaceIndex_.
        void updateIntersections(const labelList& changedFaces);

        //- Remove cells. Put exposedFaces into exposedPatchIDs.
        autoPtr<mapPolyMesh> doRemoveCells
        (
            const labelList& cellsToRemove,
            const labelList& exposedFaces,
            const labelList& exposedPatchIDs,
            removeCells& cellRemover
        );

        // Get cells which are inside any closed surface. Note that
        // all closed surfaces
        // will have already been oriented to have keepPoint outside.
        labelList getInsideCells(const word&) const;

        // Do all to remove inside cells
        autoPtr<mapPolyMesh> removeInsideCells
        (
            const string& msg,
            const label exposedPatchI
        );


        // Refinement candidate selection

            //- Mark cell for refinement (if not already marked). Return false
            // if refinelimit hit. Keeps running count (in nRefine) of cells
            // marked for refinement
            static bool markForRefine
            (
                const label markValue,
                const label nAllowRefine,
                label& cellValue,
                label& nRefine
            );

            //- Mark every cell with level of feature passing through it
            //  (or -1 if not passed through). Uses tracking.
            void markFeatureCellLevel
            (
                const pointField& keepPoints,
                labelList& maxFeatureLevel
            ) const;

            //- Calculate list of cells to refine based on intersection of
            //  features.
            label markFeatureRefinement
            (
                const pointField& keepPoints,
                const label nAllowRefine,

                labelList& refineCell,
                label& nRefine
            ) const;

            //- Mark cells for distance-to-feature based refinement.
            label markInternalDistanceToFeatureRefinement
            (
                const label nAllowRefine,
                labelList& refineCell,
                label& nRefine
            ) const;

            //- Mark cells for refinement-shells based refinement.
            label markInternalRefinement
            (
                const label nAllowRefine,
                labelList& refineCell,
                label& nRefine
            ) const;

            //- Collect faces that are intersected and whose neighbours aren't
            //  yet marked  for refinement.
            labelList getRefineCandidateFaces
            (
                const labelList& refineCell
            ) const;

            //- Mark cells for surface intersection based refinement.
            label markSurfaceRefinement
            (
                const label nAllowRefine,
                const labelList& neiLevel,
                const pointField& neiCc,
                labelList& refineCell,
                label& nRefine
            ) const;

            //- Helper: count number of normals1 that are in normals2
            label countMatches
            (
                const List<point>& normals1,
                const List<point>& normals2,
                const scalar tol = 1e-6
            ) const;

            //- Mark cells for surface curvature based refinement. Marks if
            //  local curvature > curvature or if on different regions
            //  (markDifferingRegions)
            label markSurfaceCurvatureRefinement
            (
                const scalar curvature,
                const label nAllowRefine,
                const labelList& neiLevel,
                const pointField& neiCc,
                labelList& refineCell,
                label& nRefine
            ) const;

            //- Mark cell if local a gap topology or
            bool checkProximity
            (
                const scalar planarCos,
                const label nAllowRefine,

                const label surfaceLevel,
                const vector& surfaceLocation,
                const vector& surfaceNormal,

                const label cellI,

                label& cellMaxLevel,
                vector& cellMaxLocation,
                vector& cellMaxNormal,

                labelList& refineCell,
                label& nRefine
            ) const;

            //- Mark cells for surface proximity based refinement.
            label markProximityRefinement
            (
                const scalar curvature,
                const label nAllowRefine,
                const labelList& neiLevel,
                const pointField& neiCc,

                labelList& refineCell,
                label& nRefine
            ) const;

        // Baffle handling

            //- Get faces to repatch. Returns map from face to patch.
            Map<labelPair> getZoneBafflePatches
            (
                const bool allowBoundary,
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch
            ) const;

            //- Determine patches for baffles
            void getBafflePatches
            (
                const labelList& globalToMasterPatch,
                const labelList& neiLevel,
                const pointField& neiCc,
                labelList& ownPatch,
                labelList& neiPatch
            ) const;

            //- Repatches external face or creates baffle for internal face
            //  with user specified patches (might be different for both sides).
            //  Returns label of added face.
            label createBaffle
            (
                const label faceI,
                const label ownPatch,
                const label neiPatch,
                polyTopoChange& meshMod
            ) const;

        // Problem cell handling

            //- Helper function to mark face as being on 'boundary'. Used by
            //  markFacesOnProblemCells
            void markBoundaryFace
            (
                const label faceI,
                boolList& isBoundaryFace,
                boolList& isBoundaryEdge,
                boolList& isBoundaryPoint
            ) const;

            void findNearest
            (
                const labelList& meshFaces,
                List<pointIndexHit>& nearestInfo,
                labelList& nearestSurface,
                labelList& nearestRegion,
                vectorField& nearestNormal
            ) const;

            Map<label> findEdgeConnectedProblemCells
            (
                const scalarField& perpendicularAngle,
                const labelList&
            ) const;

            bool isCollapsedFace
            (
                const pointField&,
                const pointField& neiCc,
                const scalar minFaceArea,
                const scalar maxNonOrtho,
                const label faceI
            ) const;

            bool isCollapsedCell
            (
                const pointField&,
                const scalar volFraction,
                const label cellI
            ) const;

            //- Returns list with for every internal face -1 or the patch
            //  they should be baffled into. If removeEdgeConnectedCells is set
            //  removes cells based on perpendicularAngle.
            labelList markFacesOnProblemCells
            (
                const dictionary& motionDict,
                const bool removeEdgeConnectedCells,
                const scalarField& perpendicularAngle,
                const labelList& globalToMasterPatch
            ) const;

            //- Returns list with for every face the label of the nearest
            //  patch. Any unreached face (disconnected mesh?) becomes
            //  adaptPatchIDs[0]
            labelList nearestPatch(const labelList& adaptPatchIDs) const;

            //- Returns list with for every internal face -1 or the patch
            //  they should be baffled into.
            labelList markFacesOnProblemCellsGeometric
            (
                const snapParameters& snapParams,
                const dictionary& motionDict
            ) const;


        // Baffle merging

            //- Extract those baffles (duplicate) faces that are on the edge
            //  of a baffle region. These are candidates for merging.
            List<labelPair> freeStandingBaffles
            (
                const List<labelPair>&,
                const scalar freeStandingAngle
            ) const;


        // Zone handling

            //- Finds zone per cell for cells inside closed named surfaces.
            //  (uses geometric test for insideness)
            //  Adapts namedSurfaceIndex so all faces on boundary of cellZone
            //  have corresponding faceZone.
            void findCellZoneGeometric
            (
                const pointField& neiCc,
                const labelList& closedNamedSurfaces,
                labelList& namedSurfaceIndex,
                const labelList& surfaceToCellZone,
                labelList& cellToZone
            ) const;

            //- Finds zone per cell for cells inside named surfaces that have
            //  an inside point specified.
            void findCellZoneInsideWalk
            (
                const labelList& locationSurfaces,
                const labelList& namedSurfaceIndex,
                const labelList& surfaceToCellZone,
                labelList& cellToZone
            ) const;

            //- Determines cell zone from cell region information.
            bool calcRegionToZone
            (
                const label surfZoneI,
                const label ownRegion,
                const label neiRegion,

                labelList& regionToCellZone
            ) const;

            //- Finds zone per cell. Uses topological walk with all faces
            //  marked in namedSurfaceIndex regarded as blocked.
            void findCellZoneTopo
            (
                const point& keepPoint,
                const labelList& namedSurfaceIndex,
                const labelList& surfaceToCellZone,
                labelList& cellToZone
            ) const;

            void makeConsistentFaceIndex
            (
                const labelList& cellToZone,
                labelList& namedSurfaceIndex
            ) const;

            //- Remove any loose standing cells
            void handleSnapProblems
            (
                const snapParameters& snapParams,
                const bool useTopologicalSnapDetection,
                const bool removeEdgeConnectedCells,
                const scalarField& perpendicularAngle,
                const dictionary& motionDict,
                Time& runTime,
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch
            );


            // Some patch utilities

            //- Get all faces in namedSurfaceIndex that have no cellZone on
            //  either side.
            labelList freeStandingBaffleFaces
            (
                const labelList& surfaceToCellZone,
                const labelList& namedSurfaceIndex,
                const labelList& cellToZone,
                const labelList& neiCellZone
            ) const;

            //- Determine per patch edge the number of master faces. Used
            //  to detect non-manifold situations.
            void calcPatchNumMasterFaces
            (
                const PackedBoolList& isMasterFace,
                const indirectPrimitivePatch& patch,
                labelList& nMasterFaces
            ) const;

            //- Determine per patch face the (singly-) connected zone it
            //  is in. Return overall number of zones.
            label markPatchZones
            (
                const indirectPrimitivePatch& patch,
                const labelList& nMasterFaces,
                labelList& faceToZone
            ) const;

            //- Make faces consistent.
            void consistentOrientation
            (
                const PackedBoolList& isMasterFace,
                const indirectPrimitivePatch& patch,
                const labelList& nMasterFaces,
                const labelList& faceToZone,
                const Map<label>& zoneToOrientation,
                PackedBoolList& meshFlipMap
            ) const;


        //- Disallow default bitwise copy construct
        meshRefinement(const meshRefinement&);

        //- Disallow default bitwise assignment
        void operator=(const meshRefinement&);

public:

    //- Runtime type information
    ClassName("meshRefinement");


    // Constructors

        //- Construct from components
        meshRefinement
        (
            fvMesh& mesh,
            const scalar mergeDistance,
            const bool overwrite,
            const refinementSurfaces&,
            const refinementFeatures&,
            const shellSurfaces&
        );


    // Member Functions

        // Access

            //- reference to mesh
            const fvMesh& mesh() const
            {
                return mesh_;
            }
            fvMesh& mesh()
            {
                return mesh_;
            }

            scalar mergeDistance() const
            {
                return mergeDistance_;
            }

            //- Overwrite the mesh?
            bool overwrite() const
            {
                return overwrite_;
            }

            //- (points)instance of mesh upon construction
            const word& oldInstance() const
            {
                return oldInstance_;
            }

            //- reference to surface search engines
            const refinementSurfaces& surfaces() const
            {
                return surfaces_;
            }

            //- reference to feature edge mesh
            const refinementFeatures& features() const
            {
                return features_;
            }

            //- reference to refinement shells (regions)
            const shellSurfaces& shells() const
            {
                return shells_;
            }

            //- reference to meshcutting engine
            const hexRef8& meshCutter() const
            {
                return meshCutter_;
            }

            //- per start-end edge the index of the surface hit
            const labelList& surfaceIndex() const
            {
                return surfaceIndex_;
            }

            labelList& surfaceIndex()
            {
                return surfaceIndex_;
            }

            //- Additional face data that is maintained across
            //  topo changes. Every entry is a list over all faces.
            //  Bit of a hack. Additional flag to say whether to maintain master
            //  only (false) or increase set to account for face-from-face.
            const List<Tuple2<mapType, labelList> >& userFaceData() const
            {
                return userFaceData_;
            }

            List<Tuple2<mapType, labelList> >& userFaceData()
            {
                return userFaceData_;
            }


        // Other

            //- Count number of intersections (local)
            label countHits() const;

            //- Redecompose according to cell count
            //  keepZoneFaces : find all faceZones from zoned surfaces and keep
            //                  owner and neighbour together
            //  keepBaffles   : find all baffles and keep them together
            autoPtr<mapDistributePolyMesh> balance
            (
                const bool keepZoneFaces,
                const bool keepBaffles,
                const scalarField& cellWeights,
                decompositionMethod& decomposer,
                fvMeshDistribute& distributor
            );

            //- Get faces with intersection.
            labelList intersectedFaces() const;

            //- Get points on surfaces with intersection and boundary faces.
            labelList intersectedPoints() const;

            //- Create patch from set of patches
            static autoPtr<indirectPrimitivePatch> makePatch
            (
                const polyMesh&,
                const labelList&
            );

            //- Helper function to make a pointVectorField with correct
            //  bcs for mesh movement:
            //  - adaptPatchIDs         : fixedValue
            //  - processor             : calculated (so free to move)
            //  - cyclic/wedge/symmetry : slip
            //  - other                 : slip
            static tmp<pointVectorField> makeDisplacementField
            (
                const pointMesh& pMesh,
                const labelList& adaptPatchIDs
            );

            //- Helper function: check that face zones are synced
            static void checkCoupledFaceZones(const polyMesh&);

            //- Helper: calculate edge weights (1/length)
            static void calculateEdgeWeights
            (
                const polyMesh& mesh,
                const PackedBoolList& isMasterEdge,
                const labelList& meshPoints,
                const edgeList& edges,
                scalarField& edgeWeights,
                scalarField& invSumWeight
            );

            //- Helper: weighted sum (over all subset of mesh points) by
            //  summing contribution from (master) edges
            template<class Type>
            static void weightedSum
            (
                const polyMesh& mesh,
                const PackedBoolList& isMasterEdge,
                const labelList& meshPoints,
                const edgeList& edges,
                const scalarField& edgeWeights,
                const Field<Type>& data,
                Field<Type>& sum
            );


        // Refinement

            //- Is local topology a small gap?
            bool isGap
            (
                const scalar,
                const vector&,
                const vector&,
                const vector&,
                const vector&
            ) const;

            //- Is local topology a small gap normal to the test vector
            bool isNormalGap
            (
                const scalar,
                const vector&,
                const vector&,
                const vector&,
                const vector&
            ) const;

            //- Calculate list of cells to refine.
            labelList refineCandidates
            (
                const pointField& keepPoints,
                const scalar curvature,
                const scalar planarAngle,

                const bool featureRefinement,
                const bool featureDistanceRefinement,
                const bool internalRefinement,
                const bool surfaceRefinement,
                const bool curvatureRefinement,
                const bool gapRefinement,
                const label maxGlobalCells,
                const label maxLocalCells
            ) const;

            //- Refine some cells
            autoPtr<mapPolyMesh> refine(const labelList& cellsToRefine);

            //- Refine some cells and rebalance
            autoPtr<mapDistributePolyMesh> refineAndBalance
            (
                const string& msg,
                decompositionMethod& decomposer,
                fvMeshDistribute& distributor,
                const labelList& cellsToRefine,
                const scalar maxLoadUnbalance
            );

            //- Balance before refining some cells
            autoPtr<mapDistributePolyMesh> balanceAndRefine
            (
                const string& msg,
                decompositionMethod& decomposer,
                fvMeshDistribute& distributor,
                const labelList& cellsToRefine,
                const scalar maxLoadUnbalance
            );


        // Baffle handling

            //- Split off unreachable areas of mesh.
            void baffleAndSplitMesh
            (
                const bool handleSnapProblems,

                // How to remove problem snaps
                const snapParameters& snapParams,
                const bool useTopologicalSnapDetection,
                const bool removeEdgeConnectedCells,
                const scalarField& perpendicularAngle,

                // How to handle free-standing baffles
                const bool mergeFreeStanding,
                const scalar freeStandingAngle,

                const dictionary& motionDict,
                Time& runTime,
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch,
                const point& keepPoint
            );

            //- Split off (with optional buffer layers) unreachable areas
            //  of mesh. Does not introduce baffles.
            autoPtr<mapPolyMesh> splitMesh
            (
                const label nBufferLayers,
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch,
                const point& keepPoint
            );

            //- Find boundary points that connect to more than one cell
            //  region and split them.
            autoPtr<mapPolyMesh> dupNonManifoldPoints(const localPointRegion&);

            //- Find boundary points that connect to more than one cell
            //  region and split them.
            autoPtr<mapPolyMesh> dupNonManifoldPoints();

            //- Create baffle for every internal face where ownPatch != -1.
            //  External faces get repatched according to ownPatch (neiPatch
            //  should be -1 for these)
            autoPtr<mapPolyMesh> createBaffles
            (
                const labelList& ownPatch,
                const labelList& neiPatch
            );

            //- Debug helper: check faceZones are not on processor patches
            void checkZoneFaces() const;

            //- Create baffles for faces straddling zoned surfaces. Return
            //  baffles.
            autoPtr<mapPolyMesh> createZoneBaffles
            (
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch,
                List<labelPair>&
            );

            //- Merge baffles. Gets pairs of faces.
            autoPtr<mapPolyMesh> mergeBaffles(const List<labelPair>&);

            //- Put faces/cells into zones according to surface specification.
            //  Returns null if no zone surfaces present. Region containing
            //  the keepPoint will not be put into a cellZone.
            autoPtr<mapPolyMesh> zonify
            (
                const point& keepPoint,
                const bool allowFreeStandingZoneFaces
            );


        // Other topo changes

            //- Helper:append patch to end of mesh.
            static label appendPatch
            (
                fvMesh&,
                const label insertPatchI,
                const word&,
                const dictionary&
            );

            //- Helper:add patch to mesh. Update all registered fields.
            //  Used by addMeshedPatch to add patches originating from surfaces.
            static label addPatch(fvMesh&, const word& name, const dictionary&);

            //- Add patch originating from meshing. Update meshedPatches_.
            label addMeshedPatch(const word& name, const dictionary&);

            //- Get patchIDs for patches added in addMeshedPatch.
            labelList meshedPatches() const;

            //- Select coupled faces that are not collocated
            void selectSeparatedCoupledFaces(boolList&) const;

            //- Find region point is in. Uses optional perturbation to re-test.
            static label findRegion
            (
                const polyMesh&,
                const labelList& cellRegion,
                const vector& perturbVec,
                const point& p
            );

            //- Split mesh. Keep part containing point.
            autoPtr<mapPolyMesh> splitMeshRegions
            (
                const labelList& globalToMasterPatch,
                const labelList& globalToSlavePatch,
                const point& keepPoint
            );

            //- Split faces into two
            autoPtr<mapPolyMesh> splitFaces
            (
                const labelList& splitFaces,
                const labelPairList& splits
            );

            //- Update local numbering for mesh redistribution
            void distribute(const mapDistributePolyMesh&);

            //- Update for external change to mesh. changedFaces are in new mesh
            //  face labels.
            void updateMesh
            (
                const mapPolyMesh&,
                const labelList& changedFaces
            );

            //- Helper: reorder list according to map.
            template<class T>
            static void updateList
            (
                const labelList& newToOld,
                const T& nullValue,
                List<T>& elems
            );


            // Restoring : is where other processes delete and reinsert data.

                //- Signal points/face/cells for which to store data
                void storeData
                (
                    const labelList& pointsToStore,
                    const labelList& facesToStore,
                    const labelList& cellsToStore
                );

                //- Update local numbering + undo
                //  Data to restore given as new pointlabel + stored pointlabel
                //  (i.e. what was in pointsToStore)
                void updateMesh
                (
                    const mapPolyMesh&,
                    const labelList& changedFaces,
                    const Map<label>& pointsToRestore,
                    const Map<label>& facesToRestore,
                    const Map<label>& cellsToRestore
                );

            // Merging coplanar faces and edges

                //- Merge coplanar faces. preserveFaces is != -1 for faces
                //  to be preserved
                label mergePatchFacesUndo
                (
                    const scalar minCos,
                    const scalar concaveCos,
                    const labelList& patchIDs,
                    const dictionary& motionDict,
                    const labelList& preserveFaces
                );

                autoPtr<mapPolyMesh> doRemovePoints
                (
                    removePoints& pointRemover,
                    const boolList& pointCanBeDeleted
                );

                autoPtr<mapPolyMesh> doRestorePoints
                (
                    removePoints& pointRemover,
                    const labelList& facesToRestore
                );

                labelList collectFaces
                (
                    const labelList& candidateFaces,
                    const labelHashSet& set
                ) const;

                // Pick up faces of cells of faces in set.
                labelList growFaceCellFace
                (
                    const labelHashSet& set
                ) const;

                //- Merge edges, maintain mesh quality. Return global number
                //  of edges merged
                label mergeEdgesUndo
                (
                    const scalar minCos,
                    const dictionary& motionDict
                );


        // Debug/IO

            //- Debugging: check that all faces still obey start()>end()
            void checkData();

            static void testSyncPointList
            (
                const string& msg,
                const polyMesh& mesh,
                const List<scalar>& fld
            );

            static void testSyncPointList
            (
                const string& msg,
                const polyMesh& mesh,
                const List<point>& fld
            );

            //- Compare two lists over all boundary faces
            template<class T>
            void testSyncBoundaryFaceList
            (
                const scalar mergeDistance,
                const string&,
                const UList<T>&,
                const UList<T>&
            ) const;

            //- Print list according to (collected and) sorted coordinate
            template<class T>
            static void collectAndPrint
            (
                const UList<point>& points,
                const UList<T>& data
            );

            //- Determine master point for subset of points. If coupled
            //  chooses only one
            static PackedBoolList getMasterPoints
            (
                const polyMesh& mesh,
                const labelList& meshPoints
            );

            //- Determine master edge for subset of edges. If coupled
            //  chooses only one
            static PackedBoolList getMasterEdges
            (
                const polyMesh& mesh,
                const labelList& meshEdges
            );

            //- Print some mesh stats.
            void printMeshInfo(const bool, const string&) const;

            //- Replacement for Time::timeName() : return oldInstance (if
            //  overwrite_)
            word timeName() const;

            //- Set instance of all local IOobjects
            void setInstance(const fileName&);

            //- Write mesh and all data
            bool write() const;

            //- Write refinement level as volScalarFields for postprocessing
            void dumpRefinementLevel() const;

            //- Debug: Write intersection information to OBJ format
            void dumpIntersections(const fileName& prefix) const;

            //- Do any one of above IO functions
            void write
            (
                const debugType debugFlags,
                const writeType writeFlags,
                const fileName&
            ) const;

            //- Helper: calculate average
            template<class T>
            static T gAverage
            (
                const PackedBoolList& isMasterElem,
                const UList<T>& values
            );

            //- Get/set write level
            static writeType writeLevel();
            static void writeLevel(const writeType);

            //- Get/set output level
            static outputType outputLevel();
            static void outputLevel(const outputType);


            //- Helper: convert wordList into bit pattern using provided
            //  NamedEnum
            template<class Enum>
            static int readFlags(const Enum& namedEnum, const wordList&);

};


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

} // End namespace Foam

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

#ifdef NoRepository
#   include "meshRefinementTemplates.C"
#endif

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

#endif

// ************************************************************************* //
meshRefinement.H (37,347 bytes)   

user4

2015-04-10 15:16

  ~0004598

The logic was indeed assuming that faceZones were never inside a cellZone. I've uploaded a patched src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H and meshRefinementBaffles.C. Could you see if this fixes your problem and doesn't introduce any other ones?

GRAUPS

2015-04-13 16:08

reporter   ~0004609

mattijs, thanks for the patch, I will test in the next couple days and get back to you.

GRAUPS

2015-04-16 22:01

reporter   ~0004613

mattijs, after testing your patch, I confirmed that it fixed the problem I was seeing (face normals on a faceZone inside a cellZone are now consistent). The monitor on the int faceZone now reports the correct values. I tested both in parallel and in serial and didn't observe any notable changes in the results of the simulation or the mesh that was generated. From my end, and for my uses, the patch appears good.

GRAUPS

2015-07-08 15:12

reporter   ~0005053

mattijs, has this/will this patch be incorporated into the OpenFOAM source?

GRAUPS

2015-08-04 21:51

reporter   ~0005197

Henry, Mattijs, can you please push this fix to OpenFOAM-dev? Or is there a reason that it has been left out thus far?

wyldckat

2015-10-25 20:29

updater  

wyldckat

2015-10-25 20:33

updater  

meshRefinement.H_v2 (37,347 bytes)

wyldckat

2015-10-25 20:34

updater  

wyldckat

2015-10-25 20:42

updater   ~0005492

Sorry, the attached files "meshRefinementBaffles.C_v2" and "meshRefinement.H_v2" are Mattijs original attachments adapted to the latest OpenFOAM-dev, along with, the attached case "porous_pipe_non_aligned_serial_v2.tar.gz" that is re-adapted from "porous_pipe_non_aligned_serial.tar.gz" to the latest OpenFOAM-dev. It seems to work as intended with the v2 modifications.

Problem is that I failed to notice sooner that in the OpenFOAM-history repository is another implemented solution that seems to be better that this one. I'll check it out in a few minutes.

wyldckat

2015-10-25 23:03

updater  

wyldckat

2015-10-25 23:13

updater   ~0005493

OK, using the attached case "porous_pipe_non_aligned_serial_v2.tar.gz" and the attached file "meshRefinementBaffles.C_v3" for replacing the file "src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C", it seems to work well at least for this case.

These changes are based on commit bc48d5266868407913914c70096d09d6b5585275 from the OpenFOAM-history repository: https://github.com/OpenCFD/OpenFOAM-history/commit/bc48d5266868407913914c70096d09d6b5585275 - Mattijs committed on 15 April with the message "BUG: #1479: faceZone inside cellZone".


I did my best in adapting the changes, because the code in OpenFOAM-history, namely for this library "autoMesh" has several new features and some that are probably still being worked on, which make the code considerably different from the current one in OpenFOAM-dev. The conflicting changes were made based on what Mattijs had attached here on this bug report.

Given that this v3 adaptation isn't 100% safe, this still needs considerable testing. @GRAUPS: Can you test at least with your cases that need this feature, before this is finally committed to OpenFOAM-dev?

GRAUPS

2015-10-26 16:20

reporter   ~0005494

Bruno, thanks, I really appreciate your work on this. I am unfortunately unable to test this immediately, but will in the next couple days here as my time frees up.

GRAUPS

2015-11-02 16:38

reporter   ~0005542

Bruno, sorry for the delay. I should have time today to test this and get back to you. Expect an update within the next day.

GRAUPS

2015-11-03 15:55

reporter   ~0005548

Bruno, I finally had time to test the meshRefinementBaffles.C_v3 patch that you provided. I tested it on a few different models in both serial and parallel and found no problems with it.

If you find it stable enough, can you also push the change to the newly released OpenFOAM 3.0.x as well as dev?

wyldckat

2015-11-03 16:49

updater   ~0005549

GRAUPS, many thanks for taking the time to test this with several cases!
Henry is who has the final word on integrating the v3 proposed patch.

----------

As for a clarification about what I meant in the comment #5493, namely: «The conflicting changes were made based on what Mattijs had attached here on this bug report.»

There were only two ambiguous locations where the code in OpenCFD's History repository is considerably different from the code in the Foundation's Development repository:

First location is shown in the attached image "001.png", which shows:

  1. on the left the changes provided on the bug tracker (the namely the v2 proposition that was based on the first attached code by Mattijs);
  2. on the centre is the current version in the Development repository;
  3. and on the right is the change made based on the code located on the History repository (the v3 proposition). The respective piece of code is shown here: https://github.com/OpenCFD/OpenFOAM-history/commit/bc48d5266868407913914c70096d09d6b5585275#diff-346983d8c32ae91984f813eb4d8c0412L2853

Please note that:

 * The v3 modification is completely based on the code that is part of the history repository.
 * The only detail used from the original file provided in the bug tracker by Mattijs is regarding the actual location of where the fix should be located.
 * There is an additional code modification that is not shown in the commit itself, but the "else if (ownZone == maxZone)" block was removed for v3, based on the existing code on the history repository and that this block can be considered redundant.


The second location is similar to the first location, as shown in the attached image "002.png":

 * The location shown at Github: https://github.com/OpenCFD/OpenFOAM-history/commit/bc48d5266868407913914c70096d09d6b5585275#diff-346983d8c32ae91984f813eb4d8c0412L2911
 * The code in v3 is based on the code present in the History repository.
 * The original file that Mattijs provided was only used for figuring out the actual location where the fix should be implemented.

The code originally placed on the bug report by Mattijs was not used directly in the proposed v3 fix; it was only used as a rough guideline of "where to look for the issues".

wyldckat

2015-11-03 16:50

updater  

001.png (173,840 bytes)   
001.png (173,840 bytes)   

wyldckat

2015-11-03 16:50

updater  

002.png (179,984 bytes)   
002.png (179,984 bytes)   

henry

2015-11-03 16:54

manager   ~0005550

I am setting-up OpenFOAM-3.0.x at the moment and when ready will apply these changes to it as well as OpenFOAM-dev.

henry

2015-11-04 11:58

manager   ~0005554

Thanks Bruno.

Resolved in OpenFOAM-3.0.x: commit 2fe3551252d375f33b29046d6733942e88f3a855
Resolved in OpenFOAM-dev: commit 6206280471fa1a1177711a8f6c5869e04c02baf7

henry

2015-11-08 22:39

manager   ~0005587

I have had to revert this changes as it causes problems with the creation of AMI patches.

If you run the propeller tutorials you will see

AMI: Creating addressing and weights between 20196 source faces and 20412 target faces
AMI: Patch source sum(weights) min/max/average = 0.504506, 1.00581, 0.975813
AMI: Patch target sum(weights) min/max/average = 0, 1.0012, 0.966063

and the weight of 0 causes failure. Without the patch

AMI: Creating addressing and weights between 18496 source faces and 18720 target faces
AMI: Patch source sum(weights) min/max/average = 0.984261, 1.01509, 1.00007
AMI: Patch target sum(weights) min/max/average = 0.969379, 1.00121, 0.999945

wyldckat

2015-11-09 08:56

updater   ~0005588

I was afraid and expecting something like this would happen :(
I won't be able to look into this before Friday, in case anyone wants to tackle this in the meantime.

wyldckat

2015-11-16 22:55

updater   ~0005623

My apologies for the time spent with the v3 patch that didn't work as intended.

I'm attaching the v4 proposition, which provides the following two files and respective modifications:

  - "meshRefinement.H_v4" - the change is minimal, namely only the method "freeStandingBaffleFaces()" had its argument "namedSurfaceIndex" renamed to "faceToZone".

  - "meshRefinementBaffles.C_v4" - the changes include what had already been proposed in v3, but it now includes the changes that were missing and that are from the OpenFOAM-history repository. The changes are:

    - In the method "freeStandingBaffleFaces()", it now uses the new "faceToZone" identification, instead of the "namedSurfaceIndex", making the code a bit simpler and faster in this method. This simplifies the identification of what "faceZones" are free-standing. This was copy-pasted from the OpenFOAM-history repository as-is, since it has all of the necessary pieces and it does not have any additional features.

    - In the method "zonify()", several changes were made:

      - Added the list conversion from "namedSurfaceIndex" to "faceToZone". This was one of the important details that was missing in v3.

      - The "Get coupled neighbour cellZone" block was already in v3 and is copy-pasted from the OpenFOAM-history repository.

      - The call to "freeStandingBaffleFaces" was updated accordingly.

      - Debug output added for "freeStanding.obj", which was already in v3.

      - The "Mark zones" block was slightly updated for avoiding name collision. The local variable "faceToZone" was renamed to "faceToConnectedZone".

      - The loops with the comments "Put the faces into the correct zone" and "Set owner as no-flip" have been updated accordingly to the code that in OpenFOAM-history is now in its own method, also named "zonify()". This was where the v3 patch was using v2 as a reference. This issue is now cleared up.
        In order to make this a bit more clearer, have a look at the following locations in the OpenFOAM-history repository, for the latest commit tree:

           - https://github.com/OpenCFD/OpenFOAM-history/blob/412cb0fec86a8ba12da5f9338301cf3c1dfbe0b3/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C#L4561
           - https://github.com/OpenCFD/OpenFOAM-history/blob/412cb0fec86a8ba12da5f9338301cf3c1dfbe0b3/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C#L2849

        The detail here is that the content of the method shown in the second link, is unwrapped in the main "zonify" method in the current OpenFOAM-dev: https://github.com/OpenFOAM/OpenFOAM-dev/blob/baa02e6916b9c13fb377d9299566d90e92647dc1/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C#L3293


Therefore, this v4 proposition is fully based on the OpenFOAM-history repository. Mattijs original attachment here in the bug tracker was no longer used even as a reference, since I've managed to properly understand the implementation that is currently in OpenFOAM-history.


I've tested this v4 proposition with:

  - the previously adapted case "porous_pipe_non_aligned_serial_v2.tar.gz" and the faceZones were oriented accordingly;

  - the tutorial case "incompressible/pimpleDyMFoam/propeller", for which the AMI is now properly mapped:

    - the v3 proposition gave this broken mapping summary:

      AMI: Creating addressing and weights between 20196 source faces and 20412 target faces
      AMI: Patch source sum(weights) min/max/average = 0.504506, 1.00581, 0.975813
      AMI: Patch target sum(weights) min/max/average = 0, 1.0012, 0.966063

    - all other implementations (dev/v2/v4/history) give this mapping summary:

      AMI: Creating addressing and weights between 18496 source faces and 18720 target faces
      AMI: Patch source sum(weights) min/max/average = 0.984261, 1.01509, 1.00007
      AMI: Patch target sum(weights) min/max/average = 0.969379, 1.00121, 0.999945


I still want to do more tests with the v4 proposition, namely to run the tutorial cases that use snappyHexMesh in OpenFOAM 3.0.0/x and then compare with the ones generated with the v4 proposition implemented into OpenFOAM-dev, in order to assess if the meshes are identically or nearly identically generated... and to diagnose those that aren't identical. I'm already attaching the v4 changes as a pre-emptive measure, to avoid this getting lost due to some weird reason (e.g. HD failure).

wyldckat

2015-11-16 22:55

updater  

wyldckat

2015-11-16 22:56

updater  

meshRefinement.H_v4 (37,275 bytes)

wyldckat

2015-11-21 17:04

updater  

Report summarizing any differences in meshing between 3.0.0 and dev+v4.txt (5,979 bytes)   
3.0.0 (32-bit labels) vs dev (64-bit labels):

compressible/rhoPimpleDyMFoam/annularThermalMixer
  - Same mesh features, but two files stand out as very strange:
    - constant/polyMesh/owner
    - constant/polyMesh/faces
  
  - Oddly enough, the "constant/polyMesh/faceZones" is identical, including the
    flipMaps.
  
  - But according the "log.snappyHexMesh", the meshing had several issues on
    both versions... so it's natural that there are some issues.

  - The critical difference is revealed by a full checkMesh:

        --- 3.0.0
        +++ dev
        -             rotorBlades      540      682  ok (non-closed singly connected) (-0.0375 -0.0375 0.04) (0.0375 0.0375 0.08)
        -       rotorBlades_slave      540      682  ok (non-closed singly connected) (-0.0375 -0.0375 0.04) (0.0375 0.0375 0.08)
        +             rotorBlades      540      640  ok (non-closed singly connected) (-0.0375 -0.0375 0.04) (0.0375 0.0375 0.08)
        +       rotorBlades_slave      540      640  ok (non-closed singly connected) (-0.0375 -0.0375 0.04) (0.0375 0.0375 0.08)

  - Using the filter "Clean To Grid" in ParaView, which removes duplicate
    points, revealed that 640 is the correct number of unique points.
    Therefore, the excess points is actually a bug in 3.0.0!

  - Furthermore, when using the filter "Normal Glyphs" on the surface mesh for
    the patch "rotorBlades", the normals are all well represented when the
    mesh was generated with the v4 patch on OpenFOAM-dev, which didn't happen
    with 3.0.0.


heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges
  - Identical meshes, including points on all digits (ascii, precision 6).


heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater
  - The master mesh has a virtually identical mesh, although:
    - faces, owner, neighour are identical;
    - points are similar (numerical precision related issues)
    - cellZones are identical;
    - the faceZones have different flipmaps... which is what we changed with v4.

  - "bottomAir" has a nearly identical mesh, only the points differ very little
    between them.
  
  - Other regions have a different arrangement of patches, making it a bit
    difficult to assess directly from comparing the files.

  - "constant/cellToRegion" files are identical.
  
  - Full checkMesh for the main mesh and for each region only reveals that the
    point map is different enough to give some small differences in checks.
    For example:
  
        --- 3.0.0
        +++ dev
            All angles in faces OK.
        -    Face flatness (1 = flat, 0 = butterfly) : min = 0.9799881  average = 0.9999923
        +    Face flatness (1 = flat, 0 = butterfly) : min = 0.9797262  average = 0.9999923
            All face flatness OK.
        -    Cell determinant (wellposedness) : minimum: 0.9215865 average: 10.2842
        +    Cell determinant (wellposedness) : minimum: 0.9215865 average: 10.28414
            Cell determinant check OK.
        - ***Concave cells (using face planes) found, number of cells: 944
        -  <<Writing 944 concave cells to set concaveCells
        -    Face interpolation weight : minimum: 0.3264674 average: 0.4792975
        + ***Concave cells (using face planes) found, number of cells: 933
        +  <<Writing 933 concave cells to set concaveCells
        +    Face interpolation weight : minimum: 0.3263706 average: 0.4792962
            Face interpolation weight check OK.
        -    Face volume ratio : minimum: 0.1168785 average: 0.8916654
        +    Face volume ratio : minimum: 0.1167739 average: 0.8916526
            Face volume ratio check OK.

  - Residuals and Courant Numbers all seemed within the same scale of
    numerical differences.


incompressible/pimpleDyMFoam/propeller
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


incompressible/simpleFoam/motorBike
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).

  
incompressible/simpleFoam/rotorDisk
  - Identical meshes, including points on all digits (ascii, precision 6).


incompressible/simpleFoam/turbineSiting
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


incompressible/simpleFoam/windAroundBuildings
  - Identical meshes, including points on all digits (ascii, precision 6).


lagrangian/MPPICFoam/cyclone/
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


mesh/snappyHexMesh/flange
  - Identical meshes, including points on all digits (ascii, precision 6).


multiphase/interDyMFoam/ras/DTCHull
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


multiphase/interDyMFoam/ras/mixerVesselAMI
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).

multiphase/interFoam/ras/DTCHull
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


multiphase/interPhaseChangeDyMFoam/propeller
  - Very hard to compare based on files, since it ran in binary mode by default.
  - Using "foamFormatConvert -constant" revealed the meshes to be identical
    (ascii, precision 12).


multiphase/interPhaseChangeFoam/cavitatingBullet
  - Identical meshes, including points on all digits (ascii, precision 6).

wyldckat

2015-11-21 17:19

updater   ~0005649

@Henry: The check is now completed! snappyHexMesh is operating properly, at least with all of the tutorials! Attached is the file "Report summarizing any differences in meshing between 3.0.0 and dev+v4.txt" that has the full report.

The only tutorial left out from this battery testing was "incompressible/pisoFoam/les/motorBike/", because the case generates a mesh that requires more than 6GB of system RAM, the current limit of my machine at home.

Most of the cases gave either identical meshes or very similar meshes. The only two with substantial differences were:

 - "compressible/rhoPimpleDyMFoam/annularThermalMixer" - The mesh generated with OpenFOAM-dev with the v4 proposition gave a better mesh!
    - Not only the normals on the baffles/patches "rotorBlades.*" were all now uniformly oriented;
    - but it also revealed a bug in OpenFOAM 3.0.0 where it generated duplicate points on these baffles/patches, namely 682 points on 3.0.0 versus 640 on OpenFOAM-dev with v4 proposition.

 - "heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater" - The master mesh (the one used for creating the regions) has a virtually identical mesh, although:
    - faces, owner, neighour are identical;
    - points are similar (numerical precision related issues)
    - cellZones are identical;
    - the faceZones have different flipmaps... which is what we changed with v4.


The other detail is that I compared OpenFOAM 3.0.0 with 32-bit labels versus OpenFOAM-dev with 64-bit labels (commit baa02e6916b9), and it still gave identical meshes for most of the cases! It is the desired result, but at least this confirms that 64-bit labels are not affecting the mesh generation, or at least not with most of the tutorials.


All of this to say: proposition v4 is now tested and all signs point toward it now working properly and as intended!

Reminder: there are 2 v4 files, both for folder "src/mesh/autoMesh/autoHexMesh/meshRefinement/":

  - meshRefinementBaffles.C_v4
  - meshRefinement.H_v4

GRAUPS

2015-11-21 18:09

reporter   ~0005650

@Bruno @Henry: Thanks so much for your support and hard work! I look forward to the patch being implemented in 3.0.0 and dev!

henry

2015-11-21 21:37

manager   ~0005653

Thanks Bruno.

Resolved in OpenFOAM-dev by commit ccef40cc0a08d117fc9cfc4729401d1d29769059
Resolved in OpenFOAM-3.0.x by commit e48f65e34d9e5f3704129ae0144018ecb62c1db4

Issue History

Date Modified Username Field Change
2015-01-08 16:52 GRAUPS New Issue
2015-01-08 16:52 GRAUPS File Added: non_aligned_pipes.tar.gz
2015-03-02 12:53 user4 Note Added: 0003940
2015-03-04 18:49 GRAUPS Note Added: 0003956
2015-03-04 18:50 GRAUPS File Added: porous_pipe_non_aligned_serial.tar.gz
2015-03-12 17:02 user4 Note Added: 0004096
2015-03-16 16:58 GRAUPS Note Added: 0004147
2015-03-24 00:17 liuhuafei Issue cloned: 0001603
2015-03-30 16:34 user4 File Added: normals.png
2015-03-30 16:39 user4 Note Added: 0004552
2015-03-30 18:12 GRAUPS Note Added: 0004553
2015-04-10 15:14 user4 File Added: meshRefinementBaffles.C
2015-04-10 15:15 user4 File Added: meshRefinement.H
2015-04-10 15:16 user4 Note Added: 0004598
2015-04-13 16:08 GRAUPS Note Added: 0004609
2015-04-16 22:01 GRAUPS Note Added: 0004613
2015-07-08 15:12 GRAUPS Note Added: 0005053
2015-08-04 21:51 GRAUPS Note Added: 0005197
2015-10-25 20:29 wyldckat File Added: porous_pipe_non_aligned_serial_v2.tar.gz
2015-10-25 20:33 wyldckat File Added: meshRefinement.H_v2
2015-10-25 20:34 wyldckat File Added: meshRefinementBaffles.C_v2
2015-10-25 20:42 wyldckat Note Added: 0005492
2015-10-25 23:03 wyldckat File Added: meshRefinementBaffles.C_v3
2015-10-25 23:13 wyldckat Note Added: 0005493
2015-10-26 16:20 GRAUPS Note Added: 0005494
2015-11-02 16:38 GRAUPS Note Added: 0005542
2015-11-03 15:55 GRAUPS Note Added: 0005548
2015-11-03 16:49 wyldckat Note Added: 0005549
2015-11-03 16:50 wyldckat File Added: 001.png
2015-11-03 16:50 wyldckat File Added: 002.png
2015-11-03 16:54 henry Note Added: 0005550
2015-11-04 11:58 henry Note Added: 0005554
2015-11-04 11:58 henry Status new => resolved
2015-11-04 11:58 henry Resolution open => fixed
2015-11-04 11:58 henry Assigned To => henry
2015-11-08 22:39 henry Note Added: 0005587
2015-11-08 22:39 henry Status resolved => feedback
2015-11-08 22:39 henry Resolution fixed => reopened
2015-11-09 08:56 wyldckat Note Added: 0005588
2015-11-16 22:55 wyldckat Note Added: 0005623
2015-11-16 22:55 wyldckat File Added: meshRefinementBaffles.C_v4
2015-11-16 22:56 wyldckat File Added: meshRefinement.H_v4
2015-11-21 17:04 wyldckat File Added: Report summarizing any differences in meshing between 3.0.0 and dev+v4.txt
2015-11-21 17:19 wyldckat Note Added: 0005649
2015-11-21 17:19 wyldckat Status feedback => assigned
2015-11-21 18:09 GRAUPS Note Added: 0005650
2015-11-21 21:37 henry Note Added: 0005653
2015-11-21 21:37 henry Status assigned => resolved
2015-11-21 21:37 henry Resolution reopened => fixed