From 0f2e8be91c6bea6de15e1bb7f828e7ec4eac683b Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 22 Dec 2011 11:39:54 +0000
Subject: [PATCH] ENH: mapFields: add -subtract option to subtract whilst
 mapping

---
 .../mapFields/MapConsistentVolFields.H        |  20 +-
 .../preProcessing/mapFields/MapMeshes.H       | 263 ++++++++++++++++++
 .../preProcessing/mapFields/MapVolFields.H    |  13 +-
 .../preProcessing/mapFields/mapFields.C       | 217 ++++++---------
 .../meshToMesh/meshToMesh.H                   |  45 +--
 .../meshToMesh/meshToMeshInterpolate.C        | 121 ++++----
 6 files changed, 476 insertions(+), 203 deletions(-)
 create mode 100644 applications/utilities/preProcessing/mapFields/MapMeshes.H

diff --git a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
index 2a203eac48..9015c473a2 100644
--- a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
@@ -35,12 +35,13 @@ License
 namespace Foam
 {
 
-template<class Type>
+template<class Type, class CombineOp>
 void MapConsistentVolFields
 (
     const IOobjectList& objects,
     const meshToMesh& meshToMeshInterp,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const CombineOp& cop
 )
 {
     const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@@ -84,7 +85,13 @@ void MapConsistentVolFields
             );
 
             // Interpolate field
-            meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
+            meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
 
             // Write field
             fieldTarget.write();
@@ -97,7 +104,12 @@ void MapConsistentVolFields
             GeometricField<Type, fvPatchField, volMesh> fieldTarget
             (
                 fieldTargetIOobject,
-                meshToMeshInterp.interpolate(fieldSource, mapOrder)
+                meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+                (
+                    fieldSource,
+                    mapOrder,
+                    cop
+                )
             );
 
             // Write field
diff --git a/applications/utilities/preProcessing/mapFields/MapMeshes.H b/applications/utilities/preProcessing/mapFields/MapMeshes.H
new file mode 100644
index 0000000000..e497eb0ccd
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFields/MapMeshes.H
@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "MapConsistentVolFields.H"
+#include "mapLagrangian.H"
+#include "UnMapped.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<template<class> class CombineOp>
+void MapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh::order& mapOrder
+)
+{
+    // Create the interpolation scheme
+    meshToMesh meshToMeshInterp(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
+        // ~~~~~~~~~~~~~
+        MapConsistentVolFields<scalar>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<scalar>()
+        );
+        MapConsistentVolFields<vector>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<vector>()
+        );
+        MapConsistentVolFields<sphericalTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<sphericalTensor>()
+        );
+        MapConsistentVolFields<symmTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<symmTensor>()
+        );
+        MapConsistentVolFields<tensor>
+        (
+            objects,
+            meshToMeshInterp,
+            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(meshToMeshInterp);
+}
+
+
+template<template<class> class CombineOp>
+void MapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const meshToMesh::order& mapOrder
+)
+{
+    // Create the interpolation scheme
+    meshToMesh meshToMeshInterp
+    (
+        meshSource,
+        meshTarget,
+        patchMap,
+        cuttingPatches
+    );
+
+    Info<< nl
+        << "Mapping fields for time " << meshSource.time().timeName()
+        << nl << endl;
+
+    {
+        // Search for list of source objects for this time
+        IOobjectList objects(meshSource, meshSource.time().timeName());
+
+        // Map volFields
+        // ~~~~~~~~~~~~~
+        MapVolFields<scalar>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<scalar>()
+        );
+        MapVolFields<vector>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<vector>()
+        );
+        MapVolFields<sphericalTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<sphericalTensor>()
+        );
+        MapVolFields<symmTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<symmTensor>()
+        );
+        MapVolFields<tensor>
+        (
+            objects,
+            meshToMeshInterp,
+            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(meshToMeshInterp);
+}
+
+
+template<template<class> class CombineOp>
+void MapConsistentSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh::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
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFields/MapVolFields.H b/applications/utilities/preProcessing/mapFields/MapVolFields.H
index 64547a54de..65eb44d12d 100644
--- a/applications/utilities/preProcessing/mapFields/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapVolFields.H
@@ -35,12 +35,13 @@ License
 namespace Foam
 {
 
-template<class Type>
+template<class Type, class CombineOp>
 void MapVolFields
 (
     const IOobjectList& objects,
     const meshToMesh& meshToMeshInterp,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const CombineOp& cop
 )
 {
     const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@@ -84,7 +85,13 @@ void MapVolFields
             );
 
             // Interpolate field
-            meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
+            meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
 
             // Write field
             fieldTarget.write();
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 8d70911cab..9571be6fb0 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -34,11 +34,8 @@ Description
 
 #include "fvCFD.H"
 #include "meshToMesh.H"
-#include "MapVolFields.H"
-#include "MapConsistentVolFields.H"
-#include "UnMapped.H"
 #include "processorFvPatch.H"
-#include "mapLagrangian.H"
+#include "MapMeshes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,56 +43,28 @@ void mapConsistentMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    // Create the interpolation scheme
-    meshToMesh meshToMeshInterp(meshSource, meshTarget);
-
-    Info<< nl
-        << "Consistently creating and mapping fields for time "
-        << meshSource.time().timeName() << nl << endl;
-
+    if (subtract)
     {
-        // Search for list of objects for this time
-        IOobjectList objects(meshSource, meshSource.time().timeName());
-
-        // Map volFields
-        // ~~~~~~~~~~~~~
-        MapConsistentVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<vector>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<sphericalTensor>
+        MapConsistentMesh<minusEqOp>
         (
-            objects,
-            meshToMeshInterp,
+            meshSource,
+            meshTarget,
             mapOrder
         );
-        MapConsistentVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
     }
-
+    else
     {
-        // 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);
+        MapConsistentMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
     }
-
-    mapLagrangian(meshToMeshInterp);
 }
 
 
