diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshTools.C b/src/finiteArea/faMesh/faMeshTools/faMeshTools.C
index cce1e3c6d7..102dd4046b 100644
--- a/src/finiteArea/faMesh/faMeshTools/faMeshTools.C
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshTools.C
@@ -30,6 +30,8 @@ License
 #include "faBoundaryMeshEntries.H"
 #include "areaFields.H"
 #include "edgeFields.H"
+#include "fileOperation.H"
+#include "BitOps.H"
 #include "polyMesh.H"
 #include "processorFaPatch.H"
 
@@ -66,7 +68,8 @@ void Foam::faMeshTools::forceDemandDriven(faMesh& mesh)
 }
 
 
-Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
+Foam::autoPtr<Foam::faMesh>
+Foam::faMeshTools::newMesh
 (
     const IOobject& io,
     const polyMesh& pMesh,
@@ -248,9 +251,11 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
 }
 
 
-Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
+Foam::autoPtr<Foam::faMesh>
+Foam::faMeshTools::loadOrCreateMeshImpl
 (
     const IOobject& io,
+    refPtr<fileOperation>* readHandlerPtr,  // Can be nullptr
     const polyMesh& pMesh,
     const bool decompose,
     const bool verbose
@@ -270,9 +275,11 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
     // Read and scatter master patches (without reading master mesh!)
 
     PtrList<entry> patchEntries;
-    if (Pstream::master())
+    if (UPstream::master())
     {
-        const bool oldParRun = Pstream::parRun(false);
+        const bool oldParRun = UPstream::parRun(false);
+        const label oldNumProcs = fileHandler().nProcs();
+        const int oldCache = fileOperation::cacheLevel(0);
 
         patchEntries = faBoundaryMeshEntries
         (
@@ -288,46 +295,75 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
             )
         );
 
-        Pstream::parRun(oldParRun);
+        fileOperation::cacheLevel(oldCache);
+        if (oldParRun)
+        {
+            const_cast<fileOperation&>(fileHandler()).nProcs(oldNumProcs);
+        }
+        UPstream::parRun(oldParRun);
     }
 
     // Broadcast: send patches to all
-    Pstream::broadcast(patchEntries);  // == worldComm;
+    Pstream::broadcast(patchEntries, UPstream::worldComm);
 
-    /// Info<< patchEntries << nl;
-
-    // Dummy meshes
-    // ~~~~~~~~~~~~
 
     // Check who has or needs a mesh.
-    // For 'decompose', only need mesh on master.
-    // Otherwise check for presence of the "faceLabels" file
+    bool haveLocalMesh = false;
 
-    bool haveMesh =
-    (
-        decompose
-      ? Pstream::master()
-      : fileHandler().isFile
+    if (readHandlerPtr)
+    {
+        // Non-null reference when a mesh exists on given processor
+        haveLocalMesh = (*readHandlerPtr).good();
+    }
+    else
+    {
+        // No file handler.
+        // For 'decompose', only need mesh on master.
+        // Otherwise check for presence of the "faceLabels" file
+
+        haveLocalMesh =
         (
-            fileHandler().filePath
+            decompose
+          ? UPstream::master()
+          : fileHandler().isFile
             (
-                io.time().path()/io.instance()/meshSubDir/"faceLabels"
+                fileHandler().filePath
+                (
+                    io.time().path()/io.instance()/meshSubDir/"faceLabels"
+                )
             )
-        )
+        );
+    }
+
+
+    // Globally consistent information about who has a mesh
+    boolList haveMesh
+    (
+        UPstream::allGatherValues<bool>(haveLocalMesh)
     );
 
 
-    if (!haveMesh)
+    autoPtr<faMesh> meshPtr;
+
+    if (!haveLocalMesh)
     {
-        const bool oldParRun = Pstream::parRun(false);
+        // No local mesh - need to synthesize one
+
+        const bool oldParRun = UPstream::parRun(false);
+        const label oldNumProcs = fileHandler().nProcs();
+        const int oldCache = fileOperation::cacheLevel(0);
 
         // Create dummy mesh - on procs that don't already have a mesh
-        faMesh dummyMesh
+        meshPtr.reset
         (
-            pMesh,
-            labelList(),
-            IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE)
+            new faMesh
+            (
+                pMesh,
+                labelList(),
+                IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE)
+            )
         );
