diff --git a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
index ac3b733b038228f347dbe2c8651e19a75a3a50fe..ad178b765cbe55713463e268b1f9463747318826 100644
--- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
+++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
@@ -58,11 +58,15 @@ Usage
     - changeDictionary [OPTION]
 
     \param -literalRE \n
-    Do not interpret regular expressions; treat them as any other keyword.
+    Do not interpret regular expressions or patchGroups;
+    treat them as any other keyword.
 
     \param -enableFunctionEntries \n
     By default all dictionary preprocessing of fields is disabled
 
+    \param -disablePatchGroups \n
+    By default all keys are also checked for being patchGroups
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
@@ -84,7 +88,47 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-bool merge(dictionary&, const dictionary&, const bool);
+// Extract groupPatch (= shortcut) info from boundary file info
+HashTable<wordList, word> extractPatchGroups(const dictionary& boundaryDict)
+{
+    HashTable<wordList, word> groupToPatch;
+
+    forAllConstIter(dictionary, boundaryDict, iter)
+    {
+        const word& patchName = iter().keyword();
+        const dictionary& patchDict = iter().dict();
+
+        wordList groups;
+        if (patchDict.readIfPresent("inGroups", groups))
+        {
+            forAll(groups, i)
+            {
+                HashTable<wordList, word>::iterator fndGroup = groupToPatch.find
+                (
+                    groups[i]
+                );
+                if (fndGroup == groupToPatch.end())
+                {
+                    groupToPatch.insert(groups[i], wordList(1, patchName));
+                }
+                else
+                {
+                    fndGroup().append(patchName);
+                }
+            }
+        }
+    }
+    return groupToPatch;
+}
+
+
+bool merge
+(
+    dictionary&,
+    const dictionary&,
+    const bool,
+    const HashTable<wordList, word>&
+);
 
 
 // Add thisEntry to dictionary thisDict.
@@ -93,7 +137,8 @@ bool addEntry
     dictionary& thisDict,
     entry& thisEntry,
     const entry& mergeEntry,
