diff --git a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
new file mode 100644
index 000000000..672573e8c
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MapConsistentVolFields_H
+#define MapConsistentVolFields_H
+
+#include "GeometricField.H"
+#include "meshToMesh0.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type, class CombineOp>
+void MapConsistentVolFields
+(
+    const IOobjectList& objects,
+    const meshToMesh0& meshToMesh0Interp,
+    const meshToMesh0::order& mapOrder,
+    const CombineOp& cop
+)
+{
+    const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
+    const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
+
+    word fieldClassName
+    (
+        GeometricField<Type, fvPatchField, volMesh>::typeName
+    );
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    forAllIter(IOobjectList, fields, fieldIter)
+    {
+        Info<< "    interpolating " << fieldIter()->name()
+            << endl;
+
+        // Read field
+        GeometricField<Type, fvPatchField, volMesh> fieldSource
+        (
+            *fieldIter(),
+            meshSource
+        );
+
+        IOobject fieldTargetIOobject
+        (
+            fieldIter()->name(),
+            meshTarget.time().timeName(),
+            meshTarget,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        );
+
+        if (fieldTargetIOobject.headerOk())
+        {
+            // Read fieldTarget
+            GeometricField<Type, fvPatchField, volMesh> fieldTarget
+            (
+                fieldTargetIOobject,
+                meshTarget
+            );
+
+            // Interpolate field
+            meshToMesh0Interp.interpolate
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
+
+            // Write field
+            fieldTarget.write();
+        }
+        else
+        {
+            fieldTargetIOobject.readOpt() = IOobject::NO_READ;
+
+            // Interpolate field
+            GeometricField<Type, fvPatchField, volMesh> fieldTarget
+            (
+                fieldTargetIOobject,
+                meshToMesh0Interp.interpolate
+                (
+                    fieldSource,
+                    mapOrder,
+                    cop
+                )
+            );
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFields/MapLagrangianFields.H b/applications/utilities/preProcessing/mapFields/MapLagrangianFields.H
index 07d8d8335..9244091d2 100644
--- a/applications/utilities/preProcessing/mapFields/MapLagrangianFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapLagrangianFields.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ Description
 
 #include "cloud.H"
 #include "GeometricField.H"
-#include "meshToMesh.H"
+#include "meshToMesh0.H"
 #include "IOobjectList.H"
 #include "CompactIOField.H"
 
@@ -51,18 +51,19 @@ void MapLagrangianFields
 (
     const string& cloudName,
     const IOobjectList& objects,
-    const polyMesh& meshTarget,
+    const meshToMesh0& meshToMesh0Interp,
     const labelList& addParticles
 )
 {
+    const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
+
     {
         IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
 
         forAllIter(IOobjectList, fields, fieldIter)
         {
-            const word& fieldName = fieldIter()->name();
-
-            Info<< "    mapping lagrangian field " << fieldName << endl;
+            Info<< "    mapping lagrangian field "
+                << fieldIter()->name() << endl;
 
             // Read field (does not need mesh)
             IOField<Type> fieldSource(*fieldIter());
@@ -72,7 +73,7 @@ void MapLagrangianFields
             (
                 IOobject
                 (
-                    fieldName,
+                    fieldIter()->name(),
                     meshTarget.time().timeName(),
                     cloud::prefix/cloudName,
                     meshTarget,
@@ -99,9 +100,8 @@ void MapLagrangianFields
 
         forAllIter(IOobjectList, fieldFields, fieldIter)
         {
-            const word& fieldName = fieldIter()->name();
-
-            Info<< "    mapping lagrangian fieldField " << fieldName << endl;
+            Info<< "    mapping lagrangian fieldField "
+                << fieldIter()->name() << endl;
 
             // Read field (does not need mesh)
             IOField<Field<Type> > fieldSource(*fieldIter());
@@ -112,7 +112,7 @@ void MapLagrangianFields
             (
                 IOobject
                 (
-                    fieldName,
+                    fieldIter()->name(),
                     meshTarget.time().timeName(),
                     cloud::prefix/cloudName,
                     meshTarget,
diff --git a/applications/utilities/preProcessing/mapFields/MapMeshes.H b/applications/utilities/preProcessing/mapFields/MapMeshes.H
index e75819f28..b11138da6 100644
--- a/applications/utilities/preProcessing/mapFields/MapMeshes.H
+++ b/applications/utilities/preProcessing/mapFields/MapMeshes.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,9 +27,9 @@ License
 #define MapMeshes_H
 
 #include "MapVolFields.H"
+#include "MapConsistentVolFields.H"
 #include "mapLagrangian.H"
 #include "UnMapped.H"
-#include "pointMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,62 +37,64 @@ namespace Foam
 {
 
 template<template<class> class CombineOp>
-void MapMesh
+void MapConsistentMesh
 (
-    const meshToMesh& interp,
-    const HashSet<word>& selectedFields,
-    const bool noLagrangian
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh0::order& mapOrder
 )
 {
-    {
-        const polyMesh& meshSource = interp.srcRegion();
+    // Create the interpolation scheme
+    meshToMesh0 meshToMesh0Interp(meshSource, meshTarget);
 
+    Info<< nl
+        << "Consistently creating and mapping fields for time "
+        << meshSource.time().timeName() << nl << endl;
+
+    {
         // Search for list of objects for this time
         IOobjectList objects(meshSource, meshSource.time().timeName());
 
         // Map volFields
         // ~~~~~~~~~~~~~
-        MapVolFields<scalar>
+        MapConsistentVolFields<scalar>
         (
             objects,
-            selectedFields,
-            interp,
+            meshToMesh0Interp,
+            mapOrder,
             CombineOp<scalar>()
         );
-
-        MapVolFields<vector>
+        MapConsistentVolFields<vector>
         (
             objects,
-            selectedFields,
-            interp,
+            meshToMesh0Interp,
+            mapOrder,
             CombineOp<vector>()
         );
-        MapVolFields<sphericalTensor>
+        MapConsistentVolFields<sphericalTensor>
         (
             objects,
-            selectedFields,
-            interp,
+            meshToMesh0Interp,
+            mapOrder,
             CombineOp<sphericalTensor>()
         );
-        MapVolFields<symmTensor>
+        MapConsistentVolFields<symmTensor>
         (
             objects,
-            selectedFields,
-            interp,
+            meshToMesh0Interp,
+            mapOrder,
             CombineOp<symmTensor>()
         );
-        MapVolFields<tensor>
+        MapConsistentVolFields<tensor>
         (
             objects,
-            selectedFields,
-            interp,
+            meshToMesh0Interp,
+            mapOrder,
             CombineOp<tensor>()
         );
     }
 
     {
-        const polyMesh& meshTarget = interp.tgtRegion();
-
         // Search for list of target objects for this time
         IOobjectList objects(meshTarget, meshTarget.time().timeName());
 
@@ -113,10 +115,140 @@ void MapMesh
         UnMapped<pointTensorField>(objects);
     }
 
-    if (!noLagrangian)
+    mapLagrangian(meshToMesh0Interp);
+}
+
+
+template<template<class> class CombineOp>
+void MapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const meshToMesh0::order& mapOrder
+)
+{
+    // Create the interpolation scheme
+    meshToMesh0 meshToMesh0Interp
+    (
+        meshSource,
+        meshTarget,
+        patchMap,
+        cuttingPatches
+    );
+
+    Info<< nl
+        << "Mapping fields for time " << meshSource.time().timeName()
+        << nl << endl;
+
     {
-        mapLagrangian(interp);
+        // Search for list of source objects for this time
+        IOobjectList objects(meshSource, meshSource.time().timeName());
+
+        // Map volFields
+        // ~~~~~~~~~~~~~
+        MapVolFields<scalar>
+        (
+            objects,
+            meshToMesh0Interp,
+            mapOrder,
+            CombineOp<scalar>()
+        );
+        MapVolFields<vector>
+        (
+            objects,
+            meshToMesh0Interp,
+            mapOrder,
+            CombineOp<vector>()
+        );
+        MapVolFields<sphericalTensor>
+        (
+            objects,
+            meshToMesh0Interp,
+            mapOrder,
+            CombineOp<sphericalTensor>()
+        );
+        MapVolFields<symmTensor>
+        (
+            objects,
+            meshToMesh0Interp,
+            mapOrder,
+            CombineOp<symmTensor>()
+        );
+        MapVolFields<tensor>
+        (
+            objects,
+            meshToMesh0Interp,
+            mapOrder,
+            CombineOp<tensor>()
+        );
     }
+
+    {
+        // Search for list of target objects for this time
+        IOobjectList objects(meshTarget, meshTarget.time().timeName());
+
+        // Mark surfaceFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<surfaceScalarField>(objects);
+        UnMapped<surfaceVectorField>(objects);
+        UnMapped<surfaceSphericalTensorField>(objects);
+        UnMapped<surfaceSymmTensorField>(objects);
+        UnMapped<surfaceTensorField>(objects);
+
+        // Mark pointFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<pointScalarField>(objects);
+        UnMapped<pointVectorField>(objects);
+        UnMapped<pointSphericalTensorField>(objects);
+        UnMapped<pointSymmTensorField>(objects);
+        UnMapped<pointTensorField>(objects);
+    }
+
+    mapLagrangian(meshToMesh0Interp);
+}
+
+
+template<template<class> class CombineOp>
+void MapConsistentSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh0::order& mapOrder
+)
+{
+    HashTable<word> patchMap;
+    HashTable<label> cuttingPatchTable;
+
+    forAll(meshTarget.boundary(), patchi)
+    {
+        if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
+        {
+            patchMap.insert
+            (
+                meshTarget.boundary()[patchi].name(),
+                meshTarget.boundary()[patchi].name()
+            );
+        }
+        else
+        {
+            cuttingPatchTable.insert
+            (
+                meshTarget.boundaryMesh()[patchi].name(),
+                -1
+            );
+        }
+    }
+
+    MapSubMesh<CombineOp>
+    (
+        meshSource,
+        meshTarget,
+        patchMap,
+        cuttingPatchTable.toc(),
+        mapOrder
+    );
 }
 
 
diff --git a/applications/utilities/preProcessing/mapFields/MapVolFields.H b/applications/utilities/preProcessing/mapFields/MapVolFields.H
index 791ffb4c9..c72e70962 100644
--- a/applications/utilities/preProcessing/mapFields/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapVolFields.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,11 +23,11 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MapConsistentVolFields_H
-#define MapConsistentVolFields_H
+#ifndef MapVolFields_H
+#define MapVolFields_H
 
 #include "GeometricField.H"
-#include "meshToMesh.H"
+#include "meshToMesh0.H"
 #include "IOobjectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,55 +39,62 @@ template<class Type, class CombineOp>
 void MapVolFields
 (
     const IOobjectList& objects,
-    const HashSet<word>& selectedFields,
-    const meshToMesh& interp,
+    const meshToMesh0& meshToMesh0Interp,
+    const meshToMesh0::order& mapOrder,
     const CombineOp& cop
 )
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+    const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
+    const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
 
-    const fvMesh& meshSource = static_cast<const fvMesh&>(interp.srcRegion());
-    const fvMesh& meshTarget = static_cast<const fvMesh&>(interp.tgtRegion());
+    word fieldClassName
+    (
+        GeometricField<Type, fvPatchField, volMesh>::typeName
+    );
 
-    IOobjectList fields = objects.lookupClass(fieldType::typeName);
+    IOobjectList fields = objects.lookupClass(fieldClassName);
 
     forAllIter(IOobjectList, fields, fieldIter)
     {
-        const word& fieldName = fieldIter()->name();
+        IOobject fieldTargetIOobject
+        (
+            fieldIter()->name(),
+            meshTarget.time().timeName(),
+            meshTarget,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        );
 
-        if (selectedFields.empty() || selectedFields.found(fieldName))
+        if (fieldTargetIOobject.headerOk())
         {
-            Info<< "    interpolating " << fieldName << endl;
+            Info<< "    interpolating " << fieldIter()->name()
+                << endl;
 
-            const fieldType fieldSource(*fieldIter(), meshSource);
-
-            IOobject targetIO
+            // Read field fieldSource
+            GeometricField<Type, fvPatchField, volMesh> fieldSource
             (
-                fieldName,
-                meshTarget.time().timeName(),
-                meshTarget,
-                IOobject::MUST_READ
+                *fieldIter(),
+                meshSource
             );
 
-            if (targetIO.headerOk())
-            {
-                fieldType fieldTarget(targetIO, meshTarget);
+            // Read fieldTarget
+            GeometricField<Type, fvPatchField, volMesh> fieldTarget
+            (
+                fieldTargetIOobject,
+                meshTarget
+            );
 
-                interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
+            // Interpolate field
+            meshToMesh0Interp.interpolate
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
 
-                fieldTarget.write();
-            }
-            else
-            {
-                targetIO.readOpt() = IOobject::NO_READ;
-
-                tmp<fieldType>
-                    tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
-
-                fieldType fieldTarget(targetIO, tfieldTarget);
-
-                fieldTarget.write();
-            }
+            // Write field
+            fieldTarget.write();
         }
     }
 }
diff --git a/applications/utilities/preProcessing/mapFields/UnMapped.H b/applications/utilities/preProcessing/mapFields/UnMapped.H
index 70f1e4029..bcbefb2bc 100644
--- a/applications/utilities/preProcessing/mapFields/UnMapped.H
+++ b/applications/utilities/preProcessing/mapFields/UnMapped.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/utilities/preProcessing/mapFields/createTimes.H b/applications/utilities/preProcessing/mapFields/createTimes.H
index f5fd92b60..91e6a349a 100644
--- a/applications/utilities/preProcessing/mapFields/createTimes.H
+++ b/applications/utilities/preProcessing/mapFields/createTimes.H
@@ -1,11 +1,15 @@
     Info<< "\nCreate databases as time" << endl;
 
-    HashTable<string> srcOptions(args.options());
-    srcOptions.erase("case");
-    srcOptions.insert("case", fileName(rootDirSource/caseDirSource));
+    Time runTimeSource
+    (
+        Time::controlDictName,
+        rootDirSource,
+        caseDirSource
+    );
 
-    argList argsSrc(args, srcOptions, false, false, false);
-
-    Time runTimeSource(Time::controlDictName, argsSrc);
-
-    Time runTimeTarget(Time::controlDictName, args);
+    Time runTimeTarget
+    (
+        Time::controlDictName,
+        rootDirTarget,
+        caseDirTarget
+    );
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 2109c8877..929776d9d 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -27,12 +27,14 @@ Application
 Description
     Maps volume fields from one mesh to another, reading and
     interpolating all fields present in the time directory of both cases.
+    Parallel and non-parallel cases are handled without the need to reconstruct
+    them first.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "meshToMesh.H"
-#include "processorPolyPatch.H"
+#include "meshToMesh0.H"
+#include "processorFvPatch.H"
 #include "MapMeshes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -41,33 +43,26 @@ void mapConsistentMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::interpolationMethod& mapMethod,
-    const bool subtract,
-    const HashSet<word>& selectedFields,
-    const bool noLagrangian
+    const meshToMesh0::order& mapOrder,
+    const bool subtract
 )
 {
-    Info<< nl << "Consistently creating and mapping fields for time "
-        << meshSource.time().timeName() << nl << endl;
-
-    meshToMesh interp(meshSource, meshTarget, mapMethod);
-
     if (subtract)
     {
-        MapMesh<minusEqOp>
+        MapConsistentMesh<minusEqOp>
         (
-            interp,
-            selectedFields,
-            noLagrangian
+            meshSource,
+            meshTarget,
+            mapOrder
         );
     }
     else
     {
-        MapMesh<plusEqOp>
+        MapConsistentMesh<eqOp>
         (
-            interp,
-            selectedFields,
-            noLagrangian
+            meshSource,
+            meshTarget,
+            mapOrder
         );
     }
 }
@@ -79,40 +74,59 @@ void mapSubMesh
     const fvMesh& meshTarget,
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
-    const meshToMesh::interpolationMethod& mapMethod,
-    const bool subtract,
-    const HashSet<word>& selectedFields,
-    const bool noLagrangian
+    const meshToMesh0::order& mapOrder,
+    const bool subtract
 )
 {
-    Info<< nl << "Creating and mapping fields for time "
-        << meshSource.time().timeName() << nl << endl;
-
-    meshToMesh interp
-    (
-        meshSource,
-        meshTarget,
-        mapMethod,
-        patchMap,
-        cuttingPatches
-    );
-
     if (subtract)
     {
-        MapMesh<minusEqOp>
+        MapSubMesh<minusEqOp>
         (
-            interp,
-            selectedFields,
-            noLagrangian
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
         );
     }
     else
     {
-        MapMesh<plusEqOp>
+        MapSubMesh<eqOp>
         (
-            interp,
-            selectedFields,
-            noLagrangian
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
+        );
+    }
+}
+
+
+void mapConsistentSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh0::order& mapOrder,
+    const bool subtract
+)
+{
+    if (subtract)
+    {
+        MapConsistentSubMesh<minusEqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
+    }
+    else
+    {
+        MapConsistentSubMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
         );
     }
 }
@@ -125,20 +139,30 @@ wordList addProcessorPatches
 )
 {
     // Add the processor patches to the cutting list
-    HashSet<word> cuttingPatchTable;
+    HashTable<label> cuttingPatchTable;
     forAll(cuttingPatches, i)
     {
-        cuttingPatchTable.insert(cuttingPatches[i]);
+        cuttingPatchTable.insert(cuttingPatches[i], i);
     }
 
-    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();
-
-    forAll(pbm, patchI)
+    forAll(meshTarget.boundary(), patchi)
     {
-        if (isA<processorPolyPatch>(pbm[patchI]))
+        if (isA<processorFvPatch>(meshTarget.boundary()[patchi]))
         {
-            const word& patchName = pbm[patchI].name();
-            cuttingPatchTable.insert(patchName);
+            if
+            (
+               !cuttingPatchTable.found
+                (
+                    meshTarget.boundaryMesh()[patchi].name()
+                )
+            )
+            {
+                cuttingPatchTable.insert
+                (
+                    meshTarget.boundaryMesh()[patchi].name(),
+                    -1
+                );
+            }
         }
     }
 
@@ -154,7 +178,7 @@ int main(int argc, char *argv[])
     (
         "map volume fields from one mesh to another"
     );
-
+    argList::noParallel();
     argList::validArgs.append("sourceCase");
 
     argList::addOption
@@ -176,6 +200,16 @@ int main(int argc, char *argv[])
         "specify the target region"
     );
     argList::addBoolOption
+    (
+        "parallelSource",
+        "the source is decomposed"
+    );
+    argList::addBoolOption
+    (
+        "parallelTarget",
+        "the target is decomposed"
+    );
+    argList::addBoolOption
     (
         "consistent",
         "source and target geometry and boundary conditions identical"
@@ -191,21 +225,14 @@ int main(int argc, char *argv[])
         "subtract",
         "subtract mapped source from target"
     );
-    argList::addOption
-    (
-        "fields",
-        "list",
-        "specify a list of fields to be mapped. Eg, '(U T p)' - "
-        "regular expressions not currently supported"
-    );
-    argList::addBoolOption
-    (
-        "noLagrangian",
-        "skip mapping lagrangian positions and fields"
-    );
 
     argList args(argc, argv);
 
+    if (!args.check())
+    {
+        FatalError.exit();
+    }
+
     fileName rootDirTarget(args.rootPath());
     fileName caseDirTarget(args.globalCaseName());
 
@@ -229,17 +256,35 @@ int main(int argc, char *argv[])
         Info<< "Target region: " << targetRegion << endl;
     }
 
+    const bool parallelSource = args.optionFound("parallelSource");
+    const bool parallelTarget = args.optionFound("parallelTarget");
     const bool consistent = args.optionFound("consistent");
 
-    meshToMesh::interpolationMethod mapMethod =
-        meshToMesh::imCellVolumeWeight;
-
+    meshToMesh0::order mapOrder = meshToMesh0::INTERPOLATE;
     if (args.optionFound("mapMethod"))
     {
-        mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]];
+        const word mapMethod(args["mapMethod"]);
+        if (mapMethod == "mapNearest")
+        {
+            mapOrder = meshToMesh0::MAP;
+        }
+        else if (mapMethod == "interpolate")
+        {
+            mapOrder = meshToMesh0::INTERPOLATE;
+        }
+        else if (mapMethod == "cellPointInterpolate")
+        {
+            mapOrder = meshToMesh0::CELL_POINT_INTERPOLATE;
+        }
+        else
+        {
+            FatalErrorIn(args.executable())
+                << "Unknown mapMethod " << mapMethod << ". Valid options are: "
+                << "mapNearest, interpolate and cellPointInterpolate"
+                << exit(FatalError);
+        }
 
-        Info<< "Mapping method: "
-            << meshToMesh::interpolationMethodNames_[mapMethod] << endl;
+        Info<< "Mapping method: " << mapMethod << endl;
     }
 
     const bool subtract = args.optionFound("subtract");
@@ -248,13 +293,6 @@ int main(int argc, char *argv[])
         Info<< "Subtracting mapped source field from target" << endl;
     }
 
-    HashSet<word> selectedFields;
-    if (args.optionFound("fields"))
-    {
-        args.optionLookup("fields")() >> selectedFields;
-    }
-
-    const bool noLagrangian = args.optionFound("noLagrangian");
 
     #include "createTimes.H"
 
@@ -280,58 +318,341 @@ int main(int argc, char *argv[])
         mapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;
     }
 
-    #include "setTimeIndex.H"
-
-    Info<< "\nCreate meshes\n" << endl;
-
-    fvMesh meshSource
-    (
-        IOobject
-        (
-            sourceRegion,
-            runTimeSource.timeName(),
-            runTimeSource
-        )
-    );
-
-    fvMesh meshTarget
-    (
-        IOobject
-        (
-            targetRegion,
-            runTimeTarget.timeName(),
-            runTimeTarget
-        )
-    );
-
-    Info<< "Source mesh size: " << meshSource.nCells() << tab
-        << "Target mesh size: " << meshTarget.nCells() << nl << endl;
-
-    if (consistent)
+    if (parallelSource && !parallelTarget)
     {
-        mapConsistentMesh
+        IOdictionary decompositionDict
         (
-            meshSource,
-            meshTarget,
-            mapMethod,
-            subtract,
-            selectedFields,
-            noLagrangian
+            IOobject
+            (
+                "decomposeParDict",
+                runTimeSource.system(),
+                runTimeSource,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
         );
+
+        int nProcs(readInt(decompositionDict.lookup("numberOfSubdomains")));
+
+        Info<< "Create target mesh\n" << endl;
+
+        fvMesh meshTarget
+        (
+            IOobject
+            (
+                targetRegion,
+                runTimeTarget.timeName(),
+                runTimeTarget
+            )
+        );
+
+        Info<< "Target mesh size: " << meshTarget.nCells() << endl;
+
+        for (int procI=0; procI<nProcs; procI++)
+        {
+            Info<< nl << "Source processor " << procI << endl;
+
+            Time runTimeSource
+            (
+                Time::controlDictName,
+                rootDirSource,
+                caseDirSource/fileName(word("processor") + name(procI))
+            );
+
+            #include "setTimeIndex.H"
+
+            fvMesh meshSource
+            (
+                IOobject
+                (
+                    sourceRegion,
+                    runTimeSource.timeName(),
+                    runTimeSource
+                )
+            );
+
+            Info<< "mesh size: " << meshSource.nCells() << endl;
+
+            if (consistent)
+            {
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
+            }
+            else
+            {
+                mapSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    patchMap,
+                    cuttingPatches,
+                    mapOrder,
+                    subtract
+                );
+            }
+        }
+    }
+    else if (!parallelSource && parallelTarget)
+    {
+        IOdictionary decompositionDict
+        (
+            IOobject
+            (
+                "decomposeParDict",
+                runTimeTarget.system(),
+                runTimeTarget,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
+        );
+
+        int nProcs(readInt(decompositionDict.lookup("numberOfSubdomains")));
+
+        Info<< "Create source mesh\n" << endl;
+
+        #include "setTimeIndex.H"
+
+        fvMesh meshSource
+        (
+            IOobject
+            (
+                sourceRegion,
+                runTimeSource.timeName(),
+                runTimeSource
+            )
+        );
+
+        Info<< "Source mesh size: " << meshSource.nCells() << endl;
+
+        for (int procI=0; procI<nProcs; procI++)
+        {
+            Info<< nl << "Target processor " << procI << endl;
+
+            Time runTimeTarget
+            (
+                Time::controlDictName,
+                rootDirTarget,
+                caseDirTarget/fileName(word("processor") + name(procI))
+            );
+
+            fvMesh meshTarget
+            (
+                IOobject
+                (
+                    targetRegion,
+                    runTimeTarget.timeName(),
+                    runTimeTarget
+                )
+            );
+
+            Info<< "mesh size: " << meshTarget.nCells() << endl;
+
+            if (consistent)
+            {
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
+            }
+            else
+            {
+                mapSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    patchMap,
+                    addProcessorPatches(meshTarget, cuttingPatches),
+                    mapOrder,
+                    subtract
+                );
+            }
+        }
+    }
+    else if (parallelSource && parallelTarget)
+    {
+        IOdictionary decompositionDictSource
+        (
+            IOobject
+            (
+                "decomposeParDict",
+                runTimeSource.system(),
+                runTimeSource,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
+        );
+
+        int nProcsSource
+        (
+            readInt(decompositionDictSource.lookup("numberOfSubdomains"))
+        );
+
+
+        IOdictionary decompositionDictTarget
+        (
+            IOobject
+            (
+                "decomposeParDict",
+                runTimeTarget.system(),
+                runTimeTarget,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
+        );
+
+        int nProcsTarget
+        (
+            readInt(decompositionDictTarget.lookup("numberOfSubdomains"))
+        );
+
+        List<boundBox> bbsTarget(nProcsTarget);
+        List<bool> bbsTargetSet(nProcsTarget, false);
+
+        for (int procISource=0; procISource<nProcsSource; procISource++)
+        {
+            Info<< nl << "Source processor " << procISource << endl;
+
+            Time runTimeSource
+            (
+                Time::controlDictName,
+                rootDirSource,
+                caseDirSource/fileName(word("processor") + name(procISource))
+            );
+
+            #include "setTimeIndex.H"
+
+            fvMesh meshSource
+            (
+                IOobject
+                (
+                    sourceRegion,
+                    runTimeSource.timeName(),
+                    runTimeSource
+                )
+            );
+
+            Info<< "mesh size: " << meshSource.nCells() << endl;
+
+            boundBox bbSource(meshSource.bounds());
+
+            for (int procITarget=0; procITarget<nProcsTarget; procITarget++)
+            {
+                if
+                (
+                    !bbsTargetSet[procITarget]
+                  || (
+                      bbsTargetSet[procITarget]
+                   && bbsTarget[procITarget].overlaps(bbSource)
+                     )
+                )
+                {
+                    Info<< nl << "Target processor " << procITarget << endl;
+
+                    Time runTimeTarget
+                    (
+                        Time::controlDictName,
+                        rootDirTarget,
+                        caseDirTarget/fileName(word("processor")
+                      + name(procITarget))
+                    );
+
+                    fvMesh meshTarget
+                    (
+                        IOobject
+                        (
+                            targetRegion,
+                            runTimeTarget.timeName(),
+                            runTimeTarget
+                        )
+                    );
+
+                    Info<< "mesh size: " << meshTarget.nCells() << endl;
+
+                    bbsTarget[procITarget] = meshTarget.bounds();
+                    bbsTargetSet[procITarget] = true;
+
+                    if (bbsTarget[procITarget].overlaps(bbSource))
+                    {
+                        if (consistent)
+                        {
+                            mapConsistentSubMesh
+                            (
+                                meshSource,
+                                meshTarget,
+                                mapOrder,
+                                subtract
+                            );
+                        }
+                        else
+                        {
+                            mapSubMesh
+                            (
+                                meshSource,
+                                meshTarget,
+                                patchMap,
+                                addProcessorPatches(meshTarget, cuttingPatches),
+                                mapOrder,
+                                subtract
+                            );
+                        }
+                    }
+                }
+            }
+        }
     }
     else
     {
-        mapSubMesh
+        #include "setTimeIndex.H"
+
+        Info<< "Create meshes\n" << endl;
+
+        fvMesh meshSource
         (
-            meshSource,
-            meshTarget,
-            patchMap,
-            addProcessorPatches(meshTarget, cuttingPatches),
-            mapMethod,
-            subtract,
-            selectedFields,
-            noLagrangian
+            IOobject
+            (
+                sourceRegion,
+                runTimeSource.timeName(),
+                runTimeSource
+            )
         );
+
+        fvMesh meshTarget
+        (
+            IOobject
+            (
+                targetRegion,
+                runTimeTarget.timeName(),
+                runTimeTarget
+            )
+        );
+
+        Info<< "Source mesh size: " << meshSource.nCells() << tab
+            << "Target mesh size: " << meshTarget.nCells() << nl << endl;
+
+        if (consistent)
+        {
+            mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
+        }
+        else
+        {
+            mapSubMesh
+            (
+                meshSource,
+                meshTarget,
+                patchMap,
+                cuttingPatches,
+                mapOrder,
+                subtract
+            );
+        }
     }
 
     Info<< "\nEnd\n" << endl;
diff --git a/applications/utilities/preProcessing/mapFields/mapFieldsDict b/applications/utilities/preProcessing/mapFields/mapFieldsDict
index 65e33997d..925c20d86 100644
--- a/applications/utilities/preProcessing/mapFields/mapFieldsDict
+++ b/applications/utilities/preProcessing/mapFields/mapFieldsDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
+|  \\    /   O peration     | Version:  2.2.2                                 |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -14,12 +14,6 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Specify how to map patches. There are three different options:
-// - patch exists in the source case: specify mapping (patchMap)
-// - patch should be interpolated from internal values in source case
-//   (cuttingPatches)
-// - patch should not be mapped. Default if not in patchMap or cuttingPatches
-
 // List of pairs of target/source patches for mapping
 patchMap
 (
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
index d2ed7da59..24db86343 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
@@ -80,17 +80,27 @@ static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
 }
 
 
-void mapLagrangian(const meshToMesh& interp)
+void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
 {
     // Determine which particles are in meshTarget
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    const polyMesh& meshSource = interp.srcRegion();
-    const polyMesh& meshTarget = interp.tgtRegion();
-    const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
+    // target to source cell map
+    const labelList& cellAddressing = meshToMesh0Interp.cellAddressing();
 
+    // Invert celladdressing to get source to target(s).
+    // Note: could use sparse addressing but that is too storage inefficient
+    // (Map<labelList>)
+    labelListList sourceToTargets
+    (
+        invertOneToMany(meshToMesh0Interp.fromMesh().nCells(), cellAddressing)
+    );
+
+    const fvMesh& meshSource = meshToMesh0Interp.fromMesh();
+    const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
     const pointField& targetCc = meshTarget.cellCentres();
 
+
     fileNameList cloudDirs
     (
         readDir
@@ -110,7 +120,7 @@ void mapLagrangian(const meshToMesh& interp)
             cloud::prefix/cloudDirs[cloudI]
         );
 
-        IOobject* positionsPtr = objects.lookup(word("positions"));
+        IOobject* positionsPtr = objects.lookup("positions");
 
         if (positionsPtr)
         {
@@ -158,7 +168,7 @@ void mapLagrangian(const meshToMesh& interp)
                 if (iter().cell() >= 0)
                 {
                     const labelList& targetCells =
-                        sourceToTarget[iter().cell()];
+                        sourceToTargets[iter().cell()];
 
                     // Particle probably in one of the targetcells. Try
                     // all by tracking from their cell centre to the parcel
@@ -200,7 +210,7 @@ void mapLagrangian(const meshToMesh& interp)
                 sourceParticleI++;
             }
 
-            Info<< "    after meshToMesh addressing found "
+            Info<< "    after meshToMesh0 addressing found "
                 << targetParcels.size()
                 << " parcels in target mesh." << endl;
 
@@ -249,47 +259,17 @@ void mapLagrangian(const meshToMesh& interp)
                 // ~~~~~~~~~~~~~~~~~~~~~
 
                 MapLagrangianFields<label>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
                 MapLagrangianFields<scalar>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
                 MapLagrangianFields<vector>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
                 MapLagrangianFields<sphericalTensor>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
                 MapLagrangianFields<symmTensor>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
                 MapLagrangianFields<tensor>
-                (
-                    cloudDirs[cloudI],
-                    objects,
-                    meshTarget,
-                    addParticles
-                );
+                (cloudDirs[cloudI], objects, meshToMesh0Interp, addParticles);
             }
         }
     }
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.H b/applications/utilities/preProcessing/mapFields/mapLagrangian.H
index 82e1c589c..5c4b4e4b8 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.H
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef mapLagrangian_H
 #define mapLagrangian_H
 
-#include "meshToMesh.H"
+#include "meshToMesh0.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,7 +43,7 @@ namespace Foam
 {
 
 //- Maps lagrangian positions and fields
-void mapLagrangian(const meshToMesh& interp);
+void mapLagrangian(const meshToMesh0& meshToMesh0Interp);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/utilities/preProcessing/mapFields/setTimeIndex.H b/applications/utilities/preProcessing/mapFields/setTimeIndex.H
index a468d5064..20012e3d2 100644
--- a/applications/utilities/preProcessing/mapFields/setTimeIndex.H
+++ b/applications/utilities/preProcessing/mapFields/setTimeIndex.H
@@ -1,4 +1,3 @@
-{
     instantList sourceTimes = runTimeSource.times();
     label sourceTimeIndex = runTimeSource.timeIndex();
     if (args.optionFound("sourceTime"))
@@ -30,4 +29,3 @@
     Info<< "\nSource time: " << runTimeSource.value()
         << "\nTarget time: " << runTimeTarget.value()
         << endl;
-}
diff --git a/applications/utilities/preProcessing/mapFieldsPar/Make/files b/applications/utilities/preProcessing/mapFieldsPar/Make/files
new file mode 100644
index 000000000..63742e0e3
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/Make/files
@@ -0,0 +1,4 @@
+mapLagrangian.C
+mapFieldsPar.C
+
+EXE = $(FOAM_APPBIN)/mapFieldsPar
diff --git a/applications/utilities/preProcessing/mapFieldsPar/Make/options b/applications/utilities/preProcessing/mapFieldsPar/Make/options
new file mode 100644
index 000000000..7bd964094
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/Make/options
@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lsampling \
+    -lmeshTools \
+    -llagrangian \
+    -lfiniteVolume \
+    -lgenericPatchFields
diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H b/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H
new file mode 100644
index 000000000..e04552214
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/MapLagrangianFields.H
@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+InNamespace
+    Foam
+
+Description
+    Gets the indices of (source)particles that have been appended to the
+    target cloud and maps the lagrangian fields accordingly.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MapLagrangianFields_H
+#define MapLagrangianFields_H
+
+#include "cloud.H"
+#include "GeometricField.H"
+#include "meshToMesh.H"
+#include "IOobjectList.H"
+#include "CompactIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Gets the indices of (source)particles that have been appended to the
+//  target cloud and maps the lagrangian fields accordingly.
+template<class Type>
+void MapLagrangianFields
+(
+    const string& cloudName,
+    const IOobjectList& objects,
+    const polyMesh& meshTarget,
+    const labelList& addParticles
+)
+{
+    {
+        IOobjectList fields = objects.lookupClass(IOField<Type>::typeName);
+
+        forAllIter(IOobjectList, fields, fieldIter)
+        {
+            const word& fieldName = fieldIter()->name();
+
+            Info<< "    mapping lagrangian field " << fieldName << endl;
+
+            // Read field (does not need mesh)
+            IOField<Type> fieldSource(*fieldIter());
+
+            // Map
+            IOField<Type> fieldTarget
+            (
+                IOobject
+                (
+                    fieldName,
+                    meshTarget.time().timeName(),
+                    cloud::prefix/cloudName,
+                    meshTarget,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                addParticles.size()
+            );
+
+            forAll(addParticles, i)
+            {
+                fieldTarget[i] = fieldSource[addParticles[i]];
+            }
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+
+    {
+        IOobjectList fieldFields =
+            objects.lookupClass(IOField<Field<Type> >::typeName);
+
+        forAllIter(IOobjectList, fieldFields, fieldIter)
+        {
+            const word& fieldName = fieldIter()->name();
+
+            Info<< "    mapping lagrangian fieldField " << fieldName << endl;
+
+            // Read field (does not need mesh)
+            IOField<Field<Type> > fieldSource(*fieldIter());
+
+            // Map - use CompactIOField to automatically write in
+            // compact form for binary format.
+            CompactIOField<Field<Type>, Type> fieldTarget
+            (
+                IOobject
+                (
+                    fieldName,
+                    meshTarget.time().timeName(),
+                    cloud::prefix/cloudName,
+                    meshTarget,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                addParticles.size()
+            );
+
+            forAll(addParticles, i)
+            {
+                fieldTarget[i] = fieldSource[addParticles[i]];
+            }
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+
+    {
+        IOobjectList fieldFields =
+            objects.lookupClass(CompactIOField<Field<Type>, Type>::typeName);
+
+        forAllIter(IOobjectList, fieldFields, fieldIter)
+        {
+            Info<< "    mapping lagrangian fieldField "
+                << fieldIter()->name() << endl;
+
+            // Read field (does not need mesh)
+            CompactIOField<Field<Type>, Type> fieldSource(*fieldIter());
+
+            // Map
+            CompactIOField<Field<Type>, Type> fieldTarget
+            (
+                IOobject
+                (
+                    fieldIter()->name(),
+                    meshTarget.time().timeName(),
+                    cloud::prefix/cloudName,
+                    meshTarget,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                addParticles.size()
+            );
+
+            forAll(addParticles, i)
+            {
+                fieldTarget[i] = fieldSource[addParticles[i]];
+            }
+
+            // Write field
+            fieldTarget.write();
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H b/applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H
new file mode 100644
index 000000000..83e378992
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MapMeshes_H
+#define MapMeshes_H
+
+#include "MapVolFields.H"
+#include "mapLagrangian.H"
+#include "UnMapped.H"
+#include "pointMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<template<class> class CombineOp>
+void MapMesh
+(
+    const meshToMesh& interp,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+)
+{
+    {
+        const polyMesh& meshSource = interp.srcRegion();
+
+        // Search for list of objects for this time
+        IOobjectList objects(meshSource, meshSource.time().timeName());
+
+        // Map volFields
+        // ~~~~~~~~~~~~~
+        MapVolFields<scalar>
+        (
+            objects,
+            selectedFields,
+            interp,
+            CombineOp<scalar>()
+        );
+
+        MapVolFields<vector>
+        (
+            objects,
+            selectedFields,
+            interp,
+            CombineOp<vector>()
+        );
+        MapVolFields<sphericalTensor>
+        (
+            objects,
+            selectedFields,
+            interp,
+            CombineOp<sphericalTensor>()
+        );
+        MapVolFields<symmTensor>
+        (
+            objects,
+            selectedFields,
+            interp,
+            CombineOp<symmTensor>()
+        );
+        MapVolFields<tensor>
+        (
+            objects,
+            selectedFields,
+            interp,
+            CombineOp<tensor>()
+        );
+    }
+
+    {
+        const polyMesh& meshTarget = interp.tgtRegion();
+
+        // Search for list of target objects for this time
+        IOobjectList objects(meshTarget, meshTarget.time().timeName());
+
+        // Mark surfaceFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<surfaceScalarField>(objects);
+        UnMapped<surfaceVectorField>(objects);
+        UnMapped<surfaceSphericalTensorField>(objects);
+        UnMapped<surfaceSymmTensorField>(objects);
+        UnMapped<surfaceTensorField>(objects);
+
+        // Mark pointFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<pointScalarField>(objects);
+        UnMapped<pointVectorField>(objects);
+        UnMapped<pointSphericalTensorField>(objects);
+        UnMapped<pointSymmTensorField>(objects);
+        UnMapped<pointTensorField>(objects);
+    }
+
+    if (!noLagrangian)
+    {
+        mapLagrangian(interp);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H b/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
new file mode 100644
index 000000000..ac5c18b07
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MapConsistentVolFields_H
+#define MapConsistentVolFields_H
+
+#include "GeometricField.H"
+#include "meshToMesh.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type, class CombineOp>
+void MapVolFields
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields,
+    const meshToMesh& interp,
+    const CombineOp& cop
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const fvMesh& meshSource = static_cast<const fvMesh&>(interp.srcRegion());
+    const fvMesh& meshTarget = static_cast<const fvMesh&>(interp.tgtRegion());
+
+    IOobjectList fields = objects.lookupClass(fieldType::typeName);
+
+    forAllIter(IOobjectList, fields, fieldIter)
+    {
+        const word& fieldName = fieldIter()->name();
+
+        if (selectedFields.empty() || selectedFields.found(fieldName))
+        {
+            Info<< "    interpolating " << fieldName << endl;
+
+            const fieldType fieldSource(*fieldIter(), meshSource);
+
+            IOobject targetIO
+            (
+                fieldName,
+                meshTarget.time().timeName(),
+                meshTarget,
+                IOobject::MUST_READ
+            );
+
+            if (targetIO.headerOk())
+            {
+                fieldType fieldTarget(targetIO, meshTarget);
+
+                interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
+
+                fieldTarget.write();
+            }
+            else
+            {
+                targetIO.readOpt() = IOobject::NO_READ;
+
+                tmp<fieldType>
+                    tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
+
+                fieldType fieldTarget(targetIO, tfieldTarget);
+
+                fieldTarget.write();
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/UnMapped.H b/applications/utilities/preProcessing/mapFieldsPar/UnMapped.H
new file mode 100644
index 000000000..bcbefb2bc
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/UnMapped.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef UnMapped_H
+#define UnMapped_H
+
+#include "IOobjectList.H"
+#include "OSspecific.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type>
+void UnMapped(const IOobjectList& objects)
+{
+    IOobjectList fields = objects.lookupClass(Type::typeName);
+
+    forAllConstIter(IOobjectList, fields, fieldIter)
+    {
+        mvBak(fieldIter()->objectPath(), "unmapped");
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/createTimes.H b/applications/utilities/preProcessing/mapFieldsPar/createTimes.H
new file mode 100644
index 000000000..f5fd92b60
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/createTimes.H
@@ -0,0 +1,11 @@
+    Info<< "\nCreate databases as time" << endl;
+
+    HashTable<string> srcOptions(args.options());
+    srcOptions.erase("case");
+    srcOptions.insert("case", fileName(rootDirSource/caseDirSource));
+
+    argList argsSrc(args, srcOptions, false, false, false);
+
+    Time runTimeSource(Time::controlDictName, argsSrc);
+
+    Time runTimeTarget(Time::controlDictName, args);
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsDict b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsDict
new file mode 100644
index 000000000..470250179
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsDict
@@ -0,0 +1,36 @@
+/*--------------------------------*- 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      mapFieldsParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Specify how to map patches. There are three different options:
+// - patch exists in the source case: specify mapping (patchMap)
+// - patch should be interpolated from internal values in source case
+//   (cuttingPatches)
+// - patch should not be mapped. Default if not in patchMap or cuttingPatches
+
+// List of pairs of target/source patches for mapping
+patchMap
+(
+    lid movingWall
+);
+
+// List of target patches cutting the source domain (these need to be
+// handled specially e.g. interpolated from internal values)
+cuttingPatches
+(
+    fixedWalls
+);
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
new file mode 100644
index 000000000..1ef16c7d2
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
@@ -0,0 +1,343 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+Application
+    mapFieldsPar
+
+Description
+    Maps volume fields from one mesh to another, reading and
+    interpolating all fields present in the time directory of both cases.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "meshToMesh.H"
+#include "processorPolyPatch.H"
+#include "MapMeshes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void mapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh::interpolationMethod& mapMethod,
+    const bool subtract,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+)
+{
+    Info<< nl << "Consistently creating and mapping fields for time "
+        << meshSource.time().timeName() << nl << endl;
+
+    meshToMesh interp(meshSource, meshTarget, mapMethod);
+
+    if (subtract)
+    {
+        MapMesh<minusEqOp>
+        (
+            interp,
+            selectedFields,
+            noLagrangian
+        );
+    }
+    else
+    {
+        MapMesh<plusEqOp>
+        (
+            interp,
+            selectedFields,
+            noLagrangian
+        );
+    }
+}
+
+
+void mapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const meshToMesh::interpolationMethod& mapMethod,
+    const bool subtract,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+)
+{
+    Info<< nl << "Creating and mapping fields for time "
+        << meshSource.time().timeName() << nl << endl;
+
+    meshToMesh interp
+    (
+        meshSource,
+        meshTarget,
+        mapMethod,
+        patchMap,
+        cuttingPatches
+    );
+
+    if (subtract)
+    {
+        MapMesh<minusEqOp>
+        (
+            interp,
+            selectedFields,
+            noLagrangian
+        );
+    }
+    else
+    {
+        MapMesh<plusEqOp>
+        (
+            interp,
+            selectedFields,
+            noLagrangian
+        );
+    }
+}
+
+
+wordList addProcessorPatches
+(
+    const fvMesh& meshTarget,
+    const wordList& cuttingPatches
+)
+{
+    // Add the processor patches to the cutting list
+    HashSet<word> cuttingPatchTable;
+    forAll(cuttingPatches, i)
+    {
+        cuttingPatchTable.insert(cuttingPatches[i]);
+    }
+
+    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();
+
+    forAll(pbm, patchI)
+    {
+        if (isA<processorPolyPatch>(pbm[patchI]))
+        {
+            const word& patchName = pbm[patchI].name();
+            cuttingPatchTable.insert(patchName);
+        }
+    }
+
+    return cuttingPatchTable.toc();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "map volume fields from one mesh to another"
+    );
+
+    argList::validArgs.append("sourceCase");
+
+    argList::addOption
+    (
+        "sourceTime",
+        "scalar|'latestTime'",
+        "specify the source time"
+    );
+    argList::addOption
+    (
+        "sourceRegion",
+        "word",
+        "specify the source region"
+    );
+    argList::addOption
+    (
+        "targetRegion",
+        "word",
+        "specify the target region"
+    );
+    argList::addBoolOption
+    (
+        "consistent",
+        "source and target geometry and boundary conditions identical"
+    );
+    argList::addOption
+    (
+        "mapMethod",
+        "word",
+        "specify the mapping method"
+    );
+    argList::addBoolOption
+    (
+        "subtract",
+        "subtract mapped source from target"
+    );
+    argList::addOption
+    (
+        "fields",
+        "list",
+        "specify a list of fields to be mapped. Eg, '(U T p)' - "
+        "regular expressions not currently supported"
+    );
+    argList::addBoolOption
+    (
+        "noLagrangian",
+        "skip mapping lagrangian positions and fields"
+    );
+
+    argList args(argc, argv);
+
+    fileName rootDirTarget(args.rootPath());
+    fileName caseDirTarget(args.globalCaseName());
+
+    const fileName casePath = args[1];
+    const fileName rootDirSource = casePath.path();
+    const fileName caseDirSource = casePath.name();
+
+    Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;
+    word sourceRegion = fvMesh::defaultRegion;
+    if (args.optionFound("sourceRegion"))
+    {
+        sourceRegion = args["sourceRegion"];
+        Info<< "Source region: " << sourceRegion << endl;
+    }
+
+    Info<< "Target: " << rootDirTarget << " " << caseDirTarget << endl;
+    word targetRegion = fvMesh::defaultRegion;
+    if (args.optionFound("targetRegion"))
+    {
+        targetRegion = args["targetRegion"];
+        Info<< "Target region: " << targetRegion << endl;
+    }
+
+    const bool consistent = args.optionFound("consistent");
+
+    meshToMesh::interpolationMethod mapMethod =
+        meshToMesh::imCellVolumeWeight;
+
+    if (args.optionFound("mapMethod"))
+    {
+        mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]];
+
+        Info<< "Mapping method: "
+            << meshToMesh::interpolationMethodNames_[mapMethod] << endl;
+    }
+
+    const bool subtract = args.optionFound("subtract");
+    if (subtract)
+    {
+        Info<< "Subtracting mapped source field from target" << endl;
+    }
+
+    HashSet<word> selectedFields;
+    if (args.optionFound("fields"))
+    {
+        args.optionLookup("fields")() >> selectedFields;
+    }
+
+    const bool noLagrangian = args.optionFound("noLagrangian");
+
+    #include "createTimes.H"
+
+    HashTable<word> patchMap;
+    wordList cuttingPatches;
+
+    if (!consistent)
+    {
+        IOdictionary mapFieldsParDict
+        (
+            IOobject
+            (
+                "mapFieldsParDict",
+                runTimeTarget.system(),
+                runTimeTarget,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        mapFieldsParDict.lookup("patchMap") >> patchMap;
+        mapFieldsParDict.lookup("cuttingPatches") >>  cuttingPatches;
+    }
+
+    #include "setTimeIndex.H"
+
+    Info<< "\nCreate meshes\n" << endl;
+
+    fvMesh meshSource
+    (
+        IOobject
+        (
+            sourceRegion,
+            runTimeSource.timeName(),
+            runTimeSource
+        )
+    );
+
+    fvMesh meshTarget
+    (
+        IOobject
+        (
+            targetRegion,
+            runTimeTarget.timeName(),
+            runTimeTarget
+        )
+    );
+
+    Info<< "Source mesh size: " << meshSource.nCells() << tab
+        << "Target mesh size: " << meshTarget.nCells() << nl << endl;
+
+    if (consistent)
+    {
+        mapConsistentMesh
+        (
+            meshSource,
+            meshTarget,
+            mapMethod,
+            subtract,
+            selectedFields,
+            noLagrangian
+        );
+    }
+    else
+    {
+        mapSubMesh
+        (
+            meshSource,
+            meshTarget,
+            patchMap,
+            addProcessorPatches(meshTarget, cuttingPatches),
+            mapMethod,
+            subtract,
+            selectedFields,
+            noLagrangian
+        );
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
new file mode 100644
index 000000000..d2ed7da59
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
@@ -0,0 +1,303 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "MapLagrangianFields.H"
+#include "passiveParticleCloud.H"
+#include "meshSearch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+static const scalar perturbFactor = 1e-6;
+
+
+// Special version of findCell that generates a cell guaranteed to be
+// compatible with tracking.
+static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
+{
+    label cellI = -1;
+    label tetFaceI = -1;
+    label tetPtI = -1;
+
+    const polyMesh& mesh = cloud.pMesh();
+
+    mesh.findCellFacePt(pt, cellI, tetFaceI, tetPtI);
+
+    if (cellI >= 0)
+    {
+        return cellI;
+    }
+    else
+    {
+        // See if particle on face by finding nearest face and shifting
+        // particle.
+
+        meshSearch meshSearcher
+        (
+            mesh,
+            polyMesh::FACE_PLANES    // no decomposition needed
+        );
+
+        label faceI = meshSearcher.findNearestBoundaryFace(pt);
+
+        if (faceI >= 0)
+        {
+            const point& cc = mesh.cellCentres()[mesh.faceOwner()[faceI]];
+
+            const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc;
+
+            mesh.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI);
+
+            return cellI;
+        }
+    }
+
+    return -1;
+}
+
+
+void mapLagrangian(const meshToMesh& interp)
+{
+    // Determine which particles are in meshTarget
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const polyMesh& meshSource = interp.srcRegion();
+    const polyMesh& meshTarget = interp.tgtRegion();
+    const labelListList& sourceToTarget = interp.srcToTgtCellAddr();
+
+    const pointField& targetCc = meshTarget.cellCentres();
+
+    fileNameList cloudDirs
+    (
+        readDir
+        (
+            meshSource.time().timePath()/cloud::prefix,
+            fileName::DIRECTORY
+        )
+    );
+
+    forAll(cloudDirs, cloudI)
+    {
+        // Search for list of lagrangian objects for this time
+        IOobjectList objects
+        (
+            meshSource,
+            meshSource.time().timeName(),
+            cloud::prefix/cloudDirs[cloudI]
+        );
+
+        IOobject* positionsPtr = objects.lookup(word("positions"));
+
+        if (positionsPtr)
+        {
+            Info<< nl << "    processing cloud " << cloudDirs[cloudI] << endl;
+
+            // Read positions & cell
+            passiveParticleCloud sourceParcels
+            (
+                meshSource,
+                cloudDirs[cloudI],
+                false
+            );
+            Info<< "    read " << sourceParcels.size()
+                << " parcels from source mesh." << endl;
+
+            // Construct empty target cloud
+            passiveParticleCloud targetParcels
+            (
+                meshTarget,
+                cloudDirs[cloudI],
+                IDLList<passiveParticle>()
+            );
+
+            particle::TrackingData<passiveParticleCloud> td(targetParcels);
+
+            label sourceParticleI = 0;
+
+            // Indices of source particles that get added to targetParcels
+            DynamicList<label> addParticles(sourceParcels.size());
+
+            // Unmapped particles
+            labelHashSet unmappedSource(sourceParcels.size());
+
+
+            // Initial: track from fine-mesh cell centre to particle position
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            // This requires there to be no boundary in the way.
+
+
+            forAllConstIter(Cloud<passiveParticle>, sourceParcels, iter)
+            {
+                bool foundCell = false;
+
+                // Assume that cell from read parcel is the correct one...
+                if (iter().cell() >= 0)
+                {
+                    const labelList& targetCells =
+                        sourceToTarget[iter().cell()];
+
+                    // Particle probably in one of the targetcells. Try
+                    // all by tracking from their cell centre to the parcel
+                    // position.
+
+                    forAll(targetCells, i)
+                    {
+                        // Track from its cellcentre to position to make sure.
+                        autoPtr<passiveParticle> newPtr
+                        (
+                            new passiveParticle
+                            (
+                                meshTarget,
+                                targetCc[targetCells[i]],
+                                targetCells[i]
+                            )
+                        );
+                        passiveParticle& newP = newPtr();
+
+                        label faceI = newP.track(iter().position(), td);
+
+                        if (faceI < 0 && newP.cell() >= 0)
+                        {
+                            // Hit position.
+                            foundCell = true;
+                            addParticles.append(sourceParticleI);
+                            targetParcels.addParticle(newPtr.ptr());
+                            break;
+                        }
+                    }
+                }
+
+                if (!foundCell)
+                {
+                    // Store for closer analysis
+                    unmappedSource.insert(sourceParticleI);
+                }
+
+                sourceParticleI++;
+            }
+
+            Info<< "    after meshToMesh addressing found "
+                << targetParcels.size()
+                << " parcels in target mesh." << endl;
+
+
+            // Do closer inspection for unmapped particles
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            if (unmappedSource.size())
+            {
+                sourceParticleI = 0;
+
+                forAllIter(Cloud<passiveParticle>, sourceParcels, iter)
+                {
+                    if (unmappedSource.found(sourceParticleI))
+                    {
+                        label targetCell =
+                            findCell(targetParcels, iter().position());
+
+                        if (targetCell >= 0)
+                        {
+                            unmappedSource.erase(sourceParticleI);
+                            addParticles.append(sourceParticleI);
+                            iter().cell() = targetCell;
+                            targetParcels.addParticle
+                            (
+                                sourceParcels.remove(&iter())
+                            );
+                        }
+                    }
+                    sourceParticleI++;
+                }
+            }
+            addParticles.shrink();
+
+            Info<< "    after additional mesh searching found "
+                << targetParcels.size() << " parcels in target mesh." << endl;
+
+            if (addParticles.size())
+            {
+                IOPosition<passiveParticleCloud>(targetParcels).write();
+
+                // addParticles now contains the indices of the sourceMesh
+                // particles that were appended to the target mesh.
+
+                // Map lagrangian fields
+                // ~~~~~~~~~~~~~~~~~~~~~
+
+                MapLagrangianFields<label>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+                MapLagrangianFields<scalar>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+                MapLagrangianFields<vector>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+                MapLagrangianFields<sphericalTensor>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+                MapLagrangianFields<symmTensor>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+                MapLagrangianFields<tensor>
+                (
+                    cloudDirs[cloudI],
+                    objects,
+                    meshTarget,
+                    addParticles
+                );
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.H b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.H
new file mode 100644
index 000000000..dc9cf1560
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.H
@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+InNamespace
+    Foam
+
+Description
+    Maps lagrangian positions and fields
+
+SourceFiles
+    mapLagrangian.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mapLagrangian_H
+#define mapLagrangian_H
+
+#include "meshToMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Maps lagrangian positions and fields
+void mapLagrangian(const meshToMesh& interp);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/setTimeIndex.H b/applications/utilities/preProcessing/mapFieldsPar/setTimeIndex.H
new file mode 100644
index 000000000..a468d5064
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/setTimeIndex.H
@@ -0,0 +1,33 @@
+{
+    instantList sourceTimes = runTimeSource.times();
+    label sourceTimeIndex = runTimeSource.timeIndex();
+    if (args.optionFound("sourceTime"))
+    {
+        if (args["sourceTime"] == "latestTime")
+        {
+            sourceTimeIndex = sourceTimes.size() - 1;
+        }
+        else
+        {
+            sourceTimeIndex = Time::findClosestTimeIndex
+            (
+                sourceTimes,
+                args.optionRead<scalar>("sourceTime")
+            );
+        }
+    }
+    else
+    {
+        sourceTimeIndex = Time::findClosestTimeIndex
+        (
+            sourceTimes,
+            runTimeTarget.time().value()
+        );
+    }
+
+    runTimeSource.setTime(sourceTimes[sourceTimeIndex], sourceTimeIndex);
+
+    Info<< "\nSource time: " << runTimeSource.value()
+        << "\nTarget time: " << runTimeTarget.value()
+        << endl;
+}
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 135032b10..8bbc898dd 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -55,15 +55,18 @@ graphField/writePatchGraph.C
 graphField/writeCellGraph.C
 graphField/makeGraph.C
 
-meshToMesh = meshToMeshInterpolation/meshToMesh
-$(meshToMesh)/meshToMesh.C
-$(meshToMesh)/meshToMeshParallelOps.C
-meshToMeshMethods = meshToMeshInterpolation/meshToMesh/calcMethod
+meshToMesh/meshToMesh.C
+meshToMesh/meshToMeshParallelOps.C
+meshToMeshMethods = meshToMesh/calcMethod
 $(meshToMeshMethods)/meshToMeshMethod/meshToMeshMethod.C
 $(meshToMeshMethods)/meshToMeshMethod/meshToMeshMethodNew.C
 $(meshToMeshMethods)/cellVolumeWeight/cellVolumeWeightMethod.C
 $(meshToMeshMethods)/direct/directMethod.C
 $(meshToMeshMethods)/mapNearest/mapNearestMethod.C
 
+meshToMesh0/meshToMesh0.C
+meshToMesh0/calculateMeshToMesh0Addressing.C
+meshToMesh0/calculateMeshToMesh0Weights.C
+
 
 LIB = $(FOAM_LIBBIN)/libsampling
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C b/src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
rename to src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
index b1e5141dc..863301bbc 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H b/src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
similarity index 98%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
rename to src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
index 2a92f04ca..9e8171572 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
+++ b/src/sampling/meshToMesh/calcMethod/cellVolumeWeight/cellVolumeWeightMethod.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C b/src/sampling/meshToMesh/calcMethod/direct/directMethod.C
similarity index 100%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C
rename to src/sampling/meshToMesh/calcMethod/direct/directMethod.C
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.H b/src/sampling/meshToMesh/calcMethod/direct/directMethod.H
similarity index 98%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.H
rename to src/sampling/meshToMesh/calcMethod/direct/directMethod.H
index 414f3d7b5..ddbb2c208 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.H
+++ b/src/sampling/meshToMesh/calcMethod/direct/directMethod.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C b/src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C
rename to src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C
index 2ba0d76dd..27f1610d5 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H b/src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H
similarity index 98%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H
rename to src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H
index 6df1afee0..b26bc645a 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H
+++ b/src/sampling/meshToMesh/calcMethod/mapNearest/mapNearestMethod.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
rename to src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
index a6af83e72..11b0a5d59 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H
similarity index 98%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H
rename to src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H
index e76ae9176..7397673d5 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
similarity index 95%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
rename to src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
index 373e7aa86..1e1e0a3a8 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
similarity index 97%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
rename to src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
index 9554ba9f4..6b3abcf2e 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethodNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C b/src/sampling/meshToMesh/meshToMesh.C
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
rename to src/sampling/meshToMesh/meshToMesh.C
index 91c998ce4..c03c3789d 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMesh/meshToMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMesh/meshToMesh.H
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
rename to src/sampling/meshToMesh/meshToMesh.H
index 8d5961191..718270e4f 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMesh/meshToMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H b/src/sampling/meshToMesh/meshToMeshI.H
similarity index 97%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H
rename to src/sampling/meshToMesh/meshToMeshI.H
index 103b14a0c..a3d626754 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H
+++ b/src/sampling/meshToMesh/meshToMeshI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshParallelOps.C b/src/sampling/meshToMesh/meshToMeshParallelOps.C
similarity index 99%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshParallelOps.C
rename to src/sampling/meshToMesh/meshToMeshParallelOps.C
index 01f26bcc0..ce3fcc428 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshParallelOps.C
+++ b/src/sampling/meshToMesh/meshToMeshParallelOps.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshTemplates.C b/src/sampling/meshToMesh/meshToMeshTemplates.C
similarity index 100%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshTemplates.C
rename to src/sampling/meshToMesh/meshToMeshTemplates.C
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/weightedFvPatchFieldMapper.H b/src/sampling/meshToMesh/weightedFvPatchFieldMapper.H
similarity index 100%
rename from src/sampling/meshToMeshInterpolation/meshToMesh/weightedFvPatchFieldMapper.H
rename to src/sampling/meshToMesh/weightedFvPatchFieldMapper.H
diff --git a/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
new file mode 100644
index 000000000..35b8f63f2
--- /dev/null
+++ b/src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
@@ -0,0 +1,343 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+Description
+    private member of meshToMesh0.
+    Calculates mesh to mesh addressing pattern (for each cell from one mesh
+    find the closest cell centre in the other mesh).
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshToMesh0.H"
+#include "SubField.H"
+
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::meshToMesh0::calcAddressing()
+{
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateAddressing() : "
+            << "calculating mesh-to-mesh cell addressing" << endl;
+    }
+
+    // set reference to cells
+    const cellList& fromCells = fromMesh_.cells();
+    const pointField& fromPoints = fromMesh_.points();
+
+    // In an attempt to preserve the efficiency of linear search and prevent
+    // failure, a RESCUE mechanism will be set up. Here, we shall mark all
+    // cells next to the solid boundaries. If such a cell is found as the
+    // closest, the relationship between the origin and cell will be examined.
+    // If the origin is outside the cell, a global n-squared search is
+    // triggered.
+
+    // SETTING UP RESCUE
+
+    // visit all boundaries and mark the cell next to the boundary.
+
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateAddressing() : "
+            << "Setting up rescue" << endl;
+    }
+
+    List<bool> boundaryCell(fromCells.size(), false);
+
+    // set reference to boundary
+    const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
+
+    forAll(patchesFrom, patchI)
+    {
+        // get reference to cells next to the boundary
+        const labelUList& bCells = patchesFrom[patchI].faceCells();
+
+        forAll(bCells, faceI)
+        {
+            boundaryCell[bCells[faceI]] = true;
+        }
+    }
+
+    treeBoundBox meshBb(fromPoints);
+
+    scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
+
+    treeBoundBox shiftedBb
+    (
+        meshBb.min(),
+        meshBb.max() + vector(typDim, typDim, typDim)
+    );
+
+    if (debug)
+    {
+        Info<< "\nMesh" << endl;
+        Info<< "   bounding box           : " << meshBb << endl;
+        Info<< "   bounding box (shifted) : " << shiftedBb << endl;
+        Info<< "   typical dimension      :" << shiftedBb.typDim() << endl;
+    }
+
+    indexedOctree<treeDataCell> oc
+    (
+        treeDataCell(false, fromMesh_, polyMesh::FACE_DIAG_TRIS),
+        shiftedBb,      // overall bounding box
+        8,              // maxLevel
+        10,             // leafsize
+        6.0             // duplicity
+    );
+
+    if (debug)
+    {
+        oc.print(Pout, false, 0);
+    }
+
+    cellAddresses
+    (
+        cellAddressing_,
+        toMesh_.cellCentres(),
+        fromMesh_,
+        boundaryCell,
+        oc
+    );
+
+    forAll(toMesh_.boundaryMesh(), patchi)
+    {
+        const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
+
+        if (cuttingPatches_.found(toPatch.name()))
+        {
+            boundaryAddressing_[patchi].setSize(toPatch.size());
+
+            cellAddresses
+            (
+                boundaryAddressing_[patchi],
+                toPatch.faceCentres(),
+                fromMesh_,
+                boundaryCell,
+                oc
+            );
+        }
+        else if
+        (
+            patchMap_.found(toPatch.name())
+         && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
+        )
+        {
+            const polyPatch& fromPatch = fromMesh_.boundaryMesh()
+            [
+                fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
+            ];
+
+            if (fromPatch.empty())
+            {
+                WarningIn("meshToMesh0::calcAddressing()")
+                    << "Source patch " << fromPatch.name()
+                    << " has no faces. Not performing mapping for it."
+                    << endl;
+                boundaryAddressing_[patchi] = -1;
+            }
+            else
+            {
+                treeBoundBox wallBb(fromPatch.localPoints());
+                scalar typDim =
+                    wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
+
+                treeBoundBox shiftedBb
+                (
+                    wallBb.min(),
+                    wallBb.max() + vector(typDim, typDim, typDim)
+                );
+
+                // Note: allow more levels than in meshSearch. Assume patch
+                // is not as big as all boundary faces
+                indexedOctree<treeDataFace> oc
+                (
+                    treeDataFace(false, fromPatch),
+                    shiftedBb,  // overall search domain
+                    12,         // maxLevel
+                    10,         // leafsize
+                    6.0         // duplicity
+                );
+
+                const vectorField::subField centresToBoundary =
+                    toPatch.faceCentres();
+
+                boundaryAddressing_[patchi].setSize(toPatch.size());
+
+                scalar distSqr = sqr(wallBb.mag());
+
+                forAll(toPatch, toi)
+                {
+                    boundaryAddressing_[patchi][toi] = oc.findNearest
+                    (
+                        centresToBoundary[toi],
+                        distSqr
+                    ).index();
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateAddressing() : "
+            << "finished calculating mesh-to-mesh cell addressing" << endl;
+    }
+}
+
+
+void Foam::meshToMesh0::cellAddresses
+(
+    labelList& cellAddressing_,
+    const pointField& points,
+    const fvMesh& fromMesh,
+    const List<bool>& boundaryCell,
+    const indexedOctree<treeDataCell>& oc
+) const
+{
+    // the implemented search method is a simple neighbour array search.
+    // It starts from a cell zero, searches its neighbours and finds one
+    // which is nearer to the target point than the current position.
+    // The location of the "current position" is reset to that cell and
+    // search through the neighbours continues. The search is finished
+    // when all the neighbours of the cell are farther from the target
+    // point than the current cell
+
+    // set curCell label to zero (start)
+    register label curCell = 0;
+
+    // set reference to cell to cell addressing
+    const vectorField& centresFrom = fromMesh.cellCentres();
+    const labelListList& cc = fromMesh.cellCells();
+
+    forAll(points, toI)
+    {
+        // pick up target position
+        const vector& p = points[toI];
+
+        // set the sqr-distance
+        scalar distSqr = magSqr(p - centresFrom[curCell]);
+
+        bool closer;
+
+        do
+        {
+            closer = false;
+
+            // set the current list of neighbouring cells
+            const labelList& neighbours = cc[curCell];
+
+            forAll(neighbours, nI)
+            {
+                scalar curDistSqr =
+                    magSqr(p - centresFrom[neighbours[nI]]);
+
+                // search through all the neighbours.
+                // If the cell is closer, reset current cell and distance
+                if (curDistSqr < (1 - SMALL)*distSqr)
+                {
+                    curCell = neighbours[nI];
+                    distSqr = curDistSqr;
+                    closer = true;    // a closer neighbour has been found
+                }
+            }
+        } while (closer);
+
+        cellAddressing_[toI] = -1;
+
+        // Check point is actually in the nearest cell
+        if (fromMesh.pointInCell(p, curCell))
+        {
+            cellAddressing_[toI] = curCell;
+        }
+        else
+        {
+            // If curCell is a boundary cell then the point maybe either outside
+            // the domain or in an other region of the doamin, either way use
+            // the octree search to find it.
+            if (boundaryCell[curCell])
+            {
+                cellAddressing_[toI] = oc.findInside(p);
+            }
+            else
+            {
+                // If not on the boundary search the neighbours
+                bool found = false;
+
+                // set the current list of neighbouring cells
+                const labelList& neighbours = cc[curCell];
+
+                forAll(neighbours, nI)
+                {
+                    // search through all the neighbours.
+                    // If point is in neighbour reset current cell
+                    if (fromMesh.pointInCell(p, neighbours[nI]))
+                    {
+                        cellAddressing_[toI] = neighbours[nI];
+                        found = true;
+                        break;
+                    }
+                }
+
+                if (!found)
+                {
+                    // If still not found search the neighbour-neighbours
+
+                    // set the current list of neighbouring cells
+                    const labelList& neighbours = cc[curCell];
+
+                    forAll(neighbours, nI)
+                    {
+                        // set the current list of neighbour-neighbouring cells
+                        const labelList& nn = cc[neighbours[nI]];
+
+                        forAll(nn, nI)
+                        {
+                            // search through all the neighbours.
+                            // If point is in neighbour reset current cell
+                            if (fromMesh.pointInCell(p, nn[nI]))
+                            {
+                                cellAddressing_[toI] = nn[nI];
+                                found = true;
+                                break;
+                            }
+                        }
+                        if (found) break;
+                    }
+                }
+
+                if (!found)
+                {
+                    // Still not found so us the octree
+                    cellAddressing_[toI] = oc.findInside(p);
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C b/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
new file mode 100644
index 000000000..30d3c37f8
--- /dev/null
+++ b/src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
@@ -0,0 +1,274 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "meshToMesh0.H"
+#include "tetOverlapVolume.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::meshToMesh0::calculateInverseDistanceWeights() const
+{
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateInverseDistanceWeights() : "
+            << "calculating inverse distance weighting factors" << endl;
+    }
+
+    if (inverseDistanceWeightsPtr_)
+    {
+        FatalErrorIn("meshToMesh0::calculateInverseDistanceWeights()")
+            << "weighting factors already calculated"
+            << exit(FatalError);
+    }
+
+    //- Initialise overlap volume to zero
+    V_ = 0.0;
+
+    inverseDistanceWeightsPtr_ = new scalarListList(toMesh_.nCells());
+    scalarListList& invDistCoeffs = *inverseDistanceWeightsPtr_;
+
+    // get reference to source mesh data
+    const labelListList& cc = fromMesh_.cellCells();
+    const vectorField& centreFrom = fromMesh_.C().internalField();
+    const vectorField& centreTo = toMesh_.C().internalField();
+
+    forAll(cellAddressing_, celli)
+    {
+        if (cellAddressing_[celli] != -1)
+        {
+            const vector& target = centreTo[celli];
+            scalar m = mag(target - centreFrom[cellAddressing_[celli]]);
+
+            const labelList& neighbours = cc[cellAddressing_[celli]];
+
+            // if the nearest cell is a boundary cell or there is a direct hit,
+            // pick up the value
+            label directCelli = -1;
+            if (m < directHitTol || neighbours.empty())
+            {
+                directCelli = celli;
+            }
+            else
+            {
+                forAll(neighbours, ni)
+                {
+                    scalar nm = mag(target - centreFrom[neighbours[ni]]);
+                    if (nm < directHitTol)
+                    {
+                        directCelli = neighbours[ni];
+                        break;
+                    }
+                }
+            }
+
+
+            if (directCelli != -1)
+            {
+                // Direct hit
+                invDistCoeffs[directCelli].setSize(1);
+                invDistCoeffs[directCelli][0] = 1.0;
+                V_ += fromMesh_.V()[cellAddressing_[directCelli]];
+            }
+            else
+            {
+                invDistCoeffs[celli].setSize(neighbours.size() + 1);
+
+                // The first coefficient corresponds to the centre cell.
+                // The rest is ordered in the same way as the cellCells list.
+                scalar invDist = 1.0/m;
+                invDistCoeffs[celli][0] = invDist;
+                scalar sumInvDist = invDist;
+
+                // now add the neighbours
+                forAll(neighbours, ni)
+                {
+                    invDist = 1.0/mag(target - centreFrom[neighbours[ni]]);
+                    invDistCoeffs[celli][ni + 1] = invDist;
+                    sumInvDist += invDist;
+                }
+
+                // divide by the total inverse-distance
+                forAll(invDistCoeffs[celli], i)
+                {
+                    invDistCoeffs[celli][i] /= sumInvDist;
+                }
+
+
+                V_ +=
+                    invDistCoeffs[celli][0]
+                   *fromMesh_.V()[cellAddressing_[celli]];
+                for (label i = 1; i < invDistCoeffs[celli].size(); i++)
+                {
+                    V_ +=
+                        invDistCoeffs[celli][i]*fromMesh_.V()[neighbours[i-1]];
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::meshToMesh0::calculateInverseVolumeWeights() const
+{
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateInverseVolumeWeights() : "
+            << "calculating inverse volume weighting factors" << endl;
+    }
+
+    if (inverseVolumeWeightsPtr_)
+    {
+        FatalErrorIn("meshToMesh0::calculateInverseVolumeWeights()")
+            << "weighting factors already calculated"
+            << exit(FatalError);
+    }
+
+    //- Initialise overlap volume to zero
+    V_ = 0.0;
+
+    inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
+    scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
+
+    const labelListList& cellToCell = cellToCellAddressing();
+
+    tetOverlapVolume overlapEngine;
+
+    forAll(cellToCell, celli)
+    {
+        const labelList& overlapCells = cellToCell[celli];
+
+        if (overlapCells.size() > 0)
+        {
+            invVolCoeffs[celli].setSize(overlapCells.size());
+
+            forAll(overlapCells, j)
+            {
+                label cellFrom = overlapCells[j];
+                treeBoundBox bbFromMesh
+                (
+                    pointField
+                    (
+                        fromMesh_.points(),
+                        fromMesh_.cellPoints()[cellFrom]
+                    )
+                );
+
+                scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp
+                (
+                    toMesh_,
+                    celli,
+
+                    fromMesh_,
+                    cellFrom,
+                    bbFromMesh
+                );
+                invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
+
+                V_ += v;
+            }
+        }
+    }
+}
+
+
+void Foam::meshToMesh0::calculateCellToCellAddressing() const
+{
+    if (debug)
+    {
+        Info<< "meshToMesh0::calculateCellToCellAddressing() : "
+            << "calculating cell to cell addressing" << endl;
+    }
+
+    if (cellToCellAddressingPtr_)
+    {
+        FatalErrorIn("meshToMesh0::calculateCellToCellAddressing()")
+            << "addressing already calculated"
+            << exit(FatalError);
+    }
+
+    //- Initialise overlap volume to zero
+    V_ = 0.0;
+
+    tetOverlapVolume overlapEngine;
+
+    cellToCellAddressingPtr_ = new labelListList(toMesh_.nCells());
+    labelListList& cellToCell = *cellToCellAddressingPtr_;
+
+
+    forAll(cellToCell, iTo)
+    {
+        const labelList overLapCells =
+            overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo);
+        if (overLapCells.size() > 0)
+        {
+            //Info << "To " << iTo << endl;
+            //Info << "cellToCell " << overLapCells << endl;
+
+            cellToCell[iTo].setSize(overLapCells.size());
+            forAll(overLapCells, j)
+            {
+                cellToCell[iTo][j] = overLapCells[j];
+                V_ += fromMesh_.V()[overLapCells[j]];
+            }
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const Foam::scalarListList& Foam::meshToMesh0::inverseDistanceWeights() const
+{
+    if (!inverseDistanceWeightsPtr_)
+    {
+        calculateInverseDistanceWeights();
+    }
+
+    return *inverseDistanceWeightsPtr_;
+}
+
+
+const Foam::scalarListList& Foam::meshToMesh0::inverseVolumeWeights() const
+{
+    if (!inverseVolumeWeightsPtr_)
+    {
+        calculateInverseVolumeWeights();
+    }
+
+    return *inverseVolumeWeightsPtr_;
+}
+
+
+const Foam::labelListList& Foam::meshToMesh0::cellToCellAddressing() const
+{
+    if (!cellToCellAddressingPtr_)
+    {
+        calculateCellToCellAddressing();
+    }
+
+    return *cellToCellAddressingPtr_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMesh0/meshToMesh0.C b/src/sampling/meshToMesh0/meshToMesh0.C
new file mode 100644
index 000000000..164d749cd
--- /dev/null
+++ b/src/sampling/meshToMesh0/meshToMesh0.C
@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "meshToMesh0.H"
+#include "processorFvPatch.H"
+#include "demandDrivenData.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(meshToMesh0, 0);
+}
+
+const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-5;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshToMesh0::meshToMesh0
+(
+    const fvMesh& meshFrom,
+    const fvMesh& meshTo,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatchNames
+)
+:
+    fromMesh_(meshFrom),
+    toMesh_(meshTo),
+    patchMap_(patchMap),
+    cellAddressing_(toMesh_.nCells()),
+    boundaryAddressing_(toMesh_.boundaryMesh().size()),
+    inverseDistanceWeightsPtr_(NULL),
+    inverseVolumeWeightsPtr_(NULL),
+    cellToCellAddressingPtr_(NULL),
+    V_(0.0)
+{
+    forAll(fromMesh_.boundaryMesh(), patchi)
+    {
+        fromMeshPatches_.insert
+        (
+            fromMesh_.boundaryMesh()[patchi].name(),
+            patchi
+        );
+    }
+
+    forAll(toMesh_.boundaryMesh(), patchi)
+    {
+        toMeshPatches_.insert
+        (
+            toMesh_.boundaryMesh()[patchi].name(),
+            patchi
+        );
+    }
+
+    forAll(cuttingPatchNames, i)
+    {
+        if (toMeshPatches_.found(cuttingPatchNames[i]))
+        {
+            cuttingPatches_.insert
+            (
+                cuttingPatchNames[i],
+                toMeshPatches_.find(cuttingPatchNames[i])()
+            );
+        }
+        else
+        {
+            WarningIn
+            (
+                "meshToMesh0::meshToMesh0"
+                "(const fvMesh& meshFrom, const fvMesh& meshTo,"
+                "const HashTable<word>& patchMap,"
+                "const wordList& cuttingPatchNames)"
+            )   << "Cannot find cutting-patch " << cuttingPatchNames[i]
+                << " in destination mesh" << endl;
+        }
+    }
+
+    forAll(toMesh_.boundaryMesh(), patchi)
+    {
+        // Add the processor patches in the toMesh to the cuttingPatches list
+        if (isA<processorPolyPatch>(toMesh_.boundaryMesh()[patchi]))
+        {
+            cuttingPatches_.insert
+            (
+                toMesh_.boundaryMesh()[patchi].name(),
+                patchi
+            );
+        }
+    }
+
+    calcAddressing();
+}
+
+
+Foam::meshToMesh0::meshToMesh0
+(
+    const fvMesh& meshFrom,
+    const fvMesh& meshTo
+)
+:
+    fromMesh_(meshFrom),
+    toMesh_(meshTo),
+    cellAddressing_(toMesh_.nCells()),
+    boundaryAddressing_(toMesh_.boundaryMesh().size()),
+    inverseDistanceWeightsPtr_(NULL),
+    inverseVolumeWeightsPtr_(NULL),
+    cellToCellAddressingPtr_(NULL),
+    V_(0.0)
+{
+    // check whether both meshes have got the same number
+    // of boundary patches
+    if (fromMesh_.boundary().size() != toMesh_.boundary().size())
+    {
+        FatalErrorIn
+        (
+            "meshToMesh0::meshToMesh0"
+            "(const fvMesh& meshFrom, const fvMesh& meshTo)"
+        )   << "Incompatible meshes: different number of patches, "
+            << "fromMesh = " << fromMesh_.boundary().size()
+            << ", toMesh = " << toMesh_.boundary().size()
+            << exit(FatalError);
+    }
+
+    forAll(fromMesh_.boundaryMesh(), patchi)
+    {
+        if
+        (
+            fromMesh_.boundaryMesh()[patchi].name()
+         != toMesh_.boundaryMesh()[patchi].name()
+        )
+        {
+            FatalErrorIn
+            (
+                "meshToMesh0::meshToMesh0"
+                "(const fvMesh& meshFrom, const fvMesh& meshTo)"
+            )   << "Incompatible meshes: different patch names for patch "
+                << patchi
+                << ", fromMesh = " << fromMesh_.boundary()[patchi].name()
+                << ", toMesh = " << toMesh_.boundary()[patchi].name()
+                << exit(FatalError);
+        }
+
+        if
+        (
+            fromMesh_.boundaryMesh()[patchi].type()
+         != toMesh_.boundaryMesh()[patchi].type()
+        )
+        {
+            FatalErrorIn
+            (
+                "meshToMesh0::meshToMesh0"
+                "(const fvMesh& meshFrom, const fvMesh& meshTo)"
+            )   << "Incompatible meshes: different patch types for patch "
+                << patchi
+                << ", fromMesh = " << fromMesh_.boundary()[patchi].type()
+                << ", toMesh = " << toMesh_.boundary()[patchi].type()
+                << exit(FatalError);
+        }
+
+        fromMeshPatches_.insert
+        (
+            fromMesh_.boundaryMesh()[patchi].name(),
+            patchi
+        );
+
+        toMeshPatches_.insert
+        (
+            toMesh_.boundaryMesh()[patchi].name(),
+            patchi
+        );
+
+        patchMap_.insert
+        (
+            toMesh_.boundaryMesh()[patchi].name(),
+            fromMesh_.boundaryMesh()[patchi].name()
+        );
+    }
+
+    calcAddressing();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::meshToMesh0::~meshToMesh0()
+{
+    deleteDemandDrivenData(inverseDistanceWeightsPtr_);
+    deleteDemandDrivenData(inverseVolumeWeightsPtr_);
+    deleteDemandDrivenData(cellToCellAddressingPtr_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMesh0/meshToMesh0.H b/src/sampling/meshToMesh0/meshToMesh0.H
new file mode 100644
index 000000000..d2c7ff03e
--- /dev/null
+++ b/src/sampling/meshToMesh0/meshToMesh0.H
@@ -0,0 +1,373 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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::meshToMesh0
+
+Description
+    mesh to mesh interpolation class.
+
+Note
+    This class is due to be deprecated in favour of meshToMesh0New
+
+SourceFiles
+    meshToMesh0.C
+    calculateMeshToMesh0Addressing.C
+    calculateMeshToMesh0Weights.C
+    meshToMesh0Templates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshtoMesh_H
+#define meshtoMesh_H
+
+#include "fvMesh.H"
+#include "HashTable.H"
+#include "fvPatchMapper.H"
+#include "scalarList.H"
+#include "className.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type>
+class indexedOctree;
+
+class treeDataCell;
+
+/*---------------------------------------------------------------------------*\
+                         Class meshToMesh0 Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshToMesh0
+{
+    // Private data
+
+        // mesh references
+
+        const fvMesh& fromMesh_;
+        const fvMesh& toMesh_;
+
+        //- fromMesh patch labels
+        HashTable<label> fromMeshPatches_;
+
+        //- toMesh patch labels
+        HashTable<label> toMeshPatches_;
+
+        //- Patch map
+        HashTable<word> patchMap_;
+
+        //- toMesh patch labels which cut the from-mesh
+        HashTable<label> cuttingPatches_;
+
+        //- Cell addressing
+        labelList cellAddressing_;
+
+        //- Boundary addressing
+        labelListList boundaryAddressing_;
+
+        //- Inverse-distance interpolation weights
+        mutable scalarListList* inverseDistanceWeightsPtr_;
+
+        //- Inverse-volume interpolation weights
+        mutable scalarListList* inverseVolumeWeightsPtr_;
+
+        //- Cell to cell overlap addressing
+        mutable labelListList* cellToCellAddressingPtr_;
+
+        //- Overlap volume
+        mutable scalar V_;
+
+
+    // Private Member Functions
+
+        void calcAddressing();
+
+        void cellAddresses
+        (
+            labelList& cells,
+            const pointField& points,
+            const fvMesh& fromMesh,
+            const List<bool>& boundaryCell,
+            const indexedOctree<treeDataCell>& oc
+        ) const;
+
+        void calculateInverseDistanceWeights() const;
+
+        void calculateInverseVolumeWeights() const;
+
+        void calculateCellToCellAddressing() const;
+
+        const scalarListList& inverseDistanceWeights() const;
+
+        const scalarListList& inverseVolumeWeights() const;
+
+        const labelListList& cellToCellAddressing() const;
+
+
+    // Private static data members
+
+        //- Direct hit tolerance
+        static const scalar directHitTol;
+
+
+public:
+
+    // Declare name of the class and its debug switch
+    ClassName("meshToMesh0");
+
+
+    //- Enumeration specifying required accuracy
+    enum order
+    {
+        MAP,
+        INTERPOLATE,
+        CELL_POINT_INTERPOLATE,
+        CELL_VOLUME_WEIGHT
+    };
+
+
+    // Constructors
+
+        //- Construct from the two meshes, the patch name map for the patches
+        //  to be interpolated and the names of the toMesh-patches which
+        //  cut the fromMesh
+        meshToMesh0
+        (
+            const fvMesh& fromMesh,
+            const fvMesh& toMesh,
+            const HashTable<word>& patchMap,
+            const wordList& cuttingPatchNames
+        );
+
+        //- Construct from the two meshes assuming there is an exact mapping
+        //  between the patches
+        meshToMesh0
+        (
+            const fvMesh& fromMesh,
+            const fvMesh& toMesh
+        );
+
+
+    //- Destructor
+    ~meshToMesh0();
+
+
+    //- Patch-field interpolation class
+    class patchFieldInterpolator
+    :
+        public fvPatchFieldMapper
+    {
+        const labelList& directAddressing_;
+
+
+    public:
+
+        // Constructors
+
+            //- Construct given addressing
+            patchFieldInterpolator(const labelList& addr)
+            :
+                directAddressing_(addr)
+            {}
+
+
+        //- Destructor
+        virtual ~patchFieldInterpolator()
+        {}
+
+
+        // Member Functions
+
+            label size() const
+            {
+                return directAddressing_.size();
+            }
+
+            bool direct() const
+            {
+                return true;
+            }
+
+            bool hasUnmapped() const
+            {
+                return false;
+            }
+
+            const labelList& directAddressing() const
+            {
+                return directAddressing_;
+            }
+    };
+
+
+    // Member Functions
+
+        // Access
+
+            const fvMesh& fromMesh() const
+            {
+                return fromMesh_;
+            }
+
+            const fvMesh& toMesh() const
+            {
+                return toMesh_;
+            }
+
+            //- From toMesh cells to fromMesh cells
+            const labelList& cellAddressing() const
+            {
+                return cellAddressing_;
+            }
+
+            //- Overlap volume
+            scalar V() const
+            {
+                return V_;
+            }
+
+
+        // Interpolation
+
+            //- Map field
+            template<class Type, class CombineOp>
+            void mapField
+            (
+                Field<Type>&,
+                const Field<Type>&,
+                const labelList& adr,
+                const CombineOp& cop
+            ) const;
+
+            //- Interpolate field using inverse-distance weights
+            template<class Type, class CombineOp>
+            void interpolateField
+            (
+                Field<Type>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                const labelList& adr,
+                const scalarListList& weights,
+                const CombineOp& cop
+            ) const;
+
+            //- Interpolate field using inverse-volume weights
+            template<class Type, class CombineOp>
+            void interpolateField
+            (
+                Field<Type>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                const labelListList& adr,
+                const scalarListList& weights,
+                const CombineOp& cop
+            ) const;
+
+
+            //- Interpolate field using cell-point interpolation
+            template<class Type, class CombineOp>
+            void interpolateField
+            (
+                Field<Type>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                const labelList& adr,
+                const vectorField& centres,
+                const CombineOp& cop
+            )const;
+
+
+            //- Interpolate internal volume field
+            template<class Type, class CombineOp>
+            void interpolateInternalField
+            (
+                Field<Type>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+
+            template<class Type, class CombineOp>
+            void interpolateInternalField
+            (
+                Field<Type>&,
+                const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+
+
+            //- Interpolate volume field
+            template<class Type, class CombineOp>
+            void interpolate
+            (
+                GeometricField<Type, fvPatchField, volMesh>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+
+            template<class Type, class CombineOp>
+            void interpolate
+            (
+                GeometricField<Type, fvPatchField, volMesh>&,
+                const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+
+
+            //- Interpolate volume field
+            template<class Type, class CombineOp>
+            tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+            (
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+
+            template<class Type, class CombineOp>
+            tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
+            (
+                const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "meshToMesh0Templates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/meshToMesh0/meshToMesh0Templates.C b/src/sampling/meshToMesh0/meshToMesh0Templates.C
new file mode 100644
index 000000000..9d3bff680
--- /dev/null
+++ b/src/sampling/meshToMesh0/meshToMesh0Templates.C
@@ -0,0 +1,477 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "meshToMesh0.H"
+#include "volFields.H"
+#include "interpolationCellPoint.H"
+#include "SubField.H"
+#include "mixedFvPatchField.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::mapField
+(
+    Field<Type>& toF,
+    const Field<Type>& fromVf,
+    const labelList& adr,
+    const CombineOp& cop
+) const
+{
+    // Direct mapping of nearest-cell values
+
+    forAll(toF, celli)
+    {
+        if (adr[celli] != -1)
+        {
+            cop(toF[celli], fromVf[adr[celli]]);
+        }
+    }
+
+    //toF.map(fromVf, adr);
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolateField
+(
+    Field<Type>& toF,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    const labelListList& adr,
+    const scalarListList& weights,
+    const CombineOp& cop
+) const
+{
+    // Inverse volume weighted interpolation
+    forAll(toF, celli)
+    {
+        const labelList& overlapCells = adr[celli];
+        const scalarList& w = weights[celli];
+
+        Type f = pTraits<Type>::zero;
+        forAll(overlapCells, i)
+        {
+            label fromCelli = overlapCells[i];
+            f += fromVf[fromCelli]*w[i];
+            cop(toF[celli], f);
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolateField
+(
+    Field<Type>& toF,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    const labelList& adr,
+    const scalarListList& weights,
+    const CombineOp& cop
+) const
+{
+    // Inverse distance weighted interpolation
+
+    // get reference to cellCells
+    const labelListList& cc = fromMesh_.cellCells();
+
+    forAll(toF, celli)
+    {
+        if (adr[celli] != -1)
+        {
+            const labelList& neighbours = cc[adr[celli]];
+            const scalarList& w = weights[celli];
+
+            Type f = fromVf[adr[celli]]*w[0];
+
+            for (label ni = 1; ni < w.size(); ni++)
+            {
+                f += fromVf[neighbours[ni - 1]]*w[ni];
+            }
+
+            cop(toF[celli], f);
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolateField
+(
+    Field<Type>& toF,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    const labelList& adr,
+    const vectorField& centres,
+    const CombineOp& cop
+) const
+{
+    // Cell-Point interpolation
+    interpolationCellPoint<Type> interpolator(fromVf);
+
+    forAll(toF, celli)
+    {
+        if (adr[celli] != -1)
+        {
+            cop
+            (
+                toF[celli],
+                interpolator.interpolate
+                (
+                    centres[celli],
+                    adr[celli]
+                )
+            );
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolateInternalField
+(
+    Field<Type>& toF,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    if (fromVf.mesh() != fromMesh_)
+    {
+        FatalErrorIn
+        (
+            "meshToMesh0::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh0::order, const CombineOp&) const"
+        )   << "the argument field does not correspond to the right mesh. "
+            << "Field size: " << fromVf.size()
+            << " mesh size: " << fromMesh_.nCells()
+            << exit(FatalError);
+    }
+
+    if (toF.size() != toMesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "meshToMesh0::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh0::order, const CombineOp&) const"
+        )   << "the argument field does not correspond to the right mesh. "
+            << "Field size: " << toF.size()
+            << " mesh size: " << toMesh_.nCells()
+            << exit(FatalError);
+    }
+
+    switch(ord)
+    {
+        case MAP:
+            mapField(toF, fromVf, cellAddressing_, cop);
+        break;
+
+        case INTERPOLATE:
+        {
+            interpolateField
+            (
+                toF,
+                fromVf,
+                cellAddressing_,
+                inverseDistanceWeights(),
+                cop
+            );
+            break;
+        }
+        case CELL_POINT_INTERPOLATE:
+        {
+            interpolateField
+            (
+                toF,
+                fromVf,
+                cellAddressing_,
+                toMesh_.cellCentres(),
+                cop
+            );
+
+            break;
+        }
+        case CELL_VOLUME_WEIGHT:
+        {
+            const labelListList& cellToCell = cellToCellAddressing();
+            const scalarListList& invVolWeights = inverseVolumeWeights();
+
+            interpolateField
+            (
+                toF,
+                fromVf,
+                cellToCell,
+                invVolWeights,
+                cop
+            );
+            break;
+        }
+        default:
+            FatalErrorIn
+            (
+                "meshToMesh0::interpolateInternalField(Field<Type>&, "
+                "const GeometricField<Type, fvPatchField, volMesh>&, "
+                "meshToMesh0::order, const CombineOp&) const"
+            )   << "unknown interpolation scheme " << ord
+                << exit(FatalError);
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolateInternalField
+(
+    Field<Type>& toF,
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    interpolateInternalField(toF, tfromVf(), ord, cop);
+    tfromVf.clear();
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolate
+(
+    GeometricField<Type, fvPatchField, volMesh>& toVf,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    interpolateInternalField(toVf, fromVf, ord, cop);
+
+    forAll(toMesh_.boundaryMesh(), patchi)
+    {
+        const fvPatch& toPatch = toMesh_.boundary()[patchi];
+
+        if (cuttingPatches_.found(toPatch.name()))
+        {
+            switch(ord)
+            {
+                case MAP:
+                {
+                    mapField
+                    (
+                        toVf.boundaryField()[patchi],
+                        fromVf,
+                        boundaryAddressing_[patchi],
+                        cop
+                    );
+                    break;
+                }
+
+                case INTERPOLATE:
+                {
+                    interpolateField
+                    (
+                        toVf.boundaryField()[patchi],
+                        fromVf,
+                        boundaryAddressing_[patchi],
+                        toPatch.Cf(),
+                        cop
+                    );
+                    break;
+                }
+
+                case CELL_POINT_INTERPOLATE:
+                {
+                    interpolateField
+                    (
+                        toVf.boundaryField()[patchi],
+                        fromVf,
+                        boundaryAddressing_[patchi],
+                        toPatch.Cf(),
+                        cop
+                    );
+                    break;
+                }
+                case CELL_VOLUME_WEIGHT:
+                {
+                    // Do nothing
+                    break;
+                }
+
+                default:
+                    FatalErrorIn
+                    (
+                        "meshToMesh0::interpolate("
+                        "GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const GeometricField<Type, fvPatchField, volMesh>&, "
+                        "meshToMesh0::order, const CombineOp&) const"
+                    )   << "unknown interpolation scheme " << ord
+                        << exit(FatalError);
+            }
+
+            if (isA<mixedFvPatchField<Type> >(toVf.boundaryField()[patchi]))
+            {
+                refCast<mixedFvPatchField<Type> >
+                (
+                    toVf.boundaryField()[patchi]
+                ).refValue() = toVf.boundaryField()[patchi];
+            }
+        }
+        else if
+        (
+            patchMap_.found(toPatch.name())
+         && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
+        )
+        {
+            /*
+            toVf.boundaryField()[patchi].map
+            (
+                fromVf.boundaryField()
+                [
+                    fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
+                ],
+                boundaryAddressing_[patchi]
+            );
+            */
+
+            mapField
+            (
+                toVf.boundaryField()[patchi],
+                fromVf.boundaryField()
+                [
+                    fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
+                ],
+                boundaryAddressing_[patchi],
+                cop
+            );
+        }
+    }
+}
+
+
+template<class Type, class CombineOp>
+void Foam::meshToMesh0::interpolate
+(
+    GeometricField<Type, fvPatchField, volMesh>& toVf,
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    interpolate(toVf, tfromVf(), ord, cop);
+    tfromVf.clear();
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMesh0::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    // Create and map the internal-field values
+    Field<Type> internalField(toMesh_.nCells());
+    interpolateInternalField(internalField, fromVf, ord, cop);
+
+    // check whether both meshes have got the same number
+    // of boundary patches
+    if (fromMesh_.boundary().size() != toMesh_.boundary().size())
+    {
+        FatalErrorIn
+        (
+            "meshToMesh0::interpolate"
+            "(const GeometricField<Type, fvPatchField, volMesh>&,"
+            "meshToMesh0::order, const CombineOp&) const"
+        )   << "Incompatible meshes: different number of boundaries, "
+               "only internal field may be interpolated"
+            << exit(FatalError);
+    }
+
+    // Create and map the patch field values
+    PtrList<fvPatchField<Type> > patchFields
+    (
+        boundaryAddressing_.size()
+    );
+
+    forAll(boundaryAddressing_, patchI)
+    {
+        patchFields.set
+        (
+            patchI,
+            fvPatchField<Type>::New
+            (
+                fromVf.boundaryField()[patchI],
+                toMesh_.boundary()[patchI],
+                DimensionedField<Type, volMesh>::null(),
+                patchFieldInterpolator
+                (
+                    boundaryAddressing_[patchI]
+                )
+            )
+        );
+    }
+
+
+    // Create the complete field from the pieces
+    tmp<GeometricField<Type, fvPatchField, volMesh> > ttoF
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                "interpolated(" + fromVf.name() + ')',
+                toMesh_.time().timeName(),
+                toMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            toMesh_,
+            fromVf.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+
+    return ttoF;
+}
+
+
+template<class Type, class CombineOp>
+Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::meshToMesh0::interpolate
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
+    meshToMesh0::order ord,
+    const CombineOp& cop
+) const
+{
+    tmp<GeometricField<Type, fvPatchField, volMesh> > tint =
+        interpolate(tfromVf(), ord, cop);
+    tfromVf.clear();
+
+    return tint;
+}
+
+
+// ************************************************************************* //