diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..e924e9de680bcbe388c230a109aa2d01c14b9adf
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
@@ -0,0 +1,13 @@
+regionProperties/regionProperties.C
+
+coupleManager/coupleManager.C
+
+derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
+derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
+
+
+chtMultiRegionFoam.C
+
+EXE = $(FOAM_APPBIN)/chtMultiRegionFoam
+
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..6bf54dad003cf91b4c5e8668cc5494f803a01d10
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@@ -0,0 +1,16 @@
+EXE_INC = \
+    -Ifluid \
+    -Isolid \
+    -IregionProperties \
+    -IcoupleManager \
+    -IderivedFvPatchFields/solidWallTemperatureCoupled \
+    -IderivedFvPatchFields/solidWallHeatFluxTemperatureCoupled \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModels
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..b8151196305023ca58530101a7ca4dd4c0469fd4
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    chtMultiRegionFoam
+
+Description
+    Combination of heatConductionFoam and buoyantFoam for conjugate heat
+    transfer between a solid region and fluid region
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicThermo.H"
+#include "compressible/turbulenceModel/turbulenceModel.H"
+#include "fixedGradientFvPatchFields.H"
+#include "regionProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "solveContinuityEquation.C"
+#include "solveMomentumEquation.C"
+#include "compressibleContinuityErrors.C"
+#include "solvePressureDifferenceEquation.C"
+#include "solveEnthalpyEquation.C"
+#include "compressibleCourantNo.C"
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    regionProperties rp(runTime);
+
+#   include "createFluidMeshes.H"
+#   include "createSolidMeshes.H"
+
+#   include "createFluidFields.H"
+
+#   include "createSolidFields.H"
+
+#   include "initContinuityErrs.H"
+
+#   include "readTimeControls.H"
+
+    if (fluidRegions.size())
+    {
+#       include "compressibleMultiRegionCourantNo.H"
+#       include "setInitialDeltaT.H"
+    }
+
+    while(runTime.run())
+    {
+#       include "readTimeControls.H"
+
+        if (fluidRegions.size())
+        {
+#           include "compressibleMultiRegionCourantNo.H"
+#           include "setDeltaT.H"
+        }
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        forAll(fluidRegions, i)
+        {
+            Info<< "\nSolving for fluid region "
+                << fluidRegions[i].name() << endl;
+#           include "readFluidMultiRegionPISOControls.H"
+#           include "solveFluid.H"
+        }
+
+        forAll(solidRegions, i)
+        {
+            Info<< "\nSolving for solid region "
+                << solidRegions[i].name() << endl;
+#           include "readSolidMultiRegionPISOControls.H"
+#           include "solveSolid.H"
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info << "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C
new file mode 100644
index 0000000000000000000000000000000000000000..0fbe3ea29fdda768310dad8dcf784f07c0504f2f
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.C
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "coupleManager.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::coupleManager::coupleManager(const fvPatch& patch)
+:
+    patch_(patch),
+    neighbourRegionName_("undefined-neighbourRegionName"),
+    neighbourPatchName_("undefined-neighbourPatchName"),
+    neighbourFieldName_("undefined-neighbourFieldName"),
+    localRegion_(patch_.boundaryMesh().mesh())
+{}
+
+
+Foam::coupleManager::coupleManager
+(
+    const fvPatch& patch,
+    const dictionary& dict
+)
+:
+    patch_(patch),
+    neighbourRegionName_(dict.lookup("neighbourRegionName")),
+    neighbourPatchName_(dict.lookup("neighbourPatchName")),
+    neighbourFieldName_(dict.lookup("neighbourFieldName")),
+    localRegion_(patch_.boundaryMesh().mesh())
+{}
+
+
+Foam::coupleManager::coupleManager
+(
+    const coupleManager& cm
+)
+:
+    patch_(cm.patch()),
+    neighbourRegionName_(cm.neighbourRegionName()),
+    neighbourPatchName_(cm.neighbourPatchName()),
+    neighbourFieldName_(cm.neighbourFieldName()),
+    localRegion_(patch_.boundaryMesh().mesh())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::coupleManager::~coupleManager()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::coupleManager::checkCouple() const
+{
+    Info<< "neighbourRegionName_ = " << neighbourRegionName_ << endl;
+    Info<< "neighbourPatchName_ = " << neighbourPatchName_ << endl;
+    Info<< "neighbourFieldName_ = " << neighbourFieldName_ << endl;
+
+    const fvPatch& nPatch = neighbourPatch();
+
+    if (patch_.size() != nPatch.size())
+    {
+        FatalErrorIn("Foam::coupleManager::checkCouple()")
+            << "Unequal patch sizes:" << nl
+            << "    patch name (size) = " << patch_.name()
+            << "(" << patch_.size() << ")" << nl
+            << "    neighbour patch name (size) = "
+            << nPatch.name() << "(" << patch_.size() << ")" << nl
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::coupleManager::coupleToObj() const
+{
+    const fvPatch& nPatch = neighbourPatch();
+
+    OFstream obj
+        (
+             patch_.name() + "_to_" + nPatch.name() + "_couple.obj"
+        );
+    const vectorField& c1 = patch_.Cf();
+    const vectorField& c2 = neighbourPatch().Cf();
+
+    if (c1.size() != c2.size())
+    {
+        FatalErrorIn("coupleManager::coupleToObj() const")
+            << "Coupled patches are of unequal size:" << nl
+            << "    patch0 = " << patch_.name()
+            << "(" << patch_.size() <<  ")" << nl
+            << "    patch1 = " << nPatch.name()
+            << "(" << nPatch.size() <<  ")" << nl
+            << abort(FatalError);
+    }
+
+    forAll(c1, i)
+    {
+        obj << "v " << c1[i].x() << " " << c1[i].y() << " " << c1[i].z() << nl
+            << "v " << c2[i].x() << " " << c2[i].y() << " " << c2[i].z() << nl
+            << "l " << (2*i + 1) << " " << (2*i + 2) << endl;
+    }
+}
+
+
+void Foam::coupleManager::writeEntries(Ostream& os) const
+{
+    os.writeKeyword("neighbourRegionName");
+    os << neighbourRegionName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("neighbourPatchName");
+    os << neighbourPatchName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("neighbourFieldName");
+    os << neighbourFieldName_ << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H
new file mode 100644
index 0000000000000000000000000000000000000000..a6f265faae1558f724dcb3d8a9100ca75bea30d7
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManager.H
@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    coupleManager
+
+Description
+    Object to handle the coupling of region patches. It can be queried to
+    return the neighbour information.
+
+SourceFiles
+    coupleManager.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coupleManager_H
+#define coupleManager_H
+
+#include "Ostream.H"
+#include "dictionary.H"
+#include "fvPatch.H"
+#include "fvMesh.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class coupleManager Declaration
+\*---------------------------------------------------------------------------*/
+
+class coupleManager
+{
+    // Private data
+
+        //- Reference to the local fvPatch
+        const fvPatch& patch_;
+
+        //- Name of neighbour region
+        word neighbourRegionName_;
+
+        //- Name of patch on the neighbour region
+        word neighbourPatchName_;
+
+        //- Name of field on the neighbour region
+        word neighbourFieldName_;
+
+        //- Reference to the local region
+        const fvMesh& localRegion_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+//        coupleManager(const coupleManager&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const coupleManager&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvPatch
+        coupleManager(const fvPatch& patch);
+
+        //- Construct from fvPatch and dictionary
+        coupleManager(const fvPatch& patch, const dictionary& dict);
+
+        //- Copy constructor
+        coupleManager(const coupleManager& cm);
+
+
+    // Destructor
+
+        ~coupleManager();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the local patch
+            inline const fvPatch& patch() const;
+
+            //- Return the name of the neighbour region
+            inline const word& neighbourRegionName() const;
+
+            //- Return the name of the patch on the neighbour region
+            inline const word& neighbourPatchName() const;
+
+            //- Return the name of the field on the neighbour region
+            inline const word& neighbourFieldName() const;
+
+            //- Return a reference to the neighbour mesh
+            inline const fvMesh& neighbourRegion() const;
+
+            //- Return the neighbour patch ID
+            inline const label neighbourPatchID() const;
+
+            //- Return a reference to the neighbour patch
+            inline const fvPatch& neighbourPatch() const;
+
+            //- Return a reference to the neighbour patch field
+            template<class Type>
+            inline const fvPatchField<Type>& neighbourPatchField() const;
+
+            //- Check that the couple is valid
+            void checkCouple() const;
+
+
+        // Edit
+
+            //- Return the name of the neighbour region
+            word& neighbourRegionName();
+
+            //- Return the name of the patch on the neighbour region
+            word& neighbourPatchName();
+
+            //- Return the name of the field on the neighbour region
+            word& neighbourFieldName();
+
+
+        // Write
+
+           //- Write couple to obj file
+           void coupleToObj() const;
+
+           //- Write entries for patches
+           void writeEntries(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "coupleManagerI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H
new file mode 100644
index 0000000000000000000000000000000000000000..5a252d870422421bab3c95e06e79762c0aa6a24c
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/coupleManager/coupleManagerI.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::fvPatch& Foam::coupleManager::patch() const
+{
+    return patch_;
+};
+
+
+inline const Foam::word& Foam::coupleManager::neighbourRegionName() const
+{
+    return neighbourRegionName_;
+};
+
+
+inline const Foam::word& Foam::coupleManager::neighbourPatchName() const
+{
+    return neighbourPatchName_;
+};
+
+
+inline const Foam::word& Foam::coupleManager::neighbourFieldName() const
+{
+     return neighbourFieldName_;
+};
+
+
+inline const Foam::fvMesh& Foam::coupleManager::neighbourRegion() const
+{
+    return localRegion_.objectRegistry::parent()
+        .lookupObject<fvMesh>(neighbourRegionName_);
+}
+
+
+inline const Foam::label Foam::coupleManager::neighbourPatchID() const
+{
+    return neighbourRegion().boundaryMesh().findPatchID(neighbourPatchName_);
+}
+
+
+inline const Foam::fvPatch& Foam::coupleManager::neighbourPatch() const
+{
+    return neighbourRegion().boundary()[neighbourPatchID()];
+}
+
+
+template<class Type>
+inline const Foam::fvPatchField<Type>&
+Foam::coupleManager::neighbourPatchField() const
+{
+    return neighbourPatch().lookupPatchField
+        <GeometricField<Type, fvPatchField, volMesh>, Type>
+            (neighbourFieldName_);
+}
+
+
+inline Foam::word& Foam::coupleManager::neighbourRegionName()
+{
+    return neighbourRegionName_;
+};
+
+
+inline Foam::word& Foam::coupleManager::neighbourPatchName()
+{
+    return neighbourPatchName_;
+};
+
+
+inline Foam::word& Foam::coupleManager::neighbourFieldName()
+{
+    return neighbourFieldName_;
+};
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..c1fdf3f5380db24906d8845a6070754e4d1eef94
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallHeatFluxTemperatureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    q_(p.size(), 0.0),
+    KName_("undefined-K")
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    q_(ptf.q_, mapper),
+    KName_(ptf.KName_)
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    q_("q", dict, p.size()),
+    KName_(dict.lookup("K"))
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf),
+    q_(tppsf.q_),
+    KName_(tppsf.KName_)
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF),
+    q_(tppsf.q_),
+    KName_(tppsf.KName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    q_.autoMap(m);
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
+        refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
+
+    q_.rmap(hfptf.q_, addr);
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& Kw =
+        patch().lookupPatchField<volScalarField, scalar>(KName_);
+
+    const fvPatchScalarField& Tw = *this;
+
+    operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fixedValueFvPatchScalarField::write(os);
+    q_.writeEntry("q", os);
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        solidWallHeatFluxTemperatureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..a03fe7bf39c40a1f78c8821d4a3efa9254e70956
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    solidWallHeatFluxTemperatureFvPatchScalarField
+
+Description
+    Heat flux boundary condition for temperature on solid region
+
+    Example usage:
+        myWallPatch
+        {
+            type            solidWallHeatFluxTemperature;
+            K               K;                 // Name of K field
+            q               uniform 1000;      // Heat flux / [W/m2]
+            value           300.0;             // Initial temperature / [K]
+        }
+
+
+SourceFiles
+    solidWallHeatFluxTemperatureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
+#define solidWallHeatFluxTemperatureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+      Class solidWallHeatFluxTemperatureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallHeatFluxTemperatureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Heat flux / [W/m2]
+        scalarField q_;
+
+        //- Name of thermal conductivity field
+        word KName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallHeatFluxTemperature");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // solidWallHeatFluxTemperatureFvPatchScalarField
+        // onto a new patch
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallHeatFluxTemperatureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&,
+            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 solidWallHeatFluxTemperatureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..43a859a408c446cc164d187a7cce6bee1f0e46ad
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+#include "solidWallTemperatureCoupledFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
+solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    coupleManager_(p),
+    KName_("undefined-K")
+{}
+
+
+Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
+solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
+    coupleManager_(ptf.coupleManager_),
+    KName_(ptf.KName_)
+{}
+
+
+Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
+solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    coupleManager_(p, dict),
+    KName_(dict.lookup("K"))
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        evaluate();
+    }
+}
+
+
+Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::
+solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureCoupledFvPatchScalarField& whftcsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(whftcsf, iF),
+    coupleManager_(whftcsf.coupleManager_),
+    KName_(whftcsf.KName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvPatchField<scalar>& neighbourField =
+        coupleManager_.neighbourPatchField<scalar>();
+
+    const fvPatchField<scalar>& K =
+        patch().lookupPatchField<volScalarField, scalar>(KName_);
+
+    gradient() = refCast<const solidWallTemperatureCoupledFvPatchScalarField>
+        (neighbourField).flux()/K;
+
+    fixedGradientFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallHeatFluxTemperatureCoupledFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    coupleManager_.writeEntries(os);
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..12b75c180020351050e80010b283146769b524d9
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemperatureCoupledFvPatchScalarField.H
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+
+Description
+    Fixed heat-flux boundary condition for temperature, to be used by the
+    conjugate heat transfer solver.
+
+    Example usage:
+        myInterfacePatchName
+        {
+            type                solidWallHeatFluxTemperatureCoupled;
+            neighbourRegionName fluid;
+            neighbourPatchName  fluidSolidInterface;
+            neighbourFieldName  T;
+            K                   K;
+            value               uniform 300;
+        }
+
+SourceFiles
+    solidWallHeatFluxTemperatureCoupledFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
+#define solidWallHeatFluxTemperatureCoupledFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
+#include "coupleManager.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+   Class solidWallHeatFluxTemperatureCoupledFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+:
+    public fixedGradientFvPatchScalarField
+{
+    // Private data
+
+        //- Couple manager object
+        coupleManager coupleManager_;
+
+        //- Name of thermal conductivity field
+        word KName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallHeatFluxTemperatureCoupled");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+        //  onto a new patch
+        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+                    (
+                        *this
+                    )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureCoupledFvPatchScalarField&,
+            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 solidWallHeatFluxTemperatureCoupledFvPatchScalarField
+                    (
+                        *this,
+                        iF
+                    )
+            );
+        }
+
+
+    // Member functions
+
+        //- 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/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..b1480446bd3d7b6c42689adc7dc7bb4d801475a0
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallTemperatureCoupledFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallTemperatureCoupledFvPatchScalarField::
+solidWallTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    coupleManager_(p),
+    KName_("undefined-K")
+{}
+
+
+Foam::solidWallTemperatureCoupledFvPatchScalarField::
+solidWallTemperatureCoupledFvPatchScalarField
+(
+    const solidWallTemperatureCoupledFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    coupleManager_(ptf.coupleManager_),
+    KName_(ptf.KName_)
+{}
+
+
+Foam::solidWallTemperatureCoupledFvPatchScalarField::
+solidWallTemperatureCoupledFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    coupleManager_(p, dict),
+    KName_(dict.lookup("K"))
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        evaluate();
+    }
+}
+
+
+Foam::solidWallTemperatureCoupledFvPatchScalarField::
+solidWallTemperatureCoupledFvPatchScalarField
+(
+    const solidWallTemperatureCoupledFvPatchScalarField& wtcsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(wtcsf, iF),
+    coupleManager_(wtcsf.coupleManager_),
+    KName_(wtcsf.KName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::solidWallTemperatureCoupledFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fvPatchField<scalar>& neighbourField =
+        coupleManager_.neighbourPatchField<scalar>();
+
+    operator==(neighbourField);
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallTemperatureCoupledFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    coupleManager_.writeEntries(os);
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::solidWallTemperatureCoupledFvPatchScalarField::flux() const
+{
+    const fvPatchScalarField& Kw =
+        patch().lookupPatchField<volScalarField, scalar>(KName_);
+
+    const fvPatchScalarField& Tw = *this;
+
+    return Tw.snGrad()*patch().magSf()*Kw;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    solidWallTemperatureCoupledFvPatchScalarField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..abf33a7b433a9b99346f3d78e4b35226a1923120
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.H
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    solidWallHeatFluxCoupledFvPatchScalarField
+
+Description
+    Fixed value boundary condition for temperature, to be used by the
+    conjugate heat transfer solver.
+
+    Example usage:
+        myInterfacePatchName
+        {
+            type                solidWallTemperatureCoupled;
+            neighbourRegionName fluid;
+            neighbourPatchName  fluidSolidInterface;
+            neighbourFieldName  T;
+            K                   K;
+            value               uniform 300;
+        }
+
+SourceFiles
+    solidWallTemperatureCoupledFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallTemperatureCoupledFvPatchScalarField_H
+#define solidWallTemperatureCoupledFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "coupleManager.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class solidWallTemperatureCoupledFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallTemperatureCoupledFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Couple manager object
+        coupleManager coupleManager_;
+
+        //- Name of thermal conductivity field
+        word KName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallTemperatureCoupled");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallTemperatureCoupledFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given solidWallTemperatureCoupledFvPatchScalarField
+        //  onto a new patch
+        solidWallTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallTemperatureCoupledFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallTemperatureCoupledFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallTemperatureCoupledFvPatchScalarField
+        (
+            const solidWallTemperatureCoupledFvPatchScalarField&,
+            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 solidWallTemperatureCoupledFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Flux
+        tmp<scalarField> flux() const;
+
+        //- 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/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..ba3d099a9f72730909512c6017df90d51d8818d5
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -0,0 +1,10 @@
+    tmp<fvVectorMatrix> UEqn = solveMomentumEquation
+    (
+        momentumPredictor,
+        Uf[i],
+        rhof[i],
+        phif[i],
+        pdf[i],
+        ghf[i],
+        turb[i]
+    );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C
new file mode 100644
index 0000000000000000000000000000000000000000..8101c81e9358e1d6a35e596eede27af3920f5b4d
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Continuity errors for fluid meshes
+
+\*---------------------------------------------------------------------------*/
+
+void compressibleContinuityErrors
+(
+    scalar& cumulativeContErr,
+    const volScalarField& rho,
+    const basicThermo& thermo
+)
+{
+    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
+
+    scalar sumLocalContErr =
+    (
+        fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
+    ).value();
+
+    scalar globalContErr =
+    (
+        fvc::domainIntegrate(rho - thermo.rho())/totalMass
+    ).value();
+
+    cumulativeContErr += globalContErr;
+
+    const word& regionName = rho.mesh().name();
+
+    Info<< "time step continuity errors (" << regionName << ")"
+        << ": sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
new file mode 100644
index 0000000000000000000000000000000000000000..ce1f54aca16ec8e0493a94882f83e9d4dc474d83
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers for the fluid
+    regions
+
+\*---------------------------------------------------------------------------*/
+
+scalar compressibleCourantNo
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volScalarField& rho,
+    const surfaceScalarField& phi
+)
+{
+    scalar CoNum = 0.0;
+    scalar meanCoNum = 0.0;
+
+    if (mesh.nInternalFaces())
+    {
+        surfaceScalarField SfUfbyDelta =
+            mesh.surfaceInterpolation::deltaCoeffs()
+          * mag(phi)
+          / fvc::interpolate(rho);
+
+        CoNum = max(SfUfbyDelta/mesh.magSf())
+            .value()*runTime.deltaT().value();
+
+        meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+            .value()*runTime.deltaT().value();
+    }
+
+    Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
+        << " max: " << CoNum << endl;
+
+    return CoNum;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
new file mode 100644
index 0000000000000000000000000000000000000000..68c3621ec1fdb4fa4b0c2d9a805f7d3afd086767
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -0,0 +1,15 @@
+    scalar CoNum = -GREAT;
+    forAll(fluidRegions, regionI)
+    {
+        CoNum = max
+        (
+            compressibleCourantNo
+            (
+                fluidRegions[regionI],
+                runTime,
+                rhof[regionI],
+                phif[regionI]
+            ),
+            CoNum
+        );
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..73c51deeecef39b0731350481be819a85218d013
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -0,0 +1,181 @@
+    // Initialise fluid field pointer lists
+    PtrList<basicThermo> thermof(fluidRegions.size());
+    PtrList<volScalarField> rhof(fluidRegions.size());
+    PtrList<volScalarField> Kf(fluidRegions.size());
+    PtrList<volVectorField> Uf(fluidRegions.size());
+    PtrList<surfaceScalarField> phif(fluidRegions.size());
+    PtrList<compressible::turbulenceModel> turb(fluidRegions.size());
+    PtrList<volScalarField> DpDtf(fluidRegions.size());
+    PtrList<volScalarField> ghf(fluidRegions.size());
+    PtrList<volScalarField> pdf(fluidRegions.size());
+
+    List<scalar> initialMassf(fluidRegions.size());
+
+    dimensionedScalar pRef("pRef", dimensionSet(1, -1, -2, 0, 0), 1.0E5);
+
+    // Populate fluid field pointer lists
+    forAll(fluidRegions, i)
+    {
+        Info<< "*** Reading fluid mesh thermophysical properties for region "
+            << fluidRegions[i].name() << nl << endl;
+
+        Info<< "    Adding to pdf\n" << endl;
+        pdf.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "pd",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to thermof\n" << endl;
+
+        thermof.set
+        (
+            i,
+            autoPtr<basicThermo>
+            (
+                basicThermo::New(fluidRegions[i])
+            ).ptr()
+        );
+
+        Info<< "    Adding to rhof\n" << endl;
+        rhof.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                thermof[i].rho()
+            )
+        );
+
+        Info<< "    Adding to Kf\n" << endl;
+        Kf.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                thermof[i].rho()*thermof[i].Cp()*thermof[i].alpha()
+            )
+        );
+
+        Info<< "    Adding to Uf\n" << endl;
+        Uf.set
+        (
+            i,
+            new volVectorField
+            (
+                IOobject
+                (
+                    "U",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to phif\n" << endl;
+        phif.set
+        (
+            i,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "phi",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                linearInterpolate(rhof[i]*Uf[i])
+                    & fluidRegions[i].Sf()
+            )
+        );
+
+        Info<< "    Adding to turb\n" << endl;
+        turb.set
+        (
+            i,
+            autoPtr<compressible::turbulenceModel>
+            (
+                compressible::turbulenceModel::New
+                (
+                    rhof[i],
+                    Uf[i],
+                    phif[i],
+                    thermof[i]
+                )
+            ).ptr()
+        );
+
+        Info<< "    Adding to DpDtf\n" << endl;
+        DpDtf.set
+        (
+            i,
+            new volScalarField
+            (
+                fvc::DDt
+                (
+                    surfaceScalarField
+                    (
+                        "phiU",
+                        phif[i]/fvc::interpolate(rhof[i])
+                    ),
+                    thermof[i].p()
+                )
+            )
+        );
+
+        const dictionary& environmentalProperties =
+            fluidRegions[i].lookupObject<IOdictionary>
+                ("environmentalProperties");
+        dimensionedVector g(environmentalProperties.lookup("g"));
+
+        Info<< "    Adding to ghf\n" << endl;
+        ghf.set
+        (
+            i,
+            new volScalarField
+            (
+                "gh",
+                 g & fluidRegions[i].C()
+            )
+        );
+
+        Info<< "    Updating p from pd\n" << endl;
+        thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef;
+
+
+        initialMassf[i] = fvc::domainIntegrate(rhof[i]).value();
+    }
+
+
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
new file mode 100644
index 0000000000000000000000000000000000000000..aec08349eca1f3d06d53f0f22b36d9ab19e2f3a7
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
@@ -0,0 +1,45 @@
+    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
+
+    forAll(rp.fluidRegionNames(), i)
+    {
+        Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
+            << " for time = " << runTime.timeName() << nl << endl;
+
+        fluidRegions.set
+        (
+            i,
+            new fvMesh
+            (
+                IOobject
+                (
+                    rp.fluidRegionNames()[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
+                )
+            )
+        );
+
+        // Force calculation of geometric properties to prevent it being done
+        // later in e.g. some boundary evaluation
+        //(void)fluidRegions[i].weights();
+        //(void)fluidRegions[i].deltaCoeffs();
+
+        // Attach environmental properties to each region
+        autoPtr<IOdictionary> environmentalProperties
+        (
+            new IOdictionary
+            (
+                IOobject
+                (
+                    "environmentalProperties",
+                    runTime.constant(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+        environmentalProperties.ptr()->store();
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..301ceddfdb432b980642c092d15257e87f51e71e
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
@@ -0,0 +1,9 @@
+    solveEnthalpyEquation
+    (
+        rhof[i],
+        DpDtf[i],
+        phif[i],
+        turb[i],
+        thermof[i]
+    );
+
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..9e07ab170ca34dcea523343a2fdd380976d1af35
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -0,0 +1,66 @@
+{
+    bool closedVolume = false;
+
+    //pdf[i].boundaryField() ==
+    //    thermof[i].p().boundaryField()
+    //  - rhof[i].boundaryField()*ghf[i].boundaryField()
+    //  - pRef.value();
+
+    rhof[i] = thermof[i].rho();
+
+    volScalarField rUA = 1.0/UEqn().A();
+    Uf[i] = rUA*UEqn().H();
+
+    phif[i] =
+        fvc::interpolate(rhof[i])
+       *(
+            (fvc::interpolate(Uf[i]) & fluidRegions[i].Sf())
+//          + fvc::ddtPhiCorr(rUA, rhof[i], Uf[i], phif[i])
+        )
+      - fvc::interpolate(rhof[i]*rUA*ghf[i])
+            *fvc::snGrad(rhof[i])
+            *fluidRegions[i].magSf();
+
+    // Solve pressure difference
+#   include "pdEqn.H"
+
+    // Solve continuity
+#   include "rhoEqn.H"
+
+    // Update pressure field (including bc)
+    thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef;
+    DpDtf[i] = fvc::DDt
+    (
+        surfaceScalarField("phiU", phif[i]/fvc::interpolate(rhof[i])),
+        thermof[i].p()
+    );
+
+    // Update continuity errors
+    compressibleContinuityErrors(cumulativeContErr, rhof[i], thermof[i]);
+
+    // Correct velocity field
+    Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]);
+    Uf[i].correctBoundaryConditions();
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        thermof[i].p() +=
+            (
+                dimensionedScalar
+                (
+                    "massIni",
+                    dimMass,
+                    initialMassf[i]
+                )
+              - fvc::domainIntegrate(thermof[i].psi()*thermof[i].p())
+            )
+            /fvc::domainIntegrate(thermof[i].psi());
+
+        rhof[i] = thermof[i].rho();
+    }
+
+    // Update thermal conductivity
+    Kf[i] = rhof[i]*thermof[i].Cp()*turb[i].alphaEff();
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..a3938432565cf675aa3dee0caee3a19eaafbdb14
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H
@@ -0,0 +1,14 @@
+    solvePressureDifferenceEquation
+    (
+        corr,
+        nCorr,
+        nNonOrthCorr,
+        closedVolume,
+        pdf[i],
+        pRef,
+        rhof[i],
+        thermof[i].psi(),
+        rUA,
+        ghf[i],
+        phif[i]
+    );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..75d6d96d1119f2630f6973a626959c54b7928407
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H
@@ -0,0 +1,27 @@
+    dictionary piso = fluidRegions[i].solutionDict().subDict("PISO");
+
+    int nCorr(readInt(piso.lookup("nCorrectors")));
+
+    int nNonOrthCorr = 0;
+    if (piso.found("nNonOrthogonalCorrectors"))
+    {
+        nNonOrthCorr = readInt(piso.lookup("nNonOrthogonalCorrectors"));
+    }
+
+    bool momentumPredictor = true;
+    if (piso.found("momentumPredictor"))
+    {
+        momentumPredictor = Switch(piso.lookup("momentumPredictor"));
+    }
+
+    bool transonic = false;
+    if (piso.found("transonic"))
+    {
+        transonic = Switch(piso.lookup("transonic"));
+    }
+
+    int nOuterCorr = 1;
+    if (piso.found("nOuterCorrectors"))
+    {
+        nOuterCorr = readInt(piso.lookup("nOuterCorrectors"));
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..a7e2a98c4aa5251ea9048f47d156dbbe4f628f3e
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H
@@ -0,0 +1 @@
+    solveContinuityEquation(rhof[i], phif[i]);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..161eb742e76f2e7f798ae59cbaf6c50aff0b68da
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H
@@ -0,0 +1,15 @@
+    if (adjustTimeStep)
+    {
+        if (CoNum > SMALL)
+        {
+            runTime.setDeltaT
+            (
+                min
+                (
+                    maxCo*runTime.deltaT().value()/CoNum,
+                    maxDeltaT
+                )
+            );
+        }
+    }
+
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C
new file mode 100644
index 0000000000000000000000000000000000000000..3d834a32aeda2c5f435b43c864517a30f7194fa7
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C
@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Solve continuity equation
+
+\*---------------------------------------------------------------------------*/
+
+void solveContinuityEquation
+(
+    volScalarField& rho,
+    const surfaceScalarField& phi
+)
+{
+    solve(fvm::ddt(rho) + fvc::div(phi));
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C
new file mode 100644
index 0000000000000000000000000000000000000000..9c9ba030c324abe46f988055018cb87f5c7553b0
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C
@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Solve enthalpy equation
+
+\*---------------------------------------------------------------------------*/
+
+void solveEnthalpyEquation
+(
+    const volScalarField& rho,
+    const volScalarField& DpDt,
+    const surfaceScalarField& phi,
+    const compressible::turbulenceModel& turb,
+    basicThermo& thermo
+)
+{
+    volScalarField& h = thermo.h();
+
+    tmp<fvScalarMatrix> hEqn
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turb.alphaEff(), h)
+     ==
+        DpDt
+    );
+    hEqn().relax();
+    hEqn().solve();
+
+    thermo.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T())
+        << endl;
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
new file mode 100644
index 0000000000000000000000000000000000000000..e24771256f48c6232a12e05d30438e3a42f3c3e0
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@@ -0,0 +1,15 @@
+#   include "rhoEqn.H"
+    for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
+    {
+    #   include "UEqn.H"
+
+    #   include "hEqn.H"
+
+     // --- PISO loop
+
+        for (int corr=0; corr<nCorr; corr++)
+        {
+    #       include "pEqn.H"
+        }
+    }
+    turb[i].correct();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C
new file mode 100644
index 0000000000000000000000000000000000000000..22a8b2e65ebd81578dd5631cff178dd0743107f9
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Solve momentum equation and return matrix for use in pressure equation
+
+\*---------------------------------------------------------------------------*/
+
+tmp<fvVectorMatrix> solveMomentumEquation
+(
+    const bool momentumPredictor,
+    volVectorField& U,
+    const volScalarField& rho,
+    const surfaceScalarField& phi,
+    const volScalarField& pd,
+    const volScalarField& gh,
+    const compressible::turbulenceModel& turb
+)
+{
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turb.divDevRhoReff(U)
+    );
+
+    UEqn().relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh);
+    }
+
+    return UEqn;
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C
new file mode 100644
index 0000000000000000000000000000000000000000..3dbd6c82646f444800cf728f0e3ca2d65af57450
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Solve pressure difference equation
+
+\*---------------------------------------------------------------------------*/
+
+void solvePressureDifferenceEquation
+(
+    const label corr,
+    const label nCorr,
+    const label nNonOrthCorr,
+    bool& closedVolume,
+    volScalarField& pd,
+    const dimensionedScalar& pRef,
+    const volScalarField& rho,
+    const volScalarField& psi,
+    const volScalarField& rUA,
+    const volScalarField& gh,
+    surfaceScalarField& phi
+)
+{
+    closedVolume = pd.needReference();
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pdEqn
+        (
+            fvm::ddt(psi, pd)
+          + fvc::ddt(psi)*pRef
+          + fvc::ddt(psi, rho)*gh
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, pd)
+        );
+
+        //pdEqn.solve();
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        {
+            pdEqn.solve(pd.mesh().solver(pd.name() + "Final"));
+        }
+        else
+        {
+            pdEqn.solve(pd.mesh().solver(pd.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pdEqn.flux();
+        }
+    }
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.C
new file mode 100644
index 0000000000000000000000000000000000000000..8e70c6f07ed5af2248cb1cc0144d208dd492608b
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.C
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionProperties.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionProperties::regionProperties(const Time& runTime)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "regionProperties",
+            runTime.time().constant(),
+            runTime.db(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    fluidRegionNames_(lookup("fluidRegionNames")),
+    solidRegionNames_(lookup("solidRegionNames"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionProperties::~regionProperties()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::List<Foam::word>& Foam::regionProperties::fluidRegionNames() const
+{
+    return fluidRegionNames_;
+}
+
+
+const Foam::List<Foam::word>& Foam::regionProperties::solidRegionNames() const
+{
+    return solidRegionNames_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..86e00a5476627ae17340668167f43c633293f0e8
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/regionProperties/regionProperties.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    regionProperties
+
+Description
+    Simple class to hold region information for coupled region simulations
+
+SourceFiles
+    regionProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionProperties_H
+#define regionProperties_H
+
+#include "IOdictionary.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class regionProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionProperties
+:
+    public IOdictionary
+{
+    // Private data
+
+        //- List of the fluid region names
+        List<word> fluidRegionNames_;
+
+        //- List of the solid region names
+        List<word> solidRegionNames_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        regionProperties(const regionProperties&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const regionProperties&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        regionProperties(const Time& runTime);
+
+
+    // Destructor
+
+        ~regionProperties();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const reference to the list of fluid region names
+            const List<word>& fluidRegionNames() const;
+
+            //- Return const reference to the list of solid region names
+            const List<word>& solidRegionNames() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..3361a89add381ed2f43d39b602b2b945040453cf
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
@@ -0,0 +1,91 @@
+    // Initialise solid field pointer lists
+    PtrList<volScalarField> rhos(solidRegions.size());
+    PtrList<volScalarField> cps(solidRegions.size());
+    PtrList<volScalarField> rhosCps(solidRegions.size());
+    PtrList<volScalarField> Ks(solidRegions.size());
+    PtrList<volScalarField> Ts(solidRegions.size());
+
+    // Populate solid field pointer lists
+    forAll(solidRegions, i)
+    {
+        Info<< "*** Reading solid mesh thermophysical properties for region "
+            << solidRegions[i].name() << nl << endl;
+
+        Info<< "    Adding to rhos\n" << endl;
+        rhos.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to cps\n" << endl;
+        cps.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "cp",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        rhosCps.set
+        (
+            i,
+            new volScalarField("rhosCps", rhos[i]*cps[i])
+        );
+
+        Info<< "    Adding to Ks\n" << endl;
+        Ks.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to Ts\n" << endl;
+        Ts.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "T",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
new file mode 100644
index 0000000000000000000000000000000000000000..eb50be23808e6ffdf31c17fca398bb5366f24c7e
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
@@ -0,0 +1,27 @@
+    PtrList<fvMesh> solidRegions(rp.solidRegionNames().size());
+
+    forAll(rp.solidRegionNames(), i)
+    {
+        Info<< "Create solid mesh for region " << rp.solidRegionNames()[i]
+            << " for time = " << runTime.timeName() << nl << endl;
+
+        solidRegions.set
+        (
+            i,
+            new fvMesh
+            (
+                IOobject
+                (
+                    rp.solidRegionNames()[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
+                )
+            )
+        );
+
+        // Force calculation of geometric properties to prevent it being done
+        // later in e.g. some boundary evaluation
+        //(void)solidRegions[i].weights();
+        //(void)solidRegions[i].deltaCoeffs();
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..6373e79d49c0869238a30a138aef81a7af905d05
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
@@ -0,0 +1,7 @@
+    dictionary piso = solidRegions[i].solutionDict().subDict("PISO");
+
+    int nNonOrthCorr = 0;
+    if (piso.found("nNonOrthogonalCorrectors"))
+    {
+        nNonOrthCorr = readInt(piso.lookup("nNonOrthogonalCorrectors"));
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
new file mode 100644
index 0000000000000000000000000000000000000000..40299e7b575315ebf474e119eefff142aec054a0
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -0,0 +1,9 @@
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        solve
+        (
+            fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i])
+        );
+    }
+}
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index 12e03c3c8bfacf88540761a88b499c474c10c1cc..b3046c7fbb9a55d16fbe357b67995a068b8f2364 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -46,18 +46,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createDynamicFvMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,10 +64,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readControls.H"
-#       include "CourantNo.H"
+        #include "readControls.H"
+        #include "CourantNo.H"
 
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -97,7 +96,7 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && correctPhi)
         {
-#           include "correctPhi.H"
+            #include "correctPhi.H"
         }
 
         // Make the fluxes relative to the mesh motion
@@ -108,22 +107,22 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && checkMeshCourantNo)
         {
-#           include "meshCourantNo.H"
+            #include "meshCourantNo.H"
         }
 
         twoPhaseProperties.correct();
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
 
         p = pd + rho*gh;
 
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 239779ce2e4f4085991b01c1c4badacf0bb08133..887bcbdead1a56701f879de1d94aafdf4da40463 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -45,7 +45,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -83,14 +83,29 @@
     );
 
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    surfaceScalarField ghf("gh", g & mesh.Cf());
-
-
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 89b7ecc5cda083eab519b1c49f9ca8110bb20ebb..f6931116801d86209597f07eceda5e6cbd430ae4 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -45,18 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,10 +63,10 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -75,17 +74,19 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
+
+        p = pd + rho*gh;
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/lesInterFoam/createFields.H b/applications/solvers/multiphase/lesInterFoam/createFields.H
index 666e65587146e6b66b09fa1bcd50f282f2dba8a9..c24f8bda35ce83ab699b266d7b3b4c6f87c6f54a 100644
--- a/applications/solvers/multiphase/lesInterFoam/createFields.H
+++ b/applications/solvers/multiphase/lesInterFoam/createFields.H
@@ -44,7 +44,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -87,9 +87,24 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
 
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
 
diff --git a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
index c006fc6dd6e0fd47c52c20ddc37324b650e3c6f0..acb11e1b294a2e2426def305f29f2c4626fb6f9b 100644
--- a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
+++ b/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C
@@ -26,8 +26,11 @@ Application
     lesInterFoam
 
 Description
-    Solver for 2 incompressible fluids capturing the interface.  Turbulence is
-    modelled using a runtime selectable incompressible LES model.
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.  Turbulence is modelled using a run-time
+    selectable incompressible LES model.
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,28 +41,21 @@ Description
 #include "twoPhaseMixture.H"
 #include "incompressible/LESmodel/LESmodel.H"
 
-#include "IFstream.H"
-#include "OFstream.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-
-#   include "createFields.H"
-//#   include "createAverages.H"
-
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,31 +63,29 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
         turbulence->correct();
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr < nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
-//#       include "calculateAverages.H"
+        #include "continuityErrs.H"
 
         runTime.write();
-//#       include "writeNaveragingSteps.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/multiphase/rasInterFoam/createFields.H b/applications/solvers/multiphase/rasInterFoam/createFields.H
index 1aef6db57959247d700d334baa0369f45446ca1a..988e8b3f744843a8c2cbad99520cff0c1c938821 100644
--- a/applications/solvers/multiphase/rasInterFoam/createFields.H
+++ b/applications/solvers/multiphase/rasInterFoam/createFields.H
@@ -44,7 +44,7 @@
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-    
+
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
 
@@ -88,9 +88,24 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
 
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+
     // Construct interface from gamma distribution
     interfaceProperties interface(gamma, U, twoPhaseProperties);
 
diff --git a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C b/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
index 9cf47c5f3208e057b6eb5b3b55d8c2e0c24733e4..1c78f70bd82ce4aa694919faa684f048e09e8b37 100644
--- a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
+++ b/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C
@@ -23,11 +23,14 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    lesInterFoam
+    rasInterFoam
 
 Description
-    Solver for 2 incompressible fluids capturing the interface.  Turbulence is
-    modelled using a runtime selectable incompressible RAS model.
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.  Turbulence is modelled using a run-time
+    selectable incompressible RAS model.
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,18 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-#   include "correctPhi.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,26 +63,28 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readTimeControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "gammaEqnSubCycle.H"
+        #include "gammaEqnSubCycle.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr < nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "continuityErrs.H"
+        #include "continuityErrs.H"
+
+        p = pd + rho*gh;
 
         turbulence->correct();
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
index 8dc0dc5db056b548826fa664211bb9beccfa1e46..8a487e60487faedc33383dd2b8d831d42f1f5291 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
@@ -168,10 +168,28 @@ displacementInterpolationFvMotionSolver
             label zoneI = fZones.findZoneID(zoneName);
             const faceZone& fz = fZones[zoneI];
 
-            scalarField fzCoords = fz().localPoints().component(dir);
+            scalar minCoord = VGREAT;
+            scalar maxCoord = -VGREAT;
 
-            zoneCoordinates[2*i] = gMin(fzCoords);
-            zoneCoordinates[2*i+1] = gMax(fzCoords);
+            forAll(fz().meshPoints(), localI)
+            {
+                label pointI = fz().meshPoints()[localI];
+                const scalar coord = points0_[pointI][dir];
+                minCoord = min(minCoord, coord);
+                maxCoord = max(maxCoord, coord);
+            }
+
+            zoneCoordinates[2*i] = returnReduce(minCoord, minOp<scalar>());
+            zoneCoordinates[2*i+1] = returnReduce(maxCoord, maxOp<scalar>());
+
+            if (debug)
+            {
+                Pout<< "direction " << dir << " : "
+                    << "zone " << zoneName
+                    << " ranges from coordinate " << zoneCoordinates[2*i]
+                    << " to " << zoneCoordinates[2*i+1]
+                    << endl;
+            }
         }
         zoneCoordinates.sort();
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
index 00f3fa9362a5592d6926aecfd9390f28f5a56837..876b13338341c6d07b1ae16eec590c646dbd74cf 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
@@ -36,34 +36,20 @@ License
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-//Foam::vector
-// Foam::surfaceSlipDisplacementPointPatchVectorField::projectNormal
-//(
-//    const label pointI,
-//    const point& pt
-//)
-//{
-//    vector projectNormal
-//
-//    if (wedge)
-//    {
-//        // For wedge: assume axis runs from (0 0 0) in direction axisNormal_
-//        //vector projectNormal(pt-axisPt);
-//        //projectNormal -= (axisNormal&projectNormal)*projectNormal;
-//
-//        projectNormal = pt - (axisNormal_&pt)*axisNormal_;
-//        projectNormal /= mag(projectNormal) + VSMALL;
-//    }
-//    else
-//    {
-//        projectNormal = this->patch().pointNormals()[pointI];
-//    }
-//
-//    return projectNormal;
-//}
+template<>
+const char*
+NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>::
+names[] =
+{
+    "nearest",
+    "pointNormal",
+    "fixedNormal"
+};
+
+const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>
+    surfaceSlipDisplacementPointPatchVectorField::followModeNames_;
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -77,6 +63,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(p, iF),
     surfaceNames_(),
+    projectMode_(NEAREST),
     projectDir_(vector::zero),
     wedgePlane_(-1)
 {}
@@ -92,6 +79,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(p, iF),
     surfaceNames_(dict.lookup("projectSurfaces")),
+    projectMode_(followModeNames_.read(dict.lookup("followMode"))),
     projectDir_(dict.lookup("projectDirection")),
     wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
     frozenPointsZone_(dict.lookup("frozenPointsZone"))
@@ -109,6 +97,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(p, iF),
     surfaceNames_(ppf.surfaceNames()),
+    projectMode_(ppf.projectMode()),
     projectDir_(ppf.projectDir()),
     wedgePlane_(ppf.wedgePlane()),
     frozenPointsZone_(ppf.frozenPointsZone())
@@ -123,6 +112,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(ppf),
     surfaceNames_(ppf.surfaceNames()),
+    projectMode_(ppf.projectMode()),
     projectDir_(ppf.projectDir()),
     wedgePlane_(ppf.wedgePlane()),
     frozenPointsZone_(ppf.frozenPointsZone())
@@ -138,6 +128,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(ppf, iF),
     surfaceNames_(ppf.surfaceNames()),
+    projectMode_(ppf.projectMode()),
     projectDir_(ppf.projectDir()),
     wedgePlane_(ppf.wedgePlane()),
     frozenPointsZone_(ppf.frozenPointsZone())
@@ -184,17 +175,21 @@ void surfaceSlipDisplacementPointPatchVectorField::evaluate
     // Construct large enough vector in direction of projectDir so
     // we're guaranteed to hit something.
 
-    //- Fixed projection vector:
-    //vector n = projectDir_/mag(projectDir_);
-    //vector projectVec = n*(n&(mesh.bounds().max()-mesh.bounds().min()));
+    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
+
+    // For case of fixed projection vector:
+    vector projectVec;
+    if (projectMode_ == FIXEDNORMAL)
+    {
+        vector n = projectDir_/mag(projectDir_);
+        projectVec = projectLen*n;
+    }
 
     //- Per point projection vector:
-    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
 
     const pointField& localPoints = patch().localPoints();
     const labelList& meshPoints = patch().meshPoints();
 
-    //vectorField motionU(this->patchInternalField());
     vectorField displacement(this->patchInternalField());
 
 
@@ -234,85 +229,90 @@ void surfaceSlipDisplacementPointPatchVectorField::evaluate
         }
         else
         {
-            //point start(localPoints[i] + deltaT*motionU[i]);
             point start(points0[meshPoints[i]] + displacement[i]);
 
-            vector projectVec(projectLen*patch().pointNormals()[i]);
-
             scalar offset = 0;
-            if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
-            {
-                offset = start[wedgePlane_];
-                start[wedgePlane_] = 0;
-                projectVec[wedgePlane_] = 0;
-            }
-
             pointIndexHit intersection;
 
-            // Check if already on surface
-            surfaces().findNearest(start, sqr(SMALL), intersection);
-
-            if (intersection.hit())
+            if (projectMode_ == NEAREST)
             {
-                //Pout<< "    point:" << start << " near:" << intersection.hit()
-                //    << endl;
+                surfaces().findNearest(start, sqr(projectLen), intersection);
             }
             else
             {
-                // No nearest found. Do intersection
+                // Check if already on surface
+                surfaces().findNearest(start, sqr(SMALL), intersection);
 
-                label rightSurf0, rightSurf1;
-                pointIndexHit rightHit0, rightHit1;
-                surfaces().findNearestIntersection
-                (
-                    start,
-                    start+projectVec,
-                    rightSurf0,
-                    rightHit0,
-                    rightSurf1,
-                    rightHit1
-                );
-
-                // Do intersection
-                label leftSurf0, leftSurf1;
-                pointIndexHit leftHit0, leftHit1;
-                surfaces().findNearestIntersection
-                (
-                    start,
-                    start-projectVec,
-                    leftSurf0,
-                    leftHit0,
-                    leftSurf1,
-                    leftHit1
-                );
-
-                if (rightHit0.hit())
+                if (!intersection.hit())
                 {
-                    if (leftHit0.hit())
+                    // No nearest found. Do intersection
+
+                    if (projectMode_ == POINTNORMAL)
                     {
-                        if
-                        (
-                            magSqr(rightHit0.hitPoint()-start)
-                          < magSqr(leftHit0.hitPoint()-start)
-                        )
+                        projectVec = projectLen*patch().pointNormals()[i];
+                    }
+
+                    // Knock out any wedge component
+                    if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
+                    {
+                        offset = start[wedgePlane_];
+                        start[wedgePlane_] = 0;
+                        projectVec[wedgePlane_] = 0;
+                    }
+
+                    label rightSurf0, rightSurf1;
+                    pointIndexHit rightHit0, rightHit1;
+                    surfaces().findNearestIntersection
+                    (
+                        start,
+                        start+projectVec,
+                        rightSurf0,
+                        rightHit0,
+                        rightSurf1,
+                        rightHit1
+                    );
+
+                    // Do intersection
+                    label leftSurf0, leftSurf1;
+                    pointIndexHit leftHit0, leftHit1;
+                    surfaces().findNearestIntersection
+                    (
+                        start,
+                        start-projectVec,
+                        leftSurf0,
+                        leftHit0,
+                        leftSurf1,
+                        leftHit1
+                    );
+
+                    if (rightHit0.hit())
+                    {
+                        if (leftHit0.hit())
                         {
-                            intersection = rightHit0;
+                            if
+                            (
+                                magSqr(rightHit0.hitPoint()-start)
+                              < magSqr(leftHit0.hitPoint()-start)
+                            )
+                            {
+                                intersection = rightHit0;
+                            }
+                            else
+                            {
+                                intersection = leftHit0;
+                            }
                         }
                         else
                         {
-                            intersection = leftHit0;
+                            intersection = rightHit0;
                         }
                     }
                     else
                     {
-                        intersection = rightHit0;
-                    }
-                }
-                else
-                {
-                    if (leftHit0.hit())
-                    {
-                        intersection = leftHit0;
+                        if (leftHit0.hit())
+                        {
+                            intersection = leftHit0;
+                        }
                     }
                 }
             }
@@ -327,7 +327,6 @@ void surfaceSlipDisplacementPointPatchVectorField::evaluate
                 {
                     interPt[wedgePlane_] += offset;
                 }
-                //motionU[i] = (interPt-localPoints[i])/deltaT;
                 displacement[i] = interPt-points0[meshPoints[i]];
             }
             else
@@ -337,7 +336,7 @@ void surfaceSlipDisplacementPointPatchVectorField::evaluate
                     << "  did not find any intersection between ray from "
                     << start-projectVec << " to " << start+projectVec
                     << endl;
-            }   
+            }  
         }
     }
 
@@ -356,6 +355,8 @@ void surfaceSlipDisplacementPointPatchVectorField::write(Ostream& os) const
     pointPatchVectorField::write(os);
     os.writeKeyword("projectSurfaces") << surfaceNames_
         << token::END_STATEMENT << nl;
+    os.writeKeyword("followMode") << followModeNames_[projectMode_]
+        << token::END_STATEMENT << nl;
     os.writeKeyword("projectDirection") << projectDir_
         << token::END_STATEMENT << nl;
     os.writeKeyword("wedgePlane") << wedgePlane_
diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H
index 7db1e8ef68d44ed1d37fa7f3f1dc2fe7e8bab5fa..169a7dfb08565e5095be2254087a6e8505bbedb6 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H
+++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H
@@ -27,6 +27,21 @@ Class
 
 Description
     Displacement follows a triSurface. Use in a displacement fvMotionSolver.
+    Following is either
+    - NEAREST : nearest
+    - POINTNORMAL : intersection with point normal
+    - FIXEDNORMAL : intersection with fixed vector
+
+    Optionally (intersection only) removes a component ("wedgePlane") to
+    stay in 2D.
+
+    Needs:
+    - projectSurfaces : names of triSurfaceMeshes (in constant/triSurface)
+    - followMode : see above
+    - projectDirection : if followMode = fixedNormal
+    - wedgePlane : -1 or component to knock out of intersection normal
+    - frozenPointsZone : empty or name of pointZone containing points
+                         that do not move
 
 SourceFiles
     surfaceSlipDisplacementPointPatchVectorField.C
@@ -52,11 +67,31 @@ class surfaceSlipDisplacementPointPatchVectorField
 :
     public pointPatchVectorField
 {
+
+public:
+
+    // Public data types
+
+        enum followMode
+        {
+            NEAREST,
+            POINTNORMAL,
+            FIXEDNORMAL
+        };
+
+private:
+
     // Private data
 
+        //- follow mode names
+        static const NamedEnum<followMode, 3> followModeNames_;
+
         //- names of surfaces
         const fileNameList surfaceNames_;
 
+        //- How to follow/project onto surface
+        const followMode projectMode_;
+
         //- direction to project
         const vector projectDir_;
 
@@ -72,9 +107,6 @@ class surfaceSlipDisplacementPointPatchVectorField
 
     // Private Member Functions
 
-        ////- Calculate projection direction (normalised) at pointI.
-        //vector projectNormal(const label pointI, const point& pt) const;
-
         //- Disallow default bitwise assignment
         void operator=(const surfaceSlipDisplacementPointPatchVectorField&);
 
@@ -163,6 +195,12 @@ public:
         //- Surface to follow. Demand loads surfaceNames.
         const triSurfaceMeshes& surfaces() const;
 
+        //- Mode of projection/following
+        const followMode projectMode() const
+        {
+            return projectMode_;
+        }
+
         //- Direction to project back onto surface
         const vector& projectDir() const
         {
diff --git a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
index b9da7468d6814ac377bca118b33bcee8e961d559..e476bc8348bdd477c0a5b548b0db0d85d97436fc 100644
--- a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
+++ b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
@@ -142,7 +142,7 @@ bool Foam::surfaceToCell::differingPointNormals
                     pointToNearest
                 );
 
-            if (pointTriI != cellTriI)
+            if (pointTriI != -1 && pointTriI != cellTriI)
             {
                 scalar cosAngle = normals[pointTriI] & normals[cellTriI];