-    const bool literalRE
+    const bool literalRE,
+    const HashTable<wordList, word>& shortcuts
 )
 {
     bool changed = false;
@@ -108,7 +153,8 @@ bool addEntry
             (
                 const_cast<dictionary&>(thisEntry.dict()),
                 mergeEntry.dict(),
-                literalRE
+                literalRE,
+                shortcuts
             )
         )
         {
@@ -135,7 +181,8 @@ bool merge
 (
     dictionary& thisDict,
     const dictionary& mergeDict,
-    const bool literalRE
+    const bool literalRE,
+    const HashTable<wordList, word>& shortcuts
 )
 {
     bool changed = false;
@@ -156,7 +203,7 @@ bool merge
     {
         const keyType& key = mergeIter().keyword();
 
-        if (literalRE || !key.isPattern())
+        if (literalRE || !(key.isPattern() || shortcuts.found(key)))
         {
             entry* entryPtr = thisDict.lookupEntryPtr
             (
@@ -179,7 +226,8 @@ bool merge
                         thisDict,
                        *entryPtr,
                         mergeIter(),
-                        literalRE
+                        literalRE,
+                        shortcuts
                     )
                 )
                 {
@@ -196,44 +244,70 @@ bool merge
     }
 
 
-    // Pass 2. Wildcard matches (if any) on any non-match keys.
+    // Pass 2. Wildcard or shortcut matches (if any) on any non-match keys.
 
     if (!literalRE && thisKeysSet.size() > 0)
     {
+        // Pick up remaining dictionary entries
         wordList thisKeys(thisKeysSet.toc());
 
         forAllConstIter(IDLList<entry>, mergeDict, mergeIter)
         {
             const keyType& key = mergeIter().keyword();
 
+            // List of indices into thisKeys
+            labelList matches;
+
             if (key.isPattern())
             {
-                // Find all matching entries in the original thisDict
-
-                labelList matches = findStrings(key, thisKeys);
+                // Wildcard match
+                matches = findStrings(key, thisKeys);
 
-                forAll(matches, i)
+            }
+            else if (shortcuts.size())
+            {
+                // See if patchGroups expand to valid thisKeys
+                const wordList shortcutNames = shortcuts.toc();
+                labelList indices = findStrings(key, shortcutNames);
+                forAll(indices, i)
                 {
-                    label matchI = matches[i];
+                    const word& name = shortcutNames[indices[i]];
+                    const wordList& keys = shortcuts[name];
+                    forAll(keys, j)
+                    {
+                        label index = findIndex(thisKeys, keys[j]);
+                        if (index != -1)
+                        {
+                            matches.append(index);
+                        }
+                    }
+                }
+            }
 
-                    entry& thisEntry = const_cast<entry&>
-                    (
-                        thisDict.lookupEntry(thisKeys[matchI], false, false)
-                    );
+            // Add all matches
+            forAll(matches, i)
+            {
+                const word& thisKey = thisKeys[matches[i]];
+
+                entry& thisEntry = const_cast<entry&>
+                (
+                    thisDict.lookupEntry(thisKey, false, false)
+                );
 
-                    if
+                if
+                (
+                    addEntry
                     (
-                        addEntry
-                        (
-                            thisDict,
-                            thisEntry,
-                            mergeIter(),
-                            literalRE
-                        )
+                        thisDict,
+                        thisEntry,
+                        mergeIter(),
+                        literalRE,
+                        HashTable<wordList, word>(0)    // no shortcuts
+                                                        // at deeper levels
                     )
-                    {
-                        changed = true;
-                    }
+                )
+                {
+                    changed = true;
                 }
             }
         }
@@ -273,6 +347,12 @@ int main(int argc, char *argv[])
         "enableFunctionEntries",
         "enable expansion of dictionary directives - #include, #codeStream etc"
     );
+    argList::addBoolOption
+    (
+        "disablePatchGroups",
+        "disable matching keys to patch groups"
+    );
+
     #include "addRegionOption.H"
 
     #include "setRootCase.H"
@@ -304,7 +384,6 @@ int main(int argc, char *argv[])
     }
 
     const bool literalRE = args.optionFound("literalRE");
-
     if (literalRE)
     {
         Info<< "Not interpreting any regular expressions (RE)"
@@ -328,6 +407,15 @@ int main(int argc, char *argv[])
     }
 
 
+    const bool disablePatchGroups = args.optionFound("disablePatchGroups");
+    if (disablePatchGroups)
+    {
+        Info<< "Not interpreting any keys in the changeDictionary"
+            << " as patchGroups"
+            << endl;
+    }
+
+
     fileName regionPrefix = "";
     if (regionName != fvMesh::defaultRegion)
     {
@@ -369,6 +457,70 @@ int main(int argc, char *argv[])
         << replaceDicts.toc() << endl;
 
 
+
+    // Always read boundary to get patch groups
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Info<< "Reading polyMesh/boundary file to extract patch names"
+        << endl;
+
+    // Read PtrList of dictionary as dictionary.
+    const word oldTypeName = IOPtrList<entry>::typeName;
+    const_cast<word&>(IOPtrList<entry>::typeName) = word::null;
+    IOPtrList<entry> dictList
+    (
+        IOobject
+        (
+            "boundary",
+            runTime.findInstance
+            (
+                regionPrefix/polyMesh::meshSubDir,
+                "boundary",
+                IOobject::READ_IF_PRESENT
+            ),
+            polyMesh::meshSubDir,
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+    const_cast<word&>(IOPtrList<entry>::typeName) = oldTypeName;
+    // Fake type back to what was in field
+    const_cast<word&>(dictList.type()) = dictList.headerClassName();
+
+    // Temporary convert to dictionary
+    dictionary fieldDict;
+    forAll(dictList, i)
+    {
+        fieldDict.add(dictList[i].keyword(), dictList[i].dict());
+    }
+
+    if (dictList.size())
+    {
+        Info<< "Loaded dictionary " << dictList.name()
+            << " with entries " << fieldDict.toc() << endl;
+    }
+
+    // Extract any patchGroups information (= shortcut for set of
+    // patches)
+    HashTable<wordList, word> patchGroups;
+    if (!disablePatchGroups)
+    {
+        patchGroups = extractPatchGroups(fieldDict);
+        if (patchGroups.size())
+        {
+            Info<< "Extracted patch groups:" << endl;
+            wordList groups(patchGroups.sortedToc());
+            forAll(groups, i)
+            {
+                Info<< "    group " << groups[i] << " with patches "
+                    << patchGroups[groups[i]] << endl;
+            }
+        }
+    }
+
+
     // Every replacement is a dictionary name and a keyword in this
 
     forAllConstIter(dictionary, replaceDicts, fieldIter)
@@ -384,46 +536,12 @@ int main(int argc, char *argv[])
             Info<< "Special handling of " << fieldName
                 << " as polyMesh/boundary file." << endl;
 
-            // Read PtrList of dictionary as dictionary.
-            const word oldTypeName = IOPtrList<entry>::typeName;
-            const_cast<word&>(IOPtrList<entry>::typeName) = word::null;
-            IOPtrList<entry> dictList
-            (
-                IOobject
-                (
-                    fieldName,
-                    runTime.findInstance
-                    (
-                        regionPrefix/polyMesh::meshSubDir,
-                        fieldName
-                    ),
-                    polyMesh::meshSubDir,
-                    mesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                )
-            );
-            const_cast<word&>(IOPtrList<entry>::typeName) = oldTypeName;
-            // Fake type back to what was in field
-            const_cast<word&>(dictList.type()) = dictList.headerClassName();
-
-            // Temporary convert to dictionary
-            dictionary fieldDict;
-            forAll(dictList, i)
-            {
-                fieldDict.add(dictList[i].keyword(), dictList[i].dict());
-            }
-
-            Info<< "Loaded dictionary " << fieldName
-                << " with entries " << fieldDict.toc() << endl;
-
             // Get the replacement dictionary for the field
             const dictionary& replaceDict = fieldIter().dict();
             Info<< "Merging entries from " << replaceDict.toc() << endl;
 
             // Merge the replacements in
-            merge(fieldDict, replaceDict, literalRE);
+            merge(fieldDict, replaceDict, literalRE, patchGroups);
 
             Info<< "fieldDict:" << fieldDict << endl;
 
@@ -492,7 +610,7 @@ int main(int argc, char *argv[])
             Info<< "Merging entries from " << replaceDict.toc() << endl;
 
             // Merge the replacements in
-            merge(fieldDict, replaceDict, literalRE);
+            merge(fieldDict, replaceDict, literalRE, patchGroups);
 
             Info<< "Writing modified fieldDict " << fieldName << endl;
             fieldDict.regIOobject::write();
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.H
index 56fef9354e49b159835194a49a8e38c2080f38b1..c7fdc5c644b854ec2bea35cf6b220485bc55ed4e 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,7 +104,7 @@ class codedFixedValueFvPatchField
     // Private data
 
         //- Dictionary contents for the boundary condition
-        mutable dictionary dict_;
+        const dictionary dict_;
 
         const word redirectType_;
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 65143f5d27db6a8d30eafa25fbfdab3d2749d3ed..d20552f36832c9703baa333fbd5bb4f092dddced 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -110,6 +110,7 @@ void reactingOneDim::updateQr()
     const volScalarField kappaRad_(kappaRad());
 
     // Propagate Qr through 1-D regions
+    label totalFaceId = 0;
     forAll(intCoupledPatchIDs_, i)
     {
         const label patchI = intCoupledPatchIDs_[i];
@@ -121,7 +122,7 @@ void reactingOneDim::updateQr()
         {
             const scalar Qr0 = Qrp[faceI];
             point Cf0 = Cf[faceI];
-            const labelList& cells = boundaryFaceCells_[faceI];
+            const labelList& cells = boundaryFaceCells_[totalFaceId];
             scalar kappaInt = 0.0;
             forAll(cells, k)
             {
@@ -132,6 +133,7 @@ void reactingOneDim::updateQr()
                 Qr_[cellI] = Qr0*exp(-kappaInt);
                 Cf0 = Cf1;
             }
+            totalFaceId ++;
         }
     }
 
@@ -156,6 +158,7 @@ void reactingOneDim::updatePhiGas()
         const volScalarField& HsiGas = tHsiGas();
         const volScalarField& RRiGas = tRRiGas();
 
+        label totalFaceId = 0;
         forAll(intCoupledPatchIDs_, i)
         {
             const label patchI = intCoupledPatchIDs_[i];
@@ -164,7 +167,7 @@ void reactingOneDim::updatePhiGas()
 
             forAll(phiGasp, faceI)
             {
-                const labelList& cells = boundaryFaceCells_[faceI];
+                const labelList& cells = boundaryFaceCells_[totalFaceId];
                 scalar massInt = 0.0;
                 forAllReverse(cells, k)
                 {
@@ -184,6 +187,7 @@ void reactingOneDim::updatePhiGas()
                         << " is : " << massInt
                         << " [kg/s] " << endl;
                 }
+                totalFaceId ++;
             }
         }
         tHsiGas().clear();
@@ -232,12 +236,22 @@ void reactingOneDim::solveContinuity()
         Info<< "reactingOneDim::solveContinuity()" << endl;
     }
 
-    solve
-    (
-        fvm::ddt(rho_)
-     ==
-      - solidChemistry_->RRg()
-    );
+    if (moveMesh_)
+    {
+        const scalarField mass0 = rho_*regionMesh().V();
+
+        fvScalarMatrix rhoEqn
+        (
+            fvm::ddt(rho_)
+         ==
+          - solidChemistry_->RRg()
+        );
+
+        rhoEqn.solve();
+
+        updateMesh(mass0);
+
+    }
 }
 
 
@@ -261,14 +275,15 @@ void reactingOneDim::solveSpeciesMass()
             solidChemistry_->RRs(i)
         );
 
-        if (moveMesh_)
+        if (regionMesh().moving())
         {
-            surfaceScalarField phiRhoMesh
+            surfaceScalarField phiYiRhoMesh
             (
                 fvc::interpolate(Yi*rho_)*regionMesh().phi()
             );
 
-            YiEqn -= fvc::div(phiRhoMesh);
+            YiEqn += fvc::div(phiYiRhoMesh);
+
         }
 
         YiEqn.solve(regionMesh().solver("Yi"));
@@ -303,14 +318,14 @@ void reactingOneDim::solveEnergy()
       + fvc::div(phiGas)
     );
 
-    if (moveMesh_)
+    if (regionMesh().moving())
     {
-        surfaceScalarField phiMesh
+        surfaceScalarField phihMesh
         (
             fvc::interpolate(rho_*h_)*regionMesh().phi()
         );
 
-        hEqn -= fvc::div(phiMesh);
+        hEqn += fvc::div(phihMesh);
     }
 
     hEqn.relax();
@@ -349,7 +364,18 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     pyrolysisModel(modelType, mesh),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solid()),
-    rho_(solidThermo_.rho()),
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            regionMesh().time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        solidThermo_.rho()
+    ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
     primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
