diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
index 4150cd7503..5460944288 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index c41bc9b6c7..401e1203d2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
index c41bc9b6c7..401e1203d2 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index a813a2d8b7..bcf99e1668 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index aa476b4dc8..0170d6e384 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
index dcec73b268..68d70ae718 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
@@ -1,7 +1,6 @@
 regionProperties/regionProperties.C
 
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C
 solid/solidRegionDiffNo.C
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
deleted file mode 100644
index 42a73343cc..0000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
+++ /dev/null
@@ -1,381 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 "solidWallMixedTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "directMappedPatchBase.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
-(
-    const polyMesh& nbrRegion
-) const
-{
-    const fvMesh& myRegion = patch().boundaryMesh().mesh();
-
-    const regionProperties& props =
-        myRegion.objectRegistry::parent().lookupObject<regionProperties>
-        (
-            "regionProperties"
-        );
-
-    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
-    if (myIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), myRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn
-            (
-                "solidWallMixedTemperatureCoupledFvPatchScalarField"
-                "::interfaceOwner(const polyMesh&) const"
-            )   << "Cannot find region " << myRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        myIndex = props.fluidRegionNames().size() + i;
-    }
-    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
-    if (nbrIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
-                << "Cannot find region " << nbrRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        nbrIndex = props.fluidRegionNames().size() + i;
-    }
-
-    return myIndex < nbrIndex;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    neighbourFieldName_("undefined-neighbourFieldName"),
-    KName_("undefined-K")
-{
-    this->refValue() = 0.0;
-    this->refGrad() = 0.0;
-    this->valueFraction() = 1.0;
-    this->fixesValue_ = true;
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const solidWallMixedTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mixedFvPatchScalarField(ptf, p, iF, mapper),
-    neighbourFieldName_(ptf.neighbourFieldName_),
-    KName_(ptf.KName_),
-    fixesValue_(ptf.fixesValue_)
-{}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    neighbourFieldName_(dict.lookup("neighbourFieldName")),
-    KName_(dict.lookup("K"))
-{
-    if (!isA<directMappedPatchBase>(this->patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<scalar, volMesh>& iF,\n"
-            "    const dictionary& dict\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
-            << "' not type '" << directMappedPatchBase::typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << dimensionedInternalField().name()
-            << " in file " << dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
-    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-
-    if (dict.found("refValue"))
-    {
-        // Full restart
-        refValue() = scalarField("refValue", dict, p.size());
-        refGrad() = scalarField("refGradient", dict, p.size());
-        valueFraction() = scalarField("valueFraction", dict, p.size());
-        fixesValue_ = readBool(dict.lookup("fixesValue"));
-    }
-    else
-    {
-        // Start from user entered data. Assume fixedValue.
-        refValue() = *this;
-        refGrad() = 0.0;
-        valueFraction() = 1.0;
-        fixesValue_ = true;
-    }
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const solidWallMixedTemperatureCoupledFvPatchScalarField& wtcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(wtcsf, iF),
-    neighbourFieldName_(wtcsf.neighbourFieldName_),
-    KName_(wtcsf.KName_),
-    fixesValue_(wtcsf.fixesValue_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::fvPatchScalarField&
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::K() const
-{
-    return this->patch().lookupPatchField<volScalarField, scalar>(KName_);
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    // Get the coupling information from the directMappedPatchBase
-    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
-    (
-        patch().patch()
-    );
-    const polyMesh& nbrMesh = mpp.sampleMesh();
-    // Force recalculation of mapping and schedule
-    const mapDistribute& distMap = mpp.map();
-    (void)distMap.schedule();
-
-    tmp<scalarField> intFld = patchInternalField();
-
-    if (interfaceOwner(nbrMesh))
-    {
-        // Note: other side information could be cached - it only needs
-        // to be updated the first time round the iteration (i.e. when
-        // switching regions) but unfortunately we don't have this information.
-
-        const fvPatch& nbrPatch = refCast<const fvMesh>
-        (
-            nbrMesh
-        ).boundary()[mpp.samplePolyPatch().index()];
-
-
-        // Calculate the temperature by harmonic averaging
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
-        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
-        (
-            nbrPatch.lookupPatchField<volScalarField, scalar>
-            (
-                neighbourFieldName_
-            )
-        );
-
-        // Swap to obtain full local values of neighbour internal field
-        scalarField nbrIntFld = nbrField.patchInternalField();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrIntFld
-        );
-
-        // Swap to obtain full local values of neighbour K*delta
-        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrKDelta
-        );
-
-
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-
-        // Calculate common wall temperature. Reuse *this to store common value.
-        scalarField Twall
-        (
-            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
-          / (myKDelta() + nbrKDelta)
-        );
-        // Assign to me
-        fvPatchScalarField::operator=(Twall);
-        // Distribute back and assign to neighbour
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            nbrField.size(),
-            distMap.constructMap(),     // reverse : what to send
-            distMap.subMap(),
-            Twall
-        );
-        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
-        (
-            nbrField
-        ).fvPatchScalarField::operator=(Twall);
-    }
-
-
-    // Switch between fixed value (of harmonic avg) or gradient
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    label nFixed = 0;
-
-    // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
-
-    if (debug)
-    {
-        scalar Q = gSum(K()*patch().magSf()*normalGradient());
-
-        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() :"
-            << " patch:" << patch().name()
-            << " heatFlux:" << Q
-            << " walltemperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
-            << endl;
-    }
-
-    forAll(*this, i)
-    {
-        // if outgoing flux use fixed value.
-        if (normalGradient()[i] < 0.0)
-        {
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = 0.0;   // not used by me
-            this->valueFraction()[i] = 1.0;
-            nFixed++;
-        }
-        else
-        {
-            // Fixed gradient. Make sure to have valid refValue (even though
-            // I am not using it - other boundary conditions might)
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = normalGradient()[i];
-            this->valueFraction()[i] = 0.0;
-        }
-    }
-
-    reduce(nFixed, sumOp<label>());
-
-    fixesValue_ = (nFixed > 0);
-
-    if (debug)
-    {
-        label nTotSize = returnReduce(this->size(), sumOp<label>());
-
-        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() :"
-            << " patch:" << patch().name()
-            << " out of:" << nTotSize
-            << " fixedBC:" << nFixed
-            << " gradient:" << nTotSize-nFixed << endl;
-    }
-
-    mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    mixedFvPatchScalarField::write(os);
-    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallMixedTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
deleted file mode 100644
index f4766e5dc3..0000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
+++ /dev/null
@@ -1,192 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    solidWallMixedTemperatureCoupledFvPatchScalarField
-
-Description
-    Mixed boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-    If my temperature is T1, neighbour is T2:
-
-    T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
-    
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallMixedTemperatureCoupled;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-    Needs to be on underlying directMapped(Wall)FvPatch.
-
-    Note: runs in parallel with arbitrary decomposition. Uses directMapped
-    functionality to calculate exchange.
-
-    Note: lags interface data so both sides use same data.
-    - problem: schedule to calculate average would interfere
-    with standard processor swaps.
-    - so: updateCoeffs sets both to same Twall. Only need to do
-    this for last outer iteration but don't have access to this.
-
-SourceFiles
-    solidWallMixedTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallMixedTemperatureCoupledFvPatchScalarField_H
-#define solidWallMixedTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "mixedFvPatchFields.H"
-#include "fvPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-        Class solidWallMixedTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallMixedTemperatureCoupledFvPatchScalarField
-:
-    public mixedFvPatchScalarField
-{
-    // Private data
-
-        //- Name of field on the neighbour region
-        const word neighbourFieldName_;
-        
-        //- Name of thermal conductivity field
-        const word KName_;
-
-        bool fixesValue_;
-
-
-    // Private Member Functions
-
-        //- Am I or neighbour owner of interface
-        bool interfaceOwner(const polyMesh& nbrRegion) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallMixedTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallMixedTemperatureCoupledFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallMixedTemperatureCoupledFvPatchScalarField
-                (
-                    *this,
-                    iF
-                )
-            );
-        }
-
-
-    // Member functions
-
-        //- Get corresponding K field
-        const fvPatchScalarField& K() const;
-
-        //- Return true if this patch field fixes a value.
-        //  Needed to check if a level has to be specified while solving
-        //  Poissons equations.
-        virtual bool fixesValue() const
-        {
-            return fixesValue_;
-        }
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
index 8f1f654d2b..b5db5078f1 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -1,18 +1,16 @@
     scalar CoNum = -GREAT;
-    if (fluidRegions.size())
+
+    forAll(fluidRegions, regionI)
     {
-        forAll(fluidRegions, regionI)
-        {
-            CoNum = max
+        CoNum = max
+        (
+            compressibleCourantNo
             (
-                compressibleCourantNo
-                (
-                    fluidRegions[regionI],
-                    runTime,
-                    rhoFluid[regionI],
-                    phiFluid[regionI]
-                ),
-                CoNum
-            );
-        }
+                fluidRegions[regionI],
+                runTime,
+                rhoFluid[regionI],
+                phiFluid[regionI]
+            ),
+            CoNum
+        );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 20ca225676..c2e5668a68 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -1,18 +1,16 @@
     scalar DiNum = -GREAT;
-    if (solidRegions.size())
+
+    forAll(solidRegions, regionI)
     {
-        forAll(solidRegions, regionI)
-        {
-            DiNum = max
+        DiNum = max
+        (
+            solidRegionDiffNo
             (
-                solidRegionDiffNo
-                (
-                    solidRegions[regionI],
-                    runTime,
-                    rhosCps[regionI],
-                    Ks[regionI]
-                ),
-                DiNum
-            );
-        }
+                solidRegions[regionI],
+                runTime,
+                rhosCps[regionI],
+                Ks[regionI]
+            ),
+            DiNum
+        );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
index bf48d47292..e2b1d536d2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
@@ -37,7 +37,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    fixedGradientFvPatchScalarField(p, iF),
     q_(p.size(), 0.0),
     KName_("undefined-K")
 {}
@@ -52,7 +52,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
     q_(ptf.q_, mapper),
     KName_(ptf.KName_)
 {}
@@ -66,7 +66,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict),
+    fixedGradientFvPatchScalarField(p, iF, dict),
     q_("q", dict, p.size()),
     KName_(dict.lookup("K"))
 {}
@@ -78,7 +78,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
 )
 :
