Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -270,7 +270,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G = 2*mut_*(tgradU() && dev(symm(tgradU())));
+    volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
 
     // Update espsilon and G at the wall

applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C

@@ -1019,6 +1019,7 @@ int main(int argc, char *argv[])
             }
         }
         writer.writeEnd();
+        Info<< endl;
 
 
 
@@ -1067,7 +1068,7 @@ int main(int argc, char *argv[])
             writer.writeInit
             (
                 runTime.caseName(),
-                cellVarNames,
+                allVarNames,
                 patchFileName,
                 DataFileType_Full
             );
@@ -1163,6 +1164,8 @@ int main(int argc, char *argv[])
                         << nl << endl;
                 }
             }
+            writer.writeEnd();
+            Info<< endl;
         }
 
 

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C

@@ -38,7 +38,17 @@ Foam::IPstream::IPstream
 )
 :
     Pstream(commsType, bufSize),
-    UIPstream(commsType, fromProcNo, buf_, externalBufPosition_),
+    UIPstream
+    (
+        commsType,
+        fromProcNo,
+        buf_,
+        externalBufPosition_,
+        UPstream::msgType(),        // tag
+        false,                      // do not clear buf_ if at end
+        format,
+        version
+    ),
     externalBufPosition_(0)
 {}
 

src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C

@@ -77,6 +77,21 @@ inline void Foam::UIPstream::readFromBuffer
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::UIPstream::~UIPstream()
+{
+    if (clearAtEnd_ && eof())
+    {
+        if (debug)
+        {
+            Pout<< "UIPstream::~UIPstream() : clearing externalBuf_ of size "
+                << externalBuf_.size() << endl;
+        }
+        externalBuf_.clearStorage();
+    }
+}
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::UIPstream::read(token& t)

src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.H

@@ -67,6 +67,8 @@ class UIPstream
 
         const int tag_;
 
+        const bool clearAtEnd_;
+
         label messageSize_;
 
 
@@ -96,6 +98,7 @@ public:
             DynamicList<char>& externalBuf,
             label& externalBufPosition,
             const int tag = UPstream::msgType(),
+            const bool clearAtEnd = false,   // destroy externalBuf if at end
             streamFormat format=BINARY,
             versionNumber version=currentVersion
         );
@@ -104,6 +107,11 @@ public:
         UIPstream(const int fromProcNo, PstreamBuffers&);
 
 
+    // Destructor
+
+        ~UIPstream();
+
+
     // Member functions
 
         // Inquiry

src/Pstream/dummy/UIPread.C

@@ -38,6 +38,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -48,6 +49,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     notImplemented
@@ -59,6 +61,7 @@ Foam::UIPstream::UIPstream
             "DynamicList<char>&,"
             "label&,"
             "const int tag,"
+            "const bool,"
             "streamFormat, versionNumber"
         ")"
     );

src/Pstream/mpi/UIPread.C

@@ -42,6 +42,7 @@ Foam::UIPstream::UIPstream
     DynamicList<char>& externalBuf,
     label& externalBufPosition,
     const int tag,
+    const bool clearAtEnd,
     streamFormat format,
     versionNumber version
 )
@@ -52,6 +53,7 @@ Foam::UIPstream::UIPstream
     externalBuf_(externalBuf),
     externalBufPosition_(externalBufPosition),
     tag_(tag),
