diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/alphaEqnComp.H b/applications/solvers/multiphase/multiphaseEulerFoam/alphaEqnComp.H
deleted file mode 100644
index f0ebde6332ffd8a5bc71b59d1bde0c98593ac65a..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/multiphaseEulerFoam/alphaEqnComp.H
+++ /dev/null
@@ -1,105 +0,0 @@
-surfaceScalarField alphaPhi1("alphaPhi1", phi1);
-surfaceScalarField alphaPhi2("alphaPhi2", phi2);
-
-{
-    word scheme("div(phi,alpha1)");
-    word schemer("div(phir,alpha1)");
-
-    surfaceScalarField phic("phic", phi);
-    surfaceScalarField phir("phir", phi1 - phi2);
-
-    if (g0.value() > 0.0)
-    {
-        surfaceScalarField alpha1f = fvc::interpolate(alpha1);
-        surfaceScalarField phipp = ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-        phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
-    }
-
-    for (int acorr=0; acorr<nAlphaCorr; acorr++)
-    {
-        volScalarField::DimensionedInternalField Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
-
-        volScalarField::DimensionedInternalField Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            fvc::div(phi)*min(alpha1, scalar(1))
-        );
-
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
-
-
-        fvScalarMatrix alpha1Eqn
-        (
-             fvm::ddt(alpha1)
-           + fvm::div(phic, alpha1, scheme)
-           + fvm::div(-fvc::flux(-phir, alpha2, schemer), alpha1, schemer)
-          ==
-             fvm::Sp(Sp, alpha1) + Su
-        );
-
-        if (g0.value() > 0.0)
-        {
-            ppMagf = rU1Af*fvc::interpolate
-            (
-                (1.0/(rho1*(alpha1 + scalar(0.0001))))
-               *g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
-            );
-
-            alpha1Eqn -= fvm::laplacian
-            (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
-                alpha1,
-                "laplacian(alphaPpMag,alpha1)"
-            );
-        }
-
-        alpha1Eqn.relax();
-        alpha1Eqn.solve();
-
-        //***HGW temporary boundedness-fix pending the introduction of MULES
-        alpha1 = max(min(alpha1, 1.0), 0.0);
-
-        #include "packingLimiter.H"
-
-        alphaPhi1 = alpha1Eqn.flux();
-        alphaPhi2 = phi - alphaPhi1;
-        alpha2 = scalar(1) - alpha1;
-
-        Info<< "Dispersed phase volume fraction = "
-            << alpha1.weightedAverage(mesh.V()).value()
-            << "  Min(alpha1) = " << min(alpha1).value()
-            << "  Max(alpha1) = " << max(alpha1).value()
-            << endl;
-    }
-}
-
-rho = alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
index f3ec0396843c2a4038eb16064e280d2398cb28f1..acab521e2009511a9098a879f060eb9765d96aa5 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C
@@ -46,11 +46,17 @@ Foam::dragModel::dragModel
     interfaceDict_(interfaceDict),
     phase1_(phase1),
     phase2_(phase2),
-    residualDrag_
+    residualPhaseFraction_
     (
-        "residualDrag",
-        dimensionSet(1, -3, -1, 0, 0),
-        interfaceDict.lookup("residualDrag")
+        "residualPhaseFraction",
+        dimless,
+        interfaceDict.lookup("residualPhaseFraction")
+    ),
+    residualSlip_
+    (
+        "residualSlip",
+        dimVelocity,
+        interfaceDict.lookup("residualSlip")
     )
 {}
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
index 37413d83e4267de6bcde2f4a40acc1adff5e72b2..a775d91e4566a168dc0a22dddcd3e54496b73e67 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H
@@ -57,7 +57,8 @@ protected:
         const dictionary& interfaceDict_;
         const phaseModel& phase1_;
         const phaseModel& phase2_;
-        dimensionedScalar residualDrag_;
+        dimensionedScalar residualPhaseFraction_;
+        dimensionedScalar residualSlip_;
 
 public:
 
@@ -117,9 +118,14 @@ public:
             return phase2_;
         }
 
-        const dimensionedScalar& residualDrag() const
+        const dimensionedScalar& residualPhaseFraction() const
         {
-            return residualDrag_;
+            return residualPhaseFraction_;
+        }
+
+        const dimensionedScalar& residualSlip() const
+        {
+            return residualSlip_;
         }
 
         //- the dragfunction K used in the momentum eq.
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 2403597d9cf5837bc5d33f3111f429673fc33b0c..6ae60e42dc8ca539fa7ff278a62011885ba326bf 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -30,6 +30,7 @@ License
 #include "MULES.H"
 #include "fvcSnGrad.H"
 #include "fvcFlux.H"
+#include "fvcAverage.H"
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
@@ -626,9 +627,20 @@ Foam::multiphaseSystem::dragCoeffs() const
         (
             iter.key(),
             (
-                dm.phase1()*dm.phase2()
-               *dm.K(mag(dm.phase1().U() - dm.phase2().U()))
-              + dm.residualDrag()
+                max
+                (
+                    fvc::average(dm.phase1())*fvc::average(dm.phase2()),
+                    //dm.phase1()*dm.phase2(),
+                    dm.residualPhaseFraction()
+                )
+               *dm.K
+                (
+                    max
+                    (
+                        mag(dm.phase1().U() - dm.phase2().U()),
+                        dm.residualSlip()
+                    )
+                )
             ).ptr()
         );
     }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index ad420d096f5587fd5e67c6b1251e008c3c3373bd..edafdb50f912e02a89acc6133fc0d7cf7491cd5c 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -146,6 +146,9 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
 
     decompositionMethod& decomposer = decomposerPtr_();
 
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
+
     // For each cell in the mesh has it been determined if it is fully
     // inside, outside, or overlaps the surface
     labelList volumeStatus
@@ -214,10 +217,10 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                         volumeStatus[cellI] = searchableSurface::UNKNOWN;
                     }
 
-                    cellWeights.internalField()[cellI] = max
+                    icellWeights[cellI] = max
                     (
                         1.0,
-                        cellWeights.internalField()[cellI]/8.0
+                        icellWeights[cellI]/8.0
                     );
                 }
 
@@ -405,7 +408,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
             (
                 mesh_,
                 mesh_.cellCentres(),
-                cellWeights
+                icellWeights
             );
 
             fvMeshDistribute distributor(mesh_, mergeDist_);
@@ -629,6 +632,8 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
     volScalarField& cellWeights
 ) const
 {
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
     labelHashSet cellsToRefine;
 
     // Determine/update the status of each cell
@@ -650,7 +655,7 @@ Foam::labelList Foam::backgroundMeshDecomposition::selectRefinementCells
                 (
                     cellI,
                     volumeStatus[cellI],
-                    cellWeights.internalField()[cellI]
+                    icellWeights[cellI]
                 )
             )
             {
@@ -871,15 +876,17 @@ Foam::backgroundMeshDecomposition::distribute
         mesh_.write();
     }
 
