diff --git a/src/OpenFOAM/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C
index e6f7c6c7db..57ebd84f76 100644
--- a/src/OpenFOAM/db/dictionary/dictionary.C
+++ b/src/OpenFOAM/db/dictionary/dictionary.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -222,6 +222,21 @@ Foam::dictionary::~dictionary()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::dictionary& Foam::dictionary::topDict() const
+{
+    const dictionary& p = parent();
+
+    if (&p != this && !p.name().empty())
+    {
+        return p.topDict();
+    }
+    else
+    {
+        return p;
+    }
+}
+
+
 Foam::label Foam::dictionary::startLineNumber() const
 {
     if (size())
diff --git a/src/OpenFOAM/db/dictionary/dictionary.H b/src/OpenFOAM/db/dictionary/dictionary.H
index 2eab9e2778..f7e7f8668c 100644
--- a/src/OpenFOAM/db/dictionary/dictionary.H
+++ b/src/OpenFOAM/db/dictionary/dictionary.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -250,6 +250,9 @@ public:
             return parent_;
         }
 
+        //- Return the top of the tree
+        const dictionary& topDict() const;
+
         //- Return line number of first token in dictionary
         label startLineNumber() const;
 
diff --git a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
index 1c774c352b..f855f79bc2 100644
--- a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
+++ b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
@@ -65,30 +65,12 @@ const Foam::word Foam::functionEntries::codeStream::codeTemplateC
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-const Foam::dictionary& Foam::functionEntries::codeStream::topDict
-(
-    const dictionary& dict
-)
-{
-    const dictionary& p = dict.parent();
-
-    if (&p != &dict && !p.name().empty())
-    {
-        return topDict(p);
-    }
-    else
-    {
-        return dict;
-    }
-}
-
-
 Foam::dlLibraryTable& Foam::functionEntries::codeStream::libs
 (
     const dictionary& dict
 )
 {
-    const IOdictionary& d = static_cast<const IOdictionary&>(topDict(dict));
+    const IOdictionary& d = static_cast<const IOdictionary&>(dict.topDict());
     return const_cast<Time&>(d.time()).libs();
 }
 
@@ -114,7 +96,7 @@ Foam::functionEntries::codeStream::getFunction
 
     // see if library is loaded
     void* lib = NULL;
-    if (isA<IOdictionary>(topDict(parentDict)))
+    if (isA<IOdictionary>(parentDict.topDict()))
     {
         lib = libs(parentDict).findLibrary(libPath);
     }
@@ -129,7 +111,7 @@ Foam::functionEntries::codeStream::getFunction
     // avoid compilation if possible by loading an existing library
     if (!lib)
     {
-        if (isA<IOdictionary>(topDict(parentDict)))
+        if (isA<IOdictionary>(parentDict.topDict()))
         {
             // Cached access to dl libs. Guarantees clean up upon destruction
             // of Time.
@@ -267,7 +249,7 @@ Foam::functionEntries::codeStream::getFunction
             }
         }
 