+    clearAtEnd_(clearAtEnd),
     messageSize_(0)
 {
     setOpened();
@@ -125,6 +127,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
     externalBuf_(buffers.recvBuf_[fromProcNo]),
     externalBufPosition_(buffers.recvBufPos_[fromProcNo]),
     tag_(buffers.tag_),
+    clearAtEnd_(true),
     messageSize_(0)
 {
     if (commsType() != UPstream::scheduled && !buffers.finishedSendsCalled_)

src/fvMotionSolver/Make/files

@@ -1,6 +1,7 @@
 fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/displacementFvMotionSolver/displacementFvMotionSolver.C
+fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C
 fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
 fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.C

@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementLayeredMotionFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointEdgeStructuredWalk.H"
+#include "pointFields.H"
+#include "PointEdgeWave.H"
+#include "syncTools.H"
+#include "interpolationTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementLayeredMotionFvMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        fvMotionSolver,
+        displacementLayeredMotionFvMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::displacementLayeredMotionFvMotionSolver::calcZoneMask
+(
+    const label cellZoneI, 
+    PackedBoolList& isZonePoint,
+    PackedBoolList& isZoneEdge
+) const
+{
+    if (cellZoneI == -1)
+    {
+        isZonePoint.setSize(mesh().nPoints());
+        isZonePoint = 1;
+
+        isZoneEdge.setSize(mesh().nEdges());
+        isZoneEdge = 1;
+    }
+    else
+    {
+        const cellZone& cz = mesh().cellZones()[cellZoneI];
+
+        label nPoints = 0;
+        forAll(cz, i)
+        {
+            const labelList& cPoints = mesh().cellPoints(cz[i]);
+            forAll(cPoints, cPointI)
+            {
+                if (!isZonePoint[cPoints[cPointI]])
+                {
+                    isZonePoint[cPoints[cPointI]] = 1;
+                    nPoints++;
+                }
+            }
+        }
+        syncTools::syncPointList
+        (
+            mesh(),
+            isZonePoint,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        // Mark edge inside cellZone
+        label nEdges = 0;
+        forAll(cz, i)
+        {
+            const labelList& cEdges = mesh().cellEdges(cz[i]);
+            forAll(cEdges, cEdgeI)
+            {
+                if (!isZoneEdge[cEdges[cEdgeI]])
+                {
+                    isZoneEdge[cEdges[cEdgeI]] = 1;
+                    nEdges++;
+                }
+            }
+        }
+        syncTools::syncEdgeList
+        (
+            mesh(),
+            isZoneEdge,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+        Info<< "On cellZone " << cz.name()
+            << " marked " << returnReduce(nPoints, sumOp<label>())
+            << " points and " << returnReduce(nEdges, sumOp<label>())
+            << " edges." << endl;
+    }
+}
+
+
+// Find distance to starting point
+void Foam::displacementLayeredMotionFvMotionSolver::walkStructured
+(
+    const label cellZoneI,
+    const PackedBoolList& isZonePoint,
+    const PackedBoolList& isZoneEdge,
+    const labelList& seedPoints,
+    const vectorField& seedData,
+    scalarField& distance,
+    vectorField& data
+) const
+{
+    const pointField& points = mesh().points();
+
+    List<pointEdgeStructuredWalk> seedInfo(seedPoints.size());
+
+    forAll(seedPoints, i)
+    {
+        seedInfo[i] = pointEdgeStructuredWalk
+        (
+            true,
+            points[seedPoints[i]],
+            0.0,
+            seedData[i]
+        );
+    }
+
+    // Current info on points
+    List<pointEdgeStructuredWalk> allPointInfo(mesh().nPoints());
+    // Mark points inside cellZone
+    forAll(isZonePoint, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            allPointInfo[pointI].inZone() = true;
+        }
+    }
+
+    // Current info on edges
+    List<pointEdgeStructuredWalk> allEdgeInfo(mesh().nEdges());
+    // Mark edges inside cellZone
+    forAll(isZoneEdge, edgeI)
+    {
+        if (isZoneEdge[edgeI])
+        {
+            allEdgeInfo[edgeI].inZone() = true;
+        }
+    }
+
+    // Walk
+    PointEdgeWave<pointEdgeStructuredWalk> wallCalc
+    (
+        mesh(),
+        seedPoints,
+        seedInfo,
+
+        allPointInfo,
+        allEdgeInfo,
+        mesh().globalData().nTotalPoints()  // max iterations
+    );
+
+    // Extract distance and passive data
+    forAll(allPointInfo, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            distance[pointI] = allPointInfo[pointI].dist();
+            data[pointI] = allPointInfo[pointI].data();
+        }
+    }
+}
+
+
+// Evaluate faceZone patch
+Foam::tmp<Foam::vectorField>
+Foam::displacementLayeredMotionFvMotionSolver::faceZoneEvaluate
+(
+    const faceZone& fz,
+    const labelList& meshPoints,
+    const dictionary& dict,
+    const PtrList<pointVectorField>& patchDisp,
+    const label patchI
+) const
+{
+    tmp<vectorField> tfld(new vectorField(meshPoints.size()));
+    vectorField& fld = tfld();
+
+    const word& type = dict.lookup("type");
+
+    if (type == "fixedValue")
+    {
+        fld = vectorField("value", dict, meshPoints.size());
+    }
+    else if (type == "timeVaryingUniformFixedValue")
+    {
+        interpolationTable<vector> timeSeries(dict);
+
+        fld = timeSeries(mesh().time().timeOutputValue());
+    }
+    else if (type == "slip")
+    {
+        if ((patchI % 2) != 1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+                *this
+            )   << "slip can only be used on second faceZonePatch of pair."
+                << "FaceZone:" << fz.name()
+                << exit(FatalIOError);
+        }
+        // Use field set by previous bc
+        fld = vectorField(patchDisp[patchI-1], meshPoints);
+    }
+    else if (type == "follow")
+    {
+        // Only on boundary faces - follow boundary conditions
+        fld = vectorField(pointDisplacement_, meshPoints);
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::faceZoneEvaluate(..)",
+            *this
+        )   << "Unknown faceZonePatch type " << type << " for faceZone "
+            << fz.name() << exit(FatalIOError);
+    }
+    return tfld;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::cellZoneSolve
+(
+    const label cellZoneI,
+    const dictionary& zoneDict
+)
+{
+    PackedBoolList isZonePoint(mesh().nPoints());
+    PackedBoolList isZoneEdge(mesh().nEdges());
+    calcZoneMask(cellZoneI, isZonePoint, isZoneEdge);
+
+    const dictionary& patchesDict = zoneDict.subDict("boundaryField");
+
+    if (patchesDict.size() != 2)
+    {
+        FatalIOErrorIn
+        (
+            "displacementLayeredMotionFvMotionSolver::"
+            "correctBoundaryConditions(..)",
+            *this
+        )   << "Can only handle 2 faceZones (= patches) per cellZone. "
+            << " cellZone:" << cellZoneI
+            << " patches:" << patchesDict.toc()
+            << exit(FatalIOError);
+    }
+
+    PtrList<scalarField> patchDist(patchesDict.size());
+    PtrList<pointVectorField> patchDisp(patchesDict.size());
+
+    // Allocate the fields
+    label patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::"
+                "correctBoundaryConditions(..)",
+                *this
+            )   << "Cannot find faceZone " << faceZoneName
+                << endl << "Valid zones are " << mesh().faceZones().names()
+                << exit(FatalIOError);
+        }
+
+        // Determine the points of the faceZone within the cellZone
+        const faceZone& fz = mesh().faceZones()[zoneI];
+
+        patchDist.set(patchI, new scalarField(mesh().nPoints()));
+        patchDisp.set
+        (
+            patchI,
+            new pointVectorField
+            (
+                IOobject
+                (
+                    mesh().cellZones()[cellZoneI].name() + "_" + fz.name(),
+                    mesh().time().timeName(),
+                    mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                pointDisplacement_  // to inherit the boundary conditions
+            )
+        );
+
+        patchI++;
+    }
+
+
+
+    // 'correctBoundaryConditions'
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Loops over all the faceZones and walks their boundary values
+
+    // Make sure we can pick up bc values from field
+    pointDisplacement_.correctBoundaryConditions();
+
+    patchI = 0;
+    forAllConstIter(dictionary, patchesDict, patchIter)
+    {
+        const word& faceZoneName = patchIter().keyword();
+        const dictionary& faceZoneDict = patchIter().dict();
+
+        // Determine the points of the faceZone within the cellZone
+        label zoneI = mesh().faceZones().findZoneID(faceZoneName);
+        const faceZone& fz = mesh().faceZones()[zoneI];
+        const labelList& fzMeshPoints = fz().meshPoints();
+        DynamicList<label> meshPoints(fzMeshPoints.size());
+        forAll(fzMeshPoints, i)
+        {
+            if (isZonePoint[fzMeshPoints[i]])
+            {
+                meshPoints.append(fzMeshPoints[i]);
+            }
+        }
+
+        // Get initial value for all the faceZone points
+        tmp<vectorField> tseed = faceZoneEvaluate
+        (
+            fz,
+            meshPoints,
+            faceZoneDict,
+            patchDisp,
+            patchI
+        );
+
+
+Info<< "For cellZone:" << cellZoneI
+    << " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
+    << " have patchField:"
+    << " max:" << gMax(tseed())
+    << " min:" << gMin(tseed())
+    << " avg:" << gAverage(tseed())
+    << endl;
+
+        // Set distance and transported value
+        walkStructured
+        (
+            cellZoneI,
+            isZonePoint,
+            isZoneEdge,
+
+            meshPoints,
+            tseed,
+            patchDist[patchI],
+            patchDisp[patchI]
+        );
+
+        // Implement real bc.
+        patchDisp[patchI].correctBoundaryConditions();
+
+
+//Info<< "Writing displacement for faceZone " << fz.name()
+//    << " to " << patchDisp[patchI].name() << endl;
+//patchDisp[patchI].write();
+
+//        // Copy into pointDisplacement for other fields to use
+//        forAll(isZonePoint, pointI)
+//        {
+//            if (isZonePoint[pointI])
+//            {
+//                pointDisplacement_[pointI] = patchDisp[patchI][pointI];
+//            }
+//        }
+//        pointDisplacement_.correctBoundaryConditions();
+
+
+        patchI++;
+    }
+
+
+    // Solve
+    // ~~~~~
+    // solving the interior is just interpolating
+
+//    // Get normalised distance
+//    pointScalarField distance
+//    (
+//        IOobject
+//        (
+//            "distance",
+//            mesh().time().timeName(),
+//            mesh(),
+//            IOobject::NO_READ,
+//            IOobject::NO_WRITE,
+//            false
+//        ),
+//        pointMesh::New(mesh()),
+//        dimensionedScalar("distance", dimLength, 0.0)
+//    );
+//    forAll(distance, pointI)
+//    {
+//        if (isZonePoint[pointI])
+//        {
+//            scalar d1 = patchDist[0][pointI];
+//            scalar d2 = patchDist[1][pointI];
+//            if (d1+d2 > SMALL)
+//            {
+//                scalar s = d1/(d1+d2);
+//                distance[pointI] = s;
+//            }
+//        }
+//    }
+//    Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
+//        << endl;
+//    distance.write();
+
+    // Average
+    forAll(pointDisplacement_, pointI)
+    {
+        if (isZonePoint[pointI])
+        {
+            scalar d1 = patchDist[0][pointI];
+            scalar d2 = patchDist[1][pointI];
+
+            scalar s = d1/(d1+d2+VSMALL);
+
+            pointDisplacement_[pointI] =
+                (1-s)*patchDisp[0][pointI]
+              + s*patchDisp[1][pointI];
+        }
+    }
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+displacementLayeredMotionFvMotionSolver
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    displacementFvMotionSolver(mesh, is),
+    pointDisplacement_
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            fvMesh_.time().timeName(),
+            fvMesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pointMesh::New(fvMesh_)
+    )
+{
+    pointDisplacement_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementLayeredMotionFvMotionSolver::
+~displacementLayeredMotionFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField>
+Foam::displacementLayeredMotionFvMotionSolver::curPoints() const
+{
+    tmp<pointField> tcurPoints
+    (
+        points0() + pointDisplacement_.internalField()
+    );
+
+    twoDCorrectPoints(tcurPoints());
+
+//    const pointField& pts = tcurPoints();
+//    forAll(pts, pointI)
+//    {
+//        Info<< "    from:" << mesh().points()[pointI]
+//            << " to:" << pts[pointI]
+//            << endl;
+//    }
+
+
+    return tcurPoints;
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::solve()
+{
+    const dictionary& ms = mesh().lookupObject<motionSolver>("dynamicMeshDict");
+    const dictionary& solverDict = ms.subDict(typeName + "Coeffs");
+
+    // Apply all regions (=cellZones)
+
+    const dictionary& regionDicts = solverDict.subDict("regions");
+    forAllConstIter(dictionary, regionDicts, regionIter)
+    {
+        const word& cellZoneName = regionIter().keyword();
+        const dictionary& regionDict = regionIter().dict();
+
+        label zoneI = mesh().cellZones().findZoneID(cellZoneName);
+
+        Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
+            << endl;
+
+        if (zoneI == -1)
+        {
+            FatalIOErrorIn
+            (
+                "displacementLayeredMotionFvMotionSolver::solve(..)",
+                *this
+            )   << "Cannot find cellZone " << cellZoneName
+                << endl << "Valid zones are " << mesh().cellZones().names()
+                << exit(FatalIOError);
+        }
+
+        cellZoneSolve(zoneI, regionDict);
+    }
+}
+
+
+void Foam::displacementLayeredMotionFvMotionSolver::updateMesh
+(
+    const mapPolyMesh& mpm
+)
+{
+    displacementFvMotionSolver::updateMesh(mpm);
+}
+
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionFvMotionSolver.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::displacementLayeredMotionFvMotionSolver
+
+Description
+    Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
+    structure of the mesh blocks and boundary conditions on these blocks.
+
+    Note: should not be an fvMotionSolver but just a motionSolver. Only here
+    so we can reuse displacementFvMotionSolver functionality (e.g. surface
+    following boundary conditions)
+
+
+    The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
+    per region (=cellZone) the boundary conditions on two opposing patches
+    (=faceZones). It then interpolates the boundary values using topological
+    walking so requires the cellZone to be a layered mesh.
+
+    The boundary conditions on faceZones are currently:
+
+    follow: the faceZone follows the overall mesh displacement.
+            Use this for faceZones on boundary faces (so it uses the
+            proper boundary conditions on the pointDisplacement).
+
+    fixedValue: fixed value.
+
+    timeVaryingUniformFixedValue: table-driven fixed value.
+
+    slip: the second of a pair of faceZones follows the tangential movement
+          specified by the first faceZone. (= removes the normal component).
+
+SourceFiles
+    displacementLayeredMotionFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementLayeredMotionFvMotionSolver_H
+#define displacementLayeredMotionFvMotionSolver_H
+
+#include "displacementFvMotionSolver.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+
+/*---------------------------------------------------------------------------*\
+             Class displacementLayeredMotionFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementLayeredMotionFvMotionSolver
+:
+    public displacementFvMotionSolver
+{
+    // Private data
+
+        //- Point motion field
+        mutable pointVectorField pointDisplacement_;
+
+
+    // Private Member Functions
+
+        void calcZoneMask
+        (
+            const label cellZoneI, 
+            PackedBoolList& isZonePoint,
+            PackedBoolList& isZoneEdge
+        ) const;
+
+        void walkStructured
+        (
+            const label cellZoneI,
+            const PackedBoolList& isZonePoint,
+            const PackedBoolList& isZoneEdge,
+            const labelList& seedPoints,
+            const vectorField& seedData,
+            scalarField& distance,
+            vectorField& data
+        ) const;
+
+        tmp<vectorField> faceZoneEvaluate
+        (
+            const faceZone& fz,
+            const labelList& meshPoints,
+            const dictionary& dict,
+            const PtrList<pointVectorField>& patchDisp,
+            const label patchI
+        ) const;
+
+        void cellZoneSolve
+        (
+            const label cellZoneI,
+            const dictionary& zoneDict
+        );
+
+
+        //- Disallow default bitwise copy construct
+        displacementLayeredMotionFvMotionSolver
+        (
+            const displacementLayeredMotionFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementLayeredMotionFvMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementLayeredMotion");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and data stream
+        displacementLayeredMotionFvMotionSolver
+        (
+            const polyMesh&,
+            Istream& msDataUnused
+        );
+
+
+    // Destructor
+
+        ~displacementLayeredMotionFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Return reference to the point motion displacement field
+        pointVectorField& pointDisplacement()
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return const reference to the point motion displacement field
+        const pointVectorField& pointDisplacement() const
+        {
+            return pointDisplacement_;
+        }
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointEdgeStructuredWalk.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::pointEdgeStructuredWalk& wDist
+)
+{
+    return os
+        << wDist.inZone_ << wDist.previousPoint_
+        << wDist.dist_ << wDist.data_;
+}
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::pointEdgeStructuredWalk& wDist
+)
+{
+    return is
+        >> wDist.inZone_ >> wDist.previousPoint_
+        >> wDist.dist_ >> wDist.data_;
+}
+
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::pointEdgeStructuredWalk
+
+Description
+    Determines length of string of edges walked to point.
+
+SourceFiles
+    pointEdgeStructuredWalkI.H
+    pointEdgeStructuredWalk.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointEdgeStructuredWalk_H
+#define pointEdgeStructuredWalk_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "pTraits.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointEdgeStructuredWalk Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointEdgeStructuredWalk
+{
+    // Private data
+
+        //- Is this point/edge inside the walk zone
+        bool inZone_;
+
+        //- Previuos point
+        point previousPoint_;
+
+        //- Sum of distance
+        scalar dist_;
+
+        //- Passive data
+        vector data_;
+
+    // Private Member Functions
+
+        //- Evaluate distance to point.
+        inline bool update
+        (
+            const point&,
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+        //- Combine current with w2. Update distSqr, origin if w2 has smaller
+        //  quantities and returns true.
+        inline bool update
+        (
+            const pointEdgeStructuredWalk& w2,
+            const scalar tol
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointEdgeStructuredWalk();
+
+        //- Construct from components
+        inline pointEdgeStructuredWalk
+        (
+            const bool,
+            const point&,
+            const scalar,
+            const vector&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline bool inZone() const;
+
+            inline bool& inZone();
+
+            inline const point& previousPoint() const;
+
+            inline scalar dist() const;
+
+            inline const vector& data() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether still contains original (invalid) value.
+            inline bool valid() const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            inline bool sameGeometry
+            (
+                const pointEdgeStructuredWalk&,
+                const scalar tol
+            ) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Apply rotation matrix to origin
+            inline void transform(const tensor& rotTensor);
+
+            //- Influence of edge on point
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointEdgeStructuredWalk& edgeInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            inline bool updatePoint
+            (
+                const pointEdgeStructuredWalk& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of point on edge.
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointEdgeStructuredWalk& pointInfo,
+                const scalar tol
+            );
+
+
+    // Member Operators
+
+        //Note: Used to determine whether to call update.
+        inline bool operator==(const pointEdgeStructuredWalk&) const;
+
+        inline bool operator!=(const pointEdgeStructuredWalk&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointEdgeStructuredWalk&);
+        friend Istream& operator>>(Istream&, pointEdgeStructuredWalk&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointEdgeStructuredWalkI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H

@@ -0,0 +1,302 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2.
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const point& pt,
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set. Walked from w2 to here (=pt)
+        dist_ = w2.dist_ + mag(pt-w2.previousPoint_);
+        previousPoint_ = pt;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+    
+// Update this with w2 (information on same point)
+inline bool Foam::pointEdgeStructuredWalk::update
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set so use neighbour
+        dist_ = w2.dist_;
+        previousPoint_ = w2.previousPoint_;
+        data_ = w2.data_;
+
+        return true;
+    }
+    else
+    {
+        // already nearer to pt
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk()
+:
+    inZone_(false),
+    previousPoint_(vector::max),
+    dist_(0),
+    data_(vector::zero)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointEdgeStructuredWalk::pointEdgeStructuredWalk
+(
+    const bool inZone,
+    const point& previousPoint, 
+    const scalar dist,
+    const vector& data
+)
+:
+    inZone_(inZone),
+    previousPoint_(previousPoint),
+    dist_(dist),
+    data_(data)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::inZone() const
+{
+    return inZone_;
+}
+
+
+inline bool& Foam::pointEdgeStructuredWalk::inZone()
+{
+    return inZone_;
+}
+
+
+inline const Foam::point& Foam::pointEdgeStructuredWalk::previousPoint() const
+{
+    return previousPoint_;
+}
+
+
+inline Foam::scalar Foam::pointEdgeStructuredWalk::dist() const
+{
+    return dist_;
+}
+
+
+inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
+{
+    return data_;
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::valid() const
+{
+    return previousPoint_ != vector::max;
+}
+
+
+// Checks for cyclic points
+inline bool Foam::pointEdgeStructuredWalk::sameGeometry
+(
+    const pointEdgeStructuredWalk& w2,
+    const scalar tol
+) const
+{
+    scalar diff = Foam::mag(dist() - w2.dist());
+
+    if (diff < SMALL)
+    {
+        return true;
+    }
+    else
+    {
+        if ((dist() > SMALL) && ((diff/dist()) < tol))
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    previousPoint_ -= coord;
+}
+
+
+inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
+{
+    previousPoint_ = Foam::transform(rotTensor, previousPoint_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (dist_)
+// is not affected by leaving/entering domain.
+inline void Foam::pointEdgeStructuredWalk::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    // back to absolute form
+    previousPoint_ += coord;
+}
+
+
+// Update this with information from connected edge
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointEdgeStructuredWalk& edgeInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], edgeInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// Update this with new information on same point
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        return update(mesh.points()[pointI], newPointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update this with new information on same point. No extra information.
+inline bool Foam::pointEdgeStructuredWalk::updatePoint
+(
+    const pointEdgeStructuredWalk& newPointInfo,
+    const scalar tol
+)
+{
+    return update(newPointInfo, tol);
+}    
+
+
+// Update this with information from connected point
+inline bool Foam::pointEdgeStructuredWalk::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointEdgeStructuredWalk& pointInfo,
+    const scalar tol
+)
+{
+    if (inZone_)
+    {
+        const pointField& points = mesh.points();
+
+        const edge& e = mesh.edges()[edgeI];
+
+        const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
+
+        return update(edgeMid, pointInfo, tol);
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeStructuredWalk::operator==
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return previousPoint() == rhs.previousPoint();
+}
+
+
+inline bool Foam::pointEdgeStructuredWalk::operator!=
+(
+    const Foam::pointEdgeStructuredWalk& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -30,8 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
-// #include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -244,6 +242,14 @@ void Foam::distanceSurface::createGeometry()
             distance_,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    distance_,
+        //    regularise_
+        //)
     );
 
     if (debug)

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -85,6 +85,7 @@ class distanceSurface
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -169,6 +170,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "distanceSurface.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -2054,8 +2054,8 @@ Foam::isoSurface::isoSurface
     }
 
 
-    //if (false)
+    if (false)
-    //{
+    {
         List<FixedList<label, 3> > faceEdges;
         labelList edgeFace0, edgeFace1;
         Map<labelList> edgeFacesRest;
@@ -2064,7 +2064,14 @@ Foam::isoSurface::isoSurface
         while (true)
         {
             // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+            calcAddressing
+            (
+                *this,
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
 
             // See if any dangling triangles
             boolList keepTriangles;
@@ -2108,7 +2115,7 @@ Foam::isoSurface::isoSurface
         }
 
         orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-    //}
+    }
 
 
     if (debug)

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -1398,6 +1398,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
 )
 :
     mesh_(mesh),
+    cVals_(cVals),
+    pVals_(pVals),
     iso_(iso),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
@@ -1563,6 +1565,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
+    if (false)
+    {
         List<FixedList<label, 3> > faceEdges;
         labelList edgeFace0, edgeFace1;
         Map<labelList> edgeFacesRest;
@@ -1571,7 +1575,14 @@ Foam::isoSurfaceCell::isoSurfaceCell
         while (true)
         {
             // Calculate addressing
-        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+            calcAddressing
+            (
+                *this,
+                faceEdges,
+                edgeFace0,
+                edgeFace1,
+                edgeFacesRest
+            );
 
             // See if any dangling triangles
             boolList keepTriangles;
@@ -1615,7 +1626,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
         }
 
         orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-
+    }
     //combineCellTriangles();
 }
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -85,6 +85,10 @@ class isoSurfaceCell
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        const scalarField& cVals_;
+
+        const scalarField& pVals_;
+
         //- isoSurfaceCell value
         const scalar iso_;
 
@@ -352,6 +356,14 @@ public:
             const Field<Type>& cCoords,
             const Field<Type>& pCoords
         ) const;
+
+        //- Interpolates cCoords,pCoords.
+        template <class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
 };
 
 

src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C

@@ -381,4 +381,56 @@ Foam::isoSurfaceCell::interpolate
 }
 
 
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurfaceCell::interpolate
+(
+    const Field<Type>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
+
+
+    generateTriPoints
+    (
+        cVals_,
+        pVals_,
+
+        cCoords,
+        pCoords,
+
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
+
+        triPoints,
+        triMeshCells
+    );
+
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triPoints, i)
+    {
+        label mergedPointI = triPointMergeMap_[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] = triPoints[i];
+        }
+    }
+
+    return tvalues;
+}
+
+
 // ************************************************************************* //

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -30,7 +30,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -231,6 +230,14 @@ void Foam::sampledCuttingPlane::createGeometry()
             0.0,
             regularise_
         )
+        //new isoSurfaceCell
+        //(
+        //    fvm,
+        //    cellDistance,
+        //    pointDistance_,
+        //    0.0,
+        //    regularise_
+        //)
     );
 
     if (debug)

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "isoSurface.H"
+//#include "isoSurfaceCell.H"
 #include "plane.H"
 #include "ZoneIDs.H"
 #include "fvMeshSubset.H"
@@ -89,6 +90,7 @@ class sampledCuttingPlane
         scalarField pointDistance_;
 
         //- Constructed iso surface
+        //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
         //- triangles converted to faceList
@@ -173,6 +175,7 @@ public:
         }
 
 
+        //const isoSurfaceCell& surface() const
         const isoSurface& surface() const
         {
             return isoSurfPtr_();

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C

@@ -26,7 +26,6 @@ License
 
 #include "volPointInterpolation.H"
 #include "sampledCuttingPlane.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -25,11 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledTriSurfaceMesh.H"
-#include "dictionary.H"
+#include "treeDataPoint.H"
-#include "polyMesh.H"
-#include "polyPatch.H"
-#include "volFields.H"
 #include "meshSearch.H"
+#include "Tuple2.H"
+#include "globalIndex.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -44,6 +43,25 @@ namespace Foam
         sampledTriSurfaceMesh,
         word
     );
+
+    //- Private class for finding nearest
+    //  - global index
+    //  - sqr(distance)
+    typedef Tuple2<scalar, label> nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first() < x.first())
+            {
+                x = y;
+            }
+        }
+    };
 }
 
 
@@ -66,13 +84,13 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
             "triSurface",           // local
             mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
+    pointToFace_(0)
-    faceMap_(0)
 {}
 
 
@@ -93,13 +111,13 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
             "triSurface",           // local
             mesh,                   // registry
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     ),
     needsUpdate_(true),
     cellLabels_(0),
-    pointMap_(0),
+    pointToFace_(0)
-    faceMap_(0)
 {}
 
 
@@ -128,8 +146,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
     cellLabels_.clear();
-    pointMap_.clear();
+    pointToFace_.clear();
-    faceMap_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -148,93 +165,137 @@ bool Foam::sampledTriSurfaceMesh::update()
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    cellLabels_.setSize(fc.size());
-
     meshSearch meshSearcher(mesh(), false);
 
+    const indexedOctree<treeDataPoint>& cellCentreTree =
+        meshSearcher.cellCentreTree();
+
+
+    // Global numbering for cells - only used to uniquely identify local cells.
+    globalIndex globalCells(mesh().nCells());
+    List<nearInfo> nearest(fc.size());
+    forAll(nearest, i)
+    {
+        nearest[i].first() = GREAT;
+        nearest[i].second() = labelMax;
+    }
+
+    // Search triangles using nearest on local mesh
     forAll(fc, triI)
     {
-        cellLabels_[triI] = meshSearcher.findCell
+        pointIndexHit nearInfo = cellCentreTree.findNearest
         (
             fc[triI],
-            -1,                 // seed
+            sqr(GREAT)
-            true                // use tree
         );
-    }
+        if (nearInfo.hit())
-
-
-    boolList include(surface_.size(), true);
-
-    label nNotFound = 0;
-
-    forAll(cellLabels_, triI)
         {
-        if (cellLabels_[triI] == -1)
+            nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
-        {
+            nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
-            include[triI] = false;
-            nNotFound++;
         }
     }
-    label nTotalNotFound = returnReduce(nNotFound, sumOp<label>());
 
-    if (debug)
+    // See which processor has the nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    boolList include(surface_.size(), false);
+
+    cellLabels_.setSize(fc.size());
+    cellLabels_ = -1;
+
+    label nFound = 0;
+    forAll(nearest, triI)
     {
-        Pout<< "surface:" << surface_.size()
+        if (nearest[triI].second() == labelMax)
-            << "  included:" << surface_.size()-nNotFound
+        {
-            << "  total not found" << nTotalNotFound << endl;
+            // Not found on any processor. How to map?
         }
+        else if (globalCells.isLocal(nearest[triI].second()))
+        {
+            cellLabels_[triI] = globalCells.toLocal(nearest[triI].second());
 
-
-    // Add to master all triangles are outside all meshes.
-    {
-        boolList onAnyProc(surface_.size(), false);
-        Pstream::listCombineGather(onAnyProc, orEqOp<bool>());
-
-        if (Pstream::master())
-        {
-            label nAdded = 0;
-            forAll(onAnyProc, triI)
-            {
-                if (!onAnyProc[triI])
-                {
             include[triI] = true;
-                    nAdded++;
+            nFound++;
         }
     }
 
+
     if (debug)
     {
-                Pout<< "nAdded to master:" << nAdded << endl;
+        Pout<< "Local out of faces:" << cellLabels_.size()
+            << " keeping:" << nFound << endl;
+    }
+
+    // Now subset the surface. Do not use triSurface::subsetMesh since requires
+    // original surface to be in compact numbering.
+
+    const triSurface& s = surface_;
+
+    // Compact to original triangle
+    labelList faceMap(s.size());
+    // Compact to original points
+    labelList pointMap(s.points().size());
+    // From original point to compact points
+    labelList reversePointMap(s.points().size(), -1);
+
+    {
+        label newPointI = 0;
+        label newTriI = 0;
+
+        forAll(s, triI)
+        {
+            if (include[triI])
+            {
+                faceMap[newTriI++] = triI;
+
+                const labelledTri& f = s[triI];
+                forAll(f, fp)
+                {
+                    if (reversePointMap[f[fp]] == -1)
+                    {
+                        pointMap[newPointI] = f[fp];
+                        reversePointMap[f[fp]] = newPointI++;
                     }
                 }
             }
+        }
+        faceMap.setSize(newTriI);
+        pointMap.setSize(newPointI);
+    }
 
+    // Subset cellLabels
+    cellLabels_ = UIndirectList<label>(cellLabels_, faceMap)();
 
-    // Now subset the surface
+    // Store any face per point
-    triSurface localSurface
+    pointToFace_.setSize(pointMap.size());
-    (
-        surface_.subsetMesh
-        (
-            include,
-            pointMap_,
-            faceMap_
-        )
-    );
 
-    // And convert into faces.
+    // Create faces and points for subsetted surface
     faceList& faces = this->storedFaces();
-    faces.setSize(localSurface.size());
+    faces.setSize(faceMap.size());
-    forAll(localSurface, i)
+    forAll(faceMap, i)
     {
-        faces[i] = localSurface[i].triFaceFace();
+        const triFace& f = s[faceMap[i]];
+        triFace newF
+        (
+            reversePointMap[f[0]],
+            reversePointMap[f[1]],
+            reversePointMap[f[2]]
+        );
+        faces[i] = newF.triFaceFace();
+
+        forAll(newF, fp)
+        {
+            pointToFace_[newF[fp]] = i;
+        }
     }
 
-    this->storedPoints() = localSurface.points();
+    this->storedPoints() = pointField(s.points(), pointMap);
 
     if (debug)
     {
         print(Pout);
-        Pout << endl;
+        Pout<< endl;
     }
 
     needsUpdate_ = false;

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H

@@ -29,15 +29,13 @@ Description
     A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
     of the triSurface.
 
+    It samples using the cell nearest to the triangle centre so does
+    not check the cell the centre is actually in ...
+
     In parallel every processor just operates on the part of the surface
     where the face centres are inside the mesh. It is then up to the
     caller to stitch the partial surfaces together.
 
-    No check is done for triangle centres being on multiple processors
-    (should never happen but truncation errors ...)
-
-    Any value on a triangle outside the mesh will be set to 0.
-
 SourceFiles
     sampledTriSurfaceMesh.C
 
@@ -75,14 +73,11 @@ class sampledTriSurfaceMesh
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local cell labels
+        //- From local surface triangle to mesh cell.
         labelList cellLabels_;
 
         //- From local surface back to surface_
-        labelList pointMap_;
+        labelList pointToFace_;
-
-        //- From local surface back to surface_
-        labelList faceMap_;
 
 
     // Private Member Functions

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C

@@ -36,20 +36,12 @@ Foam::sampledTriSurfaceMesh::sampleField
 ) const
 {
     // One value per face
-    tmp<Field<Type> > tvalues(new Field<Type>(faceMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(cellLabels_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(faceMap_, i)
+    forAll(cellLabels_, triI)
     {
-        label cellI = cellLabels_[faceMap_[i]];
+        values[triI] = vField[cellLabels_[triI]];
-        if (cellI != -1)
-        {
-            values[i] = vField[cellI];
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
     }
 
     return tvalues;
@@ -64,24 +56,15 @@ Foam::sampledTriSurfaceMesh::interpolateField
 ) const
 {
     // One value per vertex
-    tmp<Field<Type> > tvalues(new Field<Type>(pointMap_.size()));
+    tmp<Field<Type> > tvalues(new Field<Type>(pointToFace_.size()));
     Field<Type>& values = tvalues();
 
-    forAll(pointMap_, i)
+    forAll(pointToFace_, pointI)
     {
-        label origPointI = pointMap_[i];
+        label triI = pointToFace_[pointI];
-        label origTriI = surface_.pointFaces()[origPointI][0];
+        label cellI = cellLabels_[triI];
 
-        label cellI = cellLabels_[origTriI];
+        values[pointI] = interpolator.interpolate(points()[pointI], cellI);
-
-        if (cellI != -1)
-        {
-            values[i] = interpolator.interpolate(points()[i], cellI);
-        }
-        else
-        {
-            values[i] = pTraits<Type>::zero;
-        }
     }
 
     return tvalues;