Back again

applications/utilities/mesh/generation/extrudeMesh/Make/options

@@ -1,10 +1,14 @@
 EXE_INC = \
     -IextrudedMesh \
     -IextrudeModel/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lsurfMesh \
     -lmeshTools \
     -ldynamicMesh \
     -lextrudeModel

applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C

@@ -36,13 +36,15 @@ Description
 #include "Time.H"
 #include "dimensionedTypes.H"
 #include "IFstream.H"
-#include "faceMesh.H"
 #include "polyTopoChange.H"
 #include "polyTopoChanger.H"
 #include "edgeCollapser.H"
 #include "mathematicalConstants.H"
 #include "globalMeshData.H"
 #include "perfectInterface.H"
+#include "addPatchCellLayer.H"
+#include "fvMesh.H"
+#include "MeshedSurfaces.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -50,14 +52,148 @@ Description
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+enum ExtrudeMode
+{
+    MESH,
+    PATCH,
+    SURFACE
+};
+
+template<>
+const char* NamedEnum<ExtrudeMode, 3>::names[] =
+{
+    "mesh",
+    "patch",
+    "surface"
+};
+static const NamedEnum<ExtrudeMode, 3> ExtrudeModeNames;
+
+
+void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
+{
+    // Create dummy system/fv*
+    {
+        IOobject io
+        (
+            "fvSchemes",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        Info<< "Testing:" << io.objectPath() << endl;
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            dictionary divDict;
+            dummyDict.add("divSchemes", divDict);
+            dictionary gradDict;
+            dummyDict.add("gradSchemes", gradDict);
+            dictionary laplDict;
+            dummyDict.add("laplacianSchemes", laplDict);
+
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+    {
+        IOobject io
+        (
+            "fvSolution",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+}
+
+
+label findPatchID(const polyBoundaryMesh& patches, const word& name)
+{
+    label patchID = patches.findPatchID(name);
+
+    if (patchID == -1)
+    {
+        FatalErrorIn("findPatchID(const polyBoundaryMesh&, const word&)")
+            << "Cannot find patch " << name
+            << " in the source mesh.\n"
+            << "Valid patch names are " << patches.names()
+            << exit(FatalError);
+    }
+    return patchID;
+}
+
+
+labelList patchFaces(const polyBoundaryMesh& patches, const word& name)
+{
+    label patchID = findPatchID(patches, name);
+    const polyPatch& pp = patches[patchID];
+
+    return identity(pp.size()) + pp.start();
+}
+
+
+void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
+{
+    const labelList& reverseMap = map.reverseFaceMap();
+
+    labelList newFaceLabels(faceLabels.size());
+    label newI = 0;
+
+    forAll(faceLabels, i)
+    {
+        label oldFaceI = faceLabels[i];
+
+        if (reverseMap[oldFaceI] >= 0)
+        {
+            newFaceLabels[newI++] = reverseMap[oldFaceI];
+        }
+    }
+    newFaceLabels.setSize(newI);
+    faceLabels.transfer(newFaceLabels);
+}
+
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
+    #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTimeExtruded.H"
 
-    autoPtr<extrudedMesh> meshPtr(NULL);
+    // Get optional regionName
+    word regionName;
+    word regionDir;
+    if (args.optionReadIfPresent("region", regionName))
+    {
+        regionDir = regionName;
+        Info<< "Create mesh " << regionName << " for time = "
+            << runTimeExtruded.timeName() << nl << endl;
+    }
+    else
+    {
+        regionName = fvMesh::defaultRegion;
+        Info<< "Create mesh for time = "
+            << runTimeExtruded.timeName() << nl << endl;
+    }
+
 
     IOdictionary dict
     (
@@ -65,26 +201,44 @@ int main(int argc, char *argv[])
         (
             "extrudeProperties",
             runTimeExtruded.constant(),
+            regionDir,
             runTimeExtruded,
             IOobject::MUST_READ
         )
     );
 
+    // Point generator
     autoPtr<extrudeModel> model(extrudeModel::New(dict));
 
-    const word sourceType(dict.lookup("constructFrom"));
+    // Whether to flip normals
+    const Switch flipNormals(dict.lookup("flipNormals"));
 
-    autoPtr<faceMesh> fMesh;
+    // What to extrude
+    const ExtrudeMode mode = ExtrudeModeNames.read
+    (
+        dict.lookup("constructFrom")
+    );
 
-    if (sourceType == "patch")
+
+    // Generated mesh (one of either)
+    autoPtr<fvMesh> meshFromMesh;
+    autoPtr<polyMesh> meshFromSurface;
+
+    // Faces on front and back for stitching (in case of mergeFaces)
+    word frontPatchName;
+    labelList frontPatchFaces;
+    word backPatchName;
+    labelList backPatchFaces;
+
+    if (mode == PATCH || mode == MESH)
     {
         fileName sourceCasePath(dict.lookup("sourceCase"));
         sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
-        word patchName(dict.lookup("sourcePatch"));
+        dict.lookup("sourcePatch") >> frontPatchName;
 
-        Info<< "Extruding patch " << patchName
+        Info<< "Extruding patch " << frontPatchName
             << " on mesh " << sourceCasePath << nl
             << endl;
 
@@ -94,31 +248,183 @@ int main(int argc, char *argv[])
             sourceRootDir,
             sourceCaseDir
         );
-        #include "createPolyMesh.H"
+        #include "createMesh.H"
 
-        label patchID = mesh.boundaryMesh().findPatchID(patchName);
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        if (patchID == -1)
+
+        // Topo change container. Either copy an existing mesh or start
+        // with empty storage (number of patches only needed for checking)
+        autoPtr<polyTopoChange> meshMod
+        (
+            (
+                mode == MESH
+              ? new polyTopoChange(mesh)
+              : new polyTopoChange(patches.size())
+            )
+        );
+
+        // Extrusion engine. Either adding to existing mesh or
+        // creating separate mesh.
+        addPatchCellLayer layerExtrude(mesh, (mode == MESH));
+
+        indirectPrimitivePatch extrudePatch
+        (
+            IndirectList<face>
+            (
+                mesh.faces(),
+                patchFaces(patches, frontPatchName)
+            ),
+            mesh.points()
+        );
+
+        // Only used for addPatchCellLayer into new mesh
+        labelList exposedPatchIDs;
+        if (mode == PATCH)
         {
-            FatalErrorIn(args.executable())
-                << "Cannot find patch " << patchName
-                << " in the source mesh.\n"
-                << "Valid patch names are " << mesh.boundaryMesh().names()
-                << exit(FatalError);
+            dict.lookup("exposedPatchName") >> backPatchName;
+            exposedPatchIDs.setSize
+            (
+                extrudePatch.size(),
+                findPatchID(patches, backPatchName)
+            );
         }
 
-        const polyPatch& pp = mesh.boundaryMesh()[patchID];
-        fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints()));
 
+        pointField layer0Points(extrudePatch.nPoints());
+        pointField displacement(extrudePatch.nPoints());
+        forAll(displacement, pointI)
         {
-            fileName surfName(runTime.path()/patchName + ".sMesh");
-            Info<< "Writing patch as surfaceMesh to "
-                << surfName << nl << endl;
-            OFstream os(surfName);
-            os << fMesh() << nl;
+            const vector& patchNormal = extrudePatch.pointNormals()[pointI];
+
+            // layer0 point
+            layer0Points[pointI] = model()
+            (
+                extrudePatch.localPoints()[pointI],
+                patchNormal,
+                0
+            );
+            // layerN point
+            point extrudePt = model()
+            (
+                extrudePatch.localPoints()[pointI],
+                patchNormal,
+                model().nLayers()
+            );
+            displacement[pointI] = extrudePt - layer0Points[pointI];
         }
+
+        if (flipNormals)
+        {
+            Info<< "Flipping faces." << nl << endl;
+            displacement = -displacement;
+        }
+
+        // Check if wedge (has layer0 different from original patch points)
+        // If so move the mesh to starting position.
+        if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL)
+        {
+            Info<< "Moving mesh to layer0 points since differ from original"
+                << " points - this can happen for wedge extrusions." << nl
+                << endl;
+
+            pointField newPoints(mesh.points());
+            forAll(extrudePatch.meshPoints(), i)
+            {
+                newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i];
+            }
+            mesh.movePoints(newPoints);
+        }
+
+
+        // Layers per face
+        labelList nFaceLayers(extrudePatch.size(), model().nLayers());
+        // Layers per point
+        labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
+        // Displacement for first layer
+        vectorField firstLayerDisp = displacement*model().sumThickness(1);
+        // Expansion ratio not used.
+        scalarField ratio(extrudePatch.nPoints(), 1.0);
+
+        layerExtrude.setRefinement
+        (
+            ratio,              // expansion ratio
+            extrudePatch,       // patch faces to extrude
+            exposedPatchIDs,    // if new mesh: patches for exposed faces
+            nFaceLayers,
+            nPointLayers,
+            firstLayerDisp,
+            meshMod()
+        );
+
+        // Reset points according to extrusion model
+        forAll(layerExtrude.addedPoints(), pointI)
+        {
+            const labelList& pPoints = layerExtrude.addedPoints()[pointI];
+            forAll(pPoints, pPointI)
+            {
+                label meshPointI = pPoints[pPointI];
+
+                point modelPt
+                (
+                    model()
+                    (
+                        extrudePatch.localPoints()[pointI],
+                        extrudePatch.pointNormals()[pointI],
+                        pPointI+1       // layer
+                    )
+                );
+
+                const_cast<DynamicList<point>&>
+                (
+                    meshMod().points()
+                )[meshPointI] = modelPt;
+            }
+        }
+
+        // Store faces on front and exposed patch (if mode=patch there are
+        // only added faces so cannot used map to old faces)
+        const labelListList& layerFaces = layerExtrude.layerFaces();
+        backPatchFaces.setSize(layerFaces.size());
+        frontPatchFaces.setSize(layerFaces.size());
+        forAll(backPatchFaces, i)
+        {
+            backPatchFaces[i] = layerFaces[i][0];
+            frontPatchFaces[i] = layerFaces[i][layerFaces[i].size()-1];
+        }
+
+        // Create dummy fvSchemes, fvSolution
+        createDummyFvMeshFiles(mesh, regionDir);
+
+        // Create actual mesh from polyTopoChange container
+        autoPtr<mapPolyMesh> map = meshMod().makeMesh
+        (
+            meshFromMesh,
+            IOobject
+            (
+                regionName,
+                runTimeExtruded.constant(),
+                runTimeExtruded,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh
+        );
+
+        // Calculate face labels for front and back.
+        frontPatchFaces = renumber
+        (
+            map().reverseFaceMap(),
+            frontPatchFaces
+        );
+        backPatchFaces = renumber
+        (
+            map().reverseFaceMap(),
+            backPatchFaces
+        );
     }
-    else if (sourceType == "surface")
+    else
     {
         // Read from surface
         fileName surfName(dict.lookup("surface"));
@@ -126,52 +432,62 @@ int main(int argc, char *argv[])
         Info<< "Extruding surfaceMesh read from file " << surfName << nl
             << endl;
 
-        IFstream is(surfName);
+        MeshedSurface<face> fMesh(surfName);
 
-        fMesh.reset(new faceMesh(is));
-
-        Info<< "Read patch from file " << surfName << nl
-            << endl;
-    }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "Illegal 'constructFrom' specification. Should either be "
-            << "patch or surface." << exit(FatalError);
-    }
-
-    Switch flipNormals(dict.lookup("flipNormals"));
-
-    if (flipNormals)
-    {
-        Info<< "Flipping faces." << nl << endl;
-
-        faceList faces(fMesh().size());
-        forAll(faces, i)
+        if (flipNormals)
         {
-            faces[i] = fMesh()[i].reverseFace();
+            Info<< "Flipping faces." << nl << endl;
+            faceList& faces = const_cast<faceList&>(fMesh.faces());
+            forAll(faces, i)
+            {
+                faces[i] = fMesh[i].reverseFace();
+            }
         }
-        fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
+
+        Info<< "Extruding surface with :" << nl
+                << "    points     : " << fMesh.points().size() << nl
+                << "    faces      : " << fMesh.size() << nl
+                << "    normals[0] : " << fMesh.faceNormals()[0]
+                << nl
+                << endl;
+
+        meshFromSurface.reset
+        (
+            new extrudedMesh
+            (
+                IOobject
+                (
+                    extrudedMesh::defaultRegion,
+                    runTimeExtruded.constant(),
+                    runTimeExtruded
+                ),
+                fMesh,
+                model()
+            )
+        );
+
+
+        // Get the faces on front and back
+        frontPatchName = "originalPatch";
+        frontPatchFaces = patchFaces
+        (
+            meshFromSurface().boundaryMesh(),
+            frontPatchName
+        );
+        backPatchName = "otherSide";
+        backPatchFaces = patchFaces
+        (
+            meshFromSurface().boundaryMesh(),
+            backPatchName
+        );
     }
 
 
-    Info<< "Extruding patch with :" << nl
-            << "    points     : " << fMesh().points().size() << nl
-            << "    faces      : " << fMesh().size() << nl
-            << "    normals[0] : " << fMesh().faceNormals()[0]
-            << nl
-            << endl;
-
-    extrudedMesh mesh
+    polyMesh& mesh =
     (
-        IOobject
-        (
-            extrudedMesh::defaultRegion,
-            runTimeExtruded.constant(),
-            runTimeExtruded
-        ),
-        fMesh(),
-        model()
+        meshFromMesh.valid()
+      ? meshFromMesh()
+      : meshFromSurface()
     );
 
 
@@ -184,17 +500,6 @@ int main(int argc, char *argv[])
         << "Merge distance    : " << mergeDim << nl
         << endl;
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    const label origPatchID = patches.findPatchID("originalPatch");
-    const label otherPatchID = patches.findPatchID("otherSide");
-
-    if (origPatchID == -1 || otherPatchID == -1)
-    {
-        FatalErrorIn(args.executable())
-            << "Cannot find patch originalPatch or otherSide." << nl
-            << "Valid patches are " << patches.names() << exit(FatalError);
-    }
 
     // Collapse edges
     // ~~~~~~~~~~~~~~
@@ -237,6 +542,10 @@ int main(int argc, char *argv[])
             // Update fields
             mesh.updateMesh(map);
 
+            // Update stored data
+            updateFaceLabels(map(), frontPatchFaces);
+            updateFaceLabels(map(), backPatchFaces);
+
             // Move mesh (if inflation used)
             if (map().hasMotionPoints())
             {
@@ -252,22 +561,33 @@ int main(int argc, char *argv[])
     Switch mergeFaces(dict.lookup("mergeFaces"));
     if (mergeFaces)
     {
+        if (mode == MESH)
+        {
+            FatalErrorIn(args.executable())
+                << "Cannot stitch front and back of extrusion since"
+                << " in 'mesh' mode (extrusion appended to mesh)."
+                << exit(FatalError);
+        }
+
         Info<< "Assuming full 360 degree axisymmetric case;"
             << " stitching faces on patches "
-            << patches[origPatchID].name() << " and "
-            << patches[otherPatchID].name() << " together ..." << nl << endl;
+            << frontPatchName << " and "
+            << backPatchName << " together ..." << nl << endl;
+
+        if (frontPatchFaces.size() != backPatchFaces.size())
+        {
+            FatalErrorIn(args.executable())
+                << "Differing number of faces on front ("
+                << frontPatchFaces.size() << ") and back ("
+                << backPatchFaces.size() << ")"
+                << exit(FatalError);
+        }
+
+
 
         polyTopoChanger stitcher(mesh);
         stitcher.setSize(1);
 
-        // Make list of masterPatch faces
-        labelList isf(patches[origPatchID].size());
-
-        forAll (isf, i)
-        {
-            isf[i] = patches[origPatchID].start() + i;
-        }
-
         const word cutZoneName("originalCutFaceZone");
 
         List<faceZone*> fz
@@ -276,8 +596,8 @@ int main(int argc, char *argv[])
             new faceZone
             (
                 cutZoneName,
-                isf,
-                boolList(isf.size(), false),
+                frontPatchFaces,
+                boolList(frontPatchFaces.size(), false),
                 0,
                 mesh.faceZones()
             )
@@ -286,26 +606,58 @@ int main(int argc, char *argv[])
         mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
 
         // Add the perfect interface mesh modifier
-        stitcher.set
+        perfectInterface perfectStitcher
         (
+            "couple",
             0,
-            new perfectInterface
-            (
-                "couple",
-                0,
-                stitcher,
-                cutZoneName,
-                patches[origPatchID].name(),
-                patches[otherPatchID].name()
-            )
+            stitcher,
+            cutZoneName,
+            word::null,         // dummy patch name
+            word::null          // dummy patch name
         );
 
-        // Execute all polyMeshModifiers
-        autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
+        // Topo change container
+        polyTopoChange meshMod(mesh);
 
-        mesh.movePoints(morphMap->preMotionPoints());
+        perfectStitcher.setRefinement
+        (
+            indirectPrimitivePatch
+            (
+                IndirectList<face>
+                (
+                    mesh.faces(),
+                    frontPatchFaces
+                ),
+                mesh.points()
+            ),
+            indirectPrimitivePatch
+            (
+                IndirectList<face>
+                (
+                    mesh.faces(),
+                    backPatchFaces
+                ),
+                mesh.points()
+            ),
+            meshMod
+        );
+
+        // Construct new mesh from polyTopoChange.
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+        // Update fields
+        mesh.updateMesh(map);
+
+        // Move mesh (if inflation used)
+        if (map().hasMotionPoints())
+        {
+            mesh.movePoints(map().preMotionPoints());
+        }
     }
 
+    mesh.setInstance(runTimeExtruded.constant());
+    Info<< "Writing mesh to " << mesh.objectPath() << nl << endl;
+
     if (!mesh.write())
     {
         FatalErrorIn(args.executable()) << "Failed writing mesh"

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C

@@ -43,6 +43,7 @@ Foam::extrudeModel::extrudeModel
 )
 :
     nLayers_(readLabel(dict.lookup("nLayers"))),
+    expansionRatio_(readScalar(dict.lookup("expansionRatio"))),
     dict_(dict),
     coeffDict_(dict.subDict(modelType + "Coeffs"))
 {}
@@ -62,5 +63,28 @@ Foam::label Foam::extrudeModel::nLayers() const
 }
 
 
+Foam::scalar Foam::extrudeModel::expansionRatio() const
+{
+    return expansionRatio_;
+}
+
+
+Foam::scalar Foam::extrudeModel::sumThickness(const label layer) const
+{
+    // 1+r+r^2+ .. +r^(n-1) = (1-r^n)/(1-r)
+
+    if (mag(1.0-expansionRatio_) < SMALL)
+    {
+        return scalar(layer)/nLayers_;
+    }
+    else
+    {
+        return
+            (1.0-pow(expansionRatio_, layer))
+          / (1.0-pow(expansionRatio_, nLayers_));
+    }
+}
+
+
 // ************************************************************************* //
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H

@@ -58,6 +58,8 @@ protected:
 
         const label nLayers_;
 
+        const scalar expansionRatio_;
+
         const dictionary& dict_;
 
         const dictionary& coeffDict_;
@@ -113,9 +115,15 @@ public:
 
             label nLayers() const;
 
+            scalar expansionRatio() const;
+
 
     // Member Operators
 
+        //- Helper: calculate cumulative relative thickness for layer.
+        //  (layer=0 -> 0; layer=nLayers -> 1)
+        scalar sumThickness(const label layer) const;
+
         virtual point operator()
         (
             const point& surfacePoint,

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C

@@ -72,7 +72,8 @@ point linearNormal::operator()
     const label layer
 ) const
 {
-    scalar d = thickness_*layer/nLayers_;
+    //scalar d = thickness_*layer/nLayers_;
+    scalar d = thickness_*sumThickness(layer);
     return surfacePoint + d*surfaceNormal;
 }
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H

@@ -64,7 +64,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         linearNormal(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C

@@ -67,8 +67,7 @@ point linearRadial::operator()
     scalar rs = mag(surfacePoint);
     vector rsHat = surfacePoint/rs;
 
-    scalar delta = (R_ - rs)/nLayers_;
-    scalar r = rs + layer*delta;
+    scalar r = rs + (R_ - rs)*sumThickness(layer);
     return r*rsHat;
 }
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H