-        if (isA<IOdictionary>(topDict(parentDict)))
+        if (isA<IOdictionary>(parentDict.topDict()))
         {
             // Cached access to dl libs. Guarantees clean up upon destruction
             // of Time.
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C
index f892a0d7c8..faf2922955 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.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
@@ -26,11 +26,8 @@ License
 #include "mapPolyMesh.H"
 #include "polyMesh.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::mapPolyMesh::mapPolyMesh
 (
     const polyMesh& mesh,
@@ -115,7 +112,6 @@ Foam::mapPolyMesh::mapPolyMesh
 }
 
 
-// Construct from components and optionally reuse storage
 Foam::mapPolyMesh::mapPolyMesh
 (
     const polyMesh& mesh,
@@ -177,31 +173,32 @@ Foam::mapPolyMesh::mapPolyMesh
     oldPatchStarts_(oldPatchStarts, reUse),
     oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse)
 {
-    // Calculate old patch sizes
-    for (label patchI = 0; patchI < oldPatchStarts_.size() - 1; patchI++)
+    if (oldPatchStarts_.size() > 0)
     {
-        oldPatchSizes_[patchI] =
-            oldPatchStarts_[patchI + 1] - oldPatchStarts_[patchI];
-    }
-
-    // Set the last one by hand
-    const label lastPatchID = oldPatchStarts_.size() - 1;
-
-    oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID];
-
-    if (polyMesh::debug)
-    {
-        if (min(oldPatchSizes_) < 0)
+        // Calculate old patch sizes
+        for (label patchI = 0; patchI < oldPatchStarts_.size() - 1; patchI++)
         {
-            FatalErrorIn("mapPolyMesh::mapPolyMesh(...)")
-                << "Calculated negative old patch size.  Error in mapping data"
-                << abort(FatalError);
+            oldPatchSizes_[patchI] =
+                oldPatchStarts_[patchI + 1] - oldPatchStarts_[patchI];
+        }
+
+        // Set the last one by hand
+        const label lastPatchID = oldPatchStarts_.size() - 1;
+
+        oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID];
+
+        if (polyMesh::debug)
+        {
+            if (min(oldPatchSizes_) < 0)
+            {
+                FatalErrorIn("mapPolyMesh::mapPolyMesh(...)")
+                    << "Calculated negative old patch size."
+                    << "  Error in mapping data"
+                    << abort(FatalError);
+            }
         }
     }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
deleted file mode 100644
index e8d96e9b63..0000000000
--- a/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
+++ /dev/null
@@ -1,44 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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 "makeCombustionTypes.H"
-#include "psiCombustionModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace combustionModels
-{
-    makeCombustionTypes
-    (
-        infinitelyFastChemistry,
-        psiCombustionModel
-    );
-}
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
index b052b4de86..f99612eed6 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
@@ -58,10 +58,8 @@ void Foam::layerAdditionRemoval::checkDefinition()
 {
     if (!faceZoneID_.active())
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Master face zone named " << faceZoneID_.name()
+        FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
+            << "Master face zone named " << faceZoneID_.name()
             << " cannot be found."
             << abort(FatalError);
     }
@@ -79,13 +77,15 @@ void Foam::layerAdditionRemoval::checkDefinition()
             << abort(FatalError);
     }
 
-    if (topoChanger().mesh().faceZones()[faceZoneID_.index()].empty())
+    label nFaces = topoChanger().mesh().faceZones()[faceZoneID_.index()].size();
+
+    reduce(nFaces, sumOp<label>());
+
+    if (nFaces == 0)
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Face extrusion zone contains no faces. "
-            << " Please check your mesh definition."
+        FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
+            << "Face extrusion zone contains no faces. "
+            << "Please check your mesh definition."
             << abort(FatalError);
     }
 
@@ -120,16 +120,18 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
     const label index,
-    const polyTopoChanger& mme,
+    const polyTopoChanger& ptc,
     const word& zoneName,
     const scalar minThickness,
-    const scalar maxThickness
+    const scalar maxThickness,
+    const Switch thicknessFromVolume
 )
 :
-    polyMeshModifier(name, index, mme, true),
-    faceZoneID_(zoneName, mme.mesh().faceZones()),
+    polyMeshModifier(name, index, ptc, true),
+    faceZoneID_(zoneName, ptc.mesh().faceZones()),
     minLayerThickness_(minThickness),
     maxLayerThickness_(maxThickness),
+    thicknessFromVolume_(thicknessFromVolume),
     oldLayerThickness_(-1.0),
     pointsPairingPtr_(NULL),
     facesPairingPtr_(NULL),
@@ -145,13 +147,17 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
     const word& name,
     const dictionary& dict,
     const label index,
-    const polyTopoChanger& mme
+    const polyTopoChanger& ptc
 )
 :
