diff --git a/applications/test/Circulator/Test-Circulator.C b/applications/test/Circulator/Test-Circulator.C
index 19a21f6b79..1369ef6e59 100644
--- a/applications/test/Circulator/Test-Circulator.C
+++ b/applications/test/Circulator/Test-Circulator.C
@@ -38,94 +38,6 @@ Description
 using namespace Foam;
 
 
-// return
-//   0: no match
-//  +1: identical
-//  -1: same face, but different orientation
-label compare(const face& a, const face& b)
-{
-    // Basic rule: we assume that the sequence of labels in each list
-    // will be circular in the same order (but not necessarily in the
-    // same direction or from the same starting point).
-
-    // Trivial reject: faces are different size
-    label sizeA = a.size();
-    label sizeB = b.size();
-
-    if (sizeA != sizeB || sizeA == 0)
-    {
-        return 0;
-    }
-
-    const_circulator<face> aCirc(a);
-    const_circulator<face> bCirc(b);
-
-    // Rotate face b until its element matches the starting element of face a.
-    do
-    {
-        if (aCirc() == bCirc())
-        {
-            // Set bCirc fulcrum to its iterator and increment the iterators
-            bCirc.setFulcrumToIterator();
-            ++aCirc;
-            ++bCirc;
-
-            break;
-        }
-    } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
-
-    // Look forwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::CLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return 1;
-    }
-    else
-    {
-        // Reset the circulators back to their fulcrum
-        aCirc.setIteratorToFulcrum();
-        bCirc.setIteratorToFulcrum();
-        ++aCirc;
-        --bCirc;
-    }
-
-    // Look backwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::ANTICLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return -1;
-    }
-
-    return 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
@@ -184,40 +96,40 @@ int main(int argc, char *argv[])
 
     Info<< nl << nl << "Compare two faces: " << endl;
     face a(identity(5));
-    Info<< "Compare " << a << " and " << a << " Match = " << compare(a, a)
+    Info<< "Compare " << a << " and " << a << " Match = " << face::compare(a, a)
         << endl;
 
     face b(reverseList(a));
-    Info<< "Compare " << a << " and " << b << " Match = " << compare(a, b)
+    Info<< "Compare " << a << " and " << b << " Match = " << face::compare(a, b)
         << endl;
 
     face c(a);
     c[4] = 3;
-    Info<< "Compare " << a << " and " << c << " Match = " << compare(a, c)
+    Info<< "Compare " << a << " and " << c << " Match = " << face::compare(a, c)
         << endl;
 
     face d(rotateList(a, 2));
-    Info<< "Compare " << a << " and " << d << " Match = " << compare(a, d)
+    Info<< "Compare " << a << " and " << d << " Match = " << face::compare(a, d)
         << endl;
 
     face g(labelList(5, 1));
     face h(g);
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[0] = 2;
     h[3] = 2;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[4] = 3;
     h[4] = 3;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     face face1(identity(1));
     Info<< "Compare " << face1 << " and " << face1
-        << " Match = " << compare(face1, face1) << endl;
+        << " Match = " << face::compare(face1, face1) << endl;
 
     Info<< nl << nl << "Zero face" << nl << endl;
 
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
index 60f4cb2f43..340c33049f 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseDict
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
@@ -17,7 +17,7 @@ FoamFile
 collapseEdgesCoeffs
 {
     // Edges shorter than this absolute value will be merged
-    minimumEdgeLength   1e-8;
+    minimumEdgeLength   1e-6;
 
     // The maximum angle between two edges that share a point attached to
     // no other edges
@@ -25,7 +25,7 @@ collapseEdgesCoeffs
 
     // The amount that minimumEdgeLength will be reduced by for each
     // edge if that edge's collapse generates a poor quality face
-    reductionFactor                         0.5;
+    reductionFactor     0.5;
 }
 
 
@@ -67,14 +67,17 @@ meshQualityCoeffs
 {
     // Name of the dictionary that has the mesh quality coefficients used
     // by motionSmoother::checkMesh
-    meshQualityCoeffDict    meshQualityDict;
+    #include                    "meshQualityDict";
+
+    // Maximum number of smoothing iterations for the reductionFactors
+    maximumSmoothingIterations  2;
 
     // Maximum number of outer iterations is mesh quality checking is enabled
-    maximumIterations       10;
+    maximumIterations           10;
 
     // Maximum number of iterations deletion of a point can cause a bad face
     // to be constructed before it is forced to not be deleted
-    maxPointErrorCount      5;
+    maxPointErrorCount          5;
 }
 
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
index f1db60d53f..fe1283a73a 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    /* -DFULLDEBUG -g -O0 */ \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
index 7979fa966a..1b83599b31 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
@@ -26,12 +26,6 @@ License
 #include "patchToPoly2DMesh.H"
 #include "PatchTools.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::patchToPoly2DMesh::flipFaceOrder()
@@ -275,12 +269,9 @@ void Foam::patchToPoly2DMesh::createPolyMeshComponents()
     addPatchFacesToOwner();
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//- Construct from a primitivePatch
 Foam::patchToPoly2DMesh::patchToPoly2DMesh
 (
     const MeshedSurface<face>& patch,
@@ -321,13 +312,5 @@ void Foam::patchToPoly2DMesh::createMesh()
     }
 }
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
index ea27dc6f53..b7c81e68a2 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
@@ -51,6 +51,7 @@ class patchToPoly2DMesh
 {
     // Private data
 
+        // Reference to the meshed surface
         const MeshedSurface<face>& patch_;
 
         const wordList& patchNames_;
@@ -62,11 +63,16 @@ class patchToPoly2DMesh
         const EdgeMap<label>& mapEdgesRegion_;
 
         pointField points_;
+
         faceList faces_;
+
         labelList owner_;
+
         labelList neighbour_;
 
-        //- Description of data_
+
+    // Private Member Functions
+
         void flipFaceOrder();
 
         void createNeighbours();
@@ -79,9 +85,6 @@ class patchToPoly2DMesh
 
         void createPolyMeshComponents();
 
-
-    // Private Member Functions
-
         //- Disallow default bitwise copy construct
         patchToPoly2DMesh(const patchToPoly2DMesh&);
 
@@ -110,47 +113,47 @@ public:
     // Member Functions
 
         // Access
-        pointField& points()
-        {
-            return points_;
-        }
 
-        faceList& faces()
-        {
-            return faces_;
-        }
+            pointField& points()
+            {
+                return points_;
+            }
 
-        labelList& owner()
-        {
-            return owner_;
-        }
+            faceList& faces()
+            {
+                return faces_;
+            }
 
-        labelList& neighbour()
-        {
-            return neighbour_;
-        }
+            labelList& owner()
+            {
+                return owner_;
+            }
 
-        const wordList& patchNames() const
-        {
-            return patchNames_;
-        }
+            labelList& neighbour()
+            {
+                return neighbour_;
+            }
 
-        const labelList& patchSizes() const
-        {
-            return patchSizes_;
-        }
+            const wordList& patchNames() const
+            {
+                return patchNames_;
+            }
+
+            const labelList& patchSizes() const
+            {
+                return patchSizes_;
+            }
+
+            const labelList& patchStarts() const
+            {
+                return patchStarts_;
+            }
 
-        const labelList& patchStarts() const
-        {
-            return patchStarts_;
-        }
-        // Check
 
         // Edit
-        void createMesh();
-
-        // Write
 
+            //- Create the mesh
+            void createMesh();
 };
 
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
index 088d0594f8..0b04d91b72 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
@@ -317,12 +317,12 @@ int main(int argc, char *argv[])
     }
 
     // Take over refinement levels and write to new time directory.
-    Pout<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
+    Info<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
         << nl << endl;
 
     mesh().write();
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
index f5186c3a25..7aa4459e89 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
@@ -2,31 +2,31 @@
 #define CGAL_PSURF_RINGS_H_
 
 // This file adapted from
-//     CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h
-// Licensed under CGAL-3.7/LICENSE.FREE_USE
+//     CGAL-4.0/examples/Jet_fitting_3/PolyhedralSurf_rings.h
+// Licensed under CGAL-4.0/LICENSE.FREE_USE
 
 // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
-// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
-// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
-// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
-// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
-// (Israel).  All rights reserved.
-
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is
-// furnished to do so, subject to the following conditions:
-
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
-// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-// SOFTWARE.
-
+// Utrecht University (The Netherlands),
+// ETH Zurich (Switzerland),
+// INRIA Sophia-Antipolis (France),
+// Max-Planck-Institute Saarbruecken (Germany),
+// and Tel-Aviv University (Israel).  All rights reserved.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 #include <cassert>
 
diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
index a50d5e3995..7309af83e1 100644
--- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
+++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
@@ -25,7 +25,8 @@ Class
     Foam::dynamicIndexedOctree
 
 Description
-    Non-pointer based hierarchical recursive searching
+    Non-pointer based hierarchical recursive searching.
+    Storage is dynamic, so elements can be deleted.
 
 SourceFiles
     dynamicIndexedOctree.C
@@ -451,6 +452,7 @@ public:
             );
         }
 
+
     // Member Functions
 
         // Access
@@ -656,6 +658,7 @@ public:
 
             label removeIndex(const label nodIndex, const label index);
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or
@@ -671,6 +674,7 @@ public:
 
             void writeTreeInfo() const;
 
+
     // IOstream Operators
 
         friend Ostream& operator<< <Type>
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 1d7ee06c2e..fcf273cc35 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -337,11 +337,11 @@ int Foam::face::compare(const face& a, const face& b)
     } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
 
     // If the circulator has stopped then faces a and b do not share a matching
-    // point
-    if (!bCirc.circulate())
-    {
-        return 0;
-    }
+    // point. Doesn't work on matching, single element face.
+    //if (!bCirc.circulate())
+    //{
+    //    return 0;
+    //}
 
     // Look forwards around the faces for a match
     do
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
index c44f646025..9b184f994e 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(polyMeshFilter, 0);
+}
+
+
 Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
 {
     polyTopoChange originalMeshToNewMesh(mesh);
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
index 135a52520c..28aaf7ea06 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
@@ -25,6 +25,10 @@ Class
     Foam::polyMeshFilter
 
 Description
+    Filter the edges and faces of a polyMesh whilst satisfying the given mesh
+    quality criteria.
+
+    Works on a copy of the mesh.
 
 SourceFiles
     polyMeshFilter.C
@@ -102,10 +106,10 @@ class polyMeshFilter
         //  faces
         const scalar faceReductionFactor_;
 
-        //-
+        //- Maximum number of times a deleted point can be associated with the
+        //  creation of a bad face it is forced to be kept.
         const label maxPointErrorCount_;
 
-
         //- The minimum edge length for each edge
         scalarField minEdgeLen_;
 
@@ -195,6 +199,10 @@ class polyMeshFilter
 
 public:
 
+    //- Runtime type information
+    ClassName("polyMeshFilter");
+
+
     // Constructors
 
         //- Construct from fvMesh
@@ -225,6 +233,7 @@ public:
             //- Filter edges only.
             label filterEdges(const label nOriginalBadFaces);
 
+            //- Filter all faces that are in the face zone indirectPatchFaces
             label filterIndirectPatchFaces();
 };
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index af911c7ad7..f6c35df41e 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(edgeCollapser, 0);
+}
+
+
 Foam::label Foam::edgeCollapser::longestEdge
 (
     const face& f,
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
index 02359f39ac..dc6b1578c1 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
@@ -61,7 +61,7 @@ class face;
 class edge;
 
 /*---------------------------------------------------------------------------*\
-                           Class edgeCollapser Declaration
+                        Class edgeCollapser Declaration
 \*---------------------------------------------------------------------------*/
 
 class edgeCollapser
@@ -84,12 +84,18 @@ private:
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Controls collapse of a face to an edge
         const scalar guardFraction_;
 
+        //- Only collapse face to a point if high aspect ratio
         const scalar maxCollapseFaceToPointSideLengthCoeff_;
 
+        //- Allow a face to be collapsed to a point early, before the test
+        //  to collapse to an edge
         const Switch allowEarlyCollapseToPoint_;
 
+        //- Fraction of maxCollapseFaceToPointSideLengthCoeff_ to use when
+        //  allowEarlyCollapseToPoint_ is on
         const scalar allowEarlyCollapseCoeff_;
 
 
@@ -266,8 +272,8 @@ public:
                 const dictionary& meshQualityDict
             );
 
-            // Check mesh and mark points on faces in error
-            // Returns boolList with points in error set
+            //- Check mesh and mark points on faces in error
+            //  Returns boolList with points in error set
             static label checkMeshQuality
             (
                 const polyMesh& mesh,
@@ -300,7 +306,7 @@ public:
                 polyTopoChange& meshMod
             ) const;
 
-            // Mark (in collapseEdge) any edges to collapse
+            //- Mark (in collapseEdge) any edges to collapse
             label markSmallEdges
             (
                 const scalarField& minEdgeLen,
@@ -309,7 +315,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
-            // Mark (in collapseEdge) any edges to merge
+            //- Mark (in collapseEdge) any edges to merge
             label markMergeEdges
             (
                 const scalar maxCos,
@@ -332,6 +338,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
+            //- Marks edges in the faceZone indirectPatchFaces for collapse
             void markIndirectPatchFaces
             (
                 PackedBoolList& collapseEdge,
diff --git a/src/fvOptions/fvOptions/fvOption.C b/src/fvOptions/fvOptions/fvOption.C
index 65c36598c5..6cd5c770f4 100644
--- a/src/fvOptions/fvOptions/fvOption.C
+++ b/src/fvOptions/fvOptions/fvOption.C
@@ -86,9 +86,9 @@ void Foam::fv::option::setSelection(const dictionary& dict)
         }
         case smMapRegion:
         {
-            dict_.lookup("nbrModelName") >> nbrModelName_;
-            dict_.lookup("nbrRegionName") >> nbrRegionName_;
-            master_  = readBool(dict_.lookup("master"));
+            dict.lookup("nbrModelName") >> nbrModelName_;
+            dict.lookup("nbrRegionName") >> nbrRegionName_;
+            master_ = readBool(dict.lookup("master"));
             break;
         }
         case smAll:
@@ -131,7 +131,6 @@ void Foam::fv::option::setCellSet()
                     WarningIn("option::setCellSet()")
                         << "Unable to find owner cell for point " << points_[i]
                         << endl;
-
                 }
 
             }
@@ -167,40 +166,47 @@ void Foam::fv::option::setCellSet()
         }
         case smMapRegion:
         {
-            if (active_)
+            if (active_ && master_)
             {
                 IInfo<< "- selecting inter region mapping" << endl;
 
                 const fvMesh& nbrMesh =
                     mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
-                boundBox BB(mesh_.points(), false);
-                boundBox nbrBB(nbrMesh.points(), false);
-
-                if (nbrBB.overlaps(BB))
+                if (mesh_.name() == nbrMesh.name())
                 {
-                    // Dummy patches
-                    wordList cuttingPatches;
-                    HashTable<word> patchMap;
+                    FatalErrorIn("option::setCellIds()")
+                        << "Inter-region model selected, but local and "
+                        << "neighbour regions are the same: " << nl
+                        << "    local region: " << mesh_.name() << nl
+                        << "    secondary region: " << nbrMesh.name() << nl
+                        << exit(FatalError);
+                }
 
-                    secondaryToPrimaryInterpPtr_.reset
+                if (mesh_.bounds().overlaps(nbrMesh.bounds()))
+                {
+                    meshInterpPtr_.reset
                     (
-                        new meshToMesh
+                        new meshToMeshNew
                         (
-                            nbrMesh,
                             mesh_,
-                            patchMap,
-                            cuttingPatches
+                            nbrMesh,
+                            meshToMeshNew::interpolationMethodNames_.read
+                            (
+                                dict_.lookup("interpolationMethod")
+                            )
                         )
                     );
                 }
                 else
                 {
                     FatalErrorIn("option::setCellSet()")
-                        << "regions do not overlap " << nbrMesh.name()
-                        << " in region " << mesh_.name() << nl
+                        << "regions " << mesh_.name() << " and "
+                        << nbrMesh.name() <<  " do not intersect"
                         << exit(FatalError);
                 }
+
+                V_ = meshInterpPtr_->V();
             }
             break;
         }
@@ -257,7 +263,7 @@ Foam::fv::option::option
     selectionMode_(selectionModeTypeNames_.read(dict_.lookup("selectionMode"))),
     cellSetName_("none"),
     V_(0.0),
-    secondaryToPrimaryInterpPtr_(),
+    meshInterpPtr_(),
     nbrModelName_("none"),
     nbrRegionName_("none"),
     master_(false),
@@ -296,7 +302,7 @@ Foam::autoPtr<Foam::fv::option> Foam::fv::option::New
 {
     word modelType(coeffs.lookup("type"));
 
-    IInfo<< "Selecting source model type " << modelType << endl;
+    IInfo<< "Selecting finite volume options model type " << modelType << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);
@@ -317,12 +323,7 @@ Foam::autoPtr<Foam::fv::option> Foam::fv::option::New
 
 
 Foam::fv::option::~option()