+        faMesh& mesh = *meshPtr;
 
         // Add patches
         faPatchList patches(patchEntries.size());
@@ -361,39 +397,119 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
                         name,
                         patchDict,
                         nPatches++,
-                        dummyMesh.boundary()
+                        mesh.boundary()
                     )
                 );
             }
         }
-
         patches.resize(nPatches);
-        dummyMesh.addFaPatches(patches, false);  // No parallel comms
+        mesh.addFaPatches(patches, false);  // No parallel comms
 
-        // Bad hack, but the underlying polyMesh is NO_WRITE
-        // so it does not create the faMesh subDir for us...
-        Foam::mkDir(dummyMesh.boundary().path());
+        if (!readHandlerPtr)
+        {
+            // The 'old' way of doing things.
+            // Write the dummy mesh to disk for subsequent re-reading.
+            //
+            // This is not particularly elegant.
 
-        //Pout<< "Writing dummy mesh to " << dummyMesh.boundary().path()
-        //    << endl;
-        dummyMesh.write();
+            // Bad hack, but the underlying polyMesh is NO_WRITE
+            // so it does not create the faMesh subDir for us...
+            Foam::mkDir(mesh.boundary().path());
 
-        Pstream::parRun(oldParRun);  // Restore parallel state
+            //Pout<< "Writing dummy mesh to " << mesh.boundary().path() << nl;
+            mesh.write();
+
+            // Discard - it will be re-read later
+            meshPtr.reset(nullptr);
+        }
+
+        fileOperation::cacheLevel(oldCache);
+        if (oldParRun)
+        {
+            const_cast<fileOperation&>(fileHandler()).nProcs(oldNumProcs);
+        }
+        UPstream::parRun(oldParRun);  // Restore parallel state
+    }
+    else if (readHandlerPtr && haveLocalMesh)
+    {
+        const labelList meshProcIds(BitOps::sortedToc(haveMesh));
+
+        UPstream::communicator newCommunicator;
+        const label oldWorldComm = UPstream::commWorld();
+
+        auto& readHandler = *readHandlerPtr;
+        auto oldHandler = fileOperation::fileHandler(readHandler);
+
+        // With IO ranks the communicator of the fileOperation will
+        // only include the ranks for the current IO rank.
+        // Instead allocate a new communicator for everyone with a mesh
+
+        const auto& handlerProcIds = UPstream::procID(fileHandler().comm());
+
+        // Comparing global ranks in the communicator.
+        // Use std::equal for the List<label> vs List<int> comparison
+
+        if
+        (
+            meshProcIds.size() == handlerProcIds.size()
+         && std::equal
+            (
+                meshProcIds.cbegin(),
+                meshProcIds.cend(),
+                handlerProcIds.cbegin()
+            )
+        )
+        {
+            // Can use the handler communicator as is.
+            UPstream::commWorld(fileHandler().comm());
+        }
+        else if
+        (
+            UPstream::nProcs(fileHandler().comm())
+         != UPstream::nProcs(UPstream::worldComm)
+        )
+        {
+            // Need a new communicator for the fileHandler.
+
+            // Warning: MS-MPI currently uses MPI_Comm_create() instead of
+            // MPI_Comm_create_group() so it will block here!
+
+            newCommunicator.reset(UPstream::worldComm, meshProcIds);
+            UPstream::commWorld(newCommunicator.comm());
+        }
+
+        // Load but do not initialise
+        meshPtr = autoPtr<faMesh>::New(pMesh, labelList(), io);
+
+        readHandler = fileOperation::fileHandler(oldHandler);
+        UPstream::commWorld(oldWorldComm);
+
+        // Reset mesh communicator to the real world comm
+        meshPtr().comm() = UPstream::commWorld();
     }
 
-    // Read mesh
-    // ~~~~~~~~~
-    // Now all processors have a (possibly zero size) mesh so read in
-    // parallel
 
