diff --git a/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
index e3d255d961..a0ebeb01c5 100644
--- a/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
+++ b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
@@ -116,12 +116,16 @@ Description
     {
         const faPatch& p = patches[patchi];
 
-         Info<< "  " << "patch " << p.index()
-             << " (size: " << returnReduce(p.size(), sumOp<label>())
-             << ") name: " << p.name()
-             << nl;
+        // Report physical size (nEdges) not virtual size
+        Info<< "  " << "patch " << p.index()
+            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
+            << ") name: " << p.name()
+            << nl;
     }
 
+    Info<< "----------------" << nl
+        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
+
 
     // Geometry information
     Info<< nl;
diff --git a/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H b/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
index 5dd75c8804..084759849f 100644
--- a/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
+++ b/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -36,12 +36,17 @@ Description
     {
         const faPatch& p = patches[patchi];
 
-         Info<< "  " << "patch " << p.index()
-             << " (size: " << returnReduce(p.size(), sumOp<label>())
-             << ") name: " << p.name()
-             << nl;
+        // Report physical size (nEdges) not virtual size
+        Info<< "  " << "patch " << p.index()
+            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
+            << ") name: " << p.name()
+            << nl;
     }
 
+    Info<< "----------------" << nl
+        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
+
+
     // Geometry information
     Info<< nl;
     {
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
index 1098eb95fa..df80a98e6d 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@@ -609,32 +609,30 @@ Foam::label Foam::faBoundaryMesh::findPatchID
 
 Foam::label Foam::faBoundaryMesh::whichPatch(const label edgeIndex) const
 {
-    // Find out which patch the current face belongs to by comparing label
-    // with patch start labels.
-    // If the face is internal, return -1;
-    // if it is off the end of the list, abort
     if (edgeIndex < mesh().nInternalEdges())
     {
+        // Internal edge
         return -1;
     }
     else if (edgeIndex >= mesh().nEdges())
     {
+        // Bounds error: abort
         FatalErrorInFunction
             << "Edge " << edgeIndex
             << " out of bounds. Number of geometric edges " << mesh().nEdges()
             << abort(FatalError);
+
+        return -1;
     }
 
+    // Find patch that the edgeIndex belongs to.
+
     forAll(*this, patchi)
     {
         label start = mesh_.patchStarts()[patchi];
         label size = operator[](patchi).faPatch::size();
 
-        if
-        (
-            edgeIndex >= start
-         && edgeIndex < start + size
-        )
+        if (edgeIndex >= start && edgeIndex < start + size)
         {
             return patchi;
         }
@@ -642,7 +640,7 @@ Foam::label Foam::faBoundaryMesh::whichPatch(const label edgeIndex) const
 
     // If not in any of above, it's trouble!
     FatalErrorInFunction
-        << "error in patch search algorithm"
+        << "Error in patch search algorithm"
         << abort(FatalError);
 
     return -1;
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
index 989e1b9e04..e3945375ea 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
@@ -42,11 +42,12 @@ Author
 #ifndef Foam_faBoundaryMesh_H
 #define Foam_faBoundaryMesh_H
 
+#include "regIOobject.H"
 #include "faPatch.H"
+#include "labelPair.H"
 #include "lduInterfacePtrsList.H"
 #include "wordList.H"
 #include "pointField.H"
-#include "regIOobject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C
index 34b43b2380..67f18252ff 100644
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@@ -221,6 +221,8 @@ void Foam::faMesh::clearGeomNotAreas() const
 
     clearHalo();
     patchPtr_.reset(nullptr);
+    polyPatchFacesPtr_.reset(nullptr);
+    polyPatchIdsPtr_.reset(nullptr);
     bndConnectPtr_.reset(nullptr);
     deleteDemandDrivenData(SPtr_);
     deleteDemandDrivenData(patchStartsPtr_);
@@ -360,6 +362,8 @@ Foam::faMesh::faMesh
     curTimeIndex_(time().timeIndex()),
 
     patchPtr_(nullptr),
+    polyPatchFacesPtr_(nullptr),
+    polyPatchIdsPtr_(nullptr),
     bndConnectPtr_(nullptr),
     lduPtr_(nullptr),
 
@@ -464,6 +468,8 @@ Foam::faMesh::faMesh
     curTimeIndex_(time().timeIndex()),
 
     patchPtr_(nullptr),
+    polyPatchFacesPtr_(nullptr),
+    polyPatchIdsPtr_(nullptr),
     bndConnectPtr_(nullptr),
     lduPtr_(nullptr),
 
@@ -543,6 +549,8 @@ Foam::faMesh::faMesh
     curTimeIndex_(time().timeIndex()),
 
     patchPtr_(nullptr),
+    polyPatchFacesPtr_(nullptr),
+    polyPatchIdsPtr_(nullptr),
     bndConnectPtr_(nullptr),
     lduPtr_(nullptr),
 
@@ -706,6 +714,22 @@ const Foam::word& Foam::faMesh::regionName() const
 }
 
 
+Foam::labelList Foam::faMesh::faceCells() const
+{
+    const labelList& faceOwner = this->mesh().faceOwner();
+
+    labelList list(faceLabels_);
+
+    for (label& val : list)
+    {
+        // Transcribe from faceId to cellId (owner)
+        val = faceOwner[val];
+    }
+
+    return list;
+}
+
+
 void Foam::faMesh::removeFiles(const fileName& instanceDir) const
 {
     fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();
diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H
index 62734dfe18..7a4a17c8e3 100644
--- a/src/finiteArea/faMesh/faMesh.H
+++ b/src/finiteArea/faMesh/faMesh.H
@@ -253,6 +253,12 @@ class faMesh
         //- Primitive patch
         mutable std::unique_ptr<uindirectPrimitivePatch> patchPtr_;
 
+        //- List of poly patch/patch-face for faceLabels
+        mutable std::unique_ptr<List<labelPair>> polyPatchFacesPtr_;
+
+        //- List of polyPatch ids used by the area mesh
+        mutable std::unique_ptr<labelList> polyPatchIdsPtr_;
+
         //- List of proc/mesh-face for boundary edge neighbours
         mutable std::unique_ptr<List<labelPair>> bndConnectPtr_;
 
@@ -368,6 +374,10 @@ class faMesh
             //- Calculate patch starts in the edge list
             void calcPatchStarts() const;
 
+            //- Calculate which polyPatches, polyPatch/local-face
+            //- are related to the areaMesh
+            void calcWhichPatchFaces() const;
+
             //- Calculate edge lengths
             //  Triggers communication via calcEdgeAreaNormals
             void calcLe() const;
@@ -702,6 +712,14 @@ public:
             //- Return the underlying polyMesh face labels
             inline const labelList& faceLabels() const noexcept;
 
+            //- The area-face corresponding to the mesh-face, -1 if not found
+            inline label whichFace(const label meshFacei) const;
+
+            //- The polyPatches related to the areaMesh, in sorted order
+            inline const labelList& whichPolyPatches() const;
+
+            //- The polyPatch/local-face for each faceLabels()
+            inline const List<labelPair>& whichPatchFaces() const;
 
             //- Return parallel info
             const faGlobalMeshData& globalData() const;
@@ -762,6 +780,12 @@ public:
             tmp<vectorField> haloFaceNormals(const label patchi) const;
 
 
+    // Interfacing
+
+        //- The volume (owner) cells associated with the area-mesh
+        labelList faceCells() const;
+
+
     // Storage management
 
         //- Remove all files from mesh instance
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index c5df1546b9..7fd15f4db4 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -277,6 +277,46 @@ void Foam::faMesh::calcPatchStarts() const
 }
 
 
+void Foam::faMesh::calcWhichPatchFaces() const
+{
+    // Usually need both together
+    if (polyPatchFacesPtr_ || polyPatchIdsPtr_)
+    {
+        FatalErrorInFunction
+            << "Already allocated polyPatchFaces/polyPatchIds"
+            << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    polyPatchFacesPtr_.reset
+    (
+        new List<labelPair>(pbm.whichPatchFace(faceLabels_))
+    );
+
+    labelHashSet ids;
+
+    // Extract patch ids from (patch, facei) tuples
+    for (const labelPair& tup : *polyPatchFacesPtr_)
+    {
+        ids.insert(tup.first());
+    }
+
+    ids.erase(-1);  // Without internal faces (patchi == -1)
+
+    // parSync
+    Foam::reduce
+    (
+        ids,
+        bitOrOp<labelHashSet>(),
+        UPstream::msgType(),
+        this->comm()
+    );
+
+    polyPatchIdsPtr_.reset(new labelList(ids.sortedToc()));
+}
+
+
 void Foam::faMesh::calcLe() const
 {
     DebugInFunction
diff --git a/src/finiteArea/faMesh/faMeshI.H b/src/finiteArea/faMesh/faMeshI.H
index dfdbc00f79..3d3db7b958 100644
--- a/src/finiteArea/faMesh/faMeshI.H
+++ b/src/finiteArea/faMesh/faMeshI.H
@@ -145,6 +145,32 @@ inline Foam::uindirectPrimitivePatch& Foam::faMesh::patch()
 }
 
 
+inline Foam::label Foam::faMesh::whichFace(const label meshFacei) const
+{
+    return faceLabels_.find(meshFacei);
+}
+
+
+inline const Foam::labelList& Foam::faMesh::whichPolyPatches() const
+{
+    if (!polyPatchFacesPtr_ || !polyPatchIdsPtr_)
+    {
+        calcWhichPatchFaces();
+    }
+    return *polyPatchIdsPtr_;
+}
+
+
+inline const Foam::List<Foam::labelPair>& Foam::faMesh::whichPatchFaces() const
+{
+    if (!polyPatchFacesPtr_ || !polyPatchIdsPtr_)
+    {
+        calcWhichPatchFaces();
+    }
+    return *polyPatchFacesPtr_;
+}
+
+
 inline bool Foam::faMesh::hasInternalEdgeLabels() const noexcept
 {
     return false;
diff --git a/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.C b/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.C
index ee03325dc6..03bbb802f6 100644
--- a/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.C
+++ b/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,47 +24,472 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
 #include "volSurfaceMapping.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+void Foam::volSurfaceMapping::mapToSurface
 (
-    const GeometricBoundaryField<Type, fvPatchField, volMesh>& df
+    const GeometricBoundaryField<Type, fvPatchField, volMesh>& bfld,
+    Field<Type>& result
 ) const
 {
-    // Grab labels for all faces in faMesh
-    const labelList& faceLabels = mesh_.faceLabels();
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
 
-    auto tresult = tmp<Field<Type>>::New(faceLabels.size(), Zero);
-    auto& result = tresult.ref();
+    // FULLDEBUG: checkSize ?
+    // or simply result.resize(mesh_.nFaces());
 
-    // Get reference to volume mesh
-    const polyMesh& pMesh = mesh_();
-    const polyBoundaryMesh& bm = pMesh.boundaryMesh();
-
-    label patchID, faceID;
-
-    // Grab droplet cloud source by identifying patch and face
-    forAll(faceLabels, i)
+    forAll(patchFaces, i)
     {
-        // Escape if face is beyond active faces, eg belongs to a face zone
-        if (faceLabels[i] < pMesh.nFaces())
-        {
-            patchID = bm.whichPatch(faceLabels[i]);
-            faceID = bm[patchID].whichFace(faceLabels[i]);
+        const labelPair& patchAndFace = patchFaces[i];
 
-            result[i] = df[patchID][faceID];
+        if (patchAndFace.first() >= 0)
+        {
+            result[i] = bfld[patchAndFace];
         }
     }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const GeometricBoundaryField<Type, fvPatchField, volMesh>& bfld
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(mesh_.nFaces(), Zero);
+    mapToSurface(bfld, tresult.ref());
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vfld,
+    Field<Type>& result
+) const
+{
+    mapToSurface(vfld.boundaryField(), result);
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
+    Field<Type>& result
+) const
+{
+    mapToSurface(tvf().boundaryField(), result);
+    tvf.clear();
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vfld
+) const
+{
+    return mapToSurface(vfld.boundaryField());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+) const
+{
+    tmp<Field<Type>> tresult(mapToSurface(tvf().boundaryField()));
+    tvf.clear();
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& fld,
+    Field<Type>& result
+) const
+{
+    const auto& bfld = fld.boundaryField();
+
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
+
+    // FULLDEBUG: checkSize ?
+    // or simply result.resize(mesh_.nFaces());
+
+    forAll(patchFaces, i)
+    {
+        const labelPair& patchAndFace = patchFaces[i];
+
+        if (patchAndFace.first() >= 0)
+        {
+            // Value from boundary
+            result[i] = bfld[patchAndFace];
+        }
+        else if (patchAndFace.second() >= 0)
+        {
+            // Value from internal
+            result[i] = fld[patchAndFace.second()];
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& fld
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(mesh_.nFaces(), Zero);
+    mapToSurface(fld, tresult.ref());
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>& tsf,
+    Field<Type>& result
+) const
+{
+    mapToSurface(tsf(), result);
+    tsf.clear();
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>& tsf
+) const
+{
+    tmp<Field<Type>> tresult(mapToSurface(tsf()));
+    tsf.clear();
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const UPtrList<Field<Type>>& patchFields,
+    Field<Type>& result
+) const
+{
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
+
+    // FULLDEBUG: checkSize ?
+    // or simply result.resize(mesh_.nFaces());
+
+    forAll(patchFaces, i)
+    {
+        const label patchi = patchFaces[i].first();
+        const label facei = patchFaces[i].second();
+
+        const auto* pfld = patchFields.get(patchi);
+
+        if (pfld)
+        {
+            result[i] = (*pfld)[facei];
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToSurface
+(
+    const PtrMap<Field<Type>>& patchFields,
+    Field<Type>& result
+) const
+{
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
+
+    // FULLDEBUG: checkSize ?
+    // or simply result.resize(mesh_.nFaces());
+
+    forAll(patchFaces, i)
+    {
+        const label patchi = patchFaces[i].first();
+        const label facei = patchFaces[i].second();
+
+        const auto* pfld = patchFields.get(patchi);
+
+        if (pfld)
+        {
+            result[i] = (*pfld)[facei];
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const UPtrList<Field<Type>>& patchFields
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(mesh_.nFaces(), Zero);
+    mapToSurface(patchFields, tresult.ref());
+    return tresult;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
+(
+    const PtrMap<Field<Type>>& patchFields
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(mesh_.nFaces(), Zero);
+    mapToSurface(patchFields, tresult.ref());
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const GeometricBoundaryField<Type, fvPatchField, volMesh>& bfld,
+    Field<Type>& result
+) const
+{
+    PtrList<Field<Type>> patchFields;
+
+    // All referenced polyPatches (sorted order)
+    const labelList& patches = mesh_.whichPolyPatches();
+
+    if (!patches.empty())
+    {
+        // maxPolyPatch+1
+        patchFields.resize(patches.last()+1);
+    }
+
+    // Populate patchInternalField
+    for (const label patchi : patches)
+    {
+        patchFields.set(patchi, bfld[patchi].patchInternalField());
+    }
+
+    mapToSurface(patchFields, result);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const GeometricBoundaryField<Type, fvPatchField, volMesh>& bfld
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(mesh_.nFaces(), Zero);
+
+    mapInternalToSurface(bfld, tresult.ref());
 
     return tresult;
 }
 
 
+template<class Type>
+void Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    Field<Type>& result
+) const
+{
+    mapInternalToSurface(vf.boundaryField(), result);
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf,
+    Field<Type>& result
+) const
+{
+    mapInternalToSurface(tvf().boundaryField(), result);
+    tvf.clear();
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vfld
+) const
+{
+    return mapInternalToSurface(vfld.boundaryField());
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapInternalToSurface
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvf
+) const
+{
+    tmp<Field<Type>> tresult(mapInternalToSurface(tvf().boundaryField()));
+    tvf.clear();
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const DimensionedField<Type, areaMesh>& af,
+    GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
+
+    forAll(patchFaces, i)
+    {
+        const labelPair& patchAndFace = patchFaces[i];
+
+        if
+        (
+            patchAndFace.first() >= 0
+         && (patchAndFace.first() == destPatchi || destPatchi < 0)
+        )
+        {
+            dest[patchAndFace] = af[i];
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const tmp<DimensionedField<Type, areaMesh>>& taf,
+    GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolume(taf(), dest, destPatchi);
+    taf.clear();
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const DimensionedField<Type, areaMesh>& af,
+    GeometricField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolume(af, dest.boundaryFieldRef(), destPatchi);
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const tmp<DimensionedField<Type, areaMesh>>& taf,
+    GeometricField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolume(taf, dest.boundaryFieldRef(), destPatchi);
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>& taf,
+    GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolume(taf().internalField(), dest, destPatchi);
+    taf.clear();
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolume
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>& taf,
+    GeometricField<Type, fvPatchField, volMesh>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolume(taf().internalField(), dest.boundaryFieldRef(), destPatchi);
+    taf.clear();
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolumePatch
+(
+    const DimensionedField<Type, areaMesh>& af,
+    Field<Type>& dest,
+    const label destPatchi
+) const
+{
+    // The polyPatch/local-face for each of the faceLabels
+    const List<labelPair>& patchFaces = mesh_.whichPatchFaces();
+
+    forAll(patchFaces, i)
+    {
+        const label patchi = patchFaces[i].first();
+        const label facei = patchFaces[i].second();
+
+        if (patchi >= 0 && patchi == destPatchi)
+        {
+            dest[facei] = af[i];
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolumePatch
+(
+    const tmp<DimensionedField<Type, areaMesh>>& taf,
+    Field<Type>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolumePatch(taf(), dest, destPatchi);
+    taf.clear();
+}
+
+
+template<class Type>
+void Foam::volSurfaceMapping::mapToVolumePatch
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>& taf,
+    Field<Type>& dest,
+    const label destPatchi
+) const
+{
+    mapToVolumePatch(taf().internalField(), dest, destPatchi);
+    taf.clear();
+}
+
+
+// These are fragile/wrong. To be removed (2022-09)
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
 (
@@ -76,7 +501,7 @@ Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
     auto tresult = tmp<Field<Type>>::New(faceLabels.size(), Zero);
     auto& result = tresult.ref();
 
-    const polyMesh& pMesh = mesh_();
+    const polyMesh& pMesh = mesh_.mesh();
     const polyBoundaryMesh& bm = pMesh.boundaryMesh();
     label patchID, faceID;
 
@@ -95,86 +520,6 @@ Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapToSurface
 }
 
 
-template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::volSurfaceMapping::mapInternalToSurface
-(
-    const GeometricBoundaryField<Type, fvPatchField, volMesh>& df
-) const
-{
-    // Grab labels for all faces in faMesh
-    const labelList& faceLabels = mesh_.faceLabels();
-
-    auto tresult = tmp<Field<Type>>::New(faceLabels.size(), Zero);
-    auto& result = tresult.ref();
-
-    // Get reference to volume mesh
-    const polyMesh& pMesh = mesh_();
-    const polyBoundaryMesh& bm = pMesh.boundaryMesh();
-
-    label patchID, faceID;
-
-    // Grab droplet cloud source by identifying patch and face
-    forAll(faceLabels, i)
-    {
-        // Escape if face is beyond active faces, eg belongs to a face zone
-        if (faceLabels[i] < pMesh.nFaces())
-        {
-            patchID = bm.whichPatch(faceLabels[i]);
-            faceID = bm[patchID].whichFace(faceLabels[i]);
-
-            result[i] = df[patchID].patchInternalField()()[faceID];
-        }
-    }
-
-    return tresult;
-}
-
-
-template<class Type>
-void Foam::volSurfaceMapping::mapToVolume
-(
-    const GeometricField<Type, faPatchField, areaMesh>& af,
-    GeometricBoundaryField<Type, fvPatchField, volMesh>& bf
-) const
-{
-    // Grab labels for all faces in faMesh
-    const labelList& faceLabels = mesh_.faceLabels();
-
-    // Get reference to volume mesh
-    const polyMesh& pMesh = mesh_();
-    const polyBoundaryMesh& bm = pMesh.boundaryMesh();
-
-    label patchID, faceID;
-
-    const Field<Type>& afi = af.internalField();
-
-    forAll(faceLabels, i)
-    {
-        // Escape if face is beyond active faces, eg belongs to a face zone
-        if (faceLabels[i] < pMesh.nFaces())
-        {
-            patchID = bm.whichPatch(faceLabels[i]);
-            faceID = bm[patchID].whichFace(faceLabels[i]);
-
-            bf[patchID][faceID] = afi[i];
-        }
-    }
-}
-
-
-template<class Type>
-void Foam::volSurfaceMapping::mapToVolume
-(
-    const tmp<GeometricField<Type, faPatchField, areaMesh>>& taf,
-    GeometricBoundaryField<Type, fvPatchField, volMesh>& bf
-) const
-{
-    mapToVolume(taf(), bf);
-
-    taf.clear();
-}
-
-
 template<class Type>
 void Foam::volSurfaceMapping::mapToField
 (
@@ -182,9 +527,7 @@ void Foam::volSurfaceMapping::mapToField
     Field<Type>& f
 ) const
 {
-    const Field<Type>& afi = af.internalField();
-
-    mapToField(afi, f);
+    mapToField(af.primitiveField(), f);
 }
 
 
@@ -197,7 +540,7 @@ void Foam::volSurfaceMapping::mapToField
 {
     const labelList& faceLabels = mesh_.faceLabels();
 
-    const polyMesh& pMesh = mesh_();
+    const polyMesh& pMesh = mesh_.mesh();
     const polyBoundaryMesh& bm = pMesh.boundaryMesh();
     label patchID, faceID;
 
diff --git a/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.H b/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.H
index 3138249998..4bcd5c8177 100644
--- a/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.H
+++ b/src/finiteArea/interpolation/volSurfaceMapping/volSurfaceMapping.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,8 +41,11 @@ SourceFiles
 #ifndef Foam_volSurfaceMapping_H
 #define Foam_volSurfaceMapping_H
 
-#include "faMesh.H"
+#include "areaFields.H"
 #include "volMesh.H"
+#include "surfaceMesh.H"
+#include "PtrList.H"
+#include "PtrMap.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +54,7 @@ namespace Foam
 
 // Forward Declarations
 template<class Type> class fvPatchField;
+template<class Type> class fvsPatchField;
 
 /*---------------------------------------------------------------------------*\
                       Class volSurfaceMapping Declaration
@@ -69,7 +73,7 @@ public:
     // Constructors
 
         //- Construct from mesh
-        volSurfaceMapping(const faMesh& mesh)
+        explicit volSurfaceMapping(const faMesh& mesh)
         :
             mesh_(mesh)
         {}
@@ -85,59 +89,274 @@ public:
     ~volSurfaceMapping() = default;
 
 
-    // Member Functions
+    // Mapping to area fields
 
-        //- Map volume boundary field to surface
+        //- Map volume boundary fields as area field
+        template<class Type>
+        void mapToSurface
+        (
+            const GeometricBoundaryField<Type, fvPatchField, volMesh>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map volume boundary fields as area field
         template<class Type>
         tmp<Field<Type>> mapToSurface
         (
-            const GeometricBoundaryField<Type, fvPatchField, volMesh>& df
+            const GeometricBoundaryField<Type, fvPatchField, volMesh>&
         ) const;
 
-        //- Map vol Field to surface Field
+        //- Map volume (boundary) fields to area field
         template<class Type>
-        tmp<Field<Type>> mapToSurface(const Field<Type>& f) const;
+        void mapToSurface
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            Field<Type>& result
+        ) const;
 
-        //- Map patch internal field to surface
+        //- Map volume (boundary) fields to area field
+        template<class Type>
+        void mapToSurface
+        (
+            const tmp<GeometricField<Type, fvPatchField, volMesh>>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map volume (boundary) fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Map volume (boundary) fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const tmp<GeometricField<Type, fvPatchField, volMesh>>&
+        ) const;
+
+
+        //- Map surface fields to area field
+        template<class Type>
+        void mapToSurface
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map surface fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        ) const;
+
+        //- Map surface fields to area field
+        template<class Type>
+        void mapToSurface
+        (
+            const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map surface fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>&
+        ) const;
+
+
+        //- Map pre-calculated boundary fields to area field
+        template<class Type>
+        void mapToSurface
+        (
+            const UPtrList<Field<Type>>& patchFields,
+            Field<Type>& result
+        ) const;
+
+        //- Map pre-calculated boundary fields to area field
+        template<class Type>
+        void mapToSurface
+        (
+            const PtrMap<Field<Type>>& patchFields,
+            Field<Type>& result
+        ) const;
+
+        //- Map pre-calculated boundary fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const UPtrList<Field<Type>>& patchFields
+        ) const;
+
+        //- Map pre-calculated boundary fields to area field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface
+        (
+            const PtrMap<Field<Type>>& patchFields
+        ) const;
+
+
+        //- Map patch internal field to area field
+        template<class Type>
+        void mapInternalToSurface
+        (
+            const GeometricBoundaryField<Type, fvPatchField, volMesh>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map patch internal field to area field
         template<class Type>
         tmp<Field<Type>> mapInternalToSurface
         (
-            const GeometricBoundaryField<Type, fvPatchField, volMesh>& df
+            const GeometricBoundaryField<Type, fvPatchField, volMesh>&
         ) const;
 
-        //- Map surface field to volume boundary field
+        //- Map patch internal field to area field
+        template<class Type>
+        void mapInternalToSurface
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map patch internal field to area field
+        template<class Type>
+        void mapInternalToSurface
+        (
+            const tmp<GeometricField<Type, fvPatchField, volMesh>>&,
+            Field<Type>& result
+        ) const;
+
+        //- Map patch internal field to area field
+        template<class Type>
+        tmp<Field<Type>> mapInternalToSurface
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Map patch internal field to area field
+        template<class Type>
+        tmp<Field<Type>> mapInternalToSurface
+        (
+            const tmp<GeometricField<Type, fvPatchField, volMesh>>&
+        ) const;
+
+
+    // Mapping from area field to volume (boundary) field
+
+        //- Map area field to volume boundary field,
+        //- optionally restricted to a single destination patch
         template<class Type>
         void mapToVolume
         (
-            const GeometricField<Type, faPatchField, areaMesh>& af,
-            GeometricBoundaryField<Type, fvPatchField, volMesh>& bf
+            const DimensionedField<Type, areaMesh>&,
+            GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
         ) const;
 
-        //- Map surface tmp field to volume boundary field
+        //- Map area field to volume boundary field,
+        //- optionally restricted to a single destination patch
         template<class Type>
         void mapToVolume
+        (
+            const tmp<DimensionedField<Type, areaMesh>>&,
+            GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
+        ) const;
+
+        //- Map area field to volume boundary field,
+        //- optionally restricted to a single destination patch
+        template<class Type>
+        void mapToVolume
+        (
+            const DimensionedField<Type, areaMesh>&,
+            GeometricField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
+        ) const;
+
+        //- Map area field to volume boundary field,
+        //- optionally restricted to a single destination patch
+        template<class Type>
+        void mapToVolume
+        (
+            const tmp<DimensionedField<Type, areaMesh>>&,
+            GeometricField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
+        ) const;
+
+        //- Map tmp area field to volume boundary field,
+        //- optionally restricted to a single destination patch
+        template<class Type>
+        void mapToVolume
+        (
+            const tmp<GeometricField<Type, faPatchField, areaMesh>>&,
+            GeometricBoundaryField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
+        ) const;
+
+        //- Map tmp area field to volume boundary field,
+        //- optionally restricted to a single destination patch
+        template<class Type>
+        void mapToVolume
+        (
+            const tmp<GeometricField<Type, faPatchField, areaMesh>>&,
+            GeometricField<Type, fvPatchField, volMesh>& dest,
+            const label destPatchi = -1
+        ) const;
+
+
+    // Mapping from area field to volume patch field
+
+        //- Map area field to a volume boundary patch
+        template<class Type>
+        void mapToVolumePatch
+        (
+            const DimensionedField<Type, areaMesh>& af,
+            Field<Type>& dest,
+            const label destPatchi
+        ) const;
+
+        //- Map tmp area field to a volume boundary patch
+        template<class Type>
+        void mapToVolumePatch
+        (
+            const tmp<DimensionedField<Type, areaMesh>>& af,
+            Field<Type>& dest,
+            const label destPatchi
+        ) const;
+
+        //- Map tmp area field to a volume boundary patch
+        template<class Type>
+        void mapToVolumePatch
         (
             const tmp<GeometricField<Type, faPatchField, areaMesh>>& taf,
-            GeometricBoundaryField<Type, fvPatchField, volMesh>& bf
+            Field<Type>& dest,
+            const label destPatchi
         ) const;
 
+
+    // Deprecated (to be removed)
+
+        //- Map vol Field to surface Field
+        template<class Type>
+        tmp<Field<Type>> mapToSurface(const Field<Type>&) const;
+
         //- Map surface field to field
         //  Assumes Field faces in the same order as Boundary
         template<class Type>
         void mapToField
         (
             const GeometricField<Type, faPatchField, areaMesh>& af,
-            Field<Type>& f
+            Field<Type>&
         ) const;
 
         //- Map surface field to volume field
         //  Assumes Field faces in the same order as Boundary
         template<class Type>
-        void mapToField
-        (
-            const Field<Type>& af,
-            Field<Type>& f
-        ) const;
+        void mapToField(const Field<Type>& af, Field<Type>&) const;
 };