-{
-    if (!secondaryToPrimaryInterpPtr_.empty())
-    {
-        secondaryToPrimaryInterpPtr_.clear();
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
diff --git a/src/fvOptions/fvOptions/fvOption.H b/src/fvOptions/fvOptions/fvOption.H
index 8e51e4e3ea..3501fdfa45 100644
--- a/src/fvOptions/fvOptions/fvOption.H
+++ b/src/fvOptions/fvOptions/fvOption.H
@@ -35,7 +35,7 @@ Description
         selectionMode   cellSet;                // cellSet // points //cellZone
                                                 // mapRegion
 Note:
-    On evaluation, source expects to be added to the rhs of the equation
+    On evaluation, source/sink options are to be added to the equation rhs
 
 SourceFiles
     fvOption.C
@@ -50,7 +50,7 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "cellSet.H"
 #include "autoPtr.H"
-#include "meshToMesh.H"
+#include "meshToMeshNew.H"
 
 #include "runTimeSelectionTables.H"
 
@@ -131,8 +131,8 @@ protected:
 
         // Data for smMapRegion only
 
-            //- Mesh to mesh mapping for map option
-            autoPtr<meshToMesh> secondaryToPrimaryInterpPtr_;
+            //- Mesh to mesh interpolation object
+            autoPtr<meshToMeshNew> meshInterpPtr_;
 
             //- Name of the model in the neighbour mesh
             word nbrModelName_;
@@ -294,9 +294,8 @@ public:
             //- Return const access to the neighbour region name
             inline const word& nbrRegionName() const;
 
-            //- Return const access to the mapToMap Ptr
-            inline const autoPtr<meshToMesh>
-                secondaryToPrimaryInterpPtr() const;
+            //- Return const access to the mapToMap pointer
+            inline const meshToMeshNew& meshInterp() const;
 
             //- Return const access to the cell set
             inline const labelList& cells() const;
diff --git a/src/fvOptions/fvOptions/fvOptionI.H b/src/fvOptions/fvOptions/fvOptionI.H
index 8c3bf59233..aed4c14e55 100644
--- a/src/fvOptions/fvOptions/fvOptionI.H
+++ b/src/fvOptions/fvOptions/fvOptionI.H
@@ -138,10 +138,16 @@ inline const Foam::word& Foam::fv::option::nbrRegionName() const
 }
 
 
-inline const Foam::autoPtr<Foam::meshToMesh>
-Foam::fv::option::secondaryToPrimaryInterpPtr() const
+inline const Foam::meshToMeshNew& Foam::fv::option::meshInterp() const
 {
-    return secondaryToPrimaryInterpPtr_;
+    if (!meshInterpPtr_.valid())
+    {
+        FatalErrorIn("const meshToMeshNew& meshInterp() const")
+            << "Interpolation object not set"
+            << abort(FatalError);
+    }
+
+    return meshInterpPtr_();
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C b/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
index e1eb4e2d4a..1cafec08c4 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
@@ -57,7 +57,7 @@ Foam::fv::constantHeatTransfer::constantHeatTransfer
     htcConst_(),
     AoV_()
 {
-    if (master_)
+    if (active() && master_)
     {
         htcConst_.reset
         (
@@ -91,19 +91,7 @@ Foam::fv::constantHeatTransfer::constantHeatTransfer
             )
         );
 
-        const DimensionedField<scalar, volMesh>& htcConsti =
-            htcConst_().dimensionedInternalField();
-        const DimensionedField<scalar, volMesh>& AoVi =
-            AoV_().dimensionedInternalField();
-        dimensionedScalar interVol
-        (
-            "V",
-            dimVolume,
-            secondaryToPrimaryInterpPtr_->V()
-        );
-
-        htc_.dimensionedInternalField() = htcConsti*AoVi*interVol/mesh.V();
-        htc_.correctBoundaryConditions();
+        htc_ = htcConst_()*AoV_();
     }
 }
 
@@ -116,10 +104,9 @@ Foam::fv::constantHeatTransfer::~constantHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::constantHeatTransfer::calculateHtc()
+void Foam::fv::constantHeatTransfer::calculateHtc()
 {
-    return htc_;
+    // do nothing
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H b/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
index 52facea986..ceb4aaf15a 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
@@ -87,7 +87,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
index 3a1067bdd7..31a684afb9 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
@@ -25,7 +25,7 @@ License
 
 #include "interRegionHeatTransferModel.H"
 #include "fluidThermo.H"
-#include "fvm.H"
+#include "fvmSup.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fvcVolumeIntegrate.H"
 #include "fvOptionList.H"
@@ -41,10 +41,15 @@ namespace fv
 }
 
 
-// * * * * * * * * * * * *  Private member functions * * * * * * * * * * * //
+// * * * * * * * * * * * *  Protected member functions * * * * * * * * * * * //
 
-void Foam::fv::interRegionHeatTransferModel::check()
+void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 {
+    if (!firstIter_)
+    {
+        return;
+    }
+
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
     const optionList& fvOptions = nbrMesh.lookupObject<optionList>("fvOptions");
@@ -66,11 +71,31 @@ void Foam::fv::interRegionHeatTransferModel::check()
 
     if (!nbrModelFound)
     {
-        FatalErrorIn("interRegionHeatTransferModel::check()")
-            << "Secondary source name not found" << nbrModelName_
+        FatalErrorIn("interRegionHeatTransferModel::setNbrModel()")
+            << "Neighbour model not found" << nbrModelName_
             << " in region " << nbrMesh.name() << nl
             << exit(FatalError);
     }
+
+    firstIter_ = false;
+}
+
+
+void Foam::fv::interRegionHeatTransferModel::correct()
+{
+    if (master_)
+    {
+        if (mesh_.time().timeIndex() != timeIndex_)
+        {
+            calculateHtc();
+            timeIndex_ = mesh_.time().timeIndex();
+        }
+    }
+    else
+    {
+        nbrModel().correct();
+        interpolate(nbrModel().htc(), htc_);
+    }
 }
 
 
@@ -85,8 +110,9 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
 )
 :
     option(name, modelType, dict, mesh),
-    nbrModel_(),
+    nbrModel_(NULL),
     firstIter_(true),
+    timeIndex_(-1),
     htc_
     (
         IOobject
@@ -106,7 +132,9 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
         ),
         zeroGradientFvPatchScalarField::typeName
     ),
-    semiImplicit_(false)
+    semiImplicit_(false),
+    TName_(coeffs_.lookupOrDefault<word>("TName", "T")),
+    TNbrName_(coeffs_.lookupOrDefault<word>("TNbrName", "T"))
 {
     if (active())
     {
@@ -132,19 +160,14 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     const label fieldI
 )
 {
-    if (!secondaryToPrimaryInterpPtr_.valid())
-    {
-        return;
-    }
+    setNbrModel();
 
-    if (firstIter_)
-    {
-        check();
-        firstIter_ = false;
-    }
+    correct();
 
     const volScalarField& h = eqn.psi();
 
+    const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
+
     tmp<volScalarField> tTmapped
     (
         new volScalarField
@@ -157,8 +180,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            mesh_,
-            dimensionedScalar("T", dimTemperature, 0.0)
+            T
         )
     );
 
@@ -166,26 +188,10 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
-    const volScalarField& Tnbr = nbrMesh.lookupObject<volScalarField>("T");
+    const volScalarField& Tnbr =
+        nbrMesh.lookupObject<volScalarField>(TNbrName_);
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        Tmapped,
-        Tnbr,
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
-
-    if (!master_)
-    {
-        secondaryToPrimaryInterpPtr_->interpolateInternalField
-        (
-            htc_,
-            nbrModel_->calculateHtc(),
-            meshToMesh::CELL_VOLUME_WEIGHT,
-            eqOp<scalar>()
-        );
-    }
+    interpolate(Tnbr, Tmapped.internalField());
 
     if (debug)
     {
@@ -226,7 +232,6 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     }
     else
     {
-        const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
         eqn += htc_*(Tmapped - T);
     }
 }
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
index 2ced33cacc..736bd094fb 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
@@ -26,7 +26,12 @@ Class
 
 Description
     Base class for inter region heat exchange. The derived classes must
-    provide the heat transfer coefficient (htc)
+    provide the heat transfer coeffisine (htc) which is used as follows
+    in the energy equation.
+
+     \f[
+        -htc*Tmapped + Sp(htc, T)
+     \f]
 
 \*---------------------------------------------------------------------------*/
 
@@ -63,24 +68,67 @@ private:
         //- First iteration
         bool firstIter_;
 
-
-    // Private members
-
-        //- Check coupled interRegionHeatTransferModel
-        void check();
+        //- Time index - used for updating htc
+        label timeIndex_;
 
 
 protected:
 
     // Protected data
 
-        //- Heat transfer coefficient [W/m2/k] by area/volume [1/m]
-        // registered on the master mesh
+        //- Heat transfer coefficient [W/m2/k] times area/volume [1/m]
         volScalarField htc_;
 
         //- Flag to activate semi-implicit coupling
         bool semiImplicit_;
 
+        //- Name of temperature field; default = "T"
+        word TName_;
+
+        //- Name of neighbour temperature field; default = "T"
+        word TNbrName_;
+
+
+    // Protected member functions
+
+        //- Set the neighbour interRegionHeatTransferModel
+        void setNbrModel();
+
+        //- Correct to calculate the inter-region heat transfer coefficient
+        void correct();
+
+        //- Interpolate field with nbrModel specified
+        template<class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const interRegionHeatTransferModel& nbrModel,
+            const Field<Type>& field
+        ) const;
+
+        //- Interpolate field without nbrModel specified
+        template<class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const Field<Type>& field
+        ) const;
+
+        //- Interpolate field with nbrModel specified
+        template<class Type>
+        void interpolate
+        (
+            const interRegionHeatTransferModel& nbrModel,
+            const Field<Type>& field,
+            Field<Type>& result
+        ) const;
+
+        //- Interpolate field without nbrModel specified
+        template<class Type>
+        void interpolate
+        (
+            const Field<Type>& field,
+            Field<Type>& result
+        ) const;
+
 
 public:
 
@@ -110,16 +158,19 @@ public:
         // Access
 
             //- Return the heat transfer coefficient
-            const volScalarField& htc() const
-            {
-                return  htc_;
-            }
+            inline const volScalarField& htc() const;
+
+            //- Return const access to the neighbour model
+            inline const interRegionHeatTransferModel& nbrModel() const;
+
+            //- Return access to the neighbour model
+            inline interRegionHeatTransferModel& nbrModel();
 
             //-Source term to fvMatrix<scalar>
             virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
 
             //- Calculate heat transfer coefficient
-            virtual const tmp<volScalarField> calculateHtc() = 0;
+            virtual void calculateHtc() = 0;
 
 
         // I-O
@@ -138,6 +189,17 @@ public:
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interRegionHeatTransferModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "interRegionHeatTransferModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H
new file mode 100644
index 0000000000..6a6cbe550b
--- /dev/null
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const Foam::volScalarField&
+Foam::fv::interRegionHeatTransferModel::htc() const
+{
+    return  htc_;
+}
+
+
+inline const Foam::fv::interRegionHeatTransferModel&
+Foam::fv::interRegionHeatTransferModel::nbrModel() const
+{
+    if (nbrModel_ == NULL)
+    {
+        FatalErrorIn("const interRegionHeatTransferModel& nbrModel() const")
+            << "Neighbour model not set"
+            << abort(FatalError);
+    }
+
+    return *nbrModel_;
+}
+
+
+inline Foam::fv::interRegionHeatTransferModel&
+Foam::fv::interRegionHeatTransferModel::nbrModel()
+{
+    if (nbrModel_ == NULL)
+    {
+        FatalErrorIn("interRegionHeatTransferModel& nbrModel()")
+            << "Neighbour model not set"
+            << abort(FatalError);
+    }
+
+    return *nbrModel_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C
new file mode 100644
index 0000000000..4a5363a2d1
--- /dev/null
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const interRegionHeatTransferModel& nbrModel,
+    const Field<Type>& field
+) const
+{
+    if (master_)
+    {
+        return meshInterp().mapTgtToSrc(field);
+    }
+    else
+    {
+        return (nbrModel.meshInterp().mapSrcToTgt(field));
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const Field<Type>& field
+) const
+{
+    return interpolate(nbrModel(), field);
+}
+
+
+template<class Type>
+void Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const interRegionHeatTransferModel& nbrModel,
+    const Field<Type>& field,
+    Field<Type>& result
+) const
+{
+    if (master_)
+    {
+        meshInterp().mapTgtToSrc(field, plusEqOp<scalar>(), result);
+    }
+    else
+    {
+        nbrModel.meshInterp().mapSrcToTgt(field, plusEqOp<scalar>(), result);
+    }
+}
+
+
+template<class Type>
+void Foam::fv::interRegionHeatTransferModel::interpolate
+(
+    const Field<Type>& field,
+    Field<Type>& result
+) const
+{
+    return interpolate(nbrModel(), field, result);
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C b/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
index c49001e388..fd49c31963 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
@@ -94,6 +94,7 @@ Foam::fv::tabulatedHeatTransfer::tabulatedHeatTransfer
 :
     interRegionHeatTransferModel(name, modelType, dict, mesh),
     UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    UNbrName_(coeffs_.lookupOrDefault<word>("UNbrName", "U")),
     hTable_(),
     AoV_(),
     startTimeName_(mesh.time().timeName())
@@ -108,22 +109,16 @@ Foam::fv::tabulatedHeatTransfer::~tabulatedHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::tabulatedHeatTransfer::calculateHtc()
+void Foam::fv::tabulatedHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
 
-    const volVectorField& UNbr = nbrMesh.lookupObject<volVectorField>(UName_);
+    const volVectorField& UNbr =
+        nbrMesh.lookupObject<volVectorField>(UNbrName_);
 
-    scalarField UMagNbrMapped(htc_.internalField().size(), 0.0);
+    const scalarField UMagNbr(mag(UNbr));
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        UMagNbrMapped,
-        mag(UNbr),
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
+    const scalarField UMagNbrMapped(interpolate(UMagNbr));
 
     const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
 
@@ -134,9 +129,7 @@ Foam::fv::tabulatedHeatTransfer::calculateHtc()
         htcc[i] = hTable()(mag(U[i]), UMagNbrMapped[i]);
     }
 
-    htcc = htcc*AoV()*secondaryToPrimaryInterpPtr_->V()/mesh_.V();
-
-    return htc_;
+    htcc = htcc*AoV();
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H b/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
index c46f747a14..8feb9d68df 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
@@ -60,6 +60,9 @@ private:
         //- Name of velocity field; default = U
         word UName_;
 
+        //- Name of neighbour velocity field; default = U
+        word UNbrName_;
+
         //- 2D look up table
         autoPtr<interpolation2DTable<scalar> > hTable_;
 
@@ -101,7 +104,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C b/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
index c675e95e94..cf96fa0433 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
@@ -55,7 +55,7 @@ Foam::fv::variableHeatTransfer::variableHeatTransfer
 )
 :
     interRegionHeatTransferModel(name, modelType, dict, mesh),
-    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    UNbrName_(coeffs_.lookupOrDefault<word>("UNbrName", "U")),
     a_(0),
     b_(0),
     c_(0),
@@ -97,8 +97,7 @@ Foam::fv::variableHeatTransfer::~variableHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::variableHeatTransfer::calculateHtc()
+void Foam::fv::variableHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh =
         mesh_.time().lookupObject<fvMesh>(nbrRegionName());
@@ -109,26 +108,18 @@ Foam::fv::variableHeatTransfer::calculateHtc()
     const fluidThermo& nbrThermo =
         nbrMesh.lookupObject<fluidThermo>("thermophysicalProperties");
 
-    const volVectorField& U = nbrMesh.lookupObject<volVectorField>(UName_);
+    const volVectorField& UNbr =
+        nbrMesh.lookupObject<volVectorField>(UNbrName_);
 
-    const volScalarField Re(mag(U)*ds_*nbrThermo.rho()/nbrTurb.mut());
+    const volScalarField ReNbr(mag(UNbr)*ds_*nbrThermo.rho()/nbrTurb.mut());
 
-    const volScalarField Nu(a_*pow(Re, b_)*pow(Pr_, c_));
+    const volScalarField NuNbr(a_*pow(ReNbr, b_)*pow(Pr_, c_));
 
-    scalarField htcNbrMapped(htc_.internalField().size(), 0.0);
+    const scalarField htcNbr(NuNbr*nbrTurb.kappaEff()/ds_);
 
-    secondaryToPrimaryInterpPtr_->interpolateInternalField
-    (
-        htcNbrMapped,
-        Nu*nbrTurb.kappaEff()/ds_,
-        meshToMesh::MAP,
-        eqOp<scalar>()
-    );
+    const scalarField htcNbrMapped(interpolate(htcNbr));
 
-    htc_.internalField() =
-        htcNbrMapped*AoV_*secondaryToPrimaryInterpPtr_->V()/mesh_.V();
-
-    return htc_;
+    htc_.internalField() = htcNbrMapped*AoV_;
 }
 
 