-    /// Pout<< "Reading area mesh from " << io.objectRelPath() << endl;
-    auto meshPtr = autoPtr<faMesh>::New(pMesh, false);
-    faMesh& mesh = *meshPtr;
+    if (!meshPtr)
+    {
+        // Using the 'old' way of doing things (writing to disk and re-reading).
 
-    // Sync patches
-    // ~~~~~~~~~~~~
+        // Read mesh from disk
+        //
+        // Now all processors have a (possibly zero size) mesh so can
+        // read in parallel
 
-    if (!Pstream::master() && haveMesh)
+        /// Pout<< "Reading area mesh from " << io.objectRelPath() << endl;
+        // Load but do not initialise
+        meshPtr = autoPtr<faMesh>::New(pMesh, false);
+    }
+
+    faMesh& mesh = meshPtr();
+
+    // Check patches
+    // ~~~~~~~~~~~~~
+
+    #if 0
+    if (!UPstream::master() && haveLocalMesh)
     {
         // Check master names against mine
 
@@ -416,9 +532,8 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
             if (patchi >= patches.size())
             {
                 FatalErrorInFunction
-                    << "Non-processor patches not synchronised."
-                    << endl
-                    << "Processor " << Pstream::myProcNo()
+                    << "Non-processor patches not synchronised." << endl
+                    << "Processor " << UPstream::myProcNo()
                     << " has only " << patches.size()
                     << " patches, master has "
                     << patchi
@@ -432,12 +547,11 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
             )
             {
                 FatalErrorInFunction
-                    << "Non-processor patches not synchronised."
-                    << endl
+                    << "Non-processor patches not synchronised." << endl
                     << "Master patch " << patchi
                     << " name:" << type
                     << " type:" << type << endl
-                    << "Processor " << Pstream::myProcNo()
+                    << "Processor " << UPstream::myProcNo()
                     << " patch " << patchi
                     << " has name:" << patches[patchi].name()
                     << " type:" << patches[patchi].type()
@@ -445,6 +559,7 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
             }
         }
     }
+    #endif
 
 
     // Recreate basic geometry, globalMeshData etc.
@@ -467,4 +582,44 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
 }
 
 
+Foam::autoPtr<Foam::faMesh>
+Foam::faMeshTools::loadOrCreateMesh
+(
+    const IOobject& io,
+    const polyMesh& pMesh,
+    const bool decompose,
+    const bool verbose
+)
+{
+    return faMeshTools::loadOrCreateMeshImpl
+    (
+        io,
+        nullptr,  // fileOperation (ignore)
+        pMesh,
+        decompose,
+        verbose
+    );
+}
+
+
+Foam::autoPtr<Foam::faMesh>
+Foam::faMeshTools::loadOrCreateMesh
+(
+    const IOobject& io,
+    const polyMesh& pMesh,
+    refPtr<fileOperation>& readHandler,
+    const bool verbose
+)
+{
+    return faMeshTools::loadOrCreateMeshImpl
+    (
+        io,
+       &readHandler,
+        pMesh,
+        false,  // decompose (ignored)
+        verbose
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshTools.H b/src/finiteArea/faMesh/faMeshTools/faMeshTools.H
index 0906f5d550..ee63b1b305 100644
--- a/src/finiteArea/faMesh/faMeshTools/faMeshTools.H
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshTools.H
@@ -61,6 +61,26 @@ class IOobject;
 
 class faMeshTools
 {
+    // Private Member Functions
+
+    // Read mesh if available. Otherwise create empty mesh with same non-proc
+    // patches as proc0 mesh.
+    //  Has two modes of operation.
+    //  If the readHandler is non-nullptr, use it to decide on availability.
+    //  Requires:
+    //  - all processors to have all patches (and in same order).
+    //  - io.instance() set to facesInstance
+    static autoPtr<faMesh> loadOrCreateMeshImpl
+    (
+        const IOobject& io,
+        refPtr<fileOperation>* readHandlerPtr,  // Can be nullptr
+        const polyMesh& pMesh,
+        const bool decompose,  // Only used if readHandlerPtr == nullptr
+        const bool verbose = false
+    );
+
+
+
 public:
 
     //- Unregister the faMesh from its associated polyMesh
@@ -97,6 +117,19 @@ public:
         const bool verbose = false
     );
 
+    // Read mesh if available. Otherwise create empty mesh with same non-proc
+    // patches as proc0 mesh. Requires:
+    //  - all processors to have all patches (and in same order).
+    //  - io.instance() set to facesInstance
+    static autoPtr<faMesh> loadOrCreateMesh
+    (
+        const IOobject& io,
+        const polyMesh& pMesh,
+        //! Non-null reference if a mesh exists on given processor
+        refPtr<fileOperation>& readHandler,
+        const bool verbose = false
+    );
+
 
     //- Read decompose/reconstruct addressing
     static mapDistributePolyMesh readProcAddressing
diff --git a/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.C b/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.C
index a537c114bb..14f154017d 100644
--- a/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.C
+++ b/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,8 @@ License
 #include "pointSet.H"
 #include "faceSet.H"
 #include "cellSet.H"
+#include "fileOperation.H"
+#include "BitOps.H"
 #include "IOobjectList.H"
 #include "basicFvGeometryScheme.H"
 #include "processorPolyPatch.H"
@@ -471,7 +473,7 @@ Foam::fvMeshTools::newMesh
     // Read and scatter master patches (without reading master mesh!)
 
     PtrList<entry> patchEntries;
-    if (Pstream::master())
+    if (UPstream::master())
     {
         const bool oldParRun = Pstream::parRun(false);
 
@@ -502,7 +504,7 @@ Foam::fvMeshTools::newMesh
             )
         );
 
-        Pstream::parRun(oldParRun);
+        UPstream::parRun(oldParRun);
     }
 
     // Broadcast information to all