-    fixedValueFvPatchScalarField(tppsf),
+    fixedGradientFvPatchScalarField(tppsf),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -91,7 +91,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(tppsf, iF),
+    fixedGradientFvPatchScalarField(tppsf, iF),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -104,7 +104,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchScalarField::autoMap(m);
+    fixedGradientFvPatchScalarField::autoMap(m);
     q_.autoMap(m);
 }
 
@@ -115,7 +115,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchScalarField::rmap(ptf, addr);
+    fixedGradientFvPatchScalarField::rmap(ptf, addr);
 
     const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
         refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
@@ -131,14 +131,14 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
+    const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
+    (
+        KName_
+    );
 
-    const fvPatchScalarField& Tw = *this;
+    gradient() = q_/Kw;
 
-    operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
-
-    fixedValueFvPatchScalarField::updateCoeffs();
+    fixedGradientFvPatchScalarField::updateCoeffs();
 }
 
 
@@ -147,9 +147,10 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
     Ostream& os
 ) const
 {
-    fixedValueFvPatchScalarField::write(os);
+    fixedGradientFvPatchScalarField::write(os);
     q_.writeEntry("q", os);
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
 }
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
index 7558c4df7e..ffbcf52ed5 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
@@ -45,7 +45,7 @@ SourceFiles
 #ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
 #define solidWallHeatFluxTemperatureFvPatchScalarField_H
 