@@ -449,7 +475,18 @@ reactingOneDim::reactingOneDim
     pyrolysisModel(modelType, mesh, dict),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solid()),
-    rho_(solidThermo_.rho()),
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            regionMesh().time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        solidThermo_.rho()
+    ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
     primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
@@ -681,17 +718,15 @@ void reactingOneDim::evolveRegion()
 {
     Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
 
-    const scalarField mass0 = rho_*regionMesh().V();
-
     solidChemistry_->solve
     (
         time().value() - time().deltaTValue(),
         time().deltaTValue()
     );
 
-    solveContinuity();
+    calculateMassTransfer();
 
-    updateMesh(mass0);
+    solveContinuity();
 
     chemistrySh_ = solidChemistry_->Sh()();
 
@@ -704,9 +739,9 @@ void reactingOneDim::evolveRegion()
         solveEnergy();
     }
 
-    calculateMassTransfer();
-
     solidThermo_.correct();
+
+    rho_ = solidThermo_.rho();
 }
 
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
index aa3d7baa9b6c8c9760cf32315e876fbb7829ec82..619c9c0e9aba245e105afd17199c3530db037762 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
@@ -84,14 +84,8 @@ protected:
 
         // Reference to solid thermo properties
 
-//             //- Absorption coefficient [1/m]
-//             const volScalarField& kappaRad_;
-//
-//             //- Thermal conductivity [W/m/K]
-//             const volScalarField& kappa_;
-
             //- Density [kg/m3]