@@ -730,14 +732,14 @@ Foam::fvMeshTools::newMesh
 
 
 Foam::autoPtr<Foam::fvMesh>
-Foam::fvMeshTools::loadOrCreateMesh
+Foam::fvMeshTools::loadOrCreateMeshImpl
 (
     const IOobject& io,
+    refPtr<fileOperation>* readHandlerPtr,  // Can be nullptr
     const bool decompose,
     const bool verbose
 )
 {
-
     // Region name
     // ~~~~~~~~~~~
 
@@ -752,9 +754,11 @@ Foam::fvMeshTools::loadOrCreateMesh
     // Read and scatter master patches (without reading master mesh!)
 
     PtrList<entry> patchEntries;
-    if (Pstream::master())
+    if (UPstream::master())
     {
-        const bool oldParRun = Pstream::parRun(false);
+        const bool oldParRun = UPstream::parRun(false);
+        const label oldNumProcs = fileHandler().nProcs();
+        const int oldCache = fileOperation::cacheLevel(0);
 
         patchEntries = polyBoundaryMeshEntries
         (
@@ -770,45 +774,74 @@ Foam::fvMeshTools::loadOrCreateMesh
             )
         );
 
-        Pstream::parRun(oldParRun);
+        fileOperation::cacheLevel(oldCache);
+        if (oldParRun)
+        {
+            const_cast<fileOperation&>(fileHandler()).nProcs(oldNumProcs);
+        }
+        UPstream::parRun(oldParRun);
     }
 
     // Broadcast: send patches to all
-    Pstream::broadcast(patchEntries);  // == worldComm;
+    Pstream::broadcast(patchEntries, UPstream::worldComm);
 
 
-    // Dummy meshes
-    // ~~~~~~~~~~~~
-
     // Check who has or needs a mesh.
-    // For 'decompose', only need mesh on master.
-    // Otherwise check for presence of the "faces" file
+    bool haveLocalMesh = false;
 
-    bool haveMesh =
-    (
-        decompose
-      ? Pstream::master()
-      : fileHandler().isFile
+    if (readHandlerPtr)
+    {
+        // Non-null reference when a mesh exists on given processor
+        haveLocalMesh = (*readHandlerPtr).good();
+    }
+    else
+    {
+        // No file handler.
+        // Check for presence of the "faces" file,
+        // but for 'decompose', only need mesh on master.
+
+        haveLocalMesh =
         (
-            fileHandler().filePath
+            decompose
+          ? UPstream::master()
+          : fileHandler().isFile
             (
-                io.time().path()/io.instance()/meshSubDir/"faces"
+                fileHandler().filePath
+                (
+                    io.time().path()/io.instance()/meshSubDir/"faces"
+                )
             )
-        )
+        );
+    }
+
+    // Globally consistent information about who has a mesh
+    boolList haveMesh
+    (
+        UPstream::allGatherValues<bool>(haveLocalMesh)
     );
 
 
