diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/files b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..a798455eb7f807cacb1cdcea8d70c02d871821aa
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/files
@@ -0,0 +1,4 @@
+rhoPorousMRFPimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoPorousMRFPimpleFoam
+
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/options b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..61c1b6fe4694070e5a5a86789da6166158151879
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/options
@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I../rhoPimpleFoam \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3b4ead17e7cdfc69b84a104453844ae5daac3478
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
@@ -0,0 +1,39 @@
+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
+(
+    pZones.ddt(rho, U)
+  + fvm::div(phi, U)
+  + turbulence->divDevRhoReff(U)
+);
+
+if (oCorr == nOuterCorr-1)
+{
+    UEqn().relax(1);
+}
+else
+{
+    UEqn().relax();
+}
+
+mrfZones.addCoriolis(rho, UEqn());
+pZones.addResistance(UEqn());
+
+volScalarField rUA = 1.0/UEqn().A();
+
+if (momentumPredictor)
+{
+    if (oCorr == nOuterCorr-1)
+    {
+        solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
+    }
+    else
+    {
+        solve(UEqn() == -fvc::grad(p));
+    }
+}
+else
+{
+    U = rUA*(UEqn().H() - fvc::grad(p));
+    U.correctBoundaryConditions();
+}
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..b9a86ef9957b30923c2b3dd66606dcd8a33159e6
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H
@@ -0,0 +1,70 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar pMin
+    (
+        mesh.solutionDict().subDict("PIMPLE").lookup("pMin")
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    //dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..c74fb4d84b93c29385d419d6b858c234742e3a9d
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
@@ -0,0 +1,123 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn().A();
+U = rUA*UEqn().H();
+
+if (nCorr <= 1)
+{
+    UEqn.clear();
+}
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+    mrfZones.relativeFlux(fvc::interpolate(psi), phid);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+            && corr == nCorr-1
+            && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+        //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+//if (oCorr != nOuterCorr-1)
+{
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    rho = thermo.rho();
+    rho.relax();
+    Info<< "rho max/min : " << max(rho).value()
+        << " " << min(rho).value() << endl;
+}
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+bound(p, pMin);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+/*
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+*/
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..45d7646d042e1c3bdb6bfb24b955cb473ff1643b
--- /dev/null
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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
+    rhoPorousMRFPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    with support for porous media and MRF for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "bound.H"
+#include "MRFZones.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.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPIMPLEControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (nOuterCorr != 1)
+        {
+            p.storePrevIter();
+            rho.storePrevIter();
+        }
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
+            #include "UEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H
index 819cd0f538b6a87d0448618482caf75c52e3c793..f54526276762aeee1aacd3e87fda796a961ce529 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H
@@ -59,7 +59,17 @@
                 alpharScheme
             );
 
-        MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+        MULES::explicitSolve
+        (
+            geometricOneField(),
+            alpha1,
+            phi,
+            phiAlpha1,
+            Sp,
+            Su,
+            1,
+            0
+        );
 
         surfaceScalarField rho1f = fvc::interpolate(rho1);
         surfaceScalarField rho2f = fvc::interpolate(rho2);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
index 819cd0f538b6a87d0448618482caf75c52e3c793..dd704c0693c45b5cd7ded0bd48b378cd42802210 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H
@@ -59,7 +59,7 @@
                 alpharScheme
             );
 
-        MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
 
         surfaceScalarField rho1f = fvc::interpolate(rho1);
         surfaceScalarField rho2f = fvc::interpolate(rho2);
diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
index f9bad0a7058577487b0cb6e4fd7bb14f1e9ee24d..972d1d71218fba615638aedc827c701a7496e854 100644
--- a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
+++ b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H
@@ -70,12 +70,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha1,
             phiAlpha1BD,
             phiAlpha1,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3
@@ -107,12 +107,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha2,
             phiAlpha2BD,
             phiAlpha2,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 15a5291ee63596e75b32ed0b4459b7390d601e88..a232056b8dcead3b6b13a2868ec36415e150cbfd 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -51,8 +51,8 @@
         );
 
         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
-        //MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
-        MULES::implicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
 
         rhoPhi +=
             (runTime.deltaT()/totalDeltaT)
diff --git a/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H
new file mode 100644
index 0000000000000000000000000000000000000000..3f9e2865a7f15098caf76084205d5d030dc005f8
--- /dev/null
+++ b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::oneFieldField
+
+Description
+    A class representing the concept of a field of oneFields used to
+    avoid unnecessary manipulations for objects which are known to be one at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oneFieldField_H
+#define oneFieldField_H
+
+#include "oneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class oneField Declaration
+\*---------------------------------------------------------------------------*/
+
+class oneFieldField
+:
+    public one
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        oneFieldField()
+        {}
+
+
+    // Member Operators
+
+        inline oneField operator[](const label) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "oneFieldFieldI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldFieldI.H b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldFieldI.H
new file mode 100644
index 0000000000000000000000000000000000000000..1fe0c66c4c2c04374ba7eeaade6063d1d00a5be8
--- /dev/null
+++ b/src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldFieldI.H
@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "oneFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::oneField Foam::oneFieldField::operator[](const label) const
+{
+    return oneField();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H
new file mode 100644
index 0000000000000000000000000000000000000000..c551d39bc2b6aa0d58c7ac6f3e383a213d23058d
--- /dev/null
+++ b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::zeroFieldField
+
+Description
+    A class representing the concept of a field of zeroFields used to
+    avoid unnecessary manipulations for objects which are known to be zero at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef zeroFieldField_H
+#define zeroFieldField_H
+
+#include "zeroField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class zeroField Declaration
+\*---------------------------------------------------------------------------*/
+
+class zeroFieldField
+:
+    public zero
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        zeroFieldField()
+        {}
+
+
+    // Member Operators
+
+        inline zeroField operator[](const label) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "zeroFieldFieldI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldFieldI.H b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldFieldI.H
new file mode 100644
index 0000000000000000000000000000000000000000..2eecfae7ad871970a9afbdf49206f38ece80c15c
--- /dev/null
+++ b/src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldFieldI.H
@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "zeroFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::zeroField Foam::zeroFieldField::operator[](const label) const
+{
+    return zeroField();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/oneField/oneField.H b/src/OpenFOAM/fields/Fields/oneField/oneField.H
index 8a21131b267bfa7e440bc30af55446d0b7b13add..c4869e5a16d2dac891931dbfa5df93a548d4c124 100644
--- a/src/OpenFOAM/fields/Fields/oneField/oneField.H
+++ b/src/OpenFOAM/fields/Fields/oneField/oneField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,7 @@ class oneField
 :
     public one
 {
+
 public:
 
     // Constructors
@@ -65,10 +66,6 @@ public:
     // Member Operators
 
         inline scalar operator[](const label) const;
-
-        inline oneField field() const;
-
-        inline oneField oldTime() const;
 };
 
 
diff --git a/src/OpenFOAM/fields/Fields/oneField/oneFieldI.H b/src/OpenFOAM/fields/Fields/oneField/oneFieldI.H
index 7809caa045d93e7343f56f39312fe80f03e57ca1..1c0526f37076bd7499408fa5f80581af864e3ba2 100644
--- a/src/OpenFOAM/fields/Fields/oneField/oneFieldI.H
+++ b/src/OpenFOAM/fields/Fields/oneField/oneFieldI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,15 +33,5 @@ inline Foam::scalar Foam::oneField::operator[](const label) const
     return scalar(1);
 }
 
-inline Foam::oneField Foam::oneField::field() const
-{
-    return oneField();
-}
-
-inline Foam::oneField Foam::oneField::oldTime() const
-{
-    return oneField();
-}
-
 
 // ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
index 97f090623ccbf09d97c3b45de5bfc6f7510a662d..0281c154b6a92258eebc907f3df7ab0837b97733 100644
--- a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
+++ b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,10 +65,6 @@ public:
     // Member Operators
 
         inline scalar operator[](const label) const;
-
-        inline zeroField field() const;
-
-        inline zeroField oldTime() const;
 };
 
 
@@ -78,7 +74,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#   include "zeroFieldI.H"
+#include "zeroFieldI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H b/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
index 67d1dd0e9e0475af2938e2c082e979dd6872eaed..38072ac2e207e76b61c1b9857f9689b112fd0325 100644
--- a/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
+++ b/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
@@ -33,14 +33,4 @@ inline Foam::scalar Foam::zeroField::operator[](const label) const
     return scalar(0);
 }
 