-            volScalarField& rho_;
+            volScalarField rho_;
 
             //- List of solid components
             PtrList<volScalarField>& Ys_;
@@ -221,8 +215,8 @@ public:
 
             //- Fields
 
-                //- Return density [kg/m3]
-                virtual const volScalarField& rho() const;
+                //- Return const density [Kg/m3]
+                const volScalarField& rho() const;
 
                 //- Return const temperature [K]
                 virtual const volScalarField& T() const;
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.C b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
index 372e19a3f35c3859ec70faef5b519aa96e5b202f..60243b55cdc032d6bea6d459f42d71d10ecca6cb 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1D.C
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
@@ -117,6 +117,8 @@ void Foam::regionModels::regionModel1D::initialise()
     }
 
     boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI);
+    boundaryFaceFaces_.setSize(localPyrolysisFaceI);
+    boundaryFaceCells_.setSize(localPyrolysisFaceI);
 
     surfaceScalarField& nMagSf = nMagSfPtr_();
 
@@ -192,6 +194,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
     const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
 
+    label totalFaceId = 0;
     forAll(intCoupledPatchIDs_, localPatchI)
     {
         label patchI = intCoupledPatchIDs_[localPatchI];
@@ -200,8 +203,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
         forAll(pp, patchFaceI)
         {
-            const labelList& faces = boundaryFaceFaces_[patchFaceI];
-            const labelList& cells = boundaryFaceCells_[patchFaceI];
+            const labelList& faces = boundaryFaceFaces_[totalFaceId];
+            const labelList& cells = boundaryFaceCells_[totalFaceId];
+
             const vector n = pp.faceNormals()[patchFaceI];
             const vector sf = pp.faceAreas()[patchFaceI];
 
@@ -231,7 +235,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
                     if
                     (
-                        ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
+                        mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
                       > minDelta
                     )
                     {
@@ -242,19 +246,17 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
                 }
                 nbrCf = oldCf[i + 1] + localDelta;
             }
-
             // Modify boundary
-            const label bFaceI = boundaryFaceOppositeFace_[patchFaceI];
+            const label bFaceI = boundaryFaceOppositeFace_[totalFaceId];
             const face f = regionMesh().faces()[bFaceI];
             const label cellI = cells[cells.size() - 1];
             newDelta += (deltaV[cellI]/mag(sf))*n;
             forAll(f, pti)
             {
                 const label pointI = f[pti];
-
                 if
                 (
-                    ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
+                    mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
                   > minDelta
                 )
                 {
@@ -262,9 +264,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
                     cellMoveMap[cellI] = 1;
                 }
             }
+            totalFaceId ++;
         }
     }
-
     // Move points
     regionMesh().movePoints(newPoints);
 
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
index 1f8eaf9e1b699e2e2a4b925169ec87f8578b95be..da2d90c9faebd7acee8f6d0f6e4a2c975c809706 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
@@ -28,6 +28,7 @@ License
 #include "polyMesh.H"
 #include "polyPatch.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 
 #include "addToRunTimeSelectionTable.H"
 