@@ -150,7 +141,7 @@ bool Foam::fv::variableHeatTransfer::read(const dictionary& dict)
 {
     if (option::read(dict))
     {
-        coeffs_.readIfPresent("UName", UName_);
+        coeffs_.readIfPresent("UNbrName", UNbrName_);
 
         coeffs_.readIfPresent("a", a_);
         coeffs_.readIfPresent("b", b_);
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H b/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
index cfe0aaffae..57e42be811 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
+++ b/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
@@ -65,8 +65,8 @@ private:
 
     // Private data
 
-        //- Name of velocity field; default = U
-        word UName_;
+        //- Name of neighbour velocity field; default = U
+        word UNbrName_;
 
         //- Model constants
         scalar a_;
@@ -108,7 +108,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index d159b0df0b..401d14a383 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -201,4 +201,6 @@ regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledWallGAMGIn
 regionCoupled/GAMG/interfaceFields/regionCoupledGAMGInterfaceField/regionCoupledGAMGInterfaceField.C
 regionCoupled/GAMG/interfaceFields/regionCoupledGAMGInterfaceField/regionCoupledWallGAMGInterfaceField.C
 
+tetOverlapVolume/tetOverlapVolume.C
+
 LIB = $(FOAM_LIBBIN)/libmeshTools
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
similarity index 67%
rename from src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
rename to src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index 219d78ce76..c0f192f7de 100644
--- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,6 +130,135 @@ Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 
 // * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::tetOverlapVolume::cellCellOverlapMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB,
+    const scalar threshold
+) const
+{
+    const cell& cFacesA = meshA.cells()[cellAI];
+    const point& ccA = meshA.cellCentres()[cellAI];
+
+    const cell& cFacesB = meshB.cells()[cellBI];
+    const point& ccB = meshB.cellCentres()[cellBI];
+
+    scalar vol = 0.0;
+
+    forAll(cFacesA, cFA)
+    {
+        label faceAI = cFacesA[cFA];
+
+        const face& fA = meshA.faces()[faceAI];
+        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+        if (!pyrA.overlaps(cellBbB))
+        {
+            continue;
+        }
+
+        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+        label tetBasePtAI = 0;
+
+        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+        {
+            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+            label otherFacePtAI = fA.fcIndex(facePtAI);
+
+            label pt0I = -1;
+            label pt1I = -1;
+
+            if (ownA)
+            {
+                pt0I = fA[facePtAI];
+                pt1I = fA[otherFacePtAI];
+            }
+            else
+            {
+                pt0I = fA[otherFacePtAI];
+                pt1I = fA[facePtAI];
+            }
+
+            const tetPoints tetA
+            (
+                ccA,
+                tetBasePtA,
+                meshA.points()[pt0I],
+                meshA.points()[pt1I]
+            );
+            const treeBoundBox tetABb(tetA.bounds());
+
+
+            // Loop over tets of cellB
+            forAll(cFacesB, cFB)
+            {
+                label faceBI = cFacesB[cFB];
+
+                const face& fB = meshB.faces()[faceBI];
+                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+                if (!pyrB.overlaps(pyrA))
+                {
+                    continue;
+                }
+
+                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+                label tetBasePtBI = 0;
+
+                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+                {
+                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+                    label otherFacePtBI = fB.fcIndex(facePtBI);
+
+                    label pt0I = -1;
+                    label pt1I = -1;
+
+                    if (ownB)
+                    {
+                        pt0I = fB[facePtBI];
+                        pt1I = fB[otherFacePtBI];
+                    }
+                    else
+                    {
+                        pt0I = fB[otherFacePtBI];
+                        pt1I = fB[facePtBI];
+                    }
+
+                    const tetPoints tetB
+                    (
+                        ccB,
+                        tetBasePtB,
+                        meshB.points()[pt0I],
+                        meshB.points()[pt1I]
+                    );
+
+                    if (!tetB.bounds().overlaps(tetABb))
+                    {
+                        continue;
+                    }
+
+                    vol += tetTetOverlapVol(tetA, tetB);
+
+                    if (vol > threshold)
+                    {
+                        return true;
+                    }
+                }
+            }
+        }
+    }
+
+    return false;
+}
+
+
 Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
 (
     const primitiveMesh& meshA,
diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
similarity index 82%
rename from src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
rename to src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index e8890a3ee9..58646f2e5a 100644
--- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Class
     Foam::tetOverlapVolume
 
 Description
-    Calculates overlap volume of two tets.
+    Calculates the overlap volume of two cells using tetrahedral decomposition
 
 SourceFiles
     tetOverlapVolume.C
@@ -48,14 +48,14 @@ class polyMesh;
 class tetPoints;
 
 /*---------------------------------------------------------------------------*\
-                         Class tetOverlapVolume Declaration
+                      Class tetOverlapVolume Declaration
 \*---------------------------------------------------------------------------*/
 
 class tetOverlapVolume
 {
     // Private member functions
 
-        //- Tet Overlap Vol
+        //- Tet overlap volume
         scalar tetTetOverlapVol
         (
             const tetPoints& tetA,
@@ -94,6 +94,16 @@ public:
             const label cellBI
         ) const;
 
+        //- Return true if olverlap volume is greater than threshold
+        bool cellCellOverlapMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB,
+            const scalar threshold = 0.0
+        ) const;
 
         //- Calculates the overlap volume
         scalar cellCellOverlapVolumeMinDecomp
@@ -112,7 +122,6 @@ public:
 
 } // End namespace Foam
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
index a77459f82f..f8b689757d 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,6 +77,8 @@ Description
                 mean            on;
                 prime2Mean      on;
                 base            time;
+                window          10.0;
+                windowName      w1;
             }
             p
             {
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
index 0fdf344531..f4b03dd4b4 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,15 @@ void Foam::fieldAverage::addMeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_MEAN;
+        word meanFieldName = fieldName + EXT_MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
 
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
@@ -100,7 +108,16 @@ void Foam::fieldAverage::addPrime2MeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
+
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
         if (obr_.foundObject<fieldType2>(meanFieldName))
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
index 39bbfb4fa9..3b01b21594 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,8 @@ Foam::fieldAverageItem::fieldAverageItem()
     mean_(0),
     prime2Mean_(0),
     base_(ITER),
-    window_(-1.0)
+    window_(-1.0),
+    windowName_("")
 {}
 
 
@@ -64,7 +65,8 @@ Foam::fieldAverageItem::fieldAverageItem(const fieldAverageItem& faItem)
     mean_(faItem.mean_),
     prime2Mean_(faItem.prime2Mean_),
     base_(faItem.base_),
-    window_(faItem.window_)
+    window_(faItem.window_),
+    windowName_(faItem.windowName_)
 {}
 
 
@@ -94,6 +96,7 @@ void Foam::fieldAverageItem::operator=(const fieldAverageItem& rhs)
     prime2Mean_ = rhs.prime2Mean_;
     base_ = rhs.base_;
     window_ = rhs.window_;
+    windowName_ = rhs.windowName_;
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
index 992b5e393f..7fe21a9aef 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
         prime2Mean      on;
         base            time; // iteration
         window          200;  // optional averaging window
+        windowName      w1;   // optional window name (default = "")
     }
     \endverbatim
 
@@ -107,6 +108,9 @@ private:
         //- Averaging window - defaults to -1 for 'all iters/time'
         scalar window_;
 
+        //- Averaging window name - defaults to 'window'
+        word windowName_;
+
 
 public:
 
@@ -171,6 +175,11 @@ public:
                 return window_;
             }
 
+            const word& windowName() const
+            {
+                return windowName_;
+            }
+
 
     // Member Operators
 
@@ -190,7 +199,8 @@ public:
              && a.mean_ == b.mean_
              && a.prime2Mean_ == b.prime2Mean_
              && a.base_ == b.base_
-             && a.window_ == b.window_;
+             && a.window_ == b.window_
+             && a.windowName_ == b.windowName_;
         }
 
         friend bool operator!=
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
index 4ebad485a7..af293ff989 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,6 +46,7 @@ Foam::fieldAverageItem::fieldAverageItem(Istream& is)
     entry.lookup("prime2Mean") >> prime2Mean_;
     base_ = baseTypeNames_[entry.lookup("base")];
     window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    windowName_ = entry.lookupOrDefault<word>("windowName", "");
 }
 
 
@@ -66,6 +67,7 @@ Foam::Istream& Foam::operator>>(Istream& is, fieldAverageItem& faItem)
     entry.lookup("prime2Mean") >> faItem.prime2Mean_;
     faItem.base_ = faItem.baseTypeNames_[entry.lookup("base")];
     faItem.window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    faItem.windowName_ = entry.lookupOrDefault<word>("windowName", "");
 
     return is;
 }
@@ -90,6 +92,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fieldAverageItem& faItem)
     {
         os.writeKeyword("window") << faItem.window_
             << token::END_STATEMENT << nl;
+
+        if (faItem.windowName_ != "")
+        {
+            os.writeKeyword("windowName") << faItem.windowName_
+                << token::END_STATEMENT << nl;
+        }
     }
 
     os  << token::END_BLOCK << nl;
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index ef0415f0f2..747097a781 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -59,8 +59,8 @@ $(meshToMesh)/meshToMesh.C
 $(meshToMesh)/calculateMeshToMeshAddressing.C
 $(meshToMesh)/calculateMeshToMeshWeights.C
 
-tetOverlapVolume = meshToMeshInterpolation/tetOverlapVolume
-$(tetOverlapVolume)/tetOverlapVolume.C
-
+meshToMeshNew = meshToMeshInterpolation/meshToMeshNew
+$(meshToMeshNew)/meshToMeshNew.C
+$(meshToMeshNew)/meshToMeshNewParallelOps.C
 
 LIB = $(FOAM_LIBBIN)/libsampling
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 9344f75ba0..b99c0629cc 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -27,6 +27,9 @@ Class
 Description
     mesh to mesh interpolation class.
 