-inline Foam::zeroField Foam::zeroField::field() const
-{
-    return zeroField();
-}
-
-inline Foam::zeroField Foam::zeroField::oldTime() const
-{
-    return zeroField();
-}
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneField.H b/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneField.H
new file mode 100644
index 0000000000000000000000000000000000000000..51aa859b10710a619d690fcf51106ed6da7e1bcf
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneField.H
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::GeometricField
+
+Description
+    A class representing the concept of a GeometricField of 1 used to avoid
+    unnecessary manipulations for objects which are known to be one at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef geometricOneField_H
+#define geometricOneField_H
+
+#include "oneFieldField.H"
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class geometricOneField Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometricOneField
+:
+    public one
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        geometricOneField()
+        {}
+
+
+    // Member Operators
+
+        inline scalar operator[](const label) const;
+
+        inline oneField field() const;
+
+        inline oneField oldTime() const;
+
+        inline oneFieldField boundaryField() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "geometricOneFieldI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneFieldI.H b/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneFieldI.H
new file mode 100644
index 0000000000000000000000000000000000000000..cd9a79e2c9f4b765a8e9988728490b12e01ae0b9
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneFieldI.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "geometricOneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::geometricOneField::operator[](const label) const
+{
+    return scalar(1);
+}
+
+inline Foam::oneField Foam::geometricOneField::field() const
+{
+    return oneField();
+}
+
+inline Foam::oneField Foam::geometricOneField::oldTime() const
+{
+    return oneField();
+}
+
+inline Foam::oneFieldField Foam::geometricOneField::boundaryField() const
+{
+    return oneFieldField();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H
new file mode 100644
index 0000000000000000000000000000000000000000..1168d936ee0124da7ee8d477108a86e3089ca92e
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::GeometricField
+
+Description
+    A class representing the concept of a GeometricField of 1 used to avoid
+    unnecessary manipulations for objects which are known to be zero at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef geometricZeroField_H
+#define geometricZeroField_H
+
+#include "zeroFieldField.H"
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class geometricZeroField Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometricZeroField
+:
+    public zero
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        geometricZeroField()
+        {}
+
+
+    // Member Operators
+
+        inline scalar operator[](const label) const;
+
+        inline zeroField field() const;
+
+        inline zeroField oldTime() const;
+
+        inline zeroFieldField boundaryField() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "geometricZeroFieldI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H
new file mode 100644
index 0000000000000000000000000000000000000000..027145167e9da111c8ca7f187f995b10715556ce
--- /dev/null
+++ b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "geometricZeroField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::geometricZeroField::operator[](const label) const
+{
+    return scalar(0);
+}
+
+inline Foam::zeroField Foam::geometricZeroField::field() const
+{
+    return zeroField();
+}
+
+inline Foam::zeroField Foam::geometricZeroField::oldTime() const
+{
+    return zeroField();
+}
+
+inline Foam::zeroFieldField Foam::geometricZeroField::boundaryField() const
+{
+    return zeroFieldField();
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 3618d87ef5dfa6ea4ea03003cf5928d4d3c3c2e0..3ed30360dce138385fbe5e27c0195950e2ac4800 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ License
 #include "fvMatrices.H"
 #include "syncTools.H"
 #include "faceSet.H"
+#include "geometricOneField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -377,85 +378,33 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
 
 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 {
-    const surfaceVectorField& Cf = mesh_.Cf();
-    const surfaceVectorField& Sf = mesh_.Sf();
-
-    const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
-
-    // Internal faces
-    forAll(internalFaces_, i)
-    {
-        label facei = internalFaces_[i];
-        phi[facei] -= (Omega ^ (Cf[facei] - origin)) & Sf[facei];
-    }
-
-    // Included patches
-    forAll(includedFaces_, patchi)
-    {
-        forAll(includedFaces_[patchi], i)
-        {
-            label patchFacei = includedFaces_[patchi][i];
-
-            phi.boundaryField()[patchi][patchFacei] = 0.0;
-        }
-    }
+    relativeRhoFlux(geometricOneField(), phi);
+}
 
-    // Excluded patches
-    forAll(excludedFaces_, patchi)
-    {
-        forAll(excludedFaces_[patchi], i)
-        {
-            label patchFacei = excludedFaces_[patchi][i];
 
-            phi.boundaryField()[patchi][patchFacei] -=
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
-        }
-    }
+void Foam::MRFZone::relativeFlux
+(
+    const surfaceScalarField& rho,
+    surfaceScalarField& phi
+) const
+{
+    relativeRhoFlux(rho, phi);
 }
 
 
 void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
 {
-    const surfaceVectorField& Cf = mesh_.Cf();
-    const surfaceVectorField& Sf = mesh_.Sf();
-
-    const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
-
-    // Internal faces
-    forAll(internalFaces_, i)
-    {
-        label facei = internalFaces_[i];
-        phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
-    }
-
-    // Included patches
-    forAll(includedFaces_, patchi)
-    {
-        forAll(includedFaces_[patchi], i)
-        {
-            label patchFacei = includedFaces_[patchi][i];
-
-            phi.boundaryField()[patchi][patchFacei] +=
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
-        }
-    }
+    absoluteRhoFlux(geometricOneField(), phi);
+}
 
-    // Excluded patches
-    forAll(excludedFaces_, patchi)
-    {
-        forAll(excludedFaces_[patchi], i)
-        {
-            label patchFacei = excludedFaces_[patchi][i];
 
-            phi.boundaryField()[patchi][patchFacei] +=
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
-              & Sf.boundaryField()[patchi][patchFacei];
-        }
-    }
+void Foam::MRFZone::absoluteFlux
+(
+    const surfaceScalarField& rho,
+    surfaceScalarField& phi
+) const
+{
+    absoluteRhoFlux(rho, phi);
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index d5554a5b0810ea1fab45d5474df22397db3390d7..f36a8350e5055dab75f40fb572e7830ccebd5ac8 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,6 +97,22 @@ class MRFZone
         //- Divide faces in frame according to patch
         void setMRFFaces();
 
+        //- Make the given absolute mass/vol flux relative within the MRF region
+        template<class RhoFieldType>
+        void relativeRhoFlux
+        (
+            const RhoFieldType& rho,
+            surfaceScalarField& phi
+        ) const;
+
+        //- Make the given relative mass/vol flux absolute within the MRF region
+        template<class RhoFieldType>
+        void absoluteRhoFlux
+        (
+            const RhoFieldType& rho,
+            surfaceScalarField& phi
+        ) const;
+
         //- Disallow default bitwise copy construct
         MRFZone(const MRFZone&);
 
@@ -163,9 +179,23 @@ public:
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given absolute mass-flux relative within the MRF region
+        void relativeFlux
+        (
+            const surfaceScalarField& rho,
+            surfaceScalarField& phi
+        ) const;
+
         //- Make the given relative flux absolute within the MRF region
         void absoluteFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative mass-flux absolute within the MRF region
+        void absoluteFlux
+        (
+            const surfaceScalarField& rho,
+            surfaceScalarField& phi
+        ) const;
+
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 
@@ -186,6 +216,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "MRFZoneTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..b195f111a4457b78991afe1a2144871e4007b04c
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "MRFZone.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "geometricOneField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::MRFZone::relativeRhoFlux
+(
+    const RhoFieldType& rho,
+    surfaceScalarField& phi
+) const
+{
+    const surfaceVectorField& Cf = mesh_.Cf();
+    const surfaceVectorField& Sf = mesh_.Sf();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    // Internal faces
+    forAll(internalFaces_, i)
+    {
+        label facei = internalFaces_[i];
+        phi[facei] -= rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei];
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] = 0.0;
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] -=
+                rho.boundaryField()[patchi][patchFacei]
+               *(Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::MRFZone::absoluteRhoFlux
+(
+    const RhoFieldType& rho,
+    surfaceScalarField& phi
+) const
+{
+    const surfaceVectorField& Cf = mesh_.Cf();
+    const surfaceVectorField& Sf = mesh_.Sf();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    // Internal faces
+    forAll(internalFaces_, i)
+    {
+        label facei = internalFaces_[i];
+        phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
index 5ef8b3ee83ffbac7bdb843f8873d842aedd5208d..6bfafcca779e232275037191cabcb6b2be74dcec 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,19 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZones::relativeFlux
+(
+    const surfaceScalarField& rho,
+    surfaceScalarField& phi
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).relativeFlux(rho, phi);
+    }
+}
+
+
 void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
 {
     forAll(*this, i)
@@ -105,6 +118,19 @@ void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZones::absoluteFlux
+(
+    const surfaceScalarField& rho,
+    surfaceScalarField& phi
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).absoluteFlux(rho, phi);
+    }
+}
+
+
 void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
 {
     forAll(*this, i)
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
index 6e979903a9b803fa4b897863876b3f23765517ff..b074cd5fd54726a622d1be376139ab9a2ae189a4 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,9 +85,23 @@ public:
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given absolute mass-flux relative within the MRF region
+        void relativeFlux
+        (
+            const surfaceScalarField& rho,
+            surfaceScalarField& phi
+        ) const;
+
         //- Make the given relative flux absolute within the MRF region
         void absoluteFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative mass-flux absolute within the MRF region
+        void absoluteFlux
+        (
+            const surfaceScalarField& rho,
+            surfaceScalarField& phi
+        ) const;
+
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 };
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
index 444cd2513107f330b2d5a071462f9e3571dda3ed..92a5a822a057a90b4b10b6691a5089c942ee774e 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
+++ b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
@@ -87,7 +87,7 @@ template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> >
 Foam::PorousZones<ZoneType>::ddt
 (
-    const oneField&,
+    const geometricOneField&,
     GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
index cc55fc42646b30589cf4f32da1b40d16c3461ca4..b86bf6c573959b5f9a10e937559101c72643b35f 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
+++ b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
@@ -44,7 +44,7 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "fvMatricesFwd.H"
 #include "dimensionedScalarFwd.H"
-#include "oneField.H"
+#include "geometricOneField.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -105,7 +105,7 @@ public:
         template<class Type>
         tmp<fvMatrix<Type> > ddt
         (
-            const oneField&,
+            const geometricOneField&,
             GeometricField<Type, fvPatchField, volMesh>&
         );
 
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C
index 3251e2541eeeb5834d134dd62622477dc9c820fe..1e334d0dfc696d975d102d5d67e3bdd80b6d06fe 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C
+++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C
@@ -27,7 +27,7 @@ License
 #include "porousZone.H"
 #include "fvMesh.H"
 #include "fvMatrices.H"
-#include "oneField.H"
+#include "geometricOneField.H"
 
 // * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
 
@@ -236,7 +236,7 @@ void Foam::porousZone::addResistance(fvVectorMatrix& UEqn) const
                 Udiag,
                 cells,
                 V,
-                oneField(),
+                geometricOneField(),
                 U
             );
         }
@@ -268,7 +268,7 @@ void Foam::porousZone::addResistance(fvVectorMatrix& UEqn) const
                 Usource,
                 cells,
                 V,
-                oneField(),
+                geometricOneField(),
                 mesh_.lookupObject<volScalarField>("nu"),
                 U
             );
@@ -316,7 +316,7 @@ void Foam::porousZone::addResistance
             (
                 AU,
                 cells,
-                oneField(),
+                geometricOneField(),
                 U
             );
         }
@@ -344,7 +344,7 @@ void Foam::porousZone::addResistance
             (
                 AU,
                 cells,
-                oneField(),
+                geometricOneField(),
                 mesh_.lookupObject<volScalarField>("nu"),
                 U
             );
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
index d7610206d00b9c2781f22cc855c4750fdc23f297..b3ae5b8f37d9e4a264c85ffd9b269fb19e8c3eed 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
@@ -60,7 +60,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 ddt
 (
-    const oneField&,
+    const one&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
index 2b8124e5596a2562a4d8dded29e90ac3ab78655a..31e0e18fab6ffdd8deb52cc461243e20b74e9216 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
@@ -38,7 +38,7 @@ SourceFiles
 
 #include "volFieldsFwd.H"
 #include "fvMatrix.H"
-#include "oneField.H"
+#include "one.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,7 +60,7 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > ddt
     (
-        const oneField&,
+        const one&,
         const GeometricField<Type, fvPatchField, volMesh>&
     );
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
index 501456da3cb3cb807d90eeb683f0b5d4e4cddcba..24128f4c0d8fbc31af6b7bb23c850fd952da0c14 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
@@ -99,7 +99,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    const zeroField&,
+    const zero&,
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
@@ -115,7 +115,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    const zeroField&,
+    const zero&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
@@ -130,7 +130,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    const oneField&,
+    const one&,
     const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
@@ -143,7 +143,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    const oneField&,
+    const one&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
index c41774361f821b592f2f92bf88b1dfe73f2ab219..4025505deb3541991e0f0a0cc5d7493156199991 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
@@ -39,7 +39,8 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "surfaceFieldsFwd.H"
 #include "fvMatrix.H"
-#include "oneField.H"
+#include "zero.H"
+#include "one.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,7 +70,7 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        const zeroField&,
+        const zero&,
         const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
@@ -77,7 +78,7 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        const zeroField&,
+        const zero&,
         const GeometricField<Type, fvPatchField, volMesh>&
     );
 