@@ -274,6 +275,49 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::sample
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::sampledPatch::sample
+(
+    const surfaceScalarField& sField
+) const
+{
+    return sampleField(sField);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledPatch::sample
+(
+    const surfaceVectorField& sField
+) const
+{
+    return sampleField(sField);
+}
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledPatch::sample
+(
+    const surfaceSphericalTensorField& sField
+) const
+{
+    return sampleField(sField);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledPatch::sample
+(
+    const surfaceSymmTensorField& sField
+) const
+{
+    return sampleField(sField);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledPatch::sample
+(
+    const surfaceTensorField& sField
+) const
+{
+    return sampleField(sField);
+}
+
 Foam::tmp<Foam::scalarField> Foam::sampledPatch::interpolate
 (
     const interpolation<scalar>& interpolator
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.H b/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
index 36590a747a74e0ec6aca912fcbed2414e6c5fef7..0ec8dbd89fbc77a206ff22c0833d01fc6f382f23 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.H
@@ -87,6 +87,13 @@ class sampledPatch
             const GeometricField<Type, fvPatchField, volMesh>& vField
         ) const;
 
+        //- sample surface field on faces
+        template <class Type>
+        tmp<Field<Type> > sampleField
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+        ) const;
+
         template <class Type>
         tmp<Field<Type> >
         interpolateField(const interpolation<Type>&) const;
@@ -203,6 +210,35 @@ public:
             const volTensorField&
         ) const;
 
+        //- Surface sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const surfaceScalarField&
+        ) const;
+
+        //- Surface Sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const surfaceVectorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const surfaceSphericalTensorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const surfaceSymmTensorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const surfaceTensorField&
+        ) const;
 
         //- interpolate field on surface
         virtual tmp<scalarField> interpolate
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
index 47b31e97c2012367f51d189cd8cd9848b9ada957..8725f412bea741fd2dbc4252736f50560b1acfff 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
@@ -48,6 +48,27 @@ Foam::sampledPatch::sampleField
 }
 
 
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledPatch::sampleField
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+) const
+{
+    // One value per face
+    tmp<Field<Type> > tvalues(new Field<Type>(patchFaceLabels_.size()));
+    Field<Type>& values = tvalues();
+
+    forAll(patchFaceLabels_, i)
+    {
+        label patchI = patchIDs_[patchIndex_[i]];
+        values[i] = sField.boundaryField()[patchI][patchFaceLabels_[i]];
+    }
+
+    return tvalues;
+}
+
+
 template <class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::sampledPatch::interpolateField
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
index 5c7c5abb90c50b43d7c66f60adff1d1f7acd323a..aaca4f087c8b4c2690b5d98af62aea43568a9703 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
@@ -240,6 +240,55 @@ Foam::scalar Foam::sampledSurface::area() const
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::sampledSurface::sample
+(
+    const surfaceScalarField& sField
+) const
+{
+    notImplemented("tmp<Foam::scalarField> sampledSurface::sample");
+    return tmp<scalarField>(NULL);
+}
+
+
+Foam::tmp<Foam::vectorField> Foam::sampledSurface::sample
+(
+    const surfaceVectorField& sField
+) const
+{
+    notImplemented("tmp<Foam::vectorField> sampledSurface::sample");
+    return tmp<vectorField>(NULL);
+}
+
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledSurface::sample
+(
+    const surfaceSphericalTensorField& sField
+) const
+{
+    notImplemented("tmp<Foam::sphericalTensorField> sampledSurface::sample");
+    return tmp<sphericalTensorField>(NULL);
+}
+
+
+Foam::tmp<Foam::symmTensorField> Foam::sampledSurface::sample
+(
+    const surfaceSymmTensorField& sField
+) const
+{
+    notImplemented("tmp<Foam::symmTensorField> sampledSurface::sample");
+    return tmp<symmTensorField>(NULL);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::sampledSurface::sample
+(
+    const surfaceTensorField& sField
+) const
+{
+    notImplemented("tmp<Foam::tensorField> sampledSurface::sample");
+    return tmp<tensorField>(NULL);
+}
+
+
 Foam::tmp<Foam::Field<Foam::scalar> >
 Foam::sampledSurface::project(const Field<scalar>& field) const
 {
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index 9843143107269573d564d0223d7af14e8873ea37..7ef9654ebdc651a88cbb2addd0fb7a03234b4570 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -57,6 +57,8 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "autoPtr.H"
 #include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "surfaceMesh.H"
 #include "polyMesh.H"
 #include "coordinateSystems.H"
 #include "interpolation.H"
@@ -352,6 +354,35 @@ public:
             const volTensorField&
         ) const = 0;
 
+        //- Surface sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const surfaceScalarField&
+        ) const;
+
+        //- Surface Sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const surfaceVectorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const surfaceSphericalTensorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const surfaceSymmTensorField&
+        ) const;
+
+        //- Surface sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const surfaceTensorField&
+        ) const;
 
         //- Interpolate field on surface
         virtual tmp<scalarField> interpolate
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index a8dfaebdfcf67f3fe59e09c44d98c9733b5ca657..778d33871ce6b38438605d4a7fb3aa6faff4ca7d 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -95,12 +95,7 @@ Foam::sampledSurfaces::sampledSurfaces
     fieldSelection_(),
     interpolationScheme_(word::null),
     mergeList_(),
