Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -66,6 +66,7 @@ Description
 #include "Switch.H"
 #include "bound.H"
 #include "dynamicRefineFvMesh.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,6 +84,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -94,7 +97,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -163,10 +165,15 @@ int main(int argc, char *argv[])
 
 
         #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
+        {
             #include "UEqn.H"
 
+
             // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
                 #include "bEqn.H"
                 #include "ftEqn.H"
@@ -181,7 +188,11 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
             }
 
+            if (pimple.turbCorr())
+            {
                 turbulence->correct();
+            }
+        }
 
         runTime.write();
 

applications/test/PatchEdgeFaceWave/Make/files

@@ -0,0 +1,3 @@
+Test-PatchEdgeFaceWave.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-PatchEdgeFaceWave

applications/test/PatchEdgeFaceWave/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

applications/test/PatchEdgeFaceWave/Test-PatchEdgeFaceWave.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Test PatchEdgeFaceWave.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "PatchEdgeFaceWave.H"
+#include "patchEdgeFaceInfo.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validArgs.append("patch");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Get name of patch
+    const word patchName = args[1];
+
+    const polyPatch& patch = patches[patchName];
+
+    // Data on all edges and faces
+    List<patchEdgeFaceInfo> allEdgeInfo(patch.nEdges());
+    List<patchEdgeFaceInfo> allFaceInfo(patch.size());
+
+    // Initial seed
+    DynamicList<label> initialEdges;
+    DynamicList<patchEdgeFaceInfo> initialEdgesInfo;
+
+
+    // Just set an edge on the master
+    if (Pstream::master())
+    {
+        label edgeI = 0;
+        Info<< "Starting walk on edge " << edgeI << endl;
+
+        initialEdges.append(edgeI);
+        const edge& e = patch.edges()[edgeI];
+        initialEdgesInfo.append
+        (
+            patchEdgeFaceInfo
+            (
+                e.centre(patch.localPoints()),
+                0.0
+            )
+        );
+    }
+
+
+    // Walk
+    PatchEdgeFaceWave
+    <
+        primitivePatch,
+        patchEdgeFaceInfo
+    > calc
+    (
+        mesh,
+        patch,
+        initialEdges,
+        initialEdgesInfo,
+        allEdgeInfo,
+        allFaceInfo,
+        returnReduce(patch.nEdges(), sumOp<label>())
+    );
+
+
+    // Extract as patchField
+    volScalarField vsf
+    (
+        IOobject
+        (
+            "patchDist",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("patchDist", dimLength, 0.0)
+    );
+    scalarField pf(vsf.boundaryField()[patch.index()].size());
+    forAll(pf, faceI)
+    {
+        pf[faceI] = Foam::sqrt(allFaceInfo[faceI].distSqr());
+    }
+    vsf.boundaryField()[patch.index()] = pf;
+
+    Info<< "Writing patchDist volScalarField to " << runTime.value()
+        << endl;
+
+    vsf.write();
+
+
+    Info<< "\nEnd\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files

@@ -1,4 +1,3 @@
-patchPointEdgeCirculator.C
 createShellMesh.C
 extrudeToRegionMesh.C
 

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C

@@ -31,130 +31,381 @@ License
 #include "polyAddFace.H"
 #include "polyModifyFace.H"
 #include "polyAddCell.H"
-#include "patchPointEdgeCirculator.H"
+#include "labelPair.H"
+#include "indirectPrimitivePatch.H"
+#include "mapDistribute.H"
+#include "globalMeshData.H"
+#include "PatchTools.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::createShellMesh, 0);
 
+namespace Foam
+{
+template<>
+class minEqOp<labelPair>
+{
+public:
+    void operator()(labelPair& x, const labelPair& y) const
+    {
+        x[0] = min(x[0], y[0]);
+        x[1] = min(x[1], y[1]);
+    }
+};
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::createShellMesh::calcPointRegions
+// Synchronise edges
+void Foam::createShellMesh::syncEdges
 (
-    const primitiveFacePatch& patch,
-    const PackedBoolList& nonManifoldEdge,
-    faceList& pointRegions,
-    labelList& regionPoints
+    const globalMeshData& globalData,
+
+    const labelList& patchEdges,
+    const labelList& coupledEdges,
+    const PackedBoolList& sameEdgeOrientation,
+
+    PackedBoolList& isChangedEdge,
+    DynamicList<label>& changedEdges,
+    labelPairList& allEdgeData
 )
 {
-    pointRegions.setSize(patch.size());
-    forAll(pointRegions, faceI)
-    {
-        const face& f = patch.localFaces()[faceI];
-        pointRegions[faceI].setSize(f.size(), -1);
-    }
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+    const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation();
 
-    label nRegions = 0;
-
-    forAll(pointRegions, faceI)
-    {
-        const face& f = patch.localFaces()[faceI];
-
-        forAll(f, fp)
-        {
-            if (pointRegions[faceI][fp] == -1)
-            {
-                // Found unassigned point. Distribute current region.
-                label pointI = f[fp];
-                label edgeI = patch.faceEdges()[faceI][fp];
-
-                patchPointEdgeCirculator circ
+    // Convert patch-edge data into cpp-edge data
+    labelPairList cppEdgeData
     (
-                    patch,
-                    nonManifoldEdge,
-                    edgeI,
-                    findIndex(patch.edgeFaces()[edgeI], faceI),
-                    pointI
+        map.constructSize(),
+        labelPair(labelMax, labelMax)
     );
 
-                for
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+
+        if (isChangedEdge[patchEdgeI])
+        {
+            const labelPair& data = allEdgeData[patchEdgeI];
+
+            // Patch-edge data needs to be converted into coupled-edge data
+            // (optionally flipped) and consistent in orientation with
+            // other coupled edge (optionally flipped)
+            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
+            {
+                cppEdgeData[coupledEdgeI] = data;
+            }
+            else
+            {
+                cppEdgeData[coupledEdgeI] = labelPair(data[1], data[0]);
+            }
+        }
+    }
+
+    // Synchronise
+    globalData.syncData
     (
-                    patchPointEdgeCirculator iter = circ.begin();
-                    iter != circ.end();
-                    ++iter
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),
+        map,
+        minEqOp<labelPair>()
+    );
+
+    // Back from cpp-edge to patch-edge data
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+
+        if (cppEdgeData[coupledEdgeI] != labelPair(labelMax, labelMax))
+        {
+            const labelPair& data = cppEdgeData[coupledEdgeI];
+
+            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
+            {
+                allEdgeData[patchEdgeI] = data;
+            }
+            else
+            {
+                allEdgeData[patchEdgeI] = labelPair(data[1], data[0]);
+            }
+
+            if (!isChangedEdge[patchEdgeI])
+            {
+                changedEdges.append(patchEdgeI);
+                isChangedEdge[patchEdgeI] = true;
+            }
+        }
+    }
+}
+
+
+void Foam::createShellMesh::calcPointRegions
+(
+    const globalMeshData& globalData,
+    const primitiveFacePatch& patch,
+    const PackedBoolList& nonManifoldEdge,
+
+    faceList& pointGlobalRegions,
+    faceList& pointLocalRegions,
+    labelList& localToGlobalRegion
 )
 {
-                    label face2 = iter.faceID();
+    const indirectPrimitivePatch& cpp = globalData.coupledPatch();
 
-                    if (face2 != -1)
-                    {
-                        const face& f2 = patch.localFaces()[face2];
-                        label fp2 = findIndex(f2, pointI);
-                        label& region = pointRegions[face2][fp2];
-                        if (region != -1)
-                        {
-                            FatalErrorIn
+    // Calculate correspondence between patch and globalData.coupledPatch.
+    labelList patchEdges;
+    labelList coupledEdges;
+    PackedBoolList sameEdgeOrientation;
+    PatchTools::matchEdges
     (
-                                "createShellMesh::calcPointRegions(..)"
-                            )   << "On point " << pointI
-                                << " at:" << patch.localPoints()[pointI]
-                                << " found region:" << region
-                                << abort(FatalError);
-                        }
-                        region = nRegions;
-                    }
-                }
+        cpp,
+        patch,
 
-                nRegions++;
-            }
-        }
-    }
+        coupledEdges,
+        patchEdges,
+        sameEdgeOrientation
+    );
 
 
-    // From region back to originating point (many to one, a point might
-    // have multiple regions though)
-    regionPoints.setSize(nRegions);
-    forAll(pointRegions, faceI)
+    // Initial unique regions
+    // ~~~~~~~~~~~~~~~~~~~~~~
+    // These get merged later on across connected edges.
+
+    // 1. Count
+    label nMaxRegions = 0;
+    forAll(patch.localFaces(), faceI)
     {
         const face& f = patch.localFaces()[faceI];
+        nMaxRegions += f.size();
+    }
 
+    const globalIndex globalRegions(nMaxRegions);
+
+    // 2. Assign unique regions
+    label nRegions = 0;
+
+    pointGlobalRegions.setSize(patch.size());
+    forAll(pointGlobalRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+        labelList& pRegions = pointGlobalRegions[faceI];
+        pRegions.setSize(f.size());
+        forAll(pRegions, fp)
+        {
+            pRegions[fp] = globalRegions.toGlobal(nRegions++);
+        }
+    }
+
+
+    DynamicList<label> changedEdges(patch.nEdges());
+    labelPairList allEdgeData(patch.nEdges(), labelPair(labelMax, labelMax));
+    PackedBoolList isChangedEdge(patch.nEdges());
+
+
+    // Fill initial seed
+    // ~~~~~~~~~~~~~~~~~
+
+    forAll(patch.edgeFaces(), edgeI)
+    {
+        if (!nonManifoldEdge[edgeI])
+        {
+            // Take over value from one face only.
+            const edge& e = patch.edges()[edgeI];
+            label faceI = patch.edgeFaces()[edgeI][0];
+            const face& f = patch.localFaces()[faceI];
+
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = findIndex(f, e[1]);
+            allEdgeData[edgeI] = labelPair
+            (
+                pointGlobalRegions[faceI][fp0],
+                pointGlobalRegions[faceI][fp1]
+            );
+            if (!isChangedEdge[edgeI])
+            {
+                changedEdges.append(edgeI);
+                isChangedEdge[edgeI] = true;
+            }
+        }
+    }
+
+
+    syncEdges
+    (
+        globalData,
+
+        patchEdges,
+        coupledEdges,
+        sameEdgeOrientation,
+
+        isChangedEdge,
+        changedEdges,
+        allEdgeData
+    );
+
+
+    // Edge-Face-Edge walk across patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Across edge minimum regions win
+
+    while (true)
+    {
+        // From edge to face
+        // ~~~~~~~~~~~~~~~~~
+
+        DynamicList<label> changedFaces(patch.size());
+        PackedBoolList isChangedFace(patch.size());
+
+        forAll(changedEdges, changedI)
+        {
+            label edgeI = changedEdges[changedI];
+            const labelPair& edgeData = allEdgeData[edgeI];
+
+            const edge& e = patch.edges()[edgeI];
+            const labelList& eFaces = patch.edgeFaces()[edgeI];
+
+            forAll(eFaces, i)
+            {
+                label faceI = eFaces[i];
+                const face& f = patch.localFaces()[faceI];
+
+                // Combine edgeData with face data
+                label fp0 = findIndex(f, e[0]);
+                if (pointGlobalRegions[faceI][fp0] > edgeData[0])
+                {
+                    pointGlobalRegions[faceI][fp0] = edgeData[0];
+                    if (!isChangedFace[faceI])
+                    {
+                        isChangedFace[faceI] = true;
+                        changedFaces.append(faceI);
+                    }
+                }
+
+                label fp1 = findIndex(f, e[1]);
+                if (pointGlobalRegions[faceI][fp1] > edgeData[1])
+                {
+                    pointGlobalRegions[faceI][fp1] = edgeData[1];
+                    if (!isChangedFace[faceI])
+                    {
+                        isChangedFace[faceI] = true;
+                        changedFaces.append(faceI);
+                    }
+                }
+            }
+        }
+
+
+        label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
+        if (nChangedFaces == 0)
+        {
+            break;
+        }
+
+
+        // From face to edge
+        // ~~~~~~~~~~~~~~~~~
+
+        isChangedEdge = false;
+        changedEdges.clear();
+
+        forAll(changedFaces, i)
+        {
+            label faceI = changedFaces[i];
+            const face& f = patch.localFaces()[faceI];
+            const labelList& fEdges = patch.faceEdges()[faceI];
+
+            forAll(fEdges, fp)
+            {
+                label edgeI = fEdges[fp];
+
+                if (!nonManifoldEdge[edgeI])
+                {
+                    const edge& e = patch.edges()[edgeI];
+                    label fp0 = findIndex(f, e[0]);
+                    label region0 = pointGlobalRegions[faceI][fp0];
+                    label fp1 = findIndex(f, e[1]);
+                    label region1 = pointGlobalRegions[faceI][fp1];
+
+                    if
+                    (
+                        (allEdgeData[edgeI][0] > region0)
+                     || (allEdgeData[edgeI][1] > region1)
+                    )
+                    {
+                        allEdgeData[edgeI] = labelPair(region0, region1);
+                        if (!isChangedEdge[edgeI])
+                        {
+                            changedEdges.append(edgeI);
+                            isChangedEdge[edgeI] = true;
+                        }
+                    }
+                }
+            }
+        }
+
+        syncEdges
+        (
+            globalData,
+
+            patchEdges,
+            coupledEdges,
+            sameEdgeOrientation,
+
+            isChangedEdge,
+            changedEdges,
+            allEdgeData
+        );
+
+
+        label nChangedEdges = returnReduce(changedEdges.size(), sumOp<label>());
+        if (nChangedEdges == 0)
+        {
+            break;
+        }
+    }
+
+
+
+    // Assign local regions
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    // Calculate addressing from global region back to local region
+    pointLocalRegions.setSize(patch.size());
+    Map<label> globalToLocalRegion(globalRegions.localSize()/4);
+    DynamicList<label> dynLocalToGlobalRegion(globalToLocalRegion.size());
+    forAll(patch.localFaces(), faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+        face& pRegions = pointLocalRegions[faceI];
+        pRegions.setSize(f.size());
         forAll(f, fp)
         {
-            regionPoints[pointRegions[faceI][fp]] = f[fp];
-        }
-    }
+            label globalRegionI = pointGlobalRegions[faceI][fp];
 
+            Map<label>::iterator fnd = globalToLocalRegion.find(globalRegionI);
 
-    if (debug)
+            if (fnd != globalToLocalRegion.end())
             {
-        const labelListList& pointFaces = patch.pointFaces();
-        forAll(pointFaces, pointI)
-        {
-            label region = -1;
-            const labelList& pFaces = pointFaces[pointI];
-            forAll(pFaces, i)
-            {
-                label faceI = pFaces[i];
-                const face& f = patch.localFaces()[faceI];
-                label fp = findIndex(f, pointI);
-
-                if (region == -1)
-                {
-                    region = pointRegions[faceI][fp];
+                // Already encountered this global region. Assign same local one
+                pRegions[fp] = fnd();
             }
-                else if (region != pointRegions[faceI][fp])
+            else
             {
-                    Pout<< "Non-manifold point:" << pointI
-                        << " at " << patch.localPoints()[pointI]
-                        << " region:" << region
-                        << " otherRegion:" << pointRegions[faceI][fp]
-                        << endl;
-
-                }
+                // Region not yet seen. Create new one
+                label localRegionI = globalToLocalRegion.size();
+                pRegions[fp] = localRegionI;
+                globalToLocalRegion.insert(globalRegionI, localRegionI);
+                dynLocalToGlobalRegion.append(globalRegionI);
             }
         }
     }
+    localToGlobalRegion.transfer(dynLocalToGlobalRegion);
 }
 
 

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H

@@ -41,6 +41,7 @@ SourceFiles
 
 #include "primitiveFacePatch.H"
 #include "PackedBoolList.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,6 +51,7 @@ namespace Foam
 // Forward declaration of classes
 class mapPolyMesh;
 class polyTopoChange;
+class globalMeshData;
 
 /*---------------------------------------------------------------------------*\
                            Class createShellMesh Declaration
@@ -80,6 +82,18 @@ class createShellMesh
 
     // Private Member Functions
 
+        static void syncEdges
+        (
+            const globalMeshData&,
+            const labelList&,
+            const labelList&,
+            const PackedBoolList& sameEdgeOrientation,
+            PackedBoolList& isChangedEdge,
+            DynamicList<label>& changedEdges,
+            labelPairList& allEdgeData
+        );
+
+
         //- Disallow default bitwise copy construct
         createShellMesh(const createShellMesh&);
 
@@ -142,13 +156,22 @@ public:
 
         // Other
 
-            //- Helper: calculate point regions
+            //- Helper: calculate point regions. The point region is the
+            //  same on all faces connected to a point if they can be
+            //  reached through a face-edge-face walk without crossing
+            //  the nonManifoldEdge.
+            //  pointGlobalRegions : non-compact. Guaranteed to be the same
+            //  across processors.
+            //  pointLocalRegions : compact.
+            //  localToGlobalRegion : local to global region.
             static void calcPointRegions
             (
+                const globalMeshData& globalData,
                 const primitiveFacePatch& patch,
                 const PackedBoolList& nonManifoldEdge,
-                faceList& pointRegions,
-                labelList& regionPoints
+                faceList& pointGlobalRegions,
+                faceList& pointLocalRegions,
+                labelList& localToGlobalRegion
             );
 
             //- Play commands into polyTopoChange to create layer mesh.

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H

@@ -1,205 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::patchPointEdgeCirculator
-
-Description
-    Walks from starting edge/face around point on patch.
-
-    -# Use in-place: \n
-        \code
-            patchPointEdgeCirculator circ(..);
-
-            // Walk
-            do
-            {
-                Info<< "edge:" << circ.edgeID() << endl;
-                ++circ;
-            }
-            while (circ != circ.end());
-        \endcode
-
-    -# Use like STL iterator: \n
-        \code
-            patchPointEdgeCirculator circ(..);
-            for
-            (
-                patchPointEdgeCirculator iter = circ.begin();
-                iter != circ.end();
-                ++iter
-            )
-            {
-                Info<< "edge:" << iter.edgeID() << endl;
-            }
-        \endcode
-
-
-SourceFiles
-    patchPointEdgeCirculator.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef patchPointEdgeCirculator_H
-#define patchPointEdgeCirculator_H
-
-#include "face.H"
-#include "ListOps.H"
-#include "primitiveFacePatch.H"
-#include "PackedBoolList.H"
-#include "InfoProxy.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class patchPointEdgeCirculator Declaration
-\*---------------------------------------------------------------------------*/
-
-class patchPointEdgeCirculator
-{
-    // Static data members
-
-        //- end iterator
-        static const patchPointEdgeCirculator endConstIter;
-
-
-    // Private data
-
-        //- Patch
-        const primitiveFacePatch& patch_;
-
-        const PackedBoolList& nonManifoldEdge_;
-
-        //- Current edge
-        label edgeID_;
-
-        //- Current direction (face, expressed as index into edgeFaces()[edgeID]
-        label index_;
-
-        //- Point
-        label pointID_;
-
-        //- Starting edge so we know when to stop.
-        label startEdgeID_;
-
-
-    // Private Member Functions
-
-        //- Set to end() iterator
-        inline void setEnd();
-
-        //- Set edgeID_ to be the other edge on the face that uses the
-        //  point.
-        inline void otherEdge(const label cellI);
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        inline patchPointEdgeCirculator
-        (
-            const primitiveFacePatch&,
-            const PackedBoolList& nonManifoldEdge,
-            const label edgeID,
-            const label index,
-            const label pointID
-        );
-
-        //- Construct as copy
-        inline patchPointEdgeCirculator(const patchPointEdgeCirculator&);
-
-
-    // Member Functions
-
-        inline label edgeID() const;
-
-        inline label index() const;
-
-        inline label pointID() const;
-
-        //- Helper: get face according to the index.
-        //  Returns -1 if on end.
-        inline label faceID() const;
-
-        //- Walk back until non-manifold edge (if any) or minimum edge.
-        inline void setCanonical();
-
-
-    // Member Operators
-
-        inline void operator=(const patchPointEdgeCirculator& iter);
-
-        inline bool operator==(const patchPointEdgeCirculator& iter) const;
-
-        inline bool operator!=(const patchPointEdgeCirculator& iter) const;
-
-        //- Step to next face.
-        inline patchPointEdgeCirculator& operator++();
-
-        //- iterator set to the beginning face. For internal edges this is
-        //  the current face. For boundary edges this is the first boundary face
-        //  reached from walking back (i.e. in opposite direction to ++)
-        inline patchPointEdgeCirculator begin() const;
-        inline patchPointEdgeCirculator cbegin() const;
-
-        //- iterator set to beyond the end of the walk.
-        inline const patchPointEdgeCirculator& end() const;
-        inline const patchPointEdgeCirculator& cend() const;
-
-    // Info
-
-        //- Return info proxy.
-        //  Used to print token information to a stream
-        InfoProxy<patchPointEdgeCirculator> info() const
-        {
-            return *this;
-        }
-
-        friend Ostream& operator<<
-        (
-            Ostream&,
-            const InfoProxy<patchPointEdgeCirculator>&
-        );
-};
-
-Ostream& operator<<(Ostream&, const InfoProxy<patchPointEdgeCirculator>&);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "patchPointEdgeCirculatorI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H

@@ -1,308 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::patchPointEdgeCirculator::setEnd()
-{
-    edgeID_ = -1;
-    pointID_ = -1;
-}
-
-
-// Cross face to other edge on point
-void Foam::patchPointEdgeCirculator::otherEdge(const label faceI)
-{
-    const labelList& fEdges = patch_.faceEdges()[faceI];
-    const face& f = patch_.localFaces()[faceI];
-    label fp = findIndex(f, pointID_);
-
-    if (fEdges[fp] == edgeID_)
-    {
-        edgeID_ = fEdges[f.rcIndex(fp)];
-    }
-    else
-    {
-        // Check for now
-        if (fEdges[f.rcIndex(fp)] != edgeID_)
-        {
-            FatalErrorIn("patchPointEdgeCirculator::otherEdge(const label)")
-                << "face:" << faceI
-                << " verts:" << f
-                << " edges:" << fEdges
-                << " looking for other edge than " << edgeID_
-                << abort(FatalError);
-        }
-        edgeID_ = fEdges[fp];
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-//- Construct from components
-Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
-(
-    const primitiveFacePatch& patch,
-    const PackedBoolList& nonManifoldEdge,
-    const label edgeID,
-    const label index,
-
-    const label pointID
-)
-:
-    patch_(patch),
-    nonManifoldEdge_(nonManifoldEdge),
-    edgeID_(edgeID),
-    index_(index),
-    pointID_(pointID),
-    startEdgeID_(edgeID_)
-{
-    if (edgeID_ != -1)
-    {
-        const edge& e = patch_.edges()[edgeID_];
-
-        if (pointID_ != e[0] && pointID_ != e[1])
-        {
-            FatalErrorIn
-            (
-                "patchPointEdgeCirculator::patchPointEdgeCirculator(..)"
-            )   << "edge " << edgeID_ << " verts " << e
-                << " does not contain point " << pointID_ << abort(FatalError);
-        }
-    }
-}
-
-
-//- Construct copy
-Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
-(
-    const patchPointEdgeCirculator& circ
-)
-:
-    patch_(circ.patch_),
-    nonManifoldEdge_(circ.nonManifoldEdge_),
-    edgeID_(circ.edgeID_),
-    index_(circ.index_),
-    pointID_(circ.pointID_),
-    startEdgeID_(circ.startEdgeID_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::label Foam::patchPointEdgeCirculator::edgeID() const
-{
-    return edgeID_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::index() const
-{
-    return index_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::pointID() const
-{
-    return pointID_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::faceID() const
-{
-    if (edgeID_ != -1 && index_ != -1)
-    {
-        return patch_.edgeFaces()[edgeID_][index_];
-    }
-    else
-    {
-        return -1;
-    }
-}
-
-
-void Foam::patchPointEdgeCirculator::operator=
-(
-    const patchPointEdgeCirculator& circ
-)
-{
-    edgeID_ = circ.edgeID_;
-    index_ = circ.index_;
-    pointID_ = circ.pointID_;
-    startEdgeID_ = circ.startEdgeID_;
-}
-
-
-bool Foam::patchPointEdgeCirculator::operator==
-(
-    const patchPointEdgeCirculator& circ
-) const
-{
-    // Do just enough to have setEnd() produce an iterator equal to end().
-    // Could include the direction(i.e. index_) to make sure that two
-    // circulators around same point but in different direction are not equal.
-    return edgeID_ == circ.edgeID_ && pointID_ == circ.pointID_;
-}
-
-
-bool Foam::patchPointEdgeCirculator::operator!=
-(
-    const patchPointEdgeCirculator& circ
-) const
-{
-    return !(*this == circ);
-}
-
-
-void Foam::patchPointEdgeCirculator::setCanonical()
-{
-    if (edgeID_ == -1)
-    {
-        FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
-            << "Already reached end(). Cannot walk any further."
-            << abort(FatalError);
-    }
-
-    label minEdgeID = edgeID_;
-    label minIndex = index_;
-
-    while (true)
-    {
-        if (nonManifoldEdge_[edgeID_])
-        {
-            break;
-        }
-
-        // Step back
-        const labelList& eFaces = patch_.edgeFaces()[edgeID_];
-
-        if (eFaces.size() != 2)
-        {
-            FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
-                << "problem eFaces:" << eFaces << abort(FatalError);
-        }
-
-        label faceI = (index_ == 0 ? eFaces[1] : eFaces[0]);
-
-        // Step to other edge on face
-        otherEdge(faceI);
-
-        // Update index
-        index_ = findIndex(patch_.edgeFaces()[edgeID_], faceI);
-
-        if (edgeID_ < minEdgeID)
-        {
-            minEdgeID = edgeID_;
-            minIndex = index_;
-        }
-
-        if (edgeID_ == startEdgeID_)
-        {
-            edgeID_ = minEdgeID;
-            index_ = minIndex;
-            break;
-        }
-    }
-
-    startEdgeID_ = edgeID_;
-}
-
-
-//- Step to next edge.
-Foam::patchPointEdgeCirculator&
-Foam::patchPointEdgeCirculator::operator++()
-{
-    if (index_ == -1)
-    {
-        setEnd();
-    }
-    else
-    {
-        // Step to other edge on face
-        label faceI = patch_.edgeFaces()[edgeID_][index_];
-        otherEdge(faceI);
-
-        if (edgeID_ == startEdgeID_)
-        {
-            setEnd();
-        }
-        else if (nonManifoldEdge_[edgeID_])
-        {
-            // Reached non-manifold edge. Cannot walk further.
-            // Mark so it gets set to end next time.
-            index_ = -1;
-        }
-        else
-        {
-            const labelList& eFaces = patch_.edgeFaces()[edgeID_];
-
-            if (eFaces.size() != 2)
-            {
-                FatalErrorIn("patchPointEdgeCirculator:::operator++()")
-                    << "problem eFaces:" << eFaces << abort(FatalError);
-            }
-            // Point to the face that is not faceI
-            index_ = (eFaces[0] != faceI ? 0 : 1);
-        }
-    }
-
-    return *this;
-}
-
-
-Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::begin() const
-{
-    patchPointEdgeCirculator iter(*this);
-    iter.setCanonical();
-
-    return iter;
-}
-
-
-Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::cbegin() const
-{
-    patchPointEdgeCirculator iter(*this);
-    iter.setCanonical();
-
-    return iter;
-}
-
-
-const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::end()
- const
-{
-    return endConstIter;
-}
-
-const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::cend()
- const
-{
-    return endConstIter;
-}
-
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -36,7 +36,7 @@ Foam::label Foam::checkTopology
         label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
 
         // These are actually warnings, not errors.
-        if (nEmpty % nTotCells)
+        if (nTotCells && (nEmpty % nTotCells))
         {
             Info<< " ***Total number of faces on empty patches"
                 << " is not divisible by the number of cells in the mesh."
@@ -255,7 +255,7 @@ Foam::label Foam::checkTopology
     {
         regionSplit rs(mesh);
 
-        if (rs.nRegions() == 1)
+        if (rs.nRegions() <= 1)
         {
             Info<< "    Number of regions: " << rs.nRegions() << " (OK)."
                 << endl;

applications/utilities/mesh/manipulation/topoSet/topoSetDict

@@ -319,7 +319,7 @@ FoamFile
 //    source setToFaceZone;
 //    sourceInfo
 //    {
-//        set f0;           // name of faceSet
+//        faceSet f0;       // name of faceSet
 //    }
 //
 //    // Select based on faceSet, using cellSet to determine orientation

etc/config/settings.csh

@@ -203,10 +203,6 @@ case ThirdParty:
     switch ("$WM_COMPILER")
     case Gcc:
     case Gcc++0x:
-        set gcc_version=gcc-4.4.3
-        set gmp_version=gmp-5.0.1
-        set mpfr_version=mpfr-2.4.2
-        breaksw
     case Gcc46:
     case Gcc46++0x:
         set gcc_version=gcc-4.6.1
@@ -487,9 +483,31 @@ case QSMPI:
     breaksw
 
 case SGIMPI:
-    setenv FOAM_MPI ${MPI_ROOT##*/}
+    if ( ! $?MPI_ROOT) setenv MPI_ROOT /dummy
+
+    if ( ! -d "$MPI_ROOT" ) then
+        echo "Warning in $WM_PROJECT_DIR/etc/config/settings.csh:"
+        echo "    MPI_ROOT not a valid mpt installation directory."
+        echo "    Please set MPI_ROOT to the mpt installation directory."
+        echo "    (usually done by loading the mpt module)"
+        echo "    MPI_ROOT currently set to '$MPI_ROOT'"
+    endif
+
+    if ( "${MPI_ROOT:h}/" == $MPI_ROOT ) then
+        setenv MPI_ROOT ${MPI_ROOT:h}
+    endif
+
+    setenv FOAM_MPI ${MPI_ROOT:t}
     setenv MPI_ARCH_PATH $MPI_ROOT
 
+
+    if ($?FOAM_VERBOSE && $?prompt) then
+        echo "Using SGI MPT:"
+        echo "    MPI_ROOT : $MPI_ROOT"
+        echo "    FOAM_MPI : $FOAM_MPI"
+    endif
+
+
     _foamAddPath    $MPI_ARCH_PATH/bin
     _foamAddLib     $MPI_ARCH_PATH/lib
     breaksw

etc/config/settings.sh

@@ -222,12 +222,7 @@ fi
 case "${foamCompiler}" in
 OpenFOAM | ThirdParty)
     case "$WM_COMPILER" in
-    Gcc | Gcc++0x)
-        gcc_version=gcc-4.4.3
-        gmp_version=gmp-5.0.1
-        mpfr_version=mpfr-2.4.2
-        ;;
-    Gcc46 | Gcc46++0x)
+    Gcc | Gcc++0x | Gcc46 | Gcc46++0x)
         gcc_version=gcc-4.6.1
         gmp_version=gmp-5.0.2
         mpfr_version=mpfr-3.0.1
@@ -511,9 +506,30 @@ QSMPI)
     ;;
 
 SGIMPI)
+    lastCharID=$(( ${#MPI_ROOT} - 1 ))
+    if [ "${MPI_ROOT:$lastCharID:1}" == '/' ]
+    then
+        MPI_ROOT=${MPI_ROOT:0:$lastCharID}
+    fi
+
     export FOAM_MPI=${MPI_ROOT##*/}
     export MPI_ARCH_PATH=$MPI_ROOT
 
+    if [ ! -d "$MPI_ROOT" -o -z "$MPI_ARCH_PATH" ]
+    then
+        echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:" 1>&2
+        echo "    MPI_ROOT not a valid mpt installation directory or ending in a '/'." 1>&2
+        echo "    Please set MPI_ROOT to the mpt installation directory." 1>&2
+        echo "    MPI_ROOT currently set to '$MPI_ROOT'" 1>&2
+    fi
+
+    if [ "$FOAM_VERBOSE" -a "$PS1" ]
+    then
+        echo "Using SGI MPT:"
+        echo "    MPI_ROOT : $MPI_ROOT"
+        echo "    FOAM_MPI : $FOAM_MPI"
+    fi
+
     _foamAddPath    $MPI_ARCH_PATH/bin
     _foamAddLib     $MPI_ARCH_PATH/lib
     ;;

etc/controlDict

@@ -58,6 +58,11 @@ OptimisationSwitches
     commsType       nonBlocking; //scheduled; //blocking;
     floatTransfer   0;
     nProcsSimpleSum 0;
+
+    // Force dumping (at next timestep) upon signal (-1 to disable)
+    writeNowSignal              -1; //10;
+    // Force dumping (at next timestep) upon signal (-1 to disable) and exit
+    stopAtWriteNowSignal        -1;
 }
 
 

src/OSspecific/POSIX/Make/files

@@ -2,6 +2,8 @@ signals/sigFpe.C
 signals/sigSegv.C
 signals/sigInt.C
 signals/sigQuit.C
+signals/sigStopAtWriteNow.C
+signals/sigWriteNow.C
 regExp.C
 timer.C
 fileStat.C

src/OSspecific/POSIX/POSIX.C

@@ -123,7 +123,7 @@ bool Foam::setEnv
 }
 
 
-Foam::word Foam::hostName(bool full)
+Foam::string Foam::hostName(bool full)
 {
     char buf[128];
     ::gethostname(buf, sizeof(buf));
@@ -142,7 +142,7 @@ Foam::word Foam::hostName(bool full)
 }
 
 
-Foam::word Foam::domainName()
+Foam::string Foam::domainName()
 {
     char buf[128];
     ::gethostname(buf, sizeof(buf));
@@ -159,11 +159,11 @@ Foam::word Foam::domainName()
         }
     }
 
-    return word::null;
+    return string::null;
 }
 
 
-Foam::word Foam::userName()
+Foam::string Foam::userName()
 {
     struct passwd* pw = ::getpwuid(::getuid());
 
@@ -173,7 +173,7 @@ Foam::word Foam::userName()
     }
     else
     {
-        return word::null;
+        return string::null;
     }
 }
 
@@ -209,7 +209,7 @@ Foam::fileName Foam::home()
 }
 
 
-Foam::fileName Foam::home(const word& userName)
+Foam::fileName Foam::home(const string& userName)
 {
     struct passwd* pw;
 
@@ -1069,7 +1069,7 @@ void Foam::fdClose(const int fd)
 
 bool Foam::ping
 (
-    const word& destName,
+    const string& destName,
     const label destPort,
     const label timeOut
 )
@@ -1083,7 +1083,7 @@ bool Foam::ping
     {
         FatalErrorIn
         (
-            "Foam::ping(const word&, ...)"
+            "Foam::ping(const string&, ...)"
         )   << "gethostbyname error " << h_errno << " for host " << destName
             << abort(FatalError);
     }
@@ -1097,7 +1097,7 @@ bool Foam::ping
     {
         FatalErrorIn
         (
-            "Foam::ping(const word&, const label)"
+            "Foam::ping(const string&, const label)"
         )   << "socket error"
             << abort(FatalError);
     }
@@ -1149,7 +1149,7 @@ bool Foam::ping
 }
 
 
-bool Foam::ping(const word& hostname, const label timeOut)
+bool Foam::ping(const string& hostname, const label timeOut)
 {
     return ping(hostname, 222, timeOut) || ping(hostname, 22, timeOut);
 }

src/OSspecific/POSIX/signals/sigStopAtWriteNow.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sigStopAtWriteNow.H"
+#include "error.H"
+#include "JobInfo.H"
+#include "IOstreams.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// Signal number to catch
+int Foam::sigStopAtWriteNow::signal_
+(
+    debug::optimisationSwitch("stopAtWriteNowSignal", -1)
+);
+
+static Foam::Time const* runTimePtr_ = NULL;
+
+
+struct sigaction Foam::sigStopAtWriteNow::oldAction_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sigStopAtWriteNow::sigHandler(int)
+{
+    // Reset old handling
+    if (sigaction(signal_, &oldAction_, NULL) < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::sigStopAtWriteNow::sigHandler(int)"
+        )   << "Cannot reset " << signal_ << " trapping"
+            << abort(FatalError);
+    }
+
+    // Update jobInfo file
+    jobInfo.signalEnd();
+
+    Info<< "sigStopAtWriteNow :"
+        << " setting up write and stop at end of the next iteration"
+        << nl << endl;
+    runTimePtr_->stopAt(Time::saWriteNow);
+
+    //// Throw signal (to old handler)
+    //raise(signal_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sigStopAtWriteNow::sigStopAtWriteNow(){}
+
+
+Foam::sigStopAtWriteNow::sigStopAtWriteNow
+(
+    const bool verbose,
+    const Time& runTime
+)
+{
+    if (signal_ > 0)
+    {
+        // Check that the signal is different from the writeNowSignal
+        if (sigWriteNow::signal_ == signal_)
+        {
+            FatalErrorIn
+            (
+                "Foam::sigStopAtWriteNow::sigStopAtWriteNow"
+                "(const bool, const Time&)"
+            )   << "stopAtWriteNowSignal : " << signal_
+                << " cannot be the same as the writeNowSignal."
+                << " Please change this in the controlDict ("
+                << findEtcFile("controlDict", false) << ")."
+                << exit(FatalError);
+        }
+
+
+        // Store runTime
+        runTimePtr_ = &runTime;
+
+        struct sigaction newAction;
+        newAction.sa_handler = sigHandler;
+        newAction.sa_flags = SA_NODEFER;
+        sigemptyset(&newAction.sa_mask);
+        if (sigaction(signal_, &newAction, &oldAction_) < 0)
+        {
+            FatalErrorIn
+            (
+                "Foam::sigStopAtWriteNow::sigStopAtWriteNow"
+                "(const bool, const Time&)"
+            )   << "Cannot set " << signal_ << " trapping"
+                << abort(FatalError);
+        }
+
+        if (verbose)
+        {
+            Info<< "sigStopAtWriteNow :"
+                << " Enabling writing and stopping upon signal " << signal_
+                << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sigStopAtWriteNow::~sigStopAtWriteNow()
+{
+    // Reset old handling
+    if (signal_ > 0)
+    {
+        if (sigaction(signal_, &oldAction_, NULL) < 0)
+        {
+            FatalErrorIn
+            (
+                "Foam::sigStopAtWriteNow::~sigStopAtWriteNow()"
+            )   << "Cannot reset " << signal_ << " trapping"
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sigStopAtWriteNow::active() const
+{
+    return signal_ > 0;
+}
+
+
+// ************************************************************************* //

src/OSspecific/POSIX/signals/sigStopAtWriteNow.H

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sigStopAtWriteNow
+
+Description
+    Signal handler for interupt defined by
+        OptimisationSwitches::stopAtWriteNowSignal
+
+    Write and stop the job.
+
+SourceFiles
+    sigStopAtWriteNow.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sigStopAtWriteNow_H
+#define sigStopAtWriteNow_H
+
+#include <signal.h>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class Time;
+
+/*---------------------------------------------------------------------------*\
+                           Class sigStopAtWriteNow Declaration
+\*---------------------------------------------------------------------------*/
+
+class sigStopAtWriteNow
+{
+    // Private data
+
+        //- number of signal to use
+        static int signal_;
+
+        //- Saved old signal trapping setting
+        static struct sigaction oldAction_;
+
+    // Private Member Functions
+
+        static void sigHandler(int);
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        sigStopAtWriteNow();
+
+        //- Construct from components
+        sigStopAtWriteNow(const bool verbose, const Time& runTime);
+
+
+    //- Destructor
+    ~sigStopAtWriteNow();
+
+
+    // Member functions
+
+        //- Is active?
+        bool active() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OSspecific/POSIX/signals/sigWriteNow.C

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sigWriteNow.H"
+#include "error.H"
+#include "JobInfo.H"
+#include "IOstreams.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// Signal number to catch
+int Foam::sigWriteNow::signal_
+(
+    debug::optimisationSwitch("writeNowSignal", -1)
+);
+
+static Foam::Time* runTimePtr_ = NULL;
+
+
+struct sigaction Foam::sigWriteNow::oldAction_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sigWriteNow::sigHandler(int)
+{
+    Info<< "sigWriteNow :"
+        << " setting up write at end of the next iteration" << nl << endl;
+    runTimePtr_->writeOnce();
+
+    //// Throw signal (to old handler)
+    //raise(signal_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sigWriteNow::sigWriteNow()
+{}
+
+
+Foam::sigWriteNow::sigWriteNow(const bool verbose, Time& runTime)
+{
+    if (signal_ >= 0)
+    {
+        // Store runTime
+        runTimePtr_ = &runTime;
+
+        struct sigaction newAction;
+        newAction.sa_handler = sigHandler;
+        newAction.sa_flags = SA_NODEFER;
+        sigemptyset(&newAction.sa_mask);
+        if (sigaction(signal_, &newAction, &oldAction_) < 0)
+        {
+            FatalErrorIn
+            (
+                "Foam::sigWriteNow::sigWriteNow(const bool, const Time&)"
+            )   << "Cannot set " << signal_ << " trapping"
+                << abort(FatalError);
+        }
+
+        if (verbose)
+        {
+            Info<< "sigWriteNow :"
+                << " Enabling writing upon signal " << signal_
+                << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sigWriteNow::~sigWriteNow()
+{
+    // Reset old handling
+    if (signal_ > 0)
+    {
+        if (sigaction(signal_, &oldAction_, NULL) < 0)
+        {
+            FatalErrorIn
+            (
+                "Foam::sigWriteNow::~sigWriteNow()"
+            )   << "Cannot reset " << signal_ << " trapping"
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sigWriteNow::active() const
+{
+    return signal_ > 0;
+}
+
+
+// ************************************************************************* //

src/OSspecific/POSIX/signals/sigWriteNow.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sigWriteNow
+
+Description
+    Signal handler for interupt defined by OptimisationSwitches::writeNowSignal
+
+    Write once and continue.
+
+SourceFiles
+    sigWriteNow.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sigWriteNow_H
+#define sigWriteNow_H
+
+#include <signal.h>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class Time;
+
+/*---------------------------------------------------------------------------*\
+                           Class sigWriteNow Declaration
+\*---------------------------------------------------------------------------*/
+
+class sigWriteNow
+{
+    // Private data
+
+        //- number of signal to use
+        static int signal_;
+
+        //- Saved old signal trapping setting
+        static struct sigaction oldAction_;
+
+    // Private Member Functions
+
+        static void sigHandler(int);
+
+
+public:
+
+    friend class sigStopAtWriteNow;
+
+    // Constructors
+
+        //- Construct null
+        sigWriteNow();
+
+        //- Construct from components
+        sigWriteNow(const bool verbose, Time& runTime);
+
+
+    //- Destructor
+    ~sigWriteNow();
+
+
+    // Member functions
+
+        //- Is active?
+        bool active() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/db/Time/Time.C

@@ -81,10 +81,17 @@ void Foam::Time::adjustDeltaT()
 {
     if (writeControl_ == wcAdjustableRunTime)
     {
+        scalar interval = writeInterval_;
+        if (secondaryWriteControl_ == wcAdjustableRunTime)
+        {
+            interval = min(interval, secondaryWriteInterval_);
+        }
+
+
         scalar timeToNextWrite = max
         (
             0.0,
-            (outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_)
+            (outputTimeIndex_ + 1)*interval - (value() - startTime_)
         );
 
         scalar nSteps = timeToNextWrite/deltaT_ - SMALL;
@@ -252,8 +259,13 @@ Foam::Time::Time
     stopAt_(saEndTime),
     writeControl_(wcTimeStep),
     writeInterval_(GREAT),
+    secondaryWriteControl_(wcTimeStep),
+    secondaryWriteInterval_(labelMax),
     purgeWrite_(0),
+    writeOnce_(false),
     subCycling_(false),
+    sigWriteNow_(true, *this),
+    sigStopAtWriteNow_(true, *this),
 
     writeFormat_(IOstream::ASCII),
     writeVersion_(IOstream::currentVersion),
@@ -339,8 +351,13 @@ Foam::Time::Time
     stopAt_(saEndTime),
     writeControl_(wcTimeStep),
     writeInterval_(GREAT),
+    secondaryWriteControl_(wcTimeStep),
+    secondaryWriteInterval_(labelMax),
     purgeWrite_(0),
+    writeOnce_(false),
     subCycling_(false),
+    sigWriteNow_(true, *this),
+    sigStopAtWriteNow_(true, *this),
 
     writeFormat_(IOstream::ASCII),
     writeVersion_(IOstream::currentVersion),
@@ -429,8 +446,13 @@ Foam::Time::Time
     stopAt_(saEndTime),
     writeControl_(wcTimeStep),
     writeInterval_(GREAT),
+    secondaryWriteControl_(wcTimeStep),
+    secondaryWriteInterval_(labelMax),
     purgeWrite_(0),
+    writeOnce_(false),
     subCycling_(false),
+    sigWriteNow_(true, *this),
+    sigStopAtWriteNow_(true, *this),
 
     writeFormat_(IOstream::ASCII),
     writeVersion_(IOstream::currentVersion),
@@ -521,7 +543,10 @@ Foam::Time::Time
     stopAt_(saEndTime),
     writeControl_(wcTimeStep),
     writeInterval_(GREAT),
+    secondaryWriteControl_(wcTimeStep),
+    secondaryWriteInterval_(labelMax),
     purgeWrite_(0),
+    writeOnce_(false),
     subCycling_(false),
 
     writeFormat_(IOstream::ASCII),
@@ -944,6 +969,35 @@ Foam::Time& Foam::Time::operator++()
 
     if (!subCycling_)
     {
+        if (sigStopAtWriteNow_.active() || sigWriteNow_.active())
+        {
+            // A signal might have been sent on one processor only
+            // Reduce so all decide the same.
+
+            label flag = 0;
+            if (sigStopAtWriteNow_.active() && stopAt_ == saWriteNow)
+            {
+                flag += 1;
+            }
+            if (sigWriteNow_.active() && writeOnce_)
+            {
+                flag += 2;
+            }
+            reduce(flag, maxOp<label>());
+
+            if (flag & 1)
+            {
+                stopAt_ = saWriteNow;
+            }
+            if (flag & 2)
+            {
+                writeOnce_ = true;
+            }
+        }
+
+
+        outputTime_ = false;
+
         switch (writeControl_)
         {
             case wcTimeStep:
@@ -964,10 +1018,6 @@ Foam::Time& Foam::Time::operator++()
                     outputTime_ = true;
                     outputTimeIndex_ = outputIndex;
                 }
-                else
-                {
-                    outputTime_ = false;
-                }
             }
             break;
 
@@ -983,10 +1033,6 @@ Foam::Time& Foam::Time::operator++()
                     outputTime_ = true;
                     outputTimeIndex_ = outputIndex;
                 }
-                else
-                {
-                    outputTime_ = false;
-                }
             }
             break;
 
@@ -1002,14 +1048,69 @@ Foam::Time& Foam::Time::operator++()
                     outputTime_ = true;
                     outputTimeIndex_ = outputIndex;
                 }
-                else
+            }
+            break;
+        }
+
+
+        // Adapt for secondaryWrite controls
+        switch (secondaryWriteControl_)
         {
-                    outputTime_ = false;
+            case wcTimeStep:
+                outputTime_ =
+                    outputTime_
+                || !(timeIndex_ % label(secondaryWriteInterval_));
+            break;
+
+            case wcRunTime:
+            case wcAdjustableRunTime:
+            {
+                label outputIndex = label
+                (
+                    ((value() - startTime_) + 0.5*deltaT_)
+                  / secondaryWriteInterval_
+                );
+
+                if (outputIndex > outputTimeIndex_)
+                {
+                    outputTime_ = true;
+                    outputTimeIndex_ = outputIndex;
+                }
+            }
+            break;
+
+            case wcCpuTime:
+            {
+                label outputIndex = label
+                (
+                    returnReduce(elapsedCpuTime(), maxOp<double>())
+                  / secondaryWriteInterval_
+                );
+                if (outputIndex > outputTimeIndex_)
+                {
+                    outputTime_ = true;
+                    outputTimeIndex_ = outputIndex;
+                }
+            }
+            break;
+
+            case wcClockTime:
+            {
+                label outputIndex = label
+                (
+                    returnReduce(label(elapsedClockTime()), maxOp<label>())
+                  / secondaryWriteInterval_
+                );
+                if (outputIndex > outputTimeIndex_)
+                {
+                    outputTime_ = true;
+                    outputTimeIndex_ = outputIndex;
                 }
             }
             break;
         }
 
+
         // see if endTime needs adjustment to stop at the next run()/end() check
         if (!end())
         {
@@ -1027,6 +1128,14 @@ Foam::Time& Foam::Time::operator++()
                 endTime_ = value();
             }
         }
+
+        // Override outputTime if one-shot writing
+        if (writeOnce_)
+        {
+            outputTime_ = true;
+            writeOnce_ = false;
+        }
+
     }
 
     return *this;

src/OpenFOAM/db/Time/Time.H

@@ -52,6 +52,8 @@ SourceFiles
 #include "dlLibraryTable.H"
 #include "functionObjectList.H"
 #include "fileMonitor.H"
+#include "sigWriteNow.H"
+#include "sigStopAtWriteNow.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -130,15 +132,35 @@ protected:
 
         scalar writeInterval_;
 
+        // Additional writing
+
+            writeControls secondaryWriteControl_;
+
+            scalar secondaryWriteInterval_;
+
+
         label  purgeWrite_;
         mutable FIFOStack<word> previousOutputTimes_;
 
+        // One-shot writing
+        bool writeOnce_;
+
         //- Is the time currently being sub-cycled?
         bool subCycling_;
 
         //- If time is being sub-cycled this is the previous TimeState
         autoPtr<TimeState> prevTimeState_;
 
+
+        // Signal handlers for secondary writing
+
+            //- Enable one-shot writing upon signal
+            sigWriteNow sigWriteNow_;
+
+            //- Enable write and clean exit upon signal
+            sigStopAtWriteNow sigStopAtWriteNow_;
+
+
         //- Time directory name format
         static fmtflags format_;
 
@@ -357,12 +379,16 @@ public:
                 IOstream::compressionType
             ) const;
 
-            //- Write the objects now and continue the run
+            //- Write the objects now (not at end of iteration) and continue
+            //  the run
             bool writeNow();
 
-            //- Write the objects now and end the run
+            //- Write the objects now (not at end of iteration) and end the run
             bool writeAndEnd();
 
+            //- Write the objects once (one shot) and continue the run
+            void writeOnce();
+
 
         // Access
 

src/OpenFOAM/db/Time/TimeIO.C

@@ -58,6 +58,45 @@ void Foam::Time::readDict()
         controlDict_.lookup("writeFrequency") >> writeInterval_;
     }
 
+
+    // Additional writing
+    if (controlDict_.found("secondaryWriteControl"))
+    {
+        secondaryWriteControl_ = writeControlNames_.read
+        (
+            controlDict_.lookup("secondaryWriteControl")
+        );
+
+        if
+        (
+            controlDict_.readIfPresent
+            (
+                "secondaryWriteInterval",
+                secondaryWriteInterval_
+            )
+        )
+        {
+            if
+            (
+                secondaryWriteControl_
+             == wcTimeStep && label(secondaryWriteInterval_) < 1
+            )
+            {
+                FatalIOErrorIn("Time::readDict()", controlDict_)
+                    << "secondaryWriteInterval < 1"
+                    << " for secondaryWriteControl timeStep"
+                    << exit(FatalIOError);
+            }
+        }
+        else
+        {
+            controlDict_.lookup("secondaryWriteFrequency")
+                >> secondaryWriteInterval_;
+        }
+    }
+
+
+
     if (oldWriteInterval != writeInterval_)
     {
         switch (writeControl_)
@@ -310,4 +349,10 @@ bool Foam::Time::writeAndEnd()
 }
 
 
+void Foam::Time::writeOnce()
+{
+    writeOnce_ = true;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/db/objectRegistry/objectRegistry.C

@@ -163,11 +163,14 @@ Foam::label Foam::objectRegistry::getEvent() const
     label curEvent = event_++;
 
     if (event_ == labelMax)
+    {
+        if (objectRegistry::debug)
         {
             WarningIn("objectRegistry::getEvent() const")
                 << "Event counter has overflowed. "
                 << "Resetting counter on all dependent objects." << nl
                 << "This might cause extra evaluations." << endl;
+        }
 
         // Reset event counter
         curEvent = 1;

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C

@@ -124,11 +124,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     const DimensionedField<Type, GeoMesh>& df
 )
 :
-#   ifdef ConstructFromTmp
     regIOobject(df),
-#   else
-    regIOobject(df, true),
-#   endif
     Field<Type>(df),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_)
@@ -142,7 +138,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     bool reUse
 )
 :
-    regIOobject(df, true),
+    regIOobject(df, reUse),
     Field<Type>(df, reUse),
     mesh_(df.mesh_),
     dimensions_(df.dimensions_)
@@ -169,7 +165,7 @@ DimensionedField<Type, GeoMesh>::DimensionedField
     const tmp<DimensionedField<Type, GeoMesh> >& tdf
 )
 :
-    regIOobject(tdf(), true),
+    regIOobject(tdf(), tdf.isTmp()),
     Field<Type>
     (
         const_cast<DimensionedField<Type, GeoMesh>&>(tdf()),

src/OpenFOAM/global/argList/argList.C

@@ -726,7 +726,7 @@ Foam::argList::argList
     }
 
 
-    wordList slaveProcs;
+    stringList slaveProcs;
 
     // collect slave machine/pid
     if (parRunControl_.parRun())
@@ -734,7 +734,7 @@ Foam::argList::argList
         if (Pstream::master())
         {
             slaveProcs.setSize(Pstream::nProcs() - 1);
-            word  slaveMachine;
+            string  slaveMachine;
             label slavePid;
 
             label procI = 0;

src/OpenFOAM/include/OSspecific.H

@@ -69,13 +69,13 @@ bool setEnv(const word& name, const std::string& value, const bool overwrite);
 
 //- Return the system's host name, as per hostname(1)
 //  Optionally with the full name (as per the '-f' option)
-word hostName(const bool full=false);
+string hostName(const bool full=false);
 
 //- Return the system's domain name, as per hostname(1) with the '-d' option
-word domainName();
+string domainName();
 
 //- Return the user's login name
-word userName();
+string userName();
 
 //- Is user administrator
 bool isAdministrator();
@@ -84,7 +84,7 @@ bool isAdministrator();
 fileName home();
 
 //- Return home directory path name for a particular user
-fileName home(const word& userName);
+fileName home(const string& userName);
 
 //- Return current working directory path name
 fileName cwd();
@@ -189,10 +189,10 @@ unsigned int sleep(const unsigned int);
 void fdClose(const int);
 
 //- Check if machine is up by pinging given port
-bool ping(const word&, const label port, const label timeOut);
+bool ping(const string&, const label port, const label timeOut);
 
 //- Check if machine is up by pinging port 22 (ssh) and 222 (rsh)
-bool ping(const word&, const label timeOut=10);
+bool ping(const string&, const label timeOut=10);
 
 //- Execute the specified command
 int system(const std::string& command);

src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C

@@ -88,7 +88,7 @@ Foam::interpolationTable<Type>::interpolationTable(const fileName& fName)
     List<Tuple2<scalar, Type> >(),
     boundsHandling_(interpolationTable::WARN),
     fileName_(fName),
-    reader_(new openFoamTableReader<Type>())
+    reader_(new openFoamTableReader<Type>(dictionary()))
 {
     readTable();
 }

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -46,6 +46,20 @@ defineTypeNameAndDebug(Foam::globalMeshData, 0);
 // Geometric matching tolerance. Factor of mesh bounding box.
 const Foam::scalar Foam::globalMeshData::matchTol_ = 1E-8;
 
+namespace Foam
+{
+template<>
+class minEqOp<labelPair>
+{
+public:
+    void operator()(labelPair& x, const labelPair& y) const
+    {
+        x[0] = min(x[0], y[0]);
+        x[1] = min(x[1], y[1]);
+    }
+};
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -1063,6 +1077,128 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
 }
 
 
+void Foam::globalMeshData::calcGlobalEdgeOrientation() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
+            << " calculating edge orientation w.r.t. master edge." << endl;
+    }
+
+    const globalIndex& globalPoints = globalPointNumbering();
+
+    // 1. Determine master point
+    labelList masterPoint;
+    {
+        const mapDistribute& map = globalPointSlavesMap();
+
+        masterPoint.setSize(map.constructSize());
+        masterPoint = labelMax;
+
+        for (label pointI = 0; pointI < coupledPatch().nPoints(); pointI++)
+        {
+            masterPoint[pointI] = globalPoints.toGlobal(pointI);
+        }
+        syncData
+        (
+            masterPoint,
+            globalPointSlaves(),
+            globalPointTransformedSlaves(),
+            map,
+            minEqOp<label>()
+        );
+    }
+
+    // Now all points should know who is master by comparing their global
+    // pointID with the masterPointID. We now can use this information
+    // to find the orientation of the master edge.
+
+    {
+        const mapDistribute& map = globalEdgeSlavesMap();
+        const labelListList& slaves = globalEdgeSlaves();
+        const labelListList& transformedSlaves = globalEdgeTransformedSlaves();
+
+        // Distribute orientation of master edge (in masterPoint numbering)
+        labelPairList masterEdgeVerts(map.constructSize());
+        masterEdgeVerts = labelPair(labelMax, labelMax);
+
+        for (label edgeI = 0; edgeI < coupledPatch().nEdges(); edgeI++)
+        {
+            if
+            (
+                (
+                    slaves[edgeI].size()
+                  + transformedSlaves[edgeI].size()
+                )
+              > 0
+            )
+            {
+                // I am master. Fill in my masterPoint equivalent.
+
+                const edge& e = coupledPatch().edges()[edgeI];
+                masterEdgeVerts[edgeI] = labelPair
+                (
+                    masterPoint[e[0]],
+                    masterPoint[e[1]]
+                );
+            }
+        }
+        syncData
+        (
+            masterEdgeVerts,
+            slaves,
+            transformedSlaves,
+            map,
+            minEqOp<labelPair>()
+        );
+
+        // Now check my edges on how they relate to the master's edgeVerts
+        globalEdgeOrientationPtr_.reset
+        (
+            new PackedBoolList(coupledPatch().nEdges())
+        );
+        PackedBoolList& globalEdgeOrientation = globalEdgeOrientationPtr_();
+
+        forAll(coupledPatch().edges(), edgeI)
+        {
+            const edge& e = coupledPatch().edges()[edgeI];
+            const labelPair masterE
+            (
+                masterPoint[e[0]],
+                masterPoint[e[1]]
+            );
+
+            label stat = labelPair::compare
+            (
+                masterE,
+                masterEdgeVerts[edgeI]
+            );
+            if (stat == 0)
+            {
+                FatalErrorIn
+                (
+                    "globalMeshData::calcGlobalEdgeOrientation() const"
+                )   << "problem : my edge:" << e
+                    << " in master points:" << masterE
+                    << " v.s. masterEdgeVerts:" << masterEdgeVerts[edgeI]
+                    << exit(FatalError);
+            }
+            else
+            {
+                globalEdgeOrientation[edgeI] = (stat == 1);
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalEdgeOrientation() :"
+            << " finished calculating edge orientation."
+            << endl;
+    }
+}
+
+
 // Calculate uncoupled boundary faces (without calculating
 // primitiveMesh::pointFaces())
 void Foam::globalMeshData::calcPointBoundaryFaces
@@ -1660,6 +1796,7 @@ void Foam::globalMeshData::clearOut()
     globalEdgeNumberingPtr_.clear();
     globalEdgeSlavesPtr_.clear();
     globalEdgeTransformedSlavesPtr_.clear();
+    globalEdgeOrientationPtr_.clear();
     globalEdgeSlavesMapPtr_.clear();
 
     // Face
@@ -2095,6 +2232,16 @@ const
 }
 
 
+const Foam::PackedBoolList& Foam::globalMeshData::globalEdgeOrientation() const
+{
+    if (!globalEdgeOrientationPtr_.valid())
+    {
+        calcGlobalEdgeOrientation();
+    }
+    return globalEdgeOrientationPtr_();
+}
+
+
 const Foam::mapDistribute& Foam::globalMeshData::globalEdgeSlavesMap() const
 {
     if (!globalEdgeSlavesMapPtr_.valid())

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -95,6 +95,7 @@ class mapDistribute;
 template<class T> class EdgeMap;
 class globalIndex;
 class globalIndexAndTransform;
+class PackedBoolList;
 
 /*---------------------------------------------------------------------------*\
                            Class globalMeshData Declaration
@@ -191,6 +192,7 @@ class globalMeshData
             mutable autoPtr<globalIndex> globalEdgeNumberingPtr_;
             mutable autoPtr<labelListList> globalEdgeSlavesPtr_;
             mutable autoPtr<labelListList> globalEdgeTransformedSlavesPtr_;
+            mutable autoPtr<PackedBoolList> globalEdgeOrientationPtr_;
             mutable autoPtr<mapDistribute> globalEdgeSlavesMapPtr_;
 
 
@@ -297,6 +299,9 @@ class globalMeshData
             //- Calculate global edge addressing.
             void calcGlobalEdgeSlaves() const;
 
+            //- Calculate orientation w.r.t. edge master.
+            void calcGlobalEdgeOrientation() const;
+
 
         // Global boundary face/cell addressing
 
@@ -539,6 +544,8 @@ public:
                 const labelListList& globalEdgeSlaves() const;
                 const labelListList& globalEdgeTransformedSlaves() const;
                 const mapDistribute& globalEdgeSlavesMap() const;
+                //- Is my edge same orientation master edge
+                const PackedBoolList& globalEdgeOrientation() const;
 
             // Coupled point to boundary faces. These are uncoupled boundary
             // faces only but include empty patches.

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.C

@@ -30,6 +30,7 @@ License
 #include "PatchToolsSearch.C"
 #include "PatchToolsSortEdges.C"
 #include "PatchToolsNormals.C"
+#include "PatchToolsMatch.C"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -36,6 +36,7 @@ SourceFiles
     PatchToolsSearch.C
     PatchToolsSortEdges.C
     PatchToolsNormals.C
+    PatchToolsMatch.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -51,6 +52,7 @@ namespace Foam
 {
 
 class polyMesh;
+class PackedBoolList;
 
 /*---------------------------------------------------------------------------*\
                          Class PatchTools Declaration
@@ -169,6 +171,55 @@ public:
     );
 
 
+    //- Find corresponding points on patches sharing the same points
+    //  p1PointLabels : points on p1 that were matched
+    //  p2PointLabels : corresponding points on p2
+    template
+    <
+        class Face1,
+        template<class> class FaceList1,
+        class PointField1,
+        class PointType1,
+        class Face2,
+        template<class> class FaceList2,
+        class PointField2,
+        class PointType2
+    >
+    static void matchPoints
+    (
+        const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1,
+        const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2,
+
+        labelList& p1PointLabels,
+        labelList& p2PointLabels
+    );
+
+    //- Find corresponding edges on patches sharing the same points
+    //  p1EdgeLabels    : edges on p1 that were matched
+    //  p2EdgeLabels    : corresponding edges on p2
+    //  sameOrientation : same orientation?
+    template
+    <
+        class Face1,
+        template<class> class FaceList1,
+        class PointField1,
+        class PointType1,
+        class Face2,
+        template<class> class FaceList2,
+        class PointField2,
+        class PointType2
+    >
+    static void matchEdges
+    (
+        const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1,
+        const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2,
+
+        labelList& p1EdgeLabels,
+        labelList& p2EdgeLabels,
+        PackedBoolList& sameOrientation
+    );
+
+
     //- Return parallel consistent point normals for patches (on boundary faces)
     //  using mesh points.
     template

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsMatch.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template
+<
+    class Face1,
+    template<class> class FaceList1,
+    class PointField1,
+    class PointType1,
+    class Face2,
+    template<class> class FaceList2,
+    class PointField2,
+    class PointType2
+>
+void Foam::PatchTools::matchPoints
+(
+    const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1,
+    const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2,
+
+    labelList& p1PointLabels,
+    labelList& p2PointLabels
+)
+{
+    p1PointLabels.setSize(p1.nPoints());
+    p2PointLabels.setSize(p1.nPoints());
+
+    label nMatches = 0;
+
+    forAll(p1.meshPoints(), pointI)
+    {
+        label meshPointI = p1.meshPoints()[pointI];
+
+        Map<label>::const_iterator iter = p2.meshPointMap().find
+        (
+            meshPointI
+        );
+
+        if (iter != p2.meshPointMap().end())
+        {
+            p1PointLabels[nMatches] = pointI;
+            p2PointLabels[nMatches] = iter();
+            nMatches++;
+        }
+    }
+    p1PointLabels.setSize(nMatches);
+    p2PointLabels.setSize(nMatches);
+}
+
+
+template
+<
+    class Face1,
+    template<class> class FaceList1,
+    class PointField1,
+    class PointType1,
+    class Face2,
+    template<class> class FaceList2,
+    class PointField2,
+    class PointType2
+>
+void Foam::PatchTools::matchEdges
+(
+    const PrimitivePatch<Face1, FaceList1, PointField1, PointType1>& p1,
+    const PrimitivePatch<Face2, FaceList2, PointField2, PointType2>& p2,
+
+    labelList& p1EdgeLabels,
+    labelList& p2EdgeLabels,
+    PackedBoolList& sameOrientation
+)
+{
+    p1EdgeLabels.setSize(p1.nEdges());
+    p2EdgeLabels.setSize(p1.nEdges());
+    sameOrientation.setSize(p1.nEdges());
+    sameOrientation = 0;
+
+    label nMatches = 0;
+
+    EdgeMap<label> edgeToIndex(2*p1.nEdges());
+    forAll(p1.edges(), edgeI)
+    {
+        const edge& e = p1.edges()[edgeI];
+        const edge meshE
+        (
+            p1.meshPoints()[e[0]],
+            p1.meshPoints()[e[1]]
+        );
+        edgeToIndex.insert(meshE, edgeI);
+    }
+
+    forAll(p2.edges(), edgeI)
+    {
+        const edge& e = p2.edges()[edgeI];
+        const edge meshE(p2.meshPoints()[e[0]], p2.meshPoints()[e[1]]);
+
+        EdgeMap<label>::const_iterator iter = edgeToIndex.find(meshE);
+
+        if (iter != edgeToIndex.end())
+        {
+            p1EdgeLabels[nMatches] = iter();
+            p2EdgeLabels[nMatches] = edgeI;
+            sameOrientation[nMatches] = (meshE[0] == iter.key()[0]);
+            nMatches++;
+        }
+    }
+    p1EdgeLabels.setSize(nMatches);
+    p2EdgeLabels.setSize(nMatches);
+    sameOrientation.setSize(nMatches);
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/strings/stringOps/stringOps.C

@@ -548,7 +548,7 @@ Foam::string& Foam::stringOps::inplaceExpand
             //   ~OpenFOAM => site/user OpenFOAM configuration directory
             //   ~user     => home directory for specified user
 
-            word user;
+            string user;
             fileName file;
 
             if ((begVar = s.find('/')) != string::npos)

src/Pstream/dummy/UIPread.C

@@ -53,15 +53,42 @@ Foam::UIPstream::UIPstream
 {
     notImplemented
     (
-        "UIPstream::UIPstream"
-        "("
-            "const commsTypes,"
-            "const int fromProcNo,"
-            "DynamicList<char>&,"
-            "label&,"
-            "const int tag,"
-            "const bool,"
-            "streamFormat, versionNumber"
+        "UIPstream::UIPstream\n"
+        "(\n"
+            "const commsTypes,\n"
+            "const int,\n"
+            "DynamicList<char>&,\n"
+            "label&,\n"
+            "const int,\n"
+            "const bool,\n"
+            "streamFormat,\n"
+            "versionNumber\n"
+        ")"
+    );
+}
+
+
+Foam::UIPstream::UIPstream
+(
+    const int fromProcNo,
+    PstreamBuffers& buffers
+)
+:
+    UPstream(buffers.commsType_),
+    Istream(buffers.format_, buffers.version_),
+    fromProcNo_(fromProcNo),
+    externalBuf_(buffers.recvBuf_[fromProcNo]),
+    externalBufPosition_(buffers.recvBufPos_[fromProcNo]),
+    tag_(buffers.tag_),
+    clearAtEnd_(true),
+    messageSize_(0)
+{
+    notImplemented
+    (
+        "UIPstream::UIPstream\n"
+        "(\n"
+            "const int,\n"
+            "PstreamBuffers&\n"
         ")"
     );
 }

src/finiteVolume/Make/files

@@ -334,6 +334,7 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C
 $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
 $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
 $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
+$(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
 /*
 $(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
 $(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
@@ -388,6 +389,7 @@ $(basicSource)/basicSource/IObasicSourceList.C
 $(basicSource)/actuationDiskSource/actuationDiskSource.C
 $(basicSource)/radialActuationDiskSource/radialActuationDiskSource.C
 $(basicSource)/explicitSource/explicitSource.C
+$(basicSource)/explicitSetValue/explicitSetValue.C
 
 
 LIB = $(FOAM_LIBBIN)/libfiniteVolume

src/finiteVolume/cfdTools/general/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -48,14 +48,14 @@ void Foam::actuationDiskSource::addActuationDiskAxialInertialResistance
     E.xx() = uniDiskDir.x();
     E.yy() = uniDiskDir.y();
     E.zz() = uniDiskDir.z();
-    const vectorField U1((1.0 - a)*U);
+
     forAll(cells, i)
     {
-        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U1[cells[i]])*a/(1.0 - a);
+        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
     }
     forAll(cells, i)
     {
-        Usource[cells[i]] += ((Vcells[cells[i]]/V())*T[i]*E) & U1[cells[i]];
+        Usource[cells[i]] += ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]];
     }
 }
 

src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C

@@ -116,10 +116,15 @@ void Foam::basicSource::setCellSet()
                 label globalCellI = returnReduce(cellI, maxOp<label>());
                 if (globalCellI < 0)
                 {
-                    WarningIn("TimeActivatedExplicitSource<Type>::setCellIds()")
+                    WarningIn
+                    (
+                        "TimeActivatedExplicitSource<Type>::setCellIds()"
+                    )
                     << "Unable to find owner cell for point " << points_[i]
                     << endl;
+
                 }
+
             }
 
             cells_ = selectedCells.toc();
@@ -270,4 +275,30 @@ bool Foam::basicSource::isActive()
 }
 
 
+void Foam::basicSource::addSu(Foam::fvMatrix<vector>& Eqn)
+{
+    notImplemented
+    (
+        "Foam::basicSource addSu(Foam::fvMatrix<vector>& Eqn)"
+    );
+}
+
+
+void Foam::basicSource::addSu(Foam::fvMatrix<scalar>& Eqn)
+{
+    notImplemented
+    (
+        "Foam::basicSource addSu(Foam::fvMatrix<scalar>& Eqn)"
+    );
+}
+
+
+void Foam::basicSource::setValue(Foam::fvMatrix<scalar>& Eqn)
+{
+    notImplemented
+    (
+        "Foam::basicSource setValue(Foam::fvMatrix<scalar>& Eqn)"
+    );
+}
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.H

@@ -56,13 +56,13 @@ Description
         selectionMode   points;
         cellSet         c0;
 
-        explicitSourceCoeffs
-        {
         points            // list of points when selectionMode = points
         (
             (-0.088 0.007 -0.02)
             (-0.028 0.007 -0.02)
         );
+        explicitSourceCoeffs
+        {
             volumeMode      specific;  //absolute
             fieldData                  //field data
             {
@@ -323,10 +323,14 @@ public:
         // Evaluation
 
             //- Add source term to vector fvMatrix
-            virtual void addSu(fvMatrix<vector>& Eqn) = 0;
+            virtual void addSu(fvMatrix<vector>& Eqn);
 
             //- Add source term to scalar fvMatrix
-            virtual void addSu(fvMatrix<scalar>& Eqn) = 0;
+            virtual void addSu(fvMatrix<scalar>& Eqn);
+
+            //- Set constant value on field
+            virtual void setValue(fvMatrix<scalar>& Eq);
+
 
 
         // I-O

src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.C

@@ -90,6 +90,19 @@ void Foam::basicSourceList::addSu(fvMatrix<vector>& Eqn)
 }
 
 
+void Foam::basicSourceList::setValue(fvMatrix<scalar>& Eqn)
+{
+
+    forAll(*this, i)
+    {
+        if (this->operator[](i).isActive())
+        {
+            this->operator[](i).setValue(Eqn);
+        }
+    }
+}
+
+
 bool Foam::basicSourceList::read(const dictionary& dict)
 {
     bool allOk = true;

src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSourceList.H

@@ -92,6 +92,9 @@ public:
             //- Add source terms to vector fvMatrix
             void addSu(fvMatrix<vector>& Eq);
 
+            //- Set constant value on field
+            void setValue(fvMatrix<scalar>& Eq);
+
 
         // I-O
 

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "explicitSetValue.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(explicitSetValue, 0);
+    addToRunTimeSelectionTable
+    (
+        basicSource,
+        explicitSetValue,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::explicitSetValue::setFieldData(const dictionary& dict)
+{
+    scalarFields_.clear();
+    vectorFields_.clear();
+
+    wordList fieldTypes(dict.toc().size());
+    wordList fieldNames(dict.toc().size());
+
+    forAll(dict.toc(), i)
+    {
+        const word& fieldName = dict.toc()[i];
+        IOobject io
+        (
+            fieldName,
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+        if (io.headerOk())
+        {
+            fieldTypes[i] = io.headerClassName();
+            fieldNames[i] = dict.toc()[i];
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "explicitSetValue::setFieldData"
+            )   << "header not OK " << io.name()
+                << exit(FatalError);
+        }
+    }
+
+    addField(scalarFields_, fieldTypes, fieldNames, dict);
+    addField(vectorFields_, fieldTypes, fieldNames, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::explicitSetValue::explicitSetValue
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    basicSource(name, modelType, dict, mesh),
+    dict_(dict.subDict(modelType + "Coeffs"))
+{
+    setFieldData(dict_.subDict("fieldData"));
+}
+
+
+void Foam::explicitSetValue::setValue(fvMatrix<scalar>& Eqn)
+{
+    setFieldValue(Eqn, scalarFields_[Eqn.psi().name()]);
+}
+
+
+void Foam::explicitSetValue::setValue(fvMatrix<vector>& Eqn)
+{
+    setFieldValue(Eqn, vectorFields_[Eqn.psi().name()]);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValue.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::explicitSetValue
+
+Description
+    Explicit set values on fields.
+
+    Sources described by:
+
+    explicitSetValueCoeffs
+    {
+        fieldData            // field data - usage for multiple fields
+        {
+            k   30.7;
+            epsilon  1.5;
+        }
+    }
+
+SourceFiles
+    explicitSetValue.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef explicitSetValue_H
+#define explicitSetValue_H
+
+#include "cellSet.H"
+#include "volFieldsFwd.H"
+#include "DimensionedField.H"
+#include "basicSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class explicitSetValue Declaration
+\*---------------------------------------------------------------------------*/
+
+class explicitSetValue
+:
+    public basicSource
+{
+    // Private data
+
+        //- List of field types
+        HashTable<scalar> scalarFields_;
+        HashTable<vector> vectorFields_;
+
+        //- Set value to field
+        template<class Type>
+        void setFieldValue(fvMatrix<Type>&, const Type&) const;
+
+        //- Add field names and values to field table for types.
+        template<class Type>
+        void addField
+        (
+            HashTable<Type>& fields,
+            const wordList& fieldTypes,
+            const wordList& fieldNames,
+            const dictionary& dict
+        );
+
+
+protected:
+
+    // Protected data
+
+         //- Sub dictionary for time activated explicit sources
+        const dictionary& dict_;
+
+
+    // Protected functions
+
+        //- Set the local field data
+        void setFieldData(const dictionary& dict);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("explicitSetValue");
+
+
+    // Constructors
+
+        //- Construct from components
+        explicitSetValue
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- Return clone
+        autoPtr<explicitSetValue> clone() const
+        {
+            notImplemented
+            (
+                "autoPtr<explicitSetValue> clone() const"
+            );
+            return autoPtr<explicitSetValue>(NULL);
+        }
+
+
+    // Member Functions
+
+
+        // Edit
+
+            //- Return points
+            inline  const List<point>& points() const;
+
+
+        // Evaluation
+
+            //- Set value on vector field
+            virtual void setValue(fvMatrix<vector>& UEqn);
+
+            //- Set value on scalar field
+            virtual void setValue(fvMatrix<scalar>& UEqn);
+
+
+        // I-O
+
+            //- Write the source properties
+            virtual void writeData(Ostream&) const;
+
+            //- Read fieldData in sub-dictionary
+            virtual bool read(const dictionary& dict);
+
+            //- Ostream operator
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const explicitSetValue& source
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "explicitSetValueIO.C"
+#include "explicitSetValueI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "explicitSetValueTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueI.H

@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "explicitSetValue.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::List<Foam::point>&
+Foam::explicitSetValue::points() const
+{
+    return points_;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueIO.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "explicitSetValue.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::explicitSetValue::writeData(Ostream& os) const
+{
+    os  << indent << name_ << nl
+        << indent << token::BEGIN_BLOCK << incrIndent << nl;
+
+    if (scalarFields_.size() > 0)
+    {
+        os.writeKeyword("scalarFields") << scalarFields_
+            << token::END_STATEMENT << nl;
+    }
+
+    if (vectorFields_.size() > 0)
+    {
+        os.writeKeyword("vectorFields") << vectorFields_
+            << token::END_STATEMENT << nl;
+    }
+
+    os << decrIndent << indent << token::END_BLOCK << endl;
+}
+
+
+bool Foam::explicitSetValue::read(const dictionary& dict)
+{
+    if (basicSource::read(dict))
+    {
+        const dictionary& sourceDict = dict.subDict(name());
+        const dictionary& subDictCoeffs = sourceDict.subDict
+        (
+            typeName + "Coeffs"
+        );
+        setFieldData(subDictCoeffs.subDict("fieldData"));
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const explicitSetValue& source)
+{
+    source.writeData(os);
+    return os;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSetValue/explicitSetValueTemplates.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+template <class Type>
+void Foam::explicitSetValue::setFieldValue
+(
+    fvMatrix<Type>& Eqn,
+    const Type& value
+) const
+{
+    Type data = value;
+
+    DimensionedField<Type, volMesh> rhs
+    (
+        IOobject
+        (
+            "rhs",
+            Eqn.psi().mesh().time().timeName(),
+            Eqn.psi().mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        Eqn.psi().mesh(),
+        dimensioned<Type>
+        (
+            "zero",
+            dimless,
+            pTraits<Type>::zero
+        )
+    );
+
+    List<Type> values(this->cells().size());
+
+    forAll (values, i)
+    {
+        values[i] = data;
+    }
+
+    Eqn.setValues(this->cells(), values);
+}
+
+
+template <class Type>
+void Foam::explicitSetValue::addField
+(
+    HashTable<Type>& fields,
+    const wordList& fieldTypes,
+    const wordList& fieldNames,
+    const dictionary& fieldDataDict
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> geometricField;
+
+    forAll (fieldTypes, fieldI)
+    {
+        word fieldName = fieldNames[fieldI];
+        word fieldType = fieldTypes[fieldI];
+
+        if
+        (
+            (
+                fieldType
+             == GeometricField<Type, fvPatchField, volMesh>::typeName
+            ) &&
+            (
+                this->mesh().foundObject<geometricField>(fieldName)
+            )
+        )
+        {
+            Type fieldValue = fieldDataDict.lookupOrDefault<Type>
+            (
+                fieldName,
+                pTraits<Type>::zero
+            );
+
+            fields.insert(fieldName, fieldValue);
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/fieldSources/basicSource/explicitSource/explicitSource.C

@@ -75,7 +75,7 @@ void Foam::explicitSource::setFieldData(const dictionary& dict)
         IOobject io
         (
             fieldName,
-            this->mesh().time().timeName(0),
+            this->mesh().time().timeName(),
             this->mesh(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,

src/finiteVolume/cfdTools/general/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C

@@ -53,7 +53,6 @@ addRadialActuationDiskAxialInertialResistance
     E.xx() = uniDiskDir.x();
     E.yy() = uniDiskDir.y();
     E.zz() = uniDiskDir.z();
-    const vectorField U1((1.0 - a)*U);
 
     const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
     const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);
@@ -68,7 +67,7 @@ addRadialActuationDiskAxialInertialResistance
 
     forAll(cells, i)
     {
-        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U1[cells[i]])*a/(1.0 - a);
+        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
 
         scalar r = mag(mesh().cellCentres()[cells[i]] - avgCentre);
 
@@ -79,7 +78,7 @@ addRadialActuationDiskAxialInertialResistance
 
     forAll(cells, i)
     {
-        Usource[cells[i]] += ((Vcells[cells[i]]/V())*Tr[i]*E) & U1[cells[i]];
+        Usource[cells[i]] += ((Vcells[cells[i]]/V())*Tr[i]*E) & U[cells[i]];
     }
 
     if (debug)

src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C

@@ -110,38 +110,34 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
     const labelUList& owner = mesh_.owner();
     const labelUList& neighbour = mesh_.neighbour();
 
-    // Build the d-vectors
-    surfaceVectorField d
-    (
-        mesh_.Sf()/(mesh_.magSf()*mesh_.deltaCoeffs())
-    );
-
-    if (!mesh_.orthogonal())
-    {
-        d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
-    }
-
+    const volVectorField& C = mesh.C();
 
     // Set up temporary storage for the dd tensor (before inversion)
     symmTensorField dd(mesh_.nCells(), symmTensor::zero);
 
-    forAll(owner, faceI)
+    forAll(owner, facei)
     {
-        const symmTensor wdd(1.0/magSqr(d[faceI])*sqr(d[faceI]));
+        label own = owner[facei];
+        label nei = neighbour[facei];
 
-        dd[owner[faceI]] += wdd;
-        dd[neighbour[faceI]] += wdd;
+        vector d = C[nei] - C[own];
+
+        const symmTensor wdd(1.0/magSqr(d[facei])*sqr(d[facei]));
+
+        dd[own] += wdd;
+        dd[nei] += wdd;
     }
 
     // Visit the boundaries. Coupled boundaries are taken into account
     // in the construction of d vectors.
-    forAll(d.boundaryField(), patchI)
+    forAll(lsP.boundaryField(), patchi)
     {
-        const fvsPatchVectorField& pd = d.boundaryField()[patchI];
-
-        const fvPatch& p = pd.patch();
+        const fvPatch& p = lsP.boundaryField()[patchi].patch();
         const labelUList& faceCells = p.faceCells();
 
+        // Build the d-vectors
+        vectorField pd = p.delta();
+
         forAll(pd, patchFaceI)
         {
             dd[faceCells[patchFaceI]] +=
@@ -232,24 +228,30 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
 
 
     // Revisit all faces and calculate the lsP and lsN vectors
-    forAll(owner, faceI)
+    forAll(owner, facei)
     {
-        lsP[faceI] =
-            (1.0/magSqr(d[faceI]))*(invDd[owner[faceI]] & d[faceI]);
+        label own = owner[facei];
+        label nei = neighbour[facei];
 
-        lsN[faceI] =
-            ((-1.0)/magSqr(d[faceI]))*(invDd[neighbour[faceI]] & d[faceI]);
+        vector d = C[nei] - C[own];
+
+        lsP[facei] =
+            (1.0/magSqr(d[facei]))*(invDd[owner[facei]] & d);
+
+        lsN[facei] =
+            ((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d);
     }
 
     forAll(lsP.boundaryField(), patchI)
     {
-        const fvsPatchVectorField& pd = d.boundaryField()[patchI];
-
         fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
 
         const fvPatch& p = patchLsP.patch();
         const labelUList& faceCells = p.faceCells();
 
+        // Build the d-vectors
+        vectorField pd = p.delta();
+
         forAll(p, patchFaceI)
         {
             patchLsP[patchFaceI] =

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C

@@ -121,24 +121,12 @@ Foam::fv::fourthGrad<Type>::calcGrad
 
             const scalarField& lambdap = lambda.boundaryField()[patchi];
 
+            const fvPatch& p = fGrad.boundaryField()[patchi].patch();
+
+            const labelUList& faceCells = p.faceCells();
+
             // Build the d-vectors
-            vectorField pd
-            (
-                mesh.Sf().boundaryField()[patchi]
-              / (
-                    mesh.magSf().boundaryField()[patchi]
-                  * mesh.deltaCoeffs().boundaryField()[patchi]
-                )
-            );
-
-            if (!mesh.orthogonal())
-            {
-                pd -= mesh.correctionVectors().boundaryField()[patchi]
-                     /mesh.deltaCoeffs().boundaryField()[patchi];
-            }
-
-            const labelUList& faceCells =
-                fGrad.boundaryField()[patchi].patch().faceCells();
+            vectorField pd = p.delta();
 
             const Field<GradType> neighbourSecondfGrad
             (

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C

@@ -134,7 +134,7 @@ Foam::fv::gradScheme<Type>::grad
             regIOobject::store(tgGrad.ptr());
         }
 
-        cachePrintMessage("Retreiving", name, vsf);
+        cachePrintMessage("Retrieving", name, vsf);
         GradFieldType& gGrad = const_cast<GradFieldType&>
         (
             mesh().objectRegistry::template lookupObject<GradFieldType>(name)

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C

@@ -130,20 +130,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         const labelUList& faceCells = p.patch().faceCells();
 
         // Build the d-vectors
-        vectorField pd
-        (
-            mesh.Sf().boundaryField()[patchi]
-          / (
-                mesh.magSf().boundaryField()[patchi]
-              * mesh.deltaCoeffs().boundaryField()[patchi]
-           )
-        );
-
-        if (!mesh.orthogonal())
-        {
-            pd -= mesh.correctionVectors().boundaryField()[patchi]
-                /mesh.deltaCoeffs().boundaryField()[patchi];
-        }
+        vectorField pd = p.delta();
 
         if (p.coupled())
         {
@@ -196,21 +183,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         const labelUList& faceCells = p.faceCells();
 
         // Build the d-vectors
-        vectorField pd
-        (
-            mesh.Sf().boundaryField()[patchi]
-           /(
-               mesh.magSf().boundaryField()[patchi]
-              *mesh.deltaCoeffs().boundaryField()[patchi]
-            )
-        );
-
-        if (!mesh.orthogonal())
-        {
-            pd -= mesh.correctionVectors().boundaryField()[patchi]
-                /mesh.deltaCoeffs().boundaryField()[patchi];
-        }
-
+        vectorField pd = p.delta();
 
         if (p.coupled())
         {

src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H

@@ -211,7 +211,7 @@ public:
 
 // Add the patch constructor functions to the hash tables
 
-#define makeFvLaplacianTypeScheme(SS, Type, GType)                             \
+#define makeFvLaplacianTypeScheme(SS, GType, Type)                             \
                                                                                \
     typedef SS<Type, GType> SS##Type##GType;                                   \
     defineNamedTemplateTypeNameAndDebug(SS##Type##GType, 0);                   \
@@ -224,13 +224,20 @@ public:
 #define makeFvLaplacianScheme(SS)                                              \
                                                                                \
 makeFvLaplacianTypeScheme(SS, scalar, scalar)                                  \
-makeFvLaplacianTypeScheme(SS, scalar, symmTensor)                              \
-makeFvLaplacianTypeScheme(SS, scalar, tensor)                                  \
-makeFvLaplacianTypeScheme(SS, vector, scalar)                                  \
-makeFvLaplacianTypeScheme(SS, sphericalTensor, scalar)                         \
 makeFvLaplacianTypeScheme(SS, symmTensor, scalar)                              \
 makeFvLaplacianTypeScheme(SS, tensor, scalar)                                  \
-
+makeFvLaplacianTypeScheme(SS, scalar, vector)                                  \
+makeFvLaplacianTypeScheme(SS, symmTensor, vector)                              \
+makeFvLaplacianTypeScheme(SS, tensor, vector)                                  \
+makeFvLaplacianTypeScheme(SS, scalar, sphericalTensor)                         \
+makeFvLaplacianTypeScheme(SS, symmTensor, sphericalTensor)                     \
+makeFvLaplacianTypeScheme(SS, tensor, sphericalTensor)                         \
+makeFvLaplacianTypeScheme(SS, scalar, symmTensor)                              \
+makeFvLaplacianTypeScheme(SS, symmTensor, symmTensor)                          \
+makeFvLaplacianTypeScheme(SS, tensor, symmTensor)                              \
+makeFvLaplacianTypeScheme(SS, scalar, tensor)                                  \
+makeFvLaplacianTypeScheme(SS, symmTensor, tensor)                              \
+makeFvLaplacianTypeScheme(SS, tensor, tensor)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C

@@ -50,7 +50,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection
 
     // construct GeometricField<Type, fvsPatchField, surfaceMesh>
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf =
-        mesh.correctionVectors()
+        mesh.nonOrthCorrectionVectors()
       & linear<typename outerProduct<vector, Type>::type>(mesh).interpolate
         (
             gradScheme<Type>::New
@@ -88,7 +88,7 @@ Foam::fv::correctedSnGrad<Type>::correction
                 IOobject::NO_WRITE
             ),
             mesh,
-            vf.dimensions()*mesh.deltaCoeffs().dimensions()
+            vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
         )
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf();
@@ -98,7 +98,7 @@ Foam::fv::correctedSnGrad<Type>::correction
         ssf.replace
         (
             cmpt,
-            mesh.correctionVectors()
+            mesh.nonOrthCorrectionVectors()
           & linear
             <
                 typename

src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H

@@ -96,13 +96,13 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const
         {
-            return this->mesh().deltaCoeffs();
+            return this->mesh().nonOrthDeltaCoeffs();
         }
 
         //- Return true if this scheme uses an explicit correction
         virtual bool corrected() const
         {
-            return !this->mesh().orthogonal();
+            return true;
         }
 
         //- Return the explicit correction to the correctedSnGrad

src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.H

@@ -120,13 +120,13 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const
         {
-            return this->mesh().deltaCoeffs();
+            return this->mesh().nonOrthDeltaCoeffs();
         }
 
         //- Return true if this scheme uses an explicit correction
         virtual bool corrected() const
         {
-            return !this->mesh().orthogonal();
+            return true;
         }
 
         //- Return the explicit correction to the limitedSnGrad

src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Simple central-difference snGrad scheme without non-orthogonal correction.
+
+\*---------------------------------------------------------------------------*/
+
+#include "orthogonalSnGrad.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+orthogonalSnGrad<Type>::~orthogonalSnGrad()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+orthogonalSnGrad<Type>::correction
+(
+    const GeometricField<Type, fvPatchField, volMesh>&
+) const
+{
+    notImplemented
+    (
+        "orthogonalSnGrad<Type>::correction"
+        "(const GeometricField<Type, fvPatchField, volMesh>&)"
+    );
+    return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >(NULL);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrad.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::orthogonalSnGrad
+
+Description
+    Simple central-difference snGrad scheme without non-orthogonal correction.
+
+SourceFiles
+    orthogonalSnGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef orthogonalSnGrad_H
+#define orthogonalSnGrad_H
+
+#include "snGradScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class orthogonalSnGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class orthogonalSnGrad
+:
+    public snGradScheme<Type>
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const orthogonalSnGrad&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("uncorrected");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        orthogonalSnGrad(const fvMesh& mesh)
+        :
+            snGradScheme<Type>(mesh)
+        {}
+
+
+        //- Construct from mesh and data stream
+        orthogonalSnGrad(const fvMesh& mesh, Istream&)
+        :
+            snGradScheme<Type>(mesh)
+        {}
+
+
+    //- Destructor
+    virtual ~orthogonalSnGrad();
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors for the given field
+        virtual tmp<surfaceScalarField> deltaCoeffs
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const
+        {
+            return this->mesh().deltaCoeffs();
+        }
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return false;
+        }
+
+        //- Return the explicit correction to the orthogonalSnGrad
+        //  for the given field
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction(const GeometricField<Type, fvPatchField, volMesh>&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "orthogonalSnGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/snGradSchemes/orthogonalSnGrad/orthogonalSnGrads.C

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Simple central-difference snGrad scheme without non-orthogonal correction.
+
+\*---------------------------------------------------------------------------*/
+
+#include "orthogonalSnGrad.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    makeSnGradScheme(orthogonalSnGrad)
+}
+}
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H

@@ -109,7 +109,7 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const
         {
-            return this->mesh().deltaCoeffs();
+            return this->mesh().nonOrthDeltaCoeffs();
         }
 
         //- Return true if this scheme uses an explicit correction

src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C

@@ -247,7 +247,7 @@ Foam::label Foam::quadraticFitSnGradData::calcFit
     scalarList singVals(minSize_);
     label nSVDzeros = 0;
 
-    const scalar deltaCoeff = mesh().deltaCoeffs()[faci];
+    const scalar deltaCoeff = deltaCoeffs()[faci];
 
     bool goodFit = false;
     for (int iIt = 0; iIt < 10 && !goodFit; iIt++)

src/finiteVolume/finiteVolume/snGradSchemes/uncorrectedSnGrad/uncorrectedSnGrad.H

@@ -96,7 +96,7 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const
         {
-            return this->mesh().deltaCoeffs();
+            return this->mesh().nonOrthDeltaCoeffs();
         }
 
         //- Return true if this scheme uses an explicit correction

src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H

@@ -66,9 +66,6 @@ protected:
         //- Make patch weighting factors
         virtual void makeWeights(scalarField&) const = 0;
 
-        //- Make patch face - neighbour cell distances
-        virtual void makeDeltaCoeffs(scalarField&) const = 0;
-
 
 public:
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C

@@ -57,25 +57,6 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const
 }
 
 
-// Make patch face - neighbour cell distances
-void Foam::cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const
-{
-    //const cyclicPolyPatch& nbrPatch = cyclicPolyPatch_.neighbPatch();
-    const cyclicFvPatch& nbrPatch = neighbFvPatch();
-
-    const scalarField deltas(nf() & fvPatch::delta());
-    const scalarField nbrDeltas(nbrPatch.nf() & nbrPatch.fvPatch::delta());
-
-    forAll(deltas, facei)
-    {
-        scalar di = deltas[facei];
-        scalar dni = nbrDeltas[facei];
-
-        dc[facei] = 1.0/(di + dni);
-    }
-}
-
-
 // Return delta (P to N) vectors across coupled patch
 Foam::tmp<Foam::vectorField> Foam::cyclicFvPatch::delta() const
 {

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H

@@ -66,9 +66,6 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
-        //- Make patch face - neighbour cell distances
-        void makeDeltaCoeffs(scalarField&) const;
-
 
 public:
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C

@@ -60,27 +60,6 @@ void Foam::cyclicAMIFvPatch::makeWeights(scalarField& w) const
 }
 
 
-void Foam::cyclicAMIFvPatch::makeDeltaCoeffs(scalarField& dc) const
-{
-    const cyclicAMIFvPatch& nbrPatch = neighbFvPatch();
-
-    const scalarField deltas(nf() & fvPatch::delta());
-
-    const scalarField nbrDeltas
-    (
-        interpolate(nbrPatch.nf() & nbrPatch.fvPatch::delta())
-    );
-
-    forAll(deltas, faceI)
-    {
-        scalar di = deltas[faceI];
-        scalar dni = nbrDeltas[faceI];
-
-        dc[faceI] = 1.0/(di + dni);
-    }
-}
-
-
 Foam::tmp<Foam::vectorField> Foam::cyclicAMIFvPatch::delta() const
 {
     const vectorField patchD(fvPatch::delta());

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H

@@ -66,9 +66,6 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
-        //- Make patch face - neighbour cell distances
-        void makeDeltaCoeffs(scalarField&) const;
-
 
 public:
 

src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C

@@ -65,19 +65,6 @@ void processorFvPatch::makeWeights(scalarField& w) const
 }
 
 
-void processorFvPatch::makeDeltaCoeffs(scalarField& dc) const
-{
-    if (Pstream::parRun())
-    {
-        dc = (1.0 - weights())/(nf() & fvPatch::delta());
-    }
-    else
-    {
-        dc = 1.0/(nf() & fvPatch::delta());
-    }
-}
-
-
 tmp<vectorField> processorFvPatch::delta() const
 {
     if (Pstream::parRun())

src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H

@@ -65,9 +65,6 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
-        //- Make patch face - neighbour cell distances
-        void makeDeltaCoeffs(scalarField&) const;
-
 
 public:
 

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C

@@ -150,12 +150,6 @@ void Foam::fvPatch::makeWeights(scalarField& w) const
 }
 
 
-void Foam::fvPatch::makeDeltaCoeffs(scalarField& dc) const
-{
-    dc = 1.0/(nf() & delta());
-}
-
-
 void Foam::fvPatch::initMovePoints()
 {}
 

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H

@@ -86,9 +86,6 @@ protected:
         //- Make patch weighting factors
         virtual void makeWeights(scalarField&) const;
 
-        //- Make patch face - neighbour cell distances
-        virtual void makeDeltaCoeffs(scalarField&) const;
-
         //- Initialise the patches for moving points
         virtual void initMovePoints();
 

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

@@ -122,20 +122,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
             );
 
             // Build the d-vectors
-            vectorField pd
-            (
-                mesh.Sf().boundaryField()[patchi]
-              / (
-                    mesh.magSf().boundaryField()[patchi]
-                  * mesh.deltaCoeffs().boundaryField()[patchi]
-                )
-            );
-
-            if (!mesh.orthogonal())
-            {
-                pd -= mesh.correctionVectors().boundaryField()[patchi]
-                    /mesh.deltaCoeffs().boundaryField()[patchi];
-            }
+            vectorField pd = CDweights.boundaryField()[patchi].patch().delta();
 
             forAll(pLim, face)
             {

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

@@ -111,20 +111,7 @@ Foam::linearUpwind<Type>::correction
             );
 
             // Build the d-vectors
-            vectorField pd
-            (
-                mesh.Sf().boundaryField()[patchi]
-              / (
-                    mesh.magSf().boundaryField()[patchi]
-                  * mesh.deltaCoeffs().boundaryField()[patchi]
-                )
-            );
-
-            if (!mesh.orthogonal())
-            {
-                pd -= mesh.correctionVectors().boundaryField()[patchi]
-                    /mesh.deltaCoeffs().boundaryField()[patchi];
-            }
+            vectorField pd = Cf.boundaryField()[patchi].patch().delta();
 
             forAll(pOwner, facei)
             {

src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C

@@ -82,21 +82,17 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
     const surfaceVectorField& Sf = mesh_.Sf();
 
     const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
 
-    // Build the d-vectors
-    surfaceVectorField d(Sf/(mesh_.magSf()*mesh_.deltaCoeffs()));
-
-    if (!mesh_.orthogonal())
+    forAll(owner, facei)
     {
-        d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
-    }
+        label own = owner[facei];
+        label nei = neighbour[facei];
 
-    forAll(owner, faceI)
-    {
-        vector Cpf = Cf[faceI] - C[owner[faceI]];
+        vector d = C[nei] - C[own];
+        vector Cpf = Cf[facei] - C[own];
 
-        SkewCorrVecs[faceI] =
-            Cpf - ((Sf[faceI] & Cpf)/(Sf[faceI] & d[faceI]))*d[faceI];
+        SkewCorrVecs[facei] = Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
     }
 
 
@@ -115,7 +111,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
             const labelUList& faceCells = p.faceCells();
             const vectorField& patchFaceCentres = Cf.boundaryField()[patchI];
             const vectorField& patchSf = Sf.boundaryField()[patchI];
-            const vectorField& patchD = d.boundaryField()[patchI];
+            const vectorField patchD = p.delta();
 
             forAll(p, patchFaceI)
             {
@@ -136,7 +132,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
 
     if (Sf.internalField().size())
     {
-        skewCoeff = max(mag(SkewCorrVecs)/mag(d)).value();
+        skewCoeff = max(mag(SkewCorrVecs)*mesh_.deltaCoeffs()).value();
     }
 
     if (debug)
@@ -182,7 +178,7 @@ const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const
     if (!skew())
     {
         FatalErrorIn("skewCorrectionVectors::operator()()")
-            << "Cannot return correctionVectors; mesh is not skewed"
+            << "Cannot return skewCorrectionVectors; mesh is not skewed"
             << abort(FatalError);
     }
 

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C

@@ -31,43 +31,38 @@ Description
 #include "surfaceFields.H"
 #include "demandDrivenData.H"
 #include "coupledFvPatch.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(surfaceInterpolation, 0);
+defineTypeNameAndDebug(Foam::surfaceInterpolation, 0);
 
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void surfaceInterpolation::clearOut()
+void Foam::surfaceInterpolation::clearOut()
 {
-    deleteDemandDrivenData(weightingFactors_);
-    deleteDemandDrivenData(differenceFactors_);
-    deleteDemandDrivenData(correctionVectors_);
+    deleteDemandDrivenData(weights_);
+    deleteDemandDrivenData(deltaCoeffs_);
+    deleteDemandDrivenData(nonOrthDeltaCoeffs_);
+    deleteDemandDrivenData(nonOrthCorrectionVectors_);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
-surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm)
+Foam::surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm)
 :
     mesh_(fvm),
-    weightingFactors_(NULL),
-    differenceFactors_(NULL),
-    orthogonal_(false),
-    correctionVectors_(NULL)
+    weights_(NULL),
+    deltaCoeffs_(NULL),
+    nonOrthDeltaCoeffs_(NULL),
+    nonOrthCorrectionVectors_(NULL)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
 
-surfaceInterpolation::~surfaceInterpolation()
+Foam::surfaceInterpolation::~surfaceInterpolation()
 {
     clearOut();
 }
@@ -75,66 +70,67 @@ surfaceInterpolation::~surfaceInterpolation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const surfaceScalarField& surfaceInterpolation::weights() const
+const Foam::surfaceScalarField&
+Foam::surfaceInterpolation::weights() const
 {
-    if (!weightingFactors_)
+    if (!weights_)
     {
         makeWeights();
     }
 
-    return (*weightingFactors_);
+    return (*weights_);
 }
 
 
-const surfaceScalarField& surfaceInterpolation::deltaCoeffs() const
+const Foam::surfaceScalarField&
+Foam::surfaceInterpolation::deltaCoeffs() const
 {
-    if (!differenceFactors_)
+    if (!deltaCoeffs_)
     {
         makeDeltaCoeffs();
     }
 
-    return (*differenceFactors_);
+    return (*deltaCoeffs_);
 }
 
 
-bool surfaceInterpolation::orthogonal() const
+const Foam::surfaceScalarField&
+Foam::surfaceInterpolation::nonOrthDeltaCoeffs() const
 {
-    if (orthogonal_ == false && !correctionVectors_)
+    if (!nonOrthDeltaCoeffs_)
     {
-        makeCorrectionVectors();
+        makeNonOrthDeltaCoeffs();
     }
 
-    return orthogonal_;
+    return (*nonOrthDeltaCoeffs_);
 }
 
 
-const surfaceVectorField& surfaceInterpolation::correctionVectors() const
+const Foam::surfaceVectorField&
+Foam::surfaceInterpolation::nonOrthCorrectionVectors() const
 {
-    if (orthogonal())
+    if (!nonOrthCorrectionVectors_)
     {
-        FatalErrorIn("surfaceInterpolation::correctionVectors()")
-            << "cannot return correctionVectors; mesh is orthogonal"
-            << abort(FatalError);
+        makeNonOrthCorrectionVectors();
     }
 
-    return (*correctionVectors_);
+    return (*nonOrthCorrectionVectors_);
 }
 
 
 // Do what is neccessary if the mesh has moved
-bool surfaceInterpolation::movePoints()
+bool Foam::surfaceInterpolation::movePoints()
 {
-    deleteDemandDrivenData(weightingFactors_);
-    deleteDemandDrivenData(differenceFactors_);
-
-    orthogonal_ = false;
-    deleteDemandDrivenData(correctionVectors_);
+    deleteDemandDrivenData(weights_);
+    deleteDemandDrivenData(deltaCoeffs_);
+    deleteDemandDrivenData(nonOrthDeltaCoeffs_);
+    deleteDemandDrivenData(nonOrthCorrectionVectors_);
 
     return true;
 }
 
 
-void surfaceInterpolation::makeWeights() const
+void Foam::surfaceInterpolation::makeWeights() const
 {
     if (debug)
     {
@@ -143,20 +139,18 @@ void surfaceInterpolation::makeWeights() const
             << endl;
     }
 
-
-    weightingFactors_ = new surfaceScalarField
+    weights_ = new surfaceScalarField
     (
         IOobject
         (
-            "weightingFactors",
+            "weights",
             mesh_.pointsInstance(),
             mesh_
         ),
         mesh_,
         dimless
     );
-    surfaceScalarField& weightingFactors = *weightingFactors_;
-
+    surfaceScalarField& weights = *weights_;
 
     // Set local references to mesh data
     // (note that we should not use fvMesh sliced fields at this point yet
@@ -170,7 +164,7 @@ void surfaceInterpolation::makeWeights() const
     const vectorField& Sf = mesh_.faceAreas();
 
     // ... and reference to the internal field of the weighting factors
-    scalarField& w = weightingFactors.internalField();
+    scalarField& w = weights.internalField();
 
     forAll(owner, facei)
     {
@@ -188,11 +182,10 @@ void surfaceInterpolation::makeWeights() const
     {
         mesh_.boundary()[patchi].makeWeights
         (
-            weightingFactors.boundaryField()[patchi]
+            weights.boundaryField()[patchi]
         );
     }
 
-
     if (debug)
     {
         Info<< "surfaceInterpolation::makeWeights() : "
@@ -202,7 +195,7 @@ void surfaceInterpolation::makeWeights() const
 }
 
 
-void surfaceInterpolation::makeDeltaCoeffs() const
+void Foam::surfaceInterpolation::makeDeltaCoeffs() const
 {
     if (debug)
     {
@@ -215,18 +208,63 @@ void surfaceInterpolation::makeDeltaCoeffs() const
     // needed to make sure deltaCoeffs are calculated for parallel runs.
     weights();
 
-    differenceFactors_ = new surfaceScalarField
+    deltaCoeffs_ = new surfaceScalarField
     (
         IOobject
         (
-            "differenceFactors_",
+            "deltaCoeffs",
             mesh_.pointsInstance(),
             mesh_
         ),
         mesh_,
         dimless/dimLength
     );
-    surfaceScalarField& DeltaCoeffs = *differenceFactors_;
+    surfaceScalarField& DeltaCoeffs = *deltaCoeffs_;
+
+
+    // Set local references to mesh data
+    const volVectorField& C = mesh_.C();
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+
+    forAll(owner, facei)
+    {
+        DeltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
+    }
+
+    forAll(DeltaCoeffs.boundaryField(), patchi)
+    {
+        DeltaCoeffs.boundaryField()[patchi] =
+            1.0/mag(mesh_.boundary()[patchi].delta());
+    }
+}
+
+
+void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
+{
+    if (debug)
+    {
+        Info<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
+            << "Constructing differencing factors array for face gradient"
+            << endl;
+    }
+
+    // Force the construction of the weighting factors
+    // needed to make sure deltaCoeffs are calculated for parallel runs.
+    weights();
+
+    nonOrthDeltaCoeffs_ = new surfaceScalarField
+    (
+        IOobject
+        (
+            "nonOrthDeltaCoeffs",
+            mesh_.pointsInstance(),
+            mesh_
+        ),
+        mesh_,
+        dimless/dimLength
+    );
+    surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
 
 
     // Set local references to mesh data
@@ -242,49 +280,49 @@ void surfaceInterpolation::makeDeltaCoeffs() const
         vector unitArea = Sf[facei]/magSf[facei];
 
         // Standard cell-centre distance form
-        //DeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
+        //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
 
         // Slightly under-relaxed form
-        //DeltaCoeffs[facei] = 1.0/mag(delta);
+        //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
 
         // More under-relaxed form
-        //DeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
+        //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
 
         // Stabilised form for bad meshes
-        DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
+        nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
     }
 
-    forAll(DeltaCoeffs.boundaryField(), patchi)
+    forAll(nonOrthDeltaCoeffs.boundaryField(), patchi)
     {
-        mesh_.boundary()[patchi].makeDeltaCoeffs
-        (
-            DeltaCoeffs.boundaryField()[patchi]
-        );
+        vectorField delta = mesh_.boundary()[patchi].delta();
+
+        nonOrthDeltaCoeffs.boundaryField()[patchi] =
+            1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
     }
 }
 
 
-void surfaceInterpolation::makeCorrectionVectors() const
+void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
 {
     if (debug)
     {
-        Info<< "surfaceInterpolation::makeCorrectionVectors() : "
+        Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
             << "Constructing non-orthogonal correction vectors"
             << endl;
     }
 
-    correctionVectors_ = new surfaceVectorField
+    nonOrthCorrectionVectors_ = new surfaceVectorField
     (
         IOobject
         (
-            "correctionVectors",
+            "nonOrthCorrectionVectors",
             mesh_.pointsInstance(),
             mesh_
         ),
         mesh_,
         dimless
     );
-    surfaceVectorField& corrVecs = *correctionVectors_;
+    surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
 
     // Set local references to mesh data
     const volVectorField& C = mesh_.C();
@@ -292,14 +330,14 @@ void surfaceInterpolation::makeCorrectionVectors() const
     const labelUList& neighbour = mesh_.neighbour();
     const surfaceVectorField& Sf = mesh_.Sf();
     const surfaceScalarField& magSf = mesh_.magSf();
-    const surfaceScalarField& DeltaCoeffs = deltaCoeffs();
+    const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
 
     forAll(owner, facei)
     {
         vector unitArea = Sf[facei]/magSf[facei];
         vector delta = C[neighbour[facei]] - C[owner[facei]];
 
-        corrVecs[facei] = unitArea - delta*DeltaCoeffs[facei];
+        corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
     }
 
     // Boundary correction vectors set to zero for boundary patches
@@ -308,18 +346,18 @@ void surfaceInterpolation::makeCorrectionVectors() const
 
     forAll(corrVecs.boundaryField(), patchi)
     {
-        fvsPatchVectorField& patchcorrVecs = corrVecs.boundaryField()[patchi];
+        fvsPatchVectorField& patchCorrVecs = corrVecs.boundaryField()[patchi];
 
-        if (!patchcorrVecs.coupled())
+        if (!patchCorrVecs.coupled())
         {
-            patchcorrVecs = vector::zero;
+            patchCorrVecs = vector::zero;
         }
         else
         {
-            const fvsPatchScalarField& patchDeltaCoeffs
-                = DeltaCoeffs.boundaryField()[patchi];
+            const fvsPatchScalarField& patchNonOrthDeltaCoeffs
+                = NonOrthDeltaCoeffs.boundaryField()[patchi];
 
-            const fvPatch& p = patchcorrVecs.patch();
+            const fvPatch& p = patchCorrVecs.patch();
 
             const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
 
@@ -331,60 +369,19 @@ void surfaceInterpolation::makeCorrectionVectors() const
 
                 const vector& delta = patchDeltas[patchFacei];
 
-                patchcorrVecs[patchFacei] =
-                    unitArea - delta*patchDeltaCoeffs[patchFacei];
+                patchCorrVecs[patchFacei] =
+                    unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
             }
         }
     }
 
-    scalar NonOrthogCoeff = 0.0;
-
-    // Calculate the non-orthogonality for meshes with 1 face or more
-    if (returnReduce(magSf.size(), sumOp<label>()) > 0)
-    {
-        NonOrthogCoeff = radToDeg
-        (
-            asin
-            (
-                min
-                (
-                    (sum(magSf*mag(corrVecs))/sum(magSf)).value(),
-                    1.0
-                )
-            )
-        );
-    }
-
     if (debug)
     {
-        Info<< "surfaceInterpolation::makeCorrectionVectors() : "
-            << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
-            << endl;
-    }
-
-    //NonOrthogCoeff = 0.0;
-
-    if (NonOrthogCoeff < 0.1)
-    {
-        orthogonal_ = true;
-        deleteDemandDrivenData(correctionVectors_);
-    }
-    else
-    {
-        orthogonal_ = false;
-    }
-
-    if (debug)
-    {
-        Info<< "surfaceInterpolation::makeCorrectionVectors() : "
+        Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
             << "Finished constructing non-orthogonal correction vectors"
             << endl;
     }
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H

@@ -59,17 +59,17 @@ class surfaceInterpolation
 
         // Demand-driven data
 
-            //- Central-differencing weighting factors
-            mutable surfaceScalarField* weightingFactors_;
+            //- Linear difference weighting factors
+            mutable surfaceScalarField* weights_;
 
-            //- Face-gradient difference factors
-            mutable surfaceScalarField* differenceFactors_;
+            //- Cell-centre difference coefficients
+            mutable surfaceScalarField* deltaCoeffs_;
 
-            //- Is mesh orthogonal
-            mutable bool orthogonal_;
+            //- Non-orthogonal cell-centre difference coefficients
+            mutable surfaceScalarField* nonOrthDeltaCoeffs_;
 
             //- Non-orthogonality correction vectors
-            mutable surfaceVectorField* correctionVectors_;
+            mutable surfaceVectorField* nonOrthCorrectionVectors_;
 
 
     // Private Member Functions
@@ -80,8 +80,11 @@ class surfaceInterpolation
         //- Construct face-gradient difference factors
         void makeDeltaCoeffs() const;
 
+        //- Construct face-gradient difference factors
+        void makeNonOrthDeltaCoeffs() const;
+
         //- Construct non-orthogonality correction vectors
-        void makeCorrectionVectors() const;
+        void makeNonOrthCorrectionVectors() const;
 
 
 protected:
@@ -112,17 +115,18 @@ public:
 
     // Member functions
 
-        //- Return reference to weighting factors array
+        //- Return reference to linear difference weighting factors
         const surfaceScalarField& weights() const;
 
-        //- Return reference to difference factors array
+        //- Return reference to cell-centre difference coefficients
         const surfaceScalarField& deltaCoeffs() const;
 
-        //- Return whether mesh is orthogonal or not
-        bool orthogonal() const;
+        //- Return reference to non-orthogonal cell-centre difference
+        //  coefficients
+        const surfaceScalarField& nonOrthDeltaCoeffs() const;
 
-        //- Return reference to non-orthogonality correction vectors array
-        const surfaceVectorField& correctionVectors() const;
+        //- Return reference to non-orthogonality correction vectors
+        const surfaceVectorField& nonOrthCorrectionVectors() const;
 
         //- Do what is neccessary if the mesh has moved
         bool movePoints();

src/fvMotionSolver/Make/files

@@ -33,5 +33,6 @@ pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPat
 pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
 pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
 pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
+pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C
 
 LIB = $(FOAM_LIBBIN)/libfvMotionSolvers

src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.C

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "waveDisplacementPointPatchVectorField.H"
+#include "pointPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+waveDisplacementPointPatchVectorField::
+waveDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchField<vector>(p, iF),
+    amplitude_(vector::zero),
+    omega_(0.0),
+    waveNumber_(vector::zero)
+{}
+
+
+waveDisplacementPointPatchVectorField::
+waveDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValuePointPatchField<vector>(p, iF, dict),
+    amplitude_(dict.lookup("amplitude")),
+    omega_(readScalar(dict.lookup("omega"))),
+    waveNumber_(dict.lookupOrDefault<vector>("waveLength", vector::zero))
+{
+    if (!dict.found("value"))
+    {
+        updateCoeffs();
+    }
+}
+
+
+waveDisplacementPointPatchVectorField::
+waveDisplacementPointPatchVectorField
+(
+    const waveDisplacementPointPatchVectorField& ptf,
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const pointPatchFieldMapper& mapper
+)
+:
+    fixedValuePointPatchField<vector>(ptf, p, iF, mapper),
+    amplitude_(ptf.amplitude_),
+    omega_(ptf.omega_),
+    waveNumber_(ptf.waveNumber_)
+{}
+
+
+waveDisplacementPointPatchVectorField::
+waveDisplacementPointPatchVectorField
+(
+    const waveDisplacementPointPatchVectorField& ptf,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchField<vector>(ptf, iF),
+    amplitude_(ptf.amplitude_),
+    omega_(ptf.omega_),
+    waveNumber_(ptf.waveNumber_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void waveDisplacementPointPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const polyMesh& mesh = this->dimensionedInternalField().mesh()();
+    const Time& t = mesh.time();
+
+    const scalarField points( waveNumber_ & patch().localPoints());
+
+    Field<vector>::operator=
+    (
+        amplitude_*cos(omega_*t.value() - points)
+    );
+
+    fixedValuePointPatchField<vector>::updateCoeffs();
+}
+
+
+void waveDisplacementPointPatchVectorField::write(Ostream& os) const
+{
+    pointPatchField<vector>::write(os);
+    os.writeKeyword("amplitude")
+        << amplitude_ << token::END_STATEMENT << nl;
+    os.writeKeyword("omega")
+        << omega_ << token::END_STATEMENT << nl;
+    os.writeKeyword("waveNumber")
+        << waveNumber_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePointPatchTypeField
+(
+    pointPatchVectorField,
+    waveDisplacementPointPatchVectorField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/waveDisplacement/waveDisplacementPointPatchVectorField.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::waveDisplacementPointPatchVectorField
+
+Description
+    Foam::waveDisplacementPointPatchVectorField
+
+SourceFiles
+    waveDisplacementPointPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waveDisplacementPointPatchVectorField_H
+#define waveDisplacementPointPatchVectorField_H
+
+#include "fixedValuePointPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class waveDisplacementPointPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class waveDisplacementPointPatchVectorField
+:
+    public fixedValuePointPatchField<vector>
+{
+    // Private data
+
+        vector amplitude_;
+        scalar omega_;
+        vector waveNumber_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waveDisplacement");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        waveDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        waveDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given patchField<vector> onto a new patch
+        waveDisplacementPointPatchVectorField
+        (
+            const waveDisplacementPointPatchVectorField&,
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const pointPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<pointPatchField<vector> > clone() const
+        {
+            return autoPtr<pointPatchField<vector> >
+            (
+                new waveDisplacementPointPatchVectorField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        waveDisplacementPointPatchVectorField
+        (
+            const waveDisplacementPointPatchVectorField&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual autoPtr<pointPatchField<vector> > clone
+        (
+            const DimensionedField<vector, pointMesh>& iF
+        ) const
+        {
+            return autoPtr<pointPatchField<vector> >
+            (
+                new waveDisplacementPointPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/Make/files

@@ -30,8 +30,13 @@ meshSearch/meshSearch.C
 
 meshTools/meshTools.C
 
-PointEdgeWave/PointEdgeWaveName.C
-PointEdgeWave/pointEdgePoint.C
+pWave = PointEdgeWave
+$(pWave)/PointEdgeWaveName.C
+$(pWave)/pointEdgePoint.C
+
+patchWave = PatchEdgeFaceWave
+$(patchWave)/PatchEdgeFaceWaveName.C
+$(patchWave)/patchEdgeFaceInfo.C
 
 regionSplit/regionSplit.C
 

src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.C

@@ -0,0 +1,681 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchEdgeFaceWave.H"
+#include "polyMesh.H"
+#include "globalMeshData.H"
+#include "PatchTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::scalar Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+propagationTol_ = 0.01;
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label
+Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+dummyTrackData_ = 12345;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update info for edgeI, at position pt, with information from
+// neighbouring face.
+// Updates:
+//      - changedEdge_, changedEdges_,
+//      - statistics: nEvals_, nUnvisitedEdges_
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+bool Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+updateEdge
+(
+    const label edgeI,
+    const label neighbourFaceI,
+    const Type& neighbourInfo,
+    Type& edgeInfo
+)
+{
+    nEvals_++;
+
+    bool wasValid = edgeInfo.valid(td_);
+
+    bool propagate =
+        edgeInfo.updateEdge
+        (
+            mesh_,
+            patch_,
+            edgeI,
+            neighbourFaceI,
+            neighbourInfo,
+            propagationTol_,
+            td_
+        );
+
+    if (propagate)
+    {
+        if (!changedEdge_[edgeI])
+        {
+            changedEdge_[edgeI] = true;
+            changedEdges_.append(edgeI);
+        }
+    }
+
+    if (!wasValid && edgeInfo.valid(td_))
+    {
+        --nUnvisitedEdges_;
+    }
+
+    return propagate;
+}
+
+
+// Update info for faceI, at position pt, with information from
+// neighbouring edge.
+// Updates:
+//      - changedFace_, changedFaces_,
+//      - statistics: nEvals_, nUnvisitedFace_
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+bool Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+updateFace
+(
+    const label faceI,
+    const label neighbourEdgeI,
+    const Type& neighbourInfo,
+    Type& faceInfo
+)
+{
+    nEvals_++;
+
+    bool wasValid = faceInfo.valid(td_);
+
+    bool propagate =
+        faceInfo.updateFace
+        (
+            mesh_,
+            patch_,
+            faceI,
+            neighbourEdgeI,
+            neighbourInfo,
+            propagationTol_,
+            td_
+        );
+
+    if (propagate)
+    {
+        if (!changedFace_[faceI])
+        {
+            changedFace_[faceI] = true;
+            changedFaces_.append(faceI);
+        }
+    }
+
+    if (!wasValid && faceInfo.valid(td_))
+    {
+        --nUnvisitedFaces_;
+    }
+
+    return propagate;
+}
+
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+void Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+syncEdges()
+{
+    const globalMeshData& globalData = mesh_.globalData();
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+    const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation();
+
+    // Convert patch-edge data into cpp-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //- Construct with all data in consistent orientation
+    List<Type> cppEdgeData(map.constructSize());
+
+    forAll(patchEdges_, i)
+    {
+        label patchEdgeI = patchEdges_[i];
+        label coupledEdgeI = coupledEdges_[i];
+
+        if (changedEdge_[patchEdgeI])
+        {
+            const Type& data = allEdgeInfo_[patchEdgeI];
+
+            // Patch-edge data needs to be converted into coupled-edge data
+            // (optionally flipped) and consistent in orientation with
+            // master of coupled edge (optionally flipped)
+            bool sameOrientation =
+            (
+                sameEdgeOrientation_[i]
+             == cppOrientation[coupledEdgeI]
+            );
+
+            cppEdgeData[coupledEdgeI].updateEdge
+            (
+                mesh_,
+                patch_,
+                data,
+                sameOrientation,
+                propagationTol_,
+                td_
+            );
+        }
+    }
+
+
+    // Synchronise
+    // ~~~~~~~~~~~
+
+    globalData.syncData
+    (
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),
+        map,
+        globalData.globalTransforms(),
+        updateOp<PrimitivePatchType, Type, TrackingData>
+        (
+            mesh_,
+            patch_,
+            propagationTol_,
+            td_
+        ),
+        transformOp<PrimitivePatchType, Type, TrackingData>
+        (
+            mesh_,
+            patch_,
+            propagationTol_,
+            td_
+        )
+    );
+
+
+    // Back from cpp-edge to patch-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(patchEdges_, i)
+    {
+        label patchEdgeI = patchEdges_[i];
+        label coupledEdgeI = coupledEdges_[i];
+
+        const Type& data = cppEdgeData[coupledEdgeI];
+
+        if (data.valid(td_))
+        {
+            bool sameOrientation =
+            (
+                sameEdgeOrientation_[i]
+             == cppOrientation[coupledEdgeI]
+            );
+
+            allEdgeInfo_[patchEdgeI].updateEdge
+            (
+                mesh_,
+                patch_,
+                data,
+                sameOrientation,
+                propagationTol_,
+                td_
+            );
+
+            if (!changedEdge_[patchEdgeI])
+            {
+                changedEdges_.append(patchEdgeI);
+                changedEdge_[patchEdgeI] = true;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Iterate, propagating changedEdgesInfo across patch, until no change (or
+// maxIter reached). Initial edge values specified.
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+PatchEdgeFaceWave
+(
+    const polyMesh& mesh,
+    const PrimitivePatchType& patch,
+    const labelList& changedEdges,
+    const List<Type>& changedEdgesInfo,
+
+    UList<Type>& allEdgeInfo,
+    UList<Type>& allFaceInfo,
+
+    const label maxIter,
+    TrackingData& td
+)
+:
+    mesh_(mesh),
+    patch_(patch),
+    allEdgeInfo_(allEdgeInfo),
+    allFaceInfo_(allFaceInfo),
+    td_(td),
+    changedEdge_(patch_.nEdges()),
+    changedEdges_(patch_.size()),
+    changedFace_(patch_.size()),
+    changedFaces_(patch_.size()),
+    nEvals_(0),
+    nUnvisitedEdges_(patch_.nEdges()),
+    nUnvisitedFaces_(patch_.size())
+{
+    // Calculate addressing between patch_ and mesh.globalData().coupledPatch()
+    // for ease of synchronisation
+    PatchTools::matchEdges
+    (
+        patch_,
+        mesh_.globalData().coupledPatch(),
+
+        patchEdges_,
+        coupledEdges_,
+        sameEdgeOrientation_
+    );
+
+
+    if (allEdgeInfo_.size() != patch_.nEdges())
+    {
+        FatalErrorIn
+        (
+            "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave"
+            "(const polyMesh&, const labelList&, const List<Type>,"
+            " List<Type>&, List<Type>&, const label maxIter)"
+        )   << "size of edgeInfo work array is not equal to the number"
+            << " of edges in the patch" << endl
+            << "    edgeInfo   :" << allEdgeInfo_.size() << endl
+            << "    patch.nEdges:" << patch_.nEdges()
+            << exit(FatalError);
+    }
+    if (allFaceInfo_.size() != patch_.size())
+    {
+        FatalErrorIn
+        (
+            "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave"
+            "(const polyMesh&, const labelList&, const List<Type>,"
+            " List<Type>&, List<Type>&, const label maxIter)"
+        )   << "size of edgeInfo work array is not equal to the number"
+            << " of faces in the patch" << endl
+            << "    faceInfo   :" << allFaceInfo_.size() << endl
+            << "    patch.size:" << patch_.size()
+            << exit(FatalError);
+    }
+
+
+    // Set from initial changed edges data
+    setEdgeInfo(changedEdges, changedEdgesInfo);
+
+    if (debug)
+    {
+        Pout<< "Seed edges                : " << changedEdges_.size() << endl;
+    }
+
+    // Iterate until nothing changes
+    label iter = iterate(maxIter);
+
+    if ((maxIter > 0) && (iter >= maxIter))
+    {
+        FatalErrorIn
+        (
+            "PatchEdgeFaceWave<Type, TrackingData>::PatchEdgeFaceWave"
+            "(const polyMesh&, const labelList&, const List<Type>,"
+            " List<Type>&, List<Type>&, const label maxIter)"
+        )   << "Maximum number of iterations reached. Increase maxIter." << endl
+            << "    maxIter:" << maxIter << endl
+            << "    changedEdges:" << changedEdges_.size() << endl
+            << "    changedFaces:" << changedFaces_.size() << endl
+            << exit(FatalError);
+    }
+}
+
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+PatchEdgeFaceWave
+(
+    const polyMesh& mesh,
+    const PrimitivePatchType& patch,
+    UList<Type>& allEdgeInfo,
+    UList<Type>& allFaceInfo,
+    TrackingData& td
+)
+:
+    mesh_(mesh),
+    patch_(patch),
+    allEdgeInfo_(allEdgeInfo),
+    allFaceInfo_(allFaceInfo),
+    td_(td),
+    changedEdge_(patch_.nEdges()),
+    changedEdges_(patch_.nEdges()),
+    changedFace_(patch_.size()),
+    changedFaces_(patch_.size()),
+    nEvals_(0),
+    nUnvisitedEdges_(patch_.nEdges()),
+    nUnvisitedFaces_(patch_.size())
+{
+    // Calculate addressing between patch_ and mesh.globalData().coupledPatch()
+    // for ease of synchronisation
+    PatchTools::matchEdges
+    (
+        patch_,
+        mesh_.globalData().coupledPatch(),
+
+        patchEdges_,
+        coupledEdges_,
+        sameEdgeOrientation_
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+getUnsetEdges() const
+{
+    return nUnvisitedEdges_;
+}
+
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+getUnsetFaces() const
+{
+    return nUnvisitedFaces_;
+}
+
+
+// Copy edge information into member data
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+void Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+setEdgeInfo
+(
+    const labelList& changedEdges,
+    const List<Type>& changedEdgesInfo
+)
+{
+    forAll(changedEdges, changedEdgeI)
+    {
+        label edgeI = changedEdges[changedEdgeI];
+
+        bool wasValid = allEdgeInfo_[edgeI].valid(td_);
+
+        // Copy info for edgeI
+        allEdgeInfo_[edgeI] = changedEdgesInfo[changedEdgeI];
+
+        // Maintain count of unset edges
+        if (!wasValid && allEdgeInfo_[edgeI].valid(td_))
+        {
+            --nUnvisitedEdges_;
+        }
+
+        // Mark edgeI as changed, both on list and on edge itself.
+
+        if (!changedEdge_[edgeI])
+        {
+            changedEdge_[edgeI] = true;
+            changedEdges_.append(edgeI);
+        }
+    }
+}
+
+
+// Propagate information from face to edge. Return number of edges changed.
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+faceToEdge()
+{
+    changedEdges_.clear();
+    changedEdge_ = false;
+
+    forAll(changedFaces_, changedFaceI)
+    {
+        label faceI = changedFaces_[changedFaceI];
+
+        if (!changedFace_[faceI])
+        {
+            FatalErrorIn("PatchEdgeFaceWave<Type, TrackingData>::faceToEdge()")
+                << "face " << faceI
+                << " not marked as having been changed" << nl
+                << "This might be caused by multiple occurences of the same"
+                << " seed edge." << abort(FatalError);
+        }
+
+        const Type& neighbourWallInfo = allFaceInfo_[faceI];
+
+        // Evaluate all connected edges
+        const labelList& fEdges = patch_.faceEdges()[faceI];
+
+        forAll(fEdges, fEdgeI)
+        {
+            label edgeI = fEdges[fEdgeI];
+
+            Type& currentWallInfo = allEdgeInfo_[edgeI];
+
+            if (!currentWallInfo.equal(neighbourWallInfo, td_))
+            {
+                updateEdge
+                (
+                    edgeI,
+                    faceI,
+                    neighbourWallInfo,
+                    currentWallInfo
+                );
+            }
+        }
+    }
+
+
+    syncEdges();
+
+
+    if (debug)
+    {
+        Pout<< "Changed edges             : " << changedEdges_.size() << endl;
+    }
+
+    return returnReduce(changedEdges_.size(), sumOp<label>());
+}
+
+
+// Propagate information from edge to face. Return number of faces changed.
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+edgeToFace()
+{
+    changedFaces_.clear();
+    changedFace_ = false;
+
+    const labelListList& edgeFaces = patch_.edgeFaces();
+
+    forAll(changedEdges_, changedEdgeI)
+    {
+        label edgeI = changedEdges_[changedEdgeI];
+
+        if (!changedEdge_[edgeI])
+        {
+            FatalErrorIn("PatchEdgeFaceWave<Type, TrackingData>::edgeToFace()")
+                << "edge " << edgeI
+                << " not marked as having been changed" << nl
+                << "This might be caused by multiple occurences of the same"
+                << " seed edge." << abort(FatalError);
+        }
+
+        const Type& neighbourWallInfo = allEdgeInfo_[edgeI];
+
+        // Evaluate all connected faces
+
+        const labelList& eFaces = edgeFaces[edgeI];
+        forAll(eFaces, eFaceI)
+        {
+            label faceI = eFaces[eFaceI];
+
+            Type& currentWallInfo = allFaceInfo_[faceI];
+
+            if (!currentWallInfo.equal(neighbourWallInfo, td_))
+            {
+                updateFace
+                (
+                    faceI,
+                    edgeI,
+                    neighbourWallInfo,
+                    currentWallInfo
+                );
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Changed faces             : " << changedFaces_.size() << endl;
+    }
+
+    return returnReduce(changedFaces_.size(), sumOp<label>());
+}
+
+
+// Iterate
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData
+>
+Foam::label Foam::PatchEdgeFaceWave<PrimitivePatchType, Type, TrackingData>::
+iterate
+(
+    const label maxIter
+)
+{
+    // Make sure coupled edges contain same info
+    syncEdges();
+
+    nEvals_ = 0;
+
+    label iter = 0;
+
+    while (iter < maxIter)
+    {
+        if (debug)
+        {
+            Pout<< "Iteration " << iter << endl;
+        }
+
+        label nFaces = edgeToFace();
+
+        if (debug)
+        {
+            Pout<< "Total changed faces       : " << nFaces << endl;
+        }
+
+        if (nFaces == 0)
+        {
+            break;
+        }
+
+        label nEdges = faceToEdge();
+
+        if (debug)
+        {
+            Pout<< "Total changed edges       : " << nEdges << nl
+                << "Total evaluations         : " << nEvals_ << nl
+                << "Remaining unvisited edges : " << nUnvisitedEdges_ << nl
+                << "Remaining unvisited faces : " << nUnvisitedFaces_ << nl
+                << endl;
+        }
+
+        if (nEdges == 0)
+        {
+            break;
+        }
+
+        iter++;
+    }
+
+    return iter;
+}
+
+
+// ************************************************************************* //

src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWave.H

@@ -0,0 +1,368 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PatchEdgeFaceWave
+
+Description
+    Wave propagation of information along patch. Every iteration
+    information goes through one layer of faces. Templated on information
+    that is transferred.
+
+SourceFiles
+    PatchEdgeFaceWave.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PatchEdgeFaceWave_H
+#define PatchEdgeFaceWave_H
+
+#include "scalarField.H"
+#include "PackedBoolList.H"
+#include "PrimitivePatch.H"
+#include "vectorTensorTransform.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                    Class PatchEdgeFaceWaveName Declaration
+\*---------------------------------------------------------------------------*/
+
+TemplateName(PatchEdgeFaceWave);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class PatchEdgeFaceWave Declaration
+\*---------------------------------------------------------------------------*/
+
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData = int
+>
+class PatchEdgeFaceWave
+:
+    public PatchEdgeFaceWaveName
+{
+  // Private static data
+
+        //- Relative tolerance. Stop propagation if relative changes
+        //  less than this tolerance (responsability for checking this is
+        //  up to Type implementation)
+        static scalar propagationTol_;
+
+        //- Used as default trackdata value to satisfy default template
+        //  argument.
+        static label dummyTrackData_;
+
+
+    // Private data
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- Reference to patch
+        const PrimitivePatchType& patch_;
+
+        //- Wall information for all edges
+        UList<Type>& allEdgeInfo_;
+
+        //- Information on all patch faces
+        UList<Type>& allFaceInfo_;
+
+        //- Additional data to be passed into container
+        TrackingData& td_;
+
+        //- Has edge changed
+        PackedBoolList changedEdge_;
+
+        //- List of changed edges
+        DynamicList<label> changedEdges_;
+
+        //- Has face changed
+        PackedBoolList changedFace_;
+
+        //- List of changed faces
+        DynamicList<label> changedFaces_;
+
+        //- Number of evaluations
+        label nEvals_;
+
+        //- Number of unvisited faces/edges
+        label nUnvisitedEdges_;
+        label nUnvisitedFaces_;
+
+
+        // Addressing between edges of patch_ and globalData.coupledPatch()
+        labelList patchEdges_;
+        labelList coupledEdges_;
+        PackedBoolList sameEdgeOrientation_;
+
+
+    // Private Member Functions
+
+        //- Updates edgeInfo with information from neighbour. Updates all
+        //  statistics.
+        bool updateEdge
+        (
+            const label edgeI,
+            const label neighbourFaceI,
+            const Type& neighbourInfo,
+            Type& edgeInfo
+        );
+
+        //- Updates faceInfo with information from neighbour. Updates all
+        //  statistics.
+        bool updateFace
+        (
+            const label faceI,
+            const label neighbourEdgeI,
+            const Type& neighbourInfo,
+            Type& faceInfo
+        );
+
+        //- Update coupled edges
+        void syncEdges();
+
+        //- Disallow default bitwise copy construct
+        PatchEdgeFaceWave(const PatchEdgeFaceWave&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PatchEdgeFaceWave&);
+
+
+public:
+
+    // Static Functions
+
+        //- Access to tolerance
+        static scalar propagationTol()
+        {
+            return propagationTol_;
+        }
+
+        //- Change tolerance
+        static void setPropagationTol(const scalar tol)
+        {
+            propagationTol_ = tol;
+        }
+
+
+    // Constructors
+
+        //- Construct from patch, list of changed edges with the Type
+        //  for these edges. Gets work arrays to operate on, one of size
+        //  number of patch edges, the other number of patch faces.
+        //  Iterates until nothing changes or maxIter reached.
+        //  (maxIter can be 0)
+        PatchEdgeFaceWave
+        (
+            const polyMesh& mesh,
+            const PrimitivePatchType& patch,
+            const labelList& initialEdges,
+            const List<Type>& initialEdgesInfo,
+            UList<Type>& allEdgeInfo,
+            UList<Type>& allFaceInfo,
+            const label maxIter,
+            TrackingData& td = dummyTrackData_
+        );
+
+        //- Construct from patch. Use setEdgeInfo and iterate() to do
+        //  actual calculation
+        PatchEdgeFaceWave
+        (
+            const polyMesh& mesh,
+            const PrimitivePatchType& patch,
+            UList<Type>& allEdgeInfo,
+            UList<Type>& allFaceInfo,
+            TrackingData& td = dummyTrackData_
+        );
+
+
+    // Member Functions
+
+        //- Access allEdgeInfo
+        UList<Type>& allEdgeInfo() const
+        {
+            return allEdgeInfo_;
+        }
+
+        //- Access allFaceInfo
+        UList<Type>& allFaceInfo() const
+        {
+            return allFaceInfo_;
+        }
+
+        //- Additional data to be passed into container
+        const TrackingData& data() const
+        {
+            return td_;
+        }
+
+        //- Get number of unvisited faces, i.e. faces that were not (yet)
+        //  reached from walking across patch. This can happen from
+        //  - not enough iterations done
+        //  - a disconnected patch
+        //  - a patch without walls in it
+        label getUnsetFaces() const;
+
+        label getUnsetEdges() const;
+
+        //- Copy initial data into allEdgeInfo_
+        void setEdgeInfo
+        (
+            const labelList& changedEdges,
+            const List<Type>& changedEdgesInfo
+        );
+
+        //- Propagate from edge to face. Returns total number of faces
+        //  (over all processors) changed.
+        label edgeToFace();
+
+        //- Propagate from face to edge. Returns total number of edges
+        //  (over all processors) changed.
+        label faceToEdge();
+
+        //- Iterate until no changes or maxIter reached. Returns actual
+        //  number of iterations.
+        label iterate(const label maxIter);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Update operation
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData = int
+>
+class updateOp
+{
+    //- Additional data to be passed into container
+    const polyMesh& mesh_;
+    const PrimitivePatchType& patch_;
+    const scalar tol_;
+    TrackingData& td_;
+
+public:
+    updateOp
+    (
+        const polyMesh& mesh,
+        const PrimitivePatchType& patch,
+        const scalar tol,
+        TrackingData& td
+    )
+    :
+        mesh_(mesh),
+        patch_(patch),
+        tol_(tol),
+        td_(td)
+    {}
+
+    void operator()(Type& x, const Type& y) const
+    {
+        if (y.valid(td_))
+        {
+            x.updateEdge(mesh_, patch_, y, true, tol_, td_);
+        }
+    }
+};
+
+
+//- Transform operation
+template
+<
+    class PrimitivePatchType,
+    class Type,
+    class TrackingData = int
+>
+class transformOp
+{
+    //- Additional data to be passed into container
+    const polyMesh& mesh_;
+    const PrimitivePatchType& patch_;
+    const scalar tol_;
+    TrackingData& td_;
+
+public:
+    transformOp
+    (
+        const polyMesh& mesh,
+        const PrimitivePatchType& patch,
+        const scalar tol,
+        TrackingData& td
+    )
+    :
+        mesh_(mesh),
+        patch_(patch),
+        tol_(tol),
+        td_(td)
+    {}
+
+    void operator()
+    (
+        const vectorTensorTransform& vt,
+        const bool forward,
+        List<Type>& fld
+    ) const
+    {
+        if (forward)
+        {
+            forAll(fld, i)
+            {
+                fld[i].transform(mesh_, patch_, vt.R(), tol_, td_);
+            }
+        }
+        else
+        {
+            forAll(fld, i)
+            {
+                fld[i].transform(mesh_, patch_, vt.R().T(), tol_, td_);
+            }
+        }
+    }
+};
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PatchEdgeFaceWave.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/PatchEdgeFaceWave/PatchEdgeFaceWaveName.C

@@ -23,45 +23,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "patchPointEdgeCirculator.H"
+#include "PatchEdgeFaceWave.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::patchPointEdgeCirculator
-Foam::patchPointEdgeCirculator::endConstIter
-(
-    *reinterpret_cast<primitiveFacePatch*>(0),  // primitiveFacePatch
-    *reinterpret_cast<PackedBoolList*>(0),      // PackedBoolList
-    -1,                                         // edgeID
-    -1,                                         // index
-    -1                                          // pointID
-);
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const InfoProxy<patchPointEdgeCirculator>& ip
-)
-{
-    const patchPointEdgeCirculator& c = ip.t_;
-
-    const pointField& pts = c.patch_.localPoints();
-    const edge& e = c.patch_.edges()[c.edgeID_];
-    label faceI = c.faceID();
-
-    os  << "around point:" << c.pointID_
-        << " coord:" << pts[c.pointID_]
-        << " at edge:" << c.edgeID_
-        << " coord:" << pts[e.otherVertex(c.pointID_)]
-        << " in direction of face:" << faceI;
-
-    if (faceI != -1)
-    {
-        os  << " fc:" << c.patch_.localFaces()[faceI].centre(pts);
-    }
-    return os;
-}
-
+defineTypeNameAndDebug(Foam::PatchEdgeFaceWaveName, 0);
 
 // ************************************************************************* //

src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchEdgeFaceInfo.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::patchEdgeFaceInfo& wDist
+)
+{
+    return os << wDist.origin() << wDist.distSqr();
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::patchEdgeFaceInfo& wDist
+)
+{
+    return is >> wDist.origin_ >> wDist.distSqr_;
+}
+
+
+// ************************************************************************* //

src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfo.H

@@ -0,0 +1,212 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::patchEdgeFaceInfo
+
+Description
+
+SourceFiles
+    patchEdgeFaceInfoI.H
+    patchEdgeFaceInfo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchEdgeFaceInfo_H
+#define patchEdgeFaceInfo_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "pTraits.H"
+#include "primitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class patchEdgeFaceInfo Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchEdgeFaceInfo
+{
+    // Private data
+
+        //- position of nearest wall center
+        point origin_;
+
+        //- normal distance (squared) from point to origin
+        scalar distSqr_;
+
+
+    // Private Member Functions
+
+        //- Evaluate distance to point. Update distSqr, origin from whomever
+        //  is nearer pt. Return true if w2 is closer to point,
+        //  false otherwise.
+        template<class TrackingData>
+        inline bool update
+        (
+            const point&,
+            const patchEdgeFaceInfo& w2,
+            const scalar tol,
+            TrackingData& td
+        );
+
+        //- Combine current with w2. Update distSqr, origin if w2 has smaller
+        //  quantities and returns true.
+        template<class TrackingData>
+        inline bool update
+        (
+            const patchEdgeFaceInfo& w2,
+            const scalar tol,
+            TrackingData& td
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline patchEdgeFaceInfo();
+
+        //- Construct from origin, distance
+        inline patchEdgeFaceInfo(const point&, const scalar);
+
+        //- Construct as copy
+        inline patchEdgeFaceInfo(const patchEdgeFaceInfo&);
+
+
+    // Member Functions
+
+        // Access
+
+            inline const point& origin() const;
+
+            inline scalar distSqr() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Apply rotation matrix
+            template<class TrackingData>
+            inline void transform
+            (
+                const polyMesh& mesh,
+                const primitivePatch& patch,
+                const tensor& rotTensor,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of face on edge
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const primitivePatch& patch,
+                const label edgeI,
+                const label faceI,
+                const patchEdgeFaceInfo& faceInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- New information for edge (from e.g. coupled edge)
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const primitivePatch& patch,
+                const patchEdgeFaceInfo& edgeInfo,
+                const bool sameOrientation,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of edge on face.
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh& mesh,
+                const primitivePatch& patch,
+                const label faceI,
+                const label edgeI,
+                const patchEdgeFaceInfo& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const patchEdgeFaceInfo&, TrackingData& td) const;
+
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const patchEdgeFaceInfo&) const;
+        inline bool operator!=(const patchEdgeFaceInfo&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const patchEdgeFaceInfo&);
+        friend Istream& operator>>(Istream&, patchEdgeFaceInfo&);
+};
+
+
+//- Data associated with patchEdgeFaceInfo type are contiguous
+template<>
+inline bool contiguous<patchEdgeFaceInfo>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "patchEdgeFaceInfoI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/PatchEdgeFaceWave/patchEdgeFaceInfoI.H

@@ -0,0 +1,268 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2 if w2 nearer to pt.
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::update
+(
+    const point& pt,
+    const patchEdgeFaceInfo& w2,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    scalar dist2 = magSqr(pt - w2.origin());
+
+    if (!valid(td))
+    {
+        // current not yet set so use any value
+        distSqr_ = dist2;
+        origin_ = w2.origin();
+
+        return true;
+    }
+
+    scalar diff = distSqr_ - dist2;
+
+    if (diff < 0)
+    {
+        // already nearer to pt
+        return false;
+    }
+
+    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
+    {
+        // don't propagate small changes
+        return false;
+    }
+    else
+    {
+        // update with new values
+        distSqr_ = dist2;
+        origin_ = w2.origin();
+
+        return true;
+    }
+}
+
+
+// Update this with w2 (information on same edge)
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::update
+(
+    const patchEdgeFaceInfo& w2,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!valid(td))
+    {
+        // current not yet set so use any value
+        distSqr_ = w2.distSqr();
+        origin_ = w2.origin();
+
+        return true;
+    }
+
+    scalar diff = distSqr_ - w2.distSqr();
+
+    if (diff < 0)
+    {
+        // already nearer to pt
+        return false;
+    }
+
+    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
+    {
+        // don't propagate small changes
+        return false;
+    }
+    else
+    {
+        // update with new values
+        distSqr_ =  w2.distSqr();
+        origin_ = w2.origin();
+
+        return true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo()
+:
+    origin_(point::max),
+    distSqr_(sqr(GREAT))
+{}
+
+
+// Construct from origin, distance
+inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo
+(
+    const point& origin,
+    const scalar distSqr
+)
+:
+    origin_(origin),
+    distSqr_(distSqr)
+{}
+
+
+// Construct as copy
+inline Foam::patchEdgeFaceInfo::patchEdgeFaceInfo(const patchEdgeFaceInfo& wpt)
+:
+    origin_(wpt.origin()),
+    distSqr_(wpt.distSqr())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::point& Foam::patchEdgeFaceInfo::origin() const
+{
+    return origin_;
+}
+
+
+inline Foam::scalar Foam::patchEdgeFaceInfo::distSqr() const
+{
+    return distSqr_;
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::valid(TrackingData& td) const
+{
+    return origin_ != point::max;
+}
+
+
+template<class TrackingData>
+inline void Foam::patchEdgeFaceInfo::transform
+(
+    const polyMesh& mesh,
+    const primitivePatch& patch,
+    const tensor& rotTensor,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    origin_ = Foam::transform(rotTensor, origin_);
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::updateEdge
+(
+    const polyMesh& mesh,
+    const primitivePatch& patch,
+    const label edgeI,
+    const label faceI,
+    const patchEdgeFaceInfo& faceInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const edge& e = patch.edges()[edgeI];
+    point eMid =
+        0.5
+      * (
+            patch.points()[patch.meshPoints()[e[0]]]
+          + patch.points()[patch.meshPoints()[e[1]]]
+        );
+    return update(eMid, faceInfo, tol, td);
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::updateEdge
+(
+    const polyMesh& mesh,
+    const primitivePatch& patch,
+    const patchEdgeFaceInfo& edgeInfo,
+    const bool sameOrientation,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(edgeInfo, tol, td);
+}
+
+
+template<class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::updateFace
+(
+    const polyMesh& mesh,
+    const primitivePatch& patch,
+    const label faceI,
+    const label edgeI,
+    const patchEdgeFaceInfo& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(patch.faceCentres()[faceI], edgeInfo, tol, td);
+}
+
+
+template <class TrackingData>
+inline bool Foam::patchEdgeFaceInfo::equal
+(
+    const patchEdgeFaceInfo& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::patchEdgeFaceInfo::operator==
+(
+    const Foam::patchEdgeFaceInfo& rhs
+) const
+{
+    return origin() == rhs.origin();
+}
+
+
+inline bool Foam::patchEdgeFaceInfo::operator!=
+(
+    const Foam::patchEdgeFaceInfo& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/streamLine/controlDict

@@ -0,0 +1,99 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.0.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    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 off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+functions
+{
+    streamLines
+    {
+        type            streamLine;
+
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libfieldFunctionObjects.so");
+
+        // Output every
+        outputControl   outputTime;
+        // outputInterval 10;
+
+        setFormat       vtk; //gnuplot; //xmgr; //raw; //jplot;
+
+        // Velocity field to use for tracking.
+        UName U;
+
+        // Interpolation method. Default is cellPoint. See sampleDict.
+        //interpolationScheme pointMVC;
+
+        // Tracked forwards (+U) or backwards (-U)
+        trackForward    true;
+
+        // Names of fields to sample. Should contain above velocity field!
+        fields (p U);
+
+        // Steps particles can travel before being removed
+        lifeTime        10000;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
+
+        // Cloud name to use
+        cloudName       particleTracks;
+
+        // Seeding method. See the sampleSets in sampleDict.
+        seedSampleSet   uniform;  //cloud;//triSurfaceMeshPointSet;
+
+        uniformCoeffs
+        {
+            type        uniform;
+            axis        x;  //distance;
+
+            start       (-0.0205 0.0001 0.00001);
+            end         (-0.0205 0.0005 0.00001);
+            nPoints     100;
+        }
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -33,6 +33,7 @@ License
 #include "sampledSet.H"
 #include "globalIndex.H"
 #include "mapDistribute.H"
+#include "interpolationCellPoint.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -75,9 +76,9 @@ void Foam::streamLine::track()
 
     // Read or lookup fields
     PtrList<volScalarField> vsFlds;
-    PtrList<interpolationCellPoint<scalar> > vsInterp;
+    PtrList<interpolation<scalar> > vsInterp;
     PtrList<volVectorField> vvFlds;
-    PtrList<interpolationCellPoint<vector> > vvInterp;
+    PtrList<interpolation<vector> > vvInterp;
 
     label UIndex = -1;
 
@@ -95,13 +96,29 @@ void Foam::streamLine::track()
         vsInterp.setSize(vsFlds.size());
         forAll(vsFlds, i)
         {
-            vsInterp.set(i, new interpolationCellPoint<scalar>(vsFlds[i]));
+            vsInterp.set
+            (
+                i,
+                interpolation<scalar>::New
+                (
+                    interpolationScheme_,
+                    vsFlds[i]
+                )
+            );
         }
         ReadFields(mesh, objects, vvFlds);
         vvInterp.setSize(vvFlds.size());
         forAll(vvFlds, i)
         {
-            vvInterp.set(i, new interpolationCellPoint<vector>(vvFlds[i]));
+            vvInterp.set
+            (
+                i,
+                interpolation<vector>::New
+                (
+                    interpolationScheme_,
+                    vvFlds[i]
+                )
+            );
         }
     }
     else
@@ -146,7 +163,12 @@ void Foam::streamLine::track()
                 vsInterp.set
                 (
                     nScalar++,
-                    new interpolationCellPoint<scalar>(f)
+                    //new interpolationCellPoint<scalar>(f)
+                    interpolation<scalar>::New
+                    (
+                        interpolationScheme_,
+                        f
+                    )
                 );
             }
             else if (mesh.foundObject<volVectorField>(fields_[i]))
@@ -164,7 +186,12 @@ void Foam::streamLine::track()
                 vvInterp.set
                 (
                     nVector++,
-                    new interpolationCellPoint<vector>(f)
+                    //new interpolationCellPoint<vector>(f)
+                    interpolation<vector>::New
+                    (
+                        interpolationScheme_,
+                        f
+                    )
                 );
             }
         }
@@ -291,7 +318,22 @@ void Foam::streamLine::read(const dictionary& dict)
     {
         //dict_ = dict;
         dict.lookup("fields") >> fields_;
-        UName_ = dict.lookupOrDefault<word>("U", "U");
+        if (dict.found("UName"))
+        {
+            dict.lookup("UName") >> UName_;
+        }
+        else
+        {
+            UName_ = "U";
+            if (dict.found("U"))
+            {
+                IOWarningIn("streamLine::read(const dictionary&)", dict)
+                    << "Using deprecated entry \"U\"."
+                    << " Please use \"UName\" instead."
+                    << endl;
+                dict.lookup("U") >> UName_;
+            }
+        }
         dict.lookup("trackForward") >> trackForward_;
         dict.lookup("lifeTime") >> lifeTime_;
         if (lifeTime_ < 1)
@@ -307,6 +349,15 @@ void Foam::streamLine::read(const dictionary& dict)
             nSubCycle_ = 1;
         }
 
+        interpolationScheme_ = dict.lookupOrDefault
+        (
+            "interpolationScheme",
+            interpolationCellPoint<scalar>::typeName
+        );
+
+        //Info<< typeName << " using interpolation " << interpolationScheme_
+        //    << endl;
+
         cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine");
         dict.lookup("seedSampleSet") >> seedSet_;
 

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

@@ -86,6 +86,9 @@ class streamLine
         //- Field to transport particle with
         word UName_;
 
+        //- Interpolation scheme to use
+        word interpolationScheme_;
+
         //- Whether to use +u or -u
         bool trackForward_;
 

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

@@ -38,7 +38,7 @@ SourceFiles
 
 #include "particle.H"
 #include "autoPtr.H"
-#include "interpolationCellPoint.H"
+#include "interpolation.H"
 #include "vectorList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,8 +68,8 @@ public:
     public:
 
 
-        const PtrList<interpolationCellPoint<scalar> >& vsInterp_;
-        const PtrList<interpolationCellPoint<vector> >& vvInterp_;
+        const PtrList<interpolation<scalar> >& vsInterp_;
+        const PtrList<interpolation<vector> >& vvInterp_;
         const label UIndex_;
         const bool trackForward_;
         const label nSubCycle_;
@@ -84,8 +84,8 @@ public:
             trackingData
             (
                 Cloud<streamLineParticle>& cloud,
-                const PtrList<interpolationCellPoint<scalar> >& vsInterp,
-                const PtrList<interpolationCellPoint<vector> >& vvInterp,
+                const PtrList<interpolation<scalar> >& vsInterp,
+                const PtrList<interpolation<vector> >& vvInterp,
                 const label UIndex,
                 const bool trackForward,
                 const label nSubCycle,

src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModelCollection.C

@@ -26,6 +26,16 @@ License
 #include "pyrolysisModelCollection.H"
 #include "volFields.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTemplateTypeNameAndDebug
+(
+    Foam::IOPtrList<Foam::regionModels::pyrolysisModels::pyrolysisModel>,
+    0
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 namespace Foam
 {
 namespace regionModels
@@ -33,11 +43,6 @@ namespace regionModels
 namespace pyrolysisModels
 {
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTemplateTypeNameAndDebug(IOPtrList<pyrolysisModel>, 0);
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 pyrolysisModelCollection::pyrolysisModelCollection

src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C

@@ -49,7 +49,8 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
 template<class ThermoType>
 void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
 {
-    volScalarField Yt("Yt", Y_[0]);
+        // It changes Yt patches to "calculated"
+    volScalarField Yt("Yt", 1.0*Y_[0]);
 
     for (label n=1; n<Y_.size(); n++)
     {

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C

@@ -86,7 +86,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     z_(dict.lookup("z")),
     z0_(readScalar(dict.lookup("z0"))),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    zGround_(readScalar(dict.lookup("zGround")))
+    zGround_("zGround", dict, p.size())
 {
     if (mag(z_) < SMALL)
     {

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H

@@ -90,7 +90,7 @@ class atmBoundaryLayerInletEpsilonFvPatchScalarField
         const scalar kappa_;
 
         //- Minimum corrdinate value in z direction
-        const scalar zGround_;
+        const scalarField zGround_;
 
 
 public:

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C

@@ -115,7 +115,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     n_ /= mag(n_);
     z_ /= mag(z_);
 
-    Ustar_ = kappa_*Uref_/(log((Href_  + z0_)/min(z0_ , 0.001)));
+    Ustar_ = kappa_*Uref_/(log((Href_  + z0_)/max(z0_ , 0.001)));
 
     evaluate();
 }
@@ -150,9 +150,11 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
 
     forAll(coord, i)
     {
-        if ((coord[i] - zGround_) < Href_)
+        if ((coord[i] - zGround_[i]) < Href_)
         {
-            Un[i] = (Ustar_/kappa_)*log((coord[i] - zGround_ + z0_)/z0_);
+            Un[i] =
+                (Ustar_/kappa_)
+              * log((coord[i] - zGround_[i] + z0_)/max(z0_, 0.001));
         }
         else
         {
@@ -181,8 +183,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
         << Uref_ << token::END_STATEMENT << nl;
     os.writeKeyword("Href")
         << Href_ << token::END_STATEMENT << nl;
-    os.writeKeyword("zGround")
-        << zGround_ << token::END_STATEMENT << nl;
+    zGround_.writeEntry("zGround", os) ;
     writeEntry("value", os);
 }
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H

@@ -113,7 +113,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
         const scalar Href_;
 
         //- Minimum corrdinate value in z direction
-        const scalar zGround_;
+        const scalarField zGround_;
 
 
 public:

tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/N2

@@ -21,32 +21,27 @@ boundaryField
 {
     top
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
 
     ground
     {
-        type            zeroGradient;
+        type            calculated;
     }
 
     sides
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
 
     burner
     {
-        type            fixedValue;
-        value           uniform 0;
+        type            calculated;
     }
 
     "(region0_to.*)"
     {
-        type            zeroGradient;
+        type            calculated;
     }
 }
 

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2

@@ -23,24 +23,19 @@ boundaryField
 {
     outlet
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
     sides
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
     base
     {
-        type            zeroGradient;
+        type            calculated;
     }
     inlet
     {
-        type            fixedValue;
-        value           uniform 0;
+        type            calculated;
     }
     frontBack
     {

tutorials/combustion/fireFoam/les/smallPoolFire3D/0/N2

@@ -23,24 +23,19 @@ boundaryField
 {
     outlet
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
     sides
     {
-        type            inletOutlet;
-        inletValue      $internalField;
-        value           $internalField;
+        type            calculated;
     }
     base
     {
-        type            zeroGradient;
+        type            calculated;
     }
     inlet
     {
-        type            fixedValue;
-        value           uniform 0;
+        type            calculated;
     }
 }