-    if (!haveMesh)
+    autoPtr<fvMesh> meshPtr;
+
+    if (!haveLocalMesh)
     {
-        const bool oldParRun = Pstream::parRun(false);
+        // No local mesh - need to synthesize one
+
+        const bool oldParRun = UPstream::parRun(false);
+        const label oldNumProcs = fileHandler().nProcs();
+        const int oldCache = fileOperation::cacheLevel(0);
 
         // Create dummy mesh - on procs that don't already have a mesh
-        fvMesh dummyMesh
+        meshPtr.reset
         (
-            IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE),
-            Foam::zero{},
-            false
+            new fvMesh
+            (
+                IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE),
+                Foam::zero{},
+                false
+            )
         );
+        fvMesh& mesh = *meshPtr;
 
         // Add patches
         polyPatchList patches(patchEntries.size());
@@ -844,65 +877,142 @@ Foam::fvMeshTools::loadOrCreateMesh
                         name,
                         patchDict,
                         nPatches++,
-                        dummyMesh.boundaryMesh()
+                        mesh.boundaryMesh()
                     )
                 );
             }
         }
         patches.resize(nPatches);
-        dummyMesh.addFvPatches(patches, false);  // No parallel comms
+        mesh.addFvPatches(patches, false);  // No parallel comms
 
+        if (!readHandlerPtr)
+        {
+            // The 'old' way of doing things.
+            // Write the dummy mesh to disk for subsequent re-reading.
+            //
+            // This is not particularly elegant.
+            //
+            // Note: add some dummy zones so upon reading it does not read them
+            // from the undecomposed case. Should be done as extra argument to
+            // regIOobject::readStream?
 
-        // Add some dummy zones so upon reading it does not read them
-        // from the undecomposed case. Should be done as extra argument to
-        // regIOobject::readStream?
+            List<pointZone*> pz
+            (
+                1,
+                new pointZone("dummyZone", 0, mesh.pointZones())
+            );
+            List<faceZone*> fz
+            (
+                1,
+                new faceZone("dummyZone", 0, mesh.faceZones())
+            );
+            List<cellZone*> cz
+            (
+                1,
+                new cellZone("dummyZone", 0, mesh.cellZones())
+            );
+            mesh.addZones(pz, fz, cz);
+            mesh.pointZones().clear();
+            mesh.faceZones().clear();
+            mesh.cellZones().clear();
+            //Pout<< "Writing dummy mesh: " << mesh.polyMesh::objectPath()
+            //    << endl;
+            mesh.write();
 
-        List<pointZone*> pz
+            // Discard - it will be re-read later
+            meshPtr.reset(nullptr);
+        }
+
+        fileOperation::cacheLevel(oldCache);
+        if (oldParRun)
+        {
+            const_cast<fileOperation&>(fileHandler()).nProcs(oldNumProcs);
+        }
+        UPstream::parRun(oldParRun);  // Restore parallel state
+    }
+    else if (readHandlerPtr && haveLocalMesh)
+    {
+        const labelList meshProcIds(BitOps::sortedToc(haveMesh));
+
+        UPstream::communicator newCommunicator;
+        const label oldWorldComm = UPstream::commWorld();
+
+        auto& readHandler = *readHandlerPtr;
+        auto oldHandler = fileOperation::fileHandler(readHandler);
+
+        // With IO ranks the communicator of the fileOperation will
+        // only include the ranks for the current IO rank.
+        // Instead allocate a new communicator for everyone with a mesh
+
+        const auto& handlerProcIds = UPstream::procID(fileHandler().comm());
+
+        // Comparing global ranks in the communicator.
+        // Use std::equal for the List<label> vs List<int> comparison
+
+        if
         (
-            1,
-            new pointZone("dummyPointZone", 0, dummyMesh.pointZones())
-        );
-        List<faceZone*> fz
+            meshProcIds.size() == handlerProcIds.size()
+         && std::equal
+            (
+                meshProcIds.cbegin(),
+                meshProcIds.cend(),
+                handlerProcIds.cbegin()
+            )
+        )
+        {
+            // Can use the handler communicator as is.
+            UPstream::commWorld(fileHandler().comm());
+        }
+        else if
         (
-            1,
-            new faceZone("dummyFaceZone", 0, dummyMesh.faceZones())
-        );
-        List<cellZone*> cz
-        (
-            1,
-            new cellZone("dummyCellZone", 0, dummyMesh.cellZones())
-        );
-        dummyMesh.addZones(pz, fz, cz);
-        dummyMesh.pointZones().clear();
-        dummyMesh.faceZones().clear();
-        dummyMesh.cellZones().clear();
-        //Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
-        //    << endl;
-        dummyMesh.write();
+            UPstream::nProcs(fileHandler().comm())
+         != UPstream::nProcs(UPstream::worldComm)
+        )
+        {
+            // Need a new communicator for the fileHandler.
 
-        Pstream::parRun(oldParRun);  // Restore parallel state
+            // Warning: MS-MPI currently uses MPI_Comm_create() instead of
+            // MPI_Comm_create_group() so it will block here!
+
+            newCommunicator.reset(UPstream::worldComm, meshProcIds);
+            UPstream::commWorld(newCommunicator.comm());
+        }
+
+        // Load but do not initialise
+        meshPtr = autoPtr<fvMesh>::New(io, false);
+
+        readHandler = fileOperation::fileHandler(oldHandler);
+        UPstream::commWorld(oldWorldComm);
+
+        // Reset mesh communicator to the real world comm
+        meshPtr().polyMesh::comm() = UPstream::commWorld();
     }
 
 
-    // Read mesh
-    // ~~~~~~~~~
-    // Now all processors have a (possibly zero size) mesh so read in
-    // parallel
-
-    //Pout<< "Reading mesh from " << io.objectRelPath() << endl;
-    auto meshPtr = autoPtr<fvMesh>::New(io);
-    fvMesh& mesh = *meshPtr;
-
-    // Make sure to use a non-parallel geometry calculation method
-    fvMeshTools::setBasicGeometry(mesh);
-
-
-    // Sync patches
-    // ~~~~~~~~~~~~
-
-    if (!Pstream::master() && haveMesh)
+    if (!meshPtr)
     {
-        // Check master names against mine
+        // Using the 'old' way of doing things (writing to disk and re-reading).
+
+        // Read mesh from disk
+        //
+        // Now all processors have a (possibly zero size) mesh so can
+        // read in parallel
+
+        //Pout<< "Reading mesh from " << io.objectRelPath() << endl;
+        // Load but do not initialise
+        meshPtr = autoPtr<fvMesh>::New(io, false);
+    }
+
+    fvMesh& mesh = meshPtr();
+
+
+    // Check patches
+    // ~~~~~~~~~~~~~
+
+    #if 0
+    if (!UPstream::master() && haveLocalMesh)
+    {
+        // Check master patch entries names against local ones
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
@@ -924,12 +1034,11 @@ Foam::fvMeshTools::loadOrCreateMesh
             if (patchi >= patches.size())
             {
                 FatalErrorInFunction
-                    << "Non-processor patches not synchronised."
-                    << endl
-                    << "Processor " << Pstream::myProcNo()
+                    << "Non-processor patches not synchronised." << endl
+                    << "Processor " << UPstream::myProcNo()
                     << " has only " << patches.size()
                     << " patches, master has "
-                    << patchi
+                    << patchi << endl
                     << exit(FatalError);
             }
 
@@ -940,12 +1049,11 @@ Foam::fvMeshTools::loadOrCreateMesh
             )
             {
                 FatalErrorInFunction
-                    << "Non-processor patches not synchronised."
-                    << endl
+                    << "Non-processor patches not synchronised." << endl
                     << "Master patch " << patchi
                     << " name:" << type
                     << " type:" << type << endl
-                    << "Processor " << Pstream::myProcNo()
+                    << "Processor " << UPstream::myProcNo()
                     << " patch " << patchi
                     << " has name:" << patches[patchi].name()
                     << " type:" << patches[patchi].type()
@@ -953,92 +1061,109 @@ Foam::fvMeshTools::loadOrCreateMesh
             }
         }
     }
+    #endif
 
 
-    // Determine zones
-    // ~~~~~~~~~~~~~~~
+    // Synchronize zones
+    // ~~~~~~~~~~~~~~~~~
 
-    wordList pointZoneNames(mesh.pointZones().names());
-    wordList faceZoneNames(mesh.faceZones().names());
-    wordList cellZoneNames(mesh.cellZones().names());
-    Pstream::broadcasts
-    (
-        UPstream::worldComm,
-        pointZoneNames,
-        faceZoneNames,
-        cellZoneNames
-    );
-
-    if (!haveMesh)
     {
-        // Add the zones. Make sure to remove the old dummy ones first
-        mesh.pointZones().clear();
-        mesh.faceZones().clear();
-        mesh.cellZones().clear();
-
-        List<pointZone*> pz(pointZoneNames.size());
-        forAll(pointZoneNames, i)
-        {
-            pz[i] = new pointZone(pointZoneNames[i], i, mesh.pointZones());
-        }
-        List<faceZone*> fz(faceZoneNames.size());
-        forAll(faceZoneNames, i)
-        {
-            fz[i] = new faceZone(faceZoneNames[i], i, mesh.faceZones());
-        }
-        List<cellZone*> cz(cellZoneNames.size());
-        forAll(cellZoneNames, i)
-        {
-            cz[i] = new cellZone(cellZoneNames[i], i, mesh.cellZones());
-        }
-        mesh.addZones(pz, fz, cz);
-    }
+        wordList pointZoneNames(mesh.pointZones().names());
+        wordList faceZoneNames(mesh.faceZones().names());
+        wordList cellZoneNames(mesh.cellZones().names());
+        Pstream::broadcasts
+        (
+            UPstream::worldComm,
+            pointZoneNames,
+            faceZoneNames,
+            cellZoneNames
+        );
 
 
-    // Determine sets
-    // ~~~~~~~~~~~~~~
-
-    wordList pointSetNames;
-    wordList faceSetNames;
-    wordList cellSetNames;
-    if (Pstream::master())
-    {
-        // Read sets
-        const bool oldParRun = Pstream::parRun(false);
-        IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
-        Pstream::parRun(oldParRun);
-
-        pointSetNames = objects.sortedNames<pointSet>();
-        faceSetNames = objects.sortedNames<faceSet>();
-        cellSetNames = objects.sortedNames<cellSet>();
-    }
-    Pstream::broadcasts
-    (
-        UPstream::worldComm,
-        pointSetNames,
-        faceSetNames,
-        cellSetNames
-    );
-
-    if (!haveMesh)
-    {
-        for (const word& setName : pointSetNames)
+        if (!haveLocalMesh)
         {
-            pointSet(mesh, setName, 0).write();
-        }
-        for (const word& setName : faceSetNames)
-        {
-            faceSet(mesh, setName, 0).write();
-        }
-        for (const word& setName : cellSetNames)
-        {
-            cellSet(mesh, setName, 0).write();
+            // Add the zones. Make sure to remove the old dummy ones first
+            mesh.pointZones().clear();
+            mesh.faceZones().clear();
+            mesh.cellZones().clear();
+
+            PtrList<pointZone> pz(pointZoneNames.size());
+            forAll(pointZoneNames, i)
+            {
+                pz.emplace_set(i, pointZoneNames[i], i, mesh.pointZones());
+            }
+
+            PtrList<faceZone> fz(faceZoneNames.size());
+            forAll(faceZoneNames, i)
+            {
+                fz.emplace_set(i, faceZoneNames[i], i, mesh.faceZones());
+            }
+
+            PtrList<cellZone> cz(cellZoneNames.size());
+            forAll(cellZoneNames, i)
+            {
+                cz.emplace_set(i, cellZoneNames[i], i, mesh.cellZones());
+            }
+            mesh.addZones(std::move(pz), std::move(fz), std::move(cz));
         }
     }
 
 
-    // Force recreation of globalMeshData.
+    // Synchronize sets (on disk)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (!readHandlerPtr)
+    {
+        wordList pointSetNames;
+        wordList faceSetNames;
+        wordList cellSetNames;
+        if (UPstream::master())
+        {
+            // Read sets
+            const bool oldParRun = UPstream::parRun(false);
+
+            IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
+            UPstream::parRun(oldParRun);
+
+            pointSetNames = objects.sortedNames<pointSet>();
+            faceSetNames = objects.sortedNames<faceSet>();
+            cellSetNames = objects.sortedNames<cellSet>();
+        }
+        Pstream::broadcasts
+        (
+            UPstream::worldComm,
+            pointSetNames,
+            faceSetNames,
+            cellSetNames
+        );
+
+        if (!haveLocalMesh)
+        {
+            for (const word& setName : pointSetNames)
+            {
+                pointSet(mesh, setName, 0).write();
+            }
+            for (const word& setName : faceSetNames)
+            {
+                faceSet(mesh, setName, 0).write();
+            }
+            for (const word& setName : cellSetNames)
+            {
+                cellSet(mesh, setName, 0).write();
+            }
+        }
+    }
+
+
+    // Meshes have been done without init so now can do full initialisation
+
+    fvMeshTools::setBasicGeometry(mesh);
+    mesh.init(true);
+
+    // Force early recreation of globalMeshData etc
     mesh.globalData();
+    mesh.tetBasePtIs();
+    mesh.geometricD();
 
 
     // Do some checks.
@@ -1060,6 +1185,42 @@ Foam::fvMeshTools::loadOrCreateMesh
 }
 
 
+Foam::autoPtr<Foam::fvMesh>
+Foam::fvMeshTools::loadOrCreateMesh
+(
+    const IOobject& io,
+    const bool decompose,
+    const bool verbose
+)
+{
+    return fvMeshTools::loadOrCreateMeshImpl
+    (
+        io,
+        nullptr,  // fileOperation (ignore)
+        decompose,
+        verbose
+    );
+}
+
+
+Foam::autoPtr<Foam::fvMesh>
+Foam::fvMeshTools::loadOrCreateMesh
+(
+    const IOobject& io,
+    refPtr<fileOperation>& readHandler,
+    const bool verbose
+)
+{
+    return fvMeshTools::loadOrCreateMeshImpl
+    (
+        io,
+       &readHandler,
+        false,  // decompose (ignored)
+        verbose
+    );
+}
+
+
 void Foam::fvMeshTools::createDummyFvMeshFiles
 (
     const objectRegistry& mesh,
diff --git a/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.H b/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.H
index bdf80b80cf..2c899adb9b 100644
--- a/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.H
+++ b/src/finiteVolume/fvMesh/fvMeshTools/fvMeshTools.H
@@ -96,6 +96,21 @@ class fvMeshTools
     // Remove trailing patches
     static void trimPatches(fvMesh&, const label nPatches);
 
+    //- Read mesh if available, or create empty mesh with non-proc as per
+    //- proc0 mesh.
+    //  Has two modes of operation.
+    //  If the readHandler is non-nullptr, use it to decide on availability.
+    //  Requires:
+    //  - all processors to have all patches (and in same order).
+    //  - io.instance() set to facesInstance
+    static autoPtr<fvMesh> loadOrCreateMeshImpl
+    (
+        const IOobject& io,
+        refPtr<fileOperation>* readHandlerPtr,  // Can be nullptr
+        const bool decompose,  // Only used if readHandlerPtr == nullptr
+        const bool verbose = false
+    );
+
 
 public:
 
@@ -170,6 +185,18 @@ public:
         const bool verbose = false
     );
 
+    //- Read mesh if available, or create empty mesh with non-proc as per
+    //- proc0 mesh.
+    //  The fileHandler is specified on processors that have a processor mesh.
+    //  Generates empty mesh on other processors.
+    static autoPtr<fvMesh> loadOrCreateMesh
+    (
+        const IOobject& io,
+        //! Non-null reference if a mesh exists on given processor
+        refPtr<fileOperation>& readHandler,
+        const bool verbose = true
+    );
+
     //- Create additional fvSchemes/fvSolution files
     static void createDummyFvMeshFiles
     (