-#include "fixedValueFvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ namespace Foam
 
 class solidWallHeatFluxTemperatureFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public fixedGradientFvPatchScalarField
 {
     // Private data
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
index 1cbafe368a..7f317cc2e6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
@@ -2,7 +2,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turb.divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
index c1925c7708..688120809a 100644
--- a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevReff(U)
     );
 
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index 88ca441e98..7da3c22269 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -39,17 +39,29 @@ Description
     - mesh with cells put into cellZones (-makeCellZones)
 
     Note:
+    - cellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) cellZone.
+
+    - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
+    from the specified file. This allows one to explicitly specify the region
+    distribution and still have multiple cellZones per region.
+
+    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) cellZone.
+
+
     - Should work in parallel.
     cellZones can differ on either side of processor boundaries in which case
     the faces get moved from processor patch to directMapped patch. Not
-    ery well tested.
+    very well tested.
+
     - If a cell zone gets split into more than one region it can detect
     the largest matching region (-sloppyCellZones). This will accept any
     region that covers more than 50% of the zone. It has to be a subset
     so cannot have any cells in any other zone.
-    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
-    this if you don't mind having disconnected domains in a single region.
-    This option requires all cells to be in one (and one only) cellZone.
+
     - writes maps like decomposePar back to original mesh:
         - pointRegionAddressing : for every point in this region the point in
         the original mesh
@@ -1137,63 +1149,6 @@ EdgeMap<label> addRegionPatches
 }
 
 