@@ -61,7 +61,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         linearRadial(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C

@@ -49,7 +49,13 @@ sigmaRadial::sigmaRadial(const dictionary& dict)
     RTbyg_(readScalar(coeffDict_.lookup("RTbyg"))),
     pRef_(readScalar(coeffDict_.lookup("pRef"))),
     pStrat_(readScalar(coeffDict_.lookup("pStrat")))
-{}
+{
+    if (mag(expansionRatio() - 1.0) > SMALL)
+    {
+        WarningIn("sigmaRadial::sigmaRadial(const dictionary&)")
+            << "Ignoring expansionRatio setting." << endl;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H

@@ -63,7 +63,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         sigmaRadial(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C

@@ -88,8 +88,8 @@ point wedge::operator()
     }
     else
     {
-        //sliceAngle = angle_*(layer + 1)/nLayers_;
-        sliceAngle = angle_*layer/nLayers_;
+        //sliceAngle = angle_*layer/nLayers_;
+        sliceAngle = angle_*sumThickness(layer);
     }
 
     // Find projection onto axis (or rather decompose surfacePoint

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H

@@ -77,7 +77,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         wedge(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeProperties

@@ -14,24 +14,28 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Where to get surface from: either from surface ('surface') or
-// from (flipped) patch of existing case ('patch')
-constructFrom patch;    //surface;
+// What to extrude:
+//      patch   : from patch of another case ('sourceCase')
+//      mesh    : as above but with original case included
+//      surface : from externally read surface
 
-// If construct from (flipped) patch
-sourceCase "$FOAM_RUN/icoFoam/cavity";
+//constructFrom mesh;
+constructFrom patch;
+//constructFrom surface;
+
+// If construct from patch/mesh:
+sourceCase "../cavity";
 sourcePatch movingWall;
+// If construct from patch: patch to use for back (can be same as sourcePatch)
+exposedPatchName movingWall;
+// If construct from surface:
+surface "movingWall.stl";
+
+
 
 // Flip surface normals before usage.
 flipNormals false;
 
-// If construct from surface
-surface "movingWall.sMesh";
-
-
-// Do front and back need to be merged? Usually only makes sense for 360
-// degree wedges.
-mergeFaces true;
 
 
 //- Linear extrusion in point-normal direction
@@ -46,11 +50,13 @@ extrudeModel        wedge;
 //- Extrudes into sphere with grading according to pressure (atmospherics)
 //extrudeModel        sigmaRadial;
 
-nLayers             20;
+nLayers             10;
+
+expansionRatio      1.0;    //0.9;
 
 wedgeCoeffs
 {
-    axisPt      (0 0.1 0);
+    axisPt      (0 0.1 -0.05);
     axis        (-1 0 0);
     angle       360;  // For nLayers=1 assume symmetry so angle/2 on each side
 }
@@ -73,5 +79,10 @@ sigmaRadialCoeffs
 }
 
 
+// Do front and back need to be merged? Usually only makes sense for 360
+// degree wedges.
+mergeFaces true;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -208,8 +208,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
             // Dump halves
             {
                 OFstream str(prefix+cycPatch.name()+"_half0.obj");
-                Pout<< "Dumping cycPatch.name() half0 faces to " << str.name()
-                    << endl;
+                Pout<< "Dumping " << cycPatch.name()
+                    << " half0 faces to " << str.name() << endl;
                 meshTools::writeOBJ
                 (
                     str,
@@ -226,8 +226,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
             }
             {
                 OFstream str(prefix+cycPatch.name()+"_half1.obj");
-                Pout<< "Dumping cycPatch.name() half1 faces to " << str.name()
-                    << endl;
+                Pout<< "Dumping " << cycPatch.name()
+                    << " half1 faces to " << str.name() << endl;
                 meshTools::writeOBJ
                 (
                     str,
@@ -563,7 +563,7 @@ int main(int argc, char *argv[])
     dumpCyclicMatch("initial_", mesh);
 
     // Read patch construct info from dictionary
-    PtrList<dictionary> patchSources(dict.lookup("patches"));
+    PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
 
 
 

applications/utilities/mesh/manipulation/createPatch/createPatchDict

@@ -46,7 +46,7 @@ matchTolerance 1E-3;
 pointSync true;
 
 // Patches to create.
-patches
+patchInfo
 (
     {
         // Name of new patch

applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

@@ -105,6 +105,8 @@ int main(int argc, char *argv[])
 
     Foam::argList::validOptions.insert("overwrite", "");
 
+    Foam::argList::validOptions.insert("toleranceDict", "file with tolerances");
+
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
@@ -168,6 +170,22 @@ int main(int argc, char *argv[])
             << "If this is not the case use the -partial option" << nl << endl;
     }
 
+    // set up the tolerances for the sliding mesh
+    dictionary slidingTolerances;
+    if (args.options().found("toleranceDict")) 
+    {
+        IOdictionary toleranceFile(
+            IOobject(
+                args.options()["toleranceDict"],
+                runTime.constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE                 
+            )
+        );
+        slidingTolerances += toleranceFile;
+    }
+
     // Check for non-empty master and slave patches
     checkPatch(mesh.boundaryMesh(), masterPatchName);
     checkPatch(mesh.boundaryMesh(), slavePatchName);
@@ -320,6 +338,11 @@ int main(int argc, char *argv[])
                 true                    // couple/decouple mode
             )
         );
+        static_cast<slidingInterface&>(stitcher[0]).setTolerances
+        (
+            slidingTolerances,
+            true
+        );
     }
 
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options

@@ -5,5 +5,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian
 

applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options

@@ -6,4 +6,5 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
+    -lgenericPatchFields \
     -lconversion

applications/utilities/postProcessing/dataConversion/foamToFieldview9/Make/options

@@ -4,4 +4,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian

applications/utilities/postProcessing/dataConversion/foamToGMV/Make/options

@@ -14,4 +14,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian 

applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options

@@ -6,5 +6,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
+    -lgenericPatchFields \
     -lmeshTools
 

applications/utilities/postProcessing/dataConversion/smapToFoam/Make/options

@@ -2,4 +2,5 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lgenericPatchFields \
     -lfiniteVolume

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -22,6 +22,7 @@ FoamFile
 //      jplot
 //      gnuplot
 //      raw
+//      vtk
 setFormat raw;
 
 // Surface output format. Choice of
@@ -62,6 +63,7 @@ fields
 //      curve               specified points, not nessecary on line, uses
 //                          tracking
 //      cloud               specified points, uses findCell
+//      triSurfaceMeshPointSet  points of triSurface
 //
 // axis: how to write point coordinate. Choice of
 // - x/y/z: x/y/z coordinate only

bin/foamClearPolyMesh

@@ -130,6 +130,22 @@ for i in \
     pointLevel \
     refinementHistory \
     surfaceIndex \
+    points.gz \
+    faces.gz \
+    owner.gz \
+    neighbour.gz \
+    cells.gz \
+    boundary.gz \
+    pointZones.gz \
+    faceZones.gz \
+    cellZones.gz \
+    meshModifiers.gz \
+    parallelData.gz \
+    sets.gz \
+    cellLevel.gz \
+    pointLevel.gz \
+    refinementHistory.gz \
+    surfaceIndex.gz \
 ;
 do
     rm -rf $meshDir/$i

bin/tools/CleanFunctions

@@ -72,6 +72,7 @@ cleanCase ()
            constant/{cellToRegion,cellLevel*,pointLevel*} \
            constant/polyMesh/sets/ \
            VTK \
+           sets/streamLines \
            > /dev/null 2>&1
 
     for f in `find . -name "*Dict"`

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

@@ -117,12 +117,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking receives: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking receives: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Return next token from stream
             Istream& read(token&);
 

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

@@ -115,12 +115,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking writes: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking writes: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Write next token to stream
             Ostream& write(const token&);
 

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

@@ -264,6 +264,12 @@ public:
         //  Spawns slave processes and initialises inter-communication
         static bool init(int& argc, char**& argv);
 
+        //- Non-blocking comms: wait until all have finished.
+        static void waitRequests();
+
+        //- Non-blocking comms: has request i finished?
+        static bool finishedRequest(const label i);
+
         //- Is this a parallel run?
         static bool parRun()
         {

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C

@@ -300,8 +300,7 @@ evaluate()
         // Block for any outstanding requests
         if (Pstream::defaultCommsType == Pstream::nonBlocking)
         {
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
         }
 
         forAll(*this, patchi)

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -339,8 +339,7 @@ void Foam::mapDistribute::compact(const boolList& elemIsUsed)
 
         // Wait for all to finish
 
-        OPstream::waitRequests();
-        IPstream::waitRequests();
+        Pstream::waitRequests();
 
 
         // Compact out all submap entries that are referring to unused elements

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -287,8 +287,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -413,8 +412,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
-            IPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Collect neighbour fields
 
@@ -717,8 +715,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -842,8 +839,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
-            IPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Collect neighbour fields
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -1305,7 +1305,7 @@ bool Foam::cyclicPolyPatch::order
     {
         label baffleI = 0;
 
-        forAll(*this, faceI)
+        forAll(pp, faceI)
         {
             const face& f = pp.localFaces()[faceI];
             const labelList& pFaces = pp.pointFaces()[f[0]];

src/Pstream/dummy/IPread.C

@@ -86,17 +86,4 @@ int Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{}
-
-
-bool Foam::IPstream::finishedRequest(const label)
-{
-    notImplemented("IPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/OPwrite.C

@@ -65,17 +65,4 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{}
-
-
-bool Foam::OPstream::finishedRequest(const label)
-{
-    notImplemented("OPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/Pstream.C

@@ -61,6 +61,17 @@ void Foam::Pstream::abort()
 void Foam::reduce(scalar&, const sumOp<scalar>&)
 {}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+void Foam::Pstream::waitRequests()
+{}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    notImplemented("Pstream::finishedRequest()");
+    return false;
+}
+
 
 // ************************************************************************* //

src/Pstream/mpi/IPread.C

@@ -30,6 +30,7 @@ Description
 #include "mpi.h"
 
 #include "IPstream.H"
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,7 +38,7 @@ Description
 
 // Outstanding non-blocking operations.
 //! @cond fileScope
-Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
+//Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
 //! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
@@ -185,7 +186,7 @@ Foam::label Foam::IPstream::read
             return 0;
         }
 
-        IPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
 
         // Assume the message is completely received.
         return bufSize;
@@ -204,52 +205,6 @@ Foam::label Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{
-    if (IPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                IPstream_outstandingRequests_.size(),
-                IPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "IPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << endl;
-        }
-
-        IPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::IPstream::finishedRequest(const label i)
-{
-    if (i >= IPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "IPstream::finishedRequest(const label)"
-        )   << "There are " << IPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&IPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Make/files

@@ -1,5 +1,6 @@
 OPwrite.C
 IPread.C
 Pstream.C
+PstreamGlobals.C
 
 LIB = $(FOAM_MPI_LIBBIN)/libPstream

src/Pstream/mpi/OPwrite.C

@@ -30,13 +30,7 @@ Description
 #include "mpi.h"
 
 #include "OPstream.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// Outstanding non-blocking operations.
-//! @cond fileScope
-Foam::DynamicList<MPI_Request> OPstream_outstandingRequests_;
-//! @endcond fileScope
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
@@ -126,7 +120,7 @@ bool Foam::OPstream::write
             &request
         );
 
-        OPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
     }
     else
     {
@@ -142,52 +136,6 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{
-    if (OPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                OPstream_outstandingRequests_.size(),
-                OPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "OPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << Foam::endl;
-        }
-
-        OPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::OPstream::finishedRequest(const label i)
-{
-    if (i >= OPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "OPstream::finishedRequest(const label)"
-        )   << "There are " << OPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&OPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Pstream.C

@@ -29,6 +29,7 @@ License
 #include "Pstream.H"
 #include "PstreamReduceOps.H"
 #include "OSspecific.H"
+#include "PstreamGlobals.H"
 
 #include <cstring>
 #include <cstdlib>
@@ -130,6 +131,19 @@ void Foam::Pstream::exit(int errnum)
     delete[] buff;
 #   endif
 
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        label n = PstreamGlobals::outstandingRequests_.size();
+        PstreamGlobals::outstandingRequests_.clear();
+
+        WarningIn("Pstream::exit(int)")
+            << "There are still " << n << " outstanding MPI_Requests." << endl
+            << "This means that your code exited before doing a"
+            << " Pstream::waitRequests()." << endl
+            << "This should not happen for a normal code exit."
+            << endl;
+    }
+
     if (errnum == 0)
     {
         MPI_Finalize();
@@ -422,6 +436,57 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
 }
 
 
+void Foam::Pstream::waitRequests()
+{
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        if
+        (
+            MPI_Waitall
+            (
+                PstreamGlobals::outstandingRequests_.size(),
+                PstreamGlobals::outstandingRequests_.begin(),
+                MPI_STATUSES_IGNORE
+            )
+        )
+        {
+            FatalErrorIn
+            (
+                "Pstream::waitRequests()"
+            )   << "MPI_Waitall returned with error" << Foam::endl;
+        }
+
+        PstreamGlobals::outstandingRequests_.clear();
+    }
+}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    if (i >= PstreamGlobals::outstandingRequests_.size())
+    {
+        FatalErrorIn
+        (
+            "Pstream::finishedRequest(const label)"
+        )   << "There are " << PstreamGlobals::outstandingRequests_.size()
+            << " outstanding send requests and you are asking for i=" << i
+            << nl
+            << "Maybe you are mixing blocking/non-blocking comms?"
+            << Foam::abort(FatalError);
+    }
+
+    int flag;
+    MPI_Test
+    (
+       &PstreamGlobals::outstandingRequests_[i],
+       &flag,
+        MPI_STATUS_IGNORE
+    );
+
+    return flag != 0;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "PstreamGlobals.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// Outstanding non-blocking operations.
+//! @cond fileScope
+DynamicList<MPI_Request> PstreamGlobals::outstandingRequests_;
+//! @endcond fileScope
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Namespace
+    Foam::PstreamGlobals
+
+Description
+    Global functions and variables for working with parallel streams,
+    but principally for gamma/mpi
+
+SourceFiles
+    PstreamGlobals.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PstreamGlobals_H
+#define PstreamGlobals_H
+
+#include "mpi.h"
+
+#include "DynamicList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class PstreamGlobals Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace PstreamGlobals
+{
+
+extern DynamicList<MPI_Request> outstandingRequests_;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -2987,6 +2987,7 @@ void Foam::autoLayerDriver::addLayers
         (
             invExpansionRatio,
             pp,
+            labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
             firstDisp,          // thickness of layer nearest internal mesh

src/dynamicMesh/perfectInterface/perfectInterface.C

@@ -38,6 +38,7 @@ Description
 #include "polyModifyFace.H"
 #include "polyRemovePoint.H"
 #include "polyRemoveFace.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -61,7 +62,7 @@ const Foam::scalar Foam::perfectInterface::tol_ = 1E-3;
 
 Foam::pointField Foam::perfectInterface::calcFaceCentres
 (
-    const primitivePatch& pp
+    const indirectPrimitivePatch& pp
 )
 {
     const pointField& points = pp.points();
@@ -155,6 +156,295 @@ bool Foam::perfectInterface::changeTopology() const
 }
 
 
+void Foam::perfectInterface::setRefinement
+(
+    const indirectPrimitivePatch& pp0,
+    const indirectPrimitivePatch& pp1,
+    polyTopoChange& ref
+) const
+{
+    const polyMesh& mesh = topoChanger().mesh();
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Some aliases
+    const edgeList& edges0 = pp0.edges();
+    const pointField& pts0 = pp0.localPoints();
+    const pointField& pts1 = pp1.localPoints();    
+    const labelList& meshPts0 = pp0.meshPoints();
+    const labelList& meshPts1 = pp1.meshPoints();
+
+
+    // Get local dimension as fraction of minimum edge length
+
+    scalar minLen = GREAT;
+
+    forAll(edges0, edgeI)
+    {
+        minLen = min(minLen, edges0[edgeI].mag(pts0));
+    }
+    scalar typDim = tol_*minLen;
+
+    if (debug)
+    {
+        Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
+            << " pts0:" << pts0.size() << " pts1:" << pts1.size()
+            << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
+    }
+
+
+    // Determine pointMapping in mesh point labels. Uses geometric
+    // comparison to find correspondence between patch points.
+
+    labelList renumberPoints(mesh.points().size());
+    forAll(renumberPoints, i)
+    {
+        renumberPoints[i] = i;
+    }
+    {
+        labelList from1To0Points(pts1.size());
+
+        bool matchOk = matchPoints
+        (
+            pts1,
+            pts0,
+            scalarField(pts1.size(), typDim),   // tolerance
+            true,                               // verbose
+            from1To0Points
+        );
+
+        if (!matchOk)
+        {
+            FatalErrorIn
+            (
+                "perfectInterface::setRefinement(polyTopoChange& ref) const"
+            )   << "Points on patch sides do not match to within tolerance "
+                << typDim << exit(FatalError);
+        }
+
+        forAll(pts1, i)
+        {
+            renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
+        }
+    }
+
+
+
+    // Calculate correspondence between patch faces
+
+    labelList from0To1Faces(pp1.size());
+
+    bool matchOk = matchPoints
+    (
+        calcFaceCentres(pp0),
+        calcFaceCentres(pp1),
+        scalarField(pp0.size(), typDim),    // tolerance
+        true,                               // verbose
+        from0To1Faces
+    );
+
+    if (!matchOk)
+    {
+        FatalErrorIn
+        (
+            "perfectInterface::setRefinement(polyTopoChange& ref) const"
+        )   << "Face centres of patch sides do not match to within tolerance "
+            << typDim << exit(FatalError);
+    }
+
+
+
+    // Now
+    // - renumber faces using pts1 (except patch1 faces)
+    // - remove patch1 faces. Remember cell label on owner side.
+    // - modify patch0 faces to be internal.
+
+    // 1. Get faces to be renumbered
+    labelHashSet affectedFaces(2*pp1.size());
+    forAll(meshPts1, i)
+    {
+        label meshPointI = meshPts1[i];
+
+        if (meshPointI != renumberPoints[meshPointI])
+        {
+            const labelList& pFaces = mesh.pointFaces()[meshPointI];
+
+            forAll(pFaces, pFaceI)
+            {
+                affectedFaces.insert(pFaces[pFaceI]);
+            }
+        }
+    }
+    forAll(pp1, i)
+    {
+        affectedFaces.erase(pp1.addressing()[i]);
+    }
+    // Remove patch0 from renumbered faces. Should not be nessecary since
+    // patch0 and 1 should not share any point (if created by mergeMeshing)
+    // so affectedFaces should not contain any patch0 faces but you can
+    // never be sure what the user is doing.
+    forAll(pp0, i)
+    {
+        label faceI = pp0.addressing()[i];
+
+        if (affectedFaces.erase(faceI))
+        {
+            WarningIn
+            (
+                "perfectInterface::setRefinement(polyTopoChange&) const"
+            )   << "Found face " << faceI << " vertices "
+                << mesh.faces()[faceI] << " whose points are"
+                << " used both by master patch and slave patch" << endl;
+        }
+    }
+
+
+    // 2. Renumber (non patch0/1) faces.
+    for
+    (
+        labelHashSet::const_iterator iter = affectedFaces.begin();
+        iter != affectedFaces.end();
+        ++iter
+    )
+    {
+        label faceI = iter.key();
+
+        const face& f = mesh.faces()[faceI];
+
+        face newFace(f.size());
+
+        forAll(newFace, fp)
+        {
+            newFace[fp] = renumberPoints[f[fp]];
+        }
+
+        label nbr = -1;
+
+        label patchI = -1;
+
+        if (mesh.isInternalFace(faceI))
+        {
+            nbr = mesh.faceNeighbour()[faceI];
+        }
+        else
+        {
+            patchI = patches.whichPatch(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)];
+        }
+
+        ref.setAction
+        (
+            polyModifyFace
+            (
+                newFace,                    // modified face
+                faceI,                      // label of face being modified
+                mesh.faceOwner()[faceI],    // owner
+                nbr,                        // neighbour
+                false,                      // face flip
+                patchI,                     // patch for face
+                false,                      // remove from zone
+                zoneID,                     // zone for face
+                zoneFlip                    // face flip in zone
+            )
+        );
+    }
+
+
+    // 3. Remove patch1 points
+    forAll(meshPts1, i)
+    {
+        label meshPointI = meshPts1[i];
+
+        if (meshPointI != renumberPoints[meshPointI])
+        {
+            ref.setAction(polyRemovePoint(meshPointI));
+        }
+    }
+
+
+    // 4. Remove patch1 faces
+    forAll(pp1, i)
+    {
+        label faceI = pp1.addressing()[i];
+        ref.setAction(polyRemoveFace(faceI));
+    }
+
+
+    // 5. Modify patch0 faces for new points (not really nessecary; see
+    // comment above about patch1 and patch0 never sharing points) and
+    // becoming internal.
+    const boolList& mfFlip =
+        mesh.faceZones()[faceZoneID_.index()].flipMap();
+
+    forAll(pp0, i)
+    {
+        label faceI = pp0.addressing()[i];
+
+        const face& f = mesh.faces()[faceI];
+
+        face newFace(f.size());
+
+        forAll(newFace, fp)
+        {
+            newFace[fp] = renumberPoints[f[fp]];
+        }
+
+        label own = mesh.faceOwner()[faceI];
+
+        label pp1FaceI = pp1.addressing()[from0To1Faces[i]];
+
+        label nbr = mesh.faceOwner()[pp1FaceI];
+
+        if (own < nbr)
+        {
+            ref.setAction
+            (
+                polyModifyFace
+                (
+                    newFace,                // modified face
+                    faceI,                  // label of face being modified
+                    own,                    // owner
+                    nbr,                    // neighbour
+                    false,                  // face flip
+                    -1,                     // patch for face
+                    false,                  // remove from zone
+                    faceZoneID_.index(),    // zone for face
+                    mfFlip[i]               // face flip in zone
+                )
+            );
+        }
+        else
+        {
+            ref.setAction
+            (
+                polyModifyFace
+                (
+                    newFace.reverseFace(),  // modified face
+                    faceI,                  // label of face being modified
+                    nbr,                    // owner
+                    own,                    // neighbour
+                    true,                   // face flip
+                    -1,                     // patch for face
+                    false,                  // remove from zone
+                    faceZoneID_.index(),    // zone for face
+                    !mfFlip[i]              // face flip in zone
+                )
+            );
+        }
+    }
+}
+
+
 void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
 {
     if (debug)
@@ -176,286 +466,25 @@ void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
         const polyMesh& mesh = topoChanger().mesh();
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        const polyPatch& pp0 = patches[masterPatchID_.index()];
-        const polyPatch& pp1 = patches[slavePatchID_.index()];
-
-        // Some aliases
-        const edgeList& edges0 = pp0.edges();
-        const pointField& pts0 = pp0.localPoints();
-        const pointField& pts1 = pp1.localPoints();    
-        const labelList& meshPts0 = pp0.meshPoints();
-        const labelList& meshPts1 = pp1.meshPoints();
+        const polyPatch& patch0 = patches[masterPatchID_.index()];
+        const polyPatch& patch1 = patches[slavePatchID_.index()];
 
 
-        // Get local dimension as fraction of minimum edge length
-
-        scalar minLen = GREAT;
-
-        forAll(edges0, edgeI)
-        {
-            minLen = min(minLen, edges0[edgeI].mag(pts0));
-        }
-        scalar typDim = tol_*minLen;
-
-        if (debug)
-        {
-            Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
-                << " pts0:" << pts0.size() << " pts1:" << pts1.size()
-                << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
-        }
-
-
-        // Determine pointMapping in mesh point labels. Uses geometric
-        // comparison to find correspondence between patch points.
-
-        labelList renumberPoints(mesh.points().size());
-        forAll(renumberPoints, i)
-        {
-            renumberPoints[i] = i;
-        }
-        {
-            labelList from1To0Points(pts1.size());
-
-            bool matchOk = matchPoints
-            (
-                pts1,
-                pts0,
-                scalarField(pts1.size(), typDim),   // tolerance
-                true,                               // verbose
-                from1To0Points
-            );
-
-            if (!matchOk)
-            {
-                FatalErrorIn
-                (
-                    "perfectInterface::setRefinement(polyTopoChange& ref) const"
-                )   << "Points on patches " << pp0.name() << " and "
-                    << pp1.name() << " do not match to within tolerance "
-                    << typDim << exit(FatalError);
-            }
-
-            forAll(pts1, i)
-            {
-                renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
-            }
-        }
-
-
-
-        // Calculate correspondence between patch faces
-
-        labelList from0To1Faces(pp1.size());
-
-        bool matchOk = matchPoints
+        labelList pp0Labels(identity(patch0.size())+patch0.start());
+        indirectPrimitivePatch pp0
         (
-            calcFaceCentres(pp0),
-            calcFaceCentres(pp1),
-            scalarField(pp0.size(), typDim),    // tolerance
-            true,                               // verbose
-            from0To1Faces
+            IndirectList<face>(mesh.faces(), pp0Labels),
+            mesh.points()
         );
 
-        if (!matchOk)
-        {
-            FatalErrorIn
-            (
-                "perfectInterface::setRefinement(polyTopoChange& ref) const"
-            )   << "Face centres of patches " << pp0.name() << " and "
-                << pp1.name() << " do not match to within tolerance " << typDim
-                << exit(FatalError);
-        }
-
-
-
-        // Now
-        // - renumber faces using pts1 (except patch1 faces)
-        // - remove patch1 faces. Remember cell label on owner side.
-        // - modify patch0 faces to be internal.
-
-        // 1. Get faces to be renumbered
-        labelHashSet affectedFaces(2*pp1.size());
-        forAll(meshPts1, i)
-        {
-            label meshPointI = meshPts1[i];
-
-            if (meshPointI != renumberPoints[meshPointI])
-            {
-                const labelList& pFaces = mesh.pointFaces()[meshPointI];
-
-                forAll(pFaces, pFaceI)
-                {
-                    affectedFaces.insert(pFaces[pFaceI]);
-                }
-            }
-        }
-        forAll(pp1, i)
-        {
-            affectedFaces.erase(pp1.start() + i);
-        }
-        // Remove patch0 from renumbered faces. Should not be nessecary since
-        // patch0 and 1 should not share any point (if created by mergeMeshing)
-        // so affectedFaces should not contain any patch0 faces but you can
-        // never be sure what the user is doing.
-        forAll(pp0, i)
-        {
-            if (affectedFaces.erase(pp0.start() + i))
-            {
-                WarningIn
-                (
-                    "perfectInterface::setRefinement(polyTopoChange&) const"
-                )   << "Found face " << pp0.start() + i << " vertices "
-                    << mesh.faces()[pp0.start() + i] << " whose points are"
-                    << " used both by master patch " << pp0.name()
-                    << " and slave patch " << pp1.name()
-                    << endl;
-            }
-        }
-
-
-        // 2. Renumber (non patch0/1) faces.
-        for
+        labelList pp1Labels(identity(patch1.size())+patch1.start());
+        indirectPrimitivePatch pp1
         (
-            labelHashSet::const_iterator iter = affectedFaces.begin();
-            iter != affectedFaces.end();
-            ++iter
-        )
-        {
-            label faceI = iter.key();
+            IndirectList<face>(mesh.faces(), pp1Labels),
+            mesh.points()
+        );
 
-            const face& f = mesh.faces()[faceI];
-
-            face newFace(f.size());
-
-            forAll(newFace, fp)
-            {
-                newFace[fp] = renumberPoints[f[fp]];
-            }
-
-            label nbr = -1;
-
-            label patchI = -1;
-
-            if (mesh.isInternalFace(faceI))
-            {
-                nbr = mesh.faceNeighbour()[faceI];
-            }
-            else
-            {
-                patchI = patches.whichPatch(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)];
-            }
-
-            ref.setAction
-            (
-                polyModifyFace
-                (
-                    newFace,                    // modified face
-                    faceI,                      // label of face being modified
-                    mesh.faceOwner()[faceI],    // owner
-                    nbr,                        // neighbour
-                    false,                      // face flip
-                    patchI,                     // patch for face
-                    false,                      // remove from zone
-                    zoneID,                     // zone for face
-                    zoneFlip                    // face flip in zone
-                )
-            );
-        }
-
-
-        // 3. Remove patch1 points
-        forAll(meshPts1, i)
-        {
-            label meshPointI = meshPts1[i];
-
-            if (meshPointI != renumberPoints[meshPointI])
-            {
-                ref.setAction(polyRemovePoint(meshPointI));
-            }
-        }
-
-
-        // 4. Remove patch1 faces
-        forAll(pp1, i)
-        {
-            ref.setAction(polyRemoveFace(pp1.start() + i));
-        }
-
-
-        // 5. Modify patch0 faces for new points (not really nessecary; see
-        // comment above about patch1 and patch0 never sharing points) and
-        // becoming internal.
-        const boolList& mfFlip =
-            mesh.faceZones()[faceZoneID_.index()].flipMap();
-
-        forAll(pp0, i)
-        {
-            label faceI = pp0.start() + i;
-
-            const face& f = mesh.faces()[faceI];
-
-            face newFace(f.size());
-
-            forAll(newFace, fp)
-            {
-                newFace[fp] = renumberPoints[f[fp]];
-            }
-
-            label own = mesh.faceOwner()[faceI];
-
-            label pp1FaceI = pp1.start() + from0To1Faces[i];
-
-            label nbr = mesh.faceOwner()[pp1FaceI];
-
-            if (own < nbr)
-            {
-                ref.setAction
-                (
-                    polyModifyFace
-                    (
-                        newFace,                // modified face
-                        faceI,                  // label of face being modified
-                        own,                    // owner
-                        nbr,                    // neighbour
-                        false,                  // face flip
-                        -1,                     // patch for face
-                        false,                  // remove from zone
-                        faceZoneID_.index(),    // zone for face
-                        mfFlip[i]               // face flip in zone
-                    )
-                );
-            }
-            else
-            {
-                ref.setAction
-                (
-                    polyModifyFace
-                    (
-                        newFace.reverseFace(),  // modified face
-                        faceI,                  // label of face being modified
-                        nbr,                    // owner
-                        own,                    // neighbour
-                        false,                  // face flip
-                        -1,                     // patch for face
-                        false,                  // remove from zone
-                        faceZoneID_.index(),    // zone for face
-                        !mfFlip[i]              // face flip in zone
-                    )
-                );
-            }
-        }
+        setRefinement(pp0, pp1, ref);
     }
 }
 

src/dynamicMesh/perfectInterface/perfectInterface.H

@@ -40,6 +40,7 @@ SourceFiles
 #include "polyMeshModifier.H"
 #include "polyPatchID.H"
 #include "ZoneIDs.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,7 +76,7 @@ class perfectInterface
     // Private Member Functions
 
         //- Calculate face centres on patch
-        static pointField calcFaceCentres(const primitivePatch&);
+        static pointField calcFaceCentres(const indirectPrimitivePatch&);
 
 
         //- Disallow default bitwise copy construct
@@ -128,6 +129,17 @@ public:
         //  into the topological change
         virtual void setRefinement(polyTopoChange&) const;
 
+        //- Insert the layer addition/removal instructions
+        //  into the topological change. Uses only mesh, not any of the
+        //  patch and zone indices. Bit of a workaround - used in extruding
+        //  a mesh.
+        virtual void setRefinement
+        (
+            const indirectPrimitivePatch& pp0,
+            const indirectPrimitivePatch& pp1,
+            polyTopoChange&
+        ) const;
+
         //- Modify motion points to comply with the topological change
         virtual void modifyMotionPoints(pointField& motionPoints) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -335,17 +335,19 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     label inflateEdgeI = -1;
 
     // Check mesh faces using edge
-    forAll(meshFaces, i)
+    if (addToMesh_)
     {
-        if (mesh_.isInternalFace(meshFaces[i]))
+        forAll(meshFaces, i)
         {
-            // meshEdge uses internal faces so ok to inflate from it
-            inflateEdgeI = meshEdgeI;
-            break;
+            if (mesh_.isInternalFace(meshFaces[i]))
+            {
+                // meshEdge uses internal faces so ok to inflate from it
+                inflateEdgeI = meshEdgeI;
+                break;
+            }
         }
     }
 
-
     // Get my mesh face and its zone.
     label meshFaceI = pp.addressing()[ownFaceI];
     label zoneI = mesh_.faceZones().whichZone(meshFaceI);
@@ -504,9 +506,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
-Foam::addPatchCellLayer::addPatchCellLayer(const polyMesh& mesh)
+Foam::addPatchCellLayer::addPatchCellLayer
+(
+    const polyMesh& mesh,
+    const bool addToMesh
+)
 :
     mesh_(mesh),
+    addToMesh_(addToMesh),
     addedPoints_(0),
     layerFaces_(0)
 {}
@@ -551,6 +558,7 @@ void Foam::addPatchCellLayer::setRefinement
 (
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
+    const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
     const vectorField& firstLayerDisp,
@@ -827,6 +835,19 @@ void Foam::addPatchCellLayer::setRefinement
     }
 
 
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // Precalculated patchID for each patch face
+    labelList patchID(pp.size());
+
+    forAll(pp, patchFaceI)
+    {
+        label meshFaceI = pp.addressing()[patchFaceI];
+
+        patchID[patchFaceI] = patches.whichPatch(meshFaceI);
+    }
+
+
     // From master point (in patch point label) to added points (in mesh point
     // label)
     addedPoints_.setSize(pp.nPoints());
@@ -857,6 +878,33 @@ void Foam::addPatchCellLayer::setRefinement
     // Create new points
     //
 
+    // If creating new mesh: copy existing patch points
+    labelList copiedPatchPoints;
+    if (!addToMesh_)
+    {
+        copiedPatchPoints.setSize(firstLayerDisp.size());
+        forAll(firstLayerDisp, patchPointI)
+        {
+            if (addedPoints_[patchPointI].size())
+            {
+                label meshPointI = meshPoints[patchPointI];
+                label zoneI = mesh_.pointZones().whichZone(meshPointI);
+                copiedPatchPoints[patchPointI] = meshMod.setAction
+                (
+                    polyAddPoint
+                    (
+                        mesh_.points()[meshPointI],         // point
+                        -1,         // master point
+                        zoneI,      // zone for point
+                        true        // supports a cell
+                    )
+                );
+            }
+        }
+    }
+
+
+    // Create points for additional layers
     forAll(firstLayerDisp, patchPointI)
     {
         if (addedPoints_[patchPointI].size())
@@ -878,7 +926,7 @@ void Foam::addPatchCellLayer::setRefinement
                     polyAddPoint
                     (
                         pt,         // point
-                        meshPointI, // master point
+                        (addToMesh_ ? meshPointI : -1), // master point
                         zoneI,      // zone for point
                         true        // supports a cell
                     )
@@ -922,34 +970,15 @@ void Foam::addPatchCellLayer::setRefinement
                         -1,             // master point
                         -1,             // master edge
                         -1,             // master face
-                        mesh_.faceOwner()[meshFaceI],   // master cell id
+                        (addToMesh_ ? mesh_.faceOwner()[meshFaceI] : -1),//master
                         ownZoneI        // zone for cell
                     )
                 );
-
-                //Pout<< "For patchFace:" << patchFaceI
-                //    << " meshFace:" << pp.addressing()[patchFaceI]
-                //    << " layer:" << i << " added cell:"
-                //    << addedCells[patchFaceI][i]
-                //    << endl;
             }
         }
     }
 
 
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Precalculated patchID for each patch face
-    labelList patchID(pp.size());
-
-    forAll(pp, patchFaceI)
-    {
-        label meshFaceI = pp.addressing()[patchFaceI];
-
-        patchID[patchFaceI] = patches.whichPatch(meshFaceI);
-    }
-
-
 
     // Create faces on top of the original patch faces.
     // These faces are created from original patch faces outwards so in order
@@ -965,7 +994,6 @@ void Foam::addPatchCellLayer::setRefinement
         if (addedCells[patchFaceI].size())
         {
             layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1);
-            layerFaces_[patchFaceI][0] = meshFaceI;
 
             label zoneI = mesh_.faceZones().whichZone(meshFaceI);
 
@@ -981,7 +1009,12 @@ void Foam::addPatchCellLayer::setRefinement
                     if (addedPoints_[f[fp]].empty())
                     {
                         // Keep original point
-                        newFace[fp] = meshPoints[f[fp]];
+                        newFace[fp] =
+                        (
+                            addToMesh_
+                          ? meshPoints[f[fp]]
+                          : copiedPatchPoints[f[fp]]
+                        );
                     }
                     else
                     {
@@ -1021,18 +1054,13 @@ void Foam::addPatchCellLayer::setRefinement
                         nei,                        // neighbour
                         -1,                         // master point
                         -1,                         // master edge
-                        meshFaceI,                  // master face for addition
+                        (addToMesh_ ? meshFaceI : -1), // master face
                         false,                      // flux flip
                         patchI,                     // patch for face
                         zoneI,                      // zone for face
                         false                       // face zone flip
                     )
                 );
-                //Pout<< "Added inbetween face " << newFace
-                //    << " own:" << addedCells[patchFaceI][i]
-                //    << " nei:" << nei
-                //    << " patch:" << patchI
-                //    << endl;
             }
         }
     }