@@ -85,7 +86,7 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        const oneField&,
+        const one&,
         const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
@@ -93,7 +94,7 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        const oneField&,
+        const one&,
         const GeometricField<Type, fvPatchField, volMesh>&
     );
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.C b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
index 99e2d147467d4bdb290e704d423a9d069071e1b7..28d592c5716e8c9c54f940737a1a33828981034b 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
@@ -85,7 +85,7 @@ template<class Type>
 Foam::zeroField
 Foam::fvm::Su
 (
-    const zeroField&,
+    const zero&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
@@ -174,7 +174,7 @@ template<class Type>
 Foam::zeroField
 Foam::fvm::Sp
 (
-    const zeroField&,
+    const zero&,
     const GeometricField<Type, fvPatchField, volMesh>&
 )
 {
@@ -240,7 +240,7 @@ template<class Type>
 Foam::zeroField
 Foam::fvm::SuSp
 (
-    const zeroField&,
+    const zero&,
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.H b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
index a5662c1733beb09122ed2349d1c0bf07b3e45c74..56d3aa09e9aedb5f19dccc06c7b409fc7f1f7387 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "volFieldsFwd.H"
 #include "fvMatrix.H"
+#include "zeroField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,7 +77,7 @@ namespace fvm
         template<class Type>
         zeroField Su
         (
-            const zeroField&,
+            const zero&,
             const GeometricField<Type, fvPatchField, volMesh>&
         );
 
@@ -116,7 +117,7 @@ namespace fvm
         template<class Type>
         zeroField Sp
         (
-            const zeroField&,
+            const zero&,
             const GeometricField<Type, fvPatchField, volMesh>&
         );
 
@@ -147,7 +148,7 @@ namespace fvm
         template<class Type>
         zeroField SuSp
         (
-            const zeroField&,
+            const zero&,
             const GeometricField<Type, fvPatchField, volMesh>&
         );
 }
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 9ce8d888c4bfdedad90d8216809664424428d19f..f335eeab300fd1c1387e056511beadd89a3b931a 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1088,7 +1088,7 @@ void Foam::fvMatrix<Type>::operator-=
 template<class Type>
 void Foam::fvMatrix<Type>::operator+=
 (
-    const zeroField&
+    const zero&
 )
 {}
 
@@ -1096,7 +1096,7 @@ void Foam::fvMatrix<Type>::operator+=
 template<class Type>
 void Foam::fvMatrix<Type>::operator-=
 (
-    const zeroField&
+    const zero&
 )
 {}
 
@@ -1507,7 +1507,7 @@ template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
 (
     const fvMatrix<Type>& A,
-    const zeroField&
+    const zero&
 )
 {
     return A;
@@ -1518,7 +1518,7 @@ template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
 (
     const tmp<fvMatrix<Type> >& tA,
-    const zeroField&
+    const zero&
 )
 {
     return tA;
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index 8569862e13b7b5bafa8b32e6e49dd633534773f9..fc3efa62611728a8eebdbafc0591be38c22882af 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -44,7 +44,7 @@ SourceFiles
 #include "tmp.H"
 #include "autoPtr.H"
 #include "dimensionedTypes.H"
-#include "zeroField.H"
+#include "zero.H"
 #include "className.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -435,8 +435,8 @@ public:
         void operator+=(const dimensioned<Type>&);
         void operator-=(const dimensioned<Type>&);
 
-        void operator+=(const zeroField&);
-        void operator-=(const zeroField&);
+        void operator+=(const zero&);
+        void operator-=(const zero&);
 
         void operator*=(const DimensionedField<scalar, volMesh>&);
         void operator*=(const tmp<DimensionedField<scalar, volMesh> >&);
@@ -647,14 +647,14 @@ template<class Type>
 tmp<fvMatrix<Type> > operator==
 (
     const fvMatrix<Type>&,
-    const zeroField&
+    const zero&
 );
 
 template<class Type>
 tmp<fvMatrix<Type> > operator==
 (
     const tmp<fvMatrix<Type> >&,
-    const zeroField&
+    const zero&
 );
 
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index 9c055ca339abec6f11138a642696ebe529ea249a..dd7c35ea6e86ef65469cdcdd7a748c8780ca0745 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ License
 #include "fvcSurfaceIntegrate.H"
 #include "slicedSurfaceFields.H"
 #include "syncTools.H"
+
 #include "fvm.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,11 +50,11 @@ void Foam::MULES::explicitSolve
 {
     explicitSolve
     (
-        oneField(),
+        geometricOneField(),
         psi,
         phi,
         phiPsi,
-        zeroField(), zeroField(),
+        zero(), zero(),
         psiMax, psiMin
     );
 }
@@ -70,11 +71,11 @@ void Foam::MULES::implicitSolve
 {
     implicitSolve
     (
-        oneField(),
+        geometricOneField(),
         psi,
         phi,
         phiPsi,
-        zeroField(), zeroField(),
+        zero(), zero(),
         psiMax, psiMin
     );
 }
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index 59c4072feea1139b4aab9b3439953493b3cd0abc..8db5bf347aa41acbbf7f817fa19531c374f926e5 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,8 +46,8 @@ SourceFiles
 #include "volFields.H"
 #include "surfaceFieldsFwd.H"
 #include "primitiveFieldsFwd.H"
-#include "zeroField.H"
-#include "oneField.H"
+#include "zero.H"
+#include "geometricOneField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index bf72983e4a7e069d328aacf27bb30a670185c89c..858ff4b15022f79d04d2b4256c1581a2199179b7 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -88,8 +88,8 @@ void Foam::MULES::explicitSolve
         psi,
         phiBD,
         phiCorr,
-        Sp.field(),
-        Su.field(),
+        Sp,
+        Su,
         psiMax,
         psiMin,
         3
@@ -109,18 +109,18 @@ void Foam::MULES::explicitSolve
         psiIf =
         (
             mesh.Vsc0()*rho.oldTime()*psi0/(deltaT*mesh.Vsc())
-          + Su.field()
+          + Su
           - psiIf
-        )/(rho/deltaT - Sp.field());
+        )/(rho/deltaT - Sp);
     }
     else
     {
         psiIf =
         (
             rho.oldTime()*psi0/deltaT
-          + Su.field()
+          + Su
           - psiIf
-        )/(rho/deltaT - Sp.field());
+        )/(rho/deltaT - Sp);
     }
 
     psi.correctBoundaryConditions();
@@ -244,8 +244,8 @@ void Foam::MULES::implicitSolve
             psi,
             phiBD,
             phiCorr,
-            Sp.field(),
-            Su.field(),
+            Sp,
+            Su,
             psiMax,
             psiMin,
             nLimiterIter
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/T b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..c8760070627a095b7b08401dea7b53162fc5efd5
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/T
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/U b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..ee6588db42363bdd3e5bb01d5ffbe353e241ba0a
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/U
@@ -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       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    rotor
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    stator
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/alphat b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/alphat
new file mode 100644
index 0000000000000000000000000000000000000000..764b2c5ca06b74cd6e8427e914c3d575f05777cc
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/alphat
@@ -0,0 +1,47 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    rotor
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    stator
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    front
+    {
+        type            empty;
+    }
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/epsilon b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..4443c71b68fe074f02ba87cb4b86d3e4383dd158
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/epsilon
@@ -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      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 20;
+
+boundaryField
+{
+    rotor
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 20;
+    }
+
+    stator
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 20;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/k b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..f5ee44315fa5ddfdc436a3a9f4f12ee1fa0ee061
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/k
@@ -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      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    rotor
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0;
+    }
+
+    stator
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/mut b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/mut
new file mode 100644
index 0000000000000000000000000000000000000000..4bdc2e1afdce630c4db0d7933909854f42bf2aa6
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/mut
@@ -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      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    rotor
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+
+    stator
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/p b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..de53bf93f6d62bd02269a9004e782ccd70c0ed3c
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/0/p
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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 1e5;
+
+boundaryField
+{
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/Allrun b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..836276fe9cea7a027c4726fa021088bf33d0da6e
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/Allrun
@@ -0,0 +1,9 @@
+#!/bin/sh
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application=`getApplication`
+
+./makeMesh
+runApplication $application
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFZones b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFZones
new file mode 100644
index 0000000000000000000000000000000000000000..bff11fb77f1dbc6c68f903566b10f52ae6b6648d
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFZones
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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      MRFZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    rotor
+    {
+        // Fixed patches (by default they 'move' with the MRF zone)
+        nonRotatingPatches ();
+
+        origin    origin [0 1 0 0 0 0 0]  (0 0 0);
+        axis      axis   [0 0 0 0 0 0 0]  (0 0 1);
+        omega     omega  [0 0 -1 0 0 0 0] 1047.2;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/RASProperties b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..81b1ec9115226882e6e980b8e5f4e19cabb886d9
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/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      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/dynamicMeshDict b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..10468d3ba0fa641033d09d7e57a76f12b0667164
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/dynamicMeshDict
@@ -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      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMeshLib "libtopoChangerFvMesh.so";
+
+dynamicFvMesh   mixerFvMesh;
+
+mixerFvMeshCoeffs
+{
+    coordinateSystem
+    {
+        type            cylindrical;
+        origin          ( 0 0 0 );
+        axis            ( 0 0 1 );
+        direction       ( 1 0 0 );
+    }
+
+    rpm             10;
+
+    slider
+    {
+        inside          insideSlider;
+        outside         outsideSlider;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..ef73b0c7b75531001c82f6b5b19cefa98298e5d4
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict
@@ -0,0 +1,834 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// General macros to create 2D/extruded-2D meshes
+
+
+
+
+
+
+
+
+
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+// Hub radius
+
+
+// Impeller-tip radius
+
+
+// Baffle-tip radius
+
+
+// Tank radius
+
+
+// MRF region radius
+
+
+// Thickness of 2D slab
+
+
+// Base z
+
+
+// Top z
+
+
+// Number of cells radially between hub and impeller tip
+
+
+// Number of cells radially in each of the two regions between
+// impeller and baffle tips
+
+
+// Number of cells radially between baffle tip and tank
+
+
+// Number of cells azimuthally in each of the 8 blocks
+
+
+// Number of cells in the thickness of the slab
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+vertices
+(
+    (0.2 0 0) // Vertex r0b = 0 
+    (0.2 0 0) // Vertex r0sb = 1 
+    (0.141421356364228 -0.141421356110391 0) // Vertex r1b = 2 
+    (3.58979347393082e-10 -0.2 0) // Vertex r2b = 3 
+    (3.58979347393082e-10 -0.2 0) // Vertex r2sb = 4 
+    (-0.141421355856554 -0.141421356618065 0) // Vertex r3b = 5 
+    (-0.2 7.17958694786164e-10 0) // Vertex r4b = 6 
+    (-0.2 7.17958694786164e-10 0) // Vertex r4sb = 7 
+    (-0.141421355856554 0.141421356618065 0) // Vertex r5b = 8 
+    (3.58979347393082e-10 0.2 0) // Vertex r6b = 9 
+    (3.58979347393082e-10 0.2 0) // Vertex r6sb = 10 
+    (0.141421356364228 0.141421356110391 0) // Vertex r7b = 11 
+
+    (0.5 0 0) // Vertex rb0b = 12 
+    (0.353553390910569 -0.353553390275978 0) // Vertex rb1b = 13 
+    (8.97448368482705e-10 -0.5 0) // Vertex rb2b = 14 
+    (-0.353553389641386 -0.353553391545162 0) // Vertex rb3b = 15 
+    (-0.5 1.79489673696541e-09 0) // Vertex rb4b = 16 
+    (-0.353553389641386 0.353553391545162 0) // Vertex rb5b = 17 
+    (8.97448368482705e-10 0.5 0) // Vertex rb6b = 18 
+    (0.353553390910569 0.353553390275978 0) // Vertex rb7b = 19 
+
+    (0.6 0 0) // Vertex ri0b = 20 
+    (0.424264069092683 -0.424264068331174 0) // Vertex ri1b = 21 
+    (1.07693804217925e-09 -0.6 0) // Vertex ri2b = 22 
+    (-0.424264067569663 -0.424264069854194 0) // Vertex ri3b = 23 
+    (-0.6 2.15387608435849e-09 0) // Vertex ri4b = 24 
+    (-0.424264067569663 0.424264069854194 0) // Vertex ri5b = 25 
+    (1.07693804217925e-09 0.6 0) // Vertex ri6b = 26 
+    (0.424264069092683 0.424264068331174 0) // Vertex ri7b = 27 
+
+    (0.7 0 0) // Vertex Rb0b = 28 
+    (0.494974747274797 -0.494974746386369 0) // Vertex Rb1b = 29 
+    (1.25642771587579e-09 -0.7 0) // Vertex Rb2b = 30 
+    (-0.49497474549794 -0.494974748163226 0) // Vertex Rb3b = 31 
+    (-0.7 2.51285543175157e-09 0) // Vertex Rb4b = 32 
+    (-0.49497474549794 0.494974748163226 0) // Vertex Rb5b = 33 
+    (1.25642771587579e-09 0.7 0) // Vertex Rb6b = 34 
+    (0.494974747274797 0.494974746386369 0) // Vertex Rb7b = 35 
+
+    (1 0 0) // Vertex R0b = 36 
+    (0.707106781821139 -0.707106780551956 0) // Vertex R1b = 37 
+    (0.707106781821139 -0.707106780551956 0) // Vertex R1sb = 38 
+    (1.79489673696541e-09 -1 0) // Vertex R2b = 39 
+    (-0.707106779282772 -0.707106783090323 0) // Vertex R3b = 40 
+    (-0.707106779282772 -0.707106783090323 0) // Vertex R3sb = 41 
+    (-1 3.58979347393082e-09 0) // Vertex R4b = 42 
+    (-0.707106779282772 0.707106783090323 0) // Vertex R5b = 43 
+    (-0.707106779282772 0.707106783090323 0) // Vertex R5sb = 44 
+    (1.79489673696541e-09 1 0) // Vertex R6b = 45 
+    (0.707106781821139 0.707106780551956 0) // Vertex R7b = 46 
+    (0.707106781821139 0.707106780551956 0) // Vertex R7sb = 47 
+
+    (0.2 0 0.1) // Vertex r0t = 48 
+    (0.2 0 0.1) // Vertex r0st = 49 
+    (0.141421356364228 -0.141421356110391 0.1) // Vertex r1t = 50 
+    (3.58979347393082e-10 -0.2 0.1) // Vertex r2t = 51 
+    (3.58979347393082e-10 -0.2 0.1) // Vertex r2st = 52 
+    (-0.141421355856554 -0.141421356618065 0.1) // Vertex r3t = 53 
+    (-0.2 7.17958694786164e-10 0.1) // Vertex r4t = 54 
+    (-0.2 7.17958694786164e-10 0.1) // Vertex r4st = 55 
+    (-0.141421355856554 0.141421356618065 0.1) // Vertex r5t = 56 
+    (3.58979347393082e-10 0.2 0.1) // Vertex r6t = 57 
+    (3.58979347393082e-10 0.2 0.1) // Vertex r6st = 58 
+    (0.141421356364228 0.141421356110391 0.1) // Vertex r7t = 59 
+
+    (0.5 0 0.1) // Vertex rb0t = 60 
+    (0.353553390910569 -0.353553390275978 0.1) // Vertex rb1t = 61 
+    (8.97448368482705e-10 -0.5 0.1) // Vertex rb2t = 62 
+    (-0.353553389641386 -0.353553391545162 0.1) // Vertex rb3t = 63 
+    (-0.5 1.79489673696541e-09 0.1) // Vertex rb4t = 64 
+    (-0.353553389641386 0.353553391545162 0.1) // Vertex rb5t = 65 
+    (8.97448368482705e-10 0.5 0.1) // Vertex rb6t = 66 
+    (0.353553390910569 0.353553390275978 0.1) // Vertex rb7t = 67 
+
+    (0.6 0 0.1) // Vertex ri0t = 68 
+    (0.424264069092683 -0.424264068331174 0.1) // Vertex ri1t = 69 
+    (1.07693804217925e-09 -0.6 0.1) // Vertex ri2t = 70 
+    (-0.424264067569663 -0.424264069854194 0.1) // Vertex ri3t = 71 
+    (-0.6 2.15387608435849e-09 0.1) // Vertex ri4t = 72 
+    (-0.424264067569663 0.424264069854194 0.1) // Vertex ri5t = 73 
+    (1.07693804217925e-09 0.6 0.1) // Vertex ri6t = 74 
+    (0.424264069092683 0.424264068331174 0.1) // Vertex ri7t = 75 
+
+    (0.7 0 0.1) // Vertex Rb0t = 76 
+    (0.494974747274797 -0.494974746386369 0.1) // Vertex Rb1t = 77 
+    (1.25642771587579e-09 -0.7 0.1) // Vertex Rb2t = 78 
+    (-0.49497474549794 -0.494974748163226 0.1) // Vertex Rb3t = 79 
+    (-0.7 2.51285543175157e-09 0.1) // Vertex Rb4t = 80 
+    (-0.49497474549794 0.494974748163226 0.1) // Vertex Rb5t = 81 
+    (1.25642771587579e-09 0.7 0.1) // Vertex Rb6t = 82 
+    (0.494974747274797 0.494974746386369 0.1) // Vertex Rb7t = 83 
+
+    (1 0 0.1) // Vertex R0t = 84 
+    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1t = 85 
+    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1st = 86 
+    (1.79489673696541e-09 -1 0.1) // Vertex R2t = 87 
+    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3t = 88 
+    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3st = 89 
+    (-1 3.58979347393082e-09 0.1) // Vertex R4t = 90 
+    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5t = 91 
+    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5st = 92 
+    (1.79489673696541e-09 1 0.1) // Vertex R6t = 93 
+    (0.707106781821139 0.707106780551956 0.1) // Vertex R7t = 94 
+    (0.707106781821139 0.707106780551956 0.1) // Vertex R7st = 95 
+);
+
+blocks
+(
+    // block0
+    hex (0 2 13 12 48 50 61 60)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex (2 4 14 13 50 52 62 61)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex (3 5 15 14 51 53 63 62)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex (5 7 16 15 53 55 64 63)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex (6 8 17 16 54 56 65 64)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex (8 10 18 17 56 58 66 65)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex (9 11 19 18 57 59 67 66)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex (11 1 12 19 59 49 60 67)
+    rotor
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex (12 13 21 20 60 61 69 68)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex (13 14 22 21 61 62 70 69)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex (14 15 23 22 62 63 71 70)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex (15 16 24 23 63 64 72 71)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex (16 17 25 24 64 65 73 72)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex (17 18 26 25 65 66 74 73)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex (18 19 27 26 66 67 75 74)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex (19 12 20 27 67 60 68 75)
+    rotor
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex (20 21 29 28 68 69 77 76)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex (21 22 30 29 69 70 78 77)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex (22 23 31 30 70 71 79 78)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex (23 24 32 31 71 72 80 79)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex (24 25 33 32 72 73 81 80)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex (25 26 34 33 73 74 82 81)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex (26 27 35 34 74 75 83 82)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex (27 20 28 35 75 68 76 83)
+    stator
+    (12 4 1)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex (28 29 38 36 76 77 86 84)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex (29 30 39 37 77 78 87 85)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex (30 31 41 39 78 79 89 87)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex (31 32 42 40 79 80 90 88)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex (32 33 44 42 80 81 92 90)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex (33 34 45 43 81 82 93 91)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex (34 35 47 45 82 83 95 93)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex (35 28 36 46 83 76 84 94)
+    stator
+    (12 12 1)
+    simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc 0 2 (0.184775906536601 -0.0765366863901046 0)
+    arc 2 4 (0.0765366867217582 -0.184775906399226 0)
+    arc 3 5 (-0.0765366860584508 -0.184775906673977 0)
+    arc 5 7 (-0.18477590626185 -0.0765366870534118 0)
+    arc 6 8 (-0.18477590626185 0.0765366870534118 0)
+    arc 8 10 (-0.0765366860584508 0.184775906673977 0)
+    arc 9 11 (0.0765366867217582 0.184775906399226 0)
+    arc 11 1 (0.184775906536601 0.0765366863901046 0)
+
+    arc 12 13 (0.461939766341503 -0.191341715975262 0)
+    arc 13 14 (0.191341716804395 -0.461939765998065 0)
+    arc 14 15 (-0.191341715146127 -0.461939766684942 0)
+    arc 15 16 (-0.461939765654626 -0.19134171763353 0)
+    arc 16 17 (-0.461939765654626 0.19134171763353 0)
+    arc 17 18 (-0.191341715146127 0.461939766684942 0)
+    arc 18 19 (0.191341716804395 0.461939765998065 0)
+    arc 19 12 (0.461939766341503 0.191341715975262 0)
+
+    arc 20 21 (0.554327719609804 -0.229610059170314 0)
+    arc 21 22 (0.229610060165275 -0.554327719197677 0)
+    arc 22 23 (-0.229610058175352 -0.55432772002193 0)
+    arc 23 24 (-0.554327718785551 -0.229610061160235 0)
+    arc 24 25 (-0.554327718785551 0.229610061160235 0)
+    arc 25 26 (-0.229610058175352 0.55432772002193 0)
+    arc 26 27 (0.229610060165275 0.554327719197677 0)
+    arc 27 20 (0.554327719609804 0.229610059170314 0)
+
+    arc 28 29 (0.646715672878104 -0.267878402365366 0)
+    arc 29 30 (0.267878403526154 -0.64671567239729 0)
+    arc 30 31 (-0.267878401204578 -0.646715673358918 0)
+    arc 31 32 (-0.646715671916476 -0.267878404686941 0)
+    arc 32 33 (-0.646715671916476 0.267878404686941 0)
+    arc 33 34 (-0.267878401204578 0.646715673358918 0)
+    arc 34 35 (0.267878403526154 0.64671567239729 0)
+    arc 35 28 (0.646715672878104 0.267878402365366 0)
+
+    arc 36 38 (0.923879532683006 -0.382683431950523 0)
+    arc 37 39 (0.382683433608791 -0.923879531996129 0)
+    arc 39 41 (-0.382683430292254 -0.923879533369883 0)
+    arc 40 42 (-0.923879531309252 -0.382683435267059 0)
+    arc 42 44 (-0.923879531309252 0.382683435267059 0)
+    arc 43 45 (-0.382683430292254 0.923879533369883 0)
+    arc 45 47 (0.382683433608791 0.923879531996129 0)
+    arc 46 36 (0.923879532683006 0.382683431950523 0)
+
+    arc 48 50 (0.184775906536601 -0.0765366863901046 0.1)
+    arc 50 52 (0.0765366867217582 -0.184775906399226 0.1)
+    arc 51 53 (-0.0765366860584508 -0.184775906673977 0.1)
+    arc 53 55 (-0.18477590626185 -0.0765366870534118 0.1)
+    arc 54 56 (-0.18477590626185 0.0765366870534118 0.1)
+    arc 56 58 (-0.0765366860584508 0.184775906673977 0.1)
+    arc 57 59 (0.0765366867217582 0.184775906399226 0.1)
+    arc 59 49 (0.184775906536601 0.0765366863901046 0.1)
+
+    arc 60 61 (0.461939766341503 -0.191341715975262 0.1)
+    arc 61 62 (0.191341716804395 -0.461939765998065 0.1)
+    arc 62 63 (-0.191341715146127 -0.461939766684942 0.1)
+    arc 63 64 (-0.461939765654626 -0.19134171763353 0.1)
+    arc 64 65 (-0.461939765654626 0.19134171763353 0.1)
+    arc 65 66 (-0.191341715146127 0.461939766684942 0.1)
+    arc 66 67 (0.191341716804395 0.461939765998065 0.1)
+    arc 67 60 (0.461939766341503 0.191341715975262 0.1)
+
+    arc 68 69 (0.554327719609804 -0.229610059170314 0.1)
+    arc 69 70 (0.229610060165275 -0.554327719197677 0.1)
+    arc 70 71 (-0.229610058175352 -0.55432772002193 0.1)
+    arc 71 72 (-0.554327718785551 -0.229610061160235 0.1)
+    arc 72 73 (-0.554327718785551 0.229610061160235 0.1)
+    arc 73 74 (-0.229610058175352 0.55432772002193 0.1)
+    arc 74 75 (0.229610060165275 0.554327719197677 0.1)
+    arc 75 68 (0.554327719609804 0.229610059170314 0.1)
+
+    arc 76 77 (0.646715672878104 -0.267878402365366 0.1)
+    arc 77 78 (0.267878403526154 -0.64671567239729 0.1)
+    arc 78 79 (-0.267878401204578 -0.646715673358918 0.1)
+    arc 79 80 (-0.646715671916476 -0.267878404686941 0.1)
+    arc 80 81 (-0.646715671916476 0.267878404686941 0.1)
+    arc 81 82 (-0.267878401204578 0.646715673358918 0.1)
+    arc 82 83 (0.267878403526154 0.64671567239729 0.1)
+    arc 83 76 (0.646715672878104 0.267878402365366 0.1)
+
+    arc 84 86 (0.923879532683006 -0.382683431950523 0.1)
+    arc 85 87 (0.382683433608791 -0.923879531996129 0.1)
+    arc 87 89 (-0.382683430292254 -0.923879533369883 0.1)
+    arc 88 90 (-0.923879531309252 -0.382683435267059 0.1)
+    arc 90 92 (-0.923879531309252 0.382683435267059 0.1)
+    arc 91 93 (-0.382683430292254 0.923879533369883 0.1)
+    arc 93 95 (0.382683433608791 0.923879531996129 0.1)
+    arc 94 84 (0.923879532683006 0.382683431950523 0.1)
+);
+
+patches
+(
+    wall rotor
+    (
+        (0 2 50 48)
+        (2 4 52 50)
+        (3 5 53 51)
+        (5 7 55 53)
+        (6 8 56 54)
+        (8 10 58 56)
+        (9 11 59 57)
+        (11 1 49 59)
+
+        (0 12 60 48)
+        (1 12 60 49)
+
+        (3 14 62 51)
+        (4 14 62 52)
+
+        (6 16 64 54)
+        (7 16 64 55)
+
+        (9 18 66 57)
+        (10 18 66 58)
+    )
+
+    wall stator
+    (
+        (36 38 86 84)
+        (37 39 87 85)
+        (39 41 89 87)
+        (40 42 90 88)
+        (42 44 92 90)
+        (43 45 93 91)
+        (45 47 95 93)
+        (46 36 84 94)
+
+        (37 29 77 85)
+        (38 29 77 86)
+
+        (40 31 79 88)
+        (41 31 79 89)
+
+        (43 33 81 91)
+        (44 33 81 92)
+
+        (46 35 83 94)
+        (47 35 83 95)
+    )
+
+    empty front
+    (
+        (48 50 61 60)
+        (50 52 62 61)
+        (51 53 63 62)
+        (53 55 64 63)
+        (54 56 65 64)
+        (56 58 66 65)
+        (57 59 67 66)
+        (59 49 60 67)
+        (60 61 69 68)
+        (61 62 70 69)
+        (62 63 71 70)
+        (63 64 72 71)
+        (64 65 73 72)
+        (65 66 74 73)
+        (66 67 75 74)
+        (67 60 68 75)
+        (68 69 77 76)
+        (69 70 78 77)
+        (70 71 79 78)
+        (71 72 80 79)
+        (72 73 81 80)
+        (73 74 82 81)
+        (74 75 83 82)
+        (75 68 76 83)
+        (76 77 86 84)
+        (77 78 87 85)
+        (78 79 89 87)
+        (79 80 90 88)
+        (80 81 92 90)
+        (81 82 93 91)
+        (82 83 95 93)
+        (83 76 84 94)
+    )
+
+    empty back
+    (
+        (0 12 13 2)
+        (2 13 14 4)
+        (3 14 15 5)
+        (5 15 16 7)
+        (6 16 17 8)
+        (8 17 18 10)
+        (9 18 19 11)
+        (11 19 12 1)
+        (12 20 21 13)
+        (13 21 22 14)
+        (14 22 23 15)
+        (15 23 24 16)
+        (16 24 25 17)
+        (17 25 26 18)
+        (18 26 27 19)
+        (19 27 20 12)
+        (20 28 29 21)
+        (21 29 30 22)
+        (22 30 31 23)
+        (23 31 32 24)
+        (24 32 33 25)
+        (25 33 34 26)
+        (26 34 35 27)
+        (27 35 28 20)
+        (28 36 38 29)
+        (29 37 39 30)
+        (30 39 41 31)
+        (31 40 42 32)
+        (32 42 44 33)
+        (33 43 45 34)
+        (34 45 47 35)
+        (35 46 36 28)
+    )
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4 b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000000000000000000000000000000000..dccc042d22493ca6d40e25d7cb8ba8c76f8857aa
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,834 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// General macros to create 2D/extruded-2D meshes
+
+changecom(//)changequote([,])
+define(calc, [esyscmd(perl -e 'print ($1)')])
+define(VCOUNT, 0)
+define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])
+define(pi, 3.14159265)
+
+define(hex2D, hex ($1b $2b $3b $4b $1t $2t $3t $4t))
+define(quad2D, ($1b $2b $2t $1t))
+define(frontQuad, ($1t $2t $3t $4t))
+define(backQuad, ($1b $4b $3b $2b))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+// Hub radius
+define(r, 0.2)
+
+// Impeller-tip radius
+define(rb, 0.5)
+
+// Baffle-tip radius
+define(Rb, 0.7)
+
+// Tank radius
+define(R, 1)
+
+// MRF region radius
+define(ri, calc(0.5*(rb + Rb)))
+
+// Thickness of 2D slab
+define(z, 0.1)
+
+// Base z
+define(Zb, 0)
+
+// Top z
+define(Zt, calc(Zb + z))
+
+// Number of cells radially between hub and impeller tip
+define(Nr, 12)
+
+// Number of cells radially in each of the two regions between
+// impeller and baffle tips
+define(Ni, 4)
+
+// Number of cells radially between baffle tip and tank
+define(NR, 12)
+
+// Number of cells azimuthally in each of the 8 blocks
+define(Na, 12)
+
+// Number of cells in the thickness of the slab
+define(Nz, 1)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+define(vert, (x$1$2 y$1$2 $3))
+define(evert, (ex$1$2 ey$1$2 $3))
+
+define(a0, 0)
+define(a1, -45)
+define(a2, -90)
+define(a3, -135)
+define(a4, 180)
+define(a5, 135)
+define(a6, 90)
+define(a7, 45)
+
+define(ea0, -22.5)
+define(ea1, -67.5)
+define(ea2, -112.5)
+define(ea3, -157.5)
+define(ea4, 157.5)
+define(ea5, 112.5)
+define(ea6, 67.5)
+define(ea7, 22.5)
+
+define(ca0, calc(cos((pi/180)*a0)))
+define(ca1, calc(cos((pi/180)*a1)))
+define(ca2, calc(cos((pi/180)*a2)))
+define(ca3, calc(cos((pi/180)*a3)))
+define(ca4, calc(cos((pi/180)*a4)))
+define(ca5, calc(cos((pi/180)*a5)))
+define(ca6, calc(cos((pi/180)*a6)))
+define(ca7, calc(cos((pi/180)*a7)))
+
+define(sa0, calc(sin((pi/180)*a0)))
+define(sa1, calc(sin((pi/180)*a1)))
+define(sa2, calc(sin((pi/180)*a2)))
+define(sa3, calc(sin((pi/180)*a3)))
+define(sa4, calc(sin((pi/180)*a4)))
+define(sa5, calc(sin((pi/180)*a5)))
+define(sa6, calc(sin((pi/180)*a6)))
+define(sa7, calc(sin((pi/180)*a7)))
+
+define(cea0, calc(cos((pi/180)*ea0)))
+define(cea1, calc(cos((pi/180)*ea1)))
+define(cea2, calc(cos((pi/180)*ea2)))
+define(cea3, calc(cos((pi/180)*ea3)))
+define(cea4, calc(cos((pi/180)*ea4)))
+define(cea5, calc(cos((pi/180)*ea5)))
+define(cea6, calc(cos((pi/180)*ea6)))
+define(cea7, calc(cos((pi/180)*ea7)))
+
+define(sea0, calc(sin((pi/180)*ea0)))
+define(sea1, calc(sin((pi/180)*ea1)))
+define(sea2, calc(sin((pi/180)*ea2)))
+define(sea3, calc(sin((pi/180)*ea3)))
+define(sea4, calc(sin((pi/180)*ea4)))
+define(sea5, calc(sin((pi/180)*ea5)))
+define(sea6, calc(sin((pi/180)*ea6)))
+define(sea7, calc(sin((pi/180)*ea7)))
+
+define(x00, calc(r*ca0))
+define(x01, calc(r*ca1))
+define(x02, calc(r*ca2))
+define(x03, calc(r*ca3))
+define(x04, calc(r*ca4))
+define(x05, calc(r*ca5))
+define(x06, calc(r*ca6))
+define(x07, calc(r*ca7))
+
+define(x10, calc(rb*ca0))
+define(x11, calc(rb*ca1))
+define(x12, calc(rb*ca2))
+define(x13, calc(rb*ca3))
+define(x14, calc(rb*ca4))
+define(x15, calc(rb*ca5))
+define(x16, calc(rb*ca6))
+define(x17, calc(rb*ca7))
+
+define(x20, calc(ri*ca0))
+define(x21, calc(ri*ca1))
+define(x22, calc(ri*ca2))
+define(x23, calc(ri*ca3))
+define(x24, calc(ri*ca4))
+define(x25, calc(ri*ca5))
+define(x26, calc(ri*ca6))
+define(x27, calc(ri*ca7))
+
+define(x30, calc(Rb*ca0))
+define(x31, calc(Rb*ca1))
+define(x32, calc(Rb*ca2))
+define(x33, calc(Rb*ca3))
+define(x34, calc(Rb*ca4))
+define(x35, calc(Rb*ca5))
+define(x36, calc(Rb*ca6))
+define(x37, calc(Rb*ca7))
+
+define(x40, calc(R*ca0))
+define(x41, calc(R*ca1))
+define(x42, calc(R*ca2))
+define(x43, calc(R*ca3))
+define(x44, calc(R*ca4))
+define(x45, calc(R*ca5))
+define(x46, calc(R*ca6))
+define(x47, calc(R*ca7))
+
+define(y00, calc(r*sa0))
+define(y01, calc(r*sa1))
+define(y02, calc(r*sa2))
+define(y03, calc(r*sa3))
+define(y04, calc(r*sa4))
+define(y05, calc(r*sa5))
+define(y06, calc(r*sa6))
+define(y07, calc(r*sa7))
+
+define(y10, calc(rb*sa0))
+define(y11, calc(rb*sa1))
+define(y12, calc(rb*sa2))
+define(y13, calc(rb*sa3))
+define(y14, calc(rb*sa4))
+define(y15, calc(rb*sa5))
+define(y16, calc(rb*sa6))
+define(y17, calc(rb*sa7))
+
+define(y20, calc(ri*sa0))
+define(y21, calc(ri*sa1))
+define(y22, calc(ri*sa2))
+define(y23, calc(ri*sa3))
+define(y24, calc(ri*sa4))
+define(y25, calc(ri*sa5))
+define(y26, calc(ri*sa6))
+define(y27, calc(ri*sa7))
+
+define(y30, calc(Rb*sa0))
+define(y31, calc(Rb*sa1))
+define(y32, calc(Rb*sa2))
+define(y33, calc(Rb*sa3))
+define(y34, calc(Rb*sa4))
+define(y35, calc(Rb*sa5))
+define(y36, calc(Rb*sa6))
+define(y37, calc(Rb*sa7))
+
+define(y40, calc(R*sa0))
+define(y41, calc(R*sa1))
+define(y42, calc(R*sa2))
+define(y43, calc(R*sa3))
+define(y44, calc(R*sa4))
+define(y45, calc(R*sa5))
+define(y46, calc(R*sa6))
+define(y47, calc(R*sa7))
+
+define(ex00, calc(r*cea0))
+define(ex01, calc(r*cea1))
+define(ex02, calc(r*cea2))
+define(ex03, calc(r*cea3))
+define(ex04, calc(r*cea4))
+define(ex05, calc(r*cea5))
+define(ex06, calc(r*cea6))
+define(ex07, calc(r*cea7))
+
+define(ex10, calc(rb*cea0))
+define(ex11, calc(rb*cea1))
+define(ex12, calc(rb*cea2))
+define(ex13, calc(rb*cea3))
+define(ex14, calc(rb*cea4))
+define(ex15, calc(rb*cea5))
+define(ex16, calc(rb*cea6))
+define(ex17, calc(rb*cea7))
+
+define(ex20, calc(ri*cea0))
+define(ex21, calc(ri*cea1))
+define(ex22, calc(ri*cea2))
+define(ex23, calc(ri*cea3))
+define(ex24, calc(ri*cea4))
+define(ex25, calc(ri*cea5))
+define(ex26, calc(ri*cea6))
+define(ex27, calc(ri*cea7))
+
+define(ex30, calc(Rb*cea0))
+define(ex31, calc(Rb*cea1))
+define(ex32, calc(Rb*cea2))
+define(ex33, calc(Rb*cea3))
+define(ex34, calc(Rb*cea4))
+define(ex35, calc(Rb*cea5))
+define(ex36, calc(Rb*cea6))
+define(ex37, calc(Rb*cea7))
+
+define(ex40, calc(R*cea0))
+define(ex41, calc(R*cea1))
+define(ex42, calc(R*cea2))
+define(ex43, calc(R*cea3))
+define(ex44, calc(R*cea4))
+define(ex45, calc(R*cea5))
+define(ex46, calc(R*cea6))
+define(ex47, calc(R*cea7))
+
+define(ey00, calc(r*sea0))
+define(ey01, calc(r*sea1))
+define(ey02, calc(r*sea2))
+define(ey03, calc(r*sea3))
+define(ey04, calc(r*sea4))
+define(ey05, calc(r*sea5))
+define(ey06, calc(r*sea6))
+define(ey07, calc(r*sea7))
+
+define(ey10, calc(rb*sea0))
+define(ey11, calc(rb*sea1))
+define(ey12, calc(rb*sea2))
+define(ey13, calc(rb*sea3))
+define(ey14, calc(rb*sea4))
+define(ey15, calc(rb*sea5))
+define(ey16, calc(rb*sea6))
+define(ey17, calc(rb*sea7))
+
+define(ey20, calc(ri*sea0))
+define(ey21, calc(ri*sea1))
+define(ey22, calc(ri*sea2))
+define(ey23, calc(ri*sea3))
+define(ey24, calc(ri*sea4))
+define(ey25, calc(ri*sea5))
+define(ey26, calc(ri*sea6))
+define(ey27, calc(ri*sea7))
+
+define(ey30, calc(Rb*sea0))
+define(ey31, calc(Rb*sea1))
+define(ey32, calc(Rb*sea2))
+define(ey33, calc(Rb*sea3))
+define(ey34, calc(Rb*sea4))
+define(ey35, calc(Rb*sea5))
+define(ey36, calc(Rb*sea6))
+define(ey37, calc(Rb*sea7))
+
+define(ey40, calc(R*sea0))
+define(ey41, calc(R*sea1))
+define(ey42, calc(R*sea2))
+define(ey43, calc(R*sea3))
+define(ey44, calc(R*sea4))
+define(ey45, calc(R*sea5))
+define(ey46, calc(R*sea6))
+define(ey47, calc(R*sea7))
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+vertices
+(
+    vert(0, 0, Zb) vlabel(r0b)
+    vert(0, 0, Zb) vlabel(r0sb)
+    vert(0, 1, Zb) vlabel(r1b)
+    vert(0, 2, Zb) vlabel(r2b)
+    vert(0, 2, Zb) vlabel(r2sb)
+    vert(0, 3, Zb) vlabel(r3b)
+    vert(0, 4, Zb) vlabel(r4b)
+    vert(0, 4, Zb) vlabel(r4sb)
+    vert(0, 5, Zb) vlabel(r5b)
+    vert(0, 6, Zb) vlabel(r6b)
+    vert(0, 6, Zb) vlabel(r6sb)
+    vert(0, 7, Zb) vlabel(r7b)
+
+    vert(1, 0, Zb) vlabel(rb0b)
+    vert(1, 1, Zb) vlabel(rb1b)
+    vert(1, 2, Zb) vlabel(rb2b)
+    vert(1, 3, Zb) vlabel(rb3b)
+    vert(1, 4, Zb) vlabel(rb4b)
+    vert(1, 5, Zb) vlabel(rb5b)
+    vert(1, 6, Zb) vlabel(rb6b)
+    vert(1, 7, Zb) vlabel(rb7b)
+
+    vert(2, 0, Zb) vlabel(ri0b)
+    vert(2, 1, Zb) vlabel(ri1b)
+    vert(2, 2, Zb) vlabel(ri2b)
+    vert(2, 3, Zb) vlabel(ri3b)
+    vert(2, 4, Zb) vlabel(ri4b)
+    vert(2, 5, Zb) vlabel(ri5b)
+    vert(2, 6, Zb) vlabel(ri6b)
+    vert(2, 7, Zb) vlabel(ri7b)
+
+    vert(3, 0, Zb) vlabel(Rb0b)
+    vert(3, 1, Zb) vlabel(Rb1b)
+    vert(3, 2, Zb) vlabel(Rb2b)
+    vert(3, 3, Zb) vlabel(Rb3b)
+    vert(3, 4, Zb) vlabel(Rb4b)
+    vert(3, 5, Zb) vlabel(Rb5b)
+    vert(3, 6, Zb) vlabel(Rb6b)
+    vert(3, 7, Zb) vlabel(Rb7b)
+
+    vert(4, 0, Zb) vlabel(R0b)
+    vert(4, 1, Zb) vlabel(R1b)
+    vert(4, 1, Zb) vlabel(R1sb)
+    vert(4, 2, Zb) vlabel(R2b)
+    vert(4, 3, Zb) vlabel(R3b)
+    vert(4, 3, Zb) vlabel(R3sb)
+    vert(4, 4, Zb) vlabel(R4b)
+    vert(4, 5, Zb) vlabel(R5b)
+    vert(4, 5, Zb) vlabel(R5sb)
+    vert(4, 6, Zb) vlabel(R6b)
+    vert(4, 7, Zb) vlabel(R7b)
+    vert(4, 7, Zb) vlabel(R7sb)
+
+    vert(0, 0, Zt) vlabel(r0t)
+    vert(0, 0, Zt) vlabel(r0st)
+    vert(0, 1, Zt) vlabel(r1t)
+    vert(0, 2, Zt) vlabel(r2t)
+    vert(0, 2, Zt) vlabel(r2st)
+    vert(0, 3, Zt) vlabel(r3t)
+    vert(0, 4, Zt) vlabel(r4t)
+    vert(0, 4, Zt) vlabel(r4st)
+    vert(0, 5, Zt) vlabel(r5t)
+    vert(0, 6, Zt) vlabel(r6t)
+    vert(0, 6, Zt) vlabel(r6st)
+    vert(0, 7, Zt) vlabel(r7t)
+
+    vert(1, 0, Zt) vlabel(rb0t)
+    vert(1, 1, Zt) vlabel(rb1t)
+    vert(1, 2, Zt) vlabel(rb2t)
+    vert(1, 3, Zt) vlabel(rb3t)
+    vert(1, 4, Zt) vlabel(rb4t)
+    vert(1, 5, Zt) vlabel(rb5t)
+    vert(1, 6, Zt) vlabel(rb6t)
+    vert(1, 7, Zt) vlabel(rb7t)
+
+    vert(2, 0, Zt) vlabel(ri0t)
+    vert(2, 1, Zt) vlabel(ri1t)
+    vert(2, 2, Zt) vlabel(ri2t)
+    vert(2, 3, Zt) vlabel(ri3t)
+    vert(2, 4, Zt) vlabel(ri4t)
+    vert(2, 5, Zt) vlabel(ri5t)
+    vert(2, 6, Zt) vlabel(ri6t)
+    vert(2, 7, Zt) vlabel(ri7t)
+
+    vert(3, 0, Zt) vlabel(Rb0t)
+    vert(3, 1, Zt) vlabel(Rb1t)
+    vert(3, 2, Zt) vlabel(Rb2t)
+    vert(3, 3, Zt) vlabel(Rb3t)
+    vert(3, 4, Zt) vlabel(Rb4t)
+    vert(3, 5, Zt) vlabel(Rb5t)
+    vert(3, 6, Zt) vlabel(Rb6t)
+    vert(3, 7, Zt) vlabel(Rb7t)
+
+    vert(4, 0, Zt) vlabel(R0t)
+    vert(4, 1, Zt) vlabel(R1t)
+    vert(4, 1, Zt) vlabel(R1st)
+    vert(4, 2, Zt) vlabel(R2t)
+    vert(4, 3, Zt) vlabel(R3t)
+    vert(4, 3, Zt) vlabel(R3st)
+    vert(4, 4, Zt) vlabel(R4t)
+    vert(4, 5, Zt) vlabel(R5t)
+    vert(4, 5, Zt) vlabel(R5st)
+    vert(4, 6, Zt) vlabel(R6t)
+    vert(4, 7, Zt) vlabel(R7t)
+    vert(4, 7, Zt) vlabel(R7st)
+);
+
+blocks
+(
+    // block0
+    hex2D(r0, r1, rb1, rb0)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(r1, r2s, rb2, rb1)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(r2, r3, rb3, rb2)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(r3, r4s, rb4, rb3)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(r4, r5, rb5, rb4)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(r5, r6s, rb6, rb5)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(r6, r7, rb7, rb6)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(r7, r0s, rb0, rb7)
+    rotor
+    (Na Nr Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(rb0, rb1, ri1, ri0)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(rb1, rb2, ri2, ri1)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(rb2, rb3, ri3, ri2)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(rb3, rb4, ri4, ri3)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(rb4, rb5, ri5, ri4)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(rb5, rb6, ri6, ri5)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(rb6, rb7, ri7, ri6)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(rb7, rb0, ri0, ri7)
+    rotor
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(ri0, ri1, Rb1, Rb0)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(ri1, ri2, Rb2, Rb1)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(ri2, ri3, Rb3, Rb2)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(ri3, ri4, Rb4, Rb3)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(ri4, ri5, Rb5, Rb4)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(ri5, ri6, Rb6, Rb5)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(ri6, ri7, Rb7, Rb6)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(ri7, ri0, Rb0, Rb7)
+    stator
+    (Na Ni Nz)
+    simpleGrading (1 1 1)
+
+    // block0
+    hex2D(Rb0, Rb1, R1s, R0)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block1
+    hex2D(Rb1, Rb2, R2, R1)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block2
+    hex2D(Rb2, Rb3, R3s, R2)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block3
+    hex2D(Rb3, Rb4, R4, R3)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block4
+    hex2D(Rb4, Rb5, R5s, R4)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block5
+    hex2D(Rb5, Rb6, R6, R5)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block6
+    hex2D(Rb6, Rb7, R7s, R6)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+
+    // block7
+    hex2D(Rb7, Rb0, R0, R7)
+    stator
+    (Na NR Nz)
+    simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc r0b r1b evert(0, 0, Zb)
+    arc r1b r2sb evert(0, 1, Zb)
+    arc r2b r3b evert(0, 2, Zb)
+    arc r3b r4sb evert(0, 3, Zb)
+    arc r4b r5b evert(0, 4, Zb)
+    arc r5b r6sb evert(0, 5, Zb)
+    arc r6b r7b evert(0, 6, Zb)
+    arc r7b r0sb evert(0, 7, Zb)
+
+    arc rb0b rb1b evert(1, 0, Zb)
+    arc rb1b rb2b evert(1, 1, Zb)
+    arc rb2b rb3b evert(1, 2, Zb)
+    arc rb3b rb4b evert(1, 3, Zb)
+    arc rb4b rb5b evert(1, 4, Zb)
+    arc rb5b rb6b evert(1, 5, Zb)
+    arc rb6b rb7b evert(1, 6, Zb)
+    arc rb7b rb0b evert(1, 7, Zb)
+
+    arc ri0b ri1b evert(2, 0, Zb)
+    arc ri1b ri2b evert(2, 1, Zb)
+    arc ri2b ri3b evert(2, 2, Zb)
+    arc ri3b ri4b evert(2, 3, Zb)
+    arc ri4b ri5b evert(2, 4, Zb)
+    arc ri5b ri6b evert(2, 5, Zb)
+    arc ri6b ri7b evert(2, 6, Zb)
+    arc ri7b ri0b evert(2, 7, Zb)
+
+    arc Rb0b Rb1b evert(3, 0, Zb)
+    arc Rb1b Rb2b evert(3, 1, Zb)
+    arc Rb2b Rb3b evert(3, 2, Zb)
+    arc Rb3b Rb4b evert(3, 3, Zb)
+    arc Rb4b Rb5b evert(3, 4, Zb)
+    arc Rb5b Rb6b evert(3, 5, Zb)
+    arc Rb6b Rb7b evert(3, 6, Zb)
+    arc Rb7b Rb0b evert(3, 7, Zb)
+
+    arc R0b R1sb evert(4, 0, Zb)
+    arc R1b R2b evert(4, 1, Zb)
+    arc R2b R3sb evert(4, 2, Zb)
+    arc R3b R4b evert(4, 3, Zb)
+    arc R4b R5sb evert(4, 4, Zb)
+    arc R5b R6b evert(4, 5, Zb)
+    arc R6b R7sb evert(4, 6, Zb)
+    arc R7b R0b evert(4, 7, Zb)
+
+    arc r0t r1t evert(0, 0, Zt)
+    arc r1t r2st evert(0, 1, Zt)
+    arc r2t r3t evert(0, 2, Zt)
+    arc r3t r4st evert(0, 3, Zt)
+    arc r4t r5t evert(0, 4, Zt)
+    arc r5t r6st evert(0, 5, Zt)
+    arc r6t r7t evert(0, 6, Zt)
+    arc r7t r0st evert(0, 7, Zt)
+
+    arc rb0t rb1t evert(1, 0, Zt)
+    arc rb1t rb2t evert(1, 1, Zt)
+    arc rb2t rb3t evert(1, 2, Zt)
+    arc rb3t rb4t evert(1, 3, Zt)
+    arc rb4t rb5t evert(1, 4, Zt)
+    arc rb5t rb6t evert(1, 5, Zt)
+    arc rb6t rb7t evert(1, 6, Zt)
+    arc rb7t rb0t evert(1, 7, Zt)
+
+    arc ri0t ri1t evert(2, 0, Zt)
+    arc ri1t ri2t evert(2, 1, Zt)
+    arc ri2t ri3t evert(2, 2, Zt)
+    arc ri3t ri4t evert(2, 3, Zt)
+    arc ri4t ri5t evert(2, 4, Zt)
+    arc ri5t ri6t evert(2, 5, Zt)
+    arc ri6t ri7t evert(2, 6, Zt)
+    arc ri7t ri0t evert(2, 7, Zt)
+
+    arc Rb0t Rb1t evert(3, 0, Zt)
+    arc Rb1t Rb2t evert(3, 1, Zt)
+    arc Rb2t Rb3t evert(3, 2, Zt)
+    arc Rb3t Rb4t evert(3, 3, Zt)
+    arc Rb4t Rb5t evert(3, 4, Zt)
+    arc Rb5t Rb6t evert(3, 5, Zt)
+    arc Rb6t Rb7t evert(3, 6, Zt)
+    arc Rb7t Rb0t evert(3, 7, Zt)
+
+    arc R0t R1st evert(4, 0, Zt)
+    arc R1t R2t evert(4, 1, Zt)
+    arc R2t R3st evert(4, 2, Zt)
+    arc R3t R4t evert(4, 3, Zt)
+    arc R4t R5st evert(4, 4, Zt)
+    arc R5t R6t evert(4, 5, Zt)
+    arc R6t R7st evert(4, 6, Zt)
+    arc R7t R0t evert(4, 7, Zt)
+);
+
+patches
+(
+    wall rotor
+    (
+        quad2D(r0, r1)
+        quad2D(r1, r2s)
+        quad2D(r2, r3)
+        quad2D(r3, r4s)
+        quad2D(r4, r5)
+        quad2D(r5, r6s)
+        quad2D(r6, r7)
+        quad2D(r7, r0s)
+
+        quad2D(r0, rb0)
+        quad2D(r0s, rb0)
+
+        quad2D(r2, rb2)
+        quad2D(r2s, rb2)
+
+        quad2D(r4, rb4)
+        quad2D(r4s, rb4)
+
+        quad2D(r6, rb6)
+        quad2D(r6s, rb6)
+    )
+
+    wall stator
+    (
+        quad2D(R0, R1s)
+        quad2D(R1, R2)
+        quad2D(R2, R3s)
+        quad2D(R3, R4)
+        quad2D(R4, R5s)
+        quad2D(R5, R6)
+        quad2D(R6, R7s)
+        quad2D(R7, R0)
+
+        quad2D(R1, Rb1)
+        quad2D(R1s, Rb1)
+
+        quad2D(R3, Rb3)
+        quad2D(R3s, Rb3)
+
+        quad2D(R5, Rb5)
+        quad2D(R5s, Rb5)
+
+        quad2D(R7, Rb7)
+        quad2D(R7s, Rb7)
+    )
+
+    empty front
+    (
+        frontQuad(r0, r1, rb1, rb0)
+        frontQuad(r1, r2s, rb2, rb1)
+        frontQuad(r2, r3, rb3, rb2)
+        frontQuad(r3, r4s, rb4, rb3)
+        frontQuad(r4, r5, rb5, rb4)
+        frontQuad(r5, r6s, rb6, rb5)
+        frontQuad(r6, r7, rb7, rb6)
+        frontQuad(r7, r0s, rb0, rb7)
+        frontQuad(rb0, rb1, ri1, ri0)
+        frontQuad(rb1, rb2, ri2, ri1)
+        frontQuad(rb2, rb3, ri3, ri2)
+        frontQuad(rb3, rb4, ri4, ri3)
+        frontQuad(rb4, rb5, ri5, ri4)
+        frontQuad(rb5, rb6, ri6, ri5)
+        frontQuad(rb6, rb7, ri7, ri6)
+        frontQuad(rb7, rb0, ri0, ri7)
+        frontQuad(ri0, ri1, Rb1, Rb0)
+        frontQuad(ri1, ri2, Rb2, Rb1)
+        frontQuad(ri2, ri3, Rb3, Rb2)
+        frontQuad(ri3, ri4, Rb4, Rb3)
+        frontQuad(ri4, ri5, Rb5, Rb4)
+        frontQuad(ri5, ri6, Rb6, Rb5)
+        frontQuad(ri6, ri7, Rb7, Rb6)
+        frontQuad(ri7, ri0, Rb0, Rb7)
+        frontQuad(Rb0, Rb1, R1s, R0)
+        frontQuad(Rb1, Rb2, R2, R1)
+        frontQuad(Rb2, Rb3, R3s, R2)
+        frontQuad(Rb3, Rb4, R4, R3)
+        frontQuad(Rb4, Rb5, R5s, R4)
+        frontQuad(Rb5, Rb6, R6, R5)
+        frontQuad(Rb6, Rb7, R7s, R6)
+        frontQuad(Rb7, Rb0, R0, R7)
+    )
+
+    empty back
+    (
+        backQuad(r0, r1, rb1, rb0)
+        backQuad(r1, r2s, rb2, rb1)
+        backQuad(r2, r3, rb3, rb2)
+        backQuad(r3, r4s, rb4, rb3)
+        backQuad(r4, r5, rb5, rb4)
+        backQuad(r5, r6s, rb6, rb5)
+        backQuad(r6, r7, rb7, rb6)
+        backQuad(r7, r0s, rb0, rb7)
+        backQuad(rb0, rb1, ri1, ri0)
+        backQuad(rb1, rb2, ri2, ri1)
+        backQuad(rb2, rb3, ri3, ri2)
+        backQuad(rb3, rb4, ri4, ri3)
+        backQuad(rb4, rb5, ri5, ri4)
+        backQuad(rb5, rb6, ri6, ri5)
+        backQuad(rb6, rb7, ri7, ri6)
+        backQuad(rb7, rb0, ri0, ri7)
+        backQuad(ri0, ri1, Rb1, Rb0)
+        backQuad(ri1, ri2, Rb2, Rb1)
+        backQuad(ri2, ri3, Rb3, Rb2)
+        backQuad(ri3, ri4, Rb4, Rb3)
+        backQuad(ri4, ri5, Rb5, Rb4)
+        backQuad(ri5, ri6, Rb6, Rb5)
+        backQuad(ri6, ri7, Rb7, Rb6)
+        backQuad(ri7, ri0, Rb0, Rb7)
+        backQuad(Rb0, Rb1, R1s, R0)
+        backQuad(Rb1, Rb2, R2, R1)
+        backQuad(Rb2, Rb3, R3s, R2)
+        backQuad(Rb3, Rb4, R4, R3)
+        backQuad(Rb4, Rb5, R5s, R4)
+        backQuad(Rb5, Rb6, R6, R5)
+        backQuad(Rb6, Rb7, R7s, R6)
+        backQuad(Rb7, Rb0, R0, R7)
+    )
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..292f25b806357d9df75c7731f74dee0ec0aa3a40
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    rotor
+    {
+        type            wall;
+        nFaces          192;
+        startFace       5952;
+    }
+    stator
+    {
+        type            wall;
+        nFaces          192;
+        startFace       6144;
+    }
+    front
+    {
+        type            empty;
+        nFaces          3072;
+        startFace       6336;
+    }
+    back
+    {
+        type            empty;
+        nFaces          3072;
+        startFace       9408;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porousZones b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/porousZones
new file mode 100644
index 0000000000000000000000000000000000000000..e7986b4d348fbc45c5ca68d87a78a78043636f9f
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/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
+(
+    stator
+    {
+        coordinateSystem
+        {
+            e1  (1 0 0);
+            e2  (0 1 0);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (1e5 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/thermophysicalProperties b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..a2c84a89d38005a697cd9d8b3bcce3ada3e0db44
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/thermophysicalProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- 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      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+mixture         air 1 28.9 1007 0 1.4792e-06 116;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/transportProperties b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..36e9669391c8c0057688bb7b2bc38072dbd5ce5f
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/transportProperties
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 1e-05;
+
+CrossPowerLawCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    m               m [ 0 0 1 0 0 0 0 ] 1;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+    k               k [ 0 0 1 0 0 0 0 ] 0;
+    n               n [ 0 0 0 0 0 0 0 ] 1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/turbulenceProperties b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..f6753662e27e414f4bfb666745b7888fef8510bb
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/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  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/makeMesh b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/makeMesh
new file mode 100755
index 0000000000000000000000000000000000000000..91a75f1c48d48c5e478102002e4bfe6800858999
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/makeMesh
@@ -0,0 +1,5 @@
+#!/bin/sh
+set -x
+
+m4 < constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+blockMesh >& log.blockMesh
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/controlDict b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..b79e14b234e0191949729846152720a9687067e5
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/controlDict
@@ -0,0 +1,52 @@
+/*--------------------------------*- 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     rhoPorousMRFPimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.1;
+
+deltaT          1e-5;
+
+writeControl    adjustableRunTime;
+
+writeInterval   1e-2;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           1;
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict
new file mode 100644
index 0000000000000000000000000000000000000000..d4c9e000656215078417d03c0b5ee507e6e6e6fe
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict
@@ -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      ascii;
+    class       dictionary;
+    location    "system";
+    object      faceSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+name            rotor;
+
+action          delete;
+
+topoSetSources  ( boundaryToFace { } );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_noBoundaryFaces b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_noBoundaryFaces
new file mode 100644
index 0000000000000000000000000000000000000000..38fa0208ed993427b2efc3666064c6b779a76487
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_noBoundaryFaces
@@ -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      ascii;
+    class       dictionary;
+    location    "system";
+    object      faceSetDict_noBoundaryFaces;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+name            rotor;
+
+action          delete;
+
+topoSetSources  ( boundaryToFace { } );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_rotorFaces b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_rotorFaces
new file mode 100644
index 0000000000000000000000000000000000000000..3ccb7ad761db3fee6092b9ae449cec5f39dab7d7
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/faceSetDict_rotorFaces
@@ -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      ascii;
+    class       dictionary;
+    location    "system";
+    object      faceSetDict_rotorFaces;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+name            rotor;
+
+action          new;
+
+topoSetSources  ( cellToFace { set rotor ; option all ; } );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..2d1d7ea99d6b410ab7a3ab948f3d06de3823b2b4
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
@@ -0,0 +1,67 @@
+/*--------------------------------*- 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
+{
+    default         none;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,h)      Gauss limitedLinear 1;
+    div(phi,k)      Gauss limitedLinear 1;
+    div(phi,epsilon) Gauss limitedLinear 1;
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div(phiU,p) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..70ef25731da8e32390565e6cbafde07d701c1e98
--- /dev/null
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
@@ -0,0 +1,121 @@
+/*--------------------------------*- 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
+{
+    rho
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-6;
+        relTol          0.05;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    pFinal
+    {
+        solver          GAMG;
+        tolerance       1e-6;
+        relTol          0;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    h
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0.1;
+    }
+
+    hFinal
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    k
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0;
+    }
+
+    epsilon
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-5;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    momentumPredictor yes;
+    pMin            pMin [ 1 -1 -2 0 0 0 0 ] 1000;
+    pRefCell        0;
+    pRefValue       1e5;
+}
+
+relaxationFactors
+{
+    U               1;
+    h               1;
+    k               1;
+    epsilon         1;
+}
+
+
+// ************************************************************************* //