+    volScalarField::InternalField& icellWeights = cellWeights.internalField();
+
     while (true)
     {
         // Refine large cells if necessary
 
         label nOccupiedCells = 0;
 
-        forAll(cellWeights.internalField(), cI)
+        forAll(icellWeights, cI)
         {
-            if (cellWeights.internalField()[cI] > 1 - SMALL)
+            if (icellWeights[cI] > 1 - SMALL)
             {
                 nOccupiedCells++;
             }
@@ -910,9 +917,9 @@ Foam::backgroundMeshDecomposition::distribute
 
         labelHashSet cellsToRefine;
 
-        forAll(cellWeights, cWI)
+        forAll(icellWeights, cWI)
         {
-            if (cellWeights.internalField()[cWI] > cellWeightLimit)
+            if (icellWeights[cWI] > cellWeightLimit)
             {
                 cellsToRefine.insert(cWI);
             }
@@ -943,7 +950,7 @@ Foam::backgroundMeshDecomposition::distribute
         {
             label cellI = newCellsToRefine[nCTRI];
 
-            cellWeights.internalField()[cellI] /= 8.0;
+            icellWeights[cellI] /= 8.0;
         }
 
         // Mesh changing engine.
@@ -1075,9 +1082,9 @@ Foam::backgroundMeshDecomposition::distribute
         printMeshData(mesh_);
 
         Pout<< "    Pre distribute sum(cellWeights) "
-            << sum(cellWeights.internalField())
+            << sum(icellWeights)
             << " max(cellWeights) "
-            << max(cellWeights.internalField())
+            << max(icellWeights)
             << endl;
     }
 
@@ -1085,7 +1092,7 @@ Foam::backgroundMeshDecomposition::distribute
     (
         mesh_,
         mesh_.cellCentres(),
-        cellWeights
+        icellWeights
     );
 
     Info<< "    Redistributing background mesh cells" << endl;
@@ -1101,9 +1108,9 @@ Foam::backgroundMeshDecomposition::distribute
         printMeshData(mesh_);
 
         Pout<< "    Post distribute sum(cellWeights) "
-            << sum(cellWeights.internalField())
+            << sum(icellWeights)
             << " max(cellWeights) "
-            << max(cellWeights.internalField())
+            << max(icellWeights)
             << endl;
 
         // const_cast<Time&>(mesh_.time())++;
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 11ecb12081431a7dd61d5543644b22bb2e4956d5..5a2e599a6490a00a1891029dfe4f4ede0d57239b 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -87,8 +87,6 @@ void Foam::vtkPV3Foam::reduceMemory()
 }
 
 
-
-
 int Foam::vtkPV3Foam::setTime(int nRequest, const double requestTimes[])
 {
     Time& runTime = dbPtr_();
@@ -214,6 +212,7 @@ void Foam::vtkPV3Foam::updateMeshPartsStatus()
     }
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::vtkPV3Foam::vtkPV3Foam
@@ -378,7 +377,7 @@ void Foam::vtkPV3Foam::updateInfo()
 
     // Update mesh parts list - add Lagrangian at the bottom
     updateInfoInternalMesh(partSelection);
-    updateInfoPatches(partSelection);
+    updateInfoPatches(partSelection, enabledEntries);
     updateInfoSets(partSelection);
     updateInfoZones(partSelection);
     updateInfoLagrangian(partSelection);
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
index 7d9a481b49a232fb4d1600492d7612a62a71a5e8..b4d031b7ecc182607e2e94690420ec239110bb7f 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
@@ -347,7 +347,7 @@ class vtkPV3Foam
             void updateInfoLagrangian(vtkDataArraySelection*);
 
             //- Patch info
-            void updateInfoPatches(vtkDataArraySelection*);
+            void updateInfoPatches(vtkDataArraySelection*, stringList&);
 
             //- Set info
             void updateInfoSets(vtkDataArraySelection*);
@@ -554,18 +554,6 @@ class vtkPV3Foam
                 const labelList& faceLabels
             );
 
-            //- face set/zone field
-            template<class Type>
-            void convertFaceField
-            (
-                const GeometricField<Type, fvPatchField, volMesh>&,
-                vtkMultiBlockDataSet* output,
-                const arrayRange&,
-                const label datasetNo,
-                const fvMesh&,
-                const faceSet&
-            );
-
             //- Lagrangian fields - all types
             template<class Type>
             void convertLagrangianFields
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H
index 030c9aca05a5407dfede8c02354fc7f5530ed384..b37255a8f100f9a112865c0fc9345d998e746265 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H
@@ -110,79 +110,6 @@ void Foam::vtkPV3Foam::convertFaceField
 }
 
 