@@ -1040,35 +1068,81 @@ void Foam::addPatchCellLayer::setRefinement
     //
     // Modify old patch faces to be on the inside
     //
-    forAll(pp, patchFaceI)
+
+    if (addToMesh_)
     {
-        if (addedCells[patchFaceI].size())
+        forAll(pp, patchFaceI)
         {
-            label meshFaceI = pp.addressing()[patchFaceI];
+            if (addedCells[patchFaceI].size())
+            {
+                label meshFaceI = pp.addressing()[patchFaceI];
 
-            label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+                layerFaces_[patchFaceI][0] = meshFaceI;
 
-            meshMod.setAction
-            (
-                polyModifyFace
+                label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+
+                meshMod.setAction
                 (
-                    pp[patchFaceI],                 // modified face
-                    meshFaceI,                      // label of face
-                    mesh_.faceOwner()[meshFaceI],   // owner
-                    addedCells[patchFaceI][0],      // neighbour
-                    false,                          // face flip
-                    -1,                             // patch for face
-                    false,                          // remove from zone
-                    zoneI,                          // zone for face
-                    false                           // face flip in zone
-                )
-            );
-            //Pout<< "Modified old patch face " << meshFaceI
-            //    << " own:" << mesh_.faceOwner()[meshFaceI]
-            //    << " nei:" << addedCells[patchFaceI][0]
-            //    << endl;
+                    polyModifyFace
+                    (
+                        pp[patchFaceI],                 // modified face
+                        meshFaceI,                      // label of face
+                        mesh_.faceOwner()[meshFaceI],   // owner
+                        addedCells[patchFaceI][0],      // neighbour
+                        false,                          // face flip
+                        -1,                             // patch for face
+                        false,                          // remove from zone
+                        zoneI,                          // zone for face
+                        false                           // face flip in zone
+                    )
+                );
+            }
         }
     }
+    else
+    {
+        // If creating new mesh: reverse original faces and put them
+        // in the exposed patch ID.
+        forAll(pp, patchFaceI)
+        {
+            if (nFaceLayers[patchFaceI] > 0)
+            {
+                label meshFaceI = pp.addressing()[patchFaceI];
+                label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+                bool zoneFlip = false;
+                if (zoneI != -1)
+                {
+                    const faceZone& fz = mesh_.faceZones()[zoneI];
+                    zoneFlip = !fz.flipMap()[fz.whichFace(meshFaceI)];
+                }
+
+                // Reverse and renumber old patch face.
+                face f(pp.localFaces()[patchFaceI].reverseFace());
+                forAll(f, fp)
+                {
+                    f[fp] = copiedPatchPoints[f[fp]];
+                }
+
+                layerFaces_[patchFaceI][0] = meshMod.setAction
+                (
+                    polyAddFace
+                    (
+                        f,                          // modified face
+                        addedCells[patchFaceI][0],  // owner
+                        -1,                         // neighbour
+                        -1,                         // masterPoint
+                        -1,                         // masterEdge
+                        -1,                         // masterFace
+                        true,                       // face flip
+                        exposedPatchID[patchFaceI], // patch for face
+                        zoneI,                      // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+        }
+    }
+
 
 
     //
@@ -1255,7 +1329,16 @@ void Foam::addPatchCellLayer::setRefinement
                         forAll(stringedVerts, stringedI)
                         {
                             label v = stringedVerts[stringedI];
-                            addVertex(meshPoints[v], newFace, newFp);
+                            addVertex
+                            (
+                                (
+                                    addToMesh_
+                                  ? meshPoints[v]
+                                  : copiedPatchPoints[v]
+                                ),
+                                newFace,
+                                newFp
+                            );
                         }
                     }
                     else
@@ -1276,7 +1359,16 @@ void Foam::addPatchCellLayer::setRefinement
                             }
                             else
                             {
-                                addVertex(meshPoints[v], newFace, newFp);
+                                addVertex
+                                (
+                                    (
+                                        addToMesh_
+                                      ? meshPoints[v]
+                                      : copiedPatchPoints[v]
+                                    ),
+                                    newFace,
+                                    newFp
+                                );
                             }
                         }
                     }
@@ -1316,7 +1408,16 @@ void Foam::addPatchCellLayer::setRefinement
                         }
                         else
                         {
-                            addVertex(meshPoints[v], newFace, newFp);
+                            addVertex
+                            (
+                                (
+                                    addToMesh_
+                                  ? meshPoints[v]
+                                  : copiedPatchPoints[v]
+                                ),
+                                newFace,
+                                newFp
+                            );
                         }
                     }
 

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -26,7 +26,8 @@ Class
     Foam::addPatchCellLayer
 
 Description
-    Adds layers of cells to outside of polyPatch.
+    Adds layers of cells to outside of polyPatch. Can optionally create
+    stand-alone extruded mesh (addToMesh=false).
 
     Call setRefinement with offset vector for every patch point and number
     of layers per patch face and number of layers per patch point.
@@ -164,6 +165,9 @@ class addPatchCellLayer
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Add layers to existing mesh or create new mesh
+        const bool addToMesh_;
+
         //- For all patchpoints: list of added points (size 0 or nLayers)
         //  First point in list is one nearest to original point in patch,
         //  last one is the new point on the surface.
@@ -253,7 +257,7 @@ public:
     // Constructors
 
         //- Construct from mesh.
-        addPatchCellLayer(const polyMesh& mesh);
+        addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
 
 
     // Member Functions
@@ -291,8 +295,10 @@ public:
             //- Play commands into polyTopoChange to create layers on top
             //  of indirectPrimitivePatch (have to be outside faces).
             //  Gets displacement per patch point.
-            //  - nPointLayers : number of layers per (patch)point
-            //  - nFaceLayers : number of layers per (patch) face
+            //  - exposedPatchID : only used if creating a new mesh (addToMesh=false)
+            //    gives per pp face the patch the exposed face should get.
+            //  - nPointLayers : number of layers per (patch)point.
+            //  - nFaceLayers : number of layers per (patch) face.
             //  - firstDisplacement : displacement per point for first
             //    layer of points (i.e. nearest to original mesh). If zero
             //    do not add point.