-    formatter_(NULL),
-    scalarFields_(),
-    vectorFields_(),
-    sphericalTensorFields_(),
-    symmTensorFields_(),
-    tensorFields_()
+    formatter_(NULL)
 {
     if (Pstream::parRun())
     {
@@ -154,13 +149,6 @@ void Foam::sampledSurfaces::write()
         {
             if (debug)
             {
-                Pout<< "timeName = " << mesh_.time().timeName() << nl
-                    << "scalarFields    " << scalarFields_ << nl
-                    << "vectorFields    " << vectorFields_ << nl
-                    << "sphTensorFields " << sphericalTensorFields_ << nl
-                    << "symTensorFields " << symmTensorFields_ <<nl
-                    << "tensorFields    " << tensorFields_ <<nl;
-
                 Pout<< "Creating directory "
                     << outputPath_/mesh_.time().timeName() << nl << endl;
 
@@ -176,11 +164,19 @@ void Foam::sampledSurfaces::write()
             writeGeometry();
         }
 
-        sampleAndWrite(scalarFields_);
-        sampleAndWrite(vectorFields_);
-        sampleAndWrite(sphericalTensorFields_);
-        sampleAndWrite(symmTensorFields_);
-        sampleAndWrite(tensorFields_);
+        const IOobjectList objects(mesh_, mesh_.time().timeName());
+
+        sampleAndWrite<volScalarField>(objects);
+        sampleAndWrite<volVectorField>(objects);
+        sampleAndWrite<volSphericalTensorField>(objects);
+        sampleAndWrite<volSymmTensorField>(objects);
+        sampleAndWrite<volTensorField>(objects);
+
+        sampleAndWrite<surfaceScalarField>(objects);
+        sampleAndWrite<surfaceVectorField>(objects);
+        sampleAndWrite<surfaceSphericalTensorField>(objects);
+        sampleAndWrite<surfaceSymmTensorField>(objects);
+        sampleAndWrite<surfaceTensorField>(objects);
     }
 }
 
@@ -192,7 +188,6 @@ void Foam::sampledSurfaces::read(const dictionary& dict)
     if (surfacesFound)
     {
         dict.lookup("fields") >> fieldSelection_;
-        clearFieldGroups();
 
         dict.lookup("interpolationScheme") >> interpolationScheme_;
         const word writeType(dict.lookup("surfaceFormat"));
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index 0904c92ef0d835bce30ead29c2a9beeff08abd6d..ebff3bd94d8249db3fc4063dc4dc789953aea7e3 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -40,7 +40,9 @@ SourceFiles
 #include "sampledSurface.H"
 #include "surfaceWriter.H"
 #include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
 #include "wordReList.H"
+#include "IOobjectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,21 +62,6 @@ class sampledSurfaces
 {
     // Private classes
 
-        //- Class used for grouping field types
-        template<class Type>
-        class fieldGroup
-        :
-            public DynamicList<word>
-        {
-        public:
-
-            //- Construct null
-            fieldGroup()
-            :
-                DynamicList<word>(0)
-            {}
-        };
-
 
         //- Class used for surface merging information
         class mergeInfo
@@ -137,35 +124,43 @@ class sampledSurfaces
             //- Surface formatter
             autoPtr<surfaceWriter> formatter_;
 
-            //- Categorized scalar/vector/tensor fields
-            fieldGroup<scalar> scalarFields_;
-            fieldGroup<vector> vectorFields_;
-            fieldGroup<sphericalTensor> sphericalTensorFields_;
-            fieldGroup<symmTensor> symmTensorFields_;
-            fieldGroup<tensor> tensorFields_;
-
 
     // Private Member Functions
 
-        //- Clear old field groups
-        void clearFieldGroups();
 
-        //- Append fieldName to the appropriate group
-        label appendFieldGroup(const word& fieldName, const word& fieldType);
-
-        //- Classify field types, returns the number of fields
+        //- Return number of fields
         label classifyFields();
 
         //- Write geometry only
         void writeGeometry() const;
 
+        //- Write sampled fieldName on surface and on outputDir path
+        template<class Type>
+        void writeSurface
+        (
+            const Field<Type>& values,
+            const label surfI,
+            const word& fieldName,
+            const fileName& outputDir
+        );
+
         //- Sample and write a particular volume field
         template<class Type>
-        void sampleAndWrite(const GeometricField<Type, fvPatchField, volMesh>&);
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- Sample and write a particular surface field
+        template<class Type>
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
 
         //- Sample and write all the fields of the given type
         template<class Type>
-        void sampleAndWrite(fieldGroup<Type>&);
+        void sampleAndWrite(const IOobjectList& allObjects);
 
         //- Disallow default bitwise copy construct and assignment
         sampledSurfaces(const sampledSurfaces&);
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
index 0dda16b18011c93f400f8e7aabc2e3ee4113448e..6e67008de342e6cdd3cfaa3556e0b965e15095e5 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
@@ -30,94 +30,22 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sampledSurfaces::clearFieldGroups()
-{
-    scalarFields_.clear();
-    vectorFields_.clear();
-    sphericalTensorFields_.clear();
-    symmTensorFields_.clear();
-    tensorFields_.clear();
-}
-
-
-Foam::label Foam::sampledSurfaces::appendFieldGroup
-(
-    const word& fieldName,
-    const word& fieldType
-)
-{
-    if (fieldType == volScalarField::typeName)
-    {
-        scalarFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volVectorField::typeName)
-    {
-        vectorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volSphericalTensorField::typeName)
-    {
-        sphericalTensorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volSymmTensorField::typeName)
-    {
-        symmTensorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volTensorField::typeName)
-    {
-        tensorFields_.append(fieldName);
-        return 1;
-    }
-
-    return 0;
-}
-
-
 Foam::label Foam::sampledSurfaces::classifyFields()
 {
-    label nFields = 0;
-    clearFieldGroups();
-
-    // check files for a particular time
+     // check files for a particular time
     if (loadFromFiles_)
     {
         IOobjectList objects(mesh_, mesh_.time().timeName());
         wordList allFields = objects.sortedNames();
-
         labelList indices = findStrings(fieldSelection_, allFields);
-
-        forAll(indices, fieldI)
-        {
-            const word& fieldName = allFields[indices[fieldI]];
-
-            nFields += appendFieldGroup
-            (
-                fieldName,
-                objects.find(fieldName)()->headerClassName()
-            );
-        }
+        return indices.size();
     }
     else
     {
         wordList allFields = mesh_.sortedNames();
         labelList indices = findStrings(fieldSelection_, allFields);
-
-        forAll(indices, fieldI)
-        {
-            const word& fieldName = allFields[indices[fieldI]];
-
-            nFields += appendFieldGroup
-            (
-                fieldName,
-                mesh_.find(fieldName)()->type()
-            );
-        }
+        return indices.size();
     }
-
-    return nFields;
 }
 
 
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
index 78d7808d7af2cd09c5d80ddae203420eaa4a55e6..a56fcadf79df5a5065ef9484892f0e265ba18755 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
@@ -25,10 +25,86 @@ License
 
 #include "sampledSurfaces.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 #include "ListListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+void Foam::sampledSurfaces::writeSurface
+(
+    const Field<Type>& values,
+    const label surfI,
+    const word& fieldName,
+    const fileName& outputDir
+)
+{
+    const sampledSurface& s = operator[](surfI);
+
+    if (Pstream::parRun())
+    {
+        // Collect values from all processors
+        List<Field<Type> > gatheredValues(Pstream::nProcs());
+        gatheredValues[Pstream::myProcNo()] = values;
+        Pstream::gatherList(gatheredValues);
+
+        if (Pstream::master())
+        {
+            // Combine values into single field
+            Field<Type> allValues
+            (
+                ListListOps::combine<Field<Type> >
+                (
+                    gatheredValues,
+                    accessOp<Field<Type> >()
+                )
+            );
+
+            // Renumber (point data) to correspond to merged points
+            if (mergeList_[surfI].pointsMap.size() == allValues.size())
+            {
+                inplaceReorder(mergeList_[surfI].pointsMap, allValues);
+                allValues.setSize(mergeList_[surfI].points.size());
+            }
+
+            // Write to time directory under outputPath_
+            // skip surface without faces (eg, a failed cut-plane)
+            if (mergeList_[surfI].faces.size())
+            {
+                formatter_->write
+                (
+                    outputDir,
+                    s.name(),
+                    mergeList_[surfI].points,
+                    mergeList_[surfI].faces,
+                    fieldName,
+                    allValues,
+                    s.interpolate()
+                );
+            }
+        }
+    }
+    else
+    {
+        // Write to time directory under outputPath_
+        // skip surface without faces (eg, a failed cut-plane)
+        if (s.faces().size())
+        {
+            formatter_->write
+            (
+                outputDir,
+                s.name(),
+                s.points(),
+                s.faces(),
+                fieldName,
+                values,
+                s.interpolate()
+            );
+        }
+    }
+}
+
+
 template<class Type>
 void Foam::sampledSurfaces::sampleAndWrite
 (
@@ -65,125 +141,64 @@ void Foam::sampledSurfaces::sampleAndWrite
             values = s.sample(vField);
         }
 
-        if (Pstream::parRun())
-        {
-            // Collect values from all processors
-            List<Field<Type> > gatheredValues(Pstream::nProcs());
-            gatheredValues[Pstream::myProcNo()] = values;
-            Pstream::gatherList(gatheredValues);
-
-            if (Pstream::master())
-            {
-                // Combine values into single field
-                Field<Type> allValues
-                (
-                    ListListOps::combine<Field<Type> >
-                    (
-                        gatheredValues,
-                        accessOp<Field<Type> >()
-                    )
-                );
-
-                // Renumber (point data) to correspond to merged points
-                if (mergeList_[surfI].pointsMap.size() == allValues.size())
-                {
-                    inplaceReorder(mergeList_[surfI].pointsMap, allValues);
-                    allValues.setSize(mergeList_[surfI].points.size());
-                }
-
-                // Write to time directory under outputPath_
-                // skip surface without faces (eg, a failed cut-plane)
-                if (mergeList_[surfI].faces.size())
-                {
-                    formatter_->write
-                    (
-                        outputDir,
-                        s.name(),
-                        mergeList_[surfI].points,
-                        mergeList_[surfI].faces,
-                        fieldName,
-                        allValues,
-                        s.interpolate()
-                    );
-                }
-            }
-        }
-        else
-        {
-            // Write to time directory under outputPath_
-            // skip surface without faces (eg, a failed cut-plane)
-            if (s.faces().size())
-            {
-                formatter_->write
-                (
-                    outputDir,
-                    s.name(),
-                    s.points(),
-                    s.faces(),
-                    fieldName,
-                    values,
-                    s.interpolate()
-                );
-            }
-        }
+        writeSurface<Type>(values, surfI, fieldName, outputDir);
     }
 }
 
 
+
 template<class Type>
 void Foam::sampledSurfaces::sampleAndWrite
 (
-    fieldGroup<Type>& fields
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
 )
 {
-    if (fields.size())
+    const word& fieldName   = sField.name();
+    const fileName outputDir = outputPath_/sField.time().timeName();
+
+    forAll(*this, surfI)
     {
-        forAll(fields, fieldI)
+        const sampledSurface& s = operator[](surfI);
+        Field<Type> values = s.sample(sField);
+        writeSurface<Type>(values, surfI, fieldName, outputDir);
+    }
+}
+
+
+template<class GeoField>
+void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& allObjects)
+{
+    IOobjectList fields = allObjects.lookupClass(GeoField::typeName);
+    forAllConstIter(IOobjectList, fields, fieldIter)
+    {
+        forAll (fieldSelection_, fieldI)
         {
-            if (Pstream::master() && verbose_)
+            const wordRe field = fieldSelection_[fieldI];
+            if (field.match(fieldIter()->name()))
             {
-                Pout<< "sampleAndWrite: " << fields[fieldI] << endl;
-            }
+                if (Pstream::master() && verbose_)
+                {
+                    Pout<< "sampleAndWrite: " << field << endl;
+                }
 
-            if (loadFromFiles_)
-            {
-                sampleAndWrite
-                (
-                    GeometricField<Type, fvPatchField, volMesh>
+                if (loadFromFiles_)
+                {
+                    fieldIter()->readOpt() = IOobject::MUST_READ;
+                    sampleAndWrite
                     (
-                        IOobject
+                        GeoField
                         (
-                            fields[fieldI],
-                            mesh_.time().timeName(),
-                            mesh_,
-                            IOobject::MUST_READ,
-                            IOobject::NO_WRITE,
-                            false
-                        ),
-                        mesh_
-                    )
-                );
-            }
-            else
-            {
-                objectRegistry::const_iterator iter =
-                    mesh_.find(fields[fieldI]);
-
-                if
-                (
-                    iter != objectRegistry::end()
-                 && iter()->type()
-                 == GeometricField<Type, fvPatchField, volMesh>::typeName
-                )
+                            *fieldIter(),
+                            mesh_
+                        )
+                    );
+                }
+                else
                 {
-                   sampleAndWrite
-                   (
-                       mesh_.lookupObject
-                       <GeometricField<Type, fvPatchField, volMesh> >
-                       (
-                           fields[fieldI]
-                       )
-                   );
+                    sampleAndWrite
+                    (
+                        mesh_.lookupObject<GeoField>(fieldIter()->name())
+                    );
                 }
             }
         }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