-    polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
-    faceZoneID_(dict.lookup("faceZoneName"), mme.mesh().faceZones()),
+    polyMeshModifier(name, index, ptc, Switch(dict.lookup("active"))),
+    faceZoneID_(dict.lookup("faceZoneName"), ptc.mesh().faceZones()),
     minLayerThickness_(readScalar(dict.lookup("minLayerThickness"))),
     maxLayerThickness_(readScalar(dict.lookup("maxLayerThickness"))),
+    thicknessFromVolume_
+    (
+        dict.lookupOrDefault<Switch>("thicknessFromVolume", true)
+    ),
     oldLayerThickness_(readOldThickness(dict)),
     pointsPairingPtr_(NULL),
     facesPairingPtr_(NULL),
@@ -188,11 +194,13 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     // Layer removal:
     //     When the min thickness falls below the threshold, trigger removal.
 
-    const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
+    const polyMesh& mesh = topoChanger().mesh();
+
+    const faceZone& fz = mesh.faceZones()[faceZoneID_.index()];
     const labelList& mc = fz.masterCells();
 
-    const scalarField& V = topoChanger().mesh().cellVolumes();
-    const vectorField& S = topoChanger().mesh().faceAreas();
+    const scalarField& V = mesh.cellVolumes();
+    const vectorField& S = mesh.faceAreas();
 
     if (min(V) < -VSMALL)
     {
@@ -205,63 +213,68 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     scalar avgDelta = 0;
     scalar minDelta = GREAT;
     scalar maxDelta = 0;
+    label nDelta = 0;
 
-    forAll(fz, faceI)
+    if (thicknessFromVolume_)
     {
-        scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
-        avgDelta += curDelta;
-        minDelta = min(minDelta, curDelta);
-        maxDelta = max(maxDelta, curDelta);
+        // Thickness calculated from cell volume/face area
+        forAll(fz, faceI)
+        {
+            scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
+            avgDelta += curDelta;
+            minDelta = min(minDelta, curDelta);
+            maxDelta = max(maxDelta, curDelta);
+        }
+
+        nDelta = fz.size();
+    }
+    else
+    {
+        // Thickness calculated from edges on layer
+        const Map<label>& zoneMeshPointMap = fz().meshPointMap();
+
+        // Edges with only one point on zone
+        forAll(mc, faceI)
+        {
+            const cell& cFaces = mesh.cells()[mc[faceI]];
+            const edgeList cellEdges(cFaces.edges(mesh.faces()));
+
+            forAll(cellEdges, i)
+            {
+                const edge& e = cellEdges[i];
+
+                if (zoneMeshPointMap.found(e[0]))
+                {
+                    if (!zoneMeshPointMap.found(e[1]))
+                    {
+                        scalar curDelta = e.mag(mesh.points());
+                        avgDelta += curDelta;
+                        nDelta++;
+                        minDelta = min(minDelta, curDelta);
+                        maxDelta = max(maxDelta, curDelta);
+                    }
+                }
+                else
+                {
+                    if (zoneMeshPointMap.found(e[1]))
+                    {
+                        scalar curDelta = e.mag(mesh.points());
+                        avgDelta += curDelta;
+                        nDelta++;
+                        minDelta = min(minDelta, curDelta);
+                        maxDelta = max(maxDelta, curDelta);
+                    }
+                }
+            }
+        }
     }
 
-    avgDelta /= fz.size();
+    reduce(minDelta, minOp<scalar>());
+    reduce(maxDelta, maxOp<scalar>());
+    reduce(avgDelta, sumOp<scalar>());
+    reduce(nDelta, sumOp<label>());
 
-    ////MJ Alternative thickness determination
-    //{
-    //    // Edges on layer.
-    //    const Map<label>& zoneMeshPointMap = fz().meshPointMap();
-    //
-    //    label nDelta = 0;
-    //
-    //    // Edges with only one point on zone
-    //    const polyMesh& mesh = topoChanger().mesh();
-    //
-    //    forAll(mc, faceI)
-    //    {
-    //        const cell& cFaces = mesh.cells()[mc[faceI]];
-    //        const edgeList cellEdges(cFaces.edges(mesh.faces()));
-    //
-    //        forAll(cellEdges, i)
-    //        {
-    //            const edge& e = cellEdges[i];
-    //
-    //            if (zoneMeshPointMap.found(e[0]))
-    //            {
-    //                if (!zoneMeshPointMap.found(e[1]))
-    //                {
-    //                    scalar curDelta = e.mag(mesh.points());
-    //                    avgDelta += curDelta;
-    //                    nDelta++;
-    //                    minDelta = min(minDelta, curDelta);
-    //                    maxDelta = max(maxDelta, curDelta);
-    //                }
-    //            }
-    //            else
-    //            {
-    //                if (zoneMeshPointMap.found(e[1]))
-    //                {
-    //                    scalar curDelta = e.mag(mesh.points());
-    //                    avgDelta += curDelta;
-    //                    nDelta++;
-    //                    minDelta = min(minDelta, curDelta);
-    //                    maxDelta = max(maxDelta, curDelta);
-    //                }
-    //            }
-    //        }
-    //    }
-    //
-    //    avgDelta /= nDelta;
-    //}
+    avgDelta /= nDelta;
 
     if (debug)
     {
@@ -286,7 +299,6 @@ bool Foam::layerAdditionRemoval::changeTopology() const
         }
 
         // No topological changes allowed before first mesh motion
-        //
         oldLayerThickness_ = avgDelta;
 
         topologicalChange = false;
@@ -314,7 +326,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
                             << "Triggering layer removal" << endl;
                     }
 