@@ -105,57 +74,32 @@ void mapSubMesh
     const fvMesh& meshTarget,
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    // Create the interpolation scheme
-    meshToMesh meshToMeshInterp
-    (
-        meshSource,
-        meshTarget,
-        patchMap,
-        cuttingPatches
-    );
-
-    Info<< nl
-        << "Mapping fields for time " << meshSource.time().timeName()
-        << nl << endl;
-
+    if (subtract)
     {
-        // Search for list of source objects for this time
-        IOobjectList objects(meshSource, meshSource.time().timeName());
-
-        // Map volFields
-        // ~~~~~~~~~~~~~
-        MapVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<vector>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<sphericalTensor>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
+        MapSubMesh<minusEqOp>
+        (
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
+        );
     }
-
+    else
     {
-        // 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);
+        MapSubMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
+        );
     }
-
-    mapLagrangian(meshToMeshInterp);
 }
 
 
@@ -163,40 +107,28 @@ void mapConsistentSubMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    HashTable<word> patchMap;
-    HashTable<label> cuttingPatchTable;
-
-    forAll(meshTarget.boundary(), patchi)
+    if (subtract)
     {
-        if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
-        {
-            patchMap.insert
-            (
-                meshTarget.boundary()[patchi].name(),
-                meshTarget.boundary()[patchi].name()
-            );
-        }
-        else
-        {
-            cuttingPatchTable.insert
-            (
-                meshTarget.boundaryMesh()[patchi].name(),
-                -1
-            );
-        }
+        MapConsistentSubMesh<minusEqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
+    }
+    else
+    {
+        MapConsistentSubMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
     }
-
-    mapSubMesh
-    (
-        meshSource,
-        meshTarget,
-        patchMap,
-        cuttingPatchTable.toc(),
-        mapOrder
-    );
 }
 
 
@@ -288,6 +220,11 @@ int main(int argc, char *argv[])
         "word",
         "specify the mapping method"
     );
+    argList::addBoolOption
+    (
+        "subtract",
+        "subtract mapped source from target"
+    );
 
     argList args(argc, argv);
 
@@ -350,6 +287,13 @@ int main(int argc, char *argv[])
         Info<< "Mapping method: " << mapMethod << endl;
     }
 
+    const bool subtract = args.optionFound("subtract");
+    if (subtract)
+    {
+        Info<< "Subtracting mapped source field from target" << endl;
+    }
+
+
     #include "createTimes.H"
 
     HashTable<word> patchMap;
@@ -431,7 +375,13 @@ int main(int argc, char *argv[])
 
             if (consistent)
             {
-                mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
             }
             else
             {
@@ -441,7 +391,8 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     cuttingPatches,
-                    mapOrder
+                    mapOrder,
+                    subtract
                 );
             }
         }
@@ -503,7 +454,13 @@ int main(int argc, char *argv[])
 
             if (consistent)
             {
-                mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
             }
             else
             {
@@ -513,7 +470,8 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     addProcessorPatches(meshTarget, cuttingPatches),
-                    mapOrder
+                    mapOrder,
+                    subtract
                 );
             }
         }