index 355d1e1b36de89705c6e4bbbc83442721970e108..30c31c7471bffc10fc8a57c75bc336190538f89d 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
@@ -22,6 +22,7 @@ active true;
 infinitelyFastChemistryCoeffs
 {
     C       10;
+    semiImplicit    false;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones
index 6253de17410b3746318fb080b1d0c4421230796d..1895f851df5b662b8e39e0b663c4b2d2e358cba1 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pyrolysisZones
@@ -31,9 +31,9 @@ FoamFile
 
             radFluxName     Qr;
 
-            minimumDelta    1e-8;
+            minimumDelta    1e-12;
 
-            reactionDeltaMin 1e-8;
+            reactionDeltaMin 1e-12;
 
             moveMesh        false;
         }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
index ca98125036afe50f9646de87cbc321e80a75deea..df24e86f8dbf7d4d87a52a5ff8acc82de4f364c7 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
@@ -33,7 +33,7 @@ fvDOMCoeffs
     nPhi    3;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
     nTheta  6;          // polar angles in PI (from Z to X-Y plane)
     convergence 1e-4;   // convergence criteria for radiation iteration
-    maxIter 4;         // maximum number of iterations
+    maxIter 2;         // maximum number of iterations
 }
 
 // Number of flow iterations per radiation iteration
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
index 53c9826a9b9a59f5db2a3767088a8fb4ab919271..fb425995789bf14d15d5f89335a0ef39ec0064b9 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
@@ -16,25 +16,25 @@ FoamFile
 
 application     fireFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 
 stopAt          endTime;
 