-                    triggerRemoval_ = topoChanger().mesh().time().timeIndex();
+                    triggerRemoval_ = mesh.time().timeIndex();
 
                     // Old thickness looses meaning.
                     // Set it up to indicate layer removal
@@ -346,7 +358,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
                     << "Triggering layer addition" << endl;
             }
 
-            triggerAddition_ = topoChanger().mesh().time().timeIndex();
+            triggerAddition_ = mesh.time().timeIndex();
 
             // Old thickness looses meaning.
             // Set it up to indicate layer removal
@@ -377,8 +389,8 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         if (debug)
         {
             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer removal. " << endl;
+                << "for object " << name() << " : "
+                << "Clearing addressing after layer removal" << endl;
         }
 
         triggerRemoval_ = -1;
@@ -393,8 +405,8 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         if (debug)
         {
             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer addition. " << endl;
+                << "for object " << name() << " : "
+                << "Clearing addressing after layer addition" << endl;
         }
 
         triggerAddition_ = -1;
@@ -408,8 +420,8 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
     if (debug)
     {
         Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
-            << " for object " << name() << " : "
-            << "Clearing addressing on external request. ";
+            << "for object " << name() << " : "
+            << "Clearing addressing on external request";
 
         if (pointsPairingPtr_ || facesPairingPtr_)
         {
@@ -471,7 +483,8 @@ void Foam::layerAdditionRemoval::write(Ostream& os) const
         << faceZoneID_ << nl
         << minLayerThickness_ << nl
         << oldLayerThickness_ << nl
-        << maxLayerThickness_ << endl;
+        << maxLayerThickness_ << nl
+        << thicknessFromVolume_ << endl;
 }
 
 
@@ -486,6 +499,8 @@ void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
         << token::END_STATEMENT << nl
         << "    maxLayerThickness " << maxLayerThickness_
         << token::END_STATEMENT << nl
+        << "    thicknessFromVolume " << thicknessFromVolume_
+        << token::END_STATEMENT << nl
         << "    oldLayerThickness " << oldLayerThickness_
         << token::END_STATEMENT << nl
         << "    active " << active()
diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H
index daef765097..edb3218eae 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.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
@@ -65,6 +65,10 @@ class layerAdditionRemoval
         //- Max thickness of extrusion layer.  Triggers layer addition
         mutable scalar maxLayerThickness_;
 
+        //- Switch to calculate thickness as volume/area
+        //  If false, thickness calculated from edges
+        const bool thicknessFromVolume_;
+
         //- Layer thickness from previous step
         //  Used to decide whether to add or remove layers
         mutable scalar oldLayerThickness_;
@@ -78,7 +82,7 @@ class layerAdditionRemoval
         //- Layer removal trigger time index
         mutable label triggerRemoval_;
 
-        //- Layer addition trigger  time index
+        //- Layer addition trigger time index
         mutable label triggerAddition_;
 
 
@@ -120,6 +124,7 @@ class layerAdditionRemoval
             //- Clear addressing
             void clearAddressing() const;
 
+
         // Helpers
 
             //- Optionally read old thickness
@@ -149,10 +154,11 @@ public:
         (
             const word& name,
             const label index,
-            const polyTopoChanger& mme,
+            const polyTopoChanger& ptc,
             const word& zoneName,
             const scalar minThickness,
-            const scalar maxThickness
+            const scalar maxThickness,
+            const Switch thicknessFromVolume = true
         );
 
         //- Construct from dictionary
@@ -161,7 +167,7 @@ public:
             const word& name,
             const dictionary& dict,
             const label index,
-            const polyTopoChanger& mme
+            const polyTopoChanger& ptc
         );
 
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index 49ee1a445e..23e39f203d 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -818,22 +818,25 @@ void Foam::polyTopoChange::getFaceOrder
     patchSizes.setSize(nPatches_);
     patchSizes = 0;
 
-    patchStarts[0] = newFaceI;
-
-    for (label faceI = 0; faceI < nActiveFaces; faceI++)
+    if (nPatches_ > 0)
     {
-        if (region_[faceI] >= 0)
+        patchStarts[0] = newFaceI;
+
+        for (label faceI = 0; faceI < nActiveFaces; faceI++)
         {
-            patchSizes[region_[faceI]]++;
+            if (region_[faceI] >= 0)
+            {
+                patchSizes[region_[faceI]]++;
+            }
         }
-    }
 
-    label faceI = patchStarts[0];
+        label faceI = patchStarts[0];
 
-    forAll(patchStarts, patchI)
-    {
-        patchStarts[patchI] = faceI;
-        faceI += patchSizes[patchI];
+        forAll(patchStarts, patchI)
+        {
+            patchStarts[patchI] = faceI;
+            faceI += patchSizes[patchI];
+        }
     }
 
     //if (debug)
diff --git a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
index d2ca617bbd..e6863ad0db 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
@@ -36,7 +36,7 @@ Foam::PilchErdman<CloudType>::PilchErdman
 :
     BreakupModel<CloudType>(dict, owner, typeName),
     B1_(0.375),
-    B2_(0.236)
+    B2_(0.2274)
 {
     if (!this->defaultCoeffs(true))
     {
@@ -90,58 +90,79 @@ bool Foam::PilchErdman<CloudType>::update
     scalar& massChild
 )
 {
-    scalar semiMass = nParticle*pow3(d);
-    scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
+    // Weber number - eq (1)
+    scalar We = rhoc*sqr(Urmag)*d/sigma;
+
+    // Ohnesorge number - eq (2)
     scalar Oh = mu/sqrt(rho*d*sigma);
 
-    scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
+    // Critical Weber number - eq (5)
+    scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
 
     if (We > Wec)
     {
-        // We > 1335, wave crest stripping
+        // We > 2670, wave crest stripping - eq (12)
         scalar taubBar = 5.5;
 
-        if (We < 1335)
+        if (We < 2670)
         {
-            if (We > 175.0)
+            if (We > 351)
             {
-                // sheet stripping
-                taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
+                // sheet stripping - eq (11)
+                taubBar = 0.766*pow(We - 12.0, 0.25);
             }
-            else if (We > 22.0)
+            else if (We > 45)
             {
-                // Bag-and-stamen breakup
-                taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
+                // bag-and-stamen breakup  - eq (10)
+                taubBar = 14.1*pow(We - 12.0, 0.25);
             }
-            else if (We > 9.0)
+            else if (We > 18)
             {
-                // Bag breakup
-                taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
+                // bag breakup - eq (9)
+                taubBar = 2.45*pow(We - 12.0, 0.25);
             }
-            else if (We > 6.0)
+            else if (We > 12)
             {
-                // Vibrational breakup
-                taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+                // vibrational breakup - eq (8)
+                taubBar = 6.0*pow(We - 12.0, -0.25);
+            }
+            else
+            {
+                // no break-up
+                taubBar = GREAT;
             }
         }
 
         scalar rho12 = sqrt(rhoc/rho);
 
-        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
+        // velocity of fragmenting drop - eq (20)
+        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
+
+        // maximum stable diameter - eq (33)
         scalar Vd1 = sqr(1.0 - Vd/Urmag);
         Vd1 = max(Vd1, SMALL);
-        scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
-        scalar A = Urmag*rho12/d;
+        scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
 
-        scalar taub = taubBar/A;
+        if (d < dStable)
+        {
+            // droplet diameter already stable = no break-up
+            // - do not update d and nParticle
+            return false;
+        }
+        else
+        {
+            scalar semiMass = nParticle*pow3(d);
 
-        scalar frac = dt/taub;
+            // invert eq (3) to create a dimensional break-up time
+            scalar taub = taubBar*d/(Urmag*rho12);
 
-        // update the droplet diameter according to the rate eq. (implicitly)
-        d = (d + frac*Ds)/(1.0 + frac);
+            // update droplet diameter according to the rate eq (implicitly)
+            scalar frac = dt/taub;
+            d = (d + frac*dStable)/(1.0 + frac);
 
-        // correct the number of particles to conserve mass
-        nParticle = semiMass/pow3(d);
+            // correct the number of particles to conserve mass
+            nParticle = semiMass/pow3(d);
+        }
     }
 
     return false;
diff --git a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H
index 830c1feeb5..09d1f3fe03 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.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
@@ -25,7 +25,7 @@ Class
     Foam::PilchErdman
 
 Description
-    secondary breakup model
+    Particle secondary breakup model, based on the reference:
 
     @verbatim
     Pilch, M. and Erdman, C.A.
@@ -35,6 +35,24 @@ Description
     Int. J. Multiphase Flows 13 (1987), 741-757
     @endverbatim
 
+    The droplet fragment velocity is described by the equation:
+
+    \f[
+        V_d = V sqrt(epsilon)(B1 T + B2 T^2)
+    \f]
+
+    Where:
+        V_d     : fragment velocity
+        V       : magnitude of the relative velocity
+        epsilon : density ratio (rho_carrier/rho_droplet)
+        T       : characteristic break-up time
+        B1, B2  : model input coefficients
+
+    The authors suggest that:
+        compressible flow   : B1 = 0.75*1.0; B2 = 3*0.116
+        incompressible flow : B1 = 0.75*0.5; B2 = 3*0.0758
+
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef PilchErdman_H
diff --git a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
index 35374a7d94..e80383a303 100644
--- a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
+++ b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.H
@@ -27,6 +27,19 @@ Class
 Description
     A surface writer for the Nastran file format - both surface mesh and fields
 
+    formatOptions
+    {
+        nastran
+        {
+            // From OpenFOAM field name to Nastran field name
+            fields ((pMean PLOAD2));
+            // Optional scale
+            scale 2.0;
+            // Optional format
+            format free;    //short, long, free
+        }
+    };
+
 SourceFiles
     nastranSurfaceWriter.C
     nastranSurfaceWriterTemplates.C