diff --git a/applications/utilities/preProcessing/mapFields/Make/files b/applications/utilities/preProcessing/mapFields/Make/files
index a779c1591c..4503a5fa98 100644
--- a/applications/utilities/preProcessing/mapFields/Make/files
+++ b/applications/utilities/preProcessing/mapFields/Make/files
@@ -1,4 +1,5 @@
 mapLagrangian.C
+mapMeshes.C
 mapFields.C
 
 EXE = $(FOAM_APPBIN)/mapFields
diff --git a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
index 8dcaffa00f..8f6c36d72b 100644
--- a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,13 +35,12 @@ License
 namespace Foam
 {
 
-template<class Type, class CombineOp>
+template<class Type>
 void MapConsistentVolFields
 (
     const IOobjectList& objects,
     const meshToMesh0& meshToMesh0Interp,
-    const meshToMesh0::order& mapOrder,
-    const CombineOp& cop
+    const meshToMesh0::order& mapOrder
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
@@ -78,12 +77,11 @@ void MapConsistentVolFields
             );
 
             // Interpolate field
-            meshToMesh0Interp.interpolate
+            meshToMesh0Interp.interpolate<Type>
             (
                 fieldTarget,
                 fieldSource,
-                mapOrder,
-                cop
+                mapOrder
             );
 
             // Write field
@@ -100,8 +98,7 @@ void MapConsistentVolFields
                 meshToMesh0Interp.interpolate
                 (
                     fieldSource,
-                    mapOrder,
-                    cop
+                    mapOrder
                 )
             );
 
diff --git a/applications/utilities/preProcessing/mapFields/MapVolFields.H b/applications/utilities/preProcessing/mapFields/MapVolFields.H
index a570b36faa..61b5448a55 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     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,13 +35,12 @@ License
 namespace Foam
 {
 
-template<class Type, class CombineOp>
+template<class Type>
 void MapVolFields
 (
     const IOobjectList& objects,
     const meshToMesh0& meshToMesh0Interp,
-    const meshToMesh0::order& mapOrder,
-    const CombineOp& cop
+    const meshToMesh0::order& mapOrder
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
@@ -81,8 +80,7 @@ void MapVolFields
             (
                 fieldTarget,
                 fieldSource,
-                mapOrder,
-                cop
+                mapOrder
             );
 
             // Write field
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 83d9aa7ca9..7971b6031b 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,148 +33,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "fvMesh.H"
-#include "surfaceMesh.H"
+#include "mapMeshes.H"
 #include "decompositionMethod.H"
-#include "meshToMesh0.H"
-#include "processorFvPatch.H"
-#include "MapMeshes.H"
 
 using namespace Foam;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void mapConsistentMesh
-(
-    const fvMesh& meshSource,
-    const fvMesh& meshTarget,
-    const meshToMesh0::order& mapOrder,
-    const bool subtract
-)
-{
-    if (subtract)
-    {
-        MapConsistentMesh<minusEqOp>
-        (
-            meshSource,
-            meshTarget,
-            mapOrder
-        );
-    }
-    else
-    {
-        MapConsistentMesh<eqOp>
-        (
-            meshSource,
-            meshTarget,
-            mapOrder
-        );
-    }
-}
-
-
-void mapSubMesh
-(
-    const fvMesh& meshSource,
-    const fvMesh& meshTarget,
-    const HashTable<word>& patchMap,
-    const wordList& cuttingPatches,
-    const meshToMesh0::order& mapOrder,
-    const bool subtract
-)
-{
-    if (subtract)
-    {
-        MapSubMesh<minusEqOp>
-        (
-            meshSource,
-            meshTarget,
-            patchMap,
-            cuttingPatches,
-            mapOrder
-        );
-    }
-    else
-    {
-        MapSubMesh<eqOp>
-        (
-            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
-        );
-    }
-}
-
-
-wordList addProcessorPatches
-(
-    const fvMesh& meshTarget,
-    const wordList& cuttingPatches
-)
-{
-    // Add the processor patches to the cutting list
-    HashTable<label> cuttingPatchTable;
-    forAll(cuttingPatches, i)
-    {
-        cuttingPatchTable.insert(cuttingPatches[i], i);
-    }
-
-    forAll(meshTarget.boundary(), patchi)
-    {
-        if (isA<processorFvPatch>(meshTarget.boundary()[patchi]))
-        {
-            if
-            (
-               !cuttingPatchTable.found
-                (
-                    meshTarget.boundaryMesh()[patchi].name()
-                )
-            )
-            {
-                cuttingPatchTable.insert
-                (
-                    meshTarget.boundaryMesh()[patchi].name(),
-                    -1
-                );
-            }
-        }
-    }
-
-    return cuttingPatchTable.toc();
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -225,11 +88,6 @@ int main(int argc, char *argv[])
         "word",
         "specify the mapping method"
     );
-    argList::addBoolOption
-    (
-        "subtract",
-        "subtract mapped source from target"
-    );
 
     argList args(argc, argv);
 
@@ -292,12 +150,6 @@ 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"
 
@@ -380,8 +232,7 @@ int main(int argc, char *argv[])
                 (
                     meshSource,
                     meshTarget,
-                    mapOrder,
-                    subtract
+                    mapOrder
                 );
             }
             else
@@ -392,8 +243,7 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     cuttingPatches,
-                    mapOrder,
-                    subtract
+                    mapOrder
                 );
             }
         }
@@ -455,8 +305,7 @@ int main(int argc, char *argv[])
                 (
                     meshSource,
                     meshTarget,
-                    mapOrder,
-                    subtract
+                    mapOrder
                 );
             }
             else
@@ -467,8 +316,7 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     addProcessorPatches(meshTarget, cuttingPatches),
-                    mapOrder,
-                    subtract
+                    mapOrder
                 );
             }
         }
@@ -567,8 +415,7 @@ int main(int argc, char *argv[])
                             (
                                 meshSource,
                                 meshTarget,
-                                mapOrder,
-                                subtract
+                                mapOrder
                             );
                         }
                         else
@@ -579,8 +426,7 @@ int main(int argc, char *argv[])
                                 meshTarget,
                                 patchMap,
                                 addProcessorPatches(meshTarget, cuttingPatches),
-                                mapOrder,
-                                subtract
+                                mapOrder
                             );
                         }
                     }
@@ -621,7 +467,7 @@ int main(int argc, char *argv[])
 
         if (consistent)
         {
-            mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
+            mapConsistentMesh(meshSource, meshTarget, mapOrder);
         }
         else
         {
@@ -631,8 +477,7 @@ int main(int argc, char *argv[])
                 meshTarget,
                 patchMap,
                 cuttingPatches,
-                mapOrder,
-                subtract
+                mapOrder
             );
         }
     }
diff --git a/applications/utilities/preProcessing/mapFields/MapMeshes.H b/applications/utilities/preProcessing/mapFields/mapMeshes.C
similarity index 82%
rename from applications/utilities/preProcessing/mapFields/MapMeshes.H
rename to applications/utilities/preProcessing/mapFields/mapMeshes.C
index 164f9c5af4..6c321d3fc2 100644
--- a/applications/utilities/preProcessing/mapFields/MapMeshes.H
+++ b/applications/utilities/preProcessing/mapFields/mapMeshes.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,21 +23,17 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MapMeshes_H
-#define MapMeshes_H
-
+#include "mapMeshes.H"
+#include "surfaceMesh.H"
+#include "processorFvPatch.H"
+#include "mapLagrangian.H"
 #include "MapVolFields.H"
 #include "MapConsistentVolFields.H"
-#include "mapLagrangian.H"
 #include "UnMapped.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-template<template<class> class CombineOp>
-void MapConsistentMesh
+void Foam::mapConsistentMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
@@ -61,36 +57,31 @@ void MapConsistentMesh
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<scalar>()
+            mapOrder
         );
         MapConsistentVolFields<vector>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<vector>()
+            mapOrder
         );
         MapConsistentVolFields<sphericalTensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<sphericalTensor>()
+            mapOrder
         );
         MapConsistentVolFields<symmTensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<symmTensor>()
+            mapOrder
         );
         MapConsistentVolFields<tensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<tensor>()
+            mapOrder
         );
     }
 
@@ -119,8 +110,7 @@ void MapConsistentMesh
 }
 
 
-template<template<class> class CombineOp>
-void MapSubMesh
+void Foam::mapSubMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
@@ -152,36 +142,31 @@ void MapSubMesh
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<scalar>()
+            mapOrder
         );
         MapVolFields<vector>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<vector>()
+            mapOrder
         );
         MapVolFields<sphericalTensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<sphericalTensor>()
+            mapOrder
         );
         MapVolFields<symmTensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<symmTensor>()
+            mapOrder
         );
         MapVolFields<tensor>
         (
             objects,
             meshToMesh0Interp,
-            mapOrder,
-            CombineOp<tensor>()
+            mapOrder
         );
     }
 
@@ -210,8 +195,7 @@ void MapSubMesh
 }
 
 
-template<template<class> class CombineOp>
-void MapConsistentSubMesh
+void Foam::mapConsistentSubMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
@@ -241,7 +225,7 @@ void MapConsistentSubMesh
         }
     }
 
-    MapSubMesh<CombineOp>
+    mapSubMesh
     (
         meshSource,
         meshTarget,
@@ -252,12 +236,42 @@ void MapConsistentSubMesh
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+Foam::wordList Foam::addProcessorPatches
+(
+    const fvMesh& meshTarget,
+    const wordList& cuttingPatches
+)
+{
+    // Add the processor patches to the cutting list
+    HashTable<label> cuttingPatchTable;
+    forAll(cuttingPatches, i)
+    {
+        cuttingPatchTable.insert(cuttingPatches[i], i);
+    }
 
-} // End namespace Foam
+    forAll(meshTarget.boundary(), patchi)
+    {
+        if (isA<processorFvPatch>(meshTarget.boundary()[patchi]))
+        {
+            if
+            (
+               !cuttingPatchTable.found
+                (
+                    meshTarget.boundaryMesh()[patchi].name()
+                )
+            )
+            {
+                cuttingPatchTable.insert
+                (
+                    meshTarget.boundaryMesh()[patchi].name(),
+                    -1
+                );
+            }
+        }
+    }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    return cuttingPatchTable.toc();
+}
 
-#endif
 
 // ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFields/mapMeshes.H b/applications/utilities/preProcessing/mapFields/mapMeshes.H
new file mode 100644
index 0000000000..3a67f46f6b
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFields/mapMeshes.H
@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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 "meshToMesh0.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void mapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh0::order& mapOrder
+);
+
+void mapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const meshToMesh0::order& mapOrder
+);
+
+
+void mapConsistentSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh0::order& mapOrder
+);
+
+
+wordList addProcessorPatches
+(
+    const fvMesh& meshTarget,
+    const wordList& cuttingPatches
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/Make/files b/applications/utilities/preProcessing/mapFieldsPar/Make/files
index 63742e0e32..c2519449c4 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/Make/files
+++ b/applications/utilities/preProcessing/mapFieldsPar/Make/files
@@ -1,3 +1,4 @@
+mapMeshes.C
 mapLagrangian.C
 mapFieldsPar.C
 
diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H b/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
index fae37a73a3..79c1f080e4 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -117,13 +117,12 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void MapVolFields
 (
     const IOobjectList& objects,
     const HashSet<word>& selectedFields,
-    const meshToMesh& interp,
-    const CombineOp& cop
+    const meshToMesh& interp
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
@@ -155,7 +154,7 @@ void MapVolFields
                     << fieldName << endl;
                 fieldType fieldTarget(targetIO, meshTarget);
 
-                interp.mapSrcToTgt(fieldSource, cop, fieldTarget);
+                interp.mapSrcToTgt(fieldSource, fieldTarget);
 
                 evaluateConstraintTypes(fieldTarget);
 
@@ -168,8 +167,7 @@ void MapVolFields
 
                 targetIO.readOpt() = IOobject::NO_READ;
 
-                tmp<fieldType>
-                    tfieldTarget(interp.mapSrcToTgt(fieldSource, cop));
+                tmp<fieldType> tfieldTarget(interp.mapSrcToTgt(fieldSource));
 
                 fieldType fieldTarget(targetIO, tfieldTarget);
 
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
index 7d886d8bc9..ee9976773f 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,122 +30,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
-#include "meshToMesh.H"
+#include "argList.H"
+#include "mapMeshes.H"
 #include "cellVolumeWeightMethod.H"
-#include "processorPolyPatch.H"
-#include "MapMeshes.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void mapConsistentMesh
-(
-    const fvMesh& meshSource,
-    const fvMesh& meshTarget,
-    const word& 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 word& 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();
-}
 
+using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -187,11 +76,6 @@ int main(int argc, char *argv[])
         "word",
         "specify the mapping method"
     );
-    argList::addBoolOption
-    (
-        "subtract",
-        "subtract mapped source from target"
-    );
     argList::addOption
     (
         "fields",
@@ -242,12 +126,6 @@ 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;
-    }
-
     HashSet<word> selectedFields;
     if (args.optionFound("fields"))
     {
@@ -316,7 +194,6 @@ int main(int argc, char *argv[])
             meshSource,
             meshTarget,
             mapMethod,
-            subtract,
             selectedFields,
             noLagrangian
         );
@@ -330,7 +207,6 @@ int main(int argc, char *argv[])
             patchMap,
             cuttingPatches,
             mapMethod,
-            subtract,
             selectedFields,
             noLagrangian
         );
diff --git a/applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H b/applications/utilities/preProcessing/mapFieldsPar/mapMeshes.C
similarity index 62%
rename from applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H
rename to applications/utilities/preProcessing/mapFieldsPar/mapMeshes.C
index f33fbe7c6d..7697c01940 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/MapMeshes.H
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapMeshes.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,21 +23,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef MapMeshes_H
-#define MapMeshes_H
-
-#include "MapVolFields.H"
-#include "mapLagrangian.H"
-#include "UnMapped.H"
+#include "mapMeshes.H"
+#include "surfaceMesh.H"
 #include "pointMesh.H"
+#include "mapLagrangian.H"
+#include "MapVolFields.H"
+#include "UnMapped.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-template<template<class> class CombineOp>
-void MapMesh
+void Foam::mapMesh
 (
     const meshToMesh& interp,
     const HashSet<word>& selectedFields,
@@ -56,37 +51,32 @@ void MapMesh
         (
             objects,
             selectedFields,
-            interp,
-            CombineOp<scalar>()
+            interp
         );
 
         MapVolFields<vector>
         (
             objects,
             selectedFields,
-            interp,
-            CombineOp<vector>()
+            interp
         );
         MapVolFields<sphericalTensor>
         (
             objects,
             selectedFields,
-            interp,
-            CombineOp<sphericalTensor>()
+            interp
         );
         MapVolFields<symmTensor>
         (
             objects,
             selectedFields,
-            interp,
-            CombineOp<symmTensor>()
+            interp
         );
         MapVolFields<tensor>
         (
             objects,
             selectedFields,
-            interp,
-            CombineOp<tensor>()
+            interp
         );
     }
 
@@ -120,12 +110,87 @@ void MapMesh
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::mapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const word& mapMethod,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+)
+{
+    Info<< nl << "Consistently creating and mapping fields for time "
+        << meshSource.time().timeName() << nl << endl;
 
-} // End namespace Foam
+    meshToMesh interp(meshSource, meshTarget, mapMethod);
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    mapMesh
+    (
+        interp,
+        selectedFields,
+        noLagrangian
+    );
+}
+
+
+void Foam::mapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const word& mapMethod,
+    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
+    );
+
+    mapMesh
+    (
+        interp,
+        selectedFields,
+        noLagrangian
+    );
+}
+
+
+Foam::wordList Foam::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();
+}
 
-#endif
 
 // ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapMeshes.H b/applications/utilities/preProcessing/mapFieldsPar/mapMeshes.H
new file mode 100644
index 0000000000..edbf88474d
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapMeshes.H
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 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 "meshToMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void mapMesh
+(
+    const meshToMesh& interp,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+);
+
+
+void mapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const word& mapMethod,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+);
+
+
+void mapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const word& mapMethod,
+    const HashSet<word>& selectedFields,
+    const bool noLagrangian
+);
+
+
+wordList addProcessorPatches
+(
+    const fvMesh& meshTarget,
+    const wordList& cuttingPatches
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C b/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
index 0e594fb090..45d40dd7a3 100644
--- a/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
+++ b/src/fvModels/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,12 +81,7 @@ Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource
 
     const fvMesh& nbrMesh = mesh.time().lookupObject<fvMesh>(nbrRegionName());
 
-    meshInterp().mapTgtToSrc
-    (
-        scalarField(nbrMesh.nCells(), 1),
-        plusEqOp<scalar>(),
-        filter_
-    );
+    meshInterp().mapTgtToSrc(scalarField(nbrMesh.nCells(), 1), filter_);
 
     const word zoneName(name + ":porous");
 
diff --git a/src/fvModels/interRegion/interRegionModel/interRegionModelTemplates.C b/src/fvModels/interRegion/interRegionModel/interRegionModelTemplates.C
index f07f1b7370..b9b804fd2b 100644
--- a/src/fvModels/interRegion/interRegionModel/interRegionModelTemplates.C
+++ b/src/fvModels/interRegion/interRegionModel/interRegionModelTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,11 +67,11 @@ void Foam::fv::interRegionModel::interpolate
 {
     if (master())
     {
-        meshInterp().mapTgtToSrc(field, plusEqOp<scalar>(), result);
+        meshInterp().mapTgtToSrc(field, result);
     }
     else
     {
-        nbrModel.meshInterp().mapSrcToTgt(field, plusEqOp<scalar>(), result);
+        nbrModel.meshInterp().mapSrcToTgt(field, result);
     }
 }
 
diff --git a/src/sampling/meshToMesh/meshToMesh.C b/src/sampling/meshToMesh/meshToMesh.C
index 6e8821291f..335e238754 100644
--- a/src/sampling/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMesh/meshToMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,8 +43,7 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<scalar>& srcField,
-    Field<scalar>& tgtField,
-    const plusEqOp<scalar>& cop
+    Field<scalar>& tgtField
 ) const
 {}
 
@@ -54,8 +53,7 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<vector>& srcField,
-    Field<vector>& tgtField,
-    const plusEqOp<vector>& cop
+    Field<vector>& tgtField
 ) const
 {}
 
@@ -65,8 +63,7 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<sphericalTensor>& srcField,
-    Field<sphericalTensor>& tgtField,
-    const plusEqOp<sphericalTensor>& cop
+    Field<sphericalTensor>& tgtField
 ) const
 {}
 
@@ -76,8 +73,7 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<symmTensor>& srcField,
-    Field<symmTensor>& tgtField,
-    const plusEqOp<symmTensor>& cop
+    Field<symmTensor>& tgtField
 ) const
 {}
 
@@ -87,8 +83,7 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<tensor>& srcField,
-    Field<tensor>& tgtField,
-    const plusEqOp<tensor>& cop
+    Field<tensor>& tgtField
 ) const
 {}
 
@@ -98,8 +93,7 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<scalar>& srcField,
-    const Field<scalar>& tgtField,
-    const plusEqOp<scalar>& cop
+    const Field<scalar>& tgtField
 ) const
 {}
 
@@ -109,8 +103,7 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<vector>& srcField,
-    const Field<vector>& tgtField,
-    const plusEqOp<vector>& cop
+    const Field<vector>& tgtField
 ) const
 {}
 
@@ -120,8 +113,7 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<sphericalTensor>& srcField,
-    const Field<sphericalTensor>& tgtField,
-    const plusEqOp<sphericalTensor>& cop
+    const Field<sphericalTensor>& tgtField
 ) const
 {}
 
@@ -131,8 +123,7 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<symmTensor>& srcField,
-    const Field<symmTensor>& tgtField,
-    const plusEqOp<symmTensor>& cop
+    const Field<symmTensor>& tgtField
 ) const
 {}
 
@@ -142,8 +133,7 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<tensor>& srcField,
-    const Field<tensor>& tgtField,
-    const plusEqOp<tensor>& cop
+    const Field<tensor>& tgtField
 ) const
 {}
 
diff --git a/src/sampling/meshToMesh/meshToMesh.H b/src/sampling/meshToMesh/meshToMesh.H
index 0aada6e4ac..1ebdd4e26c 100644
--- a/src/sampling/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMesh/meshToMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -114,24 +114,22 @@ class meshToMesh
 
         //- Helper function to interpolate patch field. Template
         //  specialisations  below
-        template<class Type, class CombineOp>
+        template<class Type>
         void mapAndOpSrcToTgt
         (
             const AMIInterpolation& AMI,
             const Field<Type>& srcField,
-            Field<Type>& tgtField,
-            const CombineOp& cop
+            Field<Type>& tgtField
         ) const;
 
         //- Helper function to interpolate patch field. Template
         //  specialisations  below
-        template<class Type, class CombineOp>
+        template<class Type>
         void mapAndOpTgtToSrc
         (
             const AMIInterpolation& AMI,
             Field<Type>& srcField,
-            const Field<Type>& tgtField,
-            const CombineOp& cop
+            const Field<Type>& tgtField
         ) const;
 
         //- Return src cell IDs for the overlap region
@@ -302,34 +300,15 @@ public:
                 //- Map field from src to tgt mesh with defined operation
                 //  Values passed in via 'result' are used to initialise the
                 //  return value
-                template<class Type, class CombineOp>
+                template<class Type>
                 void mapSrcToTgt
                 (
                     const UList<Type>& srcFld,
-                    const CombineOp& cop,
                     List<Type>& result
                 ) const;
 
                 //- Return the src field mapped to the tgt mesh with a defined
                 //  operation.  Initial values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<Field<Type>> mapSrcToTgt
-                (
-                    const Field<Type>& srcFld,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a tmp field to the tgt mesh
-                //  with a defined operation
-                template<class Type, class CombineOp>
-                tmp<Field<Type>> mapSrcToTgt
-                (
-                    const tmp<Field<Type>>& tsrcFld,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a field to the tgt mesh with a
-                //  default operation (plusEqOp)
                 template<class Type>
                 tmp<Field<Type>> mapSrcToTgt
                 (
@@ -337,7 +316,7 @@ public:
                 ) const;
 
                 //- Convenience function to map a tmp field to the tgt mesh
-                //  with a default operation (plusEqOp)
+                //  with a defined operation
                 template<class Type>
                 tmp<Field<Type>> mapSrcToTgt
                 (
@@ -350,34 +329,15 @@ public:
                 //- Map field from tgt to src mesh with defined operation
                 //  Values passed in via 'result' are used to initialise the
                 //  return value
-                template<class Type, class CombineOp>
+                template<class Type>
                 void mapTgtToSrc
                 (
                     const UList<Type>& tgtFld,
-                    const CombineOp& cop,
                     List<Type>& result
                 ) const;
 
                 //- Return the tgt field mapped to the src mesh with a defined
                 //  operation.  Initial values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<Field<Type>> mapTgtToSrc
-                (
-                    const Field<Type>& tgtFld,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a tmp field to the src mesh
-                //  with a defined operation
-                template<class Type, class CombineOp>
-                tmp<Field<Type>> mapTgtToSrc
-                (
-                    const tmp<Field<Type>>& ttgtFld,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a field to the src mesh with a
-                //  default operation (plusEqOp)
                 template<class Type>
                 tmp<Field<Type>> mapTgtToSrc
                 (
@@ -385,7 +345,7 @@ public:
                 ) const;
 
                 //- Convenience function to map a tmp field to the src mesh
-                //  with a default operation (plusEqOp)
+                //  with a defined operation
                 template<class Type>
                 tmp<Field<Type>> mapTgtToSrc
                 (
@@ -398,43 +358,23 @@ public:
                 //- Interpolate a field with a defined operation.  Values
                 //  passed in via 'result' are used to initialise the return
                 //  value
-                template<class Type, class CombineOp>
+                template<class Type>
                 void mapSrcToTgt
                 (
                     const GeometricField<Type, fvPatchField, volMesh>& field,
-                    const CombineOp& cop,
                     GeometricField<Type, fvPatchField, volMesh>& result
                 ) const;
 
                 //- Interpolate a field with a defined operation.  The initial
                 //  values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<GeometricField<Type, fvPatchField, volMesh>> mapSrcToTgt
-                (
-                    const GeometricField<Type, fvPatchField, volMesh>& field,
-                    const CombineOp& cop
-                ) const;
-
-                //- Interpolate a tmp field with a defined operation.  The
-                //  initial values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<GeometricField<Type, fvPatchField, volMesh>> mapSrcToTgt
-                (
-                    const tmp<GeometricField<Type, fvPatchField, volMesh>>&
-                        tfield,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a field with a default
-                //  operation (plusEqOp)
                 template<class Type>
                 tmp<GeometricField<Type, fvPatchField, volMesh>> mapSrcToTgt
                 (
                     const GeometricField<Type, fvPatchField, volMesh>& field
                 ) const;
 
-                //- Convenience function to map a tmp field with a default
-                //  operation (plusEqOp)
+                //- Interpolate a tmp field with a defined operation.  The
+                //  initial values of the result are set to zero
                 template<class Type>
                 tmp<GeometricField<Type, fvPatchField, volMesh>> mapSrcToTgt
                 (
@@ -448,43 +388,23 @@ public:
                 //- Interpolate a field with a defined operation.  Values
                 //  passed in via 'result' are used to initialise the return
                 //  value
-                template<class Type, class CombineOp>
+                template<class Type>
                 void mapTgtToSrc
                 (
                     const GeometricField<Type, fvPatchField, volMesh>& field,
-                    const CombineOp& cop,
                     GeometricField<Type, fvPatchField, volMesh>& result
                 ) const;
 
                 //- Interpolate a field with a defined operation.  The initial
                 //  values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<GeometricField<Type, fvPatchField, volMesh>> mapTgtToSrc
-                (
-                    const GeometricField<Type, fvPatchField, volMesh>& field,
-                    const CombineOp& cop
-                ) const;
-
-                //- Interpolate a tmp field with a defined operation.  The
-                //  initial values of the result are set to zero
-                template<class Type, class CombineOp>
-                tmp<GeometricField<Type, fvPatchField, volMesh>> mapTgtToSrc
-                (
-                    const tmp<GeometricField<Type, fvPatchField, volMesh>>&
-                        tfield,
-                    const CombineOp& cop
-                ) const;
-
-                //- Convenience function to map a field with a default
-                //  operation (plusEqOp)
                 template<class Type>
                 tmp<GeometricField<Type, fvPatchField, volMesh>> mapTgtToSrc
                 (
                     const GeometricField<Type, fvPatchField, volMesh>& field
                 ) const;
 
-                //- Convenience function to map a tmp field with a default
-                //  operation (plusEqOp)
+                //- Interpolate a tmp field with a defined operation.  The
+                //  initial values of the result are set to zero
                 template<class Type>
                 tmp<GeometricField<Type, fvPatchField, volMesh>> mapTgtToSrc
                 (
@@ -509,40 +429,35 @@ void meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<scalar>& srcField,
-    Field<scalar>& tgtField,
-    const plusEqOp<scalar>& cop
+    Field<scalar>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<vector>& srcField,
-    Field<vector>& tgtField,
-    const plusEqOp<vector>& cop
+    Field<vector>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<sphericalTensor>& srcField,
-    Field<sphericalTensor>& tgtField,
-    const plusEqOp<sphericalTensor>& cop
+    Field<sphericalTensor>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<symmTensor>& srcField,
-    Field<symmTensor>& tgtField,
-    const plusEqOp<symmTensor>& cop
+    Field<symmTensor>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<tensor>& srcField,
-    Field<tensor>& tgtField,
-    const plusEqOp<tensor>& cop
+    Field<tensor>& tgtField
 ) const;
 
 
@@ -551,40 +466,35 @@ void meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<scalar>& srcField,
-    const Field<scalar>& tgtField,
-    const plusEqOp<scalar>& cop
+    const Field<scalar>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<vector>& srcField,
-    const Field<vector>& tgtField,
-    const plusEqOp<vector>& cop
+    const Field<vector>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<sphericalTensor>& srcField,
-    const Field<sphericalTensor>& tgtField,
-    const plusEqOp<sphericalTensor>& cop
+    const Field<sphericalTensor>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<symmTensor>& srcField,
-    const Field<symmTensor>& tgtField,
-    const plusEqOp<symmTensor>& cop
+    const Field<symmTensor>& tgtField
 ) const;
 template<>
 void meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<tensor>& srcField,
-    const Field<tensor>& tgtField,
-    const plusEqOp<tensor>& cop
+    const Field<tensor>& tgtField
 ) const;
 
 } // End namespace Foam
diff --git a/src/sampling/meshToMesh/meshToMeshTemplates.C b/src/sampling/meshToMesh/meshToMeshTemplates.C
index 33828e400d..3d30a96ae1 100644
--- a/src/sampling/meshToMesh/meshToMeshTemplates.C
+++ b/src/sampling/meshToMesh/meshToMeshTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,11 +77,10 @@ void Foam::meshToMesh::add
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapSrcToTgt
 (
     const UList<Type>& srcField,
-    const CombineOp& cop,
     List<Type>& result
 ) const
 {
@@ -95,8 +94,6 @@ void Foam::meshToMesh::mapSrcToTgt
             << abort(FatalError);
     }
 
-    multiplyWeightedOp<Type, CombineOp> cbop(cop);
-
     if (singleMeshProc_ == -1)
     {
         const mapDistribute& map = srcMapPtr_();
@@ -111,13 +108,10 @@ void Foam::meshToMesh::mapSrcToTgt
 
             if (srcAddress.size())
             {
-//                result[celli] = Zero;
                 result[celli] *= (1.0 - sum(srcWeight));
                 forAll(srcAddress, i)
                 {
-                    label srcI = srcAddress[i];
-                    scalar w = srcWeight[i];
-                    cbop(result[celli], celli, work[srcI], w);
+                    result[celli] += srcWeight[i]*work[srcAddress[i]];
                 }
             }
         }
@@ -131,13 +125,10 @@ void Foam::meshToMesh::mapSrcToTgt
 
             if (srcAddress.size())
             {
-//                result[celli] = Zero;
                 result[celli] *= (1.0 - sum(srcWeight));
                 forAll(srcAddress, i)
                 {
-                    label srcI = srcAddress[i];
-                    scalar w = srcWeight[i];
-                    cbop(result[celli], celli, srcField[srcI], w);
+                    result[celli] += srcWeight[i]*srcField[srcAddress[i]];
                 }
             }
         }
@@ -145,11 +136,10 @@ void Foam::meshToMesh::mapSrcToTgt
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
 (
-    const Field<Type>& srcField,
-    const CombineOp& cop
+    const Field<Type>& srcField
 ) const
 {
     tmp<Field<Type>> tresult
@@ -161,33 +151,12 @@ Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
         )
     );
 
-    mapSrcToTgt(srcField, cop, tresult.ref());
+    mapSrcToTgt(srcField, tresult.ref());
 
     return tresult;
 }
 
 
-template<class Type, class CombineOp>
-Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
-(
-    const tmp<Field<Type>>& tsrcField,
-    const CombineOp& cop
-) const
-{
-    return mapSrcToTgt(tsrcField(), cop);
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
-(
-    const Field<Type>& srcField
-) const
-{
-    return mapSrcToTgt(srcField, plusEqOp<Type>());
-}
-
-
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
 (
@@ -198,11 +167,10 @@ Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapTgtToSrc
 (
     const UList<Type>& tgtField,
-    const CombineOp& cop,
     List<Type>& result
 ) const
 {
@@ -216,8 +184,6 @@ void Foam::meshToMesh::mapTgtToSrc
             << abort(FatalError);
     }
 
-    multiplyWeightedOp<Type, CombineOp> cbop(cop);
-
     if (singleMeshProc_ == -1)
     {
         const mapDistribute& map = tgtMapPtr_();
@@ -235,9 +201,7 @@ void Foam::meshToMesh::mapTgtToSrc
                 result[celli] *= (1.0 - sum(tgtWeight));
                 forAll(tgtAddress, i)
                 {
-                    label tgtI = tgtAddress[i];
-                    scalar w = tgtWeight[i];
-                    cbop(result[celli], celli, work[tgtI], w);
+                    result[celli] += tgtWeight[i]*work[tgtAddress[i]];
                 }
             }
         }
@@ -254,9 +218,7 @@ void Foam::meshToMesh::mapTgtToSrc
                 result[celli] *= (1.0 - sum(tgtWeight));
                 forAll(tgtAddress, i)
                 {
-                    label tgtI = tgtAddress[i];
-                    scalar w = tgtWeight[i];
-                    cbop(result[celli], celli, tgtField[tgtI], w);
+                    result[celli] += tgtWeight[i]*tgtField[tgtAddress[i]];
                 }
             }
         }
@@ -264,11 +226,10 @@ void Foam::meshToMesh::mapTgtToSrc
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
 (
-    const Field<Type>& tgtField,
-    const CombineOp& cop
+    const Field<Type>& tgtField
 ) const
 {
     tmp<Field<Type>> tresult
@@ -280,50 +241,28 @@ Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
         )
     );
 
-    mapTgtToSrc(tgtField, cop, tresult.ref());
+    mapTgtToSrc(tgtField, tresult.ref());
 
     return tresult;
 }
 
 
-template<class Type, class CombineOp>
-Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
-(
-    const tmp<Field<Type>>& ttgtField,
-    const CombineOp& cop
-) const
-{
-    return mapTgtToSrc(ttgtField(), cop);
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
-(
-    const Field<Type>& tgtField
-) const
-{
-    return mapTgtToSrc(tgtField, plusEqOp<Type>());
-}
-
-
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
 (
     const tmp<Field<Type>>& ttgtField
 ) const
 {
-    return mapTgtToSrc(ttgtField(), plusEqOp<Type>());
+    return mapTgtToSrc(ttgtField());
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapAndOpSrcToTgt
 (
     const AMIInterpolation& AMI,
     const Field<Type>& srcField,
-    Field<Type>& tgtField,
-    const CombineOp& cop
+    Field<Type>& tgtField
 ) const
 {
     tgtField = pTraits<Type>::zero;
@@ -331,22 +270,20 @@ void Foam::meshToMesh::mapAndOpSrcToTgt
     AMI.interpolateToTarget
     (
         srcField,
-        multiplyWeightedOp<Type, CombineOp>(cop),
         tgtField,
         UList<Type>::null()
     );
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapSrcToTgt
 (
     const GeometricField<Type, fvPatchField, volMesh>& field,
-    const CombineOp& cop,
     GeometricField<Type, fvPatchField, volMesh>& result
 ) const
 {
-    mapSrcToTgt(field, cop, result.primitiveFieldRef());
+    mapSrcToTgt(field, result.primitiveFieldRef());
 
     const PtrList<AMIInterpolation>& AMIList = patchAMIs();
 
@@ -389,7 +326,7 @@ void Foam::meshToMesh::mapSrcToTgt
 
         // Override value to account for CombineOp (note: is dummy template
         // specialisation for plusEqOp)
-        mapAndOpSrcToTgt(AMIList[i], srcField, tgtField, cop);
+        mapAndOpSrcToTgt(AMIList[i], srcField, tgtField);
     }
 
     forAll(cuttingPatches_, i)
@@ -401,12 +338,11 @@ void Foam::meshToMesh::mapSrcToTgt
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::meshToMesh::mapSrcToTgt
 (
-    const GeometricField<Type, fvPatchField, volMesh>& field,
-    const CombineOp& cop
+    const GeometricField<Type, fvPatchField, volMesh>& field
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
@@ -485,35 +421,12 @@ Foam::meshToMesh::mapSrcToTgt
         )
     );
 
-    mapSrcToTgt(field, cop, tresult.ref());
+    mapSrcToTgt(field, tresult.ref());
 
     return tresult;
 }
 
 
-template<class Type, class CombineOp>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::meshToMesh::mapSrcToTgt
-(
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfield,
-    const CombineOp& cop
-) const
-{
-    return mapSrcToTgt(tfield(), cop);
-}
-
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::meshToMesh::mapSrcToTgt
-(
-    const GeometricField<Type, fvPatchField, volMesh>& field
-) const
-{
-    return mapSrcToTgt(field, plusEqOp<Type>());
-}
-
-
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::meshToMesh::mapSrcToTgt
@@ -521,17 +434,16 @@ Foam::meshToMesh::mapSrcToTgt
     const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfield
 ) const
 {
-    return mapSrcToTgt(tfield(), plusEqOp<Type>());
+    return mapSrcToTgt(tfield());
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapAndOpTgtToSrc
 (
     const AMIInterpolation& AMI,
     Field<Type>& srcField,
-    const Field<Type>& tgtField,
-    const CombineOp& cop
+    const Field<Type>& tgtField
 ) const
 {
     srcField = pTraits<Type>::zero;
@@ -539,22 +451,20 @@ void Foam::meshToMesh::mapAndOpTgtToSrc
     AMI.interpolateToSource
     (
         tgtField,
-        multiplyWeightedOp<Type, CombineOp>(cop),
         srcField,
         UList<Type>::null()
     );
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh::mapTgtToSrc
 (
     const GeometricField<Type, fvPatchField, volMesh>& field,
-    const CombineOp& cop,
     GeometricField<Type, fvPatchField, volMesh>& result
 ) const
 {
-    mapTgtToSrc(field, cop, result.primitiveFieldRef());
+    mapTgtToSrc(field, result.primitiveFieldRef());
 
     const PtrList<AMIInterpolation>& AMIList = patchAMIs();
 
@@ -596,7 +506,7 @@ void Foam::meshToMesh::mapTgtToSrc
 
         // Override value to account for CombineOp (could be dummy for
         // plusEqOp)
-        mapAndOpTgtToSrc(AMIList[i], srcField, tgtField, cop);
+        mapAndOpTgtToSrc(AMIList[i], srcField, tgtField);
     }
 
     forAll(cuttingPatches_, i)
@@ -608,12 +518,11 @@ void Foam::meshToMesh::mapTgtToSrc
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::meshToMesh::mapTgtToSrc
 (
-    const GeometricField<Type, fvPatchField, volMesh>& field,
-    const CombineOp& cop
+    const GeometricField<Type, fvPatchField, volMesh>& field
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
@@ -692,35 +601,12 @@ Foam::meshToMesh::mapTgtToSrc
         )
     );
 
-    mapTgtToSrc(field, cop, tresult.ref());
+    mapTgtToSrc(field, tresult.ref());
 
     return tresult;
 }
 
 
-template<class Type, class CombineOp>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::meshToMesh::mapTgtToSrc
-(
-    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfield,
-    const CombineOp& cop
-) const
-{
-    return mapTgtToSrc(tfield(), cop);
-}
-
-
-template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
-Foam::meshToMesh::mapTgtToSrc
-(
-    const GeometricField<Type, fvPatchField, volMesh>& field
-) const
-{
-    return mapTgtToSrc(field, plusEqOp<Type>());
-}
-
-
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::meshToMesh::mapTgtToSrc
@@ -728,7 +614,7 @@ Foam::meshToMesh::mapTgtToSrc
     const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfield
 ) const
 {
-    return mapTgtToSrc(tfield(), plusEqOp<Type>());
+    return mapTgtToSrc(tfield());
 }
 
 
diff --git a/src/sampling/meshToMesh0/meshToMesh0.H b/src/sampling/meshToMesh0/meshToMesh0.H
index 23e573e614..fa985563a5 100644
--- a/src/sampling/meshToMesh0/meshToMesh0.H
+++ b/src/sampling/meshToMesh0/meshToMesh0.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -254,105 +254,95 @@ public:
         // Interpolation
 
             //- Map field
-            template<class Type, class CombineOp>
+            template<class Type>
             void mapField
             (
                 Field<Type>&,
                 const Field<Type>&,
-                const labelList& adr,
-                const CombineOp& cop
+                const labelList& adr
             ) const;
 
             //- Interpolate field using inverse-distance weights
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const scalarListList& weights,
-                const CombineOp& cop
+                const scalarListList& weights
             ) const;
 
             //- Interpolate field using inverse-volume weights
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelListList& adr,
-                const scalarListList& weights,
-                const CombineOp& cop
+                const scalarListList& weights
             ) const;
 
 
             //- Interpolate field using cell-point interpolation
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const vectorField& centres,
-                const CombineOp& cop
+                const vectorField& centres
             )const;
 
 
             //- Interpolate internal volume field
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh>>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 
-            template<class Type, class CombineOp>
+            template<class Type>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh>>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type, class CombineOp>
+            template<class Type>
             tmp<GeometricField<Type, fvPatchField, volMesh>> interpolate
             (
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 
-            template<class Type, class CombineOp>
+            template<class Type>
             tmp<GeometricField<Type, fvPatchField, volMesh>> interpolate
             (
                 const tmp<GeometricField<Type, fvPatchField, volMesh>>&,
-                order=INTERPOLATE,
-                const CombineOp& cop = eqOp<Type>()
+                order=INTERPOLATE
             ) const;
 };
 
diff --git a/src/sampling/meshToMesh0/meshToMesh0Templates.C b/src/sampling/meshToMesh0/meshToMesh0Templates.C
index f77904e8bd..0fa0fd3454 100644
--- a/src/sampling/meshToMesh0/meshToMesh0Templates.C
+++ b/src/sampling/meshToMesh0/meshToMesh0Templates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,13 +31,12 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::mapField
 (
     Field<Type>& toF,
     const Field<Type>& fromVf,
-    const labelList& adr,
-    const CombineOp& cop
+    const labelList& adr
 ) const
 {
     // Direct mapping of nearest-cell values
@@ -46,7 +45,7 @@ void Foam::meshToMesh0::mapField
     {
         if (adr[celli] != -1)
         {
-            cop(toF[celli], fromVf[adr[celli]]);
+            toF[celli] = fromVf[adr[celli]];
         }
     }
 
@@ -54,14 +53,13 @@ void Foam::meshToMesh0::mapField
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelListList& adr,
-    const scalarListList& weights,
-    const CombineOp& cop
+    const scalarListList& weights
 ) const
 {
     // Inverse volume weighted interpolation
@@ -75,20 +73,19 @@ void Foam::meshToMesh0::interpolateField
         {
             label fromCelli = overlapCells[i];
             f += fromVf[fromCelli]*w[i];
-            cop(toF[celli], f);
+            toF[celli] = f;
         }
     }
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const scalarListList& weights,
-    const CombineOp& cop
+    const scalarListList& weights
 ) const
 {
     // Inverse distance weighted interpolation
@@ -110,20 +107,19 @@ void Foam::meshToMesh0::interpolateField
                 f += fromVf[neighbours[ni - 1]]*w[ni];
             }
 
-            cop(toF[celli], f);
+            toF[celli] = f;
         }
     }
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const vectorField& centres,
-    const CombineOp& cop
+    const vectorField& centres
 ) const
 {
     // Cell-Point interpolation
@@ -133,27 +129,18 @@ void Foam::meshToMesh0::interpolateField
     {
         if (adr[celli] != -1)
         {
-            cop
-            (
-                toF[celli],
-                interpolator.interpolate
-                (
-                    centres[celli],
-                    adr[celli]
-                )
-            );
+            toF[celli] = interpolator.interpolate(centres[celli], adr[celli]);
         }
     }
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolateInternalField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh0::order ord,
-    const CombineOp& cop
+    meshToMesh0::order ord
 ) const
 {
     if (fromVf.mesh() != fromMesh_)
@@ -177,7 +164,7 @@ void Foam::meshToMesh0::interpolateInternalField
     switch(ord)
     {
         case MAP:
-            mapField(toF, fromVf, cellAddressing_, cop);
+            mapField(toF, fromVf, cellAddressing_);
         break;
 
         case INTERPOLATE:
@@ -187,8 +174,7 @@ void Foam::meshToMesh0::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                inverseDistanceWeights(),
-                cop
+                inverseDistanceWeights()
             );
             break;
         }
@@ -199,8 +185,7 @@ void Foam::meshToMesh0::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                toMesh_.cellCentres(),
-                cop
+                toMesh_.cellCentres()
             );
 
             break;
@@ -215,8 +200,7 @@ void Foam::meshToMesh0::interpolateInternalField
                 toF,
                 fromVf,
                 cellToCell,
-                invVolWeights,
-                cop
+                invVolWeights
             );
             break;
         }
@@ -228,30 +212,28 @@ void Foam::meshToMesh0::interpolateInternalField
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolateInternalField
 (
     Field<Type>& toF,
     const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfromVf,
-    meshToMesh0::order ord,
-    const CombineOp& cop
+    meshToMesh0::order ord
 ) const
 {
-    interpolateInternalField(toF, tfromVf(), ord, cop);
+    interpolateInternalField(toF, tfromVf(), ord);
     tfromVf.clear();
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh0::order ord,
-    const CombineOp& cop
+    meshToMesh0::order ord
 ) const
 {
-    interpolateInternalField(toVf, fromVf, ord, cop);
+    interpolateInternalField(toVf, fromVf, ord);
 
     typename GeometricField<Type, fvPatchField, volMesh>::
         Boundary& toVfBf = toVf.boundaryFieldRef();
@@ -270,8 +252,7 @@ void Foam::meshToMesh0::interpolate
                     (
                         toVfBf[patchi],
                         fromVf,
-                        boundaryAddressing_[patchi],
-                        cop
+                        boundaryAddressing_[patchi]
                     );
                     break;
                 }
@@ -283,8 +264,7 @@ void Foam::meshToMesh0::interpolate
                         toVfBf[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf(),
-                        cop
+                        toPatch.Cf()
                     );
                     break;
                 }
@@ -296,8 +276,7 @@ void Foam::meshToMesh0::interpolate
                         toVfBf[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf(),
-                        cop
+                        toPatch.Cf()
                     );
                     break;
                 }
@@ -333,40 +312,37 @@ void Foam::meshToMesh0::interpolate
                 [
                     fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
                 ],
-                boundaryAddressing_[patchi],
-                cop
+                boundaryAddressing_[patchi]
             );
         }
     }
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 void Foam::meshToMesh0::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const tmp<GeometricField<Type, fvPatchField, volMesh>>& tfromVf,
-    meshToMesh0::order ord,
-    const CombineOp& cop
+    meshToMesh0::order ord
 ) const
 {
-    interpolate(toVf, tfromVf(), ord, cop);
+    interpolate(toVf, tfromVf(), ord);
     tfromVf.clear();
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::meshToMesh0::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh0::order ord,
-    const CombineOp& cop
+    meshToMesh0::order ord
 ) const
 {
     // Create and map the internal-field values
     Field<Type> internalField(toMesh_.nCells());
-    interpolateInternalField(internalField, fromVf, ord, cop);
+    interpolateInternalField(internalField, fromVf, ord);
 
     // check whether both meshes have got the same number
     // of boundary patches
@@ -428,17 +404,16 @@ Foam::meshToMesh0::interpolate
 }
 
 
-template<class Type, class CombineOp>
+template<class Type>
 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
+    meshToMesh0::order ord
 ) const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh>> tint =
-        interpolate(tfromVf(), ord, cop);
+        interpolate(tfromVf(), ord);
     tfromVf.clear();
 
     return tint;