@@ -305,14 +311,15 @@ public:
             //        get a cell should firstDisplacement be <> 0
             //  Note: cells get added from owner cells of patch faces
             //        (instead of e.g. from patch faces)
-           void setRefinement
-           (
-               const scalarField& expansionRatio,
-               const indirectPrimitivePatch& pp,
-               const labelList& nFaceLayers,
-               const labelList& nPointLayers,
-               const vectorField& firstLayerDisp,
-               polyTopoChange& meshMod
+            void setRefinement
+            (
+                const scalarField& expansionRatio,
+                const indirectPrimitivePatch& pp,
+                const labelList& exposedPatchID,
+                const labelList& nFaceLayers,
+                const labelList& nPointLayers,
+                const vectorField& firstLayerDisp,
+                polyTopoChange& meshMod
             );
 
 
@@ -329,6 +336,7 @@ public:
                 (
                     scalarField(pp.nPoints(), 1.0), // expansion ration
                     pp,
+                    labelList(0),
                     labelList(pp.size(), nLayers),
                     labelList(pp.nPoints(), nLayers),
                     overallDisplacement / nLayers,

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -749,6 +749,11 @@ void Foam::polyTopoChange::getFaceOrder
                 " const"
             )   << "Did not determine new position"
                 << " for face " << faceI
+                << " owner " << faceOwner_[faceI]
+                << " neighbour " << faceNeighbour_[faceI]
+                << " region " << region_[faceI] << endl
+                << "This is usually caused by not specifying a patch for"
+                << " a boundary face."
                 << abort(FatalError);
         }
     }
@@ -3176,7 +3181,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
 (
     autoPtr<fvMesh>& newMeshPtr,
     const IOobject& io,
-    const fvMesh& mesh,
+    const polyMesh& mesh,
     const bool syncParallel,
     const bool orderCells,
     const bool orderPoints

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

@@ -585,7 +585,7 @@ public:
             (
                 autoPtr<fvMesh>& newMesh,
                 const IOobject& io,
-                const fvMesh& mesh,
+                const polyMesh& mesh,
                 const bool syncParallel = true,
                 const bool orderCells = false,
                 const bool orderPoints = false

src/dynamicMesh/slidingInterface/coupleSlidingInterface.C

@@ -43,7 +43,7 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::slidingInterface::edgeCoPlanarTol_ = 0.8;
+const Foam::scalar Foam::slidingInterface::edgeCoPlanarTolDefault_ = 0.8;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/dynamicMesh/slidingInterface/slidingInterface.C

@@ -27,10 +27,8 @@ License
 #include "slidingInterface.H"
 #include "polyTopoChanger.H"
 #include "polyMesh.H"
-#include "primitiveMesh.H"
 #include "polyTopoChange.H"
 #include "addToRunTimeSelectionTable.H"
-#include "triPointRef.H"
 #include "plane.H"
 
 // Index of debug signs:
@@ -173,6 +171,14 @@ Foam::slidingInterface::slidingInterface
     attached_(false),
     projectionAlgo_(algo),
     trigger_(false),
+    pointMergeTol_(pointMergeTolDefault_),
+    edgeMergeTol_(edgeMergeTolDefault_),
+    nFacesPerSlaveEdge_(nFacesPerSlaveEdgeDefault_),
+    edgeFaceEscapeLimit_(edgeFaceEscapeLimitDefault_),
+    integralAdjTol_(integralAdjTolDefault_),
+    edgeMasterCatchFraction_(edgeMasterCatchFractionDefault_),
+    edgeCoPlanarTol_(edgeCoPlanarTolDefault_),
+    edgeEndCutoffTol_(edgeEndCutoffTolDefault_),
     cutFaceMasterPtr_(NULL),
     cutFaceSlavePtr_(NULL),
     masterFaceCellsPtr_(NULL),
@@ -280,6 +286,9 @@ Foam::slidingInterface::slidingInterface
     masterPointEdgeHitsPtr_(NULL),
     projectedSlavePointsPtr_(NULL)
 {
+    // Optionally default tolerances from dictionary
+    setTolerances(dict);
+
     checkDefinition();
 
     // If the interface is attached, the master and slave face zone addressing
@@ -686,6 +695,63 @@ const Foam::pointField& Foam::slidingInterface::pointProjection() const
     return *projectedSlavePointsPtr_;
 }
 
+void Foam::slidingInterface::setTolerances(const dictionary&dict, bool report)
+{
+    pointMergeTol_ = dict.lookupOrDefault<scalar>
+    (
+        "pointMergeTol",
+        pointMergeTol_
+    );
+    edgeMergeTol_ = dict.lookupOrDefault<scalar>
+    (
+        "edgeMergeTol",
+        edgeMergeTol_
+    );
+    nFacesPerSlaveEdge_ = dict.lookupOrDefault<scalar>
+    (
+        "nFacesPerSlaveEdge",
+        nFacesPerSlaveEdge_
+    );
+    edgeFaceEscapeLimit_ = dict.lookupOrDefault<scalar>
+    (
+        "edgeFaceEscapeLimit",
+        edgeFaceEscapeLimit_
+    );
+    integralAdjTol_ = dict.lookupOrDefault<scalar>
+    (
+        "integralAdjTol",
+        integralAdjTol_
+    );
+    edgeMasterCatchFraction_ = dict.lookupOrDefault<scalar>
+    (
+        "edgeMasterCatchFraction",
+        edgeMasterCatchFraction_
+    );
+    edgeCoPlanarTol_ = dict.lookupOrDefault<scalar>
+    (
+        "edgeCoPlanarTol",
+        edgeCoPlanarTol_
+    );
+    edgeEndCutoffTol_ = dict.lookupOrDefault<scalar>
+    (
+        "edgeEndCutoffTol",
+        edgeEndCutoffTol_
+    );
+
+    if (report)
+    {
+        Info<< "Sliding interface parameters:" << nl
+            << "pointMergeTol            : " << pointMergeTol_ << nl
+            << "edgeMergeTol             : " << edgeMergeTol_ << nl
+            << "nFacesPerSlaveEdge       : " << nFacesPerSlaveEdge_ << nl
+            << "edgeFaceEscapeLimit      : " << edgeFaceEscapeLimit_ << nl
+            << "integralAdjTol           : " << integralAdjTol_ << nl
+            << "edgeMasterCatchFraction  : " << edgeMasterCatchFraction_ << nl
+            << "edgeCoPlanarTol          : " << edgeCoPlanarTol_ << nl
+            << "edgeEndCutoffTol         : " << edgeEndCutoffTol_ << endl;
+    }
+}
+
 
 void Foam::slidingInterface::write(Ostream& os) const
 {
@@ -703,6 +769,14 @@ void Foam::slidingInterface::write(Ostream& os) const
 }
 
 
+// To write out all those tolerances
+#define WRITE_NON_DEFAULT(name) \
+    if( name ## _ != name ## Default_ )\
+    { \
+        os << "    " #name " " <<  name ## _ << token::END_STATEMENT << nl; \
+    }
+
+
 void Foam::slidingInterface::writeDict(Ostream& os) const
 {
     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
@@ -743,6 +817,15 @@ void Foam::slidingInterface::writeDict(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
+    WRITE_NON_DEFAULT(pointMergeTol)
+    WRITE_NON_DEFAULT(edgeMergeTol)
+    WRITE_NON_DEFAULT(nFacesPerSlaveEdge)
+    WRITE_NON_DEFAULT(edgeFaceEscapeLimit)
+    WRITE_NON_DEFAULT(integralAdjTol)
+    WRITE_NON_DEFAULT(edgeMasterCatchFraction)
+    WRITE_NON_DEFAULT(edgeCoPlanarTol)
+    WRITE_NON_DEFAULT(edgeEndCutoffTol)
+
     os  << token::END_BLOCK << endl;
 }
 

src/dynamicMesh/slidingInterface/slidingInterface.H

@@ -129,6 +129,32 @@ private:
         //- Trigger topological change
         mutable bool trigger_;
 
+        // Tolerances. Initialised to static ones below.
+
+            //- Point merge tolerance
+            scalar pointMergeTol_;
+
+            //- Edge merge tolerance
+            scalar edgeMergeTol_;
+
+            //- Estimated number of faces an edge goes through
+            label nFacesPerSlaveEdge_;
+
+            //- Edge-face interaction escape limit
+            label edgeFaceEscapeLimit_;
+
+            //- Integral match point adjustment tolerance
+            scalar integralAdjTol_;
+
+            //- Edge intersection master catch fraction
+            scalar edgeMasterCatchFraction_;
+
+            //- Edge intersection co-planar tolerance
+            scalar edgeCoPlanarTol_;
+
+            //- Edge end cut-off tolerance
+            scalar edgeEndCutoffTol_;
+
 
         // Private addressing data.
 
@@ -256,28 +282,28 @@ private:
     // Static data members
 
         //- Point merge tolerance
-        static const scalar pointMergeTol_;
+        static const scalar pointMergeTolDefault_;
 
         //- Edge merge tolerance
-        static const scalar edgeMergeTol_;
+        static const scalar edgeMergeTolDefault_;
 
         //- Estimated number of faces an edge goes through
-        static const label nFacesPerSlaveEdge_;
+        static const label nFacesPerSlaveEdgeDefault_;
 
         //- Edge-face interaction escape limit
-        static const label edgeFaceEscapeLimit_;
+        static const label edgeFaceEscapeLimitDefault_;
 
         //- Integral match point adjustment tolerance
-        static const scalar integralAdjTol_;
+        static const scalar integralAdjTolDefault_;
 
         //- Edge intersection master catch fraction
-        static const scalar edgeMasterCatchFraction_;
+        static const scalar edgeMasterCatchFractionDefault_;
 
         //- Edge intersection co-planar tolerance
-        static const scalar edgeCoPlanarTol_;
+        static const scalar edgeCoPlanarTolDefault_;
 
         //- Edge end cut-off tolerance
-        static const scalar edgeEndCutoffTol_;
+        static const scalar edgeEndCutoffTolDefault_;
 
 
 public:
@@ -350,6 +376,8 @@ public:
         //- Return projected points for a slave patch
         const pointField& pointProjection() const;
 
+        //- Set the tolerances from the values in a dictionary
+        void setTolerances(const dictionary&, bool report=false);
 
         //- Write
         virtual void write(Ostream&) const;

src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C

@@ -26,22 +26,19 @@ License
 
 #include "slidingInterface.H"
 #include "polyMesh.H"
-#include "primitiveMesh.H"
 #include "line.H"
-#include "triPointRef.H"
-#include "plane.H"
 #include "polyTopoChanger.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::scalar Foam::slidingInterface::pointMergeTol_ = 0.05;
-const Foam::scalar Foam::slidingInterface::edgeMergeTol_ = 0.01;
-const Foam::label Foam::slidingInterface::nFacesPerSlaveEdge_ = 5;
-const Foam::label Foam::slidingInterface::edgeFaceEscapeLimit_ = 10;
+const Foam::scalar Foam::slidingInterface::pointMergeTolDefault_ = 0.05;
+const Foam::scalar Foam::slidingInterface::edgeMergeTolDefault_ = 0.01;
+const Foam::label Foam::slidingInterface::nFacesPerSlaveEdgeDefault_ = 5;
+const Foam::label Foam::slidingInterface::edgeFaceEscapeLimitDefault_ = 10;
 
-const Foam::scalar Foam::slidingInterface::integralAdjTol_ = 0.05;
-const Foam::scalar Foam::slidingInterface::edgeMasterCatchFraction_ = 0.4;
-const Foam::scalar Foam::slidingInterface::edgeEndCutoffTol_ = 0.0001;
+const Foam::scalar Foam::slidingInterface::integralAdjTolDefault_ = 0.05;
+const Foam::scalar Foam::slidingInterface::edgeMasterCatchFractionDefault_ = 0.4;
+const Foam::scalar Foam::slidingInterface::edgeEndCutoffTolDefault_ = 0.0001;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C

@@ -45,7 +45,7 @@ Foam::linearUpwind<Type>::correction
         (
             IOobject
             (
-                "linearUpwindCorrection(" + vf.name() + ')',
+                "linearUpwind::correction(" + vf.name() + ')',
                 mesh.time().timeName(),
                 mesh,
                 IOobject::NO_READ,

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C

@@ -38,93 +38,96 @@ Foam::linearUpwindV<Type>::correction
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
-   const fvMesh& mesh = this->mesh();
+    const fvMesh& mesh = this->mesh();
 
-   tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
-   (
-       new GeometricField<Type, fvsPatchField, surfaceMesh>
-       (
-           IOobject
-           (
-               vf.name(),
-               mesh.time().timeName(),
-               mesh
-           ),
-           mesh,
-           dimensioned<Type>
-           (
-               vf.name(),
-               vf.dimensions(),
-               pTraits<Type>::zero
-           )
-       )
-   );
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "linearUpwindV::correction(" + vf.name() + ')',
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                vf.name(),
+                vf.dimensions(),
+                pTraits<Type>::zero
+            )
+        )
+    );
 
-   GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr();
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr();
 
-   const surfaceScalarField& faceFlux = this->faceFlux_;
-   const surfaceScalarField& w = mesh.weights();
+    const surfaceScalarField& faceFlux = this->faceFlux_;
+    const surfaceScalarField& w = mesh.weights();
 
-   const labelList& own = mesh.owner();
-   const labelList& nei = mesh.neighbour();
+    const labelList& own = mesh.owner();
+    const labelList& nei = mesh.neighbour();
 
-   const vectorField& C = mesh.C();
-   const vectorField& Cf = mesh.Cf();
+    const vectorField& C = mesh.C();
+    const vectorField& Cf = mesh.Cf();
 
-   GeometricField
-       <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-       gradVf = gradScheme_().grad(vf);
+    GeometricField
+        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        gradVf = gradScheme_().grad(vf);
 
-   forAll(faceFlux, facei)
-   {
-       vector maxCorr;
+    forAll(faceFlux, facei)
+    {
+        vector maxCorr;
 
-       if (faceFlux[facei] > 0.0)
-       {
-           maxCorr =
-               (1.0 - w[facei])
+        if (faceFlux[facei] > 0.0)
+        {
+            maxCorr =
+                (1.0 - w[facei])
                *(vf[nei[facei]] - vf[own[facei]]);
 
-           sfCorr[facei] =
-               (Cf[facei] - C[own[facei]]) & gradVf[own[facei]];
-       }
-       else
-       {
-           maxCorr =
-               w[facei]*(vf[own[facei]] - vf[nei[facei]]);
+            sfCorr[facei] =
+                (Cf[facei] - C[own[facei]]) & gradVf[own[facei]];
+        }
+        else
+        {
+            maxCorr =
+                w[facei]*(vf[own[facei]] - vf[nei[facei]]);
 
-           sfCorr[facei] =
+            sfCorr[facei] =
                (Cf[facei] - C[nei[facei]]) & gradVf[nei[facei]];
-       }
+        }
 
-       scalar sfCorrs = magSqr(sfCorr[facei]);
-       scalar maxCorrs = sfCorr[facei] & maxCorr;
+        scalar sfCorrs = magSqr(sfCorr[facei]);
+        scalar maxCorrs = sfCorr[facei] & maxCorr;
 
-       if (sfCorrs > 0)
-       {
-           if (maxCorrs < 0)
-           {
-               sfCorr[facei] = vector::zero;
-           }
-           else if (sfCorrs > maxCorrs)
-           {
-               sfCorr[facei] *= maxCorrs/(sfCorrs + VSMALL);
-           }
-       }
-       else if (sfCorrs < 0)
-       {
-           if (maxCorrs > 0)
-           {
-               sfCorr[facei] = vector::zero;
-           }
-           else if (sfCorrs < maxCorrs)
-           {
-               sfCorr[facei] *= maxCorrs/(sfCorrs - VSMALL);
-           }
-       }
-   }
+        if (sfCorrs > 0)
+        {
+            if (maxCorrs < 0)
+            {
+                sfCorr[facei] = vector::zero;
+            }
+            else if (sfCorrs > maxCorrs)
+            {
+                sfCorr[facei] *= maxCorrs/(sfCorrs + VSMALL);
+            }
+        }
+        else if (sfCorrs < 0)
+        {
+            if (maxCorrs > 0)
+            {
+                sfCorr[facei] = vector::zero;
+            }
+            else if (sfCorrs < maxCorrs)
+            {
+                sfCorr[facei] *= maxCorrs/(sfCorrs - VSMALL);
+            }
+        }
+    }
 
-   return tsfCorr;
+    return tsfCorr;
 }
 
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/cubic/cubic.H

@@ -131,9 +131,12 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "cubic::correction(" + vf.name() +')',
                         mesh.time().timeName(),
-                        mesh
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE,
+                        false
                     ),
                     surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc)
                 )

src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMax/localMax.H

@@ -133,7 +133,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "localMax::interpolate(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

src/finiteVolume/interpolation/surfaceInterpolation/schemes/localMin/localMin.H

@@ -133,7 +133,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "localMin::interpolate(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrected.H

@@ -144,7 +144,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "skewCorrected::skewCorrection(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

src/genericPatchFields/genericFvPatchField/genericFvPatchField.C

@@ -79,7 +79,6 @@ Foam::genericFvPatchField<Type>::genericFvPatchField
             << "    (Actual type " << actualTypeName_ << ")" << nl
             << "\n    Please add the 'value' entry to the write function "
                "of the user-defined boundary-condition\n"
-               "    or link the boundary-condition into libfoamUtil.so"
             << exit(FatalIOError);
     }
 

src/meshTools/meshSearch/meshSearch.C

@@ -428,7 +428,7 @@ Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
 :
     mesh_(mesh),
     faceDecomp_(faceDecomp),
-    cloud_(mesh_, IDLList<passiveParticle>()),
+    cloud_(mesh_, "meshSearchCloud", IDLList<passiveParticle>()),
     boundaryTreePtr_(NULL),
     cellTreePtr_(NULL),
     cellCentreTreePtr_(NULL)

src/meshTools/meshSearch/meshSearch.H

@@ -38,8 +38,7 @@ SourceFiles
 #define meshSearch_H
 
 #include "pointIndexHit.H"
-#include "Cloud.H"
-#include "passiveParticle.H"
+#include "passiveParticleCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,7 +67,7 @@ class meshSearch
         const bool faceDecomp_;
 
         //- Dummy cloud to put particles on for tracking.
-        Cloud<passiveParticle> cloud_;
+        passiveParticleCloud cloud_;
 
         //- demand driven octrees
 

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

@@ -2560,7 +2560,7 @@ Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
 
     // Triangulate
     int nTris = 0;
-    dtris2
+    int err = dtris2
     (
         pts.size(),
         geompackVertices.begin(),
@@ -2569,6 +2569,13 @@ Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
         triangle_neighbor.begin()
     );
 
+    if (err != 0)
+    {
+        FatalErrorIn("triSurfaceTools::delaunay2D(const List<vector2D>&)")
+            << "Failed dtris2 with vertices:" << pts.size()
+            << abort(FatalError);
+    }
+
     // Trim
     triangle_node.setSize(3*nTris);
     triangle_neighbor.setSize(3*nTris);

src/postProcessing/functionObjects/field/Make/files

@@ -6,4 +6,16 @@ fieldAverage/fieldAverageFunctionObject/fieldAverageFunctionObject.C
 fieldMinMax/fieldMinMax.C
 fieldMinMax/fieldMinMaxFunctionObject.C
 
+fieldValues/fieldValue/fieldValue.C
+fieldValues/face/faceSource.C
+fieldValues/face/faceSourceFunctionObject.C
+fieldValues/cell/cellSource.C
+fieldValues/cell/cellSourceFunctionObject.C
+
+streamLine/streamLine.C
+streamLine/streamLineParticle.C
+streamLine/streamLineParticleCloud.C
+streamLine/streamLineFunctionObject.C
+
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

src/postProcessing/functionObjects/field/Make/options

@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -llagrangian \
     -lsampling

src/postProcessing/functionObjects/field/fieldValues/cell/IOcellSource.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+Typedef
+    Foam::IOcellSource
+
+
+Description
+    Instance of the generic IOOutputFilter for cellSource.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOcellSource_H
+#define IOcellSource_H
+
+#include "cellSource.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOOutputFilter<cellSource> IOcellSource;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.C