-endTime         15.0;
+endTime         15;
 
 deltaT          0.03;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.5;
+writeInterval   1
 
 purgeWrite      0;
 
 writeFormat     ascii;
 
-writePrecision  6;
+writePrecision  12;
 
 writeCompression off;
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
index 9911bbedb397a9e40e00dbbe8fa96884a514623c..e683569e7a4ff14a8dce39837c09a946d0f1cecb 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
@@ -70,6 +70,7 @@ fluxRequired
 {
     default         no;
     p_rgh;
+    phiMesh;
 }
 
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
index e1d5b8db0b76425916fc4dad9c44280d248054df..f485444f2e70136dd27b16cfc4aac18533fedcea 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
@@ -32,7 +32,7 @@ solvers
         relTol          0;
     }
 
-    rhoThermo
+    rho
     {
         solver          PCG;
         preconditioner  DIC;
diff --git a/tutorials/compressible/rhoCentralFoam/forwardStep/0/Ma b/tutorials/compressible/rhoCentralFoam/forwardStep/0/Ma
deleted file mode 100644
index 061fe2b89964750c6e1964e4213be2c1cdf5825a..0000000000000000000000000000000000000000
--- a/tutorials/compressible/rhoCentralFoam/forwardStep/0/Ma
+++ /dev/null
@@ -1,264 +0,0 @@
-/*--------------------------------*- 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       volScalarField;
-    object      Ma;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 3;
-
-boundaryField
-{
-    inlet
-    {
-        type            calculated;
-        value           uniform 3;
-    }
-    outlet
-    {
-        type            calculated;
-        value           uniform 3;
-    }
-    bottom
-    {
-        type            symmetryPlane;
-    }
-    top
-    {
-        type            symmetryPlane;
-    }
-    obstacle
-    {
-        type            calculated;
-        value           nonuniform List<scalar>
-208
-(
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-3
-)
-;
-    }
-    defaultFaces
-    {
-        type            empty;
-    }
-}
-
-// ************************************************************************* //