diff --git a/applications/utilities/parallelProcessing/redistributePar/Make/files b/applications/utilities/parallelProcessing/redistributePar/Make/files
index 2dacc462a991ed7201c109c5eebf849cf13672e0..73b91011bdaf5fda0426a8c4e682dd3ade93d456 100644
--- a/applications/utilities/parallelProcessing/redistributePar/Make/files
+++ b/applications/utilities/parallelProcessing/redistributePar/Make/files
@@ -1,7 +1,10 @@
 passivePositionParticleCloud.C
+
 parLagrangianDistributor.C
 parLagrangianDistributorFields.C
+
 parPointFieldDistributor.C
+parFaFieldDistributorCache.C
 parFvFieldDistributor.C
 parFvFieldDistributorFields.C
 
diff --git a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C
index 7d789e56250ce48d15e887886595644e9263aae6..626ac35b830d91e23d9ce9f3f61d1ecb0158bb35 100644
--- a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C
+++ b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C
@@ -27,6 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "loadOrCreateMesh.H"
+#include "faMesh.H"
 #include "Pstream.H"
 #include "OSspecific.H"
 
@@ -66,6 +67,26 @@ Foam::boolList Foam::haveMeshFile
 }
 
 
+void Foam::removeProcAddressing(const faMesh& mesh)
+{
+    IOobject ioAddr
+    (
+        "procAddressing",
+        mesh.facesInstance(),
+        faMesh::meshSubDir,
+        mesh.thisDb()
+    );
+
+    for (const auto prefix : {"boundary", "edge", "face", "point"})
+    {
+        ioAddr.rename(prefix + word("ProcAddressing"));
+
+        const fileName procFile(ioAddr.objectPath());
+        Foam::rm(procFile);
+    }
+}
+
+
 void Foam::removeProcAddressing(const polyMesh& mesh)
 {
     IOobject ioAddr
diff --git a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.H b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.H
index bd468348e38234701da7b4db945d0cb931dc2033..3481e15744cabb1011b3474717d82868b8492979 100644
--- a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.H
+++ b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.H
@@ -45,6 +45,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
+class faMesh;
+
 //- Check for availability of specified mesh file (default: "faces")
 boolList haveMeshFile
 (
@@ -55,6 +58,9 @@ boolList haveMeshFile
 );
 
 
+//- Remove procAddressing
+void removeProcAddressing(const faMesh& mesh);
+
 //- Remove procAddressing
 void removeProcAddressing(const polyMesh& mesh);
 
diff --git a/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.C b/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.C
new file mode 100644
index 0000000000000000000000000000000000000000..04d3ba8857ec59dec4944023ad3ec97413a73a90
--- /dev/null
+++ b/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.C
@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "parFaFieldDistributorCache.H"
+
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "fieldsDistributor.H"
+#include "faMeshDistributor.H"
+#include "faMeshSubset.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class GeoField>
+void Foam::parFaFieldDistributorCache::redistributeAndWrite
+(
+    const faMeshDistributor& distributor,
+    PtrList<GeoField>& fields,
+    const bool isWriteProc
+)
+{
+    for (GeoField& fld : fields)
+    {
+        tmp<GeoField> tfld = distributor.distributeField(fld);
+
+        if (isWriteProc)
+        {
+            tfld().write();
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::parFaFieldDistributorCache::read
+(
+    const Time& baseRunTime,
+    const fileName& proc0CaseName,
+    const bool decompose,  // i.e. read from undecomposed case
+
+    const boolList& areaMeshOnProc,
+    const fileName& areaMeshInstance,
+    faMesh& mesh
+)
+{
+    Time& runTime = const_cast<Time&>(mesh.time());
+    const bool oldProcCase = runTime.processorCase();
+
+    autoPtr<faMeshSubset> subsetterPtr;
+
+    // Missing an area mesh somewhere?
+    if (areaMeshOnProc.found(false))
+    {
+        // A zero-sized mesh with boundaries.
+        // This is used to create zero-sized fields.
+        subsetterPtr.reset(new faMeshSubset(mesh, zero{}));
+
+        // Deregister from polyMesh ...
+        auto& obr = const_cast<objectRegistry&>
+        (
+            subsetterPtr->subMesh().thisDb()
+        );
+
+        obr.checkOut(faMesh::typeName);
+        obr.checkOut("faBoundaryMesh");
+        obr.checkOut("faSchemes");
+        obr.checkOut("faSolution");
+    }
+
+
+    // Get original objects (before incrementing time!)
+    if (Pstream::master() && decompose)
+    {
+        runTime.caseName() = baseRunTime.caseName();
+        runTime.processorCase(false);
+    }
+    IOobjectList objects(mesh.mesh(), runTime.timeName());
+    if (Pstream::master() && decompose)
+    {
+        runTime.caseName() = proc0CaseName;
+        runTime.processorCase(oldProcCase);
+    }
+
+    Info<< "From time " << runTime.timeName()
+        << " mesh:" << mesh.mesh().objectRegistry::objectRelPath()
+        << " have objects:" << objects.names() << endl;
+
+    if (Pstream::master() && decompose)
+    {
+        runTime.caseName() = baseRunTime.caseName();
+        runTime.processorCase(false);
+    }
+
+    #undef  doFieldReading
+    #define doFieldReading(Storage)                                   \
+    fieldsDistributor::readFields                                     \
+    (                                                                 \
+        areaMeshOnProc, mesh, subsetterPtr, objects, Storage,         \
+        true  /* (deregister field) */                                \
+    );
+
+    // areaFields
+    doFieldReading(scalarAreaFields_);
+    doFieldReading(vectorAreaFields_);
+    doFieldReading(sphericalTensorAreaFields_);
+    doFieldReading(symmTensorAreaFields_);
+    doFieldReading(tensorAreaFields_);
+
+    // edgeFields
+    doFieldReading(scalarEdgeFields_);
+    doFieldReading(vectorEdgeFields_);
+    doFieldReading(tensorEdgeFields_);
+    doFieldReading(sphericalTensorEdgeFields_);
+    doFieldReading(symmTensorEdgeFields_);
+    #undef doFieldReading
+}
+
+
+void Foam::parFaFieldDistributorCache::redistributeAndWrite
+(
+    const faMeshDistributor& distributor,
+    const bool isWriteProc
+)
+{
+    redistributeAndWrite(distributor, scalarAreaFields_, isWriteProc);
+    redistributeAndWrite(distributor, vectorAreaFields_, isWriteProc);
+    redistributeAndWrite(distributor, sphericalTensorAreaFields_, isWriteProc);
+    redistributeAndWrite(distributor, symmTensorAreaFields_, isWriteProc);
+    redistributeAndWrite(distributor, tensorAreaFields_, isWriteProc);
+
+    redistributeAndWrite(distributor, scalarEdgeFields_, isWriteProc);
+    redistributeAndWrite(distributor, vectorEdgeFields_, isWriteProc);
+    redistributeAndWrite(distributor, sphericalTensorEdgeFields_, isWriteProc);
+    redistributeAndWrite(distributor, symmTensorEdgeFields_, isWriteProc);
+    redistributeAndWrite(distributor, tensorEdgeFields_, isWriteProc);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.H b/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.H
new file mode 100644
index 0000000000000000000000000000000000000000..936bf9eee5e0d915cd338c9937aa0b780d92fb06
--- /dev/null
+++ b/applications/utilities/parallelProcessing/redistributePar/parFaFieldDistributorCache.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::parFaFieldDistributorCache
+
+Description
+    Simple container to manage read/write, redistribute finiteArea fields.
+
+SourceFiles
+    parFaFieldDistributorCache.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_parFaFieldDistributorCache_H
+#define Foam_parFaFieldDistributorCache_H
+
+#include "faMesh.H"
+#include "faMeshDistributor.H"
+#include "areaFieldsFwd.H"
+#include "edgeFieldsFwd.H"
+#include "PtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class parFaFieldDistributorCache Declaration
+\*---------------------------------------------------------------------------*/
+
+class parFaFieldDistributorCache
+{
+    // Private Data
+
+    #undef  declareField
+    #define declareField(Type)                                                \
+    PtrList<GeometricField<Type, faPatchField, areaMesh>> Type##AreaFields_;  \
+    PtrList<GeometricField<Type, faePatchField, edgeMesh>> Type##EdgeFields_;
+
+    declareField(scalar);
+    declareField(vector);
+    declareField(sphericalTensor);
+    declareField(symmTensor);
+    declareField(tensor);
+    #undef declareField
+
+
+    // Private Member Functions
+
+        //- Redistribute and write fields of given type (file-scope use)
+        template<class GeoField>
+        static void redistributeAndWrite
+        (
+            const faMeshDistributor& distributor,
+            PtrList<GeoField>& fields,
+            const bool isWriteProc
+        );
+
+
+public:
+
+    // Constructors
+
+        //- No copy construct
+        parFaFieldDistributorCache(const parFaFieldDistributorCache&) = delete;
+
+        //- No copy assignment
+        void operator=(const parFaFieldDistributorCache&) = delete;
+
+        //- Default construct
+        parFaFieldDistributorCache() = default;
+
+
+    // Member Functions
+
+        //- Read distributed fields
+        void read
+        (
+            const Time& baseRunTime,
+            const fileName& proc0CaseName,
+            const bool decompose,  // i.e. read from undecomposed case
+
+            const boolList& areaMeshOnProc,
+            const fileName& areaMeshInstance,
+            faMesh& mesh
+        );
+
+        void redistributeAndWrite
+        (
+            const faMeshDistributor& distributor,
+            const bool isWriteProc
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index f4d5a4cda4382584d5f6c0e46e1c0221a4e90a5e..d4f864e340dda60a666df1ab4728413846663cb5 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -100,7 +100,11 @@ Usage
 #include "meshRefinement.H"
 #include "pointFields.H"
 
-#include "readDistributedFields.H"
+#include "faMeshSubset.H"
+#include "faMeshTools.H"
+#include "faMeshDistributor.H"
+#include "parFaFieldDistributorCache.H"
+
 #include "redistributeLagrangian.H"
 
 #include "cyclicACMIFvPatch.H"
@@ -453,14 +457,12 @@ void determineDecomposition
     if (!decomposer.parallelAware())
     {
         WarningInFunction
-            << "You have selected decomposition method "
-            << decomposer.typeName
-            << " which does" << nl
-            << "not synchronise the decomposition across"
-            << " processor patches." << nl
-            << "    You might want to select a decomposition method"
-            << " which is aware of this. Continuing."
-            << endl;
+            << "You have selected decomposition method \""
+            << decomposer.type() << "\n"
+            << "    which does not synchronise decomposition across"
+               " processor patches.\n"
+               "    You might want to select a decomposition method"
+               " that is aware of this. Continuing...." << endl;
     }
 
     Time& tm = const_cast<Time&>(mesh.time());
@@ -552,16 +554,21 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
 (
     autoPtr<fileOperation>&& writeHandler,
     const Time& baseRunTime,
-    const boolList& haveMesh,
-    const fileName& meshSubDir,
+    const fileName& proc0CaseName,
+
+    // Controls
     const bool doReadFields,
     const bool decompose,       // decompose, i.e. read from undecomposed case
     const bool reconstruct,
     const bool overwrite,
-    const fileName& proc0CaseName,
+
+    // Decomposition information
     const label nDestProcs,
     const labelList& decomp,
-    const fileName& masterInstDir,
+
+    // Mesh information
+    const boolList& volMeshOnProc,
+    const fileName& volMeshInstance,
     fvMesh& mesh
 )
 {
@@ -621,24 +628,10 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         // processors
         autoPtr<fvMeshSubset> subsetterPtr;
 
-        const bool allHaveMesh = !haveMesh.found(false);
-        if (!allHaveMesh)
+        // Missing a volume mesh somewhere?
+        if (volMeshOnProc.found(false))
         {
-            // Find last non-processor patch.
-            const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-            const label nonProcI = (patches.nNonProcessor() - 1);
-
-            if (nonProcI < 0)
-            {
-                FatalErrorInFunction
-                    << "Cannot find non-processor patch on processor "
-                    << Pstream::myProcNo() << nl
-                    << " Current patches:" << patches.names()
-                    << abort(FatalError);
-            }
-
-            // Subset 0 cells, no parallel comms.
+            // A zero-sized mesh with boundaries.
             // This is used to create zero-sized fields.
             subsetterPtr.reset(new fvMeshSubset(mesh, zero{}));
         }
@@ -658,7 +651,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         }
 
         Info<< "From time " << runTime.timeName()
-            << " mesh:" << mesh.objectRegistry::objectPath()
+            << " mesh:" << mesh.objectRegistry::objectRelPath()
             << " have objects:" << objects.names() << endl;
 
         // We don't want to map the decomposition (mapping already tested when
@@ -683,7 +676,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         {                                                                     \
             fieldsDistributor::readFields                                     \
             (                                                                 \
-                haveMesh, mesh, subsetterPtr, objects, Storage                \
+                volMeshOnProc, mesh, subsetterPtr, objects, Storage           \
             );                                                                \
         }
 
@@ -716,7 +709,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         {                                                                     \
             fieldsDistributor::readFields                                     \
             (                                                                 \
-                haveMesh, oldPointMesh, subsetterPtr, objects, Storage,       \
+                volMeshOnProc, oldPointMesh, subsetterPtr, objects, Storage,  \
                 true  /* (deregister field) */                                \
             );                                                                \
             nPointFields += Storage.size();                                   \
@@ -802,7 +795,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
     }
     else
     {
-        mesh.setInstance(masterInstDir);
+        mesh.setInstance(volMeshInstance);
     }
 
 
@@ -827,7 +820,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
         if (Pstream::master())
         {
             Info<< "Setting caseName to " << baseRunTime.caseName()
-                << " to write reconstructed mesh and fields." << endl;
+                << " to write reconstructed mesh (and fields)." << endl;
             runTime.caseName() = baseRunTime.caseName();
             const bool oldProcCase(runTime.processorCase(false));
 
@@ -1050,7 +1043,12 @@ int main(int argc, char *argv[])
     (
         "Additional verbosity. (Can be used multiple times)"
     );
-
+    argList::addBoolOption
+    (
+        "no-finite-area",
+        "Suppress finiteArea mesh/field handling",
+        true  // Advanced option
+    );
 
     // Handle arguments
     // ~~~~~~~~~~~~~~~~
@@ -1104,6 +1102,7 @@ int main(int argc, char *argv[])
     const bool newTimes = args.found("newTimes");
     const int optVerbose = args.verbose();
 
+    const bool doFiniteArea = !args.found("no-finite-area");
     bool decompose = args.found("decompose");
     bool overwrite = args.found("overwrite");
 
@@ -1318,10 +1317,14 @@ int main(int argc, char *argv[])
             (
                 regionName == polyMesh::defaultRegion ? word::null : regionName
             );
-            const fileName meshSubDir(regionDir/polyMesh::meshSubDir);
 
-            Info<< nl << nl
-                << "Reconstructing mesh " << regionName << nl << endl;
+            const fileName volMeshSubDir(regionDir/polyMesh::meshSubDir);
+            const fileName areaMeshSubDir(regionDir/faMesh::meshSubDir);
+
+            Info<< nl
+                << "Reconstructing mesh " << regionDir << nl << endl;
+
+            bool areaMeshDetected = false;
 
             // Loop over all times
             forAll(timeDirs, timeI)
@@ -1332,40 +1335,71 @@ int main(int argc, char *argv[])
 
                 Info<< "Time = " << runTime.timeName() << endl << endl;
 
+                // Where meshes are
+                fileName volMeshInstance;
+                fileName areaMeshInstance;
 
-                // See where the mesh is
-                fileName facesInstance = runTime.findInstance
+                volMeshInstance = runTime.findInstance
                 (
-                    meshSubDir,
+                    volMeshSubDir,
                     "faces",
                     IOobject::READ_IF_PRESENT
                 );
-                //Pout<< "facesInstance:" << facesInstance << endl;
 
-                Pstream::scatter(facesInstance);
+                if (doFiniteArea)
+                {
+                    areaMeshInstance = runTime.findInstance
+                    (
+                        areaMeshSubDir,
+                        "faceLabels",
+                        IOobject::READ_IF_PRESENT
+                    );
+                }
 
-                // Check who has a mesh (by checking for 'faces' file)
-                const boolList haveMesh
+                Pstream::broadcasts
                 (
-                    haveMeshFile
+                    UPstream::worldComm,
+                    volMeshInstance,
+                    areaMeshInstance
+                );
+
+
+                // Check processors have meshes
+                // - check for 'faces' file (polyMesh)
+                // - check for 'faceLabels' file (faMesh)
+                boolList volMeshOnProc;
+                boolList areaMeshOnProc;
+
+                volMeshOnProc = haveMeshFile
+                (
+                    runTime,
+                    volMeshInstance/volMeshSubDir,
+                    "faces"
+                );
+
+                if (doFiniteArea)
+                {
+                    areaMeshOnProc = haveMeshFile
                     (
                         runTime,
-                        facesInstance/meshSubDir,
-                        "faces"
-                    )
-                );
+                        areaMeshInstance/areaMeshSubDir,
+                        "faceLabels"
+                    );
+
+                    areaMeshDetected = areaMeshOnProc.found(true);
+                }
 
 
                 // Addressing back to reconstructed mesh as xxxProcAddressing.
                 // - all processors have consistent faceProcAddressing
-                // - processors with no mesh don't need faceProcAddressing
+                // - processors without a mesh don't need faceProcAddressing
 
 
                 // Note: filePath searches up on processors that don't have
                 //       processor if instance = constant so explicitly check
                 //       found filename.
-                bool haveAddressing = false;
-                if (haveMesh[Pstream::myProcNo()])
+                bool haveVolAddressing = false;
+                if (volMeshOnProc[Pstream::myProcNo()])
                 {
                     // Read faces (just to know their size)
                     faceCompactIOList faces
@@ -1373,8 +1407,8 @@ int main(int argc, char *argv[])
                         IOobject
                         (
                             "faces",
-                            facesInstance,
-                            meshSubDir,
+                            volMeshInstance,
+                            volMeshSubDir,
                             runTime,
                             IOobject::MUST_READ
                         )
@@ -1386,108 +1420,312 @@ int main(int argc, char *argv[])
                         IOobject
                         (
                             "faceProcAddressing",
-                            facesInstance,
-                            meshSubDir,
+                            volMeshInstance,
+                            volMeshSubDir,
                             runTime,
                             IOobject::READ_IF_PRESENT
-                        ),
-                        labelList()
+                        )
                     );
-                    if
+
+                    haveVolAddressing =
                     (
                         faceProcAddressing.headerOk()
                      && faceProcAddressing.size() == faces.size()
-                    )
-                    {
-                        haveAddressing = true;
-                    }
+                    );
                 }
                 else
                 {
                     // Have no mesh. Don't need addressing
-                    haveAddressing = true;
+                    haveVolAddressing = true;
+                }
+
+                bool haveAreaAddressing = false;
+                if (areaMeshOnProc[Pstream::myProcNo()])
+                {
+                    // Read faces (just to know their size)
+                    labelIOList faceLabels
+                    (
+                        IOobject
+                        (
+                            "faceLabels",
+                            areaMeshInstance,
+                            areaMeshSubDir,
+                            runTime,
+                            IOobject::MUST_READ
+                        )
+                    );
+
+                    // Check faceProcAddressing
+                    labelIOList faceProcAddressing
+                    (
+                        IOobject
+                        (
+                            "faceProcAddressing",
+                            areaMeshInstance,
+                            areaMeshSubDir,
+                            runTime,
+                            IOobject::READ_IF_PRESENT
+                        )
+                    );
+
+                    haveAreaAddressing =
+                    (
+                        faceProcAddressing.headerOk()
+                     && faceProcAddressing.size() == faceLabels.size()
+                    );
+                }
+                else if (areaMeshDetected)
+                {
+                    // Have no mesh. Don't need addressing
+                    haveAreaAddressing = true;
                 }
 
 
                 // Additionally check for master faces being readable. Could
                 // do even more checks, e.g. global number of cells same
                 // as cellProcAddressing
-                bool haveUndecomposedMesh = false;
+
+                bool volMeshHaveUndecomposed = false;
+                bool areaMeshHaveUndecomposed = false;
+
                 if (Pstream::master())
                 {
                     Info<< "Checking " << baseRunTime.caseName()
-                        << " for undecomposed mesh" << endl;
+                        << " for undecomposed volume and area meshes..."
+                        << endl;
 
                     const bool oldParRun = Pstream::parRun(false);
-                    faceCompactIOList facesIO
-                    (
-                        IOobject
+
+                    // Volume
+                    {
+                        faceCompactIOList facesIO
                         (
-                            "faces",
-                            facesInstance,
-                            meshSubDir,
-                            baseRunTime,
-                            IOobject::NO_READ
-                        ),
-                        label(0)
-                    );
-                    haveUndecomposedMesh = facesIO.headerOk();
-                    Pstream::parRun(oldParRun);
+                            IOobject
+                            (
+                                "faces",
+                                volMeshInstance,
+                                volMeshSubDir,
+                                baseRunTime,
+                                IOobject::NO_READ
+                            ),
+                            label(0)
+                        );
+                        volMeshHaveUndecomposed = facesIO.headerOk();
+                    }
+
+                    // Area
+                    if (doFiniteArea)
+                    {
+                        labelIOList labelsIO
+                        (
+                            IOobject
+                            (
+                                "faceLabels",
+                                areaMeshInstance,
+                                areaMeshSubDir,
+                                baseRunTime,
+                                IOobject::NO_READ
+                            )
+                        );
+                        areaMeshHaveUndecomposed = labelsIO.headerOk();
+                    }
+
+                    Pstream::parRun(oldParRun);  // Restore parallel state
+                }
+
+                Pstream::broadcasts
+                (
+                    UPstream::worldComm,
+                    volMeshHaveUndecomposed,
+                    areaMeshHaveUndecomposed
+                );
+
+                // Report
+                {
+                    Info<< "    volume mesh ["
+                        << volMeshHaveUndecomposed << "] : "
+                        << volMeshInstance << nl
+                        << "    area   mesh ["
+                        << areaMeshHaveUndecomposed << "] : "
+                        << areaMeshInstance << nl
+                        << endl;
                 }
-                Pstream::scatter(haveUndecomposedMesh);
 
 
                 if
                 (
-                    !haveUndecomposedMesh
-                 || !returnReduce(haveAddressing, andOp<bool>())
+                    !volMeshHaveUndecomposed
+                 || !returnReduce(haveVolAddressing, andOp<bool>())
                 )
                 {
-                    Info<< "loading mesh from " << facesInstance << endl;
-                    autoPtr<fvMesh> meshPtr = fvMeshTools::loadOrCreateMesh
+                    Info<< "No undecomposed mesh. Creating from: "
+                        << volMeshInstance << endl;
+
+                    if (areaMeshHaveUndecomposed)
+                    {
+                        areaMeshHaveUndecomposed = false;
+                        Info<< "Also ignore any undecomposed area mesh"
+                            << endl;
+                    }
+
+                    autoPtr<fvMesh> volMeshPtr = fvMeshTools::loadOrCreateMesh
                     (
                         IOobject
                         (
                             regionName,
-                            facesInstance,
+                            volMeshInstance,
                             runTime,
                             Foam::IOobject::MUST_READ
                         ),
                         decompose
                     );
-                    fvMesh& mesh = meshPtr();
-
-                    // Use basic geometry calculation to avoid synchronisation
-                    // problems. See comment in routine
-                    fvMeshTools::setBasicGeometry(mesh);
-
+                    fvMeshTools::setBasicGeometry(volMeshPtr());
+                    fvMesh& mesh = volMeshPtr();
 
-                    // Determine decomposition
-                    // ~~~~~~~~~~~~~~~~~~~~~~~
 
-                    Info<< "Reconstructing mesh for time " << facesInstance
-                        << endl;
+                    Info<< nl << "Reconstructing mesh" << nl << endl;
 
-                    label nDestProcs = 1;
-                    labelList finalDecomp = labelList(mesh.nCells(), Zero);
+                    // Reconstruct (1 processor)
+                    const label nDestProcs(1);
+                    const labelList finalDecomp(mesh.nCells(), Zero);
 
                     redistributeAndWrite
                     (
                         std::move(writeHandler),
                         baseRunTime,
-                        haveMesh,
-                        meshSubDir,
+                        proc0CaseName,
+
+                        // Controls
                         false,      // do not read fields
                         false,      // do not read undecomposed case on proc0
                         true,       // write redistributed files to proc0
                         overwrite,
-                        proc0CaseName,
+
+                        // Decomposition information
                         nDestProcs,
                         finalDecomp,
-                        facesInstance,
+
+                        // For finite-volume
+                        volMeshOnProc,
+                        volMeshInstance,
                         mesh
                     );
                 }
+
+
+                // Similarly for finiteArea
+                // - may or may not have undecomposed mesh
+                // - may or may not have decomposed meshes
+
+                if
+                (
+                    areaMeshOnProc.found(true)  // ie, areaMeshDetected
+                 &&
+                    (
+                        !areaMeshHaveUndecomposed
+                     || !returnReduce(haveAreaAddressing, andOp<bool>())
+                    )
+                )
+                {
+                    Info<< "Loading area mesh from "
+                        << areaMeshInstance << endl;
+
+                    Info<< "    getting volume mesh support" << endl;
+
+                    autoPtr<fvMesh> baseMeshPtr = fvMeshTools::newMesh
+                    (
+                        IOobject
+                        (
+                            regionName,
+                            baseRunTime.timeName(),
+                            baseRunTime,
+                            IOobject::MUST_READ
+                        ),
+                        true            // read on master only
+                    );
+                    fvMeshTools::setBasicGeometry(baseMeshPtr());
+
+                    autoPtr<fvMesh> volMeshPtr = fvMeshTools::loadOrCreateMesh
+                    (
+                        IOobject
+                        (
+                            regionName,
+                            baseMeshPtr().facesInstance(),
+                            runTime,
+                            Foam::IOobject::MUST_READ
+                        ),
+                        decompose
+                    );
+                    fvMeshTools::setBasicGeometry(volMeshPtr());
+                    fvMesh& mesh = volMeshPtr();
+
+                    // Read volume proc addressing back to base mesh
+                    autoPtr<mapDistributePolyMesh> distMap
+                    (
+                        fvMeshTools::readProcAddressing(mesh, baseMeshPtr)
+                    );
+
+
+                    autoPtr<faMesh> areaMeshPtr = faMeshTools::loadOrCreateMesh
+                    (
+                        IOobject
+                        (
+                            regionName,
+                            areaMeshInstance,
+                            runTime,
+                            Foam::IOobject::MUST_READ
+                        ),
+                        mesh,  // <- The referenced polyMesh (from above)
+                        decompose
+                    );
+                    faMesh& areaMesh = areaMeshPtr();
+
+                    faMeshTools::forceDemandDriven(areaMesh);
+                    faMeshTools::unregisterMesh(areaMesh);
+
+                    autoPtr<faMesh> areaBaseMeshPtr;
+
+                    // Reconstruct using polyMesh distribute map
+                    mapDistributePolyMesh faDistMap
+                    (
+                        faMeshDistributor::distribute
+                        (
+                            areaMesh,
+                            distMap(),      // The polyMesh distMap
+                            baseMeshPtr(),  // Target polyMesh
+                            areaBaseMeshPtr
+                        )
+                    );
+
+                    faMeshTools::forceDemandDriven(areaBaseMeshPtr());
+                    faMeshTools::unregisterMesh(areaBaseMeshPtr());
+
+
+                    if (Pstream::master())
+                    {
+                        Info<< "Setting caseName to " << baseRunTime.caseName()
+                            << " to write reconstructed area mesh." << endl;
+                        runTime.caseName() = baseRunTime.caseName();
+                        const bool oldProcCase(runTime.processorCase(false));
+
+                        areaBaseMeshPtr().write();
+
+                        // Now we've written all. Reset caseName on master
+                        Info<< "Restoring caseName" << endl;
+                        runTime.caseName() = proc0CaseName;
+                        runTime.processorCase(oldProcCase);
+                    }
+
+                    // Update for the reconstructed procAddressing
+                    faMeshTools::writeProcAddressing
+                    (
+                        areaBaseMeshPtr(),  // Reconstruct location
+                        faDistMap,
+                        false,              // decompose=false
+                        std::move(writeHandler),
+                        areaMeshPtr.get()   // procMesh
+                    );
+                }
             }
 
             // Make sure all is finished writing until re-reading in pass2
@@ -1511,6 +1749,7 @@ int main(int argc, char *argv[])
             // This is a bit of tricky code and hidden inside fvMeshTools for
             // now.
             Info<< "Reading undecomposed mesh (on master)" << endl;
+
             autoPtr<fvMesh> baseMeshPtr = fvMeshTools::newMesh
             (
                 IOobject
@@ -1525,9 +1764,8 @@ int main(int argc, char *argv[])
 
             fvMeshTools::setBasicGeometry(baseMeshPtr());
 
-
             Info<< "Reading local, decomposed mesh" << endl;
-            autoPtr<fvMesh> meshPtr = fvMeshTools::loadOrCreateMesh
+            autoPtr<fvMesh> volMeshPtr = fvMeshTools::loadOrCreateMesh
             (
                 IOobject
                 (
@@ -1538,7 +1776,64 @@ int main(int argc, char *argv[])
                 ),
                 decompose
             );
-            fvMesh& mesh = meshPtr();
+            fvMesh& mesh = volMeshPtr();
+
+
+            // Similarly for finiteArea
+            autoPtr<faMesh> areaBaseMeshPtr;
+            autoPtr<faMesh> areaMeshPtr;
+            autoPtr<faMeshDistributor> faDistributor;
+            mapDistributePolyMesh areaDistMap;
+
+            if (areaMeshDetected)
+            {
+                areaBaseMeshPtr = faMeshTools::newMesh
+                (
+                    IOobject
+                    (
+                        regionName,
+                        baseRunTime.timeName(),
+                        baseRunTime,
+                        IOobject::MUST_READ
+                    ),
+                    baseMeshPtr(),
+                    true            // read on master only
+                );
+
+                areaMeshPtr = faMeshTools::loadOrCreateMesh
+                (
+                    IOobject
+                    (
+                        regionName,
+                        areaBaseMeshPtr().facesInstance(),
+                        runTime,
+                        IOobject::MUST_READ
+                    ),
+                    mesh,
+                    decompose
+                );
+
+                areaDistMap =
+                    faMeshTools::readProcAddressing
+                    (
+                        areaMeshPtr(),
+                        areaBaseMeshPtr
+                    );
+
+                faMeshTools::forceDemandDriven(areaMeshPtr());
+
+                // Create an appropriate field distributor
+                faDistributor.reset
+                (
+                    new faMeshDistributor
+                    (
+                        areaMeshPtr(),      // source
+                        areaBaseMeshPtr(),  // target
+                        areaDistMap
+                    )
+                );
+            }
+
 
             if (writeHandler && Pstream::master())
             {
@@ -1599,7 +1894,7 @@ int main(int argc, char *argv[])
                 if (newTimes && masterTimeDirSet.found(timeDirs[timeI].name()))
                 {
                     Info<< "Skipping time " << timeDirs[timeI].name()
-                        << endl << endl;
+                        << nl << endl;
                     continue;
                 }
 
@@ -1647,7 +1942,7 @@ int main(int argc, char *argv[])
                         );
                     }
 
-                    // Re-read procaddressing
+                    // Re-read procAddressing
                     distMap =
                         fvMeshTools::readProcAddressing(mesh, baseMeshPtr);
 
@@ -1681,6 +1976,17 @@ int main(int argc, char *argv[])
                     );
 
                     lagrangianDistributorPtr.reset();
+
+                    if (areaMeshPtr)
+                    {
+                        Info<< "    Discarding finite-area addressing"
+                            << " (TODO)" << nl << endl;
+
+                        areaBaseMeshPtr.reset();
+                        areaMeshPtr.reset();
+                        faDistributor.reset();
+                        areaDistMap.clear();
+                    }
                 }
 
 
@@ -1708,6 +2014,12 @@ int main(int argc, char *argv[])
                     selectedLagrangianFields
                 );
 
+                if (faDistributor)
+                {
+                    faDistributor()
+                        .distributeAllFields(objects, selectedFields);
+                }
+
                 // If there are any "uniform" directories copy them from
                 // the master processor
                 copyUniform
@@ -1772,21 +2084,16 @@ int main(int argc, char *argv[])
         forAll(regionNames, regioni)
         {
             const word& regionName = regionNames[regioni];
-            const fileName meshSubDir
+            const word& regionDir =
             (
-                regionName == polyMesh::defaultRegion
-              ? fileName(polyMesh::meshSubDir)
-              : regionNames[regioni]/polyMesh::meshSubDir
+                regionName == polyMesh::defaultRegion ? word::null : regionName
             );
+            const fileName volMeshSubDir(regionDir/polyMesh::meshSubDir);
+            const fileName areaMeshSubDir(regionDir/faMesh::meshSubDir);
 
-            if (decompose)
-            {
-                Info<< "\n\nDecomposing mesh " << regionName << nl << endl;
-            }
-            else
-            {
-                Info<< "\n\nRedistributing mesh " << regionName << nl << endl;
-            }
+            Info<< nl << nl
+                << (decompose ? "Decomposing" : "Redistributing")
+                << " mesh " << regionDir << nl << endl;
 
 
             // Get time instance directory
@@ -1795,62 +2102,191 @@ int main(int argc, char *argv[])
             // processor0. Note the changing of the processor0 casename to
             // enforce it to read/write from the undecomposed case
 
-            fileName masterInstDir;
+            fileName volMeshMasterInstance;
+            fileName areaMeshMasterInstance;
+
+            // Assume to be true
+            bool volMeshHaveUndecomposed = true;
+            bool areaMeshHaveUndecomposed = doFiniteArea;
+
             if (Pstream::master())
             {
                 if (decompose)
                 {
-                    Info<< "Setting caseName to " << baseRunTime.caseName()
-                        << " to find undecomposed mesh" << endl;
+                    Info<< "Checking undecomposed mesh in case: "
+                        << baseRunTime.caseName() << endl;
                     runTime.caseName() = baseRunTime.caseName();
                     runTime.processorCase(false);
                 }
 
                 const bool oldParRun = Pstream::parRun(false);
-                masterInstDir = runTime.findInstance
+                volMeshMasterInstance = runTime.findInstance
                 (
-                    meshSubDir,
+                    volMeshSubDir,
                     "faces",
                     IOobject::READ_IF_PRESENT
                 );
-                Pstream::parRun(oldParRun);
+
+                if (doFiniteArea)
+                {
+                    areaMeshMasterInstance = runTime.findInstance
+                    (
+                        areaMeshSubDir,
+                        "faceLabels",
+                        IOobject::READ_IF_PRESENT
+                    );
+
+                    // Note: findInstance returns "constant" even if not found,
+                    // so recheck now for a false positive.
+
+                    if ("constant" == areaMeshMasterInstance)
+                    {
+                        const boolList areaMeshOnProc
+                        (
+                            haveMeshFile
+                            (
+                                runTime,
+                                areaMeshMasterInstance/areaMeshSubDir,
+                                "faceLabels",
+                                false  // verbose=false
+                            )
+                        );
+
+                        if (areaMeshOnProc.empty() || !areaMeshOnProc[0])
+                        {
+                            areaMeshHaveUndecomposed = false;
+                        }
+                    }
+                }
+
+                Pstream::parRun(oldParRun);  // Restore parallel state
 
                 if (decompose)
                 {
-                    Info<< "Restoring caseName" << endl;
+                    Info<< "    volume mesh ["
+                        << volMeshHaveUndecomposed << "] : "
+                        << volMeshMasterInstance << nl
+                        << "    area   mesh ["
+                        << areaMeshHaveUndecomposed << "] : "
+                        << areaMeshMasterInstance << nl
+                        << nl << nl;
+
+                    // Restoring caseName
                     runTime.caseName() = proc0CaseName;
                     runTime.processorCase(oldProcCase);
                 }
             }
-            Pstream::scatter(masterInstDir);
 
-            // Check who has a polyMesh
-            const boolList haveMesh
+            Pstream::broadcasts
+            (
+                UPstream::worldComm,
+                volMeshHaveUndecomposed,
+                areaMeshHaveUndecomposed,
+                volMeshMasterInstance,
+                areaMeshMasterInstance
+            );
+
+            // Check processors have meshes
+            // - check for 'faces' file (polyMesh)
+            // - check for 'faceLabels' file (faMesh)
+            boolList volMeshOnProc;
+            boolList areaMeshOnProc;
+
+            volMeshOnProc = haveMeshFile
             (
-                haveMeshFile
+                runTime,
+                volMeshMasterInstance/volMeshSubDir,
+                "faces"
+            );
+
+            if (doFiniteArea)
+            {
+                areaMeshOnProc = haveMeshFile
                 (
                     runTime,
-                    masterInstDir/meshSubDir,
-                    "faces"
-                )
-            );
+                    areaMeshMasterInstance/areaMeshSubDir,
+                    "faceLabels"
+                );
+            }
+
+            // Prior to loadOrCreateMesh, note which meshes already exist
+            // for the current file handler.
+            // - where mesh would be written if it didn't exist already.
+            fileNameList volMeshDir(Pstream::nProcs());
+            {
+                volMeshDir[Pstream::myProcNo()] =
+                (
+                    fileHandler().objectPath
+                    (
+                        IOobject
+                        (
+                            "faces",
+                            volMeshMasterInstance/volMeshSubDir,
+                            runTime
+                        ),
+                        word::null
+                    ).path()
+                );
+
+                Pstream::allGatherList(volMeshDir);
+
+                if (optVerbose && Pstream::master())
+                {
+                    Info<< "Per processor faces dirs:" << nl
+                        << '(' << nl;
+
+                    for (const int proci : Pstream::allProcs())
+                    {
+                        Info<< "    "
+                            << runTime.relativePath(volMeshDir[proci]);
 
-            // Collect objectPath of polyMesh for the current file handler. This
-            // is where the mesh would be written if it didn't exist already.
-            fileNameList meshDir(Pstream::nProcs());
+                        if (!volMeshOnProc[proci])
+                        {
+                            Info<< " [missing]";
+                        }
+                        Info<< nl;
+                    }
+                    Info<< ')' << nl << endl;
+                }
+            }
+
+            fileNameList areaMeshDir(Pstream::nProcs());
+            if (doFiniteArea)
             {
-                const fileName fName
+                areaMeshDir[Pstream::myProcNo()] =
                 (
                     fileHandler().objectPath
                     (
-                        IOobject("faces", masterInstDir/meshSubDir, runTime),
+                        IOobject
+                        (
+                            "faceLabels",
+                            areaMeshMasterInstance/areaMeshSubDir,
+                            runTime
+                        ),
                         word::null
-                    )
+                    ).path()
                 );
-                meshDir[Pstream::myProcNo()] = fName.path();
-                Pstream::allGatherList(meshDir);
-                //Info<< "Per processor faces dirs:" << nl
-                //    << "    " << meshDir << nl << endl;
+
+                Pstream::allGatherList(areaMeshDir);
+
+                if (optVerbose && Pstream::master())
+                {
+                    Info<< "Per processor faceLabels dirs:" << nl
+                        << '(' << nl;
+
+                    for (const int proci : Pstream::allProcs())
+                    {
+                        Info<< "    "
+                            << runTime.relativePath(areaMeshDir[proci]);
+
+                        if (!areaMeshOnProc[proci])
+                        {
+                            Info<< " [missing]";
+                        }
+                        Info<< nl;
+                    }
+                    Info<< ')' << nl << endl;
+                }
             }
 
 
@@ -1865,18 +2301,55 @@ int main(int argc, char *argv[])
                 runTime.processorCase(false);
             }
 
-            autoPtr<fvMesh> meshPtr = fvMeshTools::loadOrCreateMesh
+            // Volume mesh
+            autoPtr<fvMesh> volMeshPtr = fvMeshTools::loadOrCreateMesh
             (
                 IOobject
                 (
                     regionName,
-                    masterInstDir,
+                    volMeshMasterInstance,
                     runTime,
                     Foam::IOobject::MUST_READ
                 ),
                 decompose
             );
-            fvMesh& mesh = meshPtr();
+            fvMesh& mesh = volMeshPtr();
+
+
+            // Area mesh
+
+            autoPtr<faMesh> areaMeshPtr;
+
+            // Decomposing: must have an undecomposed mesh
+            // Redistributing: have any proc mesh
+            if
+            (
+                doFiniteArea
+             &&
+                (
+                    decompose
+                  ? areaMeshHaveUndecomposed
+                  : areaMeshOnProc.found(true)
+                )
+            )
+            {
+                areaMeshPtr = faMeshTools::loadOrCreateMesh
+                (
+                    IOobject
+                    (
+                        regionName,
+                        areaMeshMasterInstance,
+                        runTime,
+                        Foam::IOobject::MUST_READ
+                    ),
+                    mesh,  // <- The referenced polyMesh (from above)
+                    decompose
+                );
+
+                faMeshTools::forceDemandDriven(*areaMeshPtr);
+                faMeshTools::unregisterMesh(*areaMeshPtr);
+            }
+
 
             if (writeHandler)
             {
@@ -1895,13 +2368,24 @@ int main(int argc, char *argv[])
                     {
                         if
                         (
-                           !haveMesh[proci]
-                         && meshDir[proci] != meshDir[myProci]
+                           !volMeshOnProc[proci]
+                         && volMeshDir[proci] != volMeshDir[myProci]
                         )
                         {
-                            Info<< "Deleting mesh dir:" << meshDir[proci]
-                                << endl;
-                            rmDir(meshDir[proci]);
+                            Info<< "Deleting mesh dir:"
+                                << volMeshDir[proci] << endl;
+                            Foam::rmDir(volMeshDir[proci]);
+                        }
+
+                        if
+                        (
+                            !areaMeshOnProc[proci]
+                         && areaMeshDir[proci] != areaMeshDir[myProci]
+                        )
+                        {
+                            Info<< "Deleting mesh dir:"
+                                << areaMeshDir[proci] << endl;
+                            Foam::rmDir(areaMeshDir[proci]);
                         }
                     }
 
@@ -1921,9 +2405,13 @@ int main(int argc, char *argv[])
             }
 
             const label nOldCells = mesh.nCells();
+
+            // const label nOldAreaFaces =
+            //     (areaMeshPtr ? areaMeshPtr().nFaces() : 0);
+            //
             //Pout<< "Loaded mesh : nCells:" << nOldCells
             //    << " nPatches:" << mesh.boundaryMesh().size() << endl;
-
+            //Pout<< "Loaded area mesh : nFaces:" << nOldAreaFaces << endl;
 
             // Determine decomposition
             // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -1943,19 +2431,47 @@ int main(int argc, char *argv[])
                 finalDecomp
             );
 
+
             if (dryrun)
             {
-                if (!Pstream::master() && !haveMesh[Pstream::myProcNo()])
+                if (!Pstream::master())
                 {
-                    // Remove dummy mesh created by loadOrCreateMesh
-                    const bool oldParRun = Pstream::parRun(false);
-                    mesh.removeFiles();
-                    rmDir(mesh.objectRegistry::objectPath());
-                    Pstream::parRun(oldParRun);  // Restore parallel state
+                    if (areaMeshPtr && !areaMeshOnProc[Pstream::myProcNo()])
+                    {
+                        // Remove dummy mesh created by loadOrCreateMesh
+                        const bool oldParRun = Pstream::parRun(false);
+                        areaMeshPtr->removeFiles();
+                        Pstream::parRun(oldParRun);  // Restore parallel state
+                    }
+
+                    if (!volMeshOnProc[Pstream::myProcNo()])
+                    {
+                        // Remove dummy mesh created by loadOrCreateMesh
+                        const bool oldParRun = Pstream::parRun(false);
+                        mesh.removeFiles();
+                        rmDir(mesh.objectRegistry::objectPath());
+                        Pstream::parRun(oldParRun);  // Restore parallel state
+                    }
                 }
                 continue;
             }
 
+            // Area fields first. Read and deregister
+            parFaFieldDistributorCache areaFields;
+            if (areaMeshPtr)
+            {
+                areaFields.read
+                (
+                    baseRunTime,
+                    proc0CaseName,
+                    decompose,
+
+                    areaMeshOnProc,
+                    areaMeshMasterInstance,
+                    (*areaMeshPtr)
+                );
+            }
+
 
             // Detect lagrangian fields
             if (Pstream::master() && decompose)
@@ -1986,16 +2502,21 @@ int main(int argc, char *argv[])
             (
                 std::move(writeHandler),
                 baseRunTime,
-                haveMesh,
-                meshSubDir,
+                proc0CaseName,
+
+                // Controls
                 true,           // read fields
                 decompose,      // decompose, i.e. read from undecomposed case
                 false,          // no reconstruction
                 overwrite,
-                proc0CaseName,
+
+                // Decomposition information
                 nDestProcs,
                 finalDecomp,
-                masterInstDir,
+
+                // For finite volume
+                volMeshOnProc,
+                volMeshMasterInstance,
                 mesh
             );
 
@@ -2011,6 +2532,100 @@ int main(int argc, char *argv[])
             );
 
 
+            // Redistribute area fields
+
+            mapDistributePolyMesh faDistMap;
+            autoPtr<faMesh> areaProcMeshPtr;
+
+            if (areaMeshPtr)
+            {
+                faDistMap = faMeshDistributor::distribute
+                (
+                    areaMeshPtr(),
+                    distMap(),
+                    areaProcMeshPtr
+                );
+
+                // Force recreation of everything that might vaguely
+                // be used by patches:
+
+                faMeshTools::forceDemandDriven(areaProcMeshPtr());
+
+
+                if (reconstruct)
+                {
+                    if (Pstream::master())
+                    {
+                        Info<< "Setting caseName to " << baseRunTime.caseName()
+                            << " to write reconstructed mesh (and fields)."
+                            << endl;
+                        runTime.caseName() = baseRunTime.caseName();
+                        const bool oldProcCase(runTime.processorCase(false));
+                        //const bool oldParRun = Pstream::parRun(false);
+
+                        areaProcMeshPtr->write();
+
+                        // Now we've written all. Reset caseName on master
+                        Info<< "Restoring caseName" << endl;
+                        runTime.caseName() = proc0CaseName;
+                        runTime.processorCase(oldProcCase);
+                    }
+                }
+                else
+                {
+                    autoPtr<fileOperation> defaultHandler;
+                    if (writeHandler)
+                    {
+                        defaultHandler = fileHandler(std::move(writeHandler));
+                    }
+
+                    IOmapDistributePolyMeshRef
+                    (
+                        IOobject
+                        (
+                            "procAddressing",
+                            areaProcMeshPtr->facesInstance(),
+                            faMesh::meshSubDir,
+                            areaProcMeshPtr->thisDb(),
+                            IOobject::NO_READ,
+                            IOobject::NO_WRITE,
+                            false
+                        ),
+                        faDistMap
+                    ).write();
+
+                    areaProcMeshPtr->write();
+
+                    if (defaultHandler)
+                    {
+                        writeHandler = fileHandler(std::move(defaultHandler));
+                    }
+
+                    if (decompose)
+                    {
+                        faMeshTools::writeProcAddressing
+                        (
+                            areaProcMeshPtr(),
+                            faDistMap,
+                            decompose,
+                            std::move(writeHandler)
+                        );
+                    }
+                }
+
+                Info<< "Written redistributed mesh to "
+                    << areaProcMeshPtr->facesInstance() << nl << endl;
+
+                faMeshDistributor distributor
+                (
+                    areaMeshPtr(),      // source
+                    areaProcMeshPtr(),  // target
+                    faDistMap
+                );
+
+                areaFields.redistributeAndWrite(distributor, true);
+            }
+
             // Copy region-specific uniform
             // (e.g. solid/uniform/cumulativeContErr)
             copyUniform
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H
index c02ed434f1a27cef9fc9587c5e3aac5a2ba4d428..0617e7e7d6a9ebc32a0c640cb8db5c66f5c6e43b 100644
--- a/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H
@@ -186,6 +186,16 @@ public:
             const bool allowUnmapped = false
         );
 
+        //- Map edge field.
+        //  Optionally allow unmapped faces not to produce a warning
+        template<class Type>
+        static tmp<GeometricField<Type, faePatchField, edgeMesh>>
+        interpolate
+        (
+            const GeometricField<Type, faePatchField, edgeMesh>&,
+            const faMesh& sMesh
+        );
+
 
     // Field Mapping
 
@@ -198,6 +208,15 @@ public:
             const GeometricField<Type, faPatchField, areaMesh>&,
             const bool allowUnmapped = false
         ) const;
+
+        //- Map edge field.
+        template<class Type>
+        tmp<GeometricField<Type, faePatchField, edgeMesh>>
+        interpolate
+        (
+            const GeometricField<Type, faePatchField, edgeMesh>&,
+            const bool allowUnmapped = false
+        ) const;
 };
 
 
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C
index c245505ca7d4dd7e18ae2c339a954b28e572271c..b6bcc5f85a11348eb716fc007212e60f580f89cc 100644
--- a/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C
@@ -144,6 +144,88 @@ Foam::faMeshSubset::interpolate
 }
 
 
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>
+>
+Foam::faMeshSubset::interpolate
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& vf,
+    const faMesh& sMesh
+)
+{
+    // 1. Create the complete field with dummy patch fields
+    PtrList<faePatchField<Type>> patchFields(sMesh.boundary().size());
+
+    forAll(patchFields, patchi)
+    {
+        patchFields.set
+        (
+            patchi,
+            faePatchField<Type>::New
+            (
+                calculatedFaePatchField<Type>::typeName,
+                sMesh.boundary()[patchi],
+                DimensionedField<Type, edgeMesh>::null()
+            )
+        );
+    }
+
+    auto tresult = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
+    (
+        IOobject
+        (
+            "subset"+vf.name(),
+            sMesh.time().timeName(),
+            sMesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sMesh,
+        vf.dimensions(),
+        Field<Type>(),
+        // Field<Type>
+        // (
+        //     vf.primitiveField(),
+        //     SubList<label>(edgeMap, sMesh.nInternalEdges())
+        // ),
+        patchFields
+    );
+    auto& result = tresult.ref();
+    result.oriented() = vf.oriented();
+
+
+    // 2. Change the faePatchFields to the correct type using a mapper
+    //  constructor (with reference to the now correct internal field)
+
+    auto& bf = result.boundaryFieldRef();
+
+    forAll(bf, patchi)
+    {
+        // Construct addressing
+        const faPatch& subPatch = sMesh.boundary()[patchi];
+
+        labelList directAddressing;
+        directFaPatchFieldMapper mapper(directAddressing);
+
+        bf.set
+        (
+            patchi,
+            faePatchField<Type>::New
+            (
+                vf.boundaryField()[patchi],
+                subPatch,
+                result(),
+                mapper
+            )
+        );
+    }
+
+    return tresult;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -166,4 +248,24 @@ Foam::faMeshSubset::interpolate
 }
 
 
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>
+>
+Foam::faMeshSubset::interpolate
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& vf,
+    const bool allowUnmapped
+) const
+{
+    if (subMeshPtr_)
+    {
+        return interpolate(vf, *subMeshPtr_);
+    }
+
+    return vf;
+}
+
+
 // ************************************************************************* //