@@ -0,0 +1,257 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "cellSource.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "IOList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace fieldValues
+    {
+        defineTypeNameAndDebug(cellSource, 0);
+    }
+
+    template<>
+    const char* NamedEnum<fieldValues::cellSource::sourceType, 1>::
+        names[] = {"cellZone"};
+
+    const NamedEnum<fieldValues::cellSource::sourceType, 1>
+        fieldValues::cellSource::sourceTypeNames_;
+
+    template<>
+    const char* NamedEnum<fieldValues::cellSource::operationType, 4>::
+        names[] = {"none", "sum", "volAverage", "volIntegrate"};
+
+    const NamedEnum<fieldValues::cellSource::operationType, 4>
+        fieldValues::cellSource::operationTypeNames_;
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fieldValues::cellSource::setCellZoneCells()
+{
+    label zoneId = mesh().cellZones().findZoneID(sourceName_);
+
+    if (zoneId < 0)
+    {
+        FatalErrorIn("cellSource::cellSource::setCellZoneCells()")
+            << "Unknown cell zone name: " << sourceName_
+            << ". Valid cell zones are: " << mesh().cellZones().names()
+            << nl << exit(FatalError);
+    }
+
+    const cellZone& cZone = mesh().cellZones()[zoneId];
+
+    cellId_.setSize(cZone.size());
+
+    label count = 0;
+    forAll(cZone, i)
+    {
+        label cellI = cZone[i];
+        cellId_[count] = cellI;
+        count++;
+    }
+
+    cellId_.setSize(count);
+
+    if (debug)
+    {
+        Info<< "Original cell zone size = " << cZone.size()
+            << ", new size = " << count << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::fieldValues::cellSource::initialise()
+{
+    switch (source_)
+    {
+        case stCellZone:
+        {
+            setCellZoneCells();
+            break;
+        }
+        default:
+        {
+            FatalErrorIn("cellSource::constructCellAddressing()")
+                << "Unknown source type. Valid source types are:"
+                << sourceTypeNames_ << nl << exit(FatalError);
+        }
+    }
+
+    Info<< type() << " " << name_ << ":" << nl
+        << "    total cells  = " << cellId_.size() << nl
+        << "    total volume = " << sum(filterField(mesh().V()))
+        << nl << endl;
+}
+
+
+void Foam::fieldValues::cellSource::makeFile()
+{
+    // Create the forces file if not already created
+    if (outputFilePtr_.empty())
+    {
+        if (debug)
+        {
+            Info<< "Creating output file." << endl;
+        }
+
+        // File update
+        if (Pstream::master())
+        {
+            fileName outputDir;
+            if (Pstream::parRun())
+            {
+                // Put in undecomposed case (Note: gives problems for
+                // distributed data running)
+                outputDir =
+                    obr_.time().path()/".."/name_/obr_.time().timeName();
+            }
+            else
+            {
+                outputDir = obr_.time().path()/name_/obr_.time().timeName();
+            }
+
+            // Create directory if does not exist
+            mkDir(outputDir);
+
+            // Open new file at start up
+            outputFilePtr_.reset(new OFstream(outputDir/(type() + ".dat")));
+
+            // Add headers to output data
+            writeFileHeader();
+        }
+    }
+}
+
+
+void Foam::fieldValues::cellSource::writeFileHeader()
+{
+    if (outputFilePtr_.valid())
+    {
+        outputFilePtr_()
+            << "# Source : " << sourceTypeNames_[source_] << " "
+            << sourceName_ <<  nl << "# Cells  : " << cellId_.size() << nl
+            << "# Time" << tab << "sum(V)";
+
+        forAll(fields_, i)
+        {
+            outputFilePtr_()
+                << tab << operationTypeNames_[operation_]
+                << "(" << fields_[i] << ")";
+        }
+
+        outputFilePtr_() << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fieldValues::cellSource::cellSource
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    fieldValue(name, obr, dict, loadFromFiles),
+    source_(sourceTypeNames_.read(dict.lookup("source"))),
+    operation_(operationTypeNames_.read(dict.lookup("operation"))),
+    cellId_(),
+    outputFilePtr_(NULL)
+{
+    initialise();
+
+    if (active_)
+    {
+        // Create the output file if not already created
+        makeFile();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fieldValues::cellSource::~cellSource()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fieldValues::cellSource::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        fieldValue::read(dict);
+        initialise();
+    }
+}
+
+
+void Foam::fieldValues::cellSource::write()
+{
+    if (active_)
+    {
+        if (log_)
+        {
+            Info<< type() << " " << name_ << " output:" << nl;
+        }
+
+        outputFilePtr_()
+            << obr_.time().value() << tab
+            << sum(filterField(mesh().V()));
+
+        forAll(fields_, i)
+        {
+            writeValues<scalar>(fields_[i]);
+            writeValues<vector>(fields_[i]);
+            writeValues<sphericalTensor>(fields_[i]);
+            writeValues<symmTensor>(fields_[i]);
+            writeValues<tensor>(fields_[i]);
+        }
+
+        outputFilePtr_()<< endl;
+
+        if (log_)
+        {
+            Info<< endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.H

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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::fieldValues::cellSource
+
+Description
+    Cell source variant of field value function object. Values of user-
+    specified fields reported for collections of cells.
+
+    cellObj1                        // Name also used to identify output folder
+    {
+        type            cellSource;
+        functionObjectLibs ("libfieldValueFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;       // log to screen?
+        valueOutput     true;       // Write values at run-time output times?
+        source          cellZone;   // Type of cell source
+        sourceName      c0;
+        operation       volAverage; // none, sum, volAverage, volIntegrate
+        fields
+        (
+            p
+            U
+        );
+    }
+
+SourceFiles
+    cellSource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellSource_H
+#define cellSource_H
+
+#include "NamedEnum.H"
+#include "fieldValue.H"
+#include "labelList.H"
+#include "OFstream.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fieldValues
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class cellSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellSource
+:
+    public fieldValue
+{
+
+public:
+
+    // Public data types
+
+        //- Source type enumeration
+        enum sourceType
+        {
+            stCellZone
+        };
+
+        //- Source type names
+        static const NamedEnum<sourceType, 1> sourceTypeNames_;
+
+
+        //- Operation type enumeration
+        enum operationType
+        {
+            opNone,
+            opSum,
+            opVolAverage,
+            opVolIntegrate
+        };
+
+        //- Operation type names
+        static const NamedEnum<operationType, 4> operationTypeNames_;
+
+
+private:
+
+    // Private member functions
+
+        //- Set cells to evaluate based on a cell zone
+        void setCellZoneCells();
+
+        //- Set cells to evaluate based on a patch
+        void setPatchCells();
+
+        //- Create the output file if not already created
+        void makeFile();
+
+
+protected:
+
+    // Protected data
+
+        //- Source type
+        sourceType source_;
+
+        //- Operation to apply to values
+        operationType operation_;
+
+        //- Local list of cell IDs
+        labelList cellId_;
+
+        //- Output file pointer
+        autoPtr<OFstream> outputFilePtr_;
+
+
+    // Protected member functions
+
+        //- Initialise, e.g. cell addressing
+        void initialise();
+
+        //- Insert field values into values list
+        template<class Type>
+        bool setFieldValues
+        (
+            const word& fieldName,
+            List<Type>& values
+        ) const;
+
+        //- Apply the 'operation' to the values
+        template<class Type>
+        Type processValues(const List<Type>& values) const;
+
+        //- Output file header information
+        virtual void writeFileHeader();
+
+
+public:
+
+    //- Run-time type information
+    TypeName("cellSource");
+
+
+    //- Construct from components
+    cellSource
+    (
+        const word& name,
+        const objectRegistry& obr,
+        const dictionary& dict,
+        const bool loadFromFiles = false
+    );
+
+
+    //- Destructor
+    virtual ~cellSource();
+
+
+    // Public member functions
+
+        // Access
+
+            //- Return the source type
+            inline const sourceType& source() const;
+
+            //- Return the local list of cell IDs
+            inline const labelList& cellId() const;
+
+
+        // Function object functions
+
+            //- Read from dictionary
+            virtual void read(const dictionary&);
+
+            //- Calculate and write
+            virtual void write();
+
+            //- Templated helper function to output field values
+            template<class Type>
+            bool writeValues(const word& fieldName);
+
+            //- Filter a field according to cellIds
+            template<class Type>
+            tmp<Field<Type> > filterField(const Field<Type>& field) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fieldValues
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "cellSourceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "cellSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.C

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "cellSourceFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug
+    (
+        cellSourceFunctionObject,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        cellSourceFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+Typedef
+    Foam::cellSourceFunctionObject
+
+Description
+    FunctionObject wrapper around cellSource to allow it to be
+    created via the functions list within controlDict.
+
+SourceFiles
+    cellSourceFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellSourceFunctionObject_H
+#define cellSourceFunctionObject_H
+
+#include "cellSource.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<fieldValues::cellSource>
+        cellSourceFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceI.H

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "cellSource.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::fieldValues::cellSource::sourceType&
+Foam::fieldValues::cellSource::source() const
+{
+    return source_;
+}
+
+
+inline const Foam::labelList&
+Foam::fieldValues::cellSource::cellId() const
+{
+    return cellId_;
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceTemplates.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "cellSource.H"
+#include "volFields.H"
+#include "IOList.H"
+#include "ListListOps.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::fieldValues::cellSource::setFieldValues
+(
+    const word& fieldName,
+    List<Type>& values
+) const
+{
+    values.setSize(cellId_.size(), pTraits<Type>::zero);
+
+    typedef GeometricField<Type, fvPatchField, volMesh> vf;
+
+    if (obr_.foundObject<vf>(fieldName))
+    {
+        const vf& field = obr_.lookupObject<vf>(fieldName);
+
+        values = UIndirectList<Type>(field, cellId_);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class Type>
+Type Foam::fieldValues::cellSource::processValues
+(
+    const List<Type>& values
+) const
+{
+    Type result = pTraits<Type>::zero;
+    switch (operation_)
+    {
+        case opSum:
+        {
+            result = sum(values);
+            break;
+        }
+        case opVolAverage:
+        {
+            tmp<scalarField> V = filterField(mesh().V());
+            result = sum(values*V())/sum(V());
+            break;
+        }
+        case opVolIntegrate:
+        {
+            result = sum(values*filterField(mesh().V()));
+            break;
+        }
+        default:
+        {
+            // Do nothing
+        }
+    }
+
+    return result;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::fieldValues::cellSource::writeValues(const word& fieldName)
+{
+    List<List<Type> > allValues(Pstream::nProcs());
+
+    bool validField =
+        setFieldValues<Type>(fieldName, allValues[Pstream::myProcNo()]);
+
+    if (validField)
+    {
+        Pstream::gatherList(allValues);
+
+        if (Pstream::master())
+        {
+            List<Type> values =
+                ListListOps::combine<List<Type> >
+                (
+                    allValues,
+                    accessOp<List<Type> >()
+                );
+
+            Type result = processValues(values);
+
+            if (valueOutput_)
+            {
+                IOList<Type>
+                (
+                    IOobject
+                    (
+                        fieldName + "_" + sourceTypeNames_[source_] + "-"
+                            + sourceName_,
+                        obr_.time().timeName(),
+                        obr_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    values
+                ).write();
+            }
+
+
+            outputFilePtr_()<< tab << result;
+
+            if (log_)
+            {
+                Info<< "    " << operationTypeNames_[operation_]
+                    << "(" << sourceName_ << ") for " << fieldName
+                    <<  " = " << result << endl;
+            }
+        }
+    }
+
+    return validField;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::fieldValues::cellSource::filterField
+(
+    const Field<Type>& field
+) const
+{
+    return tmp<Field<Type> >(new Field<Type>(field, cellId_));
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/controlDict

@@ -0,0 +1,108 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          0.005;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+functions
+(
+    faceObj1
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;
+        valueOutput     true;
+        source          patch;
+        sourceName      movingWall;
+//        source          faceZone;
+//        sourceName      f0;
+        operation       areaAverage;
+        fields
+        (
+            p
+            phi
+            U
+        );
+    }
+
+    faceObj2
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;
+        valueOutput     true;
+        source          faceZone;
+        sourceName      f0;
+        operation       sum;
+        fields
+        (
+            phi
+        );
+    }
+
+    cellObj1
+    {
+        type            cellSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;
+        valueOutput     true;
+        source          cellZone;
+        sourceName      c0;
+        operation       volAverage;
+        fields
+        (
+            p
+            U
+        );
+    }
+);
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/face/IOfaceSource.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+Typedef
+    Foam::IOfaceSource
+
+
+Description
+    Instance of the generic IOOutputFilter for faceSource.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOfaceSource_H
+#define IOfaceSource_H
+
+#include "faceSource.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOOutputFilter<faceSource> IOfaceSource;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/face/faceSource.C

@@ -0,0 +1,368 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "faceSource.H"
+#include "fvMesh.H"
+#include "cyclicPolyPatch.H"
+#include "emptyPolyPatch.H"
+#include "processorPolyPatch.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "IOList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace fieldValues
+    {
+        defineTypeNameAndDebug(faceSource, 0);
+    }
+
+    template<>
+    const char* NamedEnum<fieldValues::faceSource::sourceType, 2>::
+        names[] = {"faceZone", "patch"};
+
+    const NamedEnum<fieldValues::faceSource::sourceType, 2>
+        fieldValues::faceSource::sourceTypeNames_;
+
+    template<>
+    const char* NamedEnum<fieldValues::faceSource::operationType, 4>::
+        names[] = {"none", "sum", "areaAverage", "areaIntegrate"};
+
+    const NamedEnum<fieldValues::faceSource::operationType, 4>
+        fieldValues::faceSource::operationTypeNames_;
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fieldValues::faceSource::setFaceZoneFaces()
+{
+    label zoneId = mesh().faceZones().findZoneID(sourceName_);
+
+    if (zoneId < 0)
+    {
+        FatalErrorIn("faceSource::faceSource::setFaceZoneFaces()")
+            << "Unknown face zone name: " << sourceName_
+            << ". Valid face zones are: " << mesh().faceZones().names()
+            << nl << exit(FatalError);
+    }
+
+    const faceZone& fZone = mesh().faceZones()[zoneId];
+
+    faceId_.setSize(fZone.size());
+    facePatchId_.setSize(fZone.size());
+    flipMap_.setSize(fZone.size());
+
+    label count = 0;
+    forAll(fZone, i)
+    {
+        label faceI = fZone[i];
+
+        label faceId = -1;
+        label facePatchId = -1;
+        if (mesh().isInternalFace(faceI))
+        {
+            faceId = faceI;
+            facePatchId = -1;
+        }
+        else
+        {
+            facePatchId = mesh().boundaryMesh().whichPatch(faceI);
+            const polyPatch& pp = mesh().boundaryMesh()[facePatchId];
+            if (isA<processorPolyPatch>(pp))
+            {
+                if (refCast<const processorPolyPatch>(pp).owner())
+                {
+                    faceId = pp.whichFace(faceI);
+                }
+                else
+                {
+                    faceId = -1;
+                }
+            }
+            else if (isA<cyclicPolyPatch>(pp))
+            {
+                label patchFaceI = faceI - pp.start();
+                if (patchFaceI < pp.size()/2)
+                {
+                    faceId = patchFaceI;
+                }
+                else
+                {
+                    faceId = -1;
+                }
+            }
+            else if (!isA<emptyPolyPatch>(pp))
+            {
+                faceId = faceI - pp.start();
+            }
+            else
+            {
+                faceId = -1;
+                facePatchId = -1;
+            }
+        }
+
+        if (faceId >= 0)
+        {
+            if (fZone.flipMap()[i])
+            {
+                flipMap_[count] = -1;
+            }
+            else
+            {
+                flipMap_[count] = 1;
+            }
+            faceId_[count] = faceId;
+            facePatchId_[count] = facePatchId;
+            count++;
+        }
+    }
+
+    faceId_.setSize(count);
+    facePatchId_.setSize(count);
+    flipMap_.setSize(count);
+
+    if (debug)
+    {
+        Info<< "Original face zone size = " << fZone.size()
+            << ", new size = " << count << endl;
+    }
+}
+
+
+void Foam::fieldValues::faceSource::setPatchFaces()
+{
+    label patchId = mesh().boundaryMesh().findPatchID(sourceName_);
+
+    if (patchId < 0)
+    {
+        FatalErrorIn("faceSource::constructFaceAddressing()")
+            << "Unknown patch name: " << sourceName_
+            << ". Valid patch names are: "
+            << mesh().boundaryMesh().names() << nl
+            << exit(FatalError);
+    }
+
+    const polyPatch& pp = mesh().boundaryMesh()[patchId];
+
+    label nFaces = pp.size();
+    if (isA<cyclicPolyPatch>(pp))
+    {
+        nFaces /= 2;
+    }
+    else if (isA<emptyPolyPatch>(pp))
+    {
+        nFaces = 0;
+    }
+
+    faceId_.setSize(nFaces);
+    facePatchId_.setSize(nFaces);
+    flipMap_.setSize(nFaces);
+
+    forAll(faceId_, faceI)
+    {
+        faceId_[faceI] = faceI;
+        facePatchId_[faceI] = patchId;
+        flipMap_[faceI] = 1;
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::fieldValues::faceSource::initialise()
+{
+    switch (source_)
+    {
+        case stFaceZone:
+        {
+            setFaceZoneFaces();
+            break;
+        }
+        case stPatch:
+        {
+            setPatchFaces();
+            break;
+        }
+        default:
+        {
+            FatalErrorIn("faceSource::constructFaceAddressing()")
+                << "Unknown source type. Valid source types are:"
+                << sourceTypeNames_ << nl << exit(FatalError);
+        }
+    }
+
+    Info<< type() << " " << name_ << ":" << nl
+        << "    total faces = " << faceId_.size() << nl
+        << "    total area  = " << sum(filterField(mesh().magSf()))
+        << nl << endl;
+}
+
+
+void Foam::fieldValues::faceSource::makeFile()
+{
+    // Create the forces file if not already created
+    if (outputFilePtr_.empty())
+    {
+        if (debug)
+        {
+            Info<< "Creating output file." << endl;
+        }
+
+        // File update
+        if (Pstream::master())
+        {
+            fileName outputDir;
+            if (Pstream::parRun())
+            {
+                // Put in undecomposed case (Note: gives problems for
+                // distributed data running)
+                outputDir =
+                    obr_.time().path()/".."/name_/obr_.time().timeName();
+            }
+            else
+            {
+                outputDir = obr_.time().path()/name_/obr_.time().timeName();
+            }
+
+            // Create directory if does not exist
+            mkDir(outputDir);
+
+            // Open new file at start up
+            outputFilePtr_.reset(new OFstream(outputDir/(type() + ".dat")));
+
+            // Add headers to output data
+            writeFileHeader();
+        }
+    }
+}
+
+
+void Foam::fieldValues::faceSource::writeFileHeader()
+{
+    if (outputFilePtr_.valid())
+    {
+        outputFilePtr_()
+            << "# Source : " << sourceTypeNames_[source_] << " "
+            << sourceName_ <<  nl << "# Faces  : " << faceId_.size() << nl
+            << "# Time" << tab << "sum(magSf)";
+
+        forAll(fields_, i)
+        {
+            outputFilePtr_()
+                << tab << operationTypeNames_[operation_]
+                << "(" << fields_[i] << ")";
+        }
+
+        outputFilePtr_() << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fieldValues::faceSource::faceSource
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    fieldValue(name, obr, dict, loadFromFiles),
+    source_(sourceTypeNames_.read(dict.lookup("source"))),
+    operation_(operationTypeNames_.read(dict.lookup("operation"))),
+    faceId_(),
+    facePatchId_(),
+    flipMap_(),
+    outputFilePtr_(NULL)
+{
+    initialise();
+
+    if (active_)
+    {
+        // Create the output file if not already created
+        makeFile();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fieldValues::faceSource::~faceSource()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fieldValues::faceSource::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        fieldValue::read(dict);
+        initialise();
+    }
+}
+
+
+void Foam::fieldValues::faceSource::write()
+{
+    if (active_)
+    {
+        if (log_)
+        {
+            Info<< type() << " " << name_ << " output:" << nl;
+        }
+
+        outputFilePtr_()
+            << obr_.time().value() << tab
+            << sum(filterField(mesh().magSf()));
+
+        forAll(fields_, i)
+        {
+            writeValues<scalar>(fields_[i]);
+            writeValues<vector>(fields_[i]);
+            writeValues<sphericalTensor>(fields_[i]);
+            writeValues<symmTensor>(fields_[i]);
+            writeValues<tensor>(fields_[i]);
+        }
+
+        outputFilePtr_()<< endl;
+
+        if (log_)
+        {
+            Info<< endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/face/faceSource.H

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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::fieldValues::faceSource
+
+Description
+    Face source variant of field value function object. Values of user-
+    specified fields reported for collections of faces.
+
+    cellObj1                        // Name also used to identify output folder
+    {
+        type            cellSource;
+        functionObjectLibs ("libfieldValueFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;       // log to screen?
+        valueOutput     true;       // Write values at run-time output times?
+        source          faceZone;   // Type of face source: faceZone, patch
+        sourceName      f0;
+        operation       sum;        // none, sum, areaAverage, areaIntegrate
+        fields
+        (
+            p
+            phi
+            U
+        );
+    }
+
+SourceFiles
+    faceSource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceSource_H
+#define faceSource_H
+
+#include "NamedEnum.H"
+#include "fieldValue.H"
+#include "labelList.H"
+#include "OFstream.H"
+#include "surfaceFieldsFwd.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fieldValues
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class faceSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceSource
+:
+    public fieldValue
+{
+
+public:
+
+    // Public data types
+
+        //- Source type enumeration
+        enum sourceType
+        {
+            stFaceZone,
+            stPatch
+        };
+
+        //- Source type names
+        static const NamedEnum<sourceType, 2> sourceTypeNames_;
+
+
+        //- Operation type enumeration
+        enum operationType
+        {
+            opNone,
+            opSum,
+            opAreaAverage,
+            opAreaIntegrate
+        };
+
+        //- Operation type names
+        static const NamedEnum<operationType, 4> operationTypeNames_;
+
+
+private:
+
+    // Private member functions
+
+        //- Set faces to evaluate based on a face zone
+        void setFaceZoneFaces();
+
+        //- Set faces to evaluate based on a patch
+        void setPatchFaces();
+
+        //- Create the output file if not already created
+        void makeFile();
+
+
+protected:
+
+    // Protected data
+
+        //- Source type
+        sourceType source_;
+
+        //- Operation to apply to values
+        operationType operation_;
+
+        //- Local list of face IDs
+        labelList faceId_;
+
+        //- Local list of patch ID per face
+        labelList facePatchId_;
+
+        //- List of +1/-1 representing face flip map
+        labelList flipMap_;
+
+        //- Output file pointer
+        autoPtr<OFstream> outputFilePtr_;
+
+
+    // Protected member functions
+
+        //- Initialise, e.g. face addressing
+        void initialise();
+
+        //- Insert field values into values list
+        template<class Type>
+        bool setFieldValues
+        (
+            const word& fieldName,
+            List<Type>& values
+        ) const;
+
+        //- Apply the 'operation' to the values
+        template<class Type>
+        Type processValues(const List<Type>& values) const;
+
+        //- Output file header information
+        virtual void writeFileHeader();
+
+
+public:
+
+    //- Run-time type information
+    TypeName("faceSource");
+
+
+    //- Construct from components
+    faceSource
+    (
+        const word& name,
+        const objectRegistry& obr,
+        const dictionary& dict,
+        const bool loadFromFiles = false
+    );
+
+
+    //- Destructor
+    virtual ~faceSource();
+
+
+    // Public member functions
+
+        // Access
+
+            //- Return the source type
+            inline const sourceType& source() const;
+
+            //- Return the local list of face IDs
+            inline const labelList& faceId() const;
+
+            //- Return the local list of patch ID per face
+            inline const labelList& facePatch() const;
+
+            //- Return the list of +1/-1 representing face flip map
+            inline const labelList& flipMap() const;
+
+
+        // Function object functions
+
+            //- Read from dictionary
+            virtual void read(const dictionary&);
+
+            //- Calculate and write
+            virtual void write();
+
+            //- Templated helper function to output field values
+            template<class Type>
+            bool writeValues(const word& fieldName);
+
+            //- Filter a surface field according to faceIds
+            template<class Type>
+            tmp<Field<Type> > filterField
+            (
+                const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+            ) const;
+
+            //- Filter a volume field according to faceIds
+            template<class Type>
+            tmp<Field<Type> > filterField
+            (
+                const GeometricField<Type, fvPatchField, volMesh>& field
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fieldValues
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "faceSourceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "faceSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.C

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "faceSourceFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug
+    (
+        faceSourceFunctionObject,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        faceSourceFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+Typedef
+    Foam::faceSourceFunctionObject
+
+Description
+    FunctionObject wrapper around faceSource to allow it to be
+    created via the functions list within controlDict.
+
+SourceFiles
+    faceSourceFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceSourceFunctionObject_H
+#define faceSourceFunctionObject_H
+
+#include "faceSource.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<fieldValues::faceSource>
+        faceSourceFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/face/faceSourceI.H

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "faceSource.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::fieldValues::faceSource::sourceType&
+Foam::fieldValues::faceSource::source() const
+{
+    return source_;
+}
+
+
+inline const Foam::labelList&
+Foam::fieldValues::faceSource::faceId() const
+{
+    return faceId_;
+}
+
+
+inline const Foam::labelList&
+Foam::fieldValues::faceSource::facePatch() const
+{
+    return facePatchId_;
+}
+
+
+inline const Foam::labelList&
+Foam::fieldValues::faceSource::flipMap() const
+{
+    return flipMap_;
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/face/faceSourceTemplates.C

@@ -0,0 +1,269 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "faceSource.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "IOList.H"
+#include "ListListOps.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::fieldValues::faceSource::setFieldValues
+(
+    const word& fieldName,
+    List<Type>& values
+) const
+{
+    values.setSize(faceId_.size(), pTraits<Type>::zero);
+
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
+    typedef GeometricField<Type, fvPatchField, volMesh> vf;
+
+    if (obr_.foundObject<sf>(fieldName))
+    {
+        const sf& field = obr_.lookupObject<sf>(fieldName);
+
+        forAll(values, i)
+        {
+            label faceI = faceId_[i];
+            label patchI = facePatchId_[i];
+            if (patchI >= 0)
+            {
+                values[i] = field.boundaryField()[patchI][faceI];
+            }
+            else
+            {
+                values[i] = field[faceI];
+            }
+
+            values[i] *= flipMap_[i];
+        }
+
+        return true;
+    }
+    else if (obr_.foundObject<vf>(fieldName))
+    {
+        const vf& field = obr_.lookupObject<vf>(fieldName);
+
+        forAll(values, i)
+        {
+            label faceI = faceId_[i];
+            label patchI = facePatchId_[i];
+            if (patchI >= 0)
+            {
+                values[i] = field.boundaryField()[patchI][faceI];
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "fieldValues::faceSource::setFieldValues"
+                    "("
+                        "const word&, "
+                        "List<Type>&"
+                    ") const"
+                )   << type() << " " << name_ << ": "
+                    << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
+                    << nl
+                    << "    Unable to process internal faces for volume field "
+                    << fieldName << nl << abort(FatalError);
+            }
+
+            values[i] *= flipMap_[i];
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class Type>
+Type Foam::fieldValues::faceSource::processValues
+(
+    const List<Type>& values
+) const
+{
+    Type result = pTraits<Type>::zero;
+    switch (operation_)
+    {
+        case opSum:
+        {
+            result = sum(values);
+            break;
+        }
+        case opAreaAverage:
+        {
+            tmp<scalarField> magSf = filterField(mesh().magSf());
+            result = sum(values*magSf())/sum(magSf());
+            break;
+        }
+        case opAreaIntegrate:
+        {
+            result = sum(values*filterField(mesh().magSf()));
+            break;
+        }
+        default:
+        {
+            // Do nothing
+        }
+    }
+
+    return result;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
+{
+    List<List<Type> > allValues(Pstream::nProcs());
+
+    bool validField =
+        setFieldValues<Type>(fieldName, allValues[Pstream::myProcNo()]);
+
+    if (validField)
+    {
+        Pstream::gatherList(allValues);
+
+        if (Pstream::master())
+        {
+            List<Type> values =
+                ListListOps::combine<List<Type> >
+                (
+                    allValues,
+                    accessOp<List<Type> >()
+                );
+
+            Type result = processValues(values);
+
+            if (valueOutput_)
+            {
+                IOList<Type>
+                (
+                    IOobject
+                    (
+                        fieldName + "_" + sourceTypeNames_[source_] + "-"
+                            + sourceName_,
+                        obr_.time().timeName(),
+                        obr_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    values
+                ).write();
+            }
+
+
+            outputFilePtr_()<< tab << result;
+
+            if (log_)
+            {
+                Info<< "    " << operationTypeNames_[operation_]
+                    << "(" << sourceName_ << ") for " << fieldName
+                    <<  " = " << result << endl;
+            }
+        }
+    }
+
+    return validField;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field
+) const
+{
+    tmp<Field<Type> > tvalues(new Field<Type>(faceId_.size()));
+    Field<Type>& values = tvalues();
+
+    forAll(values, i)
+    {
+        label faceI = faceId_[i];
+        label patchI = facePatchId_[i];
+        if (patchI >= 0)
+        {
+            values[i] = field.boundaryField()[patchI][faceI];
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "fieldValues::faceSource::filterField"
+                "("
+                    "const GeometricField<Type, fvPatchField, volMesh>&"
+                ") const"
+            )   << type() << " " << name_ << ": "
+                << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
+                << nl
+                << "    Unable to process internal faces for volume field "
+                << field.name() << nl << abort(FatalError);
+        }
+
+        values[i] *= flipMap_[i];
+    }
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+) const
+{
+    tmp<Field<Type> > tvalues(new Field<Type>(faceId_.size()));
+    Field<Type>& values = tvalues();
+
+    forAll(values, i)
+    {
+        label faceI = faceId_[i];
+        label patchI = facePatchId_[i];
+        if (patchI >= 0)
+        {
+            values[i] = field.boundaryField()[patchI][faceI];
+        }
+        else
+        {
+            values[i] = field[faceI];
+        }
+
+        values[i] *= flipMap_[i];
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
+

src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "fieldValue.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fieldValue, 0);
+
+    defineTemplateTypeNameAndDebug(IOList<vector>, 0);
+    defineTemplateTypeNameAndDebug(IOList<sphericalTensor>, 0);
+    defineTemplateTypeNameAndDebug(IOList<symmTensor>, 0);
+    defineTemplateTypeNameAndDebug(IOList<tensor>, 0);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::fieldValue::updateMesh(const mapPolyMesh&)
+{
+    // Do nothing
+}
+
+
+void Foam::fieldValue::movePoints(const Field<point>&)
+{
+    // Do nothing
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fieldValue::fieldValue
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    name_(name),
+    obr_(obr),
+    active_(true),
+    log_(false),
+    sourceName_(dict.lookup("sourceName")),
+    fields_(dict.lookup("fields")),
+    valueOutput_(dict.lookup("valueOutput"))
+{
+    // Only active if obr is an fvMesh
+    if (isA<fvMesh>(obr_))
+    {
+        read(dict);
+    }
+    else
+    {
+        WarningIn
+        (
+            "fieldValue::fieldValue"
+            "("
+                "const word&, "
+                "const objectRegistry&, "
+                "const dictionary&, "
+                "const bool"
+            ")"
+        )   << "No fvMesh available, deactivating."
+            << nl << endl;
+        active_ = false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fieldValue::~fieldValue()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::fieldValue::name() const
+{
+    return name_;
+}
+
+
+const Foam::objectRegistry& Foam::fieldValue::obr() const
+{
+    return obr_;
+}
+
+
+bool Foam::fieldValue::active() const
+{
+    return active_;
+}
+
+
+const Foam::Switch& Foam::fieldValue::log() const
+{
+    return log_;
+}
+
+
+const Foam::word& Foam::fieldValue::sourceName() const
+{
+    return sourceName_;
+}
+
+
+const Foam::wordList& Foam::fieldValue::fields() const
+{
+    return fields_;
+}
+
+
+const Foam::Switch& Foam::fieldValue::valueOutput() const
+{
+    return valueOutput_;
+}
+
+
+const Foam::fvMesh& Foam::fieldValue::mesh() const
+{
+    return refCast<const fvMesh>(obr_);
+}
+
+
+void Foam::fieldValue::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        log_ = dict.lookupOrDefault<Switch>("log", false);
+        dict.lookup("fields") >> fields_;
+        dict.lookup("valueOutput") >> valueOutput_;
+    }
+}
+
+
+void Foam::fieldValue::execute()
+{
+    // Do nothing
+}
+
+
+void Foam::fieldValue::end()
+{
+    // Do nothing
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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::fieldValue
+
+Description
+    Base class for field value -based function objects.
+
+SourceFiles
+    fieldValue.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fieldValue_H
+#define fieldValue_H
+
+#include "Switch.H"
+#include "pointFieldFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class dictionary;
+class objectRegistry;
+class fvMesh;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class fieldValue Declaration
+\*---------------------------------------------------------------------------*/
+
+class fieldValue
+{
+
+protected:
+
+    // Protected data
+
+         //- Name of this fieldValue object
+        word name_;
+
+        //- Database this class is registered to
+        const objectRegistry& obr_;
+
+        //- Active flag
+        bool active_;
+
+        //- Switch to send output to Info as well as to file
+        Switch log_;
+
+        //- Name of source object
+        word sourceName_;
+
+        //- List of field names to operate on
+        wordList fields_;
+
+        //- Output field values flag
+        Switch valueOutput_;
+
+
+        // Functions to be over-ridden from IOoutputFilter class
+
+            //- Update mesh
+            virtual void updateMesh(const mapPolyMesh&);
+
+            //- Move points
+            virtual void movePoints(const Field<point>&);
+
+
+public:
+
+    //- Run-time type information
+    TypeName("fieldValue");
+
+
+    //- Construct from components
+    fieldValue
+    (
+        const word& name,
+        const objectRegistry& obr,
+        const dictionary& dict,
+        const bool loadFromFiles = false
+    );
+
+
+    //- Destructor
+    virtual ~fieldValue();
+
+
+    // Public member functions
+
+        // Access
+
+            //- Return the name of the geometric source
+            const word& name() const;
+
+            //- Return the reference to the object registry
+            const objectRegistry& obr() const;
+
+            //- Return the active flag
+            bool active() const;
+
+            //- Return the switch to send output to Info as well as to file
+            const Switch& log() const;
+
+            //- Return the source name
+            const word& sourceName() const;
+
+            //- Return the list of field names
+            const wordList& fields() const;
+
+            //- Return the output field values flag
+            const Switch& valueOutput() const;
+
+            //- Helper funvction to return the reference to the mesh
+            const fvMesh& mesh() const;
+
+
+        // Function object functions
+
+            //- Read from dictionary
+            virtual void read(const dictionary& dict);
+
+            //- Execute
+            virtual void execute();
+
+            //- Execute the at the final time-loop, currently does nothing
+            virtual void end();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLine.C

@@ -0,0 +1,622 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Pstream.H"
+#include "functionObjectList.H"
+#include "streamLine.H"
+#include "fvMesh.H"
+#include "streamLineParticleCloud.H"
+#include "ReadFields.H"
+#include "meshSearch.H"
+#include "sampledSet.H"
+#include "globalIndex.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::streamLine, 0);
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::streamLine::track()
+{
+    const Time& runTime = const_cast<Time&>(obr_.time());
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+    IDLList<streamLineParticle> initialParticles;
+    streamLineParticleCloud particles
+    (
+        mesh,
+        cloudName_,
+        initialParticles
+    );
+
+    //Pout<< "Seeding particles." << endl;
+
+    const sampledSet& seedPoints = sampledSetPtr_();
+    forAll(seedPoints, i)
+    {
+        //Pout<< "Seeded particle at " << seedPoints[i]
+        //    << " at cell:" << seedPoints.cells()[i]
+        //    << endl;
+
+        particles.addParticle
+        (
+            new streamLineParticle
+            (
+                particles,
+                seedPoints[i],
+                seedPoints.cells()[i],
+                lifeTime_               // lifetime
+            )
+        );
+    }
+    label nSeeds = returnReduce(particles.size(), sumOp<label>());
+    Info<< "Seeded " << nSeeds
+        << " particles." << endl;
+
+
+    // Read or lookup fields
+    PtrList<volScalarField> vsFlds;
+    PtrList<interpolationCellPoint<scalar> > vsInterp;
+    PtrList<volVectorField> vvFlds;
+    PtrList<interpolationCellPoint<vector> > vvInterp;
+
+
+    label UIndex = -1;
+
+    if (loadFromFiles_)
+    {
+        IOobjectList allObjects(mesh, runTime.timeName());
+
+        IOobjectList objects(2*fields_.size());
+        forAll(fields_, i)
+        {
+            objects.add(*allObjects[fields_[i]]);
+        }
+
+        ReadFields(mesh, objects, vsFlds);
+        vsInterp.setSize(vsFlds.size());
+        forAll(vsFlds, i)
+        {
+            vsInterp.set(i, new interpolationCellPoint<scalar>(vsFlds[i]));
+        }
+        ReadFields(mesh, objects, vvFlds);
+        vvInterp.setSize(vvFlds.size());
+        forAll(vvFlds, i)
+        {
+            vvInterp.set(i, new interpolationCellPoint<vector>(vvFlds[i]));
+        }
+    }
+    else
+    {
+        label nScalar = 0;
+        label nVector = 0;
+
+        forAll(fields_, i)
+        {
+            if (mesh.foundObject<volScalarField>(fields_[i]))
+            {
+                nScalar++;
+            }
+            else if (mesh.foundObject<volVectorField>(fields_[i]))
+            {
+                nVector++;
+            }
+            else
+            {
+                FatalErrorIn("streamLine::execute()")
+                    << "Cannot find field " << fields_[i] << endl
+                    << "Valid scalar fields are:"
+                    << mesh.names(volScalarField::typeName) << endl
+                    << "Valid vector fields are:"
+                    << mesh.names(volVectorField::typeName)
+                    << exit(FatalError);
+            }
+        }
+        vsInterp.setSize(nScalar);
+        nScalar = 0;
+        vvInterp.setSize(nVector);
+        nVector = 0;
+
+        forAll(fields_, i)
+        {
+            if (mesh.foundObject<volScalarField>(fields_[i]))
+            {
+                const volScalarField& f = mesh.lookupObject<volScalarField>
+                (
+                    fields_[i]
+                );
+                vsInterp.set
+                (
+                    nScalar++,
+                    new interpolationCellPoint<scalar>(f)
+                );
+            }
+            else if (mesh.foundObject<volVectorField>(fields_[i]))
+            {
+                const volVectorField& f = mesh.lookupObject<volVectorField>
+                (
+                    fields_[i]
+                );
+
+                if (f.name() == UName_)
+                {
+                    UIndex = nVector;
+                }
+
+                vvInterp.set
+                (
+                    nVector++,
+                    new interpolationCellPoint<vector>(f)
+                );
+            }
+        }
+    }
+
+    // Store the names
+    scalarNames_.setSize(vsInterp.size());
+    forAll(vsInterp, i)
+    {
+        scalarNames_[i] = vsInterp[i].psi().name();
+    }
+    vectorNames_.setSize(vvInterp.size());
+    forAll(vvInterp, i)
+    {
+        vectorNames_[i] = vvInterp[i].psi().name();
+    }
+
+
+    // Check that we know the index of U in the interpolators.
+
+    if (UIndex == -1)
+    {
+        FatalErrorIn("streamLine::execute()")
+            << "Cannot find field to move particles with : " << UName_
+            << endl
+            << "This field has to be present in the sampled fields "
+            << fields_
+            << " and in the objectRegistry." << endl
+            << exit(FatalError);
+    }
+
+
+
+    // Sampled data
+    // ~~~~~~~~~~~~
+
+    // Size to maximum expected sizes.
+    allTracks_.clear();
+    allTracks_.setCapacity(nSeeds);
+    allScalars_.setSize(vsInterp.size());
+    forAll(allScalars_, i)
+    {
+        allScalars_[i].clear();
+        allScalars_[i].setCapacity(nSeeds);
+    }
+    allVectors_.setSize(vvInterp.size());
+    forAll(allVectors_, i)
+    {
+        allVectors_[i].clear();
+        allVectors_[i].setCapacity(nSeeds);
+    }
+
+    // additional particle info
+    streamLineParticle::trackData td
+    (
+        particles,
+        vsInterp,
+        vvInterp,
+        UIndex,         // index of U in vvInterp
+        trackForward_,  // track in +u direction
+        allTracks_,
+        allScalars_,
+        allVectors_
+    );
+
+    // Track
+    //Pout<< "Tracking particles." << endl;
+    particles.move(td);
+    //Pout<< "Finished tracking particles." << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::streamLine::streamLine
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    dict_(dict),
+    name_(name),
+    obr_(obr),
+    loadFromFiles_(loadFromFiles),
+    active_(true)
+{
+    // Only active if a fvMesh is available
+    if (isA<fvMesh>(obr_))
+    {
+        read(dict_);
+    }
+    else
+    {
+        active_ = false;
+        WarningIn
+        (
+            "streamLine::streamLine\n"
+            "(\n"
+                "const word&,\n"
+                "const objectRegistry&,\n"
+                "const dictionary&,\n"
+                "const bool\n"
+            ")"
+        )   << "No fvMesh available, deactivating."
+            << nl << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::streamLine::~streamLine()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::streamLine::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        //dict_ = dict;
+        dict.lookup("fields") >> fields_;
+        UName_ = dict.lookupOrDefault<word>("U", "U");
+        dict.lookup("trackForward") >> trackForward_;
+        dict.lookup("lifeTime") >> lifeTime_;
+        cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine");
+        dict.lookup("seedSampleSet") >> seedSet_;
+
+        const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+        meshSearchPtr_.reset(new meshSearch(mesh, false));
+
+        const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
+        sampledSetPtr_ = sampledSet::New
+        (
+            seedSet_,
+            mesh,
+            meshSearchPtr_(),
+            coeffsDict
+        );
+        coeffsDict.lookup("axis") >> sampledSetAxis_;
+
+        scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
+        vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
+    }
+}
+
+
+void Foam::streamLine::execute()
+{
+//    const Time& runTime = const_cast<Time&>(obr_.time());
+//    Pout<< "**streamLine::execute : time:" << runTime.timeName() << endl;
+//    Pout<< "**streamLine::execute : time:" << runTime.timeIndex() << endl;
+//
+//    bool isOutputTime = false;
+//
+//    const functionObjectList& fobs = runTime.functionObjects();
+//
+//    forAll(fobs, i)
+//    {
+//        if (isA<streamLineFunctionObject>(fobs[i]))
+//        {
+//            const streamLineFunctionObject& fo =
+//                dynamic_cast<const streamLineFunctionObject&>(fobs[i]);
+//
+//            if (fo.name() == name_)
+//            {
+//                Pout<< "found me:" << i << endl;
+//                if (fo.outputControl().output())
+//                {
+//                    isOutputTime = true;
+//                    break;
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    if (active_ && isOutputTime)
+//    {
+//        track();
+//    }
+}
+
+
+void Foam::streamLine::end()
+{
+}
+
+
+void Foam::streamLine::write()
+{
+    if (active_)
+    {
+        const Time& runTime = const_cast<Time&>(obr_.time());
+        const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+
+        // Do all injection and tracking
+        track();
+
+
+        if (Pstream::parRun())
+        {
+            // Append slave tracks to master ones
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            globalIndex globalTrackIDs(allTracks_.size());
+
+            // Construct a distribution map to pull all to the master.
+            labelListList sendMap(Pstream::nProcs());
+            labelListList recvMap(Pstream::nProcs());
+
+            if (Pstream::master())
+            {
+                // Master: receive all. My own first, then consecutive
+                // processors.
+                label trackI = 0;
+
+                forAll(recvMap, procI)
+                {
+                    labelList& fromProc = recvMap[procI];
+                    fromProc.setSize(globalTrackIDs.localSize(procI));
+                    forAll(fromProc, i)
+                    {
+                        fromProc[i] = trackI++;
+                    }
+                }
+            }
+
+            labelList& toMaster = sendMap[0];
+            toMaster.setSize(globalTrackIDs.localSize());
+            forAll(toMaster, i)
+            {
+                toMaster[i] = i;
+            }
+
+            const mapDistribute distMap
+            (
+                globalTrackIDs.size(),
+                sendMap.xfer(),
+                recvMap.xfer()
+            );
+
+
+            // Distribute the track positions
+            mapDistribute::distribute
+            (
+                Pstream::scheduled,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                allTracks_
+            );
+
+            // Distribute the scalars
+            forAll(allScalars_, scalarI)
+            {
+                mapDistribute::distribute
+                (
+                    Pstream::scheduled,
+                    distMap.schedule(),
+                    distMap.constructSize(),
+                    distMap.subMap(),
+                    distMap.constructMap(),
+                    allScalars_[scalarI]
+                );
+            }
+            // Distribute the vectors
+            forAll(allVectors_, vectorI)
+            {
+                mapDistribute::distribute
+                (
+                    Pstream::scheduled,
+                    distMap.schedule(),
+                    distMap.constructSize(),
+                    distMap.subMap(),
+                    distMap.constructMap(),
+                    allVectors_[vectorI]
+                );
+            }
+        }
+
+
+        label n = 0;
+        forAll(allTracks_, trackI)
+        {
+            n += allTracks_[trackI].size();
+        }
+
+        Info<< "Tracks:" << allTracks_.size()
+            << "  total samples:" << n << endl;
+
+
+        // Massage into form suitable for writers
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        if (Pstream::master() && allTracks_.size())
+        {
+            // Make output directory
+
+            fileName vtkPath
+            (
+                Pstream::parRun()
+              ? runTime.path()/".."/"sets"/name()
+              : runTime.path()/"sets"/name()
+            );
+            if (mesh.name() != fvMesh::defaultRegion)
+            {
+                vtkPath = vtkPath/mesh.name();
+            }
+            vtkPath = vtkPath/mesh.time().timeName();
+
+            mkDir(vtkPath);
+
+
+            // Convert track positions
+
+            PtrList<coordSet> tracks(allTracks_.size());
+            forAll(allTracks_, trackI)
+            {
+                tracks.set
+                (
+                    trackI,
+                    new coordSet
+                    (
+                        "track" + Foam::name(trackI),
+                        sampledSetAxis_                 //"xyz"
+                    )
+                );
+                tracks[trackI].transfer(allTracks_[trackI]);
+            }
+
+
+
+            // Convert scalar values
+
+            if (allScalars_.size() > 0)
+            {
+                List<List<scalarField> > scalarValues(allScalars_.size());
+
+                forAll(allScalars_, scalarI)
+                {
+                    DynamicList<scalarList>& allTrackVals =
+                        allScalars_[scalarI];
+                    scalarValues[scalarI].setSize(allTrackVals.size());
+
+                    forAll(allTrackVals, trackI)
+                    {
+                        scalarList& trackVals = allTrackVals[trackI];
+                        scalarValues[scalarI][trackI].transfer(trackVals);
+                    }
+                }
+
+                fileName vtkFile
+                (
+                    vtkPath
+                  / scalarFormatterPtr_().getFileName
+                    (
+                        tracks[0],
+                        scalarNames_
+                    )
+                );
+
+                Info<< "Writing data to " << vtkFile.path() << endl;
+
+                scalarFormatterPtr_().write
+                (
+                    true,           // writeTracks
+                    tracks,
+                    scalarNames_,
+                    scalarValues,
+                    OFstream(vtkFile)()
+                );
+            }
+
+            // Convert vector values
+
+            if (allVectors_.size() > 0)
+            {
+                List<List<vectorField> > vectorValues(allVectors_.size());
+
+                forAll(allVectors_, vectorI)
+                {
+                    DynamicList<vectorList>& allTrackVals =
+                        allVectors_[vectorI];
+                    vectorValues[vectorI].setSize(allTrackVals.size());
+
+                    forAll(allTrackVals, trackI)
+                    {
+                        vectorList& trackVals = allTrackVals[trackI];
+                        vectorValues[vectorI][trackI].transfer(trackVals);
+                    }
+                }
+
+                fileName vtkFile
+                (
+                    vtkPath
+                  / vectorFormatterPtr_().getFileName
+                    (
+                        tracks[0],
+                        vectorNames_
+                    )
+                );
+
+                //Info<< "Writing vector data to " << vtkFile << endl;
+
+                vectorFormatterPtr_().write
+                (
+                    true,           // writeTracks
+                    tracks,
+                    vectorNames_,
+                    vectorValues,
+                    OFstream(vtkFile)()
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::streamLine::updateMesh(const mapPolyMesh&)
+{
+    read(dict_);
+}
+
+
+void Foam::streamLine::movePoints(const pointField&)
+{
+    // Moving mesh affects the search tree
+    read(dict_);
+}
+
+
+//void Foam::streamLine::readUpdate(const polyMesh::readUpdateState state)
+//{
+//    if (state != UNCHANGED)
+//    {
+//        read(dict_);
+//    }
+//}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLine.H

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::streamLine
+
+Description
+    Generation of streamlines. Samples along track of passive particle.
+
+SourceFiles
+    streamLine.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef streamLine_H
+#define streamLine_H
+
+#include "volFieldsFwd.H"
+#include "pointFieldFwd.H"
+#include "Switch.H"
+#include "DynamicList.H"
+#include "scalarList.H"
+#include "vectorList.H"
+#include "polyMesh.H"
+#include "writer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+class meshSearch;
+class sampledSet;
+
+/*---------------------------------------------------------------------------*\
+                         Class streamLine Declaration
+\*---------------------------------------------------------------------------*/
+
+class streamLine
+{
+    // Private data
+
+        //- Input dictionary
+        dictionary dict_;
+
+        //- Name of this set of field averages.
+        word name_;
+
+        //- Database this class is registered to
+        const objectRegistry& obr_;
+
+        //- Load fields from files (not from objectRegistry)
+        bool loadFromFiles_;
+
+        //- On/off switch
+        bool active_;
+
+
+        //- List of fields to sample
+        wordList fields_;
+
+        //- Field to transport particle with
+        word UName_;
+
+        //- Whether to use +u or -u
+        bool trackForward_;
+
+        //- Maximum lifetime (= number of cells) of particle
+        label lifeTime_;
+
+        //- Optional specified name of particles
+        word cloudName_;
+
+        //- Type of seed
+        word seedSet_;
+
+        //- Names of scalar fields
+        wordList scalarNames_;
+
+        //- Names of vector fields
+        wordList vectorNames_;
+
+
+
+        // Demand driven
+
+            //- Mesh searching enigne
+            autoPtr<meshSearch> meshSearchPtr_;
+
+            //- Seed set engine
+            autoPtr<sampledSet> sampledSetPtr_;
+
+            //- Axis of the sampled points to output
+            word sampledSetAxis_;
+
+            //- File output writer
+            autoPtr<writer<scalar> > scalarFormatterPtr_;
+
+            autoPtr<writer<vector> > vectorFormatterPtr_;
+
+
+        // Generated data
+
+            //- All tracks. Per particle the points it passed through
+            DynamicList<List<point> > allTracks_;
+
+            //- Per scalarField, per particle, the sampled value.
+            List<DynamicList<scalarList> > allScalars_;
+
+            //- Per scalarField, per particle, the sampled value.
+            List<DynamicList<vectorList> > allVectors_;
+
+
+
+        //- Do all seeding and tracking
+        void track();
+
+        //- Disallow default bitwise copy construct
+        streamLine(const streamLine&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const streamLine&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("streamLine");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        streamLine
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    //- Destructor
+
+        virtual ~streamLine();
+
+
+    // Member Functions
+
+        //- Return name of the set of field averages
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the field average data
+        virtual void read(const dictionary&);
+
+        //- Execute the averaging
+        virtual void execute();
+
+        //- Execute the averaging at the final time-loop, currently does nothing
+        virtual void end();
+
+        //- Calculate the field average data and write
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const pointField&);
+
+        ////- Update for changes of mesh due to readUpdate
+        //virtual void readUpdate(const polyMesh::readUpdateState state);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//#ifdef NoRepository
+//#   include "streamLineTemplates.C"
+//#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineFunctionObject.C

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "streamLineFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(streamLineFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        streamLineFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineFunctionObject.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Typedef
+    Foam::streamLineFunctionObject
+
+Description
+    FunctionObject wrapper around streamLines to allow them to be created via
+    the functions list within controlDict.
+
+SourceFiles
+    streamLineFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef streamLineFunctionObject_H
+#define streamLineFunctionObject_H
+
+#include "streamLine.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<streamLine>
+        streamLineFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C

@@ -0,0 +1,486 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "streamLineParticle.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(streamLineParticle, 0);
+    defineParticleTypeNameAndDebug(streamLineParticle, 0);
+    defineTemplateTypeNameAndDebugWithName
+    (
+        IOField<vectorField>,
+        "vectorFieldField",
+        0
+    );
+};
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::vector Foam::streamLineParticle::interpolateFields
+(
+    const streamLineParticle::trackData& td,
+    const point& position,
+    const label cellI
+)
+{
+    if (cellI == -1)
+    {
+        FatalErrorIn("streamLineParticle::interpolateFields(..)")
+            << "Cell:" << cellI << abort(FatalError);
+    }
+
+    sampledScalars_.setSize(td.vsInterp_.size());
+    forAll(td.vsInterp_, scalarI)
+    {
+        sampledScalars_[scalarI].append
+        (
+            td.vsInterp_[scalarI].interpolate
+            (
+                position,
+                cellI
+            )
+        );
+    }
+
+    sampledVectors_.setSize(td.vvInterp_.size());
+    forAll(td.vvInterp_, vectorI)
+    {
+        sampledVectors_[vectorI].append
+        (
+            td.vvInterp_[vectorI].interpolate
+            (
+                position,
+                cellI
+            )
+        );
+    }
+
+    const DynamicList<vector>& U = sampledVectors_[td.UIndex_];
+
+    return U[U.size()-1];
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct from components
+Foam::streamLineParticle::streamLineParticle
+(
+    const Cloud<streamLineParticle>& c,
+    const vector& position,
+    const label celli,
+    const label lifeTime
+)
+:
+    Particle<streamLineParticle>(c, position, celli),
+    lifeTime_(lifeTime)
+{}
+
+
+//- Construct from Istream
+Foam::streamLineParticle::streamLineParticle
+(
+    const Cloud<streamLineParticle>& c,
+    Istream& is,
+    bool readFields
+)
+:
+    Particle<streamLineParticle>(c, is, readFields)
+{
+    if (readFields)
+    {
+        //if (is.format() == IOstream::ASCII)
+        List<scalarList> sampledScalars;
+        List<vectorList> sampledVectors;
+
+        is  >> lifeTime_ >> sampledPositions_ >> sampledScalars
+            >> sampledVectors;
+
+        sampledScalars_.setSize(sampledScalars.size());
+        forAll(sampledScalars, i)
+        {
+            sampledScalars_[i].transfer(sampledScalars[i]);
+        }
+        sampledVectors_.setSize(sampledVectors.size());
+        forAll(sampledVectors, i)
+        {
+            sampledVectors_[i].transfer(sampledVectors[i]);
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "streamLineParticle::streamLineParticle"
+        "(const Cloud<streamLineParticle>&, Istream&, bool)"
+    );
+}
+
+
+// Construct copy
+Foam::streamLineParticle::streamLineParticle
+(
+    const streamLineParticle& c
+)
+:
+    Particle<streamLineParticle>(c),
+    lifeTime_(c.lifeTime_),
+    sampledPositions_(c.sampledPositions_),
+    sampledScalars_(c.sampledScalars_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    scalar deltaT = GREAT;  //cloud().pMesh().time().deltaT().value();
+    scalar tEnd = (1.0 - stepFraction())*deltaT;
+    scalar dtMax = tEnd;
+
+    while
+    (
+        td.keepParticle
+    && !td.switchProcessor
+    && lifeTime_ > 0
+    && tEnd > ROOTVSMALL
+    )
+    {
+        // set the lagrangian time-step
+        scalar dt = min(dtMax, tEnd);
+
+        // Store current position and sampled velocity.
+        --lifeTime_;
+        sampledPositions_.append(position());
+        vector U = interpolateFields(td, position(), cell());
+
+        if (!td.trackForward_)
+        {
+            U = -U;
+        }
+
+        dt *= trackToFace(position()+dt*U, td);
+
+        tEnd -= dt;
+        stepFraction() = 1.0 - tEnd/deltaT;
+    }
+
+    if (!td.keepParticle || lifeTime_ == 0)
+    {
+        if (lifeTime_ == 0)
+        {
+            if (debug)
+            {
+                Pout<< "streamLineParticle : Removing stagnant particle:"
+                    << static_cast<Particle<streamLineParticle> >(*this)
+                    << " sampled positions:" << sampledPositions_.size()
+                    << endl;
+            }
+            td.keepParticle = false;
+        }
+        else
+        {
+            // Normal exit. Store last position and fields
+            sampledPositions_.append(position());
+            interpolateFields(td, position(), cell());
+
+            if (debug)
+            {
+                Pout<< "streamLineParticle : Removing particle:"
+                    << static_cast<Particle<streamLineParticle> >(*this)
+                    << " sampled positions:" << sampledPositions_.size()
+                    << endl;
+            }
+        }
+
+        // Transfer particle data into trackData.
+        //td.allPositions_.append(sampledPositions_);
+        td.allPositions_.append(vectorList());
+        vectorList& top = td.allPositions_[td.allPositions_.size()-1];
+        top.transfer(sampledPositions_);
+
+        forAll(sampledScalars_, i)
+        {
+            //td.allScalars_[i].append(sampledScalars_[i]);
+            td.allScalars_[i].append(scalarList());
+            scalarList& top = td.allScalars_[i][td.allScalars_[i].size()-1];
+            top.transfer(sampledScalars_[i]);
+        }
+        forAll(sampledVectors_, i)
+        {
+            //td.allVectors_[i].append(sampledVectors_[i]);
+            td.allVectors_[i].append(vectorList());
+            vectorList& top = td.allVectors_[i][td.allVectors_[i].size()-1];
+            top.transfer(sampledVectors_[i]);
+        }
+    }
+
+    return td.keepParticle;
+}
+
+
+bool Foam::streamLineParticle::hitPatch
+(
+    const polyPatch&,
+    streamLineParticle::trackData& td,
+    const label patchI
+)
+{
+    // Disable generic patch interaction
+    return false;
+}
+
+
+bool Foam::streamLineParticle::hitPatch
+(
+    const polyPatch&,
+    int&,
+    const label
+)
+{
+    // Disable generic patch interaction
+    return false;
+}
+
+
+void Foam::streamLineParticle::hitWedgePatch
+(
+    const wedgePolyPatch& pp,
+    streamLineParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::streamLineParticle::hitWedgePatch
+(
+    const wedgePolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch& pp,
+    streamLineParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::streamLineParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch& pp,
+    streamLineParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::streamLineParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    streamLineParticle::trackData& td
+)
+{
+    // Switch particle
+    td.switchProcessor = true;
+}
+
+
+void Foam::streamLineParticle::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    streamLineParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::streamLineParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::hitPatch
+(
+    const polyPatch& wpp,
+    streamLineParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::streamLineParticle::hitPatch
+(
+    const polyPatch& wpp,
+    int&
+)
+{}
+
+
+void Foam::streamLineParticle::readFields
+(
+    Cloud<streamLineParticle>& c
+)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    IOField<label> lifeTime
+    (
+        c.fieldIOobject("lifeTime", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, lifeTime);
+
+    IOField<pointField> sampledPositions
+    (
+        c.fieldIOobject("sampledPositions", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, sampledPositions);
+
+//    IOField<pointField> sampleVelocity
+//    (
+//        c.fieldIOobject("sampleVelocity", IOobject::MUST_READ)
+//    );
+//    c.checkFieldIOobject(c, sampleVelocity);
+
+    label i = 0;
+    forAllIter(Cloud<streamLineParticle>, c, iter)
+    {
+        iter().lifeTime_ = lifeTime[i];
+        iter().sampledPositions_.transfer(sampledPositions[i]);
+//        iter().sampleVelocity_.transfer(sampleVelocity[i]);
+        i++;
+    }
+}
+
+
+void Foam::streamLineParticle::writeFields
+(
+    const Cloud<streamLineParticle>& c
+)
+{
+    Particle<streamLineParticle>::writeFields(c);
+
+    label np =  c.size();
+
+    IOField<label> lifeTime
+    (
+        c.fieldIOobject("lifeTime", IOobject::NO_READ),
+        np
+    );
+    IOField<pointField> sampledPositions
+    (
+        c.fieldIOobject("sampledPositions", IOobject::NO_READ),
+        np
+    );
+//    IOField<pointField> sampleVelocity
+//    (
+//        c.fieldIOobject("sampleVelocity", IOobject::NO_READ),
+//        np
+//    );
+
+    label i = 0;
+    forAllConstIter(Cloud<streamLineParticle>, c, iter)
+    {
+        lifeTime[i] = iter().lifeTime_;
+        sampledPositions[i] = iter().sampledPositions_;
+//        sampleVelocity[i] = iter().sampleVelocity_;
+        i++;
+    }
+
+    lifeTime.write();
+    sampledPositions.write();
+//    sampleVelocity.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const streamLineParticle& p)
+{
+    os  << static_cast<const Particle<streamLineParticle>&>(p)
+        << token::SPACE << p.lifeTime_
+        << token::SPACE << p.sampledPositions_
+        << token::SPACE << p.sampledScalars_
+        << token::SPACE << p.sampledVectors_;
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const streamLineParticle&)");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H

@@ -0,0 +1,300 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::streamLineParticle
+
+Description
+    Particle class that samples fields as it passes through. Used in streamline
+    calculation.
+
+SourceFiles
+    streamLineParticle.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef streamLineParticle_H
+#define streamLineParticle_H
+
+#include "Particle.H"
+#include "autoPtr.H"
+#include "interpolationCellPoint.H"
+#include "vectorList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class streamLineParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+                           Class streamLineParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class streamLineParticle
+:
+    public Particle<streamLineParticle>
+{
+
+public:
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    :
+        public Particle<streamLineParticle>::trackData
+    {
+
+    public:
+
+
+        const PtrList<interpolationCellPoint<scalar> >& vsInterp_;
+        const PtrList<interpolationCellPoint<vector> >& vvInterp_;
+        const label UIndex_;
+        const bool trackForward_;
+
+        DynamicList<vectorList>& allPositions_;
+        List<DynamicList<scalarList> >& allScalars_;
+        List<DynamicList<vectorList> >& allVectors_;
+
+
+        // Constructors
+
+            trackData
+            (
+                Cloud<streamLineParticle>& cloud,
+                const PtrList<interpolationCellPoint<scalar> >& vsInterp,
+                const PtrList<interpolationCellPoint<vector> >& vvInterp,
+                const label UIndex,
+                const bool trackForward,
+                DynamicList<List<point> >& allPositions,
+                List<DynamicList<scalarList> >& allScalars,
+                List<DynamicList<vectorList> >& allVectors
+            )
+            :
+                Particle<streamLineParticle>::trackData(cloud),
+                vsInterp_(vsInterp),
+                vvInterp_(vvInterp),
+                UIndex_(UIndex),
+                trackForward_(trackForward),
+                allPositions_(allPositions),
+                allScalars_(allScalars),
+                allVectors_(allVectors)
+            {}
+    };
+
+
+private:
+
+    // Private data
+
+        //- Lifetime of particle. Particle dies when reaches 0.
+        label lifeTime_;
+
+        //- sampled positions
+        DynamicList<point> sampledPositions_;
+
+        //- sampled scalars
+        List<DynamicList<scalar> > sampledScalars_;
+
+        //- sampled vectors
+        List<DynamicList<vector> > sampledVectors_;
+
+
+    // Private Member Functions
+
+        //- Interpolate all quantities; return interpolated velocity.
+        vector interpolateFields
+        (
+            const streamLineParticle::trackData&,
+            const point&,
+            const label cellI
+        );
+
+
+public:
+
+    //- Run-time type information
+    TypeName("streamLineParticle");
+
+
+
+    // Constructors
+
+        //- Construct from components
+        streamLineParticle
+        (
+            const Cloud<streamLineParticle>& c,
+            const vector& position,
+            const label celli,
+            const label lifeTime
+        );
+
+        //- Construct from Istream
+        streamLineParticle
+        (
+            const Cloud<streamLineParticle>& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct copy
+        streamLineParticle(const streamLineParticle& c);
+
+        //- Construct and return a clone
+        autoPtr<streamLineParticle> clone() const
+        {
+            return autoPtr<streamLineParticle>
+            (
+                new streamLineParticle(*this)
+            );
+        }
+
+
+    // Member Functions
+
+
+
+        // Tracking
+
+            //- Track all particles to their end point
+            bool move(trackData&);
+
+
+            //- Overridable function to handle the particle hitting a patch
+            //  Executed before other patch-hitting functions
+            bool hitPatch
+            (
+                const polyPatch&,
+                streamLineParticle::trackData& td,
+                const label patchI
+            );
+            bool hitPatch
+            (
+                const polyPatch&,
+                int&,
+                const label patchI
+           );
+
+            //- Overridable function to handle the particle hitting a wedge
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //  symmetryPlane
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a cyclic
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //- processorPatch
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            void hitPatch
+            (
+                const polyPatch&,
+                streamLineParticle::trackData& td
+            );
+            void hitPatch
+            (
+                const polyPatch&,
+                int&
+            );
+
+
+        // I-O
+
+            //- Read
+            static void readFields(Cloud<streamLineParticle>&);
+
+            //- Write
+            static void writeFields(const Cloud<streamLineParticle>&);
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream&, const streamLineParticle&);
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "streamLineParticleCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTemplateTypeNameAndDebug(Cloud<streamLineParticle>, 0);
+
+};
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::streamLineParticleCloud::streamLineParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName
+)
+:
+    Cloud<streamLineParticle>(mesh, cloudName, false)
+{
+    readFields();
+}
+
+
+Foam::streamLineParticleCloud::streamLineParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName,
+    const IDLList<streamLineParticle>& particles
+)
+:
+    Cloud<streamLineParticle>(mesh, cloudName, particles)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::streamLineParticleCloud::readFields()
+{
+    streamLineParticle::readFields(*this);
+}
+
+
+void Foam::streamLineParticleCloud::writeFields() const
+{
+    streamLineParticle::writeFields(*this);
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.H

@@ -23,20 +23,21 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::faceMesh
+    Foam::streamLineCloud
 
 Description
-    Storage for surface mesh i.e. points and faces.
+    A Cloud of streamLine particles
 
 SourceFiles
+    streamLineCloud.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef faceMesh_H
-#define faceMesh_H
+#ifndef streamLineCloud_H
+#define streamLineCloud_H
 
-#include "PrimitivePatch.H"
-#include "faceList.H"
+#include "Cloud.H"
+#include "streamLineParticle.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,62 +45,52 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class faceMesh Declaration
+                           Class streamLineCloud Declaration
 \*---------------------------------------------------------------------------*/
 
-class faceMesh
+class streamLineParticleCloud
 :
-    public PrimitivePatch<face, List, pointField>
+    public Cloud<streamLineParticle>
 {
-    // Private member functions
+    // Private Member Functions
 
-        PrimitivePatch<face, List, pointField> read(Istream& is)
-        {
-            pointField points(is);
-            faceList faces(is);
-            return PrimitivePatch<face, Foam::List, pointField>(faces, points);
-        }
+        //- Disallow default bitwise copy construct
+        streamLineParticleCloud(const streamLineParticleCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const streamLineParticleCloud&);
 
 
 public:
 
+    //- Type of parcel the cloud was instantiated for
+    typedef streamLineParticle parcelType;
+
     // Constructors
 
-        //- Construct from components
-        faceMesh(const faceList& faces, const pointField& points)
-        :
-            PrimitivePatch<face, Foam::List, pointField>(faces, points)
-        {}
+        //- Construct given mesh
+        streamLineParticleCloud
+        (
+            const polyMesh&,
+            const word& cloudName = "defaultCloud"
+        );
 
-        //- Construct from Istream
-        faceMesh(Istream& is)
-        :
-            PrimitivePatch<face, Foam::List, pointField>(read(is))
-        {}
+        //- Construct from mesh, cloud name, and a list of particles
+        streamLineParticleCloud
+        (
+            const polyMesh& mesh,
+            const word& cloudName,
+            const IDLList<streamLineParticle>& particles
+        );
 
 
     // Member Functions
 
-        void flip()
-        {
-            forAll(*this, i)
-            {
-                face& f = operator[](i);
-                f = f.reverseFace();
-            }
-            clearOut();
-        }
+        //- Read fields
+        virtual void readFields();
 
-    // IOstream Operators
-
-        friend Ostream& operator<<(Ostream& os, const faceMesh& fm)
-        {
-            return os
-                << fm.points()
-                << token::NL
-                << static_cast<PrimitivePatch<face, Foam::List, pointField> >
-                   (fm);
-        }
+        //- Write fields
+        virtual void writeFields() const;
 };
 
 

src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C

@@ -55,8 +55,8 @@ Foam::faceZonesIntegration::faceZonesIntegration
     obr_(obr),
     active_(true),
     log_(false),
-    faceZonesSet_(),
-    fItems_(),
+    zoneNames_(),
+    fieldNames_(),
     filePtr_(NULL)
 {
     // Check if the available mesh is an fvMesh otherise deactivate
@@ -94,9 +94,9 @@ void Foam::faceZonesIntegration::read(const dictionary& dict)
     {
         log_ = dict.lookupOrDefault<Switch>("log", false);
 
-        dict.lookup("fields") >> fItems_;
+        dict.lookup("fields") >> fieldNames_;
 
-        dict.lookup("faceZones") >> faceZonesSet_;
+        dict.lookup("faceZones") >> zoneNames_;
     }
 }
 
@@ -132,11 +132,11 @@ void Foam::faceZonesIntegration::makeFile()
             mkDir(faceZonesIntegrationDir);
 
             // Open new file at start up
-            filePtr_.resize(fItems_.size());
+            filePtr_.resize(fieldNames_.size());
 
-            forAll(fItems_, Ifields)
+            forAll(fieldNames_, fieldI)
             {
-                const word& fieldName = fItems_[Ifields];
+                const word& fieldName = fieldNames_[fieldI];
 
                 OFstream* sPtr = new OFstream
                     (
@@ -163,10 +163,9 @@ void Foam::faceZonesIntegration::writeFileHeader()
 
         os  << "#Time " << setw(w);
 
-        forAll (faceZonesSet_, zoneI)
+        forAll (zoneNames_, zoneI)
         {
-            const word name = faceZonesSet_[zoneI];
-            os  << name << setw(w);
+            os  << zoneNames_[zoneI] << setw(w);
         }
 
         os  << nl << endl;
@@ -192,9 +191,9 @@ void Foam::faceZonesIntegration::write()
     {
         makeFile();
 
-        forAll(fItems_, fieldI)
+        forAll(fieldNames_, fieldI)
         {
-            const word& fieldName = fItems_[fieldI];
+            const word& fieldName = fieldNames_[fieldI];
 
             const surfaceScalarField& sField =
                 obr_.lookupObject<surfaceScalarField>(fieldName);
@@ -203,17 +202,17 @@ void Foam::faceZonesIntegration::write()
 
             // 1. integrate over all face zones
 
-            scalarField integralVals(faceZonesSet_.size());
+            scalarField integralVals(zoneNames_.size());
 
-            forAll(faceZonesSet_, setI)
+            forAll(integralVals, zoneI)
             {
-                const word name = faceZonesSet_[setI];
+                const word& name = zoneNames_[zoneI];
 
                 label zoneID = mesh.faceZones().findZoneID(name);
 
                 const faceZone& fZone = mesh.faceZones()[zoneID];
 
-                integralVals[setI] = returnReduce
+                integralVals[zoneI] = returnReduce
                 (
                     calcIntegral(sField, fZone),
                     sumOp<scalar>()
@@ -231,15 +230,15 @@ void Foam::faceZonesIntegration::write()
 
                 os  << obr_.time().value();
 
-                forAll(integralVals, setI)
+                forAll(integralVals, zoneI)
                 {
-                    os  << ' ' << setw(w) << integralVals[setI];
+                    os  << ' ' << setw(w) << integralVals[zoneI];
 
                     if (log_)
                     {
                         Info<< "faceZonesIntegration output:" << nl
-                            << "    Integration[" << setI << "] "
-                            << integralVals[setI] << endl;
+                            << "    Integration[" << zoneI << "] "
+                            << integralVals[zoneI] << endl;
                     }
                 }
 

src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H

@@ -82,11 +82,11 @@ protected:
             //- Switch to send output to Info as well as to file
             Switch log_;
 
-            //- faceZones to integrate over
-            wordList faceZonesSet_;
+            //- List of face zone names to integrate over
+            wordList zoneNames_;
 
             //- Names of the surface fields
-            wordList fItems_;
+            wordList fieldNames_;
 
 
         //- Current open files

src/sampling/Make/files

@@ -1,16 +1,17 @@
 probes/probes.C
 probes/probesFunctionObject.C
 
-sampledSet/coordSet/coordSet.C
-sampledSet/sampledSet/sampledSet.C
 sampledSet/cloud/cloudSet.C
-sampledSet/face/faceOnlySet.C
+sampledSet/coordSet/coordSet.C
 sampledSet/curve/curveSet.C
-sampledSet/uniform/uniformSet.C
+sampledSet/face/faceOnlySet.C
 sampledSet/midPoint/midPointSet.C
 sampledSet/midPointAndFace/midPointAndFaceSet.C
+sampledSet/sampledSet/sampledSet.C
 sampledSet/sampledSets/sampledSets.C
 sampledSet/sampledSetsFunctionObject/sampledSetsFunctionObject.C
+sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
+sampledSet/uniform/uniformSet.C
 
 setWriters = sampledSet/writers
 
@@ -18,6 +19,7 @@ $(setWriters)/writers.C
 $(setWriters)/gnuplot/gnuplotSetWriterRunTime.C
 $(setWriters)/jplot/jplotSetWriterRunTime.C
 $(setWriters)/raw/rawSetWriterRunTime.C
+$(setWriters)/vtk/vtkSetWriterRunTime.C
 $(setWriters)/xmgrace/xmgraceSetWriterRunTime.C
 
 cuttingPlane/cuttingPlane.C

src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "triSurfaceMeshPointSet.H"
+#include "meshSearch.H"
+#include "DynamicList.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(triSurfaceMeshPointSet, 0);
+    addToRunTimeSelectionTable(sampledSet, triSurfaceMeshPointSet, word);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::triSurfaceMeshPointSet::calcSamples
+(
+    DynamicList<point>& samplingPts,
+    DynamicList<label>& samplingCells,
+    DynamicList<label>& samplingFaces,
+    DynamicList<label>& samplingSegments,
+    DynamicList<scalar>& samplingCurveDist
+) const
+{
+    forAll(sampleCoords_, sampleI)
+    {
+        label cellI = searchEngine().findCell(sampleCoords_[sampleI]);
+
+        if (cellI != -1)
+        {
+            samplingPts.append(sampleCoords_[sampleI]);
+            samplingCells.append(cellI);
+            samplingFaces.append(-1);
+            samplingSegments.append(0);
+            samplingCurveDist.append(1.0 * sampleI);
+        }
+    }
+}
+
+
+void Foam::triSurfaceMeshPointSet::genSamples()
+{
+    // Storage for sample points
+    DynamicList<point> samplingPts;
+    DynamicList<label> samplingCells;
+    DynamicList<label> samplingFaces;
+    DynamicList<label> samplingSegments;
+    DynamicList<scalar> samplingCurveDist;
+
+    calcSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+
+    samplingPts.shrink();
+    samplingCells.shrink();
+    samplingFaces.shrink();
+    samplingSegments.shrink();
+    samplingCurveDist.shrink();
+
+    setSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::triSurfaceMeshPointSet::triSurfaceMeshPointSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    meshSearch& searchEngine,
+    const dictionary& dict
+)
+:
+    sampledSet(name, mesh, searchEngine, dict),
+    surface_(dict.lookup("surface"))
+{
+    // Load surface.
+    if (mesh.time().foundObject<triSurfaceMesh>(surface_))
+    {
+        // Note: should use localPoints() instead of points() but assume
+        // trisurface is compact.
+        sampleCoords_ = mesh.time().lookupObject<triSurfaceMesh>
+        (
+            surface_
+        ).points();
+    }
+    else
+    {
+        sampleCoords_ = triSurfaceMesh
+        (
+            IOobject
+            (
+                surface_,
+                mesh.time().constant(),     // instance
+                "triSurface",               // local
+                mesh.time(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).points();
+    }
+
+    genSamples();
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::triSurfaceMeshPointSet::~triSurfaceMeshPointSet()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::point Foam::triSurfaceMeshPointSet::getRefPoint(const List<point>& pts)
+ const
+{
+    if (pts.size())
+    {
+        // Use first samplePt as starting point
+        return pts[0];
+    }
+    else
+    {
+        return vector::zero;
+    }
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::triSurfaceMeshPointSet
+
+Description
+    sampleSet from all points of a triSurfaceMesh.
+
+SourceFiles
+    triSurfaceMeshPointSet.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triSurfaceMeshPointSet_H
+#define triSurfaceMeshPointSet_H
+
+#include "sampledSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class triSurfaceMeshPointSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class triSurfaceMeshPointSet
+:
+    public sampledSet
+{
+    // Private data
+
+        //- Name of triSurfaceMesh
+        const word surface_;
+
+
+        //- Sampling points
+        List<point> sampleCoords_;
+
+
+    // Private Member Functions
+
+        //- Samples all points in sampleCoords.
+        void calcSamples
+        (
+            DynamicList<point>& samplingPts,
+            DynamicList<label>& samplingCells,
+            DynamicList<label>& samplingFaces,
+            DynamicList<label>& samplingSegments,
+            DynamicList<scalar>& samplingCurveDist
+        ) const;
+
+        //- Uses calcSamples to obtain samples. Copies them into *this.
+        void genSamples();
+
+public:
+
+    //- Runtime type information
+    TypeName("triSurfaceMeshPointSet");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        triSurfaceMeshPointSet
+        (
+            const word& name,
+            const polyMesh& mesh,
+            meshSearch& searchEngine,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~triSurfaceMeshPointSet();
+
+
+    // Member Functions
+
+        //- Get reference point
+        virtual point getRefPoint(const List<point>&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSet/writers/gnuplot/gnuplotSetWriter.C

@@ -69,29 +69,73 @@ void Foam::gnuplotSetWriter<Type>::write
     Ostream& os
 ) const
 {
-    os  << "set term postscript color" << endl
-        << "set output \"" << points.name() << ".ps\"" << endl
+    os  << "set term postscript color" << nl
+        << "set output \"" << points.name() << ".ps\"" << nl
         << "plot";
 
-    bool firstField = true;
-
     forAll(valueSets, i)
     {
-        if (!firstField)
+        if (i != 0)
         {
             os << ',';
         }
-        firstField = false;
 
-        os  << "'-' title \"" << valueSetNames[i] << "\" with lines";
+        os  << " \"-\" title \"" << valueSetNames[i] << "\" with lines";
     }
-    os << endl;
+    os  << nl;
 
 
     forAll(valueSets, i)
     {
-        os << endl;
         writeTable(points, *valueSets[i], os);
+        os  << "e" << nl;
+    }
+}
+
+
+template<class Type>
+void Foam::gnuplotSetWriter<Type>::write
+(
+    const bool writeTracks,
+    const PtrList<coordSet>& trackPoints,
+    const wordList& valueSetNames,
+    const List<List<Field<Type> > >& valueSets,
+    Ostream& os
+) const
+{
+    if (valueSets.size() != valueSetNames.size())
+    {
+        FatalErrorIn("gnuplotSetWriter<Type>::write(..)")
+            << "Number of variables:" << valueSetNames.size() << endl
+            << "Number of valueSets:" << valueSets.size()
+            << exit(FatalError);
+    }
+    if (trackPoints.size() > 0)
+    {
+        os  << "set term postscript color" << nl
+            << "set output \"" << trackPoints[0].name() << ".ps\"" << nl;
+
+        forAll(trackPoints, trackI)
+        {
+            os  << "plot";
+
+            forAll(valueSets, i)
+            {
+                if (i != 0)
+                {
+                    os << ',';
+                }
+
+                os  << " \"-\" title \"" << valueSetNames[i] << "\" with lines";
+            }
+            os << nl;
+
+            forAll(valueSets, i)
+            {
+                writeTable(trackPoints[trackI], valueSets[i][trackI], os);
+                os  << "e" << nl;
+            }
+        }
     }
 }
 

src/sampling/sampledSet/writers/gnuplot/gnuplotSetWriter.H

@@ -76,13 +76,22 @@ public:
             const wordList&
         ) const;
 
-        void write
+        virtual void write
         (
             const coordSet&,
             const wordList&,
             const List<const Field<Type>*>&,
             Ostream&
         ) const;
+
+        virtual void write
+        (
+            const bool writeTracks,
+            const PtrList<coordSet>&,
+            const wordList& valueSetNames,
+            const List<List<Field<Type> > >&,
+            Ostream&
+        ) const;
 };
 
 

src/sampling/sampledSet/writers/jplot/jplotSetWriter.C

@@ -39,7 +39,7 @@ Foam::Ostream& Foam::jplotSetWriter<Type>::writeHeader(Ostream& os) const
     return os
         << "# JPlot input file" << nl
         << "#" << nl << nl
-        << "# Generated by sample on " << clock::date().c_str() << endl;
+        << "# Generated by sample on " << clock::date().c_str() << nl;
 }
 
 
@@ -82,11 +82,11 @@ void Foam::jplotSetWriter<Type>::write
 ) const
 {
     os  << "# JPlot file" << nl
-        << "# column 1: " << points.name() << endl;
+        << "# column 1: " << points.name() << nl;
 
     forAll(valueSets, i)
     {
-        os  << "# column " << i + 2 << ": " << valueSetNames[i] << endl;
+        os  << "# column " << i + 2 << ": " << valueSetNames[i] << nl;
     }
 
     // Collect sets into columns

src/sampling/sampledSet/writers/jplot/jplotSetWriter.H

@@ -80,13 +80,35 @@ public:
             const wordList&
         ) const;
 
-        void write
+        virtual void write
         (
             const coordSet&,
             const wordList&,
             const List<const Field<Type>*>&,
             Ostream&
         ) const;
+
+        virtual void write
+        (
+            const bool writeTracks,
+            const PtrList<coordSet>&,
+            const wordList& valueSetNames,
+            const List<List<Field<Type> > >&,
+            Ostream&
+        ) const
+        {
+            notImplemented
+            (
+                "jplotSetWriter<Type>::write\n"
+                "(\n"
+                "    const bool,\n"
+                "    const PtrList<coordSet>&,\n"
+                "    const wordList&,\n"
+                "    const List<List<Field<Type> > >&,\n"
+                "    Ostream&\n"
+                ") const"
+            );
+        }
 };
 
 

src/sampling/sampledSet/writers/raw/rawSetWriter.C

@@ -79,4 +79,38 @@ void Foam::rawSetWriter<Type>::write
 }
 
 
+template<class Type>
+void Foam::rawSetWriter<Type>::write
+(
+    const bool writeTracks,
+    const PtrList<coordSet>& points,
+    const wordList& valueSetNames,
+    const List<List<Field<Type> > >& valueSets,
+    Ostream& os
+) const
+{
+    if (valueSets.size() != valueSetNames.size())
+    {
+        FatalErrorIn("rawSetWriter<Type>::write(..)")
+            << "Number of variables:" << valueSetNames.size() << endl
+            << "Number of valueSets:" << valueSets.size()
+            << exit(FatalError);
+    }
+
+    List<const List<Type>*> columns(valueSets.size());
+
+    forAll(points, trackI)
+    {
+        // Collect sets into columns
+        forAll(valueSets, i)
+        {
+            columns[i] = &valueSets[i][trackI];
+        }
+
+        writeTable(points[trackI], columns, os);
+        os  << nl << nl;
+    }
+}
+
+
 // ************************************************************************* //

src/sampling/sampledSet/writers/raw/rawSetWriter.H

@@ -83,6 +83,15 @@ public:
             const List<const Field<Type>*>&,
             Ostream&
         ) const;
+
+        virtual void write
+        (
+            const bool writeTracks,
+            const PtrList<coordSet>&,
+            const wordList& valueSetNames,
+            const List<List<Field<Type> > >&,
+            Ostream&
+        ) const;
 };