@@ -629,7 +587,8 @@ int main(int argc, char *argv[])
                             (
                                 meshSource,
                                 meshTarget,
-                                mapOrder
+                                mapOrder,
+                                subtract
                             );
                         }
                         else
@@ -640,7 +599,8 @@ int main(int argc, char *argv[])
                                 meshTarget,
                                 patchMap,
                                 addProcessorPatches(meshTarget, cuttingPatches),
-                                mapOrder
+                                mapOrder,
+                                subtract
                             );
                         }
                     }
@@ -679,7 +639,7 @@ int main(int argc, char *argv[])
 
         if (consistent)
         {
-            mapConsistentMesh(meshSource, meshTarget, mapOrder);
+            mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
         }
         else
         {
@@ -689,7 +649,8 @@ int main(int argc, char *argv[])
                 meshTarget,
                 patchMap,
                 cuttingPatches,
-                mapOrder
+                mapOrder,
+                subtract
             );
         }
     }
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 785829cc60..83d1ee902a 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -219,84 +219,93 @@ public:
         // Interpolation
 
             //- Map field
-            template<class Type>
+            template<class Type, class CombineOp>
             void mapField
             (
                 Field<Type>&,
                 const Field<Type>&,
-                const labelList& adr
+                const labelList& adr,
+                const CombineOp& cop
             ) const;
 
             //- Interpolate field using inverse-distance weights
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const scalarListList& weights
+                const scalarListList& weights,
+                const CombineOp& cop
             ) const;
 
             //- Interpolate field using cell-point interpolation
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const vectorField& centres
+                const vectorField& centres,
+                const CombineOp& cop
             ) const;
 
 
             //- Interpolate internal volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
             (
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
             (
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 };
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
index 991ae31826..7ed543ed3c 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
@@ -31,12 +31,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::mapField
 (
     Field<Type>& toF,
     const Field<Type>& fromVf,
-    const labelList& adr
+    const labelList& adr,
+    const CombineOp& cop
 ) const
 {
     // Direct mapping of nearest-cell values
@@ -45,7 +46,7 @@ void Foam::meshToMesh::mapField
     {
         if (adr[celli] != -1)
         {
-            toF[celli] = fromVf[adr[celli]];
+            cop(toF[celli], fromVf[adr[celli]]);
         }
     }
 
@@ -53,13 +54,14 @@ void Foam::meshToMesh::mapField
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const scalarListList& weights
+    const scalarListList& weights,
+    const CombineOp& cop
 ) const
 {
     // Inverse distance weighted interpolation
@@ -74,24 +76,27 @@ void Foam::meshToMesh::interpolateField
             const labelList& neighbours = cc[adr[celli]];
             const scalarList& w = weights[celli];
 
-            toF[celli] = fromVf[adr[celli]]*w[0];
+            Type f = fromVf[adr[celli]]*w[0];
 
             for (label ni = 1; ni < w.size(); ni++)
             {
-                toF[celli] += fromVf[neighbours[ni - 1]]*w[ni];
+                f += fromVf[neighbours[ni - 1]]*w[ni];
             }
+
+            cop(toF[celli], f);
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const vectorField& centres
+    const vectorField& centres,
+    const CombineOp& cop
 ) const
 {
     // Cell-Point interpolation
@@ -101,31 +106,36 @@ void Foam::meshToMesh::interpolateField
     {
         if (adr[celli] != -1)
         {
-            toF[celli] = interpolator.interpolate
+            cop
             (
-                centres[celli],
-                adr[celli]
+                toF[celli],
+                interpolator.interpolate
+                (
+                    centres[celli],
+                    adr[celli]
+                )
             );
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateInternalField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     if (fromVf.mesh() != fromMesh_)
     {
         FatalErrorIn
         (
-            "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-            "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-            "meshToMesh::order ord) const"
+            "meshToMesh::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh::order, const CombineOp&) const"
         )   << "the argument field does not correspond to the right mesh. "
             << "Field size: " << fromVf.size()
             << " mesh size: " << fromMesh_.nCells()
@@ -136,9 +146,9 @@ void Foam::meshToMesh::interpolateInternalField
     {
         FatalErrorIn
         (
-            "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-            "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-            "meshToMesh::order ord) const"
+            "meshToMesh::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh::order, const CombineOp&) const"
         )   << "the argument field does not correspond to the right mesh. "
             << "Field size: " << toF.size()
             << " mesh size: " << toMesh_.nCells()
@@ -148,7 +158,7 @@ void Foam::meshToMesh::interpolateInternalField
     switch(ord)
     {
         case MAP:
-            mapField(toF, fromVf, cellAddressing_);
+            mapField(toF, fromVf, cellAddressing_, cop);
         break;
 
         case INTERPOLATE:
@@ -157,7 +167,8 @@ void Foam::meshToMesh::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                inverseDistanceWeights()
+                inverseDistanceWeights(),
+                cop
             );
         break;
 
@@ -167,44 +178,47 @@ void Foam::meshToMesh::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                toMesh_.cellCentres()
+                toMesh_.cellCentres(),
+                cop
             );
         break;
 
         default:
             FatalErrorIn
             (
-                "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-                "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-                "meshToMesh::order ord) const"
+                "meshToMesh::interpolateInternalField(Field<Type>&, "
+                "const GeometricField<Type, fvPatchField, volMesh>&, "
+                "meshToMesh::order, const CombineOp&) const"
             )   << "unknown interpolation scheme " << ord
                 << exit(FatalError);
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateInternalField
 (
     Field<Type>& toF,
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolateInternalField(toF, tfromVf(), ord);
+    interpolateInternalField(toF, tfromVf(), ord, cop);
     tfromVf.clear();
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolateInternalField(toVf, fromVf, ord);
+    interpolateInternalField(toVf, fromVf, ord, cop);
 
     forAll(toMesh_.boundaryMesh(), patchi)
     {
@@ -219,7 +233,8 @@ void Foam::meshToMesh::interpolate
                     (
                         toVf.boundaryField()[patchi],
                         fromVf,
-                        boundaryAddressing_[patchi]
+                        boundaryAddressing_[patchi],
+                        cop
                     );
                 break;
 
@@ -229,7 +244,8 @@ void Foam::meshToMesh::interpolate
                         toVf.boundaryField()[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf()
+                        toPatch.Cf(),
+                        cop
                     );
                 break;
 
@@ -239,7 +255,8 @@ void Foam::meshToMesh::interpolate
                         toVf.boundaryField()[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf()
+                        toPatch.Cf(),
+                        cop
                     );
                 break;
 
@@ -247,9 +264,9 @@ void Foam::meshToMesh::interpolate
                     FatalErrorIn
                     (
                         "meshToMesh::interpolate("
-                        "GeometricField<Type, fvPatchField, volMesh>& toVf, "
-                        "const GeometricField<Type, fvPatchField, volMesh>& "
-                        "fromVf, meshToMesh::order ord) const"
+                        "GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const GeometricField<Type, fvPatchField, volMesh>&, "
+                        "meshToMesh::order, const CombineOp&) const"
                     )   << "unknown interpolation scheme " << ord
                         << exit(FatalError);
             }
@@ -286,37 +303,40 @@ void Foam::meshToMesh::interpolate
                 [
                     fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
                 ],
-                boundaryAddressing_[patchi]
+                boundaryAddressing_[patchi],
+                cop
             );
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolate(toVf, tfromVf(), ord);
+    interpolate(toVf, tfromVf(), ord, cop);
     tfromVf.clear();
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::meshToMesh::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     // Create and map the internal-field values
     Field<Type> internalField(toMesh_.nCells());
-    interpolateInternalField(internalField, fromVf, ord);
+    interpolateInternalField(internalField, fromVf, ord, cop);
 
     // check whether both meshes have got the same number
     // of boundary patches
@@ -325,8 +345,8 @@ Foam::meshToMesh::interpolate
         FatalErrorIn
         (
             "meshToMesh::interpolate"
-            "(const GeometricField<Type, fvPatchField, volMesh>& fromVf,"
-            "meshToMesh::order ord) const"
+            "(const GeometricField<Type, fvPatchField, volMesh>&,"
+            "meshToMesh::order, const CombineOp&) const"
         )   << "Incompatible meshes: different number of boundaries, "
                "only internal field may be interpolated"
             << exit(FatalError);
@@ -381,16 +401,17 @@ Foam::meshToMesh::interpolate
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::meshToMesh::interpolate
 (
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh> > tint =
-        interpolate(tfromVf(), ord);
+        interpolate(tfromVf(), ord, cop);
     tfromVf.clear();
 
     return tint;