-template<class Type>
-void Foam::vtkPV3Foam::convertFaceField
-(
-    const GeometricField<Type, fvPatchField, volMesh>& tf,
-    vtkMultiBlockDataSet* output,
-    const arrayRange& range,
-    const label datasetNo,
-    const fvMesh& mesh,
-    const faceSet& fSet
-)
-{
-    const label nComp = pTraits<Type>::nComponents;
-    const label nInternalFaces = mesh.nInternalFaces();
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeigh = mesh.faceNeighbour();
-
-    vtkFloatArray* cellData = vtkFloatArray::New();
-    cellData->SetNumberOfTuples(fSet.size());
-    cellData->SetNumberOfComponents(nComp);
-    cellData->Allocate(nComp*fSet.size());
-    cellData->SetName(tf.name().c_str());
-
-    if (debug)
-    {
-        Info<< "convert convertFaceField: "
-            << tf.name()
-            << " size = " << tf.size()
-            << " nComp=" << nComp
-            << " nTuples = " << fSet.size() <<  endl;
-    }
-
-    float vec[nComp];
-
-    // for interior faces: average owner/neighbour
-    // for boundary faces: owner
-    label faceI = 0;
-    forAllConstIter(faceSet, fSet, iter)
-    {
-        const label faceNo = iter.key();
-
-        if (faceNo < nInternalFaces)
-        {
-            Type t = 0.5*(tf[faceOwner[faceNo]] + tf[faceNeigh[faceNo]]);
-
-            for (direction d=0; d<nComp; ++d)
-            {
-                vec[d] = component(t, d);
-            }
-        }
-        else
-        {
-            const Type& t = tf[faceOwner[faceNo]];
-            for (direction d=0; d<nComp; ++d)
-            {
-                vec[d] = component(t, d);
-            }
-        }
-        vtkOpenFOAMTupleRemap<Type>(vec);
-
-        cellData->InsertTuple(faceI, vec);
-        ++faceI;
-    }
-
-
-    vtkPolyData::SafeDownCast
-    (
-        GetDataSetFromBlock(output, range, datasetNo)
-    )   ->GetCellData()
-        ->AddArray(cellData);
-
-    cellData->Delete();
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C
index deff6377ad73d9a69e48ab0bd5d9f2a760660963..642e5dbaae592ac9adef072a0101e0f09e31c686 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C
@@ -221,7 +221,8 @@ void Foam::vtkPV3Foam::updateInfoLagrangian
 
 void Foam::vtkPV3Foam::updateInfoPatches
 (
-    vtkDataArraySelection* arraySelection
+    vtkDataArraySelection* arraySelection,
+    stringList& enabledEntries
 )
 {
     if (debug)
@@ -230,12 +231,63 @@ void Foam::vtkPV3Foam::updateInfoPatches
             << " [meshPtr=" << (meshPtr_ ? "set" : "NULL") << "]" << endl;
     }
 
+
+    HashSet<string> enabledEntriesSet(enabledEntries);
+
     arrayRangePatches_.reset(arraySelection->GetNumberOfArrays());
 
     int nPatches = 0;
     if (meshPtr_)
     {
         const polyBoundaryMesh& patches = meshPtr_->boundaryMesh();
+        const HashTable<labelList, word>& groups = patches.groupPatchIDs();
+
+        const wordList allPatchNames = patches.names();
+
+        // Add patch groups
+        // ~~~~~~~~~~~~~~~~
+
+        for
+        (
+            HashTable<labelList, word>::const_iterator iter = groups.begin();
+            iter != groups.end();
+            ++iter
+        )
+        {
+            const word& groupName = iter.key();
+            const labelList& patchIDs = iter();
+
+            label nFaces = 0;
+            forAll(patchIDs, i)
+            {
+                nFaces += patches[patchIDs[i]].size();
+            }
+
+            // Valid patch if nFace > 0 - add patch to GUI list
+            if (nFaces)
+            {
+                string vtkGrpName = groupName + " - group";
+                arraySelection->AddArray(vtkGrpName.c_str());
+
+                ++nPatches;
+
+                if (enabledEntriesSet.found(vtkGrpName))
+                {
+                    forAll(patchIDs, i)
+                    {
+                        const polyPatch& pp = patches[patchIDs[i]];
+                        string vtkPatchName = pp.name() + " - patch";
+                        enabledEntriesSet.insert(vtkPatchName);
+                    }
+                    enabledEntriesSet.erase(vtkGrpName);
+                }
+            }
+        }
+
+
+        // Add patches
+        // ~~~~~~~~~~~
+
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
@@ -277,20 +329,101 @@ void Foam::vtkPV3Foam::updateInfoPatches
         {
             polyBoundaryMeshEntries patchEntries(ioObj);
 
-            // Add (non-zero) patches to the list of mesh parts
-            forAll(patchEntries, entryI)
+
+            // Read patches and determine sizes
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            wordList names(patchEntries.size());
+            labelList sizes(patchEntries.size());
+
+            forAll(patchEntries, patchI)
             {
-                label nFaces
-                (
-                    readLabel(patchEntries[entryI].dict().lookup("nFaces"))
-                );
+                const dictionary& patchDict = patchEntries[patchI].dict();
+
+                sizes[patchI] = readLabel(patchDict.lookup("nFaces"));
+                names[patchI] = patchEntries[patchI].keyword();
+            }
+
+
+            // Add (non-zero) patch groups to the list of mesh parts
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            HashTable<labelList, word> groups(patchEntries.size());
+
+            forAll(patchEntries, patchI)
+            {
+                const dictionary& patchDict = patchEntries[patchI].dict();
+
+                wordList groupNames;
+                patchDict.readIfPresent("inGroups", groupNames);
+                forAll(groupNames, groupI)
+                {
+                    HashTable<labelList, word>::iterator iter = groups.find
+                    (
+                        groupNames[groupI]
+                    );
+                    if (iter != groups.end())
+                    {
+                        iter().append(patchI);
+                    }
+                    else
+                    {
+                        groups.insert(groupNames[groupI], labelList(1, patchI));
+                    }
+                }
+            }
+
+            for
+            (
+                HashTable<labelList, word>::const_iterator iter =
+                    groups.begin();
+                iter != groups.end();
+                ++iter
+            )
+            {
+                const word& groupName = iter.key();
+                const labelList& patchIDs = iter();
+
+                label nFaces = 0;
+                forAll(patchIDs, i)
+                {
+                    nFaces += sizes[patchIDs[i]];
+                }
 
                 // Valid patch if nFace > 0 - add patch to GUI list
                 if (nFaces)
+                {
+                    string vtkGrpName = groupName + " - group";
+
+                    arraySelection->AddArray(vtkGrpName.c_str());
+
+                    ++nPatches;
+
+                    if (enabledEntriesSet.found(vtkGrpName))
+                    {
+                        forAll(patchIDs, i)
+                        {
+                            string vtkPatchName =
+                                names[patchIDs[i]] + " - patch";
+                            enabledEntriesSet.insert(vtkPatchName);
+                        }
+                        enabledEntriesSet.erase(vtkGrpName);
+                    }
+                }
+            }
+
+
+            // Add (non-zero) patches to the list of mesh parts
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            forAll(names, patchI)
+            {
+                // Valid patch if nFace > 0 - add patch to GUI list
+                if (sizes[patchI])
                 {
                     arraySelection->AddArray
                     (
-                        (patchEntries[entryI].keyword() + " - patch").c_str()
+                        (names[patchI] + " - patch").c_str()
                     );
 
                     ++nPatches;
@@ -300,6 +433,9 @@ void Foam::vtkPV3Foam::updateInfoPatches
     }
     arrayRangePatches_ += nPatches;
 
+    // Update enabled entries in case of group selection
+    enabledEntries = enabledEntriesSet.toc();
+
     if (debug)
     {
         // just for debug info
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
index a057afb51189132108ea827eecc0af90aa95f619..57e6931c12173c505ac657ab6203b66c8caadd96 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
@@ -271,7 +271,7 @@ void Foam::vtkPV3Foam::convertVolFields
                 arrayRangeFaceSets_,
                 datasetNo,
                 mesh,
-                fSet
+                fSet.toc()
             );
 
             // TODO: points
diff --git a/bin/tools/foamConfigurePaths b/bin/tools/foamConfigurePaths
index cf1093d456fa1c6813fd34cf4a18f56870380191..e027dca9fb99216c3f94199c47c00abb49ddef49 100755
--- a/bin/tools/foamConfigurePaths
+++ b/bin/tools/foamConfigurePaths
@@ -120,11 +120,17 @@ do
     -archOption | --archOption)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
         archOption="$2"
-        # replace WM_ARCH_OPTION=...
-        _inlineSed \
-            etc/bashrc \
-            '/^[^#]/s@WM_ARCH_OPTION=.*@WM_ARCH_OPTION='"$archOption@" \
-             "Replacing WM_ARCH_OPTION setting by '$archOption'"
+        current_archOption=`grep WM_ARCH_OPTION= etc/bashrc | sed "s/export WM_ARCH_OPTION=//"`
+        if [ "$archOption" != "$current_archOption" ]
+        then
+            # replace WM_ARCH_OPTION=...
+            _inlineSed \
+                etc/bashrc \
+                '/^[^#]/s@WM_ARCH_OPTION=.*@WM_ARCH_OPTION='"$archOption@" \
+                 "Replacing WM_ARCH_OPTION setting by '$archOption'"
+        else
+            echo "WM_ARCH_OPTION already set to $archOption"
+        fi
         shift 2
         ;;
     -paraviewInstall | --paraviewInstall)
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
index 6b8d4587b0644286488c7c5e13330238e948adf5..0b2a45c473897ded3db69456b1ed73b6041817dd 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C
@@ -248,15 +248,81 @@ GeometricBoundaryField
             << endl;
     }
 
-    forAll(bmesh_, patchi)
+
+    // Patch or patch-groups. (using non-wild card entries of dictionaries)
+    forAllConstIter(dictionary, dict, iter)
     {
-        if (bmesh_[patchi].type() != emptyPolyPatch::typeName)
+        if (iter().isDict())
         {
-            if
+            const labelList patchIDs = bmesh_.findIndices
             (
-                bmesh_[patchi].type() == cyclicPolyPatch::typeName
-             && !dict.found(bmesh_[patchi].name())
-            )
+                iter().keyword(),
+                true
+            );
+
+            forAll(patchIDs, i)
+            {
+                label patchi = patchIDs[i];
+                this->set
+                (
+                    patchi,
+                    PatchField<Type>::New
+                    (
+                        bmesh_[patchi],
+                        field,
+                        iter().dict()
+                    )
+                );
+            }
+        }
+    }
+
+    // Check for wildcard patch overrides
+    forAll(bmesh_, patchi)
+    {
+        if (!this->set(patchi))
+        {
+            if (bmesh_[patchi].type() == emptyPolyPatch::typeName)
+            {
+                this->set
+                (
+                    patchi,
+                    PatchField<Type>::New
+                    (
+                        emptyPolyPatch::typeName,
+                        bmesh_[patchi],
+                        field
+                    )
+                );
+            }
+            else
+            {
+                bool found = dict.found(bmesh_[patchi].name());
+
+                if (found)
+                {
+                    this->set
+                    (
+                        patchi,
+                        PatchField<Type>::New
+                        (
+                            bmesh_[patchi],
+                            field,
+                            dict.subDict(bmesh_[patchi].name())
+                        )
+                    );
+                }
+            }
+        }
+    }
+
+
+    // Check for any unset patches
+    forAll(bmesh_, patchi)
+    {
+        if (!this->set(patchi))
+        {
+            if (bmesh_[patchi].type() == cyclicPolyPatch::typeName)
             {
                 FatalIOErrorIn
                 (
@@ -274,30 +340,21 @@ GeometricBoundaryField
                     << "Run foamUpgradeCyclics to convert mesh and fields"
                     << " to split cyclics." << exit(FatalIOError);
             }
-
-            this->set
-            (
-                patchi,
-                PatchField<Type>::New
-                (
-                    bmesh_[patchi],
-                    field,
-                    dict.subDict(bmesh_[patchi].name())
-                )
-            );
-        }
-        else
-        {
-            this->set
-            (
-                patchi,
-                PatchField<Type>::New
+            else
+            {
+                FatalIOErrorIn
                 (
-                    emptyPolyPatch::typeName,
-                    bmesh_[patchi],
-                    field
-                )
-            );
+                    "GeometricField<Type, PatchField, GeoMesh>::\n"
+                    "GeometricBoundaryField::GeometricBoundaryField\n"
+                    "(\n"
+                    "    const BoundaryMesh&,\n"
+                    "    const DimensionedField<Type, GeoMesh>&,\n"
+                    "    const dictionary&\n"
+                    ")",
+                    dict
+                )   << "Cannot find patchField entry for "
+                    << bmesh_[patchi].name() << exit(FatalIOError);
+            }
         }
     }
 }
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
index c32428287a51be8bb849a766becdcdd664fd3c5e..10150f7a64b3f363a28d23244f5f3f0be914609c 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/MapGeometricFields.H
@@ -131,17 +131,16 @@ void MapGeometricFields
             );
 
             // Map the patch fields
-            forAll(field.boundaryField(), patchi)
+            typename GeometricField<Type, PatchField, GeoMesh>
+            ::GeometricBoundaryField& bfield = field.boundaryField();
+            forAll(bfield, patchi)
             {
                 // Cannot check sizes for patch fields because of
                 // empty fields in FV and because point fields get their size
                 // from the patch which has already been resized
                 //
 
-                field.boundaryField()[patchi].autoMap
-                (
-                    mapper.boundaryMap()[patchi]
-                );
+                bfield[patchi].autoMap(mapper.boundaryMap()[patchi]);
             }
 
             field.instance() = field.time().timeName();
diff --git a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
index f5baa164d0d7f4bdcfdd68c1408148a3027f0b4e..6137986165792382e8a5871dca8258a7d6adda04 100644
--- a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
+++ b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.C
@@ -32,12 +32,14 @@ Foam::patchIdentifier::patchIdentifier
 (
     const word& name,
     const label index,
-    const word& physicalType
+    const word& physicalType,
+    const wordList& inGroups
 )
 :
     name_(name),
     index_(index),
-    physicalType_(physicalType)
+    physicalType_(physicalType),
+    inGroups_(inGroups)
 {}
 
 
@@ -52,6 +54,7 @@ Foam::patchIdentifier::patchIdentifier
     index_(index)
 {
     dict.readIfPresent("physicalType", physicalType_);
+    dict.readIfPresent("inGroups", inGroups_);
 }
 
 