+Note
+    This class is due to be deprecated in favour of meshToMeshNew
+
 SourceFiles
     meshToMesh.C
     calculateMeshToMeshAddressing.C
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
new file mode 100644
index 0000000000..9ecda6cb4f
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
@@ -0,0 +1,987 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshToMeshNew.H"
+#include "OFstream.H"
+#include "Time.H"
+#include "globalIndex.H"
+#include "mergePoints.H"
+#include "treeBoundBox.H"
+#include "tetOverlapVolume.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(meshToMeshNew, 0);
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::meshToMeshNew::interpolationMethod,
+        2
+    >::names[] =
+    {
+        "map",
+        "cellVolumeWeight"
+    };
+
+    const NamedEnum<meshToMeshNew::interpolationMethod, 2>
+        meshToMeshNew::interpolationMethodNames_;
+}
+
+Foam::scalar Foam::meshToMeshNew::tolerance_ = 1e-6;
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::meshToMeshNew::writeConnectivity
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const labelListList& srcToTargetAddr
+) const
+{
+    Pout<< "Source size = " << src.nCells() << endl;
+    Pout<< "Target size = " << tgt.nCells() << endl;
+
+    word fName("addressing_" + src.name() + "_to_" + tgt.name());
+
+    if (Pstream::parRun())
+    {
+        fName = fName +  "_proc" + Foam::name(Pstream::myProcNo());
+    }
+
+    OFstream os(src.time().path()/fName + ".obj");
+
+    label vertI = 0;
+    forAll(srcToTargetAddr, i)
+    {
+        const labelList& tgtAddress = srcToTargetAddr[i];
+        forAll(tgtAddress, j)
+        {
+            label tgtI = tgtAddress[j];
+            const vector& c0 = src.cellCentres()[i];
+
+            const cell& c = tgt.cells()[tgtI];
+            const pointField pts(c.points(tgt.faces(), tgt.points()));
+            forAll(pts, j)
+            {
+                const point& p = pts[j];
+                os  << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+                vertI++;
+                os  << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
+                    << nl;
+                vertI++;
+                os  << "l " << vertI - 1 << ' ' << vertI << nl;
+            }
+        }
+    }
+}
+
+
+Foam::labelList Foam::meshToMeshNew::maskCells
+(
+    const polyMesh& src,
+    const polyMesh& tgt
+) const
+{
+    boundBox intersectBb
+    (
+        max(src.bounds().min(), tgt.bounds().min()),
+        min(src.bounds().max(), tgt.bounds().max())
+    );
+
+    intersectBb.inflate(0.01);
+
+    const cellList& srcCells = src.cells();
+    const faceList& srcFaces = src.faces();
+    const pointField& srcPts = src.points();
+
+    DynamicList<label> cells(src.size());
+    forAll(srcCells, srcI)
+    {
+        boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
+        if (intersectBb.overlaps(cellBb))
+        {
+            cells.append(srcI);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "participating source mesh cells: " << cells.size() << endl;
+    }
+
+    return cells;
+}
+
+
+bool Foam::meshToMeshNew::findInitialSeeds
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const labelList& srcCellIDs,
+    const boolList& mapFlag,
+    const label startSeedI,
+    label& srcSeedI,
+    label& tgtSeedI
+) const
+{
+    const cellList& srcCells = src.cells();
+    const faceList& srcFaces = src.faces();
+    const pointField& srcPts = src.points();
+
+    for (label i = startSeedI; i < srcCellIDs.size(); i++)
+    {
+        label srcI = srcCellIDs[i];
+
+        if (mapFlag[srcI])
+        {
+            const pointField
+                pts(srcCells[srcI].points(srcFaces, srcPts).xfer());
+
+            forAll(pts, ptI)
+            {
+                const point& pt = pts[ptI];
+                label tgtI = tgt.cellTree().findInside(pt);
+
+                if (tgtI != -1 && intersect(src, tgt, srcI, tgtI))
+                {
+                    srcSeedI = srcI;
+                    tgtSeedI = tgtI;
+
+                    return true;
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "could not find starting seed" << endl;
+    }
+
+    return false;
+}
+
+
+void Foam::meshToMeshNew::appendToDirectSeeds
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    boolList& mapFlag,
+    labelList& srcTgtSeed,
+    DynamicList<label>& srcSeeds,
+    label& srcSeedI,
+    label& tgtSeedI
+) const
+{
+    const labelList& srcNbr = src.cellCells()[srcSeedI];
+    const labelList& tgtNbr = tgt.cellCells()[tgtSeedI];
+
+    const vectorField& srcCentre = src.cellCentres();
+
+    forAll(srcNbr, i)
+    {
+        label srcI = srcNbr[i];
+
+        if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1))
+        {
+            // source cell srcI not yet mapped
+
+            // identfy if target cell exists for source cell srcI
+            bool found = false;
+            forAll(tgtNbr, j)
+            {
+                label tgtI = tgtNbr[j];
+
+                if (tgt.pointInCell(srcCentre[srcI], tgtI))
+                {
+                    // new match - append to lists
+                    found = true;
+
+                    srcTgtSeed[srcI] = tgtI;
+                    srcSeeds.append(srcI);
+
+                    break;
+                }
+            }
+
+            if (!found)
+            {
+                // no match available for source cell srcI
+                mapFlag[srcI] = false;
+            }
+        }
+    }
+
+    if (srcSeeds.size())
+    {
+        srcSeedI = srcSeeds.remove();
+        tgtSeedI = srcTgtSeed[srcSeedI];
+    }
+    else
+    {
+        srcSeedI = -1;
+        tgtSeedI = -1;
+    }
+}
+
+
+void Foam::meshToMeshNew::calcDirect
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const label srcSeedI,
+    const label tgtSeedI
+)
+{
+    // store a list of src cells already mapped
+    boolList srcSeedFlag(src.nCells(), true);
+    labelList srcTgtSeed(src.nCells(), -1);
+
+    List<DynamicList<label> > srcToTgt(src.nCells());
+    List<DynamicList<label> > tgtToSrc(tgt.nCells());
+
+    DynamicList<label> srcSeeds;
+
+    const scalarField& srcVc = src.cellVolumes();
+
+    label srcCellI = srcSeedI;
+    label tgtCellI = tgtSeedI;
+
+    do
+    {
+        // store src/tgt cell pair
+        srcToTgt[srcCellI].append(tgtCellI);
+        tgtToSrc[tgtCellI].append(srcCellI);
+
+        // mark source cell srcSeedI as matched
+        srcSeedFlag[srcCellI] = false;
+
+        // accumulate intersection volume
+        V_ += srcVc[srcCellI];
+
+        // find new source seed cell
+        appendToDirectSeeds
+        (
+            src,
+            tgt,
+            srcSeedFlag,
+            srcTgtSeed,
+            srcSeeds,
+            srcCellI,
+            tgtCellI
+        );
+
+    }
+    while (srcCellI >= 0);
+
+    // transfer addressing into persistent storage
+    forAll(srcToTgtCellAddr_, i)
+    {
+        srcToTgtCellAddr_[i].transfer(srcToTgt[i]);
+        srcToTgtCellWght_[i] = scalarList(srcToTgtCellAddr_[i].size(), 1.0);
+    }
+
+    forAll(tgtToSrcCellAddr_, i)
+    {
+        tgtToSrcCellAddr_[i].transfer(tgtToSrc[i]);
+        tgtToSrcCellWght_[i] = scalarList(tgtToSrcCellAddr_[i].size(), 1.0);
+    }
+}
+
+
+void Foam::meshToMeshNew::normaliseWeights
+(
+    const word& descriptor,
+    const scalarField& cellVolumes,
+    const labelListList& addr,
+    scalarListList& wght
+) const
+{
+    const label nCell = returnReduce(wght.size(), sumOp<label>());
+
+    if (nCell > 0)
+    {
+        scalar minW = GREAT;
+        scalar maxW = -GREAT;
+
+        forAll(wght, cellI)
+        {
+            scalarList& w = wght[cellI];
+            scalar s = sum(w);
+            scalar Vc = cellVolumes[cellI];
+
+            forAll(w, i)
+            {
+                w[i] /= Vc;
+            }
+
+            minW = min(minW, s/Vc);
+            maxW = max(maxW, s/Vc);
+        }
+
+        Info<< type() << ": " << descriptor << " weights min/max = "
+            << returnReduce(minW, minOp<scalar>()) << ", "
+            << returnReduce(maxW, maxOp<scalar>()) << endl;
+    }
+}
+
+
+void Foam::meshToMeshNew::appendNbrTgtCells
+(
+    const label tgtCellI,
+    const polyMesh& tgt,
+    const DynamicList<label>& visitedTgtCells,
+    DynamicList<label>& nbrTgtCellIDs
+) const
+{
+    const labelList& nbrCells = tgt.cellCells()[tgtCellI];
+
+    // filter out cells already visited from cell neighbours
+    forAll(nbrCells, i)
+    {
+        label nbrCellI = nbrCells[i];
+
+        if
+        (
+            (findIndex(visitedTgtCells, nbrCellI) == -1)
+         && (findIndex(nbrTgtCellIDs, nbrCellI) == -1)
+        )
+        {
+            nbrTgtCellIDs.append(nbrCellI);
+        }
+    }
+}
+
+
+void Foam::meshToMeshNew::setNextCells
+(
+    label& startSeedI,
+    label& srcCellI,
+    label& tgtCellI,
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const labelList& srcCellIDs,
+    const boolList& mapFlag,
+    const DynamicList<label>& visitedCells,
+    labelList& seedCells
+) const
+{
+    const labelList& srcNbrCells = src.cellCells()[srcCellI];
+
+    // set possible seeds for later use by querying all src cell neighbours
+    // with all visited target cells
+    bool valuesSet = false;
+    forAll(srcNbrCells, i)
+    {
+        label cellS = srcNbrCells[i];
+
+        if (mapFlag[cellS] && seedCells[cellS] == -1)
+        {
+            forAll(visitedCells, j)
+            {
+                label cellT = visitedCells[j];
+
+                if (intersect(src, tgt, cellS, cellT))
+                {
+                    seedCells[cellS] = cellT;
+
+                    if (!valuesSet)
+                    {
+                        srcCellI = cellS;
+                        tgtCellI = cellT;
+                        valuesSet = true;
+                    }
+                }
+            }
+        }
+    }
+
+    // set next src and tgt cells if not set above
+    if (valuesSet)
+    {
+        return;
+    }
+    else
+    {
+        // try to use existing seed
+        bool foundNextSeed = false;
+        for (label i = startSeedI; i < srcCellIDs.size(); i++)
+        {
+            label cellS = srcCellIDs[i];
+
+            if (mapFlag[cellS])
+            {
+                if (!foundNextSeed)
+                {
+                    startSeedI = i;
+                    foundNextSeed = true;
+                }
+
+                if (seedCells[cellS] != -1)
+                {
+                    srcCellI = cellS;
+                    tgtCellI = seedCells[cellS];
+
+                    return;
+                }
+            }
+        }
+
+        // perform new search to find match
+        if (debug)
+        {
+            Pout<< "Advancing front stalled: searching for new "
+                << "target cell" << endl;
+        }
+
+        bool restart =
+            findInitialSeeds
+            (
+                src,
+                tgt,
+                srcCellIDs,
+                mapFlag,
+                startSeedI,
+                srcCellI,
+                tgtCellI
+            );
+
+        if (restart)
+        {
+            // successfully found new starting seed-pair
+            return;
+        }
+    }
+
+    // if we have got to here, there are no more src/tgt cell intersections
+    srcCellI = -1;
+    tgtCellI = -1;
+}
+
+
+bool Foam::meshToMeshNew::intersect
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const label srcCellI,
+    const label tgtCellI
+) const
+{
+    scalar threshold = tolerance_*src.cellVolumes()[srcCellI];
+
+    tetOverlapVolume overlapEngine;
+
+    treeBoundBox bbTgtCell
+    (
+        pointField
+        (
+            tgt.points(),
+            tgt.cellPoints()[tgtCellI]
+        )
+    );
+
+    return overlapEngine.cellCellOverlapMinDecomp
+    (
+        src,
+        srcCellI,
+        tgt,
+        tgtCellI,
+        bbTgtCell,
+        threshold
+    );
+}
+
+
+Foam::scalar Foam::meshToMeshNew::interVol
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const label srcCellI,
+    const label tgtCellI
+) const
+{
+    tetOverlapVolume overlapEngine;
+
+    treeBoundBox bbTgtCell
+    (
+        pointField
+        (
+            tgt.points(),
+            tgt.cellPoints()[tgtCellI]
+        )
+    );
+
+    scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
+    (
+        src,
+        srcCellI,
+        tgt,
+        tgtCellI,
+        bbTgtCell
+    );
+
+    return vol;
+}
+
+
+void Foam::meshToMeshNew::calcIndirect
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const label srcSeedI,
+    const label tgtSeedI,
+    const labelList& srcCellIDs,
+    boolList& mapFlag,
+    label& startSeedI
+)
+{
+    label srcCellI = srcSeedI;
+    label tgtCellI = tgtSeedI;
+
+    List<DynamicList<label> > srcToTgtAddr(src.nCells());
+    List<DynamicList<scalar> > srcToTgtWght(src.nCells());
+
+    List<DynamicList<label> > tgtToSrcAddr(tgt.nCells());
+    List<DynamicList<scalar> > tgtToSrcWght(tgt.nCells());
+
+    // list of tgt cell neighbour cells
+    DynamicList<label> nbrTgtCells(10);
+
+    // list of tgt cells currently visited for srcCellI to avoid multiple hits
+    DynamicList<label> visitedTgtCells(10);
+
+    // list to keep track of tgt cells used to seed src cells
+    labelList seedCells(src.nCells(), -1);
+    seedCells[srcCellI] = tgtCellI;
+
+    const scalarField& srcVol = src.cellVolumes();
+
+    do
+    {
+        nbrTgtCells.clear();
+        visitedTgtCells.clear();
+
+        // append initial target cell and neighbours
+        nbrTgtCells.append(tgtCellI);
+        appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+
+        do
+        {
+            tgtCellI = nbrTgtCells.remove();
+            visitedTgtCells.append(tgtCellI);
+
+            scalar vol = interVol(src, tgt, srcCellI, tgtCellI);
+
+            // accumulate addressing and weights for valid intersection
+            if (vol/srcVol[srcCellI] > tolerance_)
+            {
+                // store src/tgt cell pair
+                srcToTgtAddr[srcCellI].append(tgtCellI);
+                srcToTgtWght[srcCellI].append(vol);
+
+                tgtToSrcAddr[tgtCellI].append(srcCellI);
+                tgtToSrcWght[tgtCellI].append(vol);
+
+                appendNbrTgtCells(tgtCellI, tgt, visitedTgtCells, nbrTgtCells);
+
+                // accumulate intersection volume
+                V_ += vol;
+            }
+        }
+        while (!nbrTgtCells.empty());
+
+        mapFlag[srcCellI] = false;
+
+        // find new source seed cell
+        setNextCells
+        (
+            startSeedI,
+            srcCellI,
+            tgtCellI,
+            src,
+            tgt,
+            srcCellIDs,
+            mapFlag,
+            visitedTgtCells,
+            seedCells
+        );
+    }
+    while (srcCellI != -1);
+
+    // transfer addressing into persistent storage
+    forAll(srcToTgtCellAddr_, i)
+    {
+        srcToTgtCellAddr_[i].transfer(srcToTgtAddr[i]);
+        srcToTgtCellWght_[i].transfer(srcToTgtWght[i]);
+    }
+
+    forAll(tgtToSrcCellAddr_, i)
+    {
+        tgtToSrcCellAddr_[i].transfer(tgtToSrcAddr[i]);
+        tgtToSrcCellWght_[i].transfer(tgtToSrcWght[i]);
+    }
+}
+
+
+void Foam::meshToMeshNew::calcAddressing
+(
+    const polyMesh& src,
+    const polyMesh& tgt
+)
+{
+    srcToTgtCellAddr_.setSize(src.nCells());
+    srcToTgtCellWght_.setSize(src.nCells());
+
+    tgtToSrcCellAddr_.setSize(tgt.nCells());
+    tgtToSrcCellWght_.setSize(tgt.nCells());
+
+    if (!src.nCells() || !tgt.nCells())
+    {
+        if (debug)
+        {
+            Pout<< "mesh interpolation: cells not on processor: Source cells = "
+                << src.nCells() << ", target cells = " << tgt.nCells()
+                << endl;
+        }
+    }
+
+    if (!src.nCells())
+    {
+        return;
+    }
+    else if (!tgt.nCells())
+    {
+        if (debug)
+        {
+            Pout<< "mesh interpolation: hhave " << src.nCells() << " source "
+                << " cells but no target cells" << endl;
+        }
+
+        return;
+    }
+
+    // (potentially) participating source mesh cells
+    const labelList srcCellIDs = maskCells(src, tgt);
+
+    // list to keep track of whether src cell can be mapped
+    boolList mapFlag(src.nCells(), false);
+    UIndirectList<bool>(mapFlag, srcCellIDs) = true;
+
+    // find initial point in tgt mesh
+    label srcSeedI = -1;
+    label tgtSeedI = -1;
+    label startSeedI = 0;
+
+    bool startWalk =
+        findInitialSeeds
+        (
+            src,
+            tgt,
+            srcCellIDs,
+            mapFlag,
+            startSeedI,
+            srcSeedI,
+            tgtSeedI
+        );
+
+    if (!startWalk)
+    {
+        // if meshes are collocated, after inflating the source mesh bounding
+        // box tgt mesh cells may be transferred, but may still not overlap
+        // with the source mesh
+        return;
+    }
+
+
+    switch (method_)
+    {
+        case imMap:
+        {
+            calcDirect(src, tgt, srcSeedI, tgtSeedI);
+            break;
+        }
+        case imCellVolumeWeight:
+        {
+            calcIndirect
+            (
+                src,
+                tgt,
+                srcSeedI,
+                tgtSeedI,
+                srcCellIDs,
+                mapFlag,
+                startSeedI
+            );
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::meshToMeshNew::calcAddressing"
+                "("
+                    "const polyMesh&, "
+                    "const polyMesh&"
+                ")"
+            )
+                << "Unknown interpolation method"
+                << abort(FatalError);
+        }
+    }
+
+
+    if (debug)
+    {
+        writeConnectivity(src, tgt, srcToTgtCellAddr_);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshToMeshNew::meshToMeshNew
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const interpolationMethod& method
+)
+:
+    srcRegionName_(src.name()),
+    tgtRegionName_(tgt.name()),
+    srcToTgtCellAddr_(),
+    tgtToSrcCellAddr_(),
+    srcToTgtCellWght_(),
+    tgtToSrcCellWght_(),
+    method_(method),
+    V_(0.0),
+    singleMeshProc_(-1),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    Info<< "Creating mesh-to-mesh addressing for " << src.name()
+        << " and " << tgt.name() << " regions" << endl;
+
+    singleMeshProc_ = calcDistribution(src, tgt);
+
+    if (singleMeshProc_ == -1)
+    {
+        // create global indexing for src and tgt meshes
+        globalIndex globalSrcCells(src.nCells());
+        globalIndex globalTgtCells(tgt.nCells());
+
+        // Create processor map of overlapping cells. This map gets
+        // (possibly remote) cells from the tgt mesh such that they (together)
+        // cover all of the src mesh
+        autoPtr<mapDistribute> mapPtr = calcProcMap(src, tgt);
+        const mapDistribute& map = mapPtr();
+
+        pointField newTgtPoints;
+        faceList newTgtFaces;
+        labelList newTgtFaceOwners;
+        labelList newTgtFaceNeighbours;
+        labelList newTgtCellIDs;
+
+        distributeAndMergeCells
+        (
+            map,
+            tgt,
+            globalTgtCells,
+            newTgtPoints,
+            newTgtFaces,
+            newTgtFaceOwners,
+            newTgtFaceNeighbours,
+            newTgtCellIDs
+        );
+
+
+        // create a new target mesh
+        polyMesh newTgt
+        (
+            IOobject
+            (
+                "newTgt::" + Foam::name(Pstream::myProcNo()),
+                tgt.time().timeName(),
+                tgt.time(),
+                IOobject::NO_READ
+            ),
+            xferMove(newTgtPoints),
+            xferMove(newTgtFaces),
+            xferMove(newTgtFaceOwners),
+            xferMove(newTgtFaceNeighbours),
+            false                                   // no parallel comms
+        );
+
+        // create some dummy patch info
+        List<polyPatch*> patches(1);
+        patches[0] = new polyPatch
+        (
+            "defaultFaces",
+            newTgt.nFaces() - newTgt.nInternalFaces(),
+            newTgt.nInternalFaces(),
+            0,
+            newTgt.boundaryMesh(),
+            word::null
+        );
+
+        newTgt.addPatches(patches);
+
+        // force calculation of tet-base points used for point-in-cell
+        (void)newTgt.tetBasePtIs();
+
+        // force construction of cell tree
+//        (void)newTgt.cellTree();
+
+        if (debug)
+        {
+            Pout<< "Created newTgt mesh:" << nl
+                << " old cells = " << tgt.nCells()
+                << ", new cells = " << newTgt.nCells() << nl
+                << " old faces = " << tgt.nFaces()
+                << ", new faces = " << newTgt.nFaces() << endl;
+
+            if (debug > 1)
+            {
+                Pout<< "Writing newTgt mesh: " << newTgt.name() << endl;
+                newTgt.write();
+            }
+        }
+
+        calcAddressing(src, newTgt);
+
+        // per source cell the target cell address in newTgt mesh
+        forAll(srcToTgtCellAddr_, i)
+        {
+            labelList& addressing = srcToTgtCellAddr_[i];
+            forAll(addressing, addrI)
+            {
+                addressing[addrI] = newTgtCellIDs[addressing[addrI]];
+            }
+        }
+
+        // convert target addresses in newTgtMesh into global cell numbering
+        forAll(tgtToSrcCellAddr_, i)
+        {
+            labelList& addressing = tgtToSrcCellAddr_[i];
+            forAll(addressing, addrI)
+            {
+                addressing[addrI] = globalSrcCells.toGlobal(addressing[addrI]);
+            }
+        }
+
+        // set up as a reverse distribute
+        mapDistribute::distribute
+        (
+            Pstream::nonBlocking,
+            List<labelPair>(),
+            tgt.nCells(),
+            map.constructMap(),
+            map.subMap(),
+            tgtToSrcCellAddr_,
+            ListPlusEqOp<label>(),
+            labelList()
+        );
+
+        // set up as a reverse distribute
+        mapDistribute::distribute
+        (
+            Pstream::nonBlocking,
+            List<labelPair>(),
+            tgt.nCells(),
+            map.constructMap(),
+            map.subMap(),
+            tgtToSrcCellWght_,
+            ListPlusEqOp<scalar>(),
+            scalarList()
+        );
+
+        // weights normalisation
+        normaliseWeights
+        (
+            "source",
+            src.cellVolumes(),
+            srcToTgtCellAddr_,
+            srcToTgtCellWght_
+        );
+
+        normaliseWeights
+        (
+            "target",
+            tgt.cellVolumes(),
+            tgtToSrcCellAddr_,
+            tgtToSrcCellWght_
+        );
+
+        // cache maps and reset addresses
+        List<Map<label> > cMap;
+        srcMapPtr_.reset
+        (
+            new mapDistribute(globalSrcCells, tgtToSrcCellAddr_, cMap)
+        );
+        tgtMapPtr_.reset
+        (
+            new mapDistribute(globalTgtCells, srcToTgtCellAddr_, cMap)
+        );
+
+        // collect volume intersection contributions
+        reduce(V_, sumOp<scalar>());
+    }
+    else
+    {
+        calcAddressing(src, tgt);
+
+        normaliseWeights
+        (
+            "source",
+            src.cellVolumes(),
+            srcToTgtCellAddr_,
+            srcToTgtCellWght_
+        );
+
+        normaliseWeights
+        (
+            "target",
+            tgt.cellVolumes(),
+            tgtToSrcCellAddr_,
+            tgtToSrcCellWght_
+        );
+    }
+
+    Info<< "    Overlap volume: " << V_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::meshToMeshNew::~meshToMeshNew()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
new file mode 100644
index 0000000000..85f0c1332e
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
@@ -0,0 +1,511 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::meshToMeshNew
+
+Description
+    Class to calculate the cell-addressing between two overlapping meshes
+
+SourceFiles
+    meshToMeshNew.C
+    meshToMeshNewTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshToMeshNew_H
+#define meshToMeshNew_H
+
+#include "polyMesh.H"
+#include "boundBox.H"
+#include "mapDistribute.H"
+#include "volFieldsFwd.H"
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class meshToMeshNew Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshToMeshNew
+{
+public:
+
+    // Public data types
+
+        //- Enumeration specifying required accuracy
+        enum interpolationMethod
+        {
+            imMap,
+            imCellVolumeWeight
+        };
+
+        static const NamedEnum<interpolationMethod, 2>
+            interpolationMethodNames_;
+
+private:
+
+    // Private data
+
+        //- Name of source mesh region
+        const word srcRegionName_;
+
+        //- Name of target mesh region
+        const word tgtRegionName_;
+
+        //- Source to target cell addressing
+        labelListList srcToTgtCellAddr_;
+
+        //- Target to source cell addressing
+        labelListList tgtToSrcCellAddr_;
+
+        //- Source to target cell interplation weights
+        scalarListList srcToTgtCellWght_;
+
+        //- Target to source cell interpolation weights
+        scalarListList tgtToSrcCellWght_;
+
+        //- Interpolation method
+        interpolationMethod method_;
+
+        //- Cell total volume in overlap region [m3]
+        scalar V_;
+
+        //- Index of processor that holds all of both sides. -1 in all other
+        //  cases
+        label singleMeshProc_;
+
+        //- Source map pointer - parallel running only
+        autoPtr<mapDistribute> srcMapPtr_;
+
+        //- Target map pointer - parallel running only
+        autoPtr<mapDistribute> tgtMapPtr_;
+
+        //- Tolerance used in volume overlap calculations
+        static scalar tolerance_;
+
+
+    // Private Member Functions
+
+        //- Helper function to add a constant offset to a list
+        template<class Type>
+        void add(UList<Type>& fld, const label offset) const;
+
+        //- Write the connectivity (debugging)
+        void writeConnectivity
+        (
+            const polyMesh& src,
+            const polyMesh& tgt,
+            const labelListList& srcToTargetAddr
+        ) const;
+
+        //- Return src cell IDs for the overlap region
+        labelList maskCells(const polyMesh& src, const polyMesh& tgt) const;
+
+        //- Find indices of overlapping cells in src and tgt meshes - returns
+        //  true if found a matching pair
+        bool findInitialSeeds
+        (
+            const polyMesh& src,
+            const polyMesh& tgt,
+            const labelList& srcCellIDs,
+            const boolList& mapFlag,
+            const label startSeedI,
+            label& srcSeedI,
+            label& tgtSeedI
+        ) const;
+
+
+        // Direct (one-to-one) mapping
+
+            //- Append to list of src mesh seed indices
+            void appendToDirectSeeds
+            (
+                const polyMesh& src,
+                const polyMesh& tgt,
+                boolList& mapFlag,
+                labelList& srcTgtSeed,
+                DynamicList<label>& srcSeeds,
+                label& srcSeedI,
+                label& tgtSeedI
+            ) const;
+
+            //- Main driver routine for direct mapping
+            void calcDirect
+            (
+                const polyMesh& src,
+                const polyMesh& tgt,
+                const label srcSeedI,
+                const label tgtSeedI
+            );
+
+
+        // Indirect (non-conformal) mapping
+
+            //- Normalise the interpolation weights
+            void normaliseWeights
+            (
+                const word& descriptor,
+                const scalarField& cellVolumes,
+                const labelListList& addr,
+                scalarListList& wght
+            ) const;
+
+            //- Append target cell neihgbour cells to cellIDs list
+            void appendNbrTgtCells
+            (
+                const label tgtCellI,
+                const polyMesh& tgt,
+                const DynamicList<label>& visitedTgtCells,
+                DynamicList<label>& nbrTgtCellIDs
+            ) const;
+
+            //- Set the next cells in the advancing front algorithm
+            void setNextCells
+            (
+                label& startSeedI,
+                label& srcCellI,
+                label& tgtCellI,
+                const polyMesh& src,
+                const polyMesh& tgt,
+                const labelList& srcCellIDs,
+                const boolList& mapFlag,
+                const DynamicList<label>& visitedCells,
+                labelList& seedCells
+            ) const;
+
+            //- Return the true if cells intersect
+            bool intersect
+            (
+                const polyMesh& src,
+                const polyMesh& tgt,
+                const label srcCellI,
+                const label tgtCellI
+            ) const;
+
+            //- Return the intersection volume between two cells
+            scalar interVol
+            (
+                const polyMesh& src,
+                const polyMesh& tgt,
+                const label srcCellI,
+                const label tgtCellI
+            ) const;
+
+            //- Main driver routine for indirect mapping
+            void calcIndirect
+            (
+                const polyMesh& src,
+                const polyMesh& tgt,
+                const label srcSeedI,
+                const label tgtSeedI,
+                const labelList& srcCellIDs,
+                boolList& mapFlag,
+                label& startSeedI
+            );
+
+
+        //- Calculate the addressing between overalping regions of src and tgt
+        //  meshes - main driver function
+        void calcAddressing(const polyMesh& src, const polyMesh& tgt);
+
+
+        // Parallel operations
+
+            //- Determine whether the meshes are split across multiple pocessors
+            label calcDistribution
+            (
+                const polyMesh& src,
+                const polyMesh& tgt
+            ) const;
+
+            //- Determine which processor bounding-boxes overlap
+            label calcOverlappingProcs
+            (
+                const List<boundBox>& procBb,
+                const boundBox& bb,
+                boolList& overlaps
+            ) const;
+
+            //- Calculate the mapping between processors
+            autoPtr<mapDistribute> calcProcMap
+            (
+                const polyMesh& src,
+                const polyMesh& tgt
+            ) const;
+
+            //- Distribute mesh info from 'my' processor to others
+            void distributeCells
+            (
+                const mapDistribute& map,
+                const polyMesh& tgtMesh,
+                const globalIndex& globalI,
+                List<pointField>& points,
+                List<label>& nInternalFaces,
+                List<faceList>& faces,
+                List<labelList>& faceOwner,
+                List<labelList>& faceNeighbour,
+                List<labelList>& cellIDs,
+                List<labelList>& nbrProcIDs,
+                List<labelList>& procLocalFaceIDs
+            ) const;
+
+            //- Collect pieces of tgt mesh from other procssors and restructure
+            void distributeAndMergeCells
+            (
+                const mapDistribute& map,
+                const polyMesh& tgt,
+                const globalIndex& globalI,
+                pointField& tgtPoints,
+                faceList& tgtFaces,
+                labelList& tgtFaceOwners,
+                labelList& tgtFaceNeighbours,
+                labelList& tgtCellIDs
+            ) const;
+
+
+        //- Disallow default bitwise copy construct
+        meshToMeshNew(const meshToMeshNew&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const meshToMeshNew&);
+
+
+public:
+
+    //- Run-time type information
+    TypeName("meshToMeshNew");
+
+
+    //- Construct from source and target meshes
+    meshToMeshNew
+    (
+        const polyMesh& src,
+        const polyMesh& tgt,
+        const interpolationMethod& method
+    );
+
+
+    //- Destructor
+    virtual ~meshToMeshNew();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the source to target cell addressing
+            inline const labelListList& srcToTgtCellAddr() const;
+
+            //- Return const access to the target to source cell addressing
+            inline const labelListList& tgtToSrcCellAddr() const;
+
+            //- Return const access to the source to target cell weights
+            inline const scalarListList& srcToTgtCellWght() const;
+
+            //- Return const access to the target to source cell weights
+            inline const scalarListList& tgtToSrcCellWght() const;
+
+            //- Return const access to the overlap volume
+            inline scalar V() const;
+
+
+        // Evaluation
+
+            // Source-to-target field mapping
+
+                //- Map field from src to tgt mesh with defined operation
+                //  Values passed in via 'result' are used to initialise the
+                //  return value
+                template<class Type, class CombineOp>
+                void mapSrcToTgt
+                (
+                    const UList<Type>& srcFld,
+                    const CombineOp& cop,
+                    List<Type>& result
+                ) const;
+
+                //- Return the src field mapped to the tgt mesh with a defined
+                //  operation.  Initial values of the result are set to zero
+                template<class Type, class CombineOp>
+                tmp<Field<Type> > mapSrcToTgt
+                (
+                    const Field<Type>& srcFld,
+                    const CombineOp& cop
+                ) const;
+
+                //- Convenience function to map a tmp field to the tgt mesh
+                //  with a defined operation
+                template<class Type, class CombineOp>
+                tmp<Field<Type> > mapSrcToTgt
+                (
+                    const tmp<Field<Type> >& tsrcFld,
+                    const CombineOp& cop
+                ) const;
+
+                //- Convenience function to map a field to the tgt mesh with a
+                //  default operation (plusEqOp)
+                template<class Type>
+                tmp<Field<Type> > mapSrcToTgt
+                (
+                    const Field<Type>& srcFld
+                ) const;
+
+                //- Convenience function to map a tmp field to the tgt mesh
+                //  with a default operation (plusEqOp)
+                template<class Type>
+                tmp<Field<Type> > mapSrcToTgt
+                (
+                    const tmp<Field<Type> >& tsrcFld
+                ) const;
+
+
+            // Target-to-source field mapping
+
+                //- Map field from tgt to src mesh with defined operation
+                //  Values passed in via 'result' are used to initialise the
+                //  return value
+                template<class Type, class CombineOp>
+                void mapTgtToSrc
+                (
+                    const UList<Type>& tgtFld,
+                    const CombineOp& cop,
+                    List<Type>& result
+                ) const;
+
+                //- Return the tgt field mapped to the src mesh with a defined
+                //  operation.  Initial values of the result are set to zero
+                template<class Type, class CombineOp>
+                tmp<Field<Type> > mapTgtToSrc
+                (
+                    const Field<Type>& tgtFld,
+                    const CombineOp& cop
+                ) const;
+
+                //- Convenience function to map a tmp field to the src mesh
+                //  with a defined operation
+                template<class Type, class CombineOp>
+                tmp<Field<Type> > mapTgtToSrc
+                (
+                    const tmp<Field<Type> >& ttgtFld,
+                    const CombineOp& cop
+                ) const;
+
+                //- Convenience function to map a field to the src mesh with a
+                //  default operation (plusEqOp)
+                template<class Type>
+                tmp<Field<Type> > mapTgtToSrc
+                (
+                    const Field<Type>& tgtFld
+                ) const;
+
+                //- Convenience function to map a tmp field to the src mesh
+                //  with a default operation (plusEqOp)
+                template<class Type>
+                tmp<Field<Type> > mapTgtToSrc
+                (
+                    const tmp<Field<Type> >& ttgtFld
+                ) const;
+
+
+            // Volume field mapping
+
+                //- Interpolare a field with a defined operation.  Values
+                //  passed in via 'result' are used to initialise the return
+                //  value.  Optionally interpolate patch values
+                template<class Type, class CombineOp>
+                void interpolate
+                (
+                    const GeometricField<Type, fvPatchField, volMesh>& field,
+                    const CombineOp& cop,
+                    GeometricField<Type, fvPatchField, volMesh>& result,
+                    const bool interpPatches = true
+                ) const;
+
+                //- Interpolare a field with a defined operation.  The initial
+                //  values of the result are set to zero.  Optionally
+                //  interpolate patch values
+                template<class Type, class CombineOp>
+                tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+                (
+                    const GeometricField<Type, fvPatchField, volMesh>& field,
+                    const CombineOp& cop,
+                    const bool interpPatches = true
+                ) const;
+
+                //- Interpolare a tmp field with a defined operation.  The
+                //  initial values of the result are set to zero.  Optionally
+                //  interpolate patch values
+                template<class Type, class CombineOp>
+                tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+                (
+                    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
+                        tfield,
+                    const CombineOp& cop,
+                    const bool interpPatches = true
+                ) const;
+
+                //- Convenience function to map a field with a default
+                //  operation (plusEqOp).  Optionally interpolate patch values
+                template<class Type>
+                tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+                (
+                    const GeometricField<Type, fvPatchField, volMesh>& field,
+                    const bool interpPatches = true
+                ) const;
+
+                //- Convenience function to map a tmp field with a default
+                //  operation (plusEqOp).  Optionally interpolate patch values
+                template<class Type>
+                tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+                (
+                    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
+                        tfield,
+                    const bool interpPatches = true
+                ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "meshToMeshNewI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "meshToMeshNewTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H
new file mode 100644
index 0000000000..41068e8c45
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewI.H
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshToMeshNew.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline const Foam::labelListList&
+Foam::meshToMeshNew::srcToTgtCellAddr() const
+{
+    return srcToTgtCellAddr_;
+}
+
+
+inline const Foam::labelListList&
+Foam::meshToMeshNew::tgtToSrcCellAddr() const
+{
+    return tgtToSrcCellAddr_;
+}
+
+
+inline const Foam::scalarListList&
+Foam::meshToMeshNew::srcToTgtCellWght() const
+{
+    return srcToTgtCellWght_;
+}
+
+
+inline const Foam::scalarListList&
+Foam::meshToMeshNew::tgtToSrcCellWght() const
+{
+    return tgtToSrcCellWght_;
+}
+
+
+inline Foam::scalar Foam::meshToMeshNew::V() const
+{
+    return V_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C
new file mode 100644
index 0000000000..3cb9d0597d
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewParallelOps.C
@@ -0,0 +1,884 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshToMeshNew.H"
+#include "OFstream.H"
+#include "Time.H"
+#include "globalIndex.H"
+#include "mergePoints.H"
+#include "processorPolyPatch.H"
+#include "SubField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::meshToMeshNew::calcDistribution
+(
+    const polyMesh& src,
+    const polyMesh& tgt
+) const
+{
+    label procI = 0;
+
+    if (Pstream::parRun())
+    {
+        List<label> cellsPresentOnProc(Pstream::nProcs(), 0);
+        if ((src.nCells() > 0) || (tgt.nCells() > 0))
+        {
+            cellsPresentOnProc[Pstream::myProcNo()] = 1;
+        }
+        else
+        {
+            cellsPresentOnProc[Pstream::myProcNo()] = 0;
+        }
+
+        Pstream::gatherList(cellsPresentOnProc);
+        Pstream::scatterList(cellsPresentOnProc);
+
+        label nHaveCells = sum(cellsPresentOnProc);
+
+        if (nHaveCells > 1)
+        {
+            procI = -1;
+            if (debug)
+            {
+                Info<< "meshToMeshNew::calcDistribution: "
+                    << "Meshes split across multiple processors" << endl;
+            }
+        }
+        else if (nHaveCells == 1)
+        {
+            procI = findIndex(cellsPresentOnProc, 1);
+            if (debug)
+            {
+                Info<< "meshToMeshNew::calcDistribution: "
+                    << "Meshes local to processor" << procI << endl;
+            }
+        }
+    }
+
+    return procI;
+}
+
+
+Foam::label Foam::meshToMeshNew::calcOverlappingProcs
+(
+    const List<boundBox>& procBb,
+    const boundBox& bb,
+    boolList& overlaps
+) const
+{
+    overlaps = false;
+
+    label nOverlaps = 0;
+
+    forAll(procBb, procI)
+    {
+        const boundBox& bbp = procBb[procI];
+
+        if (bbp.overlaps(bb))
+        {
+            overlaps[procI] = true;
+            nOverlaps++;
+        }
+    }
+
+    return nOverlaps;
+}
+
+
+Foam::autoPtr<Foam::mapDistribute> Foam::meshToMeshNew::calcProcMap
+(
+    const polyMesh& src,
+    const polyMesh& tgt
+) const
+{
+    // get decomposition of cells on src mesh
+    List<boundBox> procBb(Pstream::nProcs());
+
+    if (src.nCells() > 0)
+    {
+        // bounding box for my mesh - do not parallel reduce
+        procBb[Pstream::myProcNo()] = boundBox(src.points(), false);
+
+        // slightly increase size of bounding boxes to allow for cases where
+        // bounding boxes are perfectly alligned
+        procBb[Pstream::myProcNo()].inflate(0.01);
+    }
+    else
+    {
+        procBb[Pstream::myProcNo()] = boundBox();
+    }
+
+
+    Pstream::gatherList(procBb);
+    Pstream::scatterList(procBb);
+
+
+    if (debug)
+    {
+        Info<< "Determining extent of src mesh per processor:" << nl
+            << "\tproc\tbb" << endl;
+        forAll(procBb, procI)
+        {
+            Info<< '\t' << procI << '\t' << procBb[procI] << endl;
+        }
+    }
+
+
+    // determine which cells of tgt mesh overlaps src mesh per proc
+    const cellList& cells = tgt.cells();
+    const faceList& faces = tgt.faces();
+    const pointField& points = tgt.points();
+
+    labelListList sendMap;
+
+    {
+        // per processor indices into all segments to send
+        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
+        label iniSize = floor(tgt.nCells()/Pstream::nProcs());
+
+        forAll(dynSendMap, procI)
+        {
+            dynSendMap[procI].setCapacity(iniSize);
+        }
+
+        // work array - whether src processor bb overlaps the tgt cell bounds
+        boolList procBbOverlaps(Pstream::nProcs());
+        forAll(cells, cellI)
+        {
+            const cell& c = cells[cellI];
+
+            // determine bounding box of tgt cell
+            boundBox cellBb(point::max, point::min);
+            forAll(c, faceI)
+            {
+                const face& f = faces[c[faceI]];
+                forAll(f, fp)
+                {
+                    cellBb.min() = min(cellBb.min(), points[f[fp]]);
+                    cellBb.max() = max(cellBb.max(), points[f[fp]]);
+                }
+            }
+
+            // find the overlapping tgt cells on each src processor
+            (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
+
+            forAll(procBbOverlaps, procI)
+            {
+                if (procBbOverlaps[procI])
+                {
+                    dynSendMap[procI].append(cellI);
+                }
+            }
+        }
+
+        // convert dynamicList to labelList
+        sendMap.setSize(Pstream::nProcs());
+        forAll(sendMap, procI)
+        {
+            sendMap[procI].transfer(dynSendMap[procI]);
+        }
+    }
+
+    // debug printing
+    if (debug)
+    {
+        Pout<< "Of my " << cells.size() << " target cells I need to send to:"
+            << nl << "\tproc\tcells" << endl;
+        forAll(sendMap, procI)
+        {
+            Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
+        }
+    }
+
+
+    // send over how many tgt cells I need to receive from each processor
+    labelListList sendSizes(Pstream::nProcs());
+    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(sendMap, procI)
+    {
+        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+    }
+    Pstream::gatherList(sendSizes);
+    Pstream::scatterList(sendSizes);
+
+
+    // determine order of receiving
+    labelListList constructMap(Pstream::nProcs());
+
+    label segmentI = 0;
+    forAll(constructMap, procI)
+    {
+        // what I need to receive is what other processor is sending to me
+        label nRecv = sendSizes[procI][Pstream::myProcNo()];
+        constructMap[procI].setSize(nRecv);
+
+        for (label i = 0; i < nRecv; i++)
+        {
+            constructMap[procI][i] = segmentI++;
+        }
+    }
+
+    autoPtr<mapDistribute> mapPtr
+    (
+        new mapDistribute
+        (
+            segmentI,       // size after construction
+            sendMap.xfer(),
+            constructMap.xfer()
+        )
+    );
+
+    return mapPtr;
+}
+
+
+void Foam::meshToMeshNew::distributeCells
+(
+    const mapDistribute& map,
+    const polyMesh& tgtMesh,
+    const globalIndex& globalI,
+    List<pointField>& points,
+    List<label>& nInternalFaces,
+    List<faceList>& faces,
+    List<labelList>& faceOwner,
+    List<labelList>& faceNeighbour,
+    List<labelList>& cellIDs,
+    List<labelList>& nbrProcIDs,
+    List<labelList>& procLocalFaceIDs
+) const
+{
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    points.setSize(Pstream::nProcs());
+    nInternalFaces.setSize(Pstream::nProcs(), 0);
+    faces.setSize(Pstream::nProcs());
+    faceOwner.setSize(Pstream::nProcs());
+    faceNeighbour.setSize(Pstream::nProcs());
+    cellIDs.setSize(Pstream::nProcs());
+
+    nbrProcIDs.setSize(Pstream::nProcs());;
+    procLocalFaceIDs.setSize(Pstream::nProcs());;
+
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& sendElems = map.subMap()[domain];
+
+        if (sendElems.size())
+        {
+            // reverse cell map
+            labelList reverseCellMap(tgtMesh.nCells(), -1);
+            forAll(sendElems, subCellI)
+            {
+                reverseCellMap[sendElems[subCellI]] = subCellI;
+            }
+
+            DynamicList<face> subFaces(tgtMesh.nFaces());
+            DynamicList<label> subFaceOwner(tgtMesh.nFaces());
+            DynamicList<label> subFaceNeighbour(tgtMesh.nFaces());
+
+            DynamicList<label> subNbrProcIDs(tgtMesh.nFaces());
+            DynamicList<label> subProcLocalFaceIDs(tgtMesh.nFaces());
+
+            label nInternal = 0;
+
+            // internal faces
+            forAll(tgtMesh.faceNeighbour(), faceI)
+            {
+                label own = tgtMesh.faceOwner()[faceI];
+                label nbr = tgtMesh.faceNeighbour()[faceI];
+                label subOwn = reverseCellMap[own];
+                label subNbr = reverseCellMap[nbr];
+
+                if (subOwn != -1 && subNbr != -1)
+                {
+                    nInternal++;
+
+                    if (subOwn < subNbr)
+                    {
+                        subFaces.append(tgtMesh.faces()[faceI]);
+                        subFaceOwner.append(subOwn);
+                        subFaceNeighbour.append(subNbr);
+                        subNbrProcIDs.append(-1);
+                        subProcLocalFaceIDs.append(-1);
+                    }
+                    else
+                    {
+                        subFaces.append(tgtMesh.faces()[faceI].reverseFace());
+                        subFaceOwner.append(subNbr);
+                        subFaceNeighbour.append(subOwn);
+                        subNbrProcIDs.append(-1);
+                        subProcLocalFaceIDs.append(-1);
+                    }
+                }
+            }
+
+            // boundary faces for new region
+            forAll(tgtMesh.faceNeighbour(), faceI)
+            {
+                label own = tgtMesh.faceOwner()[faceI];
+                label nbr = tgtMesh.faceNeighbour()[faceI];
+                label subOwn = reverseCellMap[own];
+                label subNbr = reverseCellMap[nbr];
+
+                if (subOwn != -1 && subNbr == -1)
+                {
+                    subFaces.append(tgtMesh.faces()[faceI]);
+                    subFaceOwner.append(subOwn);
+                    subFaceNeighbour.append(subNbr);
+                    subNbrProcIDs.append(-1);
+                    subProcLocalFaceIDs.append(-1);
+                }
+                else if (subOwn == -1 && subNbr != -1)
+                {
+                    subFaces.append(tgtMesh.faces()[faceI].reverseFace());
+                    subFaceOwner.append(subNbr);
+                    subFaceNeighbour.append(subOwn);
+                    subNbrProcIDs.append(-1);
+                    subProcLocalFaceIDs.append(-1);
+                }
+            }
+
+            // boundary faces of existing region
+            forAll(tgtMesh.boundaryMesh(), patchI)
+            {
+                const polyPatch& pp = tgtMesh.boundaryMesh()[patchI];
+
+                label nbrProcI = -1;
+
+                // store info for faces on processor patches
+                if (isA<processorPolyPatch>(pp))
+                {
+                    const processorPolyPatch& ppp =
+                        dynamic_cast<const processorPolyPatch&>(pp);
+
+                    nbrProcI = ppp.neighbProcNo();
+                }
+
+                forAll(pp, i)
+                {
+                    label faceI = pp.start() + i;
+                    label own = tgtMesh.faceOwner()[faceI];
+
+                    if (reverseCellMap[own] != -1)
+                    {
+                        subFaces.append(tgtMesh.faces()[faceI]);
+                        subFaceOwner.append(reverseCellMap[own]);
+                        subFaceNeighbour.append(-1);
+                        subNbrProcIDs.append(nbrProcI);
+                        subProcLocalFaceIDs.append(i);
+                    }
+                }
+            }
+
+            // reverse point map
+            labelList reversePointMap(tgtMesh.nPoints(), -1);
+            DynamicList<point> subPoints(tgtMesh.nPoints());
+            forAll(subFaces, subFaceI)
+            {
+                face& f = subFaces[subFaceI];
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+                    if (reversePointMap[pointI] == -1)
+                    {
+                        reversePointMap[pointI] = subPoints.size();
+                        subPoints.append(tgtMesh.points()[pointI]);
+                    }
+
+                    f[fp] = reversePointMap[pointI];
+                }
+            }
+
+            // tgt cells into global numbering
+            labelList globalElems(sendElems.size());
+            forAll(sendElems, i)
+            {
+                if (debug)
+                {
+                    Pout<< "tgtProc:" << Pstream::myProcNo()
+                        << " sending tgt cell " << sendElems[i]
+                        << "[" << globalI.toGlobal(sendElems[i]) << "]"
+                        << " to srcProc " << domain << endl;
+                }
+
+                globalElems[i] = globalI.toGlobal(sendElems[i]);
+            }
+
+            // pass data
+            if (domain == Pstream::myProcNo())
+            {
+                // allocate my own data
+                points[Pstream::myProcNo()] = subPoints;
+                nInternalFaces[Pstream::myProcNo()] = nInternal;
+                faces[Pstream::myProcNo()] = subFaces;
+                faceOwner[Pstream::myProcNo()] = subFaceOwner;
+                faceNeighbour[Pstream::myProcNo()] = subFaceNeighbour;
+                cellIDs[Pstream::myProcNo()] = globalElems;
+                nbrProcIDs[Pstream::myProcNo()] = subNbrProcIDs;
+                procLocalFaceIDs[Pstream::myProcNo()] = subProcLocalFaceIDs;
+            }
+            else
+            {
+                // send data to other processor domains
+                UOPstream toDomain(domain, pBufs);
+
+                toDomain
+                    << subPoints
+                    << nInternal
+                    << subFaces
+                    << subFaceOwner
+                    << subFaceNeighbour
+                    << globalElems
+                    << subNbrProcIDs
+                    << subProcLocalFaceIDs;
+            }
+        }
+    }
+
+    // Start receiving
+    pBufs.finishedSends();
+
+    // Consume
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& recvElems = map.constructMap()[domain];
+
+        if (domain != Pstream::myProcNo() && recvElems.size())
+        {
+            UIPstream str(domain, pBufs);
+
+            str >> points[domain]
+                >> nInternalFaces[domain]
+                >> faces[domain]
+                >> faceOwner[domain]
+                >> faceNeighbour[domain]
+                >> cellIDs[domain]
+                >> nbrProcIDs[domain]
+                >> procLocalFaceIDs[domain];
+        }
+
+        if (debug)
+        {
+            Pout<< "Target mesh send sizes[" << domain << "]"
+                << ": points="<< points[domain].size()
+                << ", faces=" << faces[domain].size()
+                << ", nInternalFaces=" << nInternalFaces[domain]
+                << ", faceOwn=" << faceOwner[domain].size()
+                << ", faceNbr=" << faceNeighbour[domain].size()
+                << ", cellIDs=" << cellIDs[domain].size() << endl;
+        }
+    }
+}
+
+
+void Foam::meshToMeshNew::distributeAndMergeCells
+(
+    const mapDistribute& map,
+    const polyMesh& tgt,
+    const globalIndex& globalI,
+    pointField& tgtPoints,
+    faceList& tgtFaces,
+    labelList& tgtFaceOwners,
+    labelList& tgtFaceNeighbours,
+    labelList& tgtCellIDs
+) const
+{
+    // Exchange per-processor data
+    List<pointField> allPoints;
+    List<label> allNInternalFaces;
+    List<faceList> allFaces;
+    List<labelList> allFaceOwners;
+    List<labelList> allFaceNeighbours;
+    List<labelList> allTgtCellIDs;
+
+    // Per target mesh face the neighbouring proc and index in
+    // processor patch (all -1 for normal boundary face)
+    List<labelList> allNbrProcIDs;
+    List<labelList> allProcLocalFaceIDs;
+
+    distributeCells
+    (
+        map,
+        tgt,
+        globalI,
+        allPoints,
+        allNInternalFaces,
+        allFaces,
+        allFaceOwners,
+        allFaceNeighbours,
+        allTgtCellIDs,
+        allNbrProcIDs,
+        allProcLocalFaceIDs
+    );
+
+    // Convert lists into format that can be used to generate a valid polyMesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    //
+    // Points and cells are collected into single flat lists:
+    // - i.e. proc0, proc1 ... procN
+    //
+    // Faces need to be sorted after collection to that internal faces are
+    // contiguous, followed by all boundary faces
+    //
+    // Processor patch faces between included cells on neighbouring processors
+    // are converted into internal faces
+    //
+    // Face list structure:
+    // - Per processor:
+    //   - internal faces
+    //   - processor faces that have been converted into internal faces
+    // - Followed by all boundary faces
+    //   - from 'normal' boundary faces
+    //   - from singularly-sided processor patch faces
+
+
+    // Number of internal+coupled faces
+    labelList allNIntCoupledFaces(allNInternalFaces);
+
+    // Starting offset for points
+    label nPoints = 0;
+    labelList pointOffset(Pstream::nProcs(), 0);
+    forAll(allPoints, procI)
+    {
+        pointOffset[procI] = nPoints;
+        nPoints += allPoints[procI].size();
+    }
+
+    // Starting offset for cells
+    label nCells = 0;
+    labelList cellOffset(Pstream::nProcs(), 0);
+    forAll(allTgtCellIDs, procI)
+    {
+        cellOffset[procI] = nCells;
+        nCells += allTgtCellIDs[procI].size();
+    }
+
+    // Count any coupled faces
+    typedef FixedList<label, 3> label3;
+    typedef HashTable<label, label3, label3::Hash<> > procCoupleInfo;
+    procCoupleInfo procFaceToGlobalCell;
+
+    forAll(allNbrProcIDs, procI)
+    {
+        const labelList& nbrProcI = allNbrProcIDs[procI];
+        const labelList& localFaceI = allProcLocalFaceIDs[procI];
+
+        forAll(nbrProcI, i)
+        {
+            if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+            {
+                label3 key;
+                key[0] = min(procI, nbrProcI[i]);
+                key[1] = max(procI, nbrProcI[i]);
+                key[2] = localFaceI[i];
+
+                procCoupleInfo::const_iterator fnd =
+                    procFaceToGlobalCell.find(key);
+
+                if (fnd == procFaceToGlobalCell.end())
+                {
+                    procFaceToGlobalCell.insert(key, -1);
+                }
+                else
+                {
+                    if (debug)
+                    {
+                        Pout<< "Additional internal face between procs:"
+                            << key[0] << " and " << key[1]
+                            << " across local face " << key[2] << endl;
+                    }
+
+                    allNIntCoupledFaces[procI]++;
+                }
+            }
+        }
+    }
+
+
+    // Starting offset for internal faces
+    label nIntFaces = 0;
+    label nFacesTotal = 0;
+    labelList internalFaceOffset(Pstream::nProcs(), 0);
+    forAll(allNIntCoupledFaces, procI)
+    {
+        label nCoupledFaces =
+            allNIntCoupledFaces[procI] - allNInternalFaces[procI];
+
+        internalFaceOffset[procI] = nIntFaces;
+        nIntFaces += allNIntCoupledFaces[procI];
+        nFacesTotal += allFaceOwners[procI].size() - nCoupledFaces;
+    }
+
+    tgtPoints.setSize(nPoints);
+    tgtFaces.setSize(nFacesTotal);
+    tgtFaceOwners.setSize(nFacesTotal);
+    tgtFaceNeighbours.setSize(nFacesTotal);
+    tgtCellIDs.setSize(nCells);
+
+    // Insert points
+    forAll(allPoints, procI)
+    {
+        const pointField& pts = allPoints[procI];
+        SubList<point>(tgtPoints, pts.size(), pointOffset[procI]).assign(pts);
+    }
+
+    // Insert cellIDs
+    forAll(allTgtCellIDs, procI)
+    {
+        const labelList& cellIDs = allTgtCellIDs[procI];
+        SubList<label>(tgtCellIDs, cellIDs.size(), cellOffset[procI]).assign
+        (
+            cellIDs
+        );
+    }
+
+
+    // Insert internal faces (from internal faces)
+    forAll(allFaces, procI)
+    {
+        const faceList& fcs = allFaces[procI];
+        const labelList& faceOs = allFaceOwners[procI];
+        const labelList& faceNs = allFaceNeighbours[procI];
+
+        SubList<face> slice
+        (
+            tgtFaces,
+            allNInternalFaces[procI],
+            internalFaceOffset[procI]
+        );
+        slice.assign(SubList<face>(fcs, allNInternalFaces[procI]));
+        forAll(slice, i)
+        {
+            add(slice[i], pointOffset[procI]);
+        }
+
+        SubField<label> ownSlice
+        (
+            tgtFaceOwners,
+            allNInternalFaces[procI],
+            internalFaceOffset[procI]
+        );
+        ownSlice.assign(SubField<label>(faceOs, allNInternalFaces[procI]));
+        add(ownSlice, cellOffset[procI]);
+
+        SubField<label> nbrSlice
+        (
+            tgtFaceNeighbours,
+            allNInternalFaces[procI],
+            internalFaceOffset[procI]
+        );
+        nbrSlice.assign(SubField<label>(faceNs, allNInternalFaces[procI]));
+        add(nbrSlice, cellOffset[procI]);
+
+        internalFaceOffset[procI] += allNInternalFaces[procI];
+    }
+
+
+    // Insert internal faces (from coupled face-pairs)
+    forAll(allNbrProcIDs, procI)
+    {
+        const labelList& nbrProcI = allNbrProcIDs[procI];
+        const labelList& localFaceI = allProcLocalFaceIDs[procI];
+        const labelList& faceOs = allFaceOwners[procI];
+        const faceList& fcs = allFaces[procI];
+
+        forAll(nbrProcI, i)
+        {
+            if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+            {
+                label3 key;
+                key[0] = min(procI, nbrProcI[i]);
+                key[1] = max(procI, nbrProcI[i]);
+                key[2] = localFaceI[i];
+
+                procCoupleInfo::iterator fnd = procFaceToGlobalCell.find(key);
+
+                if (fnd != procFaceToGlobalCell.end())
+                {
+                    label tgtFaceI = fnd();
+                    if (tgtFaceI == -1)
+                    {
+                        // on first visit store the new cell on this side
+                        fnd() = cellOffset[procI] + faceOs[i];
+                    }
+                    else
+                    {
+                        // get owner and neighbour in new cell numbering
+                        label newOwn = cellOffset[procI] + faceOs[i];
+                        label newNbr = fnd();
+                        label tgtFaceI = internalFaceOffset[procI]++;
+
+                        if (debug)
+                        {
+                            Pout<< "    proc " << procI
+                                << "\tinserting face:" << tgtFaceI
+                                << " connection between owner " << newOwn
+                                << " and neighbour " << newNbr
+                                << endl;
+                        }
+
+                        if (newOwn < newNbr)
+                        {
+                            // we have correct orientation
+                            tgtFaces[tgtFaceI] = fcs[i];
+                            tgtFaceOwners[tgtFaceI] = newOwn;
+                            tgtFaceNeighbours[tgtFaceI] = newNbr;
+                        }
+                        else
+                        {
+                            // reverse orientation
+                            tgtFaces[tgtFaceI] = fcs[i].reverseFace();
+                            tgtFaceOwners[tgtFaceI] = newNbr;
+                            tgtFaceNeighbours[tgtFaceI] = newOwn;
+                        }
+
+                        add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+                        // mark with unique value
+                        fnd() = -2;
+                    }
+                }
+            }
+        }
+    }
+
+
+    forAll(allNbrProcIDs, procI)
+    {
+        const labelList& nbrProcI = allNbrProcIDs[procI];
+        const labelList& localFaceI = allProcLocalFaceIDs[procI];
+        const labelList& faceOs = allFaceOwners[procI];
+        const labelList& faceNs = allFaceNeighbours[procI];
+        const faceList& fcs = allFaces[procI];
+
+        forAll(nbrProcI, i)
+        {
+            // coupled boundary face
+            if (nbrProcI[i] != -1 && localFaceI[i] != -1)
+            {
+                label3 key;
+                key[0] = min(procI, nbrProcI[i]);
+                key[1] = max(procI, nbrProcI[i]);
+                key[2] = localFaceI[i];
+
+                label tgtFaceI = procFaceToGlobalCell[key];
+
+                if (tgtFaceI == -1)
+                {
+                    FatalErrorIn
+                    (
+                        "void Foam::meshToMeshNew::"
+                        "distributeAndMergeCells"
+                        "("
+                            "const mapDistribute&, "
+                            "const polyMesh&, "
+                            "const globalIndex&, "
+                            "pointField&, "
+                            "faceList&, "
+                            "labelList&, "
+                            "labelList&, "
+                            "labelList&"
+                        ") const"
+                    )
+                        << "Unvisited " << key
+                        << abort(FatalError);
+                }
+                else if (tgtFaceI != -2)
+                {
+                    label newOwn = cellOffset[procI] + faceOs[i];
+                    label tgtFaceI = nIntFaces++;
+
+                    if (debug)
+                    {
+                        Pout<< "    proc " << procI
+                            << "\tinserting boundary face:" << tgtFaceI
+                            << " from coupled face " << key
+                            << endl;
+                    }
+
+                    tgtFaces[tgtFaceI] = fcs[i];
+                    add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+                    tgtFaceOwners[tgtFaceI] = newOwn;
+                    tgtFaceNeighbours[tgtFaceI] = -1;
+                }
+            }
+            // normal boundary face
+            else
+            {
+                label own = faceOs[i];
+                label nbr = faceNs[i];
+                if ((own != -1) && (nbr == -1))
+                {
+                    label newOwn = cellOffset[procI] + faceOs[i];
+                    label tgtFaceI = nIntFaces++;
+
+                    tgtFaces[tgtFaceI] = fcs[i];
+                    add(tgtFaces[tgtFaceI], pointOffset[procI]);
+
+                    tgtFaceOwners[tgtFaceI] = newOwn;
+                    tgtFaceNeighbours[tgtFaceI] = -1;
+                }
+            }
+        }
+    }
+
+
+    if (debug)
+    {
+        // only merging points in debug mode
+
+        labelList oldToNew;
+        pointField newTgtPoints;
+        bool hasMerged = mergePoints
+        (
+            tgtPoints,
+            SMALL,
+            false,
+            oldToNew,
+            newTgtPoints
+        );
+
+        if (hasMerged)
+        {
+            if (debug)
+            {
+                Pout<< "Merged from " << tgtPoints.size()
+                    << " down to " << newTgtPoints.size() << " points" << endl;
+            }
+
+            tgtPoints.transfer(newTgtPoints);
+            forAll(tgtFaces, i)
+            {
+                inplaceRenumber(oldToNew, tgtFaces[i]);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
new file mode 100644
index 0000000000..69bbc69667
--- /dev/null
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
@@ -0,0 +1,542 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMesh.H"
+#include "volFields.H"
+//#include "ops.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- Helper class for list
+    template<class Type>
+    class ListPlusEqOp
+    {
+        public:
+        void operator()(List<Type>& x, const List<Type> y) const
+        {
+            if (y.size())
+            {
+                if (x.size())
+                {
+                    label sz = x.size();
+                    x.setSize(sz + y.size());
+                    forAll(y, i)
+                    {
+                        x[sz++] = y[i];
+                    }
+                }
+                else
+                {
+                    x = y;
+                }
+            }
+        }
+    };
+
+    //- Combine operator for maps/interpolations
+    template<class Type, class CombineOp>
+    class combineBinaryOp
+    {
+        const CombineOp& cop_;
+
+        public:
+
+            combineBinaryOp(const CombineOp& cop)
+            :
+                cop_(cop)
+            {}
+
+            void operator()
+            (
+                Type& x,
+                const label faceI,
+                const Type& y,
+                const scalar weight
+            ) const
+            {
+                cop_(x, weight*y);
+            }
+    };
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::meshToMeshNew::add
+(
+    UList<Type>& fld,
+    const label offset
+) const
+{
+    forAll(fld, i)
+    {
+        fld[i] += offset;
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::mapSrcToTgt
+(
+    const UList<Type>& srcField,
+    const CombineOp& cop,
+    List<Type>& result
+) const
+{
+    if (result.size() != tgtToSrcCellAddr_.size())
+    {
+        FatalErrorIn
+        (
+            "void Foam::meshToMeshNew::mapSrcToTgt"
+            "("
+                "const UList<Type>&, "
+                "const CombineOp&, "
+                "List<Type>&"
+            ") const"
+        )   << "Supplied field size is not equal to target mesh size" << nl
+            << "    source mesh   = " << srcToTgtCellAddr_.size() << nl
+            << "    target mesh   = " << tgtToSrcCellAddr_.size() << nl
+            << "    supplied field = " << result.size()
+            << abort(FatalError);
+    }
+
+    combineBinaryOp<Type, CombineOp> cbop(cop);
+
+    if (singleMeshProc_ == -1)
+    {
+        const mapDistribute& map = srcMapPtr_();
+
+        List<Type> work(srcField);
+        map.distribute(work);
+
+        forAll(result, cellI)
+        {
+            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
+            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
+
+            if (srcAddress.size())
+            {
+//                result[cellI] = pTraits<Type>::zero;
+                result[cellI] *= (1.0 - sum(srcWeight));
+                forAll(srcAddress, i)
+                {
+                    label srcI = srcAddress[i];
+                    scalar w = srcWeight[i];
+                    cbop(result[cellI], cellI, work[srcI], w);
+                }
+            }
+        }
+    }
+    else
+    {
+        forAll(result, cellI)
+        {
+            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
+            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
+
+            if (srcAddress.size())
+            {
+//                result[cellI] = pTraits<Type>::zero;
+                result[cellI] *= (1.0 - sum(srcWeight));
+                forAll(srcAddress, i)
+                {
+                    label srcI = srcAddress[i];
+                    scalar w = srcWeight[i];
+                    cbop(result[cellI], cellI, srcField[srcI], w);
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+    const Field<Type>& srcField,
+    const CombineOp& cop
+) const
+{
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            tgtToSrcCellAddr_.size(),
+            pTraits<Type>::zero
+        )
+    );
+
+    mapSrcToTgt(srcField, cop, tresult());
+
+    return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+    const tmp<Field<Type> >& tsrcField,
+    const CombineOp& cop
+) const
+{
+    return mapSrcToTgt(tsrcField(), cop);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+    const Field<Type>& srcField
+) const
+{
+    return mapSrcToTgt(srcField, plusEqOp<Type>());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapSrcToTgt
+(
+    const tmp<Field<Type> >& tsrcField
+) const
+{
+    return mapSrcToTgt(tsrcField());
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::mapTgtToSrc
+(
+    const UList<Type>& tgtField,
+    const CombineOp& cop,
+    List<Type>& result
+) const
+{
+    if (result.size() != srcToTgtCellAddr_.size())
+    {
+        FatalErrorIn
+        (
+            "void Foam::meshToMeshNew::mapTgtToSrc"
+            "("
+                "const UList<Type>&, "
+                "const CombineOp&, "
+                "List<Type>&"
+            ") const"
+        )   << "Supplied field size is not equal to source mesh size" << nl
+            << "    source mesh   = " << srcToTgtCellAddr_.size() << nl
+            << "    target mesh   = " << tgtToSrcCellAddr_.size() << nl
+            << "    supplied field = " << result.size()
+            << abort(FatalError);
+    }
+
+    combineBinaryOp<Type, CombineOp> cbop(cop);
+
+    if (singleMeshProc_ == -1)
+    {
+        const mapDistribute& map = tgtMapPtr_();
+
+        List<Type> work(tgtField);
+        map.distribute(work);
+
+        forAll(result, cellI)
+        {
+            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
+            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
+
+            if (tgtAddress.size())
+            {
+//                result[cellI] = pTraits<Type>::zero;
+                result[cellI] *= (1.0 - sum(tgtWeight));
+                forAll(tgtAddress, i)
+                {
+                    label tgtI = tgtAddress[i];
+                    scalar w = tgtWeight[i];
+                    cbop(result[cellI], cellI, work[tgtI], w);
+                }
+            }
+        }
+    }
+    else
+    {
+        forAll(result, cellI)
+        {
+            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
+            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
+
+            if (tgtAddress.size())
+            {
+//                result[cellI] = pTraits<Type>::zero;
+                result[cellI] *= (1.0 - sum(tgtWeight));
+                forAll(tgtAddress, i)
+                {
+                    label tgtI = tgtAddress[i];
+                    scalar w = tgtWeight[i];
+                    cbop(result[cellI], cellI, tgtField[tgtI], w);
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+    const Field<Type>& tgtField,
+    const CombineOp& cop
+) const
+{
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            srcToTgtCellAddr_.size(),
+            pTraits<Type>::zero
+        )
+    );
+
+    mapTgtToSrc(tgtField, cop, tresult());
+
+    return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+    const tmp<Field<Type> >& ttgtField,
+    const CombineOp& cop
+) const
+{
+    return mapTgtToSrc(ttgtField(), cop);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+    const Field<Type>& tgtField
+) const
+{
+    return mapTgtToSrc(tgtField, plusEqOp<Type>());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::meshToMeshNew::mapTgtToSrc
+(
+    const tmp<Field<Type> >& ttgtField
+) const
+{
+    return mapTgtToSrc(ttgtField(), plusEqOp<Type>());
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMeshNew::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const CombineOp& cop,
+    GeometricField<Type, fvPatchField, volMesh>& result,
+    const bool interpPatches
+) const
+{
+    const fvMesh& mesh = field.mesh();
+
+    if (mesh.name() == srcRegionName_)
+    {
+        mapSrcToTgt(field, cop, result.internalField());
+    }
+    else if (mesh.name() == tgtRegionName_)
+    {
+        mapTgtToSrc(field, cop, result.internalField());
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::meshToMeshNew::interpolate"
+            "("
+                "const GeometricField<Type, fvPatchField, volMesh>&, "
+                "const CombineOp&, "
+                "GeometricField<Type, fvPatchField, volMesh>&, "
+                "const bool"
+            ") const"
+        )
+            << "Supplied field " << field.name() << " did not originate from "
+            << "either the source or target meshes used to create this "
+            << "interpolation object"
+            << abort(FatalError);
+    }
+
+    if (interpPatches)
+    {
+        switch (method_)
+        {
+            case imMap:
+            {
+                result.boundaryField() == field.boundaryField();
+                break;
+            }
+            default:
+            {
+                notImplemented
+                (
+                    "void Foam::meshToMeshNew::interpolate"
+                    "("
+                        "const GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const CombineOp&, "
+                        "GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const bool"
+                    ") const - non-conformal patches"
+                )
+
+                // do something...
+            }
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const CombineOp& cop,
+    const bool interpPatches
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const fvMesh& mesh = field.mesh();
+
+    tmp<fieldType> tresult;
+
+    if (mesh.name() == srcRegionName_)
+    {
+        const fvMesh& tgtMesh =
+            mesh.time().lookupObject<fvMesh>(tgtRegionName_);
+
+        tresult =
+            new fieldType
+            (
+                IOobject
+                (
+                    type() + "::interpolate(" + field.name() + ")",
+                    tgtMesh.time().timeName(),
+                    tgtMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                tgtMesh,
+                dimensioned<Type>
+                (
+                    "zero",
+                    field.dimensions(),
+                    pTraits<Type>::zero
+                )
+            );
+
+         interpolate(field, cop, tresult(), interpPatches);
+    }
+    else if (mesh.name() == tgtRegionName_)
+    {
+        const fvMesh& srcMesh =
+            mesh.time().lookupObject<fvMesh>(srcRegionName_);
+
+        tresult =
+            new fieldType
+            (
+                IOobject
+                (
+                    type() + "::interpolate(" + field.name() + ")",
+                    srcMesh.time().timeName(),
+                    srcMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                srcMesh,
+                dimensioned<Type>
+                (
+                    "zero",
+                    field.dimensions(),
+                    pTraits<Type>::zero
+                )
+            );
+
+         interpolate(field, cop, tresult(), interpPatches);
+    }
+
+    return tresult;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
+    const CombineOp& cop,
+    const bool interpPatches
+) const
+{
+    return
+        interpolate
+        (
+            tfield(),
+            combineBinaryOp<Type, CombineOp>(cop),
+            interpPatches
+        );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const bool interpPatches
+) const
+{
+    return interpolate(field, plusEqOp<Type>(), interpPatches);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMeshNew::interpolate
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
+    const bool interpPatches
+) const
+{
+    return interpolate(tfield(), plusEqOp<Type>(), interpPatches);
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 3edfd4ce38..08dd1e3764 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -312,7 +312,30 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
     {
         const volScalarField& dQ =
             mesh_.lookupObject<volScalarField>("dQ");
-        E().internalField() = EhrrCoeff_*dQ;
+
+        if (dQ.dimensions() == dimEnergy/dimTime)
+        {
+            E().internalField() = EhrrCoeff_*dQ/mesh_.V();
+        }
+        else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
+        {
+            E().internalField() = EhrrCoeff_*dQ;
+        }
+        else
+        {
+            if (debug)
+            {
+                WarningIn
+                (
+                    "tmp<volScalarField>"
+                    "radiation::greyMeanAbsorptionEmission::ECont"
+                    "("
+                        "const label"
+                    ") const"
+                )
+                    << "Incompatible dimensions for dQ field" << endl;
+            }
+        }
     }
 
     return E;
diff --git a/tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvOptions b/tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvOptions
index 5b2d93ccfa..246601aa16 100644
--- a/tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvOptions
+++ b/tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvOptions b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvOptions
index f0a9fed5d8..4f0ffe9c1b 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvOptions
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvOptions b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvOptions
index 9f8e69357a..d75a43c565 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvOptions
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/T b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/T
new file mode 100644
index 0000000000..fcb0067f35
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/T
@@ -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       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 300;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+    blades
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/U b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/U
new file mode 100644
index 0000000000..a5466e3636
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/U
@@ -0,0 +1,48 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 0 5);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+    blades
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/epsilon b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/epsilon
new file mode 100644
index 0000000000..46cfcab3a2
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/epsilon
@@ -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       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.54;
+
+boundaryField
+{
+    walls
+    {
+        type            compressible::epsilonWallFunction;
+        value           $internalField;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    blades
+    {
+        type            compressible::epsilonWallFunction;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/k b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/k
new file mode 100644
index 0000000000..6a2ab6f940
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/k
@@ -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       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.375;
+
+boundaryField
+{
+    walls
+    {
+        type            compressible::kqRWallFunction;
+        value           $internalField;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    blades
+    {
+        type            compressible::kqRWallFunction;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p
new file mode 100644
index 0000000000..f20dd774ef
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p
@@ -0,0 +1,47 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    walls
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    blades
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p_rgh b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p_rgh
new file mode 100644
index 0000000000..7561126b38
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/air/p_rgh
@@ -0,0 +1,49 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    walls
+    {
+        type            buoyantPressure;
+        gradient        uniform 0;
+        value           uniform 100000;
+    }
+    inlet
+    {
+        type            buoyantPressure;
+        gradient        uniform 0;
+        value           uniform 100000;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+    blades
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/AoV b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/AoV
new file mode 100644
index 0000000000..0791bb96be
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/AoV
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      AoV;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 -1 0 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    ".*"
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/T b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/T
new file mode 100644
index 0000000000..bd86a29732
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/T
@@ -0,0 +1,42 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 350;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 400;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 350;
+        value           uniform 350;
+    }
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U
new file mode 100644
index 0000000000..23454ad6de
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- 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       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.01 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0.01 0 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+    fixedWalls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/htcConst b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/htcConst
new file mode 100644
index 0000000000..4f3040e85a
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/htcConst
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      htcConst;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 -1 0 0 0];
+
+internalField   uniform 10;
+
+boundaryField
+{
+    ".*"
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p
new file mode 100644
index 0000000000..db0ae810f9
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p
@@ -0,0 +1,42 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    fixedWalls
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p_rgh b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p_rgh
new file mode 100644
index 0000000000..6b4bc21fb3
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/p_rgh
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    inlet
+    {
+        type            buoyantPressure;
+        gradient        uniform 0;
+        value           uniform 100000;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+    fixedWalls
+    {
+        type            buoyantPressure;
+        gradient        uniform 0;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allclean b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allclean
new file mode 100755
index 0000000000..a915451689
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allclean
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+foamClearPolyMesh -region air
+foamClearPolyMesh -region porous
+
+rm -f *.OpenFOAM
+
+rm -rf 0
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-parallel b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-parallel
new file mode 100755
index 0000000000..6c26ff7801
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-parallel
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+./Allrun.pre
+
+runApplication decomposePar -region air
+mv log.decomposePar log.decomposePar.air
+
+runApplication decomposePar -region porous
+mv log.decomposePar log.decomposePar.porous
+
+runParallel $(getApplication) 4
+
+runApplication reconstructPar -latestTime -region air
+mv log.reconstructPar log.reconstructPar.air
+
+runApplication reconstructPar -latestTime -region porous
+mv log.reconstructPar log.reconstructPar.porous
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-serial b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-serial
new file mode 100755
index 0000000000..c01a61be81
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-serial
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+cd ${0%/*} || exit 1    # run from this directory
+
+./Allrun.pre
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication $(getApplication)
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun.pre b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun.pre
new file mode 100755
index 0000000000..99db2af1b6
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun.pre
@@ -0,0 +1,31 @@
+#!/bin/bash
+
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+#create meshes
+runApplication blockMesh -region air
+mv log.blockMesh log.blockMesh.air
+runApplication blockMesh -region porous
+mv log.blockMesh log.blockMesh.porous
+
+# create rotor blades in air region
+runApplication topoSet -region air -dict system/topoSetDict.1
+mv log.topoSet log.topoSet.air.1
+
+runApplication createBaffles -region air -overwrite
+
+# create rotor zone in air region for MRF
+runApplication topoSet -region air -dict system/topoSetDict.2
+mv log.topoSet log.topoSet.air.2
+
+rm -rf constant/air/polyMesh/sets
+
+# create dummy files for post-processing
+paraFoam -touch -region porous
+paraFoam -touch -region air
+
+
+cp -rf 0.org 0
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/RASProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/RASProperties
new file mode 100644
index 0000000000..a4937b503a
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/g b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/g
new file mode 100644
index 0000000000..6343b983a8
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/g
@@ -0,0 +1,23 @@
+/*--------------------------------*- 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       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+//value           ( 0 0 9.81);
+value           ( 0 0 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..d0439108c7
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict
@@ -0,0 +1,149 @@
+/*--------------------------------*- 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       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.05;
+
+vertices
+(
+    (0 0 0)             // 0
+    (10 0 0)            // 1
+    (10 10 0)           // 2
+    (0 10 0)            // 3
+
+    (2.879 2.879 0)     // 4
+    (7.121 2.879 0)     // 5
+    (7.121 7.121 0)     // 6
+    (2.879 7.121 0)     // 7
+
+    (4.293 4.293 0)     // 8
+    (5.707 4.293 0)     // 9
+    (5.707 5.707 0)     // 10
+    (4.293 5.707 0)     // 11
+
+
+    (0 0 10)            // 12
+    (10 0 10)           // 13
+    (10 10 10)          // 14
+    (0 10 10)           // 15
+
+    (2.879 2.879 10)    // 16
+    (7.121 2.879 10)    // 17
+    (7.121 7.121 10)    // 18
+    (2.879 7.121 10)    // 19
+
+    (4.293 4.293 10)    // 20
+    (5.707 4.293 10)    // 20
+    (5.707 5.707 10)    // 22
+    (4.293 5.707 10)    // 23
+);
+
+blocks
+(
+    hex (0 1  5  4 12 13 17 16) (30 30 30) simpleGrading (1 1 1)
+    hex (1 2  6  5 13 14 18 17) (30 30 30) simpleGrading (1 1 1)
+    hex (2 3  7  6 14 15 19 18) (30 30 30) simpleGrading (1 1 1)
+    hex (3 0  4  7 15 12 16 19) (30 30 30) simpleGrading (1 1 1)
+    hex (4 5  9  8 16 17 21 20) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (5 6 10  9 17 18 22 21) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (6 7 11 10 18 19 23 22) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (7 4  8 11 19 16 20 23) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (8 9 10 11 20 21 22 23) innerCylinder (30 30 30) simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc 4 5 (5 2 0)
+    arc 5 6 (8 5 0)
+    arc 6 7 (5 8 0)
+    arc 7 4 (2 5 0)
+
+    arc 8 9 (5 4 0)
+    arc 9 10 (6 5 0)
+    arc 10 11 (5 6 0)
+    arc 11 8 (4 5 0)
+
+    arc 16 17 (5 2 10)
+    arc 17 18 (8 5 10)
+    arc 18 19 (5 8 10)
+    arc 19 16 (2 5 10)
+
+    arc 20 21 (5 4 10)
+    arc 21 22 (6 5 10)
+    arc 22 23 (5 6 10)
+    arc 23 20 (4 5 10)
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 5 1)
+            (1 5 6 2)
+            (2 6 7 3)
+            (3 7 4 0)
+
+            (12 13 17 16)
+            (13 14 18 17)
+            (14 15 19 18)
+            (15 12 16 19)
+
+            (0 12 13 1)
+            (1 13 14 2)
+            (2 14 15 3)
+            (3 15 12 0)
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (4 8 9 5)
+            (5 9 10 6)
+            (6 10 11 7)
+            (7 11 8 4)
+            (8 11 10 9)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (16 17 21 20)
+            (17 18 22 21)
+            (18 19 23 22)
+            (19 16 20 23)
+            (20 21 22 23)
+        );
+    }
+
+    blades
+    {
+        type wall;
+        faces
+        ();
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/radiationProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/radiationProperties
new file mode 100644
index 0000000000..aa0a16678f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/radiationProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/thermophysicalProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/thermophysicalProperties
new file mode 100644
index 0000000000..b3c0a58998
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/thermophysicalProperties
@@ -0,0 +1,55 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       polynomial;
+    thermo          hPolynomial;
+    equationOfState icoPolynomial;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    // coefficients for air
+
+    specie
+    {
+        nMoles          1;
+        molWeight       28.85;
+    }
+    equationOfState
+    {
+        rhoCoeffs<8>    ( 4.0097 -0.016954 3.3057e-05 -3.0042e-08 1.0286e-11 0 0 0 );
+    }
+    thermodynamics
+    {
+        Hf              0;
+        Sf              0;
+        CpCoeffs<8>     ( 948.76 0.39171 -0.00095999 1.393e-06 -6.2029e-10 0 0 0 );
+    }
+    transport
+    {
+        muCoeffs<8>     ( 1.5061e-06 6.16e-08 -1.819e-11 0 0 0 0 0 );
+        kappaCoeffs<8>  ( 0.0025219 8.506e-05 -1.312e-08 0 0 0 0 0 );
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/turbulenceProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/turbulenceProperties
new file mode 100644
index 0000000000..3721a46a2e
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/RASProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/RASProperties
new file mode 100644
index 0000000000..1cc42b8af8
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/RASProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- 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       dictionary;
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel laminar;
+
+turbulence      off;
+
+printCoeffs     off;
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/g b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/g
new file mode 100644
index 0000000000..6343b983a8
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/g
@@ -0,0 +1,23 @@
+/*--------------------------------*- 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       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+//value           ( 0 0 9.81);
+value           ( 0 0 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..eae25e8926
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict
@@ -0,0 +1,76 @@
+/*--------------------------------*- 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       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.05;
+
+vertices
+(
+    (-2 4 4)     // 0
+    (12 4 4)     // 1
+    (12 6 4)     // 2
+    (-2 6 4)     // 3
+
+    (-2 4 6)     // 4
+    (12 4 6)     // 5
+    (12 6 6)     // 6
+    (-2 6 6)     // 7
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (40 15 15) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (1 2 6 5)
+        );
+    }
+    fixedWalls
+    {
+        type wall;
+        faces
+        (
+            (0 1 2 3)
+            (4 5 6 7)
+            (3 2 6 7)
+            (0 1 5 4)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/thermophysicalProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/thermophysicalProperties
new file mode 100644
index 0000000000..41e355a66b
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/thermophysicalProperties
@@ -0,0 +1,55 @@
+/*--------------------------------*- 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       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       polynomial;
+    thermo          hPolynomial;
+    equationOfState icoPolynomial;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+
+mixture
+{
+    // coefficients for water
+
+    specie
+    {
+        nMoles          1;
+        molWeight       18;
+    }
+    equationOfState
+    {
+        rhoCoeffs<8>    ( 1000 0 0 0 0 0 0 0 );
+    }
+    thermodynamics
+    {
+        Hf              0;
+        Sf              0;
+        CpCoeffs<8>     ( 4183 0 0 0 0 0 0 0 );
+    }
+    transport
+    {
+        muCoeffs<8>     ( 0.001 0 0 0 0 0 0 0 );
+        kappaCoeffs<8>  ( 0.58  0 0 0 0 0 0 0 );
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/turbulenceProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/turbulenceProperties
new file mode 100644
index 0000000000..e63bbc5081
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/turbulenceProperties
@@ -0,0 +1,19 @@
+/*--------------------------------*- 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       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/regionProperties b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/regionProperties
new file mode 100644
index 0000000000..749edd93d5
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/regionProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      regionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+regions
+(
+    fluid   (air porous)
+    solid   ()
+);
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/createBafflesDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/createBafflesDict
new file mode 100644
index 0000000000..818b3b2320
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/createBafflesDict
@@ -0,0 +1,52 @@
+/*--------------------------------*- 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       dictionary;
+    object      createBafflesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Whether to convert internal faces only (so leave boundary faces intact).
+// This is only relevant if your face selection type can pick up boundary
+// faces.
+internalFacesOnly true;
+
+
+// Baffles to create.
+baffles
+{
+    baffleFaces
+    {
+        //- Use predefined faceZone to select faces and orientation.
+        type        faceZone;
+        zoneName    rotorBlades;
+
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+                name            blades;
+                type            patch;
+            }
+            slave
+            {
+                //- Slave side patch
+                name            blades;
+                type            patch;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/decomposeParDict
new file mode 100644
index 0000000000..645728dc50
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions
new file mode 100644
index 0000000000..44aba3100f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions
@@ -0,0 +1,50 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+airToporous
+{
+    type            constantHeatTransfer;
+    active          on;
+    selectionMode   mapRegion;
+    interpolationMethod cellVolumeWeight;
+    nbrModelName    porousToair;
+    nbrRegionName   porous;
+    master          false;
+
+    constantHeatTransferCoeffs
+    {
+        fieldNames      (h);
+        semiImplicit    no;
+    }
+}
+
+MRF1
+{
+    type            MRFSource;
+    active          true;
+    selectionMode   cellZone;
+    cellZone        rotor;
+
+    MRFSourceCoeffs
+    {
+        origin    (0.25 0.25 0.25);
+        axis      (0 0 1);
+        omega     5.305; // 500 rpm
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSchemes
new file mode 100644
index 0000000000..a39397c7c3
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSchemes
@@ -0,0 +1,63 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(U)         cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,e)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited 0.333;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited 0.333;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution
new file mode 100644
index 0000000000..36d404271f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution
@@ -0,0 +1,70 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p_rgh
+    {
+        solver           GAMG;
+        tolerance        1e-7;
+        relTol           0.01;
+
+        smoother         DIC;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          100;
+    }
+
+    "(U|h|e|k|epsilon)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-6;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    rhoMax          rhoMax [ 1 -3 0 0 0 ] 2;
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 1;
+}
+
+
+relaxationFactors
+{
+    fields
+    {
+        rho             1;
+        p_rgh           0.7;
+    }
+    equations
+    {
+        U               0.3;
+        "(h|e)"         0.3;
+        k               0.3;
+        epsilon         0.3;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.1 b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.1
new file mode 100644
index 0000000000..36dc18f95f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.1
@@ -0,0 +1,120 @@
+/*--------------------------------*- 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       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // rotor faces
+    {
+        name    rotorCells;
+        type    cellSet;
+        action  new;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name    cylinder;
+        }
+    }
+    {
+        name    rotorFaces;
+        type    faceSet;
+        action  new;
+        source  cellToFace;
+        sourceInfo
+        {
+            set     rotorCells;
+            option  all;
+        }
+    }
+    {
+        name    rotorFaces;
+        type    faceSet;
+        action  subset;
+        source  boxToFace;
+        sourceInfo
+        {
+            box (-100 -100 0.1) (100 100 0.15);
+        }
+    }
+
+    {
+        name    rotorBlades;
+        type    faceSet;
+        action  new;
+        source  faceToFace;
+        sourceInfo
+        {
+            set     rotorFaces;
+        }
+    }
+
+    {
+        name    rotorBlades;
+        type    faceSet;
+        action  subset;
+        source  boxToFace;
+        sourceInfo
+        {
+            box (-100 0.249 -100) (100 0.251 100);
+        }
+    }
+
+    {
+        name    rotorBlades2;
+        type    faceSet;
+        action  new;
+        source  faceToFace;
+        sourceInfo
+        {
+            set     rotorFaces;
+        }
+    }
+
+    {
+        name    rotorBlades2;
+        type    faceSet;
+        action  subset;
+        source  boxToFace;
+        sourceInfo
+        {
+            box (0.249 -100 -100) (0.251 100 100);
+        }
+    }
+
+    {
+        name    rotorBlades;
+        type    faceSet;
+        action  add;
+        source  faceToFace;
+        sourceInfo
+        {
+            set     rotorBlades2;
+        }
+    }
+    {
+        name    rotorBlades;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        sourceInfo
+        {
+            faceSet rotorBlades;
+        }
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.2 b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.2
new file mode 100644
index 0000000000..f9f33b1095
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/topoSetDict.2
@@ -0,0 +1,64 @@
+/*--------------------------------*- 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       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // rotor cell zone
+    {
+        name    rotorCells;
+        type    cellSet;
+        action  new;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name    cylinder;
+        }
+    }
+    {
+        name    rotorCells;
+        type    cellSet;
+        action  add;
+        source  zoneToCell;
+        sourceInfo
+        {
+            name    innerCylinder;
+        }
+    }
+    {
+        name    rotorCells;
+        type    cellSet;
+        action  subset;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-100 -100 0.1) (100 100 0.15);
+        }
+    }
+    {
+        name    rotor;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set rotorCells;
+        }
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/controlDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/controlDict
new file mode 100644
index 0000000000..afc00df3f2
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/controlDict
@@ -0,0 +1,49 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     chtMultiRegionSimpleFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/decomposeParDict
new file mode 100644
index 0000000000..645728dc50
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSchemes
new file mode 100644
index 0000000000..e7d321e959
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSchemes
@@ -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       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSolution
new file mode 100644
index 0000000000..cde9e49189
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/fvSolution
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/decomposeParDict
new file mode 100644
index 0000000000..e7e490bf74
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions
new file mode 100644
index 0000000000..64561c37f6
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions
@@ -0,0 +1,35 @@
+/*--------------------------------*- 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       dictionary;
+    location    "system";
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+porousToair
+{
+    type            constantHeatTransfer;
+    active          on;
+    selectionMode   mapRegion;
+    interpolationMethod cellVolumeWeight;
+    nbrModelName    airToporous;
+    nbrRegionName   air;
+    master          true;
+
+    constantHeatTransferCoeffs
+    {
+        fieldNames      (h);
+        semiImplicit    no;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSchemes
new file mode 100644
index 0000000000..f787aa1bb4
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSchemes
@@ -0,0 +1,61 @@
+/*--------------------------------*- 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       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,e)      bounded Gauss upwind;
+    div(phi,Ekp)    bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution
new file mode 100644
index 0000000000..fd6da4016c
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution
@@ -0,0 +1,68 @@
+/*--------------------------------*- 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       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p_rgh
+    {
+        solver           GAMG;
+        tolerance        1e-7;
+        relTol           0.01;
+
+        smoother         DIC;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          100;
+    }
+
+    "(U|h|e|k|epsilon)"
+    {
+        solver           PBiCG;
+        preconditioner   DILU;
+        tolerance        1e-7;
+        relTol           0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    rhoMax          rhoMax [ 1 -3 0 0 0 ] 1100;
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 900;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        rho             1;
+        p_rgh           0.7;
+    }
+    equations
+    {
+        U               0.3;
+        "(h|e)"         0.3;
+        k               0.3;
+        epsilon         0.3;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/channel395/system/fvOptions b/tutorials/incompressible/pimpleFoam/channel395/system/fvOptions
index 97bf730df9..681d7e39ff 100644
--- a/tutorials/incompressible/pimpleFoam/channel395/system/fvOptions
+++ b/tutorials/incompressible/pimpleFoam/channel395/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/simpleFoam/mixerVessel2D/system/fvOptions b/tutorials/incompressible/simpleFoam/mixerVessel2D/system/fvOptions
index 594e280853..f55734732b 100644
--- a/tutorials/incompressible/simpleFoam/mixerVessel2D/system/fvOptions
+++ b/tutorials/incompressible/simpleFoam/mixerVessel2D/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvOptions b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvOptions
index a09b67dc9f..61e7340025 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvOptions
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvOptions b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvOptions
index ee4efec416..82e5c6b4b0 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvOptions
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvOptions b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvOptions
index ee4efec416..82e5c6b4b0 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvOptions
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvOptions b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvOptions
index ef873bfa25..f3194e5954 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvOptions
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/system/fvOptions b/tutorials/lagrangian/reactingParcelFoam/filter/system/fvOptions
index 82c4f9ed19..c7dd4992aa 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/system/fvOptions
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/fvOptions b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/fvOptions
index ee4efec416..82e5c6b4b0 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/fvOptions
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/fvOptions
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "system";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //