Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C

@@ -46,8 +46,8 @@ int main(int argc, char *argv[])
     );
 
     argList::noBanner();
-    argList::noParallel();
     argList::addBoolOption("times", "list available times");
+    argList::addBoolOption("latestTime", "list last time");
     argList::addBoolOption
     (
         "keywords",
@@ -75,6 +75,15 @@ int main(int argc, char *argv[])
             Info<< times[i].name() << endl;
         }
     }
+    else if (args.optionFound("latestTime"))
+    {
+        instantList times
+        (
+            Foam::Time::findTimes(args.rootPath()/args.caseName())
+        );
+
+        Info<< times.last().name() << endl;
+    }
 
     if (args.optionFound("dict"))
     {

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -99,6 +99,40 @@ Usage
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+const labelIOList& procAddressing
+(
+    const PtrList<fvMesh>& procMeshList,
+    const label procI,
+    const word& name,
+    PtrList<labelIOList>& procAddressingList
+)
+{
+    const fvMesh& procMesh = procMeshList[procI];
+
+    if (!procAddressingList.set(procI))
+    {
+        procAddressingList.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    name,
+                    procMesh.facesInstance(),
+                    procMesh.meshSubDir,
+                    procMesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+    }
+    return procAddressingList[procI];
+}
+
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -805,74 +839,29 @@ int main(int argc, char *argv[])
                 }
                 const fvMesh& procMesh = procMeshList[procI];
 
+                const labelIOList& faceProcAddressing = procAddressing
+                (
+                    procMeshList,
+                    procI,
+                    "faceProcAddressing",
+                    faceProcAddressingList
+                );
 
-                if (!faceProcAddressingList.set(procI))
+                const labelIOList& cellProcAddressing = procAddressing
-                {
+                (
-                    faceProcAddressingList.set
+                    procMeshList,
-                    (
+                    procI,
-                        procI,
+                    "cellProcAddressing",
-                        new labelIOList
+                    cellProcAddressingList
-                        (
+                );
-                            IOobject
-                            (
-                                "faceProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& faceProcAddressing =
-                    faceProcAddressingList[procI];
 
-
+                const labelIOList& boundaryProcAddressing = procAddressing
-                if (!cellProcAddressingList.set(procI))
+                (
-                {
+                    procMeshList,
-                    cellProcAddressingList.set
+                    procI,
-                    (
+                    "boundaryProcAddressing",
-                        procI,
+                    boundaryProcAddressingList
-                        new labelIOList
+                );
-                        (
-                            IOobject
-                            (
-                                "cellProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& cellProcAddressing =
-                    cellProcAddressingList[procI];
-
-
-                if (!boundaryProcAddressingList.set(procI))
-                {
-                    boundaryProcAddressingList.set
-                    (
-                        procI,
-                        new labelIOList
-                        (
-                            IOobject
-                            (
-                                "boundaryProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& boundaryProcAddressing =
-                    boundaryProcAddressingList[procI];
 
 
                 // FV fields
@@ -959,27 +948,13 @@ int main(int argc, char *argv[])
                  || pointTensorFields.size()
                 )
                 {
-                    if (!pointProcAddressingList.set(procI))
+                    const labelIOList& pointProcAddressing = procAddressing
-                    {
+                    (
-                        pointProcAddressingList.set
+                        procMeshList,
-                        (
+                        procI,
-                            procI,
+                        "pointProcAddressing",
-                            new labelIOList
+                        pointProcAddressingList
-                            (
+                    );
-                                IOobject
-                                (
-                                    "pointProcAddressing",
-                                    procMesh.facesInstance(),
-                                    procMesh.meshSubDir,
-                                    procMesh,
-                                    IOobject::MUST_READ,
-                                    IOobject::NO_WRITE
-                                )
-                            )
-                        );
-                    }
-                    const labelIOList& pointProcAddressing =
-                        pointProcAddressingList[procI];
 
                     const pointMesh& procPMesh = pointMesh::New(procMesh);
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -38,6 +38,7 @@ License
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
+#include "uniformDimensionedFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -195,6 +196,60 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
     }
 
 
+    autoPtr<labelIOList> cellLevelPtr;
+    {
+        IOobject io
+        (
+            "cellLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            cellLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<labelIOList> pointLevelPtr;
+    {
+        IOobject io
+        (
+            "pointLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            pointLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<uniformDimensionedScalarField> level0EdgePtr;
+    {
+        IOobject io
+        (
+            "level0Edge",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            level0EdgePtr.reset(new uniformDimensionedScalarField(io));
+        }
+    }
+
+
+
     label maxProcCells = 0;
     label totProcFaces = 0;
     label maxProcPatches = 0;
@@ -767,8 +822,28 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
 
         procMesh.write();
 
+        // Write points if pointsInstance differing from facesInstance
+        if (facesInstancePointsPtr_.valid())
+        {
+            pointIOField pointsInstancePoints
+            (
+                IOobject
+                (
+                    "points",
+                    pointsInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                xferMove(procPoints)
+            );
+            pointsInstancePoints.write();
+        }
 
 
+        // Decompose any sets
         if (decomposeSets)
         {
             forAll(cellSets, i)
@@ -813,25 +888,67 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
         }
 
 
-        // Write points if pointsInstance differing from facesInstance
+        // hexRef8 data
-        if (facesInstancePointsPtr_.valid())
+        if (cellLevelPtr.valid())
         {
-            pointIOField pointsInstancePoints
+            labelIOList
             (
                 IOobject
                 (
-                    "points",
+                    cellLevelPtr().name(),
-                    pointsInstance(),
+                    facesInstance(),
                     polyMesh::meshSubDir,
                     procMesh,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE,
+                    IOobject::AUTO_WRITE
-                    false
                 ),
-                xferMove(procPoints)
+                UIndirectList<label>
-            );
+                (
-            pointsInstancePoints.write();
+                    cellLevelPtr(),
+                    procCellAddressing_[procI]
+                )()
+            ).write();
         }
+        if (pointLevelPtr.valid())
+        {
+            labelIOList
+            (
+                IOobject
+                (
+                    pointLevelPtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                UIndirectList<label>
+                (
+                    pointLevelPtr(),
+                    procPointAddressing_[procI]
+                )()
+            ).write();
+        }
+        if (level0EdgePtr.valid())
+        {
+            uniformDimensionedScalarField
+            (
+                IOobject
+                (
+                    level0EdgePtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                level0EdgePtr()
+            ).write();
+        }
+
+
+
+        // Statistics
 
         Info<< endl
             << "Processor " << procI << nl

src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,7 +55,9 @@ timeVaryingUniformFixedValuePointPatchField
 :
     fixedValuePointPatchField<Type>(ptf, p, iF, mapper),
     timeSeries_(ptf.timeSeries_)
-{}
+{
+    updateCoeffs();
+}
 
 
 template<class Type>

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

@@ -251,27 +251,35 @@ public:
         {
         public:
 
-            template<class T>
+            template<class Type>
             void operator()
             (
                 const vectorTensorTransform& vt,
                 const bool forward,
-                List<T>& fld
+                List<Type>& fld
             ) const
             {
-                if (forward)
+                const tensor T(forward ? vt.R() : vt.R().T());
+                transformList(T, fld);
+            }
+
+            template<class Type>
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<Type> >& flds
+            ) const
+            {
+                forAll(flds, i)
                 {
-                    transformList(vt.R(), fld);
+                    operator()(vt, forward, flds[i]);
-                }
-                else
-                {
-                    transformList(vt.R().T(), fld);
                 }
             }
 
             //- Transform patch-based field
-            template<class T>
+            template<class Type>
-            void operator()(const coupledPolyPatch& cpp, UList<T>& fld) const
+            void operator()(const coupledPolyPatch& cpp, UList<Type>& fld) const
             {
                 if (!cpp.parallel())
                 {
@@ -280,8 +288,8 @@ public:
             }
 
             //- Transform sparse field
-            template<class T, template<class> class Container>
+            template<class Type, template<class> class Container>
-            void operator()(const coupledPolyPatch& cpp, Container<T>& map)
+            void operator()(const coupledPolyPatch& cpp, Container<Type>& map)
             const
             {
                 if (!cpp.parallel())
@@ -313,6 +321,18 @@ public:
                     fld = vt.invTransformPosition(pfld);
                 }
             }
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<point> >& flds
+            ) const
+            {
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
+                }
+            }
             //- Transform patch-based field
             void operator()(const coupledPolyPatch& cpp, pointField& fld) const
             {

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

@@ -28,43 +28,6 @@ License
 #include "indirectPrimitivePatch.H"
 #include "globalMeshData.H"
 
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    //- Transformation
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template
@@ -141,7 +104,7 @@ Foam::PatchTools::pointNormals
     (
         transforms,
         pointFaceNormals,
-        listTransform()
+        mapDistribute::transform()
     );
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -219,7 +219,7 @@ Foam::label Foam::removePoints::countPointUsage
 
     pointCanBeDeleted.setSize(mesh_.nPoints());
     pointCanBeDeleted = false;
-    label nDeleted = 0;
+    //label nDeleted = 0;
 
     forAll(edge0, pointI)
     {
@@ -243,14 +243,14 @@ Foam::label Foam::removePoints::countPointUsage
             if ((e0Vec & e1Vec) > minCos)
             {
                 pointCanBeDeleted[pointI] = true;
-                nDeleted++;
+                //nDeleted++;
             }
         }
         else if (edge0[pointI] == -1)
         {
             // point not used at all
             pointCanBeDeleted[pointI] = true;
-            nDeleted++;
+            //nDeleted++;
         }
     }
     edge0.clear();
@@ -300,6 +300,15 @@ Foam::label Foam::removePoints::countPointUsage
         true                // null value
     );
 
+    label nDeleted = 0;
+    forAll(pointCanBeDeleted, pointI)
+    {
+        if (pointCanBeDeleted[pointI])
+        {
+            nDeleted++;
+        }
+    }
+
     return returnReduce(nDeleted, sumOp<label>());
 }
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -33,6 +33,7 @@ Description
 #include "fvMesh.H"
 #include "Time.H"
 #include "OFstream.H"
+#include "OBJstream.H"
 #include "mapPolyMesh.H"
 #include "pointEdgePoint.H"
 #include "PointEdgeWave.H"
@@ -1193,12 +1194,29 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
     const pointField avgCc(avgCellCentres(mesh, pp));
 
-    // Construct rays from localPoints to beyond cell centre
+    // Construct rays through localPoints to beyond cell centre
+    pointField start(pp.nPoints());
     pointField end(pp.nPoints());
     forAll(localPoints, pointI)
     {
         const point& pt = localPoints[pointI];
-        end[pointI] = pt + 2*(avgCc[pointI]-pt);
+        const vector d = 2*(avgCc[pointI]-pt);
+        start[pointI] = pt - d;
+        end[pointI] = pt + d;
+    }
+
+
+    autoPtr<OBJstream> gapStr;
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        gapStr.reset
+        (
+            new OBJstream
+            (
+                mesh.time().path()
+              / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
+            )
+        );
     }
 
 
@@ -1226,7 +1244,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
         surfaces.findNearestIntersection
         (
             unzonedSurfaces,
-            localPoints,
+            start,
             end,
 
             surface1,
@@ -1248,38 +1266,67 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
             bool override = false;
 
-            if (hit1[pointI].hit())
+            //if (hit1[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit1[pointI].hitPoint(),
+            //            normal1[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit1[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+            //if (hit2[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit2[pointI].hitPoint(),
+            //            normal2[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit2[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+
+            if (hit1[pointI].hit() && hit2[pointI].hit())
             {
                 if
                 (
                     meshRefiner_.isGap
                     (
                         planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
                         hit1[pointI].hitPoint(),
-                        normal1[pointI]
+                        normal1[pointI],
-                    )
-                )
-                {
-                    disp[pointI] = hit1[pointI].hitPoint()-pt;
-                    override = true;
-                }
-            }
-            if (hit2[pointI].hit())
-            {
-                if
-                (
-                    meshRefiner_.isGap
-                    (
-                        planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
                         hit2[pointI].hitPoint(),
                         normal2[pointI]
                     )
                 )
                 {
+                    // TBD: check if the attraction (to nearest) would attract
+                    // good enough and not override attraction
+
+                    if (gapStr.valid())
+                    {
+                        const point& intPt = hit2[pointI].hitPoint();
+                        gapStr().write(linePointRef(pt, intPt));
+                    }
+
+                    // Choose hit2 : nearest to end point (so inside the domain)
                     disp[pointI] = hit2[pointI].hitPoint()-pt;
                     override = true;
                 }
@@ -1337,7 +1384,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
             surfaces.findNearestIntersection
             (
                 surfacesToTest,
-                pointField(localPoints, zonePointIndices),
+                pointField(start, zonePointIndices),
                 pointField(end, zonePointIndices),
 
                 surface1,
@@ -1363,38 +1410,63 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
                 bool override = false;
 
-                if (hit1[i].hit())
+                //if (hit1[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit1[i].hitPoint(),
+                //            normal1[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit1[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+                //if (hit2[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit2[i].hitPoint(),
+                //            normal2[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit2[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+
+                if (hit1[i].hit() && hit2[i].hit())
                 {
                     if
                     (
                         meshRefiner_.isGap
                         (
                             planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
                             hit1[i].hitPoint(),
-                            normal1[i]
+                            normal1[i],
-                        )
-                    )
-                    {
-                        disp[pointI] = hit1[i].hitPoint()-pt;
-                        override = true;
-                    }
-                }
-                if (hit2[i].hit())
-                {
-                    if
-                    (
-                        meshRefiner_.isGap
-                        (
-                            planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
                             hit2[i].hitPoint(),
                             normal2[i]
                         )
                     )
                     {
+                        if (gapStr.valid())
+                        {
+                            const point& intPt = hit2[i].hitPoint();
+                            gapStr().write(linePointRef(pt, intPt));
+                        }
+
                         disp[pointI] = hit2[i].hitPoint()-pt;
                         override = true;
                     }

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -623,6 +623,149 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 }
 
 
+// Split faces
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitFaces
+(
+    const labelList& splitFaces,
+    const labelPairList& splits
+)
+{
+    polyTopoChange meshMod(mesh_);
+
+    forAll(splitFaces, i)
+    {
+        label faceI = splitFaces[i];
+        const face& f = mesh_.faces()[faceI];
+
+
+        // Split as start and end index in face
+        const labelPair& split = splits[i];
+
+        label nVerts = split[1]-split[0];
+        if (nVerts < 0)
+        {
+            nVerts += f.size();
+        }
+        nVerts += 1;
+
+
+        // Split into f0, f1
+        face f0(nVerts);
+
+        label fp = split[0];
+        forAll(f0, i)
+        {
+            f0[i] = f[fp];
+            fp = f.fcIndex(fp);
+        }
+
+        face f1(f.size()-f0.size()+2);
+        fp = split[1];
+        forAll(f1, i)
+        {
+            f1[i] = f[fp];
+            fp = f.fcIndex(fp);
+        }
+
+
+        // Determine face properties
+        label own = mesh_.faceOwner()[faceI];
+        label nei = -1;
+        label patchI = -1;
+        if (faceI >= mesh_.nInternalFaces())
+        {
+            patchI = mesh_.boundaryMesh().whichPatch(faceI);
+        }
+        else
+        {
+            nei = mesh_.faceNeighbour()[faceI];
+        }
+
+        label zoneI = mesh_.faceZones().whichZone(faceI);
+        bool zoneFlip = false;
+        if (zoneI != -1)
+        {
+            const faceZone& fz = mesh_.faceZones()[zoneI];
+            zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
+        }
+
+
+Pout<< "face:" << faceI << " verts:" << f
+    << " split into f0:" << f0
+    << " f1:" << f1 << endl;
+
+        // Change/add faces
+        meshMod.modifyFace
+        (
+            f0,                         // modified face
+            faceI,                      // label of face
+            own,                        // owner
+            nei,                        // neighbour
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+        meshMod.addFace
+        (
+            f1,                         // modified face
+            own,                        // owner
+            nei,                        // neighbour
+            -1,                         // master point
+            -1,                         // master edge
+            faceI,                      // master face
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+    }
+
+
+
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+    // Update fields
+    mesh_.updateMesh(map);
+
+    // Move mesh (since morphing might not do this)
+    if (map().hasMotionPoints())
+    {
+        mesh_.movePoints(map().preMotionPoints());
+    }
+    else
+    {
+        // Delete mesh volumes. No other way to do this?
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
+
+    // Update local mesh data
+    const labelList& oldToNew = map().reverseFaceMap();
+    labelList newSplitFaces(renumber(oldToNew, splitFaces));
+    // Add added faces (every splitFaces becomes two faces)
+    label sz = newSplitFaces.size();
+    newSplitFaces.setSize(2*sz);
+    forAll(map().faceMap(), faceI)
+    {
+        label oldFaceI = map().faceMap()[faceI];
+        if (oldToNew[oldFaceI] != faceI)
+        {
+            // Added face
+            newSplitFaces[sz++] = faceI;
+        }
+    }
+
+    updateMesh(map, newSplitFaces);
+
+    return map;
+}
+
+
 //// Determine for multi-processor regions the lowest numbered cell on the
 //// lowest numbered processor.
 //void Foam::meshRefinement::getCoupledRegionMaster

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -866,6 +866,13 @@ public:
                 const point& keepPoint
             );
 
+            //- Split faces into two
+            autoPtr<mapPolyMesh> splitFaces
+            (
+                const labelList& splitFaces,
+                const labelPairList& splits
+            );
+
             //- Update local numbering for mesh redistribution
             void distribute(const mapDistributePolyMesh&);
 
@@ -1004,6 +1011,26 @@ public:
             //- Do any one of above IO functions. flag is combination of
             //  writeFlag values.
             void write(const label flag, const fileName&) const;
+
+            //- Helper: calculate average
+            template<class T>
+            static T gAverage
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterElem,
+                const UList<T>& values
+            );
+
+            //- Helper: calculate average over selected elements
+            template<class T>
+            static T gAverage
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterElem,
+                const labelList& meshPoints,
+                const UList<T>& values
+            );
+
 };
 
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -267,7 +267,12 @@ void Foam::meshRefinement::markFeatureCellLevel
 
 
     // Find all start cells of features. Is done by tracking from keepPoint.
-    Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
+    Cloud<trackedParticle> startPointCloud
+    (
+        mesh_,
+        "startPointCloud",
+        IDLList<trackedParticle>()
+    );
 
 
     // Features are identical on all processors. Number them so we know
@@ -315,7 +320,7 @@ void Foam::meshRefinement::markFeatureCellLevel
                     }
 
                     // Non-manifold point. Create particle.
-                    cloud.addParticle
+                    startPointCloud.addParticle
                     (
                         new trackedParticle
                         (
@@ -359,7 +364,7 @@ void Foam::meshRefinement::markFeatureCellLevel
                     }
 
                     // Non-manifold point. Create particle.
-                    cloud.addParticle
+                    startPointCloud.addParticle
                     (
                         new trackedParticle
                         (
@@ -384,12 +389,13 @@ void Foam::meshRefinement::markFeatureCellLevel
     maxFeatureLevel = labelList(mesh_.nCells(), -1);
 
     // Database to pass into trackedParticle::move
-    trackedParticle::trackingData td(cloud, maxFeatureLevel);
+    trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
 
     // Track all particles to their end position (= starting feature point)
     // Note that the particle might have started on a different processor
     // so this will transfer across nicely until we can start tracking proper.
-    cloud.move(td, GREAT);
+    startPointCloud.move(td, GREAT);
+
 
     // Reset level
     maxFeatureLevel = -1;
@@ -403,8 +409,64 @@ void Foam::meshRefinement::markFeatureCellLevel
         featureEdgeVisited[featI] = 0u;
     }
 
+
+    Cloud<trackedParticle> cloud
+    (
+        mesh_,
+        "featureCloud",
+        IDLList<trackedParticle>()
+    );
+
+    forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
+    {
+        const trackedParticle& startTp = iter();
+
+        label featI = startTp.i();
+        label pointI = startTp.j();
+
+        const featureEdgeMesh& featureMesh = features_[featI];
+        const labelList& pEdges = featureMesh.pointEdges()[pointI];
+
+        // Now shoot particles down all pEdges.
+        forAll(pEdges, pEdgeI)
+        {
+            label edgeI = pEdges[pEdgeI];
+
+            if (featureEdgeVisited[featI].set(edgeI, 1u))
+            {
+                // Unvisited edge. Make the particle go to the other point
+                // on the edge.
+
+                const edge& e = featureMesh.edges()[edgeI];
+                label otherPointI = e.otherVertex(pointI);
+
+                trackedParticle* tp(new trackedParticle(startTp));
+                tp->end() = featureMesh.points()[otherPointI];
+                tp->j() = otherPointI;
+
+                if (debug&meshRefinement::FEATURESEEDS)
+                {
+                    Pout<< "Adding particle for point:" << pointI
+                        << " coord:" << tp->position()
+                        << " feature:" << featI
+                        << " to track to:" << tp->end()
+                        << endl;
+                }
+
+                cloud.addParticle(tp);
+            }
+        }
+    }
+
+    startPointCloud.clear();
+
+
     while (true)
     {
+        // Track all particles to their end position.
+        cloud.move(td, GREAT);
+
+
         label nParticles = 0;
 
         // Make particle follow edge.
@@ -460,10 +522,23 @@ void Foam::meshRefinement::markFeatureCellLevel
         {
             break;
         }
-
-        // Track all particles to their end position.
-        cloud.move(td, GREAT);
     }
+
+
+
+    //if (debug&meshRefinement::FEATURESEEDS)
+    //{
+    //    forAll(maxFeatureLevel, cellI)
+    //    {
+    //        if (maxFeatureLevel[cellI] != -1)
+    //        {
+    //            Pout<< "Feature went through cell:" << cellI
+    //                << " coord:" << mesh_.cellCentres()[cellI]
+    //                << " leve:" << maxFeatureLevel[cellI]
+    //                << endl;
+    //        }
+    //    }
+    //}
 }
 
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C

@@ -57,6 +57,108 @@ template<class T> void meshRefinement::updateList
 }
 
 
+template<class T>
+T meshRefinement::gAverage
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterElem,
+    const UList<T>& values
+)
+{
+    if (values.size() != isMasterElem.size())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::gAverage\n"
+            "(\n"
+            "    const polyMesh&,\n"
+            "    const PackedBoolList& isMasterElem,\n"
+            "    const UList<T>& values\n"
+            ")\n"
+        )   << "Number of elements in list " << values.size()
+            << " does not correspond to number of elements in isMasterElem "
+            << isMasterElem.size()
+            << exit(FatalError);
+    }
+
+    T sum = pTraits<T>::zero;
+    label n = 0;
+
+    forAll(values, i)
+    {
+        if (isMasterElem[i])
+        {
+            sum += values[i];
+            n++;
+        }
+    }
+
+    reduce(sum, sumOp<T>());
+    reduce(n, sumOp<label>());
+
+    if (n > 0)
+    {
+        return sum/n;
+    }
+    else
+    {
+        return pTraits<T>::max;
+    }
+}
+
+
+template<class T>
+T meshRefinement::gAverage
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterElem,
+    const labelList& meshPoints,
+    const UList<T>& values
+)
+{
+    if (values.size() != meshPoints.size())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::gAverage\n"
+            "(\n"
+            "    const polyMesh&,\n"
+            "    const labelList&,\n"
+            "    const PackedBoolList& isMasterElem,\n"
+            "    const UList<T>& values\n"
+            ")\n"
+        )   << "Number of elements in list " << values.size()
+            << " does not correspond to number of elements in meshPoints "
+            << meshPoints.size()
+            << exit(FatalError);
+    }
+
+    T sum = pTraits<T>::zero;
+    label n = 0;
+
+    forAll(values, i)
+    {
+        if (isMasterElem[meshPoints[i]])
+        {
+            sum += values[i];
+            n++;
+        }
+    }
+
+    reduce(sum, sumOp<T>());
+    reduce(n, sumOp<label>());
+
+    if (n > 0)
+    {
+        return sum/n;
+    }
+    else
+    {
+        return pTraits<T>::max;
+    }
+}
+
+
 // Compare two lists over all boundary faces
 template<class T>
 void meshRefinement::testSyncBoundaryFaceList

src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C

@@ -169,7 +169,7 @@ void Foam::trackedParticle::hitProcessorPatch
     trackingData& td
 )
 {
-    // Remove particle
+    // Move to different processor
     td.switchProcessor = true;
 }
 

src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,12 +163,24 @@ public:
             return end_;
         }
 
+        //- transported label
+        label i() const
+        {
+            return i_;
+        }
+
         //- transported label
         label& i()
         {
             return i_;
         }
 
+        //- transported label
+        label j() const
+        {
+            return j_;
+        }
+
         //- transported label
         label& j()
         {

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -135,145 +135,12 @@ Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const word& surfaceName,
-    const samplingSource sampleSource
-)
-:
-    sampledSurface(name, mesh),
-    surface_
-    (
-        IOobject
-        (
-            surfaceName,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                   // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(sampleSource),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    sampledSurface(name, mesh, dict),
-    surface_
-    (
-        IOobject
-        (
-            dict.lookup("surface"),
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                   // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(samplingSourceNames_[dict.lookup("source")]),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName
-)
-:
-    sampledSurface(name, mesh),
-    surface_
-    (
-        IOobject
-        (
-            name,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                  // registry
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        surface
-    ),
-    sampleSource_(samplingSourceNames_[sampleSourceName]),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledTriSurfaceMesh::needsUpdate() const
 {
-    return needsUpdate_;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::expire()
-{
-    // already marked as expired
-    if (needsUpdate_)
-    {
-        return false;
-    }
-
-    sampledSurface::clearGeom();
-    MeshStorage::clear();
-
-    boundaryTreePtr_.clear();
-    sampleElements_.clear();
-    samplePoints_.clear();
-
-    needsUpdate_ = true;
-    return true;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::update()
-{
-    if (!needsUpdate_)
-    {
-        return false;
-    }
-
-
     // Find the cells the triangles of the surface are in.
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
-
-
     List<nearInfo> nearest(fc.size());
 
     // Global numbering for cells/faces - only used to uniquely identify local
@@ -611,6 +478,157 @@ bool Foam::sampledTriSurfaceMesh::update()
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const word& surfaceName,
+    const samplingSource sampleSource
+)
+:
+    sampledSurface(name, mesh),
+    surface_
+    (
+        IOobject
+        (
+            surfaceName,
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(sampleSource),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    surface_
+    (
+        IOobject
+        (
+            dict.lookup("surface"),
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(samplingSourceNames_[dict.lookup("source")]),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const triSurface& surface,
+    const word& sampleSourceName
+)
+:
+    sampledSurface(name, mesh),
+    surface_
+    (
+        IOobject
+        (
+            name,
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                  // registry
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        surface
+    ),
+    sampleSource_(samplingSourceNames_[sampleSourceName]),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledTriSurfaceMesh::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::expire()
+{
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+    MeshStorage::clear();
+
+    boundaryTreePtr_.clear();
+    sampleElements_.clear();
+    samplePoints_.clear();
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    // Mesh search engine, no triangulation of faces.
+    meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
+
+    return update(meshSearcher);
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update(const treeBoundBox& bb)
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    // Mesh search engine on subset, no triangulation of faces.
+    meshSearch meshSearcher(mesh(), bb, polyMesh::FACEPLANES);
+
+    return update(meshSearcher);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
 (
     const volScalarField& vField

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H

@@ -72,6 +72,7 @@ namespace Foam
 {
 
 class treeDataFace;
+class meshSearch;
 
 /*---------------------------------------------------------------------------*\
                        Class sampledTriSurfaceMesh Declaration
@@ -136,6 +137,8 @@ private:
         tmp<Field<Type> >
         interpolateField(const interpolation<Type>&) const;
 
+        bool update(const meshSearch& meshSearcher);
+
 public:
 
     //- Runtime type information
@@ -189,6 +192,10 @@ public:
         //  Do nothing (and return false) if no update was needed
         virtual bool update();
 
+        //- Update the surface using a bound box to limit the searching.
+        //  For direct use, i.e. not through sample.
+        //  Do nothing (and return false) if no update was needed
+        bool update(const treeBoundBox&);
 
         //- Points of surface
         virtual const pointField& points() const

src/thermophysicalModels/laminarFlameSpeed/Make/files

@@ -3,5 +3,6 @@ laminarFlameSpeed/laminarFlameSpeedNew.C
 constant/constant.C
 Gulders/Gulders.C
 GuldersEGR/GuldersEGR.C
+RaviPetersen/RaviPetersen.C
 
 LIB = $(FOAM_LIBBIN)/liblaminarFlameSpeedModels

src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C

@@ -0,0 +1,365 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "RaviPetersen.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+    defineTypeNameAndDebug(RaviPetersen, 0);
+
+    addToRunTimeSelectionTable
+    (
+        laminarFlameSpeed,
+        RaviPetersen,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::RaviPetersen
+(
+    const dictionary& dict,
+    const psiuReactionThermo& ct
+)
+:
+    laminarFlameSpeed(dict, ct),
+    coeffsDict_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
+    pPoints_(coeffsDict_.lookup("pPoints")),
+    EqRPoints_(coeffsDict_.lookup("EqRPoints")),
+    alpha_(coeffsDict_.lookup("alpha")),
+    beta_(coeffsDict_.lookup("beta")),
+    TRef_(readScalar(coeffsDict_.lookup("TRef")))
+{
+    checkPointsMonotonicity("equivalenceRatio", EqRPoints_);
+    checkPointsMonotonicity("pressure", pPoints_);
+    checkCoefficientArrayShape("alpha", alpha_);
+    checkCoefficientArrayShape("beta", beta_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::~RaviPetersen()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity
+(
+    const word& name,
+    const List<scalar>& x
+) const
+{
+    for (label i = 1; i < x.size(); i ++)
+    {
+        if (x[i] <= x[i-1])
+        {
+            FatalIOErrorIn
+            (
+                "laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity"
+                "("
+                    "const word&, "
+                    "const List<scalar>&"
+                ") const",
+                coeffsDict_
+            )   << "Data points for the " << name
+                << " do not increase monotonically" << endl
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape
+(
+    const word& name,
+    const List<List<List<scalar> > >& x
+) const
+{
+    bool ok = true;
+
+    ok &= x.size() == EqRPoints_.size() - 1;
+
+    forAll(x, i)
+    {
+        ok &= x[i].size() == pPoints_.size();
+
+        forAll(x[i], j)
+        {
+            ok &= x[i][j].size() == x[i][0].size();
+        }
+    }
+
+    if (!ok)
+    {
+        FatalIOErrorIn
+        (
+            "laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape"
+            "("
+                "const word&, "
+                "const List<List<List<scalar> > >&"
+            ") const",
+            coeffsDict_
+        )   << "Inconsistent size of " << name << " coefficients array" << endl
+            << exit(FatalIOError);
+    }
+}
+
+
+inline bool Foam::laminarFlameSpeedModels::RaviPetersen::interval
+(
+    const List<scalar>& xPoints,
+    const scalar x,
+    label& xIndex,
+    scalar& xXi,
+    scalar& xLim
+) const
+{
+    if (x < xPoints.first())
+    {
+        xIndex = 0;
+        xXi = 0.0;
+        xLim = xPoints.first();
+        return false;
+    }
+
+    else if (x > xPoints.last())
+    {
+        xIndex = xPoints.size() - 2;
+        xXi = 1.0;
+        xLim = xPoints.last();
+        return false;
+    }
+
+    for (xIndex = 0; x > xPoints[xIndex+1]; xIndex ++)
+    {
+        // increment xIndex until xPoints[xIndex] < x < xPoints[xIndex+1]
+    }
+
+    xXi = (x - xPoints[xIndex])/(xPoints[xIndex+1] - xPoints[xIndex]);
+    xLim = x;
+
+    return true;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::polynomial
+(
+    const List<scalar>& coeffs,
+    const scalar x
+) const
+{
+    scalar xPow = 1.0;
+    scalar y = 0.0;
+    forAll(coeffs, i)
+    {
+        y += coeffs[i]*xPow;
+        xPow *= x;
+    }
+    return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::dPolynomial
+(
+    const List<scalar>& coeffs,
+    const scalar x
+) const
+{
+    scalar xPow = 1.0;
+    scalar y = 0.0;
+    for (label i = 1; i < coeffs.size(); i++)
+    {
+        y += i*coeffs[i]*xPow;
+        xPow *= x;
+    }
+    return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::THatPowB
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar Tu
+) const
+{
+    return pow
+    (
+        Tu/TRef_,                                                     
+        polynomial(beta_[EqRIndex][pIndex],EqR)                       
+    );
+}
+
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationInRange
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar Tu
+) const
+{
+    // standard correlation
+    return
+        polynomial(alpha_[EqRIndex][pIndex],EqR)
+       *THatPowB(EqRIndex, pIndex, EqR, Tu);
+}
+ 
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationOutOfRange
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar EqRLim,
+    const scalar Tu
+) const
+{
+    scalar A = polynomial(alpha_[EqRIndex][pIndex], EqRLim);
+    scalar dA = dPolynomial(alpha_[EqRIndex][pIndex], EqRLim);
+    scalar dB = dPolynomial(beta_[EqRIndex][pIndex], EqRLim);
+    scalar TB = THatPowB(EqRIndex, pIndex, EqRLim, Tu);
+
+    // linear extrapolation from the bounds of the correlation
+    return max(TB*(A + (dA + A*log(Tu/TRef_)*dB)*(EqR - EqRLim)), 0.0);
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::speed
+(
+    const scalar EqR,
+    const scalar p,
+    const scalar Tu
+) const
+{
+    scalar Su = 0, s;
+
+    label EqRIndex, pIndex;
+    scalar EqRXi, pXi;
+    scalar EqRLim, pLim;
+    bool EqRInRange;
+
+    EqRInRange = interval(EqRPoints_, EqR, EqRIndex, EqRXi, EqRLim);
+
+    interval(pPoints_, p, pIndex, pXi, pLim);
+
+    for (label pI = 0; pI < 2; pI ++)
+    {
+        if (EqRInRange)
+        {
+            s = correlationInRange(EqRIndex, pIndex + pI, EqR, Tu);
+        }
+        else
+        {
+            s = correlationOutOfRange(EqRIndex, pIndex + pI, EqR, EqRLim, Tu);
+        }
+
+        Su += (1 - pXi)*s;
+        pXi = 1 - pXi;
+    }
+
+    return Su;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::laminarFlameSpeedModels::RaviPetersen::operator()() const
+{
+    const volScalarField& p = psiuReactionThermo_.p();
+    const volScalarField& Tu = psiuReactionThermo_.Tu();
+
+    volScalarField EqR
+    (
+        IOobject
+        (
+            "EqR",
+            p.time().timeName(),
+            p.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        p.mesh(),
+        dimensionedScalar("EqR", dimless, 0.0)
+    );
+
+    if (psiuReactionThermo_.composition().contains("ft"))
+    {
+        const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
+
+        EqR =
+            dimensionedScalar
+            (
+                psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
+            )*ft/max(1 - ft, SMALL);
+    }
+    else
+    {
+        EqR = equivalenceRatio_;
+    }
+
+    tmp<volScalarField> tSu0
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar("Su0", dimVelocity, 0.0)
+        )
+    );
+
+    volScalarField& Su0 = tSu0();
+
+    forAll(Su0, cellI)
+    {
+        Su0[cellI] = speed(EqR[cellI], p[cellI], Tu[cellI]);
+    }
+
+    return tSu0;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.H

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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::laminarFlameSpeedModels::RaviPetersen
+
+Description
+    Laminar flame speed obtained from Ravi and Petersen's correlation.
+
+    The correlation for the laminar flame speed \f$Su\f$ is of the following
+    form:
+    \f[
+        Su = \left( \sum \alpha_i \phi^i \right) 
+        \left( \frac{T}{T_{ref}} \right)^{\left( \sum \beta_j \phi^j \right)} 
+    \f]
+
+    Where \f$\phi\f$ is the equivalence ratio, and \f$\alpha\f$ and \f$\beta\f$
+    are polynomial coefficients given for a number of pressure and equivalence
+    ratio points.
+
+SourceFiles
+    RaviPetersen.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RaviPetersen_H
+#define RaviPetersen_H
+
+#include "laminarFlameSpeed.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class RaviPetersen Declaration
+\*---------------------------------------------------------------------------*/
+
+class RaviPetersen
+:
+    public laminarFlameSpeed
+{
+    // Private Data
+
+        dictionary coeffsDict_;
+
+        //- Correlation pressure values
+        List<scalar> pPoints_;
+
+        //- Correlation equivalence ratios
+        List<scalar> EqRPoints_;
+
+        //- Correlation alpha coefficients
+        List<List<List<scalar> > > alpha_;
+
+        //- Correlation beta coefficients
+        List<List<List<scalar> > > beta_;
+
+        //- Reference temperature
+        scalar TRef_;
+
+
+    // Private Member Functions
+
+        //- Check that input points are ordered
+        void checkPointsMonotonicity
+        (
+            const word& name,
+            const List<scalar>& x
+        ) const;
+
+        //- Check that the coefficient arrays are of the correct shape
+        void checkCoefficientArrayShape
+        (
+            const word& name,
+            const List<List<List<scalar> > >& x
+        ) const;
+
+        //- Find and interpolate a value in the data point arrays
+        inline bool interval
+        (
+            const List<scalar>& xPoints,
+            const scalar x,
+            label& xIndex,
+            scalar& xXi,
+            scalar& xLim
+        ) const;
+
+        //- Evaluate a polynomial
+        inline scalar polynomial
+        (
+            const List<scalar>& coeffs,
+            const scalar x
+        ) const;
+
+        //- Evaluate a polynomial differential
+        inline scalar dPolynomial
+        (
+            const List<scalar>& coeffs,
+            const scalar x
+        ) const;
+
+        //- Calculate normalised temperature to the power of the B polynomial
+        inline scalar THatPowB
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar Tu
+        ) const;
+
+        //- Return the flame speed within the correlation range
+        inline scalar correlationInRange
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar Tu
+        ) const;
+
+        //- Extrapolate the flame speed correlation outside its range
+        inline scalar correlationOutOfRange
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar EqRLim,
+            const scalar Tu
+        ) const;
+
+        //- Return the laminar flame speed [m/s]
+        inline scalar speed
+        (
+            const scalar EqR,
+            const scalar p,
+            const scalar Tu
+        ) const;
+
+
+        //- Construct as copy (not implemented)
+        RaviPetersen(const RaviPetersen&);
+        void operator=(const RaviPetersen&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("RaviPetersen");
+
+    // Constructors
+
+        //- Construct from dictionary and psiuReactionThermo
+        RaviPetersen
+        (
+            const dictionary&,
+            const psiuReactionThermo&
+        );
+
+
+    //- Destructor
+    virtual ~RaviPetersen();
+
+
+    // Member functions
+
+        //- Return the laminar flame speed [m/s]
+        tmp<volScalarField> operator()() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End laminarFlameSpeedModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/combustion/XiFoam/ras/Allclean

@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 keepCases="moriyoshiHomogeneous"
-loseCases="moriyoshiHomogeneousPart2"
+loseCases="moriyoshiHomogeneousPart2 moriyoshiHomogeneousHydrogen"
 
 for caseName in $keepCases
 do

tutorials/combustion/XiFoam/ras/Allrun

@@ -6,21 +6,25 @@ cd ${0%/*} || exit 1    # run from this directory
 
 setControlDict()
 {
-    controlDict="system/controlDict"
-
     sed \
-    -e s/"\(deltaT[ \t]*\) 5e-06;"/"\1 1e-05;"/g \
+    -e "s/\(deltaT[ \t]*\) 5e-06;/\1 1e-05;/g" \
-    -e s/"\(endTime[ \t]*\) 0.005;"/"\1 0.015;"/g \
+    -e "s/\(endTime[ \t]*\) 0.005;/\1 0.015;/g" \
-    -e s/"\(writeInterval[ \t]*\) 10;"/"\1 50;"/g \
+    -e "s/\(writeInterval[ \t]*\) 10;/\1 50;/g" \
-    $controlDict > temp.$$
+    -i system/controlDict
-    mv temp.$$ $controlDict
 }
 
+setCombustionProperties()
+{
+    sed \
+    -e "s/\(laminarFlameSpeedCorrelation[ \t]*\) Gulders;/\1 RaviPetersen;/g" \
+    -e "s/\(fuel[ \t]*\) Propane;/\1 HydrogenInAir;/g" \
+    -i constant/combustionProperties
+}
 
 # Do moriyoshiHomogeneous
 ( cd moriyoshiHomogeneous && foamRunTutorials )
 
-# Clone case
+# Clone case for second phase
 cloneCase moriyoshiHomogeneous moriyoshiHomogeneousPart2
 
 # Modify and execute
@@ -32,4 +36,19 @@ cloneCase moriyoshiHomogeneous moriyoshiHomogeneousPart2
     runApplication `getApplication`
 )
 
+# Clone case for hydrogen
+cloneCase moriyoshiHomogeneous moriyoshiHomogeneousHydrogen
+
+# Modify and execute
+(
+    cd moriyoshiHomogeneousHydrogen || exit
+
+    setCombustionProperties
+    mv constant/thermophysicalProperties \
+       constant/thermophysicalProperties.propane
+    mv constant/thermophysicalProperties.hydrogen \
+       constant/thermophysicalProperties
+    runApplication `getApplication`
+)
+
 # ----------------------------------------------------------------- end-of-file

tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/constant/combustionProperties

@@ -68,6 +68,36 @@ GuldersCoeffs
     }
 }
 
+RaviPetersenCoeffs
+{
+    HydrogenInAir
+    {
+        TRef            320;
+        pPoints         ( 1.0e05 5.0e05 1.0e06 2.0e06 3.0e06 );
+        EqRPoints       ( 0.5 2.0 5.0 );
+        alpha           ( ( (-0.03   -2.347  9.984  -6.734  1.361)
+                            ( 1.61   -9.708 19.026 -11.117  2.098)
+                            ( 2.329 -12.287 21.317 -11.973  2.207)
+                            ( 2.593 -12.813 20.815 -11.471  2.095)
+                            ( 2.728 -13.164 20.794 -11.418  2.086) )
+                          ( ( 3.558   0.162 -0.247   0.0253 0    )
+                            ( 4.818  -0.872 -0.053   0.0138 0    )
+                            ( 3.789  -0.312 -0.208   0.028  0    )
+                            ( 4.925  -1.841  0.211  -0.0059 0    )
+                            ( 4.505  -1.906  0.259  -0.0105 0    ) ) );
+        beta            ( ( ( 5.07  -6.42  3.87  -0.767)
+                            ( 5.52  -6.73  3.88  -0.728)
+                            ( 5.76  -6.92  3.92  -0.715)
+                            ( 6.02  -7.44  4.37  -0.825)
+                            ( 7.84 -11.55  7.14  -1.399) )
+                          ( ( 1.405  0.053 0.022  0    )
+                            ( 1.091  0.317 0      0    )
+                            ( 1.64  -0.03  0.07   0    )
+                            ( 0.84   0.56  0      0    )
+                            ( 0.81   0.64  0      0    ) ) );
+    }
+}
+
 ignite          yes;
 
 ignitionSites

tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/constant/thermophysicalProperties.hydrogen

@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heheuPsiThermo;
+    mixture         homogeneousMixture;
+    transport       sutherland;
+    thermo          janaf;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          absoluteEnthalpy;
+}
+
+stoichiometricAirFuelMassRatio stoichiometricAirFuelMassRatio [ 0 0 0 0 0 0 0 ] 34.074;
+
+reactants
+{
+    specie
+    {
+        nMoles          24.8095;
+        molWeight       16.0243;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 3.02082 0.00104314 -2.88613e-07 4.20369e-11 -2.37182e-15 -902.964 2.3064 );
+        lowCpCoeffs     ( 2.99138 0.00343493 -8.43792e-06 9.57755e-09 -3.75097e-12 -987.16 1.95123 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+products
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       17.9973;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 2.879 0.00161934 -4.61257e-07 6.41382e-11 -3.3855e-15 -8023.54 4.11691 );
+        lowCpCoeffs     ( 3.3506 0.00176018 -4.28718e-06 5.63372e-09 -2.35948e-12 -8211.42 1.36387 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/Allclean

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm *.obj
+rm -r constant/extendedFeatureEdgeMesh
+rm constant/triSurface/boundaryAndFaceZones.eMesh
+rm constant/polyMesh/boundary
+
+cleanCase
+
+# ----------------------------------------------------------------- end-of-file

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/Allrun.pre

@@ -19,9 +19,11 @@ runApplication collapseEdges -collapseFaces -latestTime -overwrite
 
 runApplication checkMesh -allTopology -allGeometry -latestTime
 
+latestTime=`foamInfoExec -latestTime`
+
 # Move the mesh into polyMesh
 \rm -r constant/polyMesh
-\mv 100/polyMesh constant
+\mv "${latestTime}"/polyMesh constant
 
 # Clean up intermediate meshes
 \rm -r [1-9]*

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/constant/porosityProperties

@@ -33,7 +33,7 @@ porosity1
             coordinateRotation
             {
                 type    axesRotation;
-                e1      (0.70710678 0.70710678 0);
+                e1      (1 0 0);    //(0.70710678 0.70710678 0);
                 e2      (0 0 1);
             }
         }

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/system/controlDict

@@ -14,13 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//DebugSwitches
+application     porousSimpleFoam;
-//{
-//    regIOobject           2;
-//    cellShapeControlMesh  0;
-//}
-
-application     foamyHexMesh;
 
 startFrom       startTime;
 

tutorials/incompressible/porousSimpleFoam/straightDuctImplicit/system/meshDict.geometry

@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-angledDuct.stl
+straightDuct.stl
 {
     type triSurfaceMesh;
     name angledDuct;