diff --git a/applications/solvers/combustion/dieselEngineFoam/createFields.H b/applications/solvers/combustion/dieselEngineFoam/createFields.H
index 9d9229cc3c1dbed27219d4e53ff49af24f78767d..6987608006944b90c490be22152d5b893dc648b0 100644
--- a/applications/solvers/combustion/dieselEngineFoam/createFields.H
+++ b/applications/solvers/combustion/dieselEngineFoam/createFields.H
@@ -13,6 +13,14 @@ PtrList<volScalarField>& Y = composition.Y();
 
 word inertSpecie(thermo.lookup("inertSpecie"));
 
+if (!composition.contains(inertSpecie))
+{
+    FatalErrorIn(args.executable())
+        << "Specified inert specie '" << inertSpecie << "' not found in "
+        << "species list. Available species:" << composition.species()
+        << exit(FatalError);
+}
+
 volScalarField rho
 (
     IOobject
diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/files b/applications/solvers/incompressible/porousSimpleFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..e72f7ed78c75d9aec66277a58f9807aee0073d75
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/Make/files
@@ -0,0 +1,3 @@
+porousSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/porousSimpleFoam
diff --git a/applications/solvers/incompressible/porousSimpleFoam/Make/options b/applications/solvers/incompressible/porousSimpleFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..f490133ba67a4764a473d9a150b2b06f839a3790
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I../simpleFoam \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..c1925c7708c8fb821515092efcc420b22f816ad1
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
@@ -0,0 +1,47 @@
+    // Construct the Momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
+      + turbulence->divDevReff(U)
+    );
+
+    UEqn().relax();
+
+    // Include the porous media resistance and solve the momentum equation
+    // either implicit in the tensorial resistance or transport using by
+    // including the spherical part of the resistance in the momentum diagonal
+
+    tmp<volScalarField> trAU;
+    tmp<volTensorField> trTU;
+
+    if (pressureImplicitPorosity)
+    {
+        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
+        pZones.addResistance(UEqn(), tTU());
+        trTU = inv(tTU());
+        trTU().rename("rAU");
+
+        volVectorField gradp = fvc::grad(p);
+
+        for (int UCorr=0; UCorr<nUCorr; UCorr++)
+        {
+            U = trTU() & (UEqn().H() - gradp);
+        }
+        U.correctBoundaryConditions();
+    }
+    else
+    {
+        pZones.addResistance(UEqn());
+
+        eqnResidual = solve
+        (
+            UEqn() == -fvc::grad(p)
+        ). initialResidual();
+
+        maxResidual = max(eqnResidual, maxResidual);
+
+        trAU = 1.0/UEqn().A();
+        trAU().rename("rAU");
+    }
diff --git a/applications/solvers/incompressible/porousSimpleFoam/createFields.H b/applications/solvers/incompressible/porousSimpleFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..04d57d0571472ea1ee00c9a8ec1e24171eefe640
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/createFields.H
@@ -0,0 +1,64 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info << "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    int nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
+        {
+            nUCorr = readInt
+            (
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
+            );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
diff --git a/applications/solvers/incompressible/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..e6797cd8603b8032e1c901583b6ddb15992b7ddc
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/pEqn.H
@@ -0,0 +1,59 @@
+if (pressureImplicitPorosity)
+{
+    U = trTU()&UEqn().H();
+}
+else
+{
+    U = trAU()*UEqn().H();
+}
+
+UEqn.clear();
+phi = fvc::interpolate(U) & mesh.Sf();
+adjustPhi(phi, U, p);
+
+for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+    tmp<fvScalarMatrix> tpEqn;
+
+    if (pressureImplicitPorosity)
+    {
+        tpEqn = (fvm::laplacian(trTU(), p) == fvc::div(phi));
+    }
+    else
+    {
+        tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phi));
+    }
+
+    tpEqn().setReference(pRefCell, pRefValue);
+    // retain the residual from the first iteration
+    if (nonOrth == 0)
+    {
+        eqnResidual = tpEqn().solve().initialResidual();
+        maxResidual = max(eqnResidual, maxResidual);
+    }
+    else
+    {
+        tpEqn().solve();
+    }
+
+    if (nonOrth == nNonOrthCorr)
+    {
+        phi -= tpEqn().flux();
+    }
+}
+
+#include "continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+if (pressureImplicitPorosity)
+{
+    U -= trTU()&fvc::grad(p);
+}
+else
+{
+    U -= trAU()*fvc::grad(p);
+}
+
+U.correctBoundaryConditions();
diff --git a/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..b78f1688fe9bacd942529d1edc031e7fce981eff
--- /dev/null
+++ b/applications/solvers/incompressible/porousSimpleFoam/porousSimpleFoam.C
@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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
+    porousSimpleFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow with
+    implicit or explicit porosity treatment
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "porousZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readSIMPLEControls.H"
+        #include "initConvergenceCheck.H"
+
+        p.storePrevIter();
+
+        // Pressure-velocity SIMPLE corrector
+        {
+            #include "UEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        #include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index 5c84cfdb3f5eae620a53dd50b8a0d32e6abb1fac..abe03f60c0693786e07694888078dd8ae0e08a28 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -14,7 +14,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for (label i=0; i<Y.size(); i++)
+    forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)
         {
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index bba39fe8627d330b08eab8f586835fbe37449e39..14d6f6f2298d26096ed33aead5eb30fac51fe886 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -13,6 +13,14 @@
 
     word inertSpecie(thermo.lookup("inertSpecie"));
 
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
index e802fdfe0d85929101f5e6723080b6666dfa2fb4..4e9c29815aaa0a96f1a017be68befcff7b154b59 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
@@ -15,7 +15,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for (label i=0; i<Y.size(); i++)
+    forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)
         {
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index a3054db0f5853192df84abadd77f6a3e3fed6687..3d6c5500ea9a3b3460eefb90c0abc4ead007f969 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
@@ -13,6 +13,14 @@
 
     word inertSpecie(thermo.lookup("inertSpecie"));
 
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
index 5cf44cb09fef5500b0813c099abc57c433cee9cc..c687f2035ba3fb4bb6d48efcd631d052eb189fce 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@@ -14,7 +14,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for (label i=0; i<Y.size(); i++)
+    forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)
         {
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index ec820fa5e364dcc3a56781d08bf9ac2dc7925ad0..359599a61f4f1bf68362720d9d9e9e006a9fb2cc 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -13,6 +13,14 @@
 
     word inertSpecie(thermo.lookup("inertSpecie"));
 
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
diff --git a/applications/solvers/multiphase/interMixingFoam/Make/files b/applications/solvers/multiphase/interMixingFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..488cd77e56dca265db6dab1b98933df8fe9ae72a
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/Make/files
@@ -0,0 +1,6 @@
+incompressibleThreePhaseMixture/threePhaseMixture.C
+threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+interMixingFoam.C
+
+EXE = $(FOAM_APPBIN)/interMixingFoam
+
diff --git a/applications/solvers/multiphase/interMixingFoam/Make/options b/applications/solvers/multiphase/interMixingFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..d8e4da2313e4343d0c6200c7fe9c7e8d2c090d48
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/Make/options
@@ -0,0 +1,17 @@
+INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
+
+EXE_INC = \
+    -I$(INTERFOAM) \
+    -IincompressibleThreePhaseMixture \
+    -IthreePhaseInterfaceProperties \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
new file mode 100644
index 0000000000000000000000000000000000000000..f9bad0a7058577487b0cb6e4fd7bb14f1e9ee24d
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
@@ -0,0 +1,164 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phir
+    (
+        IOobject
+        (
+            "phir",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
+    );
+
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        // Create the limiter to be used for all phase-fractions
+        scalarField allLambda(mesh.nFaces(), 1.0);
+
+        // Split the limiter into a surfaceScalarField
+        slicedSurfaceScalarField lambda
+        (
+            IOobject
+            (
+                "lambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allLambda
+        );
+
+
+        // Create the complete convection flux for alpha1
+        surfaceScalarField phiAlpha1 =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha3, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha1
+        surfaceScalarField phiAlpha1BD =
+            upwind<scalar>(mesh, phi).flux(alpha1);
+
+        // Calculate the flux correction for alpha1
+        phiAlpha1 -= phiAlpha1BD;
+
+        // Calculate the limiter for alpha1
+        MULES::limiter
+        (
+            allLambda,
+            oneField(),
+            alpha1,
+            phiAlpha1BD,
+            phiAlpha1,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Create the complete flux for alpha2
+        surfaceScalarField phiAlpha2 =
+            fvc::flux
+            (
+                phi,
+                alpha2,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(phir, alpha1, alpharScheme),
+                alpha2,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha2
+        surfaceScalarField phiAlpha2BD =
+            upwind<scalar>(mesh, phi).flux(alpha2);
+
+        // Calculate the flux correction for alpha2
+        phiAlpha2 -= phiAlpha2BD;
+
+        // Further limit the limiter for alpha2
+        MULES::limiter
+        (
+            allLambda,
+            oneField(),
+            alpha2,
+            phiAlpha2BD,
+            phiAlpha2,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Construct the limited fluxes
+        phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
+        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
+
+        // Solve for alpha1
+        solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
+
+        // Create the diffusion coefficients for alpha2<->alpha3
+        volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
+        volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
+
+        // Add the diffusive flux for alpha3->alpha2
+        phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
+
+        // Solve for alpha2
+        fvScalarMatrix alpha2Eqn
+        (
+            fvm::ddt(alpha2)
+          + fvc::div(phiAlpha2)
+          - fvm::laplacian(Dc23 + Dc32, alpha2)
+        );
+        alpha2Eqn.solve();
+
+        // Construct the complete mass flux
+        rhoPhi =
+              phiAlpha1*(rho1 - rho3)
+            + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+            + phi*rho3;
+
+        alpha3 = 1.0 - alpha1 - alpha2;
+    }
+
+    Info<< "Air phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha2.weightedAverage(mesh.V()).value()
+        << "  Min(alpha2) = " << min(alpha2).value()
+        << "  Max(alpha2) = " << max(alpha2).value()
+        << endl;
+}
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..765087a183baf3db62ec9b4dcf382d667b6209f3
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H
@@ -0,0 +1,43 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqns.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqns.H"
+}
+
+interface.correct();
+
+{
+    volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
+
+    //solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+    //Info<< "density error = "
+    //    << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl;
+
+    rho == rhoNew;
+}
diff --git a/applications/solvers/multiphase/interMixingFoam/createFields.H b/applications/solvers/multiphase/interMixingFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..9615e6c18645286de8cbdf87ba5c1a127b251236
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/createFields.H
@@ -0,0 +1,132 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            "alpha1",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field alpha2\n" << endl;
+    volScalarField alpha2
+    (
+        IOobject
+        (
+            "alpha2",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field alpha3\n" << endl;
+    volScalarField alpha3
+    (
+        IOobject
+        (
+            "alpha3",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    alpha3 == 1.0 - alpha1 - alpha2;
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    threePhaseMixture threePhaseProperties(U, phi);
+
+    const dimensionedScalar& rho1 = threePhaseProperties.rho1();
+    const dimensionedScalar& rho2 = threePhaseProperties.rho2();
+    const dimensionedScalar& rho3 = threePhaseProperties.rho3();
+
+    dimensionedScalar D23(threePhaseProperties.lookup("D23"));
+
+    // Need to store rho for ddt(rho, U)
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
+        alpha1.boundaryField().types()
+    );
+    rho.oldTime();
+
+
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rho*phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        rho1*phi
+    );
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    // Construct interface from alpha distribution
+    threePhaseInterfaceProperties interface(threePhaseProperties);
+
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
+    );
diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C
new file mode 100644
index 0000000000000000000000000000000000000000..f27e17923f952e5fafec44cf6536241071ad9743
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+//- Calculate and return the laminar viscosity
+void Foam::threePhaseMixture::calcNu()
+{
+    // Average kinematic viscosity calculated from dynamic viscosity
+    nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseMixture::threePhaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    transportModel(U, phi),
+
+    phase1Name_("phase1"),
+    phase2Name_("phase2"),
+    phase3Name_("phase3"),
+
+    nuModel1_
+    (
+        viscosityModel::New
+        (
+            "nu1",
+            subDict(phase1Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel2_
+    (
+        viscosityModel::New
+        (
+            "nu2",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel3_
+    (
+        viscosityModel::New
+        (
+            "nu3",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+
+    rho1_(nuModel1_->viscosityProperties().lookup("rho")),
+    rho2_(nuModel2_->viscosityProperties().lookup("rho")),
+    rho3_(nuModel3_->viscosityProperties().lookup("rho")),
+
+    U_(U),
+    phi_(phi),
+
+    alpha1_(U_.db().lookupObject<const volScalarField> ("alpha1")),
+    alpha2_(U_.db().lookupObject<const volScalarField> ("alpha2")),
+    alpha3_(U_.db().lookupObject<const volScalarField> ("alpha3")),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            U_.time().timeName(),
+            U_.db()
+        ),
+        U_.mesh(),
+        dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    calcNu();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            "mu",
+            alpha1_*rho1_*nuModel1_->nu()
+          + alpha2_*rho2_*nuModel2_->nu()
+          + alpha3_*rho3_*nuModel3_->nu()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "mu",
+            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+          + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+          + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "nu",
+            (
+                alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+              + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+              + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+            )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
+        )
+    );
+}
+
+
+bool Foam::threePhaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        if
+        (
+            nuModel1_().read(*this)
+         && nuModel2_().read(*this)
+         && nuModel3_().read(*this)
+        )
+        {
+            nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
+            nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
+            nuModel3_->viscosityProperties().lookup("rho") >> rho3_;
+
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H
new file mode 100644
index 0000000000000000000000000000000000000000..2712b5d0cbfc18500f47c9b8db5661a5c0721fa6
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H
@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+Description
+
+SourceFiles
+    threePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseMixture_H
+#define threePhaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class threePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseMixture
+:
+    public transportModel
+{
+    // Private data
+
+        word phase1Name_;
+        word phase2Name_;
+        word phase3Name_;
+
+        autoPtr<viscosityModel> nuModel1_;
+        autoPtr<viscosityModel> nuModel2_;
+        autoPtr<viscosityModel> nuModel3_;
+
+        dimensionedScalar rho1_;
+        dimensionedScalar rho2_;
+        dimensionedScalar rho3_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        const volScalarField& alpha1_;
+        const volScalarField& alpha2_;
+        const volScalarField& alpha3_;
+
+        volScalarField nu_;
+
+
+    // Private Member Functions
+
+        //- Calculate and return the laminar viscosity
+        void calcNu();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        threePhaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    // Destructor
+
+        ~threePhaseMixture()
+        {}
+
+
+    // Member Functions
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel1() const
+        {
+            return nuModel1_();
+        }
+
+        //- Return const-access to phase2 viscosityModel
+        const viscosityModel& nuModel2() const
+        {
+            return nuModel2_();
+        }
+
+        //- Return const-access to phase3 viscosityModel
+        const viscosityModel& nuModel3() const
+        {
+            return nuModel3_();
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho1() const
+        {
+            return rho1_;
+        }
+
+        //- Return const-access to phase2 density
+        const dimensionedScalar& rho2() const
+        {
+            return rho2_;
+        };
+
+        //- Return const-access to phase3 density
+        const dimensionedScalar& rho3() const
+        {
+            return rho3_;
+        };
+
+        const volScalarField& alpha1() const
+        {
+            return alpha1_;
+        }
+
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        const volScalarField& alpha3() const
+        {
+            return alpha3_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const
+        {
+            return nu_;
+        }
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Correct the laminar viscosity
+        void correct()
+        {
+            calcNu();
+        }
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..a18df8d12a59ee7d3914109374467a5a10e14a80
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    interMixingFoam
+
+Description
+    Solver for 3 incompressible fluids, two of which are miscible,
+    using a VOF method to capture the interface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "threePhaseMixture.H"
+#include "threePhaseInterfaceProperties.H"
+#include "turbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "correctPhi.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        threePhaseProperties.correct();
+
+        #include "alphaEqnsSubCycle.H"
+
+        #define twoPhaseProperties threePhaseProperties
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        #include "continuityErrs.H"
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "\n end \n";
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
new file mode 100644
index 0000000000000000000000000000000000000000..1ab10bd74348be3faafce4e5db2b66c57c5f7a05
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseInterfaceProperties.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "mathematicalConstants.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::threePhaseInterfaceProperties::correctContactAngle
+(
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& alpha1 =
+        mixture_.alpha1().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha2 =
+        mixture_.alpha2().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha3 =
+        mixture_.alpha3().boundaryField();
+    const volVectorField::GeometricBoundaryField& U =
+        mixture_.U().boundaryField();
+
+    const fvMesh& mesh = mixture_.U().mesh();
+    const fvBoundaryMesh& boundary = mesh.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& a2cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha2[patchi]);
+
+            const alphaContactAngleFvPatchScalarField& a3cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha3[patchi]);
+
+            scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
+            scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
+
+            scalarField sumTwoPhaseAlpha =
+                twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
+
+            twoPhaseAlpha2 /= sumTwoPhaseAlpha;
+            twoPhaseAlpha3 /= sumTwoPhaseAlpha;
+
+            fvsPatchVectorField& nHatp = nHatb[patchi];
+
+            scalarField theta =
+                convertToRad
+               *(
+                   twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+                 + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
+               );
+
+            vectorField nf = boundary[patchi].nf();
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12 = nHatp & nf;
+
+            scalarField b1 = cos(theta);
+
+            scalarField b2(nHatp.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det = 1.0 - a12*a12;
+
+            scalarField a = (b1 - a12*b2)/det;
+            scalarField b = (b2 - a12*b1)/det;
+
+            nHatp = a*nf + b*nHatp;
+
+            nHatp /= (mag(nHatp) + deltaN_.value());
+        }
+    }
+}
+
+
+void Foam::threePhaseInterfaceProperties::calculateK()
+{
+    const volScalarField& alpha1 = mixture_.alpha1();
+
+    const fvMesh& mesh = alpha1.mesh();
+    const surfaceVectorField& Sf = mesh.Sf();
+
+    // Cell gradient of alpha
+    volVectorField gradAlpha = fvc::grad(alpha1);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+
+    // Face unit interface normal
+    surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
+    correctContactAngle(nHatfv.boundaryField());
+
+    // Face unit interface normal flux
+    nHatf_ = nHatfv & Sf;
+
+    // Simple expression for curvature
+    K_ = -fvc::div(nHatf_);
+
+    // Complex expression for curvature.
+    // Correction is formally zero but numerically non-zero.
+    //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
+    //nHat.boundaryField() = nHatfv.boundaryField();
+    //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
+(
+    const threePhaseMixture& mixture
+)
+:
+    mixture_(mixture),
+    cAlpha_
+    (
+        readScalar
+        (
+            mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
+        )
+    ),
+    sigma12_(mixture.lookup("sigma12")),
+    sigma13_(mixture.lookup("sigma13")),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
+    ),
+
+    nHatf_
+    (
+        (
+            fvc::interpolate(fvc::grad(mixture.alpha1()))
+           /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
+        ) & mixture.alpha1().mesh().Sf()
+    ),
+
+    K_
+    (
+        IOobject
+        (
+            "K",
+            mixture.alpha1().time().timeName(),
+            mixture.alpha1().mesh()
+        ),
+        -fvc::div(nHatf_)
+    )
+{
+    calculateK();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..3a402feeedf37ebda540220c8773b76b4cb4f6aa
--- /dev/null
+++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+SourceFiles
+    threePhaseInterfaceProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseInterfaceProperties_H
+#define threePhaseInterfaceProperties_H
+
+#include "threePhaseMixture.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class threePhaseInterfaceProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseInterfaceProperties
+{
+    // Private data
+
+        const threePhaseMixture& mixture_;
+
+        //- Compression coefficient
+        scalar cAlpha_;
+
+        //- Surface tension 1-2
+        dimensionedScalar sigma12_;
+
+        //- Surface tension 1-3
+        dimensionedScalar sigma13_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        surfaceScalarField nHatf_;
+        volScalarField K_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        threePhaseInterfaceProperties(const threePhaseInterfaceProperties&);
+        void operator=(const threePhaseInterfaceProperties&);
+
+        //- Correction for the boundary condition on the unit normal nHat on
+        //  walls to produce the correct contact dynamic angle
+        //  calculated from the component of U parallel to the wall
+        void correctContactAngle
+        (
+            surfaceVectorField::GeometricBoundaryField& nHat
+        ) const;
+
+        //- Re-calculate the interface curvature
+        void calculateK();
+
+
+public:
+
+    //- Conversion factor for degrees into radians
+    static const scalar convertToRad;
+
+
+    // Constructors
+
+        //- Construct from volume fraction field alpha and IOdictionary
+        threePhaseInterfaceProperties(const threePhaseMixture& mixture);
+
+
+    // Member Functions
+
+        scalar cAlpha() const
+        {
+            return cAlpha_;
+        }
+
+        const dimensionedScalar& deltaN() const
+        {
+            return deltaN_;
+        }
+
+        const surfaceScalarField& nHatf() const
+        {
+            return nHatf_;
+        }
+
+        const volScalarField& K() const
+        {
+            return K_;
+        }
+
+        tmp<volScalarField> sigma() const
+        {
+            volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0));
+            volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0));
+
+            return
+                (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
+               /(limitedAlpha2 + limitedAlpha3 + SMALL);
+        }
+
+        tmp<volScalarField> sigmaK() const
+        {
+            return sigma()*K_;
+        }
+
+        void correct()
+        {
+            calculateK();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index cadaed0d3683cd8e16725c601992d4df953b903d..6729cab4c2e9932765ee6bc60bcd9034cd19b160 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -121,7 +121,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         return;
     }
 
-    // Far field gas molar fractions
+    // Far field carrier  molar fractions
     scalarField Xinf(Y_.size());
 
     forAll(Xinf, i)
@@ -135,7 +135,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
     // Molar fraction of far field species at particle surface
     const scalar Xsff = 1.0 - min(sum(Cs)*specie::RR*this->T_/pc_, 1.0);
 
-    // Surface gas total molar concentration
+    // Surface carrier total molar concentration
     const scalar CsTot = pc_/(specie::RR*this->T_);
 
     // Surface carrier composition (molar fraction)
@@ -171,10 +171,10 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
             cbrt(td.cloud().mcCarrierThermo().speciesData()[i].W());
 
         rhos += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].W();
-        cps += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].Cp(T);
         mus += Ys[i]*sqrtW*td.cloud().mcCarrierThermo().speciesData()[i].mu(T);
         kappa +=
             Ys[i]*cbrtW*td.cloud().mcCarrierThermo().speciesData()[i].kappa(T);
+        cps += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].Cp(T);
 
         sumYiSqrtW += Ys[i]*sqrtW;
         sumYiCbrtW += Ys[i]*cbrtW;
@@ -417,7 +417,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalarField& YComponents,
     scalarField& dMassPC,
     scalar& Sh,
-    scalar& dhsTrans,
+    scalar& dhsTrans,               // TODO: not used
     scalar& N,
     scalar& NCpW,
     scalarField& Cs
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 7851841a7486b3ea0f4ecdc775a86385dae995a5..a87342c6fe1e78e9b4ec9f02110c6977260cfd93 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -233,7 +233,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
     {
-        return  T + dt*Sh/(this->volume(d)*rho*cp);
+        return max(T + dt*Sh/(this->volume(d)*rho*cp), td.constProps().TMin());
     }
 
     const scalar As = this->areaS(d);
@@ -256,9 +256,11 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     IntegrationScheme<scalar>::integrationResult Tres =
         td.cloud().TIntegrator().integrate(T, dt, ap, bp);
 
-    dhsTrans += dt*htc*As*(Tres.average() - Tc_);
+    scalar Tnew = max(Tres.value(), td.constProps().TMin());
 
-    return Tres.value();
+    dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
+
+    return Tnew;
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index dc3b9ff61b2661400ea7c4919e5ca7f183359f20..bd2521eff9d6159c9aea8c34797a487954bfe90c 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -118,7 +118,7 @@ protected:
             scalar volumeTotal_;
 
             //- Total mass to inject [kg]
-            const scalar massTotal_;
+            scalar massTotal_;
 
             //- Total mass injected to date [kg]
             scalar massInjected_;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
index b79a8a1af1052c65079478035175df424ae2e24b..4e88b90bfebf896aab6bdad356bda3f6dfd301e3 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C
@@ -124,11 +124,12 @@ Foam::PatchInjection<CloudType>::PatchInjection
 
     label patchSize = cellOwners_.size();
     label totalPatchSize = patchSize;
-    reduce(totalPatchSize, sumOp<scalar>());
-    fraction_ = patchSize/totalPatchSize;
+    reduce(totalPatchSize, sumOp<label>());
+    fraction_ = scalar(patchSize)/totalPatchSize;
 
-    // Set total volume to inject
+    // Set total volume/mass to inject
     this->volumeTotal_ = fraction_*volumeFlowRate_().integrate(0.0, duration_);
+    this->massTotal_ *= fraction_;
 }
 
 
@@ -165,10 +166,19 @@ void Foam::PatchInjection<CloudType>::setPositionAndCell
     label& cellOwner
 )
 {
-    label cellI = this->owner().rndGen().integer(0, cellOwners_.size() - 1);
+    if (cellOwners_.size() > 0)
+    {
+        label cellI = this->owner().rndGen().integer(0, cellOwners_.size() - 1);
 
-    cellOwner = cellOwners_[cellI];
-    position = this->owner().mesh().C()[cellOwner];
+        cellOwner = cellOwners_[cellI];
+        position = this->owner().mesh().C()[cellOwner];
+    }
+    else
+    {
+        cellOwner = -1;
+        // dummy position
+        position = pTraits<vector>::max;
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
index 44d8915497bc096eedd2d22f3f74e94c8b476108..7c62bfda6eedf9e11c408c455038146cde18bf50 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
@@ -71,15 +71,12 @@ bool Foam::Rebound<CloudType>::correct
     nw /= mag(nw);
 
     scalar Un = U & nw;
-    vector Ut = U - Un*nw;
 
     if (Un > 0.0)
     {
         U -= UFactor_*2.0*Un*nw;
     }
 
-    U -= Ut;
-
     return true;
 }
 
diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files
index e950e4d4d5d50e162e74db3535bfb5764db6e0bc..6195b98a41140e7ae2da1b28c53aecf338ae50b0 100644
--- a/src/postProcessing/functionObjects/field/Make/files
+++ b/src/postProcessing/functionObjects/field/Make/files
@@ -7,10 +7,10 @@ fieldMinMax/fieldMinMax.C
 fieldMinMax/fieldMinMaxFunctionObject.C
 
 fieldValues/fieldValue/fieldValue.C
-fieldValues/face/faceSource.C
-fieldValues/face/faceSourceFunctionObject.C
-fieldValues/cell/cellSource.C
-fieldValues/cell/cellSourceFunctionObject.C
+fieldValues/faceSource/faceSource.C
+fieldValues/faceSource/faceSourceFunctionObject.C
+fieldValues/cellSource/cellSource.C
+fieldValues/cellSource/cellSourceFunctionObject.C
 
 streamLine/streamLine.C
 streamLine/streamLineParticle.C
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/IOcellSource.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/IOcellSource.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/IOcellSource.H
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/IOcellSource.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.C
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSource.H
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceFunctionObject.C
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.C
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceFunctionObject.C
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceFunctionObject.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceFunctionObject.H
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceFunctionObject.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceI.H b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceI.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceI.H
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceI.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/cell/cellSourceTemplates.C
rename to src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/IOfaceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/IOfaceSource.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/face/IOfaceSource.H
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/IOfaceSource.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
similarity index 82%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSource.C
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
index 76c1ac9515c7ee3d849cd333e2d32687bf32ca83..710e346ed5dc236baba9ab65fb11a6df968bd94e 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/face/faceSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
@@ -50,10 +50,13 @@ namespace Foam
         fieldValues::faceSource::sourceTypeNames_;
 
     template<>
-    const char* NamedEnum<fieldValues::faceSource::operationType, 4>::
-        names[] = {"none", "sum", "areaAverage", "areaIntegrate"};
+    const char* NamedEnum<fieldValues::faceSource::operationType, 5>::
+        names[] =
+        {
+            "none", "sum", "areaAverage", "areaIntegrate", "weightedAverage"
+        };
 
-    const NamedEnum<fieldValues::faceSource::operationType, 4>
+    const NamedEnum<fieldValues::faceSource::operationType, 5>
         fieldValues::faceSource::operationTypeNames_;
 
 }
@@ -68,7 +71,9 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
     if (zoneId < 0)
     {
         FatalErrorIn("faceSource::faceSource::setFaceZoneFaces()")
-            << "Unknown face zone name: " << sourceName_
+            << type() << " " << name_ << ": "
+            << sourceTypeNames_[source_] << "(" << sourceName_ << "):" << nl
+            << "    Unknown face zone name: " << sourceName_
             << ". Valid face zones are: " << mesh().faceZones().names()
             << nl << exit(FatalError);
     }
@@ -164,7 +169,9 @@ void Foam::fieldValues::faceSource::setPatchFaces()
     if (patchId < 0)
     {
         FatalErrorIn("faceSource::constructFaceAddressing()")
-            << "Unknown patch name: " << sourceName_
+            << type() << " " << name_ << ": "
+            << sourceTypeNames_[source_] << "(" << sourceName_ << "):" << nl
+            << "    Unknown patch name: " << sourceName_
             << ". Valid patch names are: "
             << mesh().boundaryMesh().names() << nl
             << exit(FatalError);
@@ -197,7 +204,7 @@ void Foam::fieldValues::faceSource::setPatchFaces()
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::fieldValues::faceSource::initialise()
+void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
 {
     switch (source_)
     {
@@ -214,15 +221,40 @@ void Foam::fieldValues::faceSource::initialise()
         default:
         {
             FatalErrorIn("faceSource::constructFaceAddressing()")
-                << "Unknown source type. Valid source types are:"
+                << type() << " " << name_ << ": "
+                << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
+                << nl << "    Unknown source type. Valid source types are:"
                 << sourceTypeNames_ << nl << exit(FatalError);
         }
     }
 
     Info<< type() << " " << name_ << ":" << nl
-        << "    total faces = " << faceId_.size() << nl
-        << "    total area  = " << sum(filterField(mesh().magSf()))
-        << nl << endl;
+        << "    total faces  = " << faceId_.size() << nl
+        << "    total area   = " << sum(filterField(mesh().magSf())) << nl;
+
+    if (operation_ == opWeightedAverage)
+    {
+        dict.lookup("weightField") >> weightFieldName_;
+        if
+        (
+            obr().foundObject<volScalarField>(weightFieldName_)
+         || obr().foundObject<surfaceScalarField>(weightFieldName_)
+        )
+        {
+            Info<< "    weight field = " << weightFieldName_;
+        }
+        else
+        {
+            FatalErrorIn("faceSource::constructFaceAddressing()")
+                << type() << " " << name_ << ": "
+                << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
+                << nl << "    Weight field " << weightFieldName_
+                << " must be either a " << volScalarField::typeName << " or "
+                << surfaceScalarField::typeName << nl << exit(FatalError);
+        }
+    }
+
+    Info<< nl << endl;
 }
 
 
@@ -302,12 +334,13 @@ Foam::fieldValues::faceSource::faceSource
     faceId_(),
     facePatchId_(),
     flipMap_(),
-    outputFilePtr_(NULL)
+    outputFilePtr_(NULL),
+    weightFieldName_("undefinedWeightedFieldName")
 {
-    initialise();
-
     if (active_)
     {
+        initialise(dict);
+
         // Create the output file if not already created
         makeFile();
     }
@@ -327,7 +360,7 @@ void Foam::fieldValues::faceSource::read(const dictionary& dict)
     if (active_)
     {
         fieldValue::read(dict);
-        initialise();
+        initialise(dict);
     }
 }
 
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
similarity index 93%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSource.H
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
index fa89fcdd6c58db55621872a90e47f3444fe9c8d3..3509d733cd0f61388dfb649f0c5de68735fc32c4 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/face/faceSource.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
@@ -39,7 +39,7 @@ Description
         valueOutput     true;       // Write values at run-time output times?
         source          faceZone;   // Type of face source: faceZone, patch
         sourceName      f0;
-        operation       sum;        // none, sum, areaAverage, areaIntegrate
+        operation       sum;
         fields
         (
             p
@@ -48,6 +48,13 @@ Description
         );
     }
 
+    where operation is one of:
+      - none
+      - sum
+      - areaAverage
+      - areaIntegrate
+      - weightedAverage
+
 SourceFiles
     faceSource.C
 
@@ -100,11 +107,12 @@ public:
             opNone,
             opSum,
             opAreaAverage,
-            opAreaIntegrate
+            opAreaIntegrate,
+            opWeightedAverage
         };
 
         //- Operation type names
-        static const NamedEnum<operationType, 4> operationTypeNames_;
+        static const NamedEnum<operationType, 5> operationTypeNames_;
 
 
 private:
@@ -143,11 +151,14 @@ protected:
         //- Output file pointer
         autoPtr<OFstream> outputFilePtr_;
 
+        //- Weight field name - only used for opWeightedAverage mode
+        word weightFieldName_;
+
 
     // Protected member functions
 
         //- Initialise, e.g. face addressing
-        void initialise();
+        void initialise(const dictionary& dict);
 
         //- Insert field values into values list
         template<class Type>
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceFunctionObject.C
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.C
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceFunctionObject.C
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceFunctionObject.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSourceFunctionObject.H
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceFunctionObject.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSourceI.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H
similarity index 100%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSourceI.H
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceI.H
diff --git a/src/postProcessing/functionObjects/field/fieldValues/face/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
similarity index 82%
rename from src/postProcessing/functionObjects/field/fieldValues/face/faceSourceTemplates.C
rename to src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
index ced2ee82eed1841f7dfd3c02fa8a2d99895348a2..d7609573a72cb18d4b7cbeb6c0cee65e67459654 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/face/faceSourceTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
@@ -129,6 +129,47 @@ Type Foam::fieldValues::faceSource::processValues
             result = sum(values*filterField(mesh().magSf()));
             break;
         }
+        case opWeightedAverage:
+        {
+            if (mesh().foundObject<volScalarField>(weightFieldName_))
+            {
+                tmp<scalarField> wField =
+                    filterField
+                    (
+                        mesh().lookupObject<volScalarField>(weightFieldName_)
+                    );
+               result = sum(values*wField())/sum(wField());
+            }
+            else if (mesh().foundObject<surfaceScalarField>(weightFieldName_))
+            {
+                tmp<scalarField> wField =
+                    filterField
+                    (
+                        mesh().lookupObject<surfaceScalarField>
+                        (
+                            weightFieldName_
+                        )
+                    );
+               result = sum(values*wField())/sum(wField());
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "fieldValues::faceSource::processValues"
+                    "("
+                        "List<Type>&"
+                    ") const"
+                )   << type() << " " << name_ << ": "
+                    << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
+                    << nl
+                    << "    Weight field " << weightFieldName_
+                    << " must be either a " << volScalarField::typeName
+                    << " or " << surfaceScalarField::typeName << nl
+                    << abort(FatalError);
+            }
+            break;
+        }
         default:
         {
             // Do nothing
diff --git a/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H b/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H
index cdedaae131a9f65b5a03971f02a193bf876b4d2d..3d862b1f6be9303ec7ba21fb5d524943ad93ba04 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H
@@ -137,7 +137,7 @@ public:
             //- Return the output field values flag
             const Switch& valueOutput() const;
 
-            //- Helper funvction to return the reference to the mesh
+            //- Helper function to return the reference to the mesh
             const fvMesh& mesh() const;
 
 
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index f575e970b4aed986ace1d0977970b381b4673d64..21e2171ec060931aa865e224ad69e953b835d473 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -63,10 +63,7 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     cpPolynomial_(pt.cpPolynomial_),
     dhPolynomial_(pt.dhPolynomial_),
     dsPolynomial_(pt.dsPolynomial_)
-{
-    // Offset dh poly so that it is relative to the enthalpy at Tstd
-    dhPolynomial_[0] -= dhPolynomial_.evaluate(specie::Tstd);
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index 217604e1335b74050e526d42c194e6ef012bffad..5d68ecf6710789f42fb4c2c1489bfda036e743a1 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -56,7 +56,7 @@ inline Foam::janafThermo<equationOfState>::janafThermo
 template<class equationOfState>
 inline void Foam::janafThermo<equationOfState>::checkT(const scalar T) const
 {
-    if (T <  Tlow_ || T > Thigh_)
+    if (T < Tlow_ || T > Thigh_)
     {
         FatalErrorIn
         (
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0 b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0
new file mode 120000
index 0000000000000000000000000000000000000000..f1c4a884b51ae4171513426c3b7f5f585d0d889e
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/0
@@ -0,0 +1 @@
+../angledDuctImplicit/0
\ No newline at end of file
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant
new file mode 120000
index 0000000000000000000000000000000000000000..28205c782b36471b118c5be2948a56e345961ad4
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/constant
@@ -0,0 +1 @@
+../angledDuctImplicit/constant
\ No newline at end of file
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..9ffe4063cf457efe9e2a167ad38617dc42f6dd8b
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPorousSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         200;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..a70c35e1bf81ab94f439453dfac9b9349f09bcae
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(U)         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,U)      Gauss upwind;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,k)      Gauss upwind;
+}
+
+laplacianSchemes
+{
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(rAU,p) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(1,p)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..b7d4dbd573caabf8c8e753563360f72f2c29af0b
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSolution
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration off;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    p               0.3;
+    U               0.7;
+    k               0.9;
+    epsilon         0.9;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..e1df94398773b013813d4f9f37769c8587889a1c
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 293;
+
+boundaryField
+{
+
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..88f509425865e94093e76cca468b5839e437ab77
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/U
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    front
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    back
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    wall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    porosityWall
+    {
+        type            slip;
+        value           uniform (0 0 0);
+    }
+    inlet
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.1;
+        value           uniform (0 0 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..d43279e3d7b721b6783cb6f10ccc0c884fde38c0
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/epsilon
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    front
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    back
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    wall
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    porosityWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 200;
+    }
+
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..8a7e20d154faf52c2bf200cd7c0bf3e66f42d817
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/k
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    front
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    back
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    wall
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    porosityWall
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..de1c9477df5f7ccb984f523c27282804fe92337f
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/nut
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    back
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    wall
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    porosityWall
+    {
+        type            nutWallFunction;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..4794b84b189f7020d3f1bc09e5d43bbf914d7b45
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/0/p
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..48cf724e269cffbb00c1dc8bcf837a8a6af5543d
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0438819b2509cf399923c3a2ce92ea0999735aea
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict
@@ -0,0 +1,123 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+                             
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -150  0  -25 )  // pt 0 (in1b) 
+    ( -150 35.35533906  -25 ) // pt 1 (in2b) 
+    ( -150  0  25 )  // pt 2 (in1f) 
+    ( -150 35.35533906  25 ) // pt 3 (in2f) 
+
+    // join inlet->outlet
+    (  0 0  -25 )    // pt 4 (join1b) 
+    ( -35.35533906   35.35533906  -25 ) // pt 5 (join2b) 
+    (  0 0  25 )    // pt 6 (join1f) 
+    ( -35.35533906   35.35533906  25 ) // pt 7 (join2f) 
+
+    // porosity ends ->outlet
+    ( 70.71067812 70.71067812  -25 )  // pt 8 (poro1b) 
+    ( 35.35533906 106.06601718  -25 )  // pt 9 (poro2b) 
+    ( 70.71067812 70.71067812  25 )  // pt 10 (poro1f) 
+    ( 35.35533906 106.06601718  25 )  // pt 11 (poro2f) 
+
+    // outlet
+    ( 141.42135624 141.42135624 -25 ) // pt 12 (out1b) 
+    ( 106.06601718 176.7766953 -25 ) // pt 13 (out2b) 
+    ( 141.42135624 141.42135624 25 ) // pt 14 (out1f) 
+    ( 106.06601718 176.7766953 25 ) // pt 15 (out2f) 
+);
+
+blocks
+(
+    // inlet block
+    hex (0 4 5 1 2 6 7 3)
+    inlet ( 15 20 20 ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex (4 8 9 5 6 10 11 7)
+    porosity ( 20 20 20 ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex (8 12 13 9 10 14 15 11)
+    outlet ( 20 20 20 )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    (2 6 7 3)
+    // outlet block
+    (10 14 15 11)
+    )
+
+    wall back
+    (
+    // inlet block
+    (1 5 4 0)
+    // outlet block
+    (9 13 12 8)
+    )
+
+    wall wall
+    (
+    // inlet block
+    (2 0 4 6)
+    (7 5 1 3)
+    // outlet block
+    (10 8 12 14)
+    (15 13 9 11)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    (6 10 11 7)
+    // porosity block
+    (5 9 8 4)
+    // porosity block
+    (6 4 8 10)
+    (11 9 5 7)
+    )
+
+    patch inlet
+    (
+    (3 1 0 2)
+    )
+
+    patch outlet
+    (
+    (15 13 12 14)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4 b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000000000000000000000000000000000..6d6d0669392b2396160799c4747d2e9f6c2c54d2
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,165 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    `format'      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+dnl> -----------------------------------------------------------------
+dnl> <STANDARD DEFINTIONS>
+dnl>
+changecom(//)changequote([,]) dnl>
+define(calc, [esyscmd(perl -e 'print ($1)')]) dnl>
+define(VCOUNT, 0)  dnl>
+define(vlabel, [[// ]pt VCOUNT ($1) define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])  dnl>
+dnl>
+define(hex2D, hex ($1b $2b $3b $4b $1f $2f $3f $4f)) dnl>
+define(quad2D, ($1f $1b $2b $2f))  dnl>
+define(frontQuad, ($1f $2f $3f $4f)) dnl>
+define(backQuad, ($4b $3b $2b $1b)) dnl>
+dnl>
+dnl> </STANDARD DEFINTIONS>
+dnl> -----------------------------------------------------------------
+dnl>
+define(ncells, 20) dnl>
+define(ninlet, 15) dnl>
+define(nporo, 20) dnl>
+define(noutlet, 20) dnl>
+dnl>
+define(x0,0) dnl>
+define(y0,0) dnl>
+define(y0,0) dnl>
+define(Cos,0.7071067812) dnl>   == cos(45)
+define(Sin,0.7071067812) dnl>   == sin(45)
+dnl>
+define(width,50) dnl>
+define(zBack,calc(-width/2)) dnl>
+define(zFront,calc(width/2)) dnl>
+define(leninlet,150)dnl>
+define(lenporo,100)dnl>
+define(lenoutlet,100)dnl>
+dnl>
+define(xhyp,calc(Sin*width)) dnl>
+define(yhyp,calc(Cos*width)) dnl>
+define(xinlet,leninlet)dnl>
+define(xporo,calc(Cos*lenporo)) dnl>
+define(yporo,calc(Sin*lenporo)) dnl>
+define(xoutlet,calc(xporo + Cos*lenoutlet)) dnl>
+define(youtlet,calc(yporo + Sin*lenoutlet)) dnl>
+dnl>
+
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -xinlet  y0  zBack )  vlabel(in1b)
+    ( -xinlet yhyp  zBack ) vlabel(in2b)
+    ( -xinlet  y0  zFront )  vlabel(in1f)
+    ( -xinlet yhyp  zFront ) vlabel(in2f)
+
+    // join inlet->outlet
+    (  x0 y0  zBack )    vlabel(join1b)
+    ( -xhyp   yhyp  zBack ) vlabel(join2b)
+    (  x0 y0  zFront )    vlabel(join1f)
+    ( -xhyp   yhyp  zFront ) vlabel(join2f)
+
+    // porosity ends ->outlet
+    ( xporo yporo  zBack )  vlabel(poro1b)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zBack )  vlabel(poro2b)
+    ( xporo yporo  zFront )  vlabel(poro1f)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zFront )  vlabel(poro2f)
+
+    // outlet
+    ( xoutlet youtlet zBack ) vlabel(out1b)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zBack ) vlabel(out2b)
+    ( xoutlet youtlet zFront ) vlabel(out1f)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zFront ) vlabel(out2f)
+);
+
+blocks
+(
+    // inlet block
+    hex2D(in1, join1, join2, in2)
+    inlet ( ninlet ncells ncells ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex2D(join1, poro1, poro2, join2)
+    porosity ( nporo ncells ncells ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex2D(poro1, out1, out2, poro2)
+    outlet ( noutlet ncells ncells )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    frontQuad(in1, join1, join2, in2)
+    // outlet block
+    frontQuad(poro1, out1, out2, poro2)
+    )
+
+    wall back
+    (
+    // inlet block
+    backQuad(in1, join1, join2, in2)
+    // outlet block
+    backQuad(poro1, out1, out2, poro2)
+    )
+
+    wall wall
+    (
+    // inlet block
+    quad2D(in1, join1)
+    quad2D(join2, in2)
+    // outlet block
+    quad2D(poro1, out1)
+    quad2D(out2, poro2)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    frontQuad(join1, poro1, poro2, join2)
+    // porosity block
+    backQuad(join1, poro1, poro2, join2)
+    // porosity block
+    quad2D(join1, poro1)
+    quad2D(poro2, join2)
+    )
+
+    patch inlet
+    (
+    quad2D(in2, in1)
+    )
+
+    patch outlet
+    (
+    quad2D(out2, out1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..948cb99e4312759b00458320d02423017267fd81
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    front
+    {
+        type            wall;
+        nFaces          700;
+        startFace       63400;
+    }
+    back
+    {
+        type            wall;
+        nFaces          700;
+        startFace       64100;
+    }
+    wall
+    {
+        type            wall;
+        nFaces          1400;
+        startFace       64800;
+    }
+    porosityWall
+    {
+        type            wall;
+        nFaces          1600;
+        startFace       66200;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       67800;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       68200;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
new file mode 100644
index 0000000000000000000000000000000000000000..634799837eaf7009528df07d027deafc7c6cc1f4
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/porousZones
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      porousZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    porosity
+    {
+        coordinateSystem
+        {
+            e1  (0.70710678 0.70710678 0);
+            e2  (0 0 1);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..7f3bedaae59eebf8d32dd99ac2f9c8ff8adf77a6
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/transportProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1.5e-05;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..88d62b98cc2648330cce6426b6043397caed194d
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPorousSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..a70c35e1bf81ab94f439453dfac9b9349f09bcae
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -0,0 +1,64 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(U)         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,U)      Gauss upwind;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,k)      Gauss upwind;
+}
+
+laplacianSchemes
+{
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian(rAU,p) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(1,p)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..f5f9feafa8b313dea68bfc6599515ae1f6b7b6ec
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration off;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    "(k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         2;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nUCorrectors    2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    p               0.3;
+    U               0.7;
+    k               0.9;
+    epsilon         0.9;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/constant/polyMesh/boundary b/tutorials/incompressible/simpleFoam/motorBike/constant/polyMesh/boundary
index b79c6ad432aaa5f1f3da3b5a4b9561d698e8cdb5..3c2a064fff04bf911ac962c8708db1b970fd7998 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/constant/polyMesh/boundary
+++ b/tutorials/incompressible/simpleFoam/motorBike/constant/polyMesh/boundary
@@ -21,433 +21,433 @@ FoamFile
     {
         type            patch;
         nFaces          320;
-        startFace       1016586;
+        startFace       1016673;
     }
     inlet
     {
         type            patch;
         nFaces          64;
-        startFace       1016906;
+        startFace       1016993;
     }
     outlet
     {
         type            patch;
         nFaces          64;
-        startFace       1016970;
+        startFace       1017057;
     }
     lowerWall
     {
         type            wall;
         nFaces          5330;
-        startFace       1017034;
+        startFace       1017121;
     }
     upperWall
     {
         type            patch;
         nFaces          160;
-        startFace       1022364;
+        startFace       1022451;
     }
     motorBike_frt-fairing:001%1
     {
         type            wall;
         nFaces          6626;
-        startFace       1022524;
+        startFace       1022611;
     }
     motorBike_windshield:002%2
     {
         type            wall;
         nFaces          50;
-        startFace       1029150;
+        startFace       1029237;
     }
     motorBike_rr-wh-rim:005%5
     {
         type            wall;
         nFaces          181;
-        startFace       1029200;
+        startFace       1029287;
     }
     motorBike_rr-wh-rim:010%10
     {
         type            wall;
         nFaces          340;
-        startFace       1029381;
+        startFace       1029468;
     }
     motorBike_fr-wh-rim:011%11
     {
         type            wall;
         nFaces          474;
-        startFace       1029721;
+        startFace       1029808;
     }
     motorBike_fr-wh-brake-disk:012%12
     {
         type            wall;
         nFaces          54;
-        startFace       1030195;
+        startFace       1030282;
     }
     motorBike_frame:016-shadow%13
     {
         type            wall;
         nFaces          131;
-        startFace       1030249;
+        startFace       1030336;
     }
     motorBike_rear-susp:014%14
     {
         type            wall;
         nFaces          1073;
-        startFace       1030380;
+        startFace       1030467;
     }
     motorBike_rear-susp:014-shadow%15
     {
         type            wall;
         nFaces          159;
-        startFace       1031453;
+        startFace       1031540;
     }
     motorBike_frame:016%16
     {
         type            wall;
         nFaces          20;
-        startFace       1031612;
+        startFace       1031699;
     }
     motorBike_rr-wh-rim:005-shadow%17
     {
         type            wall;
         nFaces          25;
-        startFace       1031632;
+        startFace       1031719;
     }
     motorBike_rr-wh-chain-hub:022%22
     {
         type            wall;
         nFaces          141;
-        startFace       1031657;
+        startFace       1031744;
     }
     motorBike_rearseat%24
     {
         type            wall;
         nFaces          432;
-        startFace       1031798;
+        startFace       1031885;
     }
     motorBike_frt-fairing%25
     {
         type            wall;
         nFaces          626;
-        startFace       1032230;
+        startFace       1032317;
     }
     motorBike_windshield%26
     {
         type            wall;
         nFaces          428;
-        startFace       1032856;
+        startFace       1032943;
     }
     motorBike_headlights%27
     {
         type            wall;
         nFaces          161;
-        startFace       1033284;
+        startFace       1033371;
     }
     motorBike_driversseat%28
     {
         type            wall;
         nFaces          367;
-        startFace       1033445;
+        startFace       1033532;
     }
     motorBike_rear-body%29
     {
         type            wall;
         nFaces          2076;
-        startFace       1033812;
+        startFace       1033899;
     }
     motorBike_fuel-tank%30
     {
         type            wall;
         nFaces          912;
-        startFace       1035888;
+        startFace       1035975;
     }
     motorBike_exhaust%31
     {
         type            wall;
         nFaces          2391;
-        startFace       1036800;
+        startFace       1036887;
     }
     motorBike_rr-wh-rim%32
     {
         type            wall;
         nFaces          1430;
-        startFace       1039191;
+        startFace       1039278;
     }
     motorBike_fr-mud-guard%33
     {
         type            wall;
         nFaces          767;
-        startFace       1040621;
+        startFace       1040708;
     }
     motorBike_fr-wh-rim%34
     {
         type            wall;
         nFaces          592;
-        startFace       1041388;
+        startFace       1041475;
     }
     motorBike_fr-wh-brake-disk%35
     {
         type            wall;
         nFaces          533;
-        startFace       1041980;
+        startFace       1042067;
     }
     motorBike_fr-brake-caliper%36
     {
         type            wall;
         nFaces          164;
-        startFace       1042513;
+        startFace       1042600;
     }
     motorBike_fr-wh-tyre%37
     {
         type            wall;
         nFaces          1118;
-        startFace       1042677;
+        startFace       1042764;
     }
     motorBike_hbars%38
     {
         type            wall;
         nFaces          535;
-        startFace       1043795;
+        startFace       1043882;
     }
     motorBike_fr-forks%39
     {
         type            wall;
         nFaces          1144;
-        startFace       1044330;
+        startFace       1044417;
     }
     motorBike_chain%40
     {
         type            wall;
         nFaces          474;
-        startFace       1045474;
+        startFace       1045561;
     }
     motorBike_rr-wh-tyre%41
     {
         type            wall;
         nFaces          1785;
-        startFace       1045948;
+        startFace       1046035;
     }
     motorBike_square-dial%42
     {
         type            wall;
         nFaces          6;
-        startFace       1047733;
+        startFace       1047820;
     }
     motorBike_round-dial%43
     {
         type            wall;
         nFaces          18;
-        startFace       1047739;
+        startFace       1047826;
     }
     motorBike_dial-holder%44
     {
         type            wall;
         nFaces          87;
-        startFace       1047757;
+        startFace       1047844;
     }
     motorBike_rear-susp%45
     {
         type            wall;
         nFaces          1787;
-        startFace       1047844;
+        startFace       1047931;
     }
     motorBike_rear-brake-lights%46
     {
         type            wall;
         nFaces          54;
-        startFace       1049631;
+        startFace       1049718;
     }
     motorBike_rear-light-bracket%47
     {
         type            wall;
         nFaces          163;
-        startFace       1049685;
+        startFace       1049772;
     }
     motorBike_frame%48
     {
         type            wall;
         nFaces          2040;
-        startFace       1049848;
+        startFace       1049935;
     }
     motorBike_rear-mud-guard%49
     {
         type            wall;
         nFaces          804;
-        startFace       1051888;
+        startFace       1051975;
     }
     motorBike_rear-susp-spring-damp%50
     {
         type            wall;
         nFaces          125;
-        startFace       1052692;
+        startFace       1052779;
     }
     motorBike_fairing-inner-plate%51
     {
         type            wall;
         nFaces          446;
-        startFace       1052817;
+        startFace       1052904;
     }
     motorBike_clutch-housing%52
     {
         type            wall;
         nFaces          966;
-        startFace       1053263;
+        startFace       1053350;
     }
     motorBike_radiator%53
     {
         type            wall;
         nFaces          48;
-        startFace       1054229;
+        startFace       1054316;
     }
     motorBike_water-pipe%54
     {
         type            wall;
         nFaces          103;
-        startFace       1054277;
+        startFace       1054364;
     }
     motorBike_water-pump%55
     {
         type            wall;
         nFaces          74;
-        startFace       1054380;
+        startFace       1054467;
     }
     motorBike_engine%56
     {
         type            wall;
         nFaces          2384;
-        startFace       1054454;
+        startFace       1054541;
     }
     motorBike_rear-shock-link%57
     {
         type            wall;
         nFaces          29;
-        startFace       1056838;
+        startFace       1056925;
     }
     motorBike_rear-brake-fluid-pot-bracket%58
     {
         type            wall;
         nFaces          59;
-        startFace       1056867;
+        startFace       1056954;
     }
     motorBike_rear-brake-fluid-pot%59
     {
         type            wall;
         nFaces          53;
-        startFace       1056926;
+        startFace       1057013;
     }
     motorBike_footpeg%60
     {
         type            wall;
         nFaces          87;
-        startFace       1056979;
+        startFace       1057066;
     }
     motorBike_rr-wh-chain-hub%61
     {
         type            wall;
         nFaces          145;
-        startFace       1057066;
+        startFace       1057153;
     }
     motorBike_rear-brake-caliper%62
     {
         type            wall;
         nFaces          142;
-        startFace       1057211;
+        startFace       1057298;
     }
     motorBike_rider-helmet%65
     {
         type            wall;
         nFaces          583;
-        startFace       1057353;
+        startFace       1057440;
     }
     motorBike_rider-visor%66
     {
         type            wall;
         nFaces          95;
-        startFace       1057936;
+        startFace       1058023;
     }
     motorBike_rider-boots%67
     {
         type            wall;
         nFaces          1025;
-        startFace       1058031;
+        startFace       1058118;
     }
     motorBike_rider-gloves%68
     {
         type            wall;
         nFaces          320;
-        startFace       1059056;
+        startFace       1059143;
     }
     motorBike_rider-body%69
     {
         type            wall;
         nFaces          4555;
-        startFace       1059376;
+        startFace       1059463;
     }
     motorBike_frame:0%70
     {
         type            wall;
         nFaces          37;
-        startFace       1063931;
+        startFace       1064018;
     }
     motorBike_frt-fairing:001-shadow%74
     {
         type            wall;
         nFaces          1274;
-        startFace       1063968;
+        startFace       1064055;
     }
     motorBike_windshield-shadow%75
     {
         type            wall;
         nFaces          101;
-        startFace       1065242;
+        startFace       1065329;
     }
     motorBike_fr-mud-guard-shadow%81
     {
         type            wall;
         nFaces          129;
-        startFace       1065343;
+        startFace       1065430;
     }
     motorBike_fr-wh-brake-disk-shadow%83
     {
         type            wall;
         nFaces          77;
-        startFace       1065472;
+        startFace       1065559;
     }
     motorBike_rear-mud-guard-shadow%84
     {
         type            wall;
         nFaces          138;
-        startFace       1065549;
+        startFace       1065636;
     }
     motorBike_rear-susp-spring-damp-shadow%85
     {
         type            wall;
         nFaces          15;
-        startFace       1065687;
+        startFace       1065774;
     }
     motorBike_radiator-shadow%86
     {
         type            wall;
         nFaces          12;
-        startFace       1065702;
+        startFace       1065789;
     }
     motorBike_rear-shock-link-shadow%87
     {
         type            wall;
         nFaces          7;
-        startFace       1065714;
+        startFace       1065801;
     }
     motorBike_rear-brake-fluid-pot-bracket-shadow%88
     {
         type            wall;
         nFaces          6;
-        startFace       1065721;
+        startFace       1065808;
     }
     motorBike_rr-wh-chain-hub-shadow%89
     {
         type            wall;
         nFaces          24;
-        startFace       1065727;
+        startFace       1065814;
     }
 )
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
index 50814119794de05fc809e924e6aac45e5589974c..72974933474f1f84d9bca41ef9f42a67048fd18b 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
@@ -24,13 +24,13 @@ gradSchemes
     default         Gauss linear;
     grad(p)         Gauss linear;
     grad(U)         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
+    //grad(U)         cellLimited Gauss linear 1;
 }
 
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss linearUpwindV Gauss linear;
+    div(phi,U)      Gauss linearUpwind grad(U);
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
     div((nuEff*dev(grad(U).T()))) Gauss linear;
@@ -39,8 +39,8 @@ divSchemes
 laplacianSchemes
 {
     default         Gauss linear corrected;
-//    default         Gauss linear limited 0.5;
-//    default         Gauss linear limited 0.333;
+    //default         Gauss linear limited 0.5;
+    //default         Gauss linear limited 0.333;
 }
 
 interpolationSchemes
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
index 3b0a3e35b679d7ba421255aaa4b43dfc003b8fe0..e99098b5b9c395c01a65884467579249c4e1dfa4 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
@@ -71,4 +71,9 @@ relaxationFactors
     omega           0.7;
 }
 
+cache
+{
+    grad(U);
+}
+
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/U b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/U
new file mode 100644
index 0000000000000000000000000000000000000000..2eee81b3a0d74c066ef53838efb6298e41d2bc8a
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/U
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+    floatingObject
+    {
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/alpha1 b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/alpha1
new file mode 100644
index 0000000000000000000000000000000000000000..8b3768458dcd903da0e526c0c8a44df4f5ab6329
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/alpha1
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha1;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    floatingObject
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/epsilon b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..849d3713a5f662ec9a6cb5aadcc61521fde2fab5
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.1;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.1;
+        value           uniform 0.1;
+    }
+    floatingObject
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/k b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/k
new file mode 100644
index 0000000000000000000000000000000000000000..84563be30c3703daa0ebae0b1e307bbeab826a75
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/k
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.1;
+        value           uniform 0.1;
+    }
+    floatingObject
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/nut b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/nut
new file mode 100644
index 0000000000000000000000000000000000000000..a86c3e35e008a370381b517395b721d5b80247f3
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/nut
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    atmosphere
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    floatingObject
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/p b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/p
new file mode 100644
index 0000000000000000000000000000000000000000..b5127205565f8a4d0e734bf5f989b385e51ff768
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/p
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      pd;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 0;
+        U               U;
+        phi             phi;
+        rho             rho;
+        psi             none;
+        gamma           1;
+        value           uniform 0;
+    }
+    floatingObject
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement
new file mode 100644
index 0000000000000000000000000000000000000000..711e09e22d868ecb68ee12141b909f673409fab7
--- /dev/null
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       pointVectorField;
+    location    "0.01";
+    object      pointDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    stationaryWalls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    floatingObject
+    {
+        type            sixDoFRigidBodyDisplacement;
+        centreOfMass    (0.5 0.5 0.5);
+        momentOfInertia (0.08622222 0.8622222 0.144);
+        mass            9.6;
+        rhoInf          1;
+        Q               (1 0 0 0 1 0 0 0 1);
+        v               (0 0 0);
+        a               (0 0 0);
+        pi              (0 0 0);
+        tau             (0 0 0);
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/Allrun b/tutorials/multiphase/interDyMFoam/ras/floatingObject/Allrun
index 1bece85bd98222c0358d809de9c9d3033aeda6c4..7c7ab52940262234486363281c2195f97e86b06e 100755
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/Allrun
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/Allrun
@@ -20,7 +20,7 @@ makeMeshByCellSet()
 runApplication blockMesh
 makeMeshByCellSet 1 2
 runApplication subsetMesh -overwrite c0 -patch floatingObject
-cp -r 0.orig 0 > /dev/null 2>&1
+cp -r 0.org 0 > /dev/null 2>&1
 runApplication setFields
 runApplication $application
 
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..26852141c2105c1e29071bcfc24a91c920937a02
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha1.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha1.org
new file mode 100644
index 0000000000000000000000000000000000000000..f0bdefa8cddb1ae29fe7f3c675625c288f9a509e
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha1.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha1;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+
+    rightWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org
new file mode 100644
index 0000000000000000000000000000000000000000..122e71dacedafedd83ea01b45380bc03d789950a
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha2.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+
+    rightWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha3.org b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha3.org
new file mode 100644
index 0000000000000000000000000000000000000000..f3d69ec08266a13c1c1f9ce2a80b2f1bdb6655cd
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/alpha3.org
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha3;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+
+    rightWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/p b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..b77bfe911838f761d5818a97593e1a4652949d22
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/0/p
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    rightWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    lowerWall
+    {
+        type            buoyantPressure;
+        value           uniform 0;
+    }
+
+    atmosphere
+    {
+        type            totalPressure;
+        p0              uniform 0;
+        U               U;
+        phi             phi;
+        rho             rho;
+        psi             none;
+        gamma           1;
+        value           uniform 0;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..cff4dbd23191f6a4ca343bf6898b0762ed0cce6d
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allclean
@@ -0,0 +1,4 @@
+#!/bin/sh
+
+foamCleanTutorials cases
+rm -rf 0/alpha[1-3].gz
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..15f6fafd90d19159e7c83d7192668c297d0e6798
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/Allrun
@@ -0,0 +1,10 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+cp 0/alpha1.org 0/alpha1
+cp 0/alpha2.org 0/alpha2
+cp 0/alpha3.org 0/alpha3
+runApplication setFields
+runApplication interMixingFoam
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..8fa634496563aaf871d8be56a7d7a1159758692c
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/dynamicMeshDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   staticFvMesh;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..27d4d324881b52b2a37d45b4cfbed46ed94d8051
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..a2553926d290a80076f1bf827324fd256d933b53
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/blockMeshDict
@@ -0,0 +1,92 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.146;
+
+vertices        
+(
+    (0 0 0)
+    (2 0 0)
+    (2.16438 0 0)
+    (4 0 0)
+    (0 0.32876 0)
+    (2 0.32876 0)
+    (2.16438 0.32876 0)
+    (4 0.32876 0)
+    (0 4 0)
+    (2 4 0)
+    (2.16438 4 0)
+    (4 4 0)
+    (0 0 0.1)
+    (2 0 0.1)
+    (2.16438 0 0.1)
+    (4 0 0.1)
+    (0 0.32876 0.1)
+    (2 0.32876 0.1)
+    (2.16438 0.32876 0.1)
+    (4 0.32876 0.1)
+    (0 4 0.1)
+    (2 4 0.1)
+    (2.16438 4 0.1)
+    (4 4 0.1)
+);
+
+blocks          
+(
+    hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1)
+    hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1)
+    hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1)
+    hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1)
+    hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1)
+);
+
+edges           
+(
+);
+
+patches         
+(
+    wall leftWall 
+    (
+        (0 12 16 4)
+        (4 16 20 8)
+    )
+    wall rightWall 
+    (
+        (7 19 15 3)
+        (11 23 19 7)
+    )
+    wall lowerWall 
+    (
+        (0 1 13 12)
+        (1 5 17 13)
+        (5 6 18 17)
+        (2 14 18 6)
+        (2 3 15 14)
+    )
+    patch atmosphere 
+    (
+        (8 20 21 9)
+        (9 21 22 10)
+        (10 22 23 11)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..cfd3f8c4d431b8e5a2ab2b5035b4d4f6747f6b36
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/polyMesh/boundary
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    leftWall
+    {
+        type            wall;
+        nFaces          50;
+        startFace       4432;
+    }
+    rightWall
+    {
+        type            wall;
+        nFaces          50;
+        startFace       4482;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          62;
+        startFace       4532;
+    }
+    atmosphere
+    {
+        type            patch;
+        nFaces          46;
+        startFace       4594;
+    }
+    defaultFaces
+    {
+        type            empty;
+        nFaces          4536;
+        startFace       4640;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..3cada9d4de2381a5996e8ef87393adffd7f10119
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/transportProperties
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Air
+phase1
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0] 1.48e-05;
+    rho             rho [1 -3 0 0 0 0 0] 1;
+}
+
+// Other Liquid
+phase2
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0]  1e-6;
+    rho             rho [1 -3 0 0 0 0 0] 1010;
+}
+
+// Water
+phase3
+{
+    transportModel  Newtonian;
+    nu              nu [0 2 -1 0 0 0 0]  1e-6;
+    rho             rho [1 -3 0 0 0 0 0] 1000;
+}
+
+// Surface tension coefficients
+sigma12           sigma12 [1 0 -2 0 0 0 0] 0.05;
+sigma13           sigma13 [1 0 -2 0 0 0 0] 0.04;
+
+// Diffusivity between miscible phases
+D23               D23   [0 2 -1 0 0 0 0]  3e-09;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..9cfc50a3d927937cf2ce96c64d3c9d6bd3144e2c
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/controlDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..a6f556508b1043f8de1653b692df0c37a7050d82
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.05;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..f1482468d421307958abf604cda46138a6feee91
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/decomposeParDict
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+metisCoeffs
+{
+    processorWeights ( 1 1 1 1 );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..bc9a0699584078169b5ede5cde5d4a94eac8c86d
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(rho*phi,U)  Gauss limitedLinearV 1;
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss interfaceCompression;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+    pcorr;
+    "alpha.";
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..f2543e14e1d63795f0d15e7deb675f27b1affb02
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSolution
@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha."
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+        nSweeps         1;
+    }
+
+    pcorr
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-10;
+        relTol          0;
+    }
+
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-07;
+        relTol          0.05;
+    }
+
+    pFinal
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-07;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    momentumPredictor no;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    nAlphaCorr      1;
+    nAlphaSubCycles 2;
+    cAlpha          1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..51ad3170f569e195dd9a8613a81f53657a277b86
--- /dev/null
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/setFieldsDict
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha1 0
+    volScalarFieldValue alpha2 1
+    volScalarFieldValue alpha3 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 0
+            volScalarFieldValue alpha2 0
+            volScalarFieldValue alpha3 1
+        );
+    }
+    boxToCell
+    {
+        box (0.1461 0.05 -1) (1 1 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 1
+            volScalarFieldValue alpha2 0
+            volScalarFieldValue alpha3 0
+        );
+    }
+);
+
+
+// ************************************************************************* //