@@ -63,7 +66,8 @@ Foam::patchIdentifier::patchIdentifier
 :
     name_(p.name_),
     index_(index),
-    physicalType_(p.physicalType_)
+    physicalType_(p.physicalType_),
+    inGroups_(p.inGroups_)
 {}
 
 
@@ -82,6 +86,11 @@ void Foam::patchIdentifier::write(Ostream& os) const
         os.writeKeyword("physicalType") << physicalType_
             << token::END_STATEMENT << nl;
     }
+    if (inGroups_.size())
+    {
+        os.writeKeyword("inGroups") << inGroups_
+            << token::END_STATEMENT << nl;
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
index bd355f175600d3450db161e2d722ea04ee04bfff..cb7d60b7b5b153ba230ccb1b3fe98645ec5d34bd 100644
--- a/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
+++ b/src/OpenFOAM/meshes/Identifiers/patch/patchIdentifier.H
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef patchIdentifier_H
 #define patchIdentifier_H
 
-#include "word.H"
+#include "wordList.H"
 #include "label.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,6 +68,8 @@ class patchIdentifier
         //- Optional physical type
         mutable word physicalType_;
 
+        //- Optional groups patch belongs to
+        wordList inGroups_;
 
 public:
 
@@ -78,7 +80,8 @@ public:
         (
             const word& name,
             const label index,
-            const word& physicalType = word::null
+            const word& physicalType = word::null,
+            const wordList& inGroups = wordList()
         );
 
         //- Construct from dictionary
@@ -139,6 +142,19 @@ public:
             return index_;
         }
 
+        //- Return the optional groups patch belongs to
+        const wordList& inGroups() const
+        {
+            return inGroups_;
+        }
+
+        //- Return the optional groups patch belongs to for modification
+        wordList& inGroups()
+        {
+            return inGroups_;
+        }
+
+
         //- Write patchIdentifier as a dictionary
         void write(Ostream&) const;
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
index 113a9719a025f3c1eefe479063a2b95e3863384e..c8ec20fff2c03ac53301d391839364f9851e7369 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
@@ -58,6 +58,16 @@ Foam::pointBoundaryMesh::pointBoundaryMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::labelList Foam::pointBoundaryMesh::findIndices
+(
+    const keyType& key,
+    const bool usePatchGroups
+) const
+{
+    return mesh()().boundaryMesh().findIndices(key, usePatchGroups);
+}
+
+
 void Foam::pointBoundaryMesh::calcGeometry()
 {
     PstreamBuffers pBufs(Pstream::defaultCommsType);
diff --git a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
index 0e6e2b9016e6b3767206dac551959b875fc01edd..afbc52c93f91d856be9f8262ca47147a3e05610d 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
@@ -96,6 +96,9 @@ public:
             return mesh_;
         }
 
+        //- Find patch indices given a name
+        labelList findIndices(const keyType&, const bool useGroups) const;
+
         //- Correct polyBoundaryMesh after moving points
         void movePoints(const pointField&);
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index 94a1597d96f2a1d648e3f7c4220b97a82df8be3a..ab13b344440ac5b6fe3648b4559dd4ba104eb103 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -151,6 +151,7 @@ void Foam::polyBoundaryMesh::clearAddressing()
 {
     neighbourEdgesPtr_.clear();
     patchIDPtr_.clear();
+    groupPatchIDsPtr_.clear();
 
     forAll(*this, patchI)
     {
@@ -369,6 +370,54 @@ const Foam::labelList& Foam::polyBoundaryMesh::patchID() const
 }
 
 
+const Foam::HashTable<Foam::labelList, Foam::word>&
+Foam::polyBoundaryMesh::groupPatchIDs() const
+{
+    if (!groupPatchIDsPtr_.valid())
+    {
+        groupPatchIDsPtr_.reset(new HashTable<labelList, word>(10));
+        HashTable<labelList, word>& groupPatchIDs = groupPatchIDsPtr_();
+
+        const polyBoundaryMesh& bm = *this;
+
+        forAll(bm, patchI)
+        {
+            const wordList& groups = bm[patchI].inGroups();
+
+            forAll(groups, i)
+            {
+                const word& name = groups[i];
+
+                if (findPatchID(name) != -1)
+                {
+                    WarningIn("polyBoundaryMesh::groupPatchIDs() const")
+                        << "Patch " << bm[patchI].name()
+                        << " specifies a group " << name
+                        << " which is also a patch name."
+                        << " This might give problems later on." << endl;
+                }
+
+
+                HashTable<labelList, word>::iterator iter = groupPatchIDs.find
+                (
+                    name
+                );
+
+                if (iter != groupPatchIDs.end())
+                {
+                    iter().append(patchI);
+                }
+                else
+                {
+                    groupPatchIDs.insert(name, labelList(1, patchI));
+                }
+            }
+        }
+    }
+    return groupPatchIDsPtr_();
+}
+
+
 Foam::wordList Foam::polyBoundaryMesh::names() const
 {
     const polyPatchList& patches = *this;
@@ -414,28 +463,74 @@ Foam::wordList Foam::polyBoundaryMesh::physicalTypes() const
 }
 
 
-Foam::labelList Foam::polyBoundaryMesh::findIndices(const keyType& key) const
+Foam::labelList Foam::polyBoundaryMesh::findIndices
+(
+    const keyType& key,
+    const bool usePatchGroups
+) const
 {
-    labelList indices;
+    DynamicList<label> indices;
 
     if (!key.empty())
     {
         if (key.isPattern())
         {
             indices = findStrings(key, this->names());
+
+            if (usePatchGroups && groupPatchIDs().size())
+            {
+                labelHashSet indexSet(indices);
+
+                const wordList allGroupNames = groupPatchIDs().toc();
+                labelList groupIDs = findStrings(key, allGroupNames);
+                forAll(groupIDs, i)
+                {
+                    const word& grpName = allGroupNames[groupIDs[i]];
+                    const labelList& patchIDs = groupPatchIDs()[grpName];
+                    forAll(patchIDs, j)
+                    {
+                        if (indexSet.insert(patchIDs[j]))
+                        {
+                            indices.append(patchIDs[j]);
+                        }
+                    }
+                }
+            }
         }
         else
         {
-            indices.setSize(this->size());
-            label nFound = 0;
+            // Literal string. Special version of above to avoid
+            // unnecessary memory allocations
+
+            indices.setCapacity(1);
             forAll(*this, i)
             {
                 if (key == operator[](i).name())
                 {
-                    indices[nFound++] = i;
+                    indices.append(i);
+                    break;
+                }
+            }
+
+            if (usePatchGroups && groupPatchIDs().size())
+            {
+                const HashTable<labelList, word>::const_iterator iter =
+                    groupPatchIDs().find(key);
+
+                if (iter != groupPatchIDs().end())
+                {
+                    labelHashSet indexSet(indices);
+
+                    const labelList& patchIDs = iter();
+                    forAll(patchIDs, j)
+                    {
+                        if (indexSet.insert(patchIDs[j]))
+                        {
+                            indices.append(patchIDs[j]);
+                        }
+                    }
                 }
             }
-            indices.setSize(nFound);
         }
     }
 
@@ -552,7 +647,8 @@ Foam::label Foam::polyBoundaryMesh::whichPatch(const label faceIndex) const
 Foam::labelHashSet Foam::polyBoundaryMesh::patchSet
 (
     const UList<wordRe>& patchNames,
-    const bool warnNotFound
+    const bool warnNotFound,
+    const bool usePatchGroups
 ) const
 {
     const wordList allPatchNames(this->names());
@@ -560,23 +656,57 @@ Foam::labelHashSet Foam::polyBoundaryMesh::patchSet
 
     forAll(patchNames, i)
     {
+         const word& patchName = patchNames[i];
+
         // Treat the given patch names as wild-cards and search the set
         // of all patch names for matches
-        labelList patchIDs = findStrings(patchNames[i], allPatchNames);
+        labelList patchIDs = findStrings(patchName, allPatchNames);
 
-        if (patchIDs.empty() && warnNotFound)
+        forAll(patchIDs, j)
         {
-            WarningIn
-            (
-                "polyBoundaryMesh::patchSet"
-                "(const wordReList&, const bool) const"
-            )   << "Cannot find any patch names matching " << patchNames[i]
-                << endl;
+            ids.insert(patchIDs[j]);
         }
 
-        forAll(patchIDs, j)
+        if (patchIDs.empty())
         {
-            ids.insert(patchIDs[j]);
+            if (usePatchGroups)
+            {
+                const wordList allGroupNames = groupPatchIDs().toc();
+
+                // Regard as group name
+                labelList groupIDs = findStrings(patchName, allGroupNames);
+
+                forAll(groupIDs, i)
+                {
+                    const word& name = allGroupNames[groupIDs[i]];
+                    const labelList& extraPatchIDs = groupPatchIDs()[name];
+
+                    forAll(extraPatchIDs, extraI)
+                    {
+                        ids.insert(extraPatchIDs[extraI]);
+                    }
+                }
+
+                if (groupIDs.empty() && warnNotFound)
+                {
+                    WarningIn
+                    (
+                        "polyBoundaryMesh::patchSet"
+                        "(const wordReList&, const bool, const bool) const"
+                    )   << "Cannot find any patch or group names matching "
+                        << patchName
+                        << endl;
+                }
+            }
+            else if (warnNotFound)
+            {
+                WarningIn
+                (
+                    "polyBoundaryMesh::patchSet"
+                    "(const wordReList&, const bool, const bool) const"
+                )   << "Cannot find any patch names matching " << patchName
+                    << endl;
+            }
         }
     }
 
@@ -778,6 +908,7 @@ void Foam::polyBoundaryMesh::updateMesh()
 {
     neighbourEdgesPtr_.clear();
     patchIDPtr_.clear();
+    groupPatchIDsPtr_.clear();
 
     PstreamBuffers pBufs(Pstream::defaultCommsType);
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 8cc837dc78612b3b56c0430b062085883224a779..211d9148e941d6a6a799f149265dcecaadea0dae 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -70,6 +70,8 @@ class polyBoundaryMesh
 
         mutable autoPtr<labelList> patchIDPtr_;
 
+        mutable autoPtr<HashTable<labelList, word> > groupPatchIDsPtr_;
+
         //- Edges of neighbouring patches
         mutable autoPtr<List<labelPairList> > neighbourEdgesPtr_;
 
@@ -153,8 +155,12 @@ public:
         //- Return a list of physical types
         wordList physicalTypes() const;
 
-        //- Return patch indices for all matches
-        labelList findIndices(const keyType&) const;
+        //- Return patch indices for all matches. Optionally matches patchGroups
+        labelList findIndices
+        (
+            const keyType&,
+            const bool usePatchGroups = false
+        ) const;
 
         //- Return patch index for the first match, return -1 if not found
         label findIndex(const keyType&) const;
@@ -168,12 +174,17 @@ public:
         //- Per boundary face label the patch index
         const labelList& patchID() const;
 
+        //- Per patch group the patch indices
+        const HashTable<labelList, word>& groupPatchIDs() const;
+
         //- Return the set of patch IDs corresponding to the given names
-        //  By default warns if given names are not found.
+        //  By default warns if given names are not found. Optionally
+        //  matches to patchGroups as well as patchNames
         labelHashSet patchSet
         (
             const UList<wordRe>& patchNames,
-            const bool warnNotFound = true
+            const bool warnNotFound = true,
+            const bool usePatchGroups = false
         ) const;
 
         //- Check whether all procs have all patches and in same order. Return
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index cbc223c5833fad72c22a4c5294b522aa948d3675..960e35ddd4ca620744ee5733e4843b9f539c139b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -703,11 +703,11 @@ void Foam::polyMesh::resetPrimitives
     {
         boundary_[patchI] = polyPatch
         (
-            boundary_[patchI].name(),
-            patchSizes[patchI],
-            patchStarts[patchI],
+            boundary_[patchI],
+            boundary_,
             patchI,
-            boundary_
+            patchSizes[patchI],
+            patchStarts[patchI]
         );
     }
 
diff --git a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
index 44f1e79ca8c9ef9d5adafbb2d7ec9dcbabfb902c..782f8456c4abfd6ed12b0a356dfa3219a2e7bb0d 100644
--- a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
+++ b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
@@ -104,16 +104,20 @@ Foam::tmp<Foam::volScalarField> Foam::consumptionSpeed::omega0Sigma
 
     volScalarField& omega0 = tomega0();
 
-    forAll(omega0, celli)
+    volScalarField::InternalField& iomega0 = omega0.internalField();
+
+    forAll(iomega0, celli)
     {
-        omega0[celli] = omega0Sigma(sigma[celli], 1.0);
+        iomega0[celli] = omega0Sigma(sigma[celli], 1.0);
     }
 
-    forAll(omega0.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& bomega0 = omega0.boundaryField();
+
+    forAll(bomega0, patchi)
     {
-        forAll(omega0.boundaryField()[patchi], facei)
+        forAll(bomega0[patchi], facei)
         {
-            omega0.boundaryField()[patchi][facei] =
+            bomega0[patchi][facei] =
                 omega0Sigma
                 (
                     sigma.boundaryField()[patchi][facei],
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 0af332a336de54520688e2108e15e03c61dd7b7d..3b295f46dd8401afe885ccacaa9a39e7c340db22 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -360,9 +360,11 @@ Foam::dynamicRefineFvMesh::refine
             }
 
             // Recalculate new boundary faces.
-            forAll(phi.boundaryField(), patchI)
+            surfaceScalarField::GeometricBoundaryField& bphi =
+                phi.boundaryField();
+            forAll(bphi, patchI)
             {
-                fvsPatchScalarField& patchPhi = phi.boundaryField()[patchI];
+                fvsPatchScalarField& patchPhi = bphi[patchI];
                 const fvsPatchScalarField& patchPhiU =
                     phiU.boundaryField()[patchI];
 
@@ -404,8 +406,7 @@ Foam::dynamicRefineFvMesh::refine
                     const fvsPatchScalarField& patchPhiU =
                         phiU.boundaryField()[patchI];
 
-                    fvsPatchScalarField& patchPhi =
-                        phi.boundaryField()[patchI];
+                    fvsPatchScalarField& patchPhi = bphi[patchI];
 
                     patchPhi[i] = patchPhiU[i];
                 }
@@ -549,6 +550,9 @@ Foam::dynamicRefineFvMesh::unrefine
             }
 
             surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
+            surfaceScalarField::GeometricBoundaryField& bphi =
+                phi.boundaryField();
+
             const surfaceScalarField phiU
             (
                 fvc::interpolate
@@ -582,10 +586,7 @@ Foam::dynamicRefineFvMesh::unrefine
 
                             const fvsPatchScalarField& patchPhiU =
                                 phiU.boundaryField()[patchI];
-
-                            fvsPatchScalarField& patchPhi =
-                                phi.boundaryField()[patchI];
-
+                            fvsPatchScalarField& patchPhi = bphi[patchI];
                             patchPhi[i] = patchPhiU[i];
                         }
                     }
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
index f841e2b74dcb5aa27710b7fdc8710d712f0864de..bb01e2fa833baa6ab36826620b9f1e27cf24279f 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
@@ -81,6 +81,9 @@ void Foam::fvMeshAdder::MapVolField
     // Patch fields from old mesh
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    typename GeometricField<Type, fvPatchField, volMesh>::
+    GeometricBoundaryField& bfld = fld.boundaryField();
+
     {
         const labelList& oldPatchMap = meshMap.oldPatchMap();
         const labelList& oldPatchStarts = meshMap.oldPatchStarts();
@@ -121,18 +124,18 @@ void Foam::fvMeshAdder::MapVolField
 
 
         // Sort deleted ones last so is now in newPatch ordering
-        fld.boundaryField().reorder(oldToNew);
+        bfld.reorder(oldToNew);
         // Extend to covers all patches
-        fld.boundaryField().setSize(mesh.boundaryMesh().size());
+        bfld.setSize(mesh.boundaryMesh().size());
         // Delete unused patches
         for
         (
             label newPatchI = nUsedPatches;
-            newPatchI < fld.boundaryField().size();
+            newPatchI < bfld.size();
             newPatchI++
         )
         {
-            fld.boundaryField().set(newPatchI, NULL);
+            bfld.set(newPatchI, NULL);
         }
 
 
@@ -167,12 +170,12 @@ void Foam::fvMeshAdder::MapVolField
                 //   value
                 // - hope that field mapping allows aliasing since old and new
                 //   are same memory!
-                fld.boundaryField().set
+                bfld.set
                 (
                     newPatchI,
                     fvPatchField<Type>::New
                     (
-                        fld.boundaryField()[newPatchI], // old field
+                        bfld[newPatchI],                // old field
                         mesh.boundary()[newPatchI],     // new fvPatch
                         fld.dimensionedInternalField(), // new internal field
                         patchMapper                     // mapper (new to old)
@@ -201,7 +204,7 @@ void Foam::fvMeshAdder::MapVolField
                 const polyPatch& oldPatch =
                     fldToAdd.mesh().boundaryMesh()[patchI];
 
-                if (!fld.boundaryField()(newPatchI))
+                if (!bfld(newPatchI))
                 {
                     // First occurrence of newPatchI. Map from existing
                     // patchField
@@ -221,7 +224,7 @@ void Foam::fvMeshAdder::MapVolField
 
                     directFvPatchFieldMapper patchMapper(newToAdded);
 
-                    fld.boundaryField().set
+                    bfld.set
                     (
                         newPatchI,
                         fvPatchField<Type>::New
@@ -255,7 +258,7 @@ void Foam::fvMeshAdder::MapVolField
                     const fvPatchField<Type>& addedFld =
                         fldToAdd.boundaryField()[patchI];
 
-                    fvPatchField<Type>& newFld = fld.boundaryField()[newPatchI];
+                    fvPatchField<Type>& newFld = bfld[newPatchI];
 
                     forAll(newFld, i)
                     {
@@ -357,6 +360,9 @@ void Foam::fvMeshAdder::MapSurfaceField
     const fvMesh& mesh = fld.mesh();
     const labelList& oldPatchStarts = meshMap.oldPatchStarts();
 
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+    GeometricBoundaryField& bfld = fld.boundaryField();
+
     // Internal field
     // ~~~~~~~~~~~~~~
 
@@ -375,9 +381,9 @@ void Foam::fvMeshAdder::MapSurfaceField
         // Faces that were boundary faces but are not anymore.
         // Use owner value (so lowest numbered cell, i.e. from 'old' not 'added'
         // mesh)
-        forAll(fld.boundaryField(), patchI)
+        forAll(bfld, patchI)
         {
-            const fvsPatchField<Type>& pf = fld.boundaryField()[patchI];
+            const fvsPatchField<Type>& pf = bfld[patchI];
 
             label start = oldPatchStarts[patchI];
 
@@ -436,18 +442,18 @@ void Foam::fvMeshAdder::MapSurfaceField
 
 
         // Sort deleted ones last so is now in newPatch ordering
-        fld.boundaryField().reorder(oldToNew);
+        bfld.reorder(oldToNew);
         // Extend to covers all patches
-        fld.boundaryField().setSize(mesh.boundaryMesh().size());
+        bfld.setSize(mesh.boundaryMesh().size());
         // Delete unused patches
         for
         (
             label newPatchI = nUsedPatches;
-            newPatchI < fld.boundaryField().size();
+            newPatchI < bfld.size();
             newPatchI++
         )
         {
-            fld.boundaryField().set(newPatchI, NULL);
+            bfld.set(newPatchI, NULL);
         }
 
 
@@ -478,16 +484,16 @@ void Foam::fvMeshAdder::MapSurfaceField
                 // Create new patchField with same type as existing one.
                 // Note:
                 // - boundaryField already in new order so access with newPatchI
-                // - fld.boundaryField()[newPatchI] both used for type and old
+                // - bfld[newPatchI] both used for type and old
                 //   value
                 // - hope that field mapping allows aliasing since old and new
                 //   are same memory!
-                fld.boundaryField().set
+                bfld.set
                 (
                     newPatchI,
                     fvsPatchField<Type>::New
                     (
-                        fld.boundaryField()[newPatchI], // old field
+                        bfld[newPatchI],                // old field
                         mesh.boundary()[newPatchI],     // new fvPatch
                         fld.dimensionedInternalField(), // new internal field
                         patchMapper                     // mapper (new to old)
@@ -516,7 +522,7 @@ void Foam::fvMeshAdder::MapSurfaceField
                 const polyPatch& oldPatch =
                     fldToAdd.mesh().boundaryMesh()[patchI];
 
-                if (!fld.boundaryField()(newPatchI))
+                if (!bfld(newPatchI))
                 {
                     // First occurrence of newPatchI. Map from existing
                     // patchField
@@ -536,7 +542,7 @@ void Foam::fvMeshAdder::MapSurfaceField
 
                     directFvPatchFieldMapper patchMapper(newToAdded);
 
-                    fld.boundaryField().set
+                    bfld.set
                     (
                         newPatchI,
                         fvsPatchField<Type>::New
@@ -570,8 +576,7 @@ void Foam::fvMeshAdder::MapSurfaceField
                     const fvsPatchField<Type>& addedFld =
                         fldToAdd.boundaryField()[patchI];
 
-                    fvsPatchField<Type>& newFld =
-                        fld.boundaryField()[newPatchI];
+                    fvsPatchField<Type>& newFld = bfld[newPatchI];
 
                     forAll(newFld, i)
                     {
diff --git a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
index 7c72beaf99d5a7c4fea84146b544b018de48d2f3..206e7c4e95e5da9ce5e12b583d9a8e8175f12330 100644
--- a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
+++ b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
@@ -47,10 +47,12 @@ bool Foam::adjustPhi
         scalar fixedMassOut = 0.0;
         scalar adjustableMassOut = 0.0;
 
-        forAll(phi.boundaryField(), patchi)
+        surfaceScalarField::GeometricBoundaryField& bphi = phi.boundaryField();
+
+        forAll(bphi, patchi)
         {
             const fvPatchVectorField& Up = U.boundaryField()[patchi];
-            const fvsPatchScalarField& phip = phi.boundaryField()[patchi];
+            const fvsPatchScalarField& phip = bphi[patchi];
 
             if (!isA<processorFvsPatchScalarField>(phip))
             {
@@ -123,10 +125,10 @@ bool Foam::adjustPhi
                 << exit(FatalError);
         }
 
-        forAll(phi.boundaryField(), patchi)
+        forAll(bphi, patchi)
         {
             const fvPatchVectorField& Up = U.boundaryField()[patchi];
-            fvsPatchScalarField& phip = phi.boundaryField()[patchi];
+            fvsPatchScalarField& phip = bphi[patchi];
 
             if (!isA<processorFvsPatchScalarField>(phip))
             {
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index e92de99aec8f7d15ad4e83c945a3a479d84bc063..d6ef3989cdd5d0b7ebf570d88c9e981fb109c5fb 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -75,7 +75,10 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
     }
 
-    forAll(corDeltaT.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& bcorDeltaT =
+        corDeltaT.boundaryField();
+
+    forAll(bcorDeltaT, patchi)
     {
         const fvsPatchScalarField& pcofrDeltaT =
             cofrDeltaT.boundaryField()[patchi];
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
index 91c67d823a7e4cbd4be3ec22d63d7115cd9ffd1a..2756f56de5696de3642269ac3d22276f4e9d9d83 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcAverage.C
@@ -73,9 +73,12 @@ average
         surfaceSum(mesh.magSf()*ssf)/surfaceSum(mesh.magSf())
     )().internalField();
 
-    forAll(av.boundaryField(), patchi)
+    typename GeometricField<Type, fvPatchField, volMesh>::
+    GeometricBoundaryField& bav = av.boundaryField();
+
+    forAll(bav, patchi)
     {
-        av.boundaryField()[patchi] = ssf.boundaryField()[patchi];
+        bav[patchi] = ssf.boundaryField()[patchi];
     }
 
     av.correctBoundaryConditions();
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
index 4bf1a4dd6b50be8dcc9b24fd2910f1ea749c0aca..ebec7bb7d5ac45099c4d767f6626f3fcce89927e 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
@@ -130,9 +130,11 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
 
     // Visit the boundaries. Coupled boundaries are taken into account
     // in the construction of d vectors.
-    forAll(lsP.boundaryField(), patchi)
+    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
+
+    forAll(blsP, patchi)
     {
-        const fvPatch& p = lsP.boundaryField()[patchi].patch();
+        const fvPatch& p = blsP[patchi].patch();
         const labelUList& faceCells = p.faceCells();
 
         // Build the d-vectors
@@ -242,9 +244,9 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
             ((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d);
     }
 
-    forAll(lsP.boundaryField(), patchI)
+    forAll(blsP, patchI)
     {
-        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
+        fvsPatchVectorField& patchLsP = blsP[patchI];
 
         const fvPatch& p = patchLsP.patch();
         const labelUList& faceCells = p.faceCells();
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
index 27f5b72fcfd483336cc71e2c4549a92f3da24b38..86b4bc0ad314e4380f51b18aa1199258e2a7e93c 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
@@ -121,7 +121,9 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
     }
 
 
-    forAll(lsP.boundaryField(), patchi)
+    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
+
+    forAll(blsP, patchi)
     {
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
@@ -172,9 +174,9 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
     }
 
-    forAll(lsP.boundaryField(), patchi)
+    forAll(blsP, patchi)
     {
-        fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi];
+        fvsPatchVectorField& patchLsP = blsP[patchi];
 
         const fvsPatchScalarField& pw = w.boundaryField()[patchi];
         const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index 8cf43b1df9433a2263406b3dfced9f1f1e74b340..6f5b58c3e5d564a5b566250c76ba684de5065224 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -125,13 +125,6 @@ void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
                 }
             }
         }
-        else
-        {
-            for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
-            {
-                phiPsiCorrs[phasei][facei] = 0;
-            }
-        }
     }
 }
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 55478e59aa4ab74e9b7c540c62473eef53a65cae..98ed53686089aab9add585aabec71ef18273370e 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -668,9 +668,12 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
         limitSum(phiPsiCorrsInternal);
     }
 
-    forAll(phiPsiCorrs[0].boundaryField(), patchi)
+    surfaceScalarField::GeometricBoundaryField& bfld =
+        phiPsiCorrs[0].boundaryField();
+
+    forAll(bfld, patchi)
     {
-        if (phiPsiCorrs[0].boundaryField()[patchi].coupled())
+        if (bfld[patchi].coupled())
         {
             UPtrList<scalarField> phiPsiCorrsPatch(phiPsiCorrs.size());
             forAll(phiPsiCorrs, phasei)
diff --git a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C
index 5f01380527a3ff54ff0dbf0d74ffe11c45b0b36d..035ab4261534749a5263c3fbea961547be6b1091 100644
--- a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C
+++ b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.C
@@ -23,8 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvMesh.H"
 #include "fvBoundaryMesh.H"
+#include "fvMesh.H"
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -87,6 +87,16 @@ Foam::label Foam::fvBoundaryMesh::findPatchID(const word& patchName) const
 }
 
 
+Foam::labelList Foam::fvBoundaryMesh::findIndices
+(
+    const keyType& key,
+    const bool usePatchGroups
+) const
+{
+    return mesh().boundaryMesh().findIndices(key, usePatchGroups);
+}
+
+
 void Foam::fvBoundaryMesh::movePoints()
 {
     forAll(*this, patchI)
diff --git a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
index 8270fc9348ec30355dd04cd951f45b4ecd62af0f..df9616855150870848d5fc40798ecc1884c9956f 100644
--- a/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
+++ b/src/finiteVolume/fvMesh/fvBoundaryMesh/fvBoundaryMesh.H
@@ -109,6 +109,9 @@ public:
             //- Find patch index given a name
             label findPatchID(const word& patchName) const;
 
+            //- Find patch indices given a name
+            labelList findIndices(const keyType&, const bool useGroups) const;
+
 
         //- Correct patches after moving points
         void movePoints();
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
index 215d6bbef8fc28f3ae36ce0d91ed8e5b39f03b04..a19febe07b637ad81759ec08044641b6cb002772 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H
@@ -76,18 +76,18 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalVelocity.internalField(), tV)
+    volVectorField::InternalField& itotalVelocity =
+        totalVelocity.internalField();
+
+    forAll(itotalVelocity, tV)
     {
         if (totalMass_sum[tV] > VSMALL)
         {
-            totalVelocity.internalField()[tV] =
-                totalMomentum_sum[tV]
-                /totalMass_sum[tV];
+            itotalVelocity[tV] = totalMomentum_sum[tV]/totalMass_sum[tV];
         }
         else
         {
-            totalVelocity.internalField()[tV] =
-            vector::zero;
+            itotalVelocity[tV] = vector::zero;
         }
     }
 
@@ -150,17 +150,20 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalTemperature.internalField(), tT)
+    volScalarField::InternalField& itotalTemperature =
+        totalTemperature.internalField();
+
+    forAll(itotalTemperature, tT)
     {
         if (totalN_sum[tT] > 0)
         {
-            totalTemperature.internalField()[tT] =
+            itotalTemperature[tT] =
                 totalTemperatureVTerms_sum[tT]
                 /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
         }
         else
         {
-            totalTemperature.internalField()[tT] = 0.0;
+            itotalTemperature[tT] = 0.0;
         }
     }
 
@@ -206,17 +209,19 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalMeanKE.internalField(), tMKE)
+    volScalarField::InternalField& itotalMeanKE = totalMeanKE.internalField();
+
+    forAll(itotalMeanKE, tMKE)
     {
         if (totalN_sum[tMKE] > 0)
         {
-            totalMeanKE.internalField()[tMKE] =
+            itotalMeanKE[tMKE] =
                 totalKE_sum[tMKE]
                 /totalN_sum[tMKE];
         }
         else
         {
-            totalMeanKE.internalField()[tMKE] = 0.0;
+            itotalMeanKE[tMKE] = 0.0;
         }
     }
 }
diff --git a/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H b/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
index 5a195aa9389a1a1cd1365f5b1bb0449a5aaf20e4..a19febe07b637ad81759ec08044641b6cb002772 100644
--- a/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
+++ b/src/lagrangian/molecularDynamics/old/mdTools/averageMDFields.H
@@ -76,18 +76,18 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalVelocity.internalField(), tV)
+    volVectorField::InternalField& itotalVelocity =
+        totalVelocity.internalField();
+
+    forAll(itotalVelocity, tV)
     {
         if (totalMass_sum[tV] > VSMALL)
         {
-            totalVelocity.internalField()[tV] =
-                totalMomentum_sum[tV]
-                /totalMass_sum[tV];
+            itotalVelocity[tV] = totalMomentum_sum[tV]/totalMass_sum[tV];
         }
         else
         {
-            totalVelocity.internalField()[tV] =
-            vector::zero;
+            itotalVelocity[tV] = vector::zero;
         }
     }
 
@@ -150,17 +150,20 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalTemperature.internalField(), tT)
+    volScalarField::InternalField& itotalTemperature =
+        totalTemperature.internalField();
+
+    forAll(itotalTemperature, tT)
     {
         if (totalN_sum[tT] > 0)
         {
-            totalTemperature.internalField()[tT] =
+            itotalTemperature[tT] =
                 totalTemperatureVTerms_sum[tT]
                 /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
         }
         else
         {
-            totalTemperature.internalField()[tT] = 0.0;
+            itotalTemperature[tT] = 0.0;
         }
     }
 
@@ -206,18 +209,19 @@ if (runTime.outputTime())
         }
     }
 
-    forAll(totalMeanKE.internalField(), tMKE)
+    volScalarField::InternalField& itotalMeanKE = totalMeanKE.internalField();
+
+    forAll(itotalMeanKE, tMKE)
     {
         if (totalN_sum[tMKE] > 0)
         {
-            totalMeanKE.internalField()[tMKE] =
+            itotalMeanKE[tMKE] =
                 totalKE_sum[tMKE]
                 /totalN_sum[tMKE];
         }
         else
         {
-            totalMeanKE.internalField()[tMKE] = 0.0;
+            itotalMeanKE[tMKE] = 0.0;
         }
     }
 }
-
diff --git a/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C b/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
index 912564714a53adcf690f29232ddd2723c412d719..024f6cfcedb8df8391a9742028f9e4cecabbe24d 100644
--- a/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
+++ b/src/thermophysicalModels/basicSolidThermo/solidMixtureThermo/solidMixtureThermo/solidMixtureThermo.C
@@ -39,16 +39,20 @@ void Foam::solidMixtureThermo<MixtureType>::calculate()
     scalarField& sigmaSCells = sigmaS_.internalField();
     scalarField& emissivityCells = emissivity_.internalField();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
-        rhoCells[celli] = MixtureType::rho(T_[celli], celli);
-        kappaCells[celli] = MixtureType::kappa(T_[celli], celli);
-        sigmaSCells[celli] = MixtureType::sigmaS(T_[celli], celli);
-        KCells[celli] = MixtureType::K(T_[celli], celli);
-        emissivityCells[celli] = MixtureType::emissivity(T_[celli], celli);
+        rhoCells[celli] = MixtureType::rho(iT[celli], celli);
+        kappaCells[celli] = MixtureType::kappa(iT[celli], celli);
+        sigmaSCells[celli] = MixtureType::sigmaS(iT[celli], celli);
+        KCells[celli] = MixtureType::K(iT[celli], celli);
+        emissivityCells[celli] = MixtureType::emissivity(iT[celli], celli);
     }
 
-    forAll(T_.boundaryField(), patchI)
+    const volScalarField::GeometricBoundaryField& bT = T_.boundaryField();
+
+    forAll(bT, patchI)
     {
         rho_.boundaryField()[patchI] == this->rho(patchI)();
         K_.boundaryField()[patchI] == this->K(patchI)();
@@ -158,14 +162,18 @@ Foam::solidMixtureThermo<MixtureType>::Cp() const
     );
     volScalarField& Cp = tCp();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         Cp[celli] = MixtureType::Cp(T_[celli], celli);
     }
 
-    forAll(Cp.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bCp = Cp.boundaryField();
+
+    forAll(bCp, patchI)
     {
-        Cp.boundaryField()[patchI] == this->Cp(patchI)();
+        bCp[patchI] == this->Cp(patchI)();
     }
 
     return tCp;
@@ -194,14 +202,18 @@ Foam::solidMixtureThermo<MixtureType>::hs() const
     );
     volScalarField& hs = ths();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         hs[celli] = MixtureType::hs(T_[celli], celli);
     }
 
-    forAll(hs.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bHs = hs.boundaryField();
+
+    forAll(bHs, patchI)
     {
-        hs.boundaryField()[patchI] == this->hs(patchI)();
+        bHs[patchI] == this->hs(patchI)();
     }
 
     return ths;
@@ -230,14 +242,18 @@ Foam::solidMixtureThermo<MixtureType>::Hf() const
     );
     volScalarField& hf = thF();
 
-    forAll(T_.internalField(), celli)
+    const volScalarField::InternalField& iT = T_.internalField();
+
+    forAll(iT, celli)
     {
         hf[celli] = MixtureType::hf(T_[celli], celli);
     }
 
-    forAll(hf.boundaryField(), patchI)
+    volScalarField::GeometricBoundaryField& bhf = hf.boundaryField();
+
+    forAll(bhf, patchI)
     {
-        hf.boundaryField()[patchI] == this->Hf(patchI)();
+        bhf[patchI] == this->Hf(patchI)();
     }
 
     return thF;
diff --git a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
index 1c5cc36d09d590df6c749f91bc92036c9f8ed761..6c68214b63f230de8cd9a02b92a3fe2e6302f8d1 100644
--- a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
@@ -760,11 +760,11 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
         )
     );
 
-    volScalarField& gasHs = tHs();
+    volScalarField::InternalField& gasHs = tHs().internalField();
 
     const GasThermo& mixture = gasThermo_[index];
 
-    forAll(gasHs.internalField(), cellI)
+    forAll(gasHs, cellI)
     {
         gasHs[cellI] = mixture.Hs(T[cellI]);
     }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
index 03db056e8b26e67fb249df01a89d3715cba994bf..2883377ba97ece4f643fab08878e8db7bad8802c 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes
@@ -32,12 +32,12 @@ divSchemes
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,Yi_hs)  Gauss multivariateSelection
     {
-        O2              limitedLinear01 1;
-        N2              limitedLinear01 1;
-        C3H8            limitedLinear01 1;
-        H2O             limitedLinear01 1;
-        CO2             limitedLinear01 1;
-        hs              limitedLinear 1;
+        O2              linearUpwind grad(O2);
+        N2              linearUpwind grad(N2);
+        C3H8            linearUpwind grad(C3H8);
+        H2O             linearUpwind grad(H2O);
+        CO2             linearUpwind grad(CO2);
+        hs              linearUpwind grad(hs);
     };
     div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,K)          Gauss limitedLinear 1;
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
index 058c08def36d93640ad3ebab96aff71302a74d33..ce35cbafd7f8f50e49222d6ed42831fefccde7ef 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
@@ -15,61 +15,117 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-7
+9
 (
     floor
     {
         type            wall;
         nFaces          50;
-        startFace       3896;
+        startFace       3882;
     }
     ceiling
     {
         type            wall;
         nFaces          50;
-        startFace       3946;
+        startFace       3932;
     }
     inlet
     {
         type            patch;
         nFaces          40;
-        startFace       3996;
+        startFace       3982;
     }
     outlet
     {
         type            patch;
         nFaces          40;
-        startFace       4036;
+        startFace       4022;
     }
     fixedWalls
     {
         type            empty;
         nFaces          4000;
-        startFace       4076;
+        startFace       4062;
     }
     baffle1Wall_0
     {
-        type            directMappedWall;
+        type            mappedWall;
         nFaces          14;
-        startFace       8076;
+        startFace       8062;
         sampleMode      nearestPatchFace;
         sampleRegion    region0;
         samplePatch     baffle1Wall_1;
         offsetMode      uniform;
-        offset          ( 0 0 0 );
-        faces           ( );
+        offset          (0 0 0);
     }
     baffle1Wall_1
     {
-        type            directMappedWall;
+        type            mappedWall;
         nFaces          14;
-        startFace       8090;
+        startFace       8076;
         sampleMode      nearestPatchFace;
         sampleRegion    region0;
         samplePatch     baffle1Wall_0;
         offsetMode      uniform;
-        offset          ( 0 0 0 );
-        faces           ( );
+        offset          (0 0 0);
+    }
+    region0_to_baffleRegion_baffleFaces2_top
+    {
+        type            mappedWall;
+        nFaces          14;
+        startFace       8090;
+        sampleMode      nearestPatchFace;
+        sampleRegion    baffleRegion;
+        samplePatch     region0_to_baffleRegion_baffleFaces2_top;
+        offsetMode      nonuniform;
+        offsets         
+14
+(
+(0.02 8.67362e-19 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -6.93889e-18 -0)
+(0.02 -6.93889e-18 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -0 -0)
+(0.02 -2.77556e-17 -0)
+)
+;
+    }
+    region0_to_baffleRegion_baffleFaces2_bottom
+    {
+        type            mappedWall;
+        nFaces          14;
+        startFace       8104;
+        sampleMode      nearestPatchFace;
+        sampleRegion    baffleRegion;
+        samplePatch     region0_to_baffleRegion_baffleFaces2_bottom;
+        offsetMode      nonuniform;
+        offsets         
+14
+(
+(1.11022e-16 8.67362e-19 -0)
+(1.11022e-16 -3.46945e-18 -0)
+(-1.11022e-16 -0 -0)
+(-1.11022e-16 -6.93889e-18 1.11022e-16)
+(-1.11022e-16 -1.38778e-17 -0)
+(-1.11022e-16 1.38778e-17 -0)
+(-1.11022e-16 -0 -0)
+(-0 -0 -0)
+(-1.11022e-16 -0 -0)
+(-1.11022e-16 -0 -0)
+(-0 -0 -0)
+(-0 -2.77556e-17 -0)
+(-0 -0 -0)
+(-0 -2.77556e-17 -0)
+)
+;
     }
 )
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index ca4ad66e3319d604128f1fac908b57acc69e10ee..68819f01c42149e557daa31e094cc4fb84192297 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -113,6 +113,14 @@ castellatedMeshControls
         {
             // Surface-wise min and max refinement level
             level (5 6);
+
+            // Optional specification of patch type (default is wall). No
+            // constraint types (cyclic, symmetry) etc. are allowed.
+            patchInfo
+            {
+                type wall;
+                inGroups (motorBike);
+            }
         }
     }
 
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary b/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
index 7b72e85f595d014f15070eef4db4e25d5088980a..7abdd81f4fef7404693356fc268d5daa1c6fff5d 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
@@ -15,37 +15,43 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-5
+6
 (
     outlet
     {
         type            patch;
-        nFaces          600;
-        startFace       51900;
+        nFaces          922;
+        startFace       364825;
     }
     sides
     {
         type            patch;
-        nFaces          1200;
-        startFace       52500;
+        nFaces          1834;
+        startFace       365747;
     }
     inlet
     {
         type            patch;
-        nFaces          600;
-        startFace       53700;
+        nFaces          923;
+        startFace       367581;
     }
     ground
     {
         type            wall;
-        nFaces          900;
-        startFace       54300;
+        nFaces          0;
+        startFace       368504;
     }
     top
     {
         type            patch;
         nFaces          900;
-        startFace       55200;
+        startFace       368504;
+    }
+    terrain_patch0
+    {
+        type            wall;
+        nFaces          16037;
+        startFace       369404;
     }
 )
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/transportProperties b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/transportProperties
index 1f390a2c972b62b7f99a536c9fe82b081ffacdda..22545b07217af1c36a8388fb31ebf106c6e7cb27 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/transportProperties
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/constant/transportProperties
@@ -71,16 +71,19 @@ drag
         air
         {
             type SchillerNaumann;
-            residualDrag 0;
+            residualPhaseFraction 0;
+            residualSlip 0;
         }
 
         water
         {
             type SchillerNaumann;
-            residualDrag 0;
+            residualPhaseFraction 0;
+            residualSlip 0;
         }
 
-        residualDrag 1;
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
     }
 );
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/constant/transportProperties b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/constant/transportProperties
index d0fdbfd12bce6c60b99b697cc89c1612f5ffdb38..e59dbd8d6fad3c2fc58faa3a466461138fb9a9bd 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/constant/transportProperties
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/constant/transportProperties
@@ -79,9 +79,9 @@ sigmas
     (air water)     0.07
     (air oil)       0.07
     (air mercury)   0.07
-    (water oil)     0.07
-    (water mercury) 0.07
-    (oil mercury)   0.07
+    (water oil)     0
+    (water mercury) 0
+    (oil mercury)   0
 );
 
 interfaceCompression
@@ -109,19 +109,22 @@ drag
     (air water)
     {
         type SchillerNaumann;
-        residualDrag 100;
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
     }
 
     (air oil)
     {
         type SchillerNaumann;
-        residualDrag 100;
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
     }
 
     (air mercury)
     {
         type SchillerNaumann;
-        residualDrag 100;
+        residualPhaseFraction 1e-3;
+        residualSlip 1e-3;
     }
 );