-//// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
-//// if no zone found, zone otherwise
-//label findCorrespondingSubZone
-//(
-//    const cellZoneMesh& cellZones,
-//    const labelList& existingZoneID,
-//    const labelList& cellRegion,
-//    const label regionI
-//)
-//{
-//    // Zone corresponding to region. No corresponding zone.
-//    label zoneI = labelMax;
-//
-//    labelList regionCells = findIndices(cellRegion, regionI);
-//
-//    if (regionCells.empty())
-//    {
-//        // My local portion is empty. Maps to any empty cellZone. Mark with
-//        // special value which can get overwritten by other processors.
-//        zoneI = -1;
-//    }
-//    else
-//    {
-//        // Get zone for first element.
-//        zoneI = existingZoneID[regionCells[0]];
-//
-//        if (zoneI == -1)
-//        {
-//            zoneI = labelMax;
-//        }
-//        else
-//        {
-//            // 1. All regionCells in zoneI?
-//            forAll(regionCells, i)
-//            {
-//                if (existingZoneID[regionCells[i]] != zoneI)
-//                {
-//                    zoneI = labelMax;
-//                    break;
-//                }
-//            }
-//        }
-//    }
-//
-//    // Determine same zone over all processors.
-//    reduce(zoneI, maxOp<label>());
-//
-//    if (zoneI == labelMax)
-//    {
-//        // Cells in region that are not in zoneI
-//        zoneI = -1;
-//    }
-//
-//    return zoneI;
-//}
-
-
 // Find region that covers most of cell zone
 label findCorrespondingRegion
 (
@@ -1314,62 +1269,20 @@ label findCorrespondingRegion
 //}
 
 
-// Main program:
-
-int main(int argc, char *argv[])
+// Get zone per cell
+// - non-unique zoning
+// - coupled zones
+void getZoneID
+(
+    const polyMesh& mesh,
+    const cellZoneMesh& cellZones,
+    labelList& zoneID,
+    labelList& neiZoneID
+)
 {
-#   include "addOverwriteOption.H"
-    argList::addBoolOption("cellZones");
-    argList::addBoolOption("cellZonesOnly");
-    argList::addOption("blockedFaces", "faceSet");
-    argList::addBoolOption("makeCellZones");
-    argList::addBoolOption("largestOnly");
-    argList::addOption("insidePoint", "point");
-    argList::addBoolOption("detectOnly");
-    argList::addBoolOption("sloppyCellZones");
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-    runTime.functionObjects().off();
-#   include "createMesh.H"
-    const word oldInstance = mesh.pointsInstance();
-
-    word blockedFacesName;
-    if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
-    {
-        Info<< "Reading blocked internal faces from faceSet "
-            << blockedFacesName << nl << endl;
-    }
-
-    const bool makeCellZones    = args.optionFound("makeCellZones");
-    const bool largestOnly      = args.optionFound("largestOnly");
-    const bool insidePoint      = args.optionFound("insidePoint");
-    const bool useCellZones     = args.optionFound("cellZones");
-    const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
-    const bool overwrite        = args.optionFound("overwrite");
-    const bool detectOnly       = args.optionFound("detectOnly");
-    const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
-
-    if (insidePoint && largestOnly)
-    {
-        FatalErrorIn(args.executable())
-            << "You cannot specify both -largestOnly"
-            << " (keep region with most cells)"
-            << " and -insidePoint (keep region containing point)"
-            << exit(FatalError);
-    }
-
-
-    const cellZoneMesh& cellZones = mesh.cellZones();
-
-
-    // Collect zone per cell
-    // ~~~~~~~~~~~~~~~~~~~~~
-    // - non-unique zoning
-    // - coupled zones
-
     // Existing zoneID
-    labelList zoneID(mesh.nCells(), -1);
+    zoneID.setSize(mesh.nCells());
+    zoneID = -1;
 
     forAll(cellZones, zoneI)
     {
@@ -1384,7 +1297,7 @@ int main(int argc, char *argv[])
             }
             else
             {
-                FatalErrorIn(args.executable())
+                FatalErrorIn("getZoneID(..)")
                     << "Cell " << cellI << " with cell centre "
                     << mesh.cellCentres()[cellI]
                     << " is multiple zones. This is not allowed." << endl
@@ -1396,175 +1309,42 @@ int main(int argc, char *argv[])
     }
 
     // Neighbour zoneID.
-    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
 
     forAll(neiZoneID, i)
     {
         neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
     }
     syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
+}
 
 
-    // Determine connected regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+void matchRegions
+(
+    const bool sloppyCellZones,
+    const polyMesh& mesh,
 
-    // Mark additional faces that are blocked
-    boolList blockedFace;
+    const label nCellRegions,
+    const labelList& cellRegion,
 
-    // Read from faceSet
-    if (blockedFacesName.size())
-    {
-        faceSet blockedFaceSet(mesh, blockedFacesName);
-        Info<< "Read " << returnReduce(blockedFaceSet.size(), sumOp<label>())
-            << " blocked faces from set " << blockedFacesName << nl << endl;
+    labelList& regionToZone,
+    wordList& regionNames,
+    labelList& zoneToRegion
+)
+{
+    const cellZoneMesh& cellZones = mesh.cellZones();
 
-        blockedFace.setSize(mesh.nFaces(), false);
-
-        forAllConstIter(faceSet, blockedFaceSet, iter)
-        {
-            blockedFace[iter.key()] = true;
-        }
-    }
-
-    // Imply from differing cellZones
-    if (useCellZones)
-    {
-        blockedFace.setSize(mesh.nFaces(), false);
-
-        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-        {
-            label own = mesh.faceOwner()[faceI];
-            label nei = mesh.faceNeighbour()[faceI];
-
-            if (zoneID[own] != zoneID[nei])
-            {
-                blockedFace[faceI] = true;
-            }
-        }
-
-        // Different cellZones on either side of processor patch.
-        forAll(neiZoneID, i)
-        {
-            label faceI = i+mesh.nInternalFaces();
-
-            if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
-            {
-                blockedFace[faceI] = true;
-            }
-        }
-    }
-
-
-    // Determine per cell the region it belongs to
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // cellRegion is the labelList with the region per cell.
-    labelList cellRegion;
-    label nCellRegions = 0;
-    if (useCellZonesOnly)
-    {
-        label unzonedCellI = findIndex(zoneID, -1);
-        if (unzonedCellI != -1)
-        {
-            FatalErrorIn(args.executable())
-                << "For the cellZonesOnly option all cells "
-                << "have to be in a cellZone." << endl
-                << "Cell " << unzonedCellI
-                << " at" << mesh.cellCentres()[unzonedCellI]
-                << " is not in a cellZone. There might be more unzoned cells."
-                << exit(FatalError);
-        }
-        cellRegion = zoneID;
-        nCellRegions = gMax(cellRegion)+1;
-    }
-    else
-    {
-        // Do a topological walk to determine regions
-        regionSplit regions(mesh, blockedFace);
-        nCellRegions = regions.nRegions();
-        cellRegion.transfer(regions);
-    }
-
-    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
-
-
-    // Write to manual decomposition option
-    {
-        labelIOList cellToRegion
-        (
-            IOobject
-            (
-                "cellToRegion",
-                mesh.facesInstance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            cellRegion
-        );
-        cellToRegion.write();
-
-        Info<< "Writing region per cell file (for manual decomposition) to "
-            << cellToRegion.objectPath() << nl << endl;
-    }
-    // Write for postprocessing
-    {
-        volScalarField cellToRegion
-        (
-            IOobject
-            (
-                "cellToRegion",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, 0),
-            zeroGradientFvPatchScalarField::typeName
-        );
-        forAll(cellRegion, cellI)
-        {
-            cellToRegion[cellI] = cellRegion[cellI];
-        }
-        cellToRegion.write();
-
-        Info<< "Writing region per cell as volScalarField to "
-            << cellToRegion.objectPath() << nl << endl;
-    }
-
-
-    // Sizes per region
-    // ~~~~~~~~~~~~~~~~
-
-    labelList regionSizes(nCellRegions, 0);
-
-    forAll(cellRegion, cellI)
-    {
-        regionSizes[cellRegion[cellI]]++;
-    }
-    forAll(regionSizes, regionI)
-    {
-        reduce(regionSizes[regionI], sumOp<label>());
-    }
-
-    Info<< "Region\tCells" << nl
-        << "------\t-----" << endl;
-
-    forAll(regionSizes, regionI)
-    {
-        Info<< regionI << '\t' << regionSizes[regionI] << nl;
-    }
-    Info<< endl;
+    regionToZone.setSize(nCellRegions, -1);
+    regionNames.setSize(nCellRegions);
+    zoneToRegion.setSize(cellZones.size(), -1);
 
+    // Get current per cell zoneID
+    labelList zoneID(mesh.nCells(), -1);
+    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    getZoneID(mesh, cellZones, zoneID, neiZoneID);
 
     // Sizes per cellzone
-    // ~~~~~~~~~~~~~~~~~~
-
     labelList zoneSizes(cellZones.size(), 0);
-    if (useCellZones || makeCellZones || sloppyCellZones)
     {
         List<wordList> zoneNames(Pstream::nProcs());
         zoneNames[Pstream::myProcNo()] = cellZones.names();
@@ -1575,7 +1355,7 @@ int main(int argc, char *argv[])
         {
             if (zoneNames[procI] != zoneNames[0])
             {
-                FatalErrorIn(args.executable())
+                FatalErrorIn("matchRegions(..)")
                     << "cellZones not synchronised across processors." << endl
                     << "Master has cellZones " << zoneNames[0] << endl
                     << "Processor " << procI
@@ -1595,16 +1375,6 @@ int main(int argc, char *argv[])
     }
 
 
-    // Whether region corresponds to a cellzone
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Region per zone
-    labelList regionToZone(nCellRegions, -1);
-    // Name of region
-    wordList regionNames(nCellRegions);
-    // Zone to region
-    labelList zoneToRegion(cellZones.size(), -1);
-
     if (sloppyCellZones)
     {
         Info<< "Trying to match regions to existing cell zones;"
@@ -1624,7 +1394,7 @@ int main(int argc, char *argv[])
             if (regionI != -1)
             {
                 Info<< "Sloppily matched region " << regionI
-                    << " size " << regionSizes[regionI]
+                    //<< " size " << regionSizes[regionI]
                     << " to zone " << zoneI << " size " << zoneSizes[zoneI]
                     << endl;
                 zoneToRegion[zoneI] = regionI;
@@ -1664,6 +1434,330 @@ int main(int argc, char *argv[])
             regionNames[regionI] = "domain" + Foam::name(regionI);
         }
     }
+}
+
+
+void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
+{
+    // Write to manual decomposition option
+    {
+        labelIOList cellToRegion
+        (
+            IOobject
+            (
+                "cellToRegion",
+                mesh.facesInstance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            cellRegion
+        );
+        cellToRegion.write();
+
+        Info<< "Writing region per cell file (for manual decomposition) to "
+            << cellToRegion.objectPath() << nl << endl;
+    }
+    // Write for postprocessing
+    {
+        volScalarField cellToRegion
+        (
+            IOobject
+            (
+                "cellToRegion",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        forAll(cellRegion, cellI)
+        {
+            cellToRegion[cellI] = cellRegion[cellI];
+        }
+        cellToRegion.write();
+
+        Info<< "Writing region per cell as volScalarField to "
+            << cellToRegion.objectPath() << nl << endl;
+    }
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addOverwriteOption.H"
+    argList::addBoolOption("cellZones");
+    argList::addBoolOption("cellZonesOnly");
+    argList::addOption("cellZonesFileOnly", "cellZonesName");
+    argList::addOption("blockedFaces", "faceSet");
+    argList::addBoolOption("makeCellZones");
+    argList::addBoolOption("largestOnly");
+    argList::addOption("insidePoint", "point");
+    argList::addBoolOption("detectOnly");
+    argList::addBoolOption("sloppyCellZones");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+    runTime.functionObjects().off();
+#   include "createMesh.H"
+    const word oldInstance = mesh.pointsInstance();
+
+    word blockedFacesName;
+    if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
+    {
+        Info<< "Reading blocked internal faces from faceSet "
+            << blockedFacesName << nl << endl;
+    }
+
+    const bool makeCellZones    = args.optionFound("makeCellZones");
+    const bool largestOnly      = args.optionFound("largestOnly");
+    const bool insidePoint      = args.optionFound("insidePoint");
+    const bool useCellZones     = args.optionFound("cellZones");
+    const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
+    const bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
+    const bool overwrite        = args.optionFound("overwrite");
+    const bool detectOnly       = args.optionFound("detectOnly");
+    const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
+    if
+    (
+        (useCellZonesOnly || useCellZonesFile)
+     && (
+            blockedFacesName != word::null
+         || useCellZones
+        )
+    )
+    {
+        FatalErrorIn(args.executable())
+            << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
+            << " (which specify complete split)"
+            << " in combination with -blockedFaces or -cellZones"
+            << " (which imply a split based on topology)"
+            << exit(FatalError);
+    }
+
+
+
+    if (insidePoint && largestOnly)
+    {
+        FatalErrorIn(args.executable())
+            << "You cannot specify both -largestOnly"
+            << " (keep region with most cells)"
+            << " and -insidePoint (keep region containing point)"
+            << exit(FatalError);
+    }
+
+
+    const cellZoneMesh& cellZones = mesh.cellZones();
+
+    // Existing zoneID
+    labelList zoneID(mesh.nCells(), -1);
+    // Neighbour zoneID.
+    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    getZoneID(mesh, cellZones, zoneID, neiZoneID);
+
+
+
+    // Determine per cell the region it belongs to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // cellRegion is the labelList with the region per cell.
+    labelList cellRegion;
+    // Region per zone
+    labelList regionToZone;
+    // Name of region
+    wordList regionNames;
+    // Zone to region
+    labelList zoneToRegion;
+
+    label nCellRegions = 0;
+    if (useCellZonesOnly)
+    {
+        Info<< "Using current cellZones to split mesh into regions."
+            << " This requires all"
+            << " cells to be in one and only one cellZone." << nl << endl;
+
+        label unzonedCellI = findIndex(zoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = zoneID;
+        nCellRegions = gMax(cellRegion)+1;
+        regionToZone.setSize(nCellRegions);
+        regionNames.setSize(nCellRegions);
+        zoneToRegion.setSize(cellZones.size(), -1);
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
+        {
+            regionToZone[regionI] = regionI;
+            zoneToRegion[regionI] = regionI;
+            regionNames[regionI] = cellZones[regionI].name();
+        }
+    }
+    else if (useCellZonesFile)
+    {
+        const word zoneFile = args.option("cellZonesFileOnly");
+        Info<< "Reading split from cellZones file " << zoneFile << endl
+            << "This requires all"
+            << " cells to be in one and only one cellZone." << nl << endl;
+
+        cellZoneMesh newCellZones
+        (
+            IOobject
+            (
+                zoneFile,
+                mesh.facesInstance(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh
+        );
+
+        labelList newZoneID(mesh.nCells(), -1);
+        labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+        getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
+
+        label unzonedCellI = findIndex(newZoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesFileOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = newZoneID;
+        nCellRegions = gMax(cellRegion)+1;
+        zoneToRegion.setSize(newCellZones.size(), -1);
+        regionToZone.setSize(nCellRegions);
+        regionNames.setSize(nCellRegions);
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
+        {
+            regionToZone[regionI] = regionI;
+            zoneToRegion[regionI] = regionI;
+            regionNames[regionI] = newCellZones[regionI].name();
+        }
+    }
+    else
+    {
+        // Determine connected regions
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Mark additional faces that are blocked
+        boolList blockedFace;
+
+        // Read from faceSet
+        if (blockedFacesName.size())
+        {
+            faceSet blockedFaceSet(mesh, blockedFacesName);
+            Info<< "Read "
+                << returnReduce(blockedFaceSet.size(), sumOp<label>())
+                << " blocked faces from set " << blockedFacesName << nl << endl;
+
+            blockedFace.setSize(mesh.nFaces(), false);
+
+            forAllConstIter(faceSet, blockedFaceSet, iter)
+            {
+                blockedFace[iter.key()] = true;
+            }
+        }
+
+        // Imply from differing cellZones
+        if (useCellZones)
+        {
+            blockedFace.setSize(mesh.nFaces(), false);
+
+            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+            {
+                label own = mesh.faceOwner()[faceI];
+                label nei = mesh.faceNeighbour()[faceI];
+
+                if (zoneID[own] != zoneID[nei])
+                {
+                    blockedFace[faceI] = true;
+                }
+            }
+
+            // Different cellZones on either side of processor patch.
+            forAll(neiZoneID, i)
+            {
+                label faceI = i+mesh.nInternalFaces();
+
+                if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
+                {
+                    blockedFace[faceI] = true;
+                }
+            }
+        }
+
+        // Do a topological walk to determine regions
+        regionSplit regions(mesh, blockedFace);
+        nCellRegions = regions.nRegions();
+        cellRegion.transfer(regions);
+
+        // Make up region names. If possible match them to existing zones.
+        matchRegions
+        (
+            sloppyCellZones,
+            mesh,
+            nCellRegions,
+            cellRegion,
+
+            regionToZone,
+            regionNames,
+            zoneToRegion
+        );
+    }
+
+    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
+
+
+    // Write decomposition to file
+    writeCellToRegion(mesh, cellRegion);
+
+
+
+    // Sizes per region
+    // ~~~~~~~~~~~~~~~~
+
+    labelList regionSizes(nCellRegions, 0);
+
+    forAll(cellRegion, cellI)
+    {
+        regionSizes[cellRegion[cellI]]++;
+    }
+    forAll(regionSizes, regionI)
+    {
+        reduce(regionSizes[regionI], sumOp<label>());
+    }
+
+    Info<< "Region\tCells" << nl
+        << "------\t-----" << endl;
+
+    forAll(regionSizes, regionI)
+    {
+        Info<< regionI << '\t' << regionSizes[regionI] << nl;
+    }
+    Info<< endl;
+
 
 
     // Print region to zone
@@ -1676,16 +1770,6 @@ int main(int argc, char *argv[])
     }
     Info<< endl;
 
-    //// Print zone to region
-    //Info<< "Zone\tName\tRegion" << nl
-    //    << "----\t----\t------" << endl;
-    //forAll(zoneToRegion, zoneI)
-    //{
-    //    Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
-    //        << zoneToRegion[zoneI] << nl;
-    //}
-    //Info<< endl;
-
 
 
     // Since we're going to mess with patches make sure all non-processor ones
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index ab7193a682..0d47c7db9d 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -70,6 +70,7 @@ namespace Foam
 
 class mapPolyMesh;
 class globalIndex;
+class PstreamBuffers;
 
 /*---------------------------------------------------------------------------*\
                            Class mapDistribute Declaration
@@ -318,6 +319,13 @@ public:
                 }
             }
 
+            //- Do all sends using PstreamBuffers
+            template<class T>
+            void send(PstreamBuffers&, const List<T>&) const;
+            //- Do all receives using PstreamBuffers
+            template<class T>
+            void receive(PstreamBuffers&, List<T>&) const;
+
             //- Correct for topo change.
             void updateMesh(const mapPolyMesh&)
             {
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index e219bc95ab..8797519d4d 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -185,7 +185,7 @@ void Foam::mapDistribute::distribute
     {
         if (!contiguous<T>())
         {
-            PstreamBuffers pBuffs(Pstream::nonBlocking);
+            PstreamBuffers pBufs(Pstream::nonBlocking);
 
             // Stream data into buffer
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -195,13 +195,13 @@ void Foam::mapDistribute::distribute
                 if (domain != Pstream::myProcNo() && map.size())
                 {
                     // Put data into send buffer
-                    UOPstream toDomain(domain, pBuffs);
+                    UOPstream toDomain(domain, pBufs);
                     toDomain << UIndirectList<T>(field, map);
                 }
             }
 
             // Start receiving
-            pBuffs.finishedSends();
+            pBufs.finishedSends();
 
             {
                 // Set up 'send' to myself
@@ -231,7 +231,7 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    UIPstream str(domain, pBuffs);
+                    UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
                     if (recvField.size() != map.size())
@@ -551,7 +551,7 @@ void Foam::mapDistribute::distribute
     {
         if (!contiguous<T>())
         {
-            PstreamBuffers pBuffs(Pstream::nonBlocking);
+            PstreamBuffers pBufs(Pstream::nonBlocking);
 
             // Stream data into buffer
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -561,13 +561,13 @@ void Foam::mapDistribute::distribute
                 if (domain != Pstream::myProcNo() && map.size())
                 {
                     // Put data into send buffer
-                    UOPstream toDomain(domain, pBuffs);
+                    UOPstream toDomain(domain, pBufs);
                     toDomain << UIndirectList<T>(field, map);
                 }
             }
 
             // Start receiving
-            pBuffs.finishedSends();
+            pBufs.finishedSends();
 
             {
                 // Set up 'send' to myself
@@ -597,7 +597,7 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    UIPstream str(domain, pBuffs);
+                    UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
                     if (recvField.size() != map.size())
@@ -757,4 +757,66 @@ void Foam::mapDistribute::distribute
 }
 
 
+template<class T>
+void Foam::mapDistribute::send(PstreamBuffers& pBufs, const List<T>& field)
+const
+{
+    // Stream data into buffer
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = subMap_[domain];
+
+        if (map.size())
+        {
+            // Put data into send buffer
+            UOPstream toDomain(domain, pBufs);
+            toDomain << UIndirectList<T>(field, map);
+        }
+    }
+
+    // Start sending and receiving but do not block.
+    pBufs.finishedSends(false);
+}
+
+
+template<class T>
+void Foam::mapDistribute::receive(PstreamBuffers& pBufs, List<T>& field) const
+{
+    // Consume
+    field.setSize(constructSize_);
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        if (map.size())
+        {
+            UIPstream str(domain, pBufs);
+            List<T> recvField(str);
+
+            if (recvField.size() != map.size())
+            {
+                FatalErrorIn
+                (
+                    "template<class T>\n"
+                    "void mapDistribute::receive\n"
+                    "(\n"
+                    "    PstreamBuffers&,\n"
+                    "    List<T>&\n"
+                    ")\n"
+                )   << "Expected from processor " << domain
+                    << " " << map.size() << " but received "
+                    << recvField.size() << " elements."
+                    << abort(FatalError);
+            }
+
+            forAll(map, i)
+            {
+                field[map[i]] = recvField[i];
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 7dcfceaf10..f28e9affd0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -2189,8 +2189,8 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         if (debug || report)
         {
-            Info<< " ***Concave cells found, number of cells: "
-                << nConcaveCells << endl;
+            Info<< " ***Concave cells (using face planes) found,"
+                << " number of cells: " << nConcaveCells << endl;
         }
 
         return true;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index e2a2fee10a..0d9cecd039 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1234,7 +1234,7 @@ void Foam::meshRefinement::findCellZoneTopo
         {
             label surfI = namedSurfaceIndex[faceI];
 
-            if (surfI != -1 && surfaceToCellZone[surfI] != -1)
+            if (surfI != -1)
             {
                 // Calculate region to zone from cellRegions on either side
                 // of internal face.
@@ -1286,7 +1286,7 @@ void Foam::meshRefinement::findCellZoneTopo
 
                     label surfI = namedSurfaceIndex[faceI];
 
-                    if (surfI != -1 && surfaceToCellZone[surfI] != -1)
+                    if (surfI != -1)
                     {
                         bool changedCell = calcRegionToZone
                         (
@@ -2327,9 +2327,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Set using walking
     // ~~~~~~~~~~~~~~~~~
 
-    if (!allowFreeStandingZoneFaces)
+    //if (!allowFreeStandingZoneFaces)
     {
-        Info<< "Walking to assign cellZones." << nl << endl;
+        Info<< "Walking from location-in-mesh " << keepPoint
+            << " to assign cellZones "
+            << "- crossing a faceZone face changes cellZone" << nl << endl;
 
         // Topological walk
         findCellZoneTopo
@@ -2520,6 +2522,32 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         mesh_.setInstance(oldInstance());
     }
 
+    // Print some stats (note: zones are synchronised)
+    if (mesh_.cellZones().size() > 0)
+    {
+        Info<< "CellZones:" << endl;
+        forAll(mesh_.cellZones(), zoneI)
+        {
+            const cellZone& cz = mesh_.cellZones()[zoneI];
+            Info<< "    " << cz.name()
+                << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
+                << endl;
+        }
+        Info<< endl;
+    }
+    if (mesh_.faceZones().size() > 0)
+    {
+        Info<< "FaceZones:" << endl;
+        forAll(mesh_.faceZones(), zoneI)
+        {
+            const faceZone& fz = mesh_.faceZones()[zoneI];
+            Info<< "    " << fz.name()
+                << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
+                << endl;
+        }
+        Info<< endl;
+    }
+
     // None of the faces has changed, only the zones. Still...
     updateMesh(map, labelList());
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index c6aae8ccae..ac01e052cd 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -77,11 +77,22 @@ Foam::refinementSurfaces::refinementSurfaces
         if (dict.found("faceZone"))
         {
             dict.lookup("faceZone") >> faceZoneNames_[surfI];
+            bool hasSide = dict.readIfPresent("zoneInside", zoneInside_[surfI]);
             if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
             {
-                dict.lookup("zoneInside") >> zoneInside_[surfI];
+                if (hasSide && !allGeometry_[surfaces_[surfI]].hasVolumeType())
+                {
+                    IOWarningIn
+                    (
+                        "refinementSurfaces::refinementSurfaces(..)",
+                        dict
+                    )   << "Unused entry zoneInside for faceZone "
+                        << faceZoneNames_[surfI]
+                        << " since surface " << names_[surfI]
+                        << " is not closed." << endl;
+                }
             }
-            else if (dict.found("zoneInside"))
+            else if (hasSide)
             {
                 IOWarningIn("refinementSurfaces::refinementSurfaces(..)", dict)
                     << "Unused entry zoneInside for faceZone "
@@ -324,11 +335,30 @@ Foam::refinementSurfaces::refinementSurfaces
             if (dict.found("faceZone"))
             {
                 dict.lookup("faceZone") >> faceZoneNames_[surfI];
+                bool hasSide = dict.readIfPresent
+                (
+                    "zoneInside",
+                    zoneInside_[surfI]
+                );
                 if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
                 {
-                    dict.lookup("zoneInside") >> zoneInside_[surfI];
+                    if
+                    (
+                        hasSide
+                    && !allGeometry_[surfaces_[surfI]].hasVolumeType()
+                    )
+                    {
+                        IOWarningIn
+                        (
+                            "refinementSurfaces::refinementSurfaces(..)",
+                            dict
+                        )   << "Unused entry zoneInside for faceZone "
+                            << faceZoneNames_[surfI]
+                            << " since surface " << names_[surfI]
+                            << " is not closed." << endl;
+                    }
                 }
-                else if (dict.found("zoneInside"))
+                else if (hasSide)
                 {
                     IOWarningIn
                     (
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index c727a74d65..0941f03528 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -64,7 +64,6 @@ int main(int argc, char *argv[])
             tmp<fvVectorMatrix> UEqn
             (
                 fvm::div(phi, U)
-              - fvm::Sp(fvc::div(phi), U)
               + turbulence->divDevReff(U)
             );
             mrfZones.addCoriolis(UEqn());