diff --git a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
index 2aaede0be8..bb34931e9a 100644
--- a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
+++ b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,8 @@ Original Authors
 #include "edgeFields.H"
 #include "processorFaPatch.H"
 #include "foamVtkIndPatchWriter.H"
+#include "foamVtkLineWriter.H"
+#include "faMeshTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,11 +68,27 @@ int main(int argc, char *argv[])
         "Write mesh as a vtp (vtk) file for display or debugging"
     );
 
+    argList::addOption
+    (
+        "geometryOrder",
+        "N",
+        "Test different geometry order - experimental!!",
+        true  // Advanced option
+    );
+
     #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedPolyMesh.H"
 
+    int geometryOrder(1);
+    if (args.readIfPresent("geometryOrder", geometryOrder))
+    {
+        Info<< "Setting faMesh::geometryOrder = " << geometryOrder << nl
+            << "(experimental)" << nl << endl;
+
+        faMesh::geometryOrder(geometryOrder);
+    }
 
     // Create
     faMesh aMesh(mesh);
diff --git a/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H b/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
index 1df8b4d312..c6edf85a3e 100644
--- a/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
+++ b/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -16,7 +16,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    // finiteArea
+    // finiteArea - faces
     vtk::uindirectPatchWriter writer
     (
         aMesh.patch(),
@@ -29,9 +29,23 @@ Description
 
     writer.writeGeometry();
 
+    globalIndex procAddr(aMesh.nFaces());
+    labelList cellIDs;
+
+    if (Pstream::master())
+    {
+        cellIDs.resize(procAddr.totalSize());
+        for (const labelRange& range : procAddr.ranges())
+        {
+            auto slice = cellIDs.slice(range);
+            slice = identity(range);
+        }
+    }
+
     // CellData
-    writer.beginCellData(3);
+    writer.beginCellData(4);
     writer.writeProcIDs();
+    writer.write("cellID", cellIDs);
     writer.write("area", aMesh.S().field());
     writer.write("normal", aMesh.faceAreaNormals());
 
@@ -43,5 +57,79 @@ Description
         << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
 }
 
+{
+    // finiteArea - edges
+    vtk::lineWriter writer
+    (
+        aMesh.points(),
+        aMesh.edges(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea-edges"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // CellData
+    writer.beginCellData(4);
+    writer.writeProcIDs();
+    {
+        // Use primitive patch order
+        Field<scalar> fld
+        (
+            faMeshTools::flattenEdgeField(aMesh.magLe(), true)
+        );
+        writer.write("magLe", fld);
+    }
+
+    // PointData
+    writer.beginPointData(1);
+    writer.write("normal", aMesh.pointAreaNormals());
+
+    Info<< nl
+        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+}
+
+{
+    // finiteArea - edgeCentres
+    // (no other convenient way to display vectors on the edges)
+    vtk::lineWriter writer
+    (
+        aMesh.edgeCentres(),
+        edgeList::null(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea-edgesCentres"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // PointData
+    writer.beginPointData(4);
+    {
+        // Use primitive patch order
+        Field<vector> fld
+        (
+            faMeshTools::flattenEdgeField(aMesh.Le(), true)
+        );
+        writer.write("Le", fld);
+    }
+    {
+        // Use primitive patch order
+        Field<vector> fld
+        (
+            faMeshTools::flattenEdgeField(aMesh.edgeAreaNormals(), true)
+        );
+        writer.write("normal", fld);
+    }
+
+    Info<< nl
+        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+}
+
 
 // ************************************************************************* //
diff --git a/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
index 256d2b1371..e3d255d961 100644
--- a/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
+++ b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -20,24 +20,7 @@ Description
     const label nNonProcessor = patches.nNonProcessor();
     const label nPatches = patches.size();
 
-    Info<< "----------------" << nl
-        << "Mesh Information" << nl
-        << "----------------" << nl
-        << "  " << "boundingBox: " << boundBox(aMesh.points()) << nl;
-
-    Info<< "  Number of points: "
-        << returnReduce(aMesh.nPoints(), sumOp<label>()) << nl
-        << "  Number of faces: "
-        << returnReduce(aMesh.nFaces(), sumOp<label>()) << nl;
-
-    Info<< "  Number of edges: "
-        << returnReduce(aMesh.nEdges(), sumOp<label>()) << nl
-        << "  Number of internal edges: "
-        << returnReduce(aMesh.nInternalEdges(), sumOp<label>()) << nl;
-
-
-    label nProcEdges = 0;
-
+    label nLocalProcEdges = 0;
     if (Pstream::parRun())
     {
         for (const faPatch& fap : patches)
@@ -46,22 +29,84 @@ Description
 
             if (cpp)
             {
-                nProcEdges += fap.nEdges();
+                nLocalProcEdges += fap.nEdges();
             }
         }
     }
 
-    const label nBoundEdges = aMesh.nBoundaryEdges() - nProcEdges;
+    const labelList nFaces
+    (
+        UPstream::listGatherValues<label>(aMesh.nFaces())
+    );
+    const labelList nPoints
+    (
+        UPstream::listGatherValues<label>(aMesh.nPoints())
+    );
+    const labelList nEdges
+    (
+        UPstream::listGatherValues<label>(aMesh.nEdges())
+    );
+    const labelList nIntEdges
+    (
+        UPstream::listGatherValues<label>(aMesh.nInternalEdges())
+    );
 
-    Info<< "  Number of boundary edges: "
-        << returnReduce(nBoundEdges - nProcEdges, sumOp<label>()) << nl;
+    // The "real" (non-processor) boundary edges
+    const labelList nBndEdges
+    (
+        UPstream::listGatherValues<label>
+        (
+            aMesh.nBoundaryEdges() - nLocalProcEdges
+        )
+    );
+    const labelList nProcEdges
+    (
+        UPstream::listGatherValues<label>(nLocalProcEdges)
+    );
 
-    if (Pstream::parRun())
+
+    // Format output as
+    //  Number of faces: ...
+    //         per-proc: (...)
+
+    const auto reporter =
+        [&](const char* tag, const labelList& list)
+        {
+            Info<< "  Number of " << tag << ": " << sum(list) << nl;
+            if (Pstream::parRun())
+            {
+                int padding = static_cast<int>
+                (
+                    // strlen("  Number of ") - strlen("per-proc")
+                    (12 - 8)
+                  + strlen(tag)
+                );
+
+                do { Info<< ' '; } while (--padding > 0);
+
+                Info<< "per-proc: " << flatOutput(list) << nl;
+            }
+        };
+
+
+    Info<< "----------------" << nl
+        << "Mesh Information" << nl
+        << "----------------" << nl
+        << "  " << "boundingBox: " << boundBox(aMesh.points()) << nl;
+
+    if (Pstream::master())
     {
-        Info<< "  Number of processor edges: "
-            << returnReduce(nProcEdges, sumOp<label>()) << nl;
-    }
+        reporter("faces", nFaces);
+        reporter("points", nPoints);
+        reporter("edges", nEdges);
+        reporter("internal edges", nIntEdges);
+        reporter("boundary edges", nBndEdges);
 
+        if (Pstream::parRun())
+        {
+            reporter("processor edges", nProcEdges);
+        }
+    }
 
     Info<< "----------------" << nl
         << "Patches" << nl
diff --git a/src/finiteArea/Make/files b/src/finiteArea/Make/files
index 99d971d4db..7db25ed3f3 100644
--- a/src/finiteArea/Make/files
+++ b/src/finiteArea/Make/files
@@ -9,6 +9,10 @@ faMesh/faMeshBoundaryHalo.C
 faMesh/faBoundaryMesh/faBoundaryMesh.C
 faMesh/faBoundaryMesh/faBoundaryMeshEntries.C
 
+faMesh/faMeshSubset/faMeshSubset.C
+faMesh/faMeshTools/faMeshTools.C
+faMesh/faMeshTools/faMeshToolsProcAddr.C
+
 faPatches = faMesh/faPatches
 $(faPatches)/faPatch/faPatch.C
 $(faPatches)/faPatch/faPatchData.C
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
index 5e3af56280..d6787152de 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@@ -235,6 +235,22 @@ void Foam::faBoundaryMesh::calcGeometry()
 }
 
 
+Foam::UPtrList<const Foam::labelUList>
+Foam::faBoundaryMesh::edgeLabels() const
+{
+    const faPatchList& patches = *this;
+
+    UPtrList<const labelUList> list(patches.size());
+
+    forAll(list, patchi)
+    {
+        list.set(patchi, &patches[patchi].edgeLabels());
+    }
+
+    return list;
+}
+
+
 Foam::UPtrList<const Foam::labelUList>
 Foam::faBoundaryMesh::edgeFaces() const
 {
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
index 1f17979b57..d6cf8d2a75 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
@@ -132,6 +132,9 @@ public:
             return mesh_;
         }
 
+        //- Return a list of edgeLabels for each patch
+        UPtrList<const labelUList> edgeLabels() const;
+
         //- Return a list of edgeFaces for each patch
         UPtrList<const labelUList> edgeFaces() const;
 
diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C
index 32e767ee26..32de35d34a 100644
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@@ -337,10 +337,12 @@ Foam::faMesh::faMesh
         *this
     ),
     comm_(Pstream::worldComm),
+    curTimeIndex_(time().timeIndex()),
+
     patchPtr_(nullptr),
     bndConnectPtr_(nullptr),
     lduPtr_(nullptr),
-    curTimeIndex_(time().timeIndex()),
+
     SPtr_(nullptr),
     S0Ptr_(nullptr),
     S00Ptr_(nullptr),
@@ -439,10 +441,12 @@ Foam::faMesh::faMesh
         label(0)
     ),
     comm_(Pstream::worldComm),
+    curTimeIndex_(time().timeIndex()),
+
     patchPtr_(nullptr),
     bndConnectPtr_(nullptr),
     lduPtr_(nullptr),
-    curTimeIndex_(time().timeIndex()),
+
     SPtr_(nullptr),
     S0Ptr_(nullptr),
     S00Ptr_(nullptr),
diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H
index c6de146f0f..4d59b4ac79 100644
--- a/src/finiteArea/faMesh/faMesh.H
+++ b/src/finiteArea/faMesh/faMesh.H
@@ -209,37 +209,43 @@ class faMesh
         faBoundaryMesh boundary_;
 
 
-        // Primitive mesh data
+    // Primitive mesh data
 
-            //- Edges, addressing into local point list
-            edgeList edges_;
+        //- Edges, addressing into local point list
+        edgeList edges_;
 
-            //- Edge owner
-            labelList edgeOwner_;
+        //- Edge owner
+        labelList edgeOwner_;
 
-            //- Edge neighbour
-            labelList edgeNeighbour_;
+        //- Edge neighbour
+        labelList edgeNeighbour_;
 
 
-        // Primitive size data
+    // Primitive size data
 
-            //- Total number of points
-            mutable label nPoints_;
+        //- Total number of points
+        mutable label nPoints_;
 
-            //- Total number of edges
-            mutable label nEdges_;
+        //- Total number of edges
+        mutable label nEdges_;
 
-            //- Number of internal edges
-            mutable label nInternalEdges_;
+        //- Number of internal edges
+        mutable label nInternalEdges_;
 
-            //- Number of faces
-            mutable label nFaces_;
+        //- Number of faces
+        mutable label nFaces_;
 
 
-        // Communication support
+    // Communication support, updating
 
-            //- Communicator used for parallel communication
-            label comm_;
+        //- Communicator used for parallel communication
+        label comm_;
+
+        //- Current time index for motion
+        //  Note.  The whole mechanism will be replaced once the
+        //  dimensionedField is created and the dimensionedField
+        //  will take care of the old-time levels.
+        mutable label curTimeIndex_;
 
 
     // Demand-driven data
@@ -253,11 +259,8 @@ class faMesh
         //- Ldu addressing data
         mutable faMeshLduAddressing* lduPtr_;
 
-        //- Current time index for motion
-        //  Note.  The whole mechanism will be replaced once the
-        //  dimensionedField is created and the dimensionedField
-        //  will take care of the old-time levels.
-        mutable label curTimeIndex_;
+
+    // Geometric Data
 
         //- Face areas
         mutable DimensionedField<scalar, areaMesh>* SPtr_;
@@ -545,6 +548,22 @@ public:
 
     // Static Functions
 
+        //- Return the current geometry treatment (0: primitive, 1: standard)
+        //  A zero level is with restricted neighbour information
+        static int geometryOrder() noexcept
+        {
+            return geometryOrder_;
+        }
+
+        //- Set the preferred geometry treatment
+        //  \return the previous value
+        static int geometryOrder(const int order) noexcept
+        {
+            int old(geometryOrder_);
+            geometryOrder_ = order;
+            return old;
+        }
+
         //- Read construction from polyMesh if all files are available
         static autoPtr<faMesh> TryNew(const polyMesh& pMesh);
 
@@ -571,73 +590,79 @@ public:
         bool init(const bool doInit);
 
 
-        // Database
+    // Database
 
-            //- Return access to polyMesh
-            inline const polyMesh& mesh() const;
+        //- Return access to polyMesh
+        inline const polyMesh& mesh() const;
 
-            //- Interface to referenced polyMesh (similar to GeoMesh)
-            const polyMesh& operator()() const { return mesh(); }
+        //- Interface to referenced polyMesh (similar to GeoMesh)
+        const polyMesh& operator()() const { return mesh(); }
 
-            //- Return the local mesh directory (dbDir()/meshSubDir)
-            fileName meshDir() const;
+        //- Return the local mesh directory (dbDir()/meshSubDir)
+        fileName meshDir() const;
 
-            //- Return reference to time
-            const Time& time() const;
+        //- Return reference to time
+        const Time& time() const;
 
-            //- Return the current instance directory for points
-            //  Used in the construction of geometric mesh data dependent
-            //  on points
-            const fileName& pointsInstance() const;
+        //- Return the current instance directory for points
+        //  Used in the construction of geometric mesh data dependent
+        //  on points
+        const fileName& pointsInstance() const;
 
-            //- Return the current instance directory for faces
-            const fileName& facesInstance() const;
+        //- Return the current instance directory for faces
+        const fileName& facesInstance() const;
 
 
-        // Communication support
+    // Communication support
 
-            //- Return communicator used for parallel communication
-            inline label comm() const noexcept;
+        //- Return communicator used for parallel communication
+        inline label comm() const noexcept;
 
-            //- Return communicator used for parallel communication
-            inline label& comm() noexcept;
+        //- Return communicator used for parallel communication
+        inline label& comm() noexcept;
 
 
-            // Mesh size parameters
+    // Access: Mesh size parameters
 
-                //- Number of local mesh points
-                inline label nPoints() const noexcept;
+        //- Number of local mesh points
+        inline label nPoints() const noexcept;
 
-                //- Number of local mesh edges
-                inline label nEdges() const noexcept;
+        //- Number of local mesh edges
+        inline label nEdges() const noexcept;
 
-                //- Number of internal faces
-                inline label nInternalEdges() const noexcept;
+        //- Number of internal faces
+        inline label nInternalEdges() const noexcept;
 
-                //- Number of boundary edges (== nEdges - nInternalEdges)
-                inline label nBoundaryEdges() const noexcept;
+        //- Number of boundary edges (== nEdges - nInternalEdges)
+        inline label nBoundaryEdges() const noexcept;
 
-                //- Number of patch faces
-                inline label nFaces() const noexcept;
+        //- Number of patch faces
+        inline label nFaces() const noexcept;
 
 
-            // Primitive mesh data
+    // Access: Primitive mesh data
 
-                //- Return local patch points
-                inline const pointField& points() const;
+        //- Return local points
+        inline const pointField& points() const;
 
-                //- Return local patch edges with reordered boundary
-                inline const edgeList& edges() const noexcept;
+        //- Return local edges with reordered boundary
+        inline const edgeList& edges() const noexcept;
 
-                //- Return local patch faces
-                inline const faceList& faces() const;
+        //- Sub-list of local internal edges
+        inline const edgeList::subList internalEdges() const;
 
-                //- Edge owner addressing
-                inline const labelList& edgeOwner() const noexcept;
+        //- Return local faces
+        inline const faceList& faces() const;
 
-                //- Edge neighbour addressing
-                inline const labelList& edgeNeighbour() const noexcept;
+        //- Edge owner addressing
+        inline const labelList& edgeOwner() const noexcept;
 
+        //- Edge neighbour addressing
+        inline const labelList& edgeNeighbour() const noexcept;
+
+        //- True if the internalEdges use an ordering that does not
+        //- correspond 1-to-1 with the patch internalEdges.
+        inline bool hasInternalEdgeLabels() const noexcept;
 
 
         // Access
@@ -649,7 +674,7 @@ public:
             virtual const objectRegistry& thisDb() const;
 
             //- Name function is needed to disambiguate those inherited
-            //  from base classes
+            //- from base classes
             const word& name() const
             {
                 return thisDb().name();
@@ -688,7 +713,7 @@ public:
             }
 
             //- True if given edge label is internal to the mesh
-            bool isInternalEdge(const label edgeIndex) const
+            bool isInternalEdge(const label edgeIndex) const noexcept
             {
                 return (edgeIndex < nInternalEdges_);
             }
@@ -817,6 +842,9 @@ public:
             boolList& correctPatchPointNormals() const;
 
 
+        // Storage management
+
+
         //- Write mesh
         virtual bool write(const bool valid = true) const;
 
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index c6af8c0499..fe9016edee 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -37,6 +37,7 @@ License
 #include "scalarMatrices.H"
 #include "processorFaPatchFields.H"
 #include "emptyFaPatchFields.H"
+#include "triPointRef.H"
 
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
@@ -102,6 +103,121 @@ static inline vector areaInvDistSqrWeightedNormalDualEdge
     );
 }
 
+
+// Calculate transform tensor with reference vector (unitAxis1)
+// and direction vector (axis2).
+//
+// This is nearly identical to the meshTools axesRotation
+// with an E3_E1 transformation with the following exceptions:
+//
+// - axis1 (e3 == unitAxis1): is already normalized (unit vector)
+// - axis2 (e1 == dirn): no difference
+// - transformation is row-vectors, not column-vectors
+static inline tensor rotation_e3e1
+(
+    const vector& unitAxis1,
+    vector dirn
+)
+{
+    dirn.removeCollinear(unitAxis1);
+    dirn.normalise();
+
+    // Set row vectors
+    return tensor
+    (
+        dirn,
+        (unitAxis1^dirn),
+        unitAxis1
+    );
+}
+
+
+// Simple area-weighted normal calculation for the specified edge vector
+// and its owner/neighbour face centres (internal edges).
+//
+// Uses four triangles since adjacent faces may be non-planar
+// and/or the edge centre is skewed from the face centres.
+
+/*---------------------------------------------------------------------------*\
+ *          *           |
+ *         /|\          | Triangles (0,1) are on the owner side.
+ *        / | \         |
+ *       /  |  \        | Triangles (2,3) are on the neighbour side.
+ *      /  1|3  \       |
+ * own *----|----* nbr  | Boundary edges are the same, but without a neighbour
+ *      \  0|2  /       |
+ *       \  |  /        |
+ *        \ | /         |
+ *         \|/          |
+ *          *           |
+\*---------------------------------------------------------------------------*/
+
+static inline vector calcEdgeNormalFromFace
+(
+    const linePointRef& edgeVec,
+    const point& ownCentre,
+    const point& neiCentre
+)
+{
+    const scalar magEdge(edgeVec.mag());
+
+    if (magEdge < ROOTVSMALL)
+    {
+        return Zero;
+    }
+
+    const vector edgeCtr = edgeVec.centre();
+
+    vector edgeNorm
+    (
+        // owner
+        triPointRef(edgeCtr, ownCentre, edgeVec.first()).areaNormal()
+      + triPointRef(edgeCtr, edgeVec.last(), ownCentre).areaNormal()
+        // neighbour
+      + triPointRef(edgeCtr, edgeVec.first(), neiCentre).areaNormal()
+      + triPointRef(edgeCtr, neiCentre, edgeVec.last()).areaNormal()
+    );
+
+    // Requires a unit-vector (already tested its mag)
+    edgeNorm.removeCollinear(edgeVec.vec()/magEdge);
+    edgeNorm.normalise();
+
+    // Scale with the original edge length
+    return magEdge*edgeNorm;
+}
+
+
+// As above, but for boundary edgess (no neighbour)
+static inline vector calcEdgeNormalFromFace
+(
+    const linePointRef& edgeVec,
+    const point& ownCentre
+)
+{
+    const scalar magEdge(edgeVec.mag());
+
+    if (magEdge < ROOTVSMALL)
+    {
+        return Zero;
+    }
+
+    const vector edgeCtr = edgeVec.centre();
+
+    vector edgeNorm
+    (
+        // owner
+        triPointRef(edgeCtr, ownCentre, edgeVec.first()).areaNormal()
+      + triPointRef(edgeCtr, edgeVec.last(), ownCentre).areaNormal()
+    );
+
+    // Requires a unit-vector (already tested its mag)
+    edgeNorm.removeCollinear(edgeVec.vec()/magEdge);
+    edgeNorm.normalise();
+
+    // Scale with the original edge length
+    return magEdge*edgeNorm;
+}
+
 } // End namespace Foam
 
 
@@ -189,71 +305,187 @@ void Foam::faMesh::calcLe() const
 
     edgeVectorField& Le = *LePtr_;
 
+    const pointField& localPoints = points();
 
-    const pointField& pPoints = points();
-    const edgeList& pEdges = edges();
+    if (faMesh::geometryOrder() < 1)
+    {
+        // Simple (primitive) edge normal calculations.
+        // These are primarly designed to avoid any communication
+        // but are thus necessarily inconsistent across processor boundaries!
+
+        // Reasonable to assume that the volume mesh already has faceCentres
+        // eg, used magSf somewhere.
+        // Can use these instead of triggering our calcAreaCentres().
+
+        WarningInFunction
+            << "Using geometryOrder < 1 : "
+               "simplified edge area-normals for Le() calculation"
+            << endl;
+
+        UIndirectList<vector> fCentres(mesh().faceCentres(), faceLabels());
+
+        // Flat addressing
+        vectorField edgeNormals(nEdges_);
+
+        // Internal (edge normals)
+        for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
+        {
+            edgeNormals[edgei] =
+                calcEdgeNormalFromFace
+                (
+                    edges_[edgei].line(localPoints),
+                    fCentres[owner()[edgei]],
+                    fCentres[neighbour()[edgei]]
+                );
+        }
+
+        // Boundary (edge normals). Like above, but only has owner
+        for (label edgei = nInternalEdges_; edgei < nEdges_; ++edgei)
+        {
+            edgeNormals[edgei] =
+                calcEdgeNormalFromFace
+                (
+                    edges_[edgei].line(localPoints),
+                    fCentres[owner()[edgei]]
+                );
+        }
+
+
+        // Now use these edge normals for calculating Le
+
+        // Internal (edge vector)
+        {
+            vectorField& internalFld = Le.ref();
+            for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
+            {
+                vector& leVector = internalFld[edgei];
+
+                vector edgeVec = edges_[edgei].vec(localPoints);
+                const scalar magEdge(mag(edgeVec));
+
+                if (magEdge < ROOTVSMALL)
+                {
+                    // Too small
+                    leVector = Zero;
+                    continue;
+                }
+
+                const vector edgeCtr = edges_[edgei].centre(localPoints);
+                const vector& edgeNorm = edgeNormals[edgei];
+                const vector& ownCentre = fCentres[owner()[edgei]];
+
+                leVector = magEdge*normalised(edgeVec ^ edgeNorm);
+                leVector *= sign(leVector & (edgeCtr - ownCentre));
+            }
+        }
+
+        // Boundary (edge vector)
+
+        forAll(boundary(), patchi)
+        {
+            const labelUList& bndEdgeFaces = boundary()[patchi].edgeFaces();
+
+            const edgeList::subList bndEdges =
+                boundary()[patchi].patchSlice(edges_);
+
+            vectorField& patchLe = Le.boundaryFieldRef()[patchi];
+
+            forAll(patchLe, bndEdgei)
+            {
+                vector& leVector = patchLe[bndEdgei];
+                const label meshEdgei(boundary()[patchi].start() + bndEdgei);
+
+                vector edgeVec = bndEdges[bndEdgei].vec(localPoints);
+                const scalar magEdge(mag(edgeVec));
+
+                if (magEdge < ROOTVSMALL)
+                {
+                    // Too small
+                    leVector = Zero;
+                    continue;
+                }
+
+                const vector edgeCtr = bndEdges[bndEdgei].centre(localPoints);
+                const vector& edgeNorm = edgeNormals[meshEdgei];
+                const vector& ownCentre = fCentres[bndEdgeFaces[bndEdgei]];
+
+                leVector = magEdge*normalised(edgeVec ^ edgeNorm);
+                leVector *= sign(leVector & (edgeCtr - ownCentre));
+            }
+        }
+
+        // Done
+        return;
+    }
+
+
+    // Longer forms.
+    // Using edgeAreaNormals, which uses pointAreaNormals (communication!)
 
     const edgeVectorField& eCentres = edgeCentres();
     const areaVectorField& fCentres = areaCentres();
-
     const edgeVectorField& edgeNormals = edgeAreaNormals();
 
-    vectorField& leInternal = Le.ref();
-    const vectorField& edgeNormalsInternal = edgeNormals.internalField();
-    const vectorField& fCentresInternal = fCentres.internalField();
-    const vectorField& eCentresInternal = eCentres.internalField();
-    const scalarField& magLeInternal = magLe().internalField();
+    // Internal (edge vector)
 
-    forAll(leInternal, edgeI)
     {
-        leInternal[edgeI] =
-            pEdges[edgeI].vec(pPoints) ^ edgeNormalsInternal[edgeI];
+        vectorField& internalFld = Le.ref();
+        for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
+        {
+            vector& leVector = internalFld[edgei];
 
-        leInternal[edgeI] *=
-          - sign
-            (
-                leInternal[edgeI] &
-                (
-                    fCentresInternal[owner()[edgeI]]
-                  - eCentresInternal[edgeI]
-                )
-            );
+            vector edgeVec = edges_[edgei].vec(localPoints);
+            const scalar magEdge(mag(edgeVec));
 
-        leInternal[edgeI] *=
-            magLeInternal[edgeI]/mag(leInternal[edgeI]);
+            if (magEdge < ROOTVSMALL)
+            {
+                // Too small
+                leVector = Zero;
+                continue;
+            }
+
+            const vector& edgeCtr = eCentres[edgei];
+            const vector& edgeNorm = edgeNormals[edgei];
+            const vector& ownCentre = fCentres[owner()[edgei]];
+
+            leVector = magEdge*normalised(edgeVec ^ edgeNorm);
+            leVector *= sign(leVector & (edgeCtr - ownCentre));
+        }
     }
 
-    forAll(boundary(), patchI)
+    forAll(boundary(), patchi)
     {
-        const labelUList& bndEdgeFaces = boundary()[patchI].edgeFaces();
+        const labelUList& bndEdgeFaces = boundary()[patchi].edgeFaces();
 
         const edgeList::subList bndEdges =
-            boundary()[patchI].patchSlice(pEdges);
+            boundary()[patchi].patchSlice(edges_);
 
         const vectorField& bndEdgeNormals =
-            edgeNormals.boundaryField()[patchI];
+            edgeNormals.boundaryField()[patchi];
 
-        vectorField& patchLe = Le.boundaryFieldRef()[patchI];
-        const vectorField& patchECentres = eCentres.boundaryField()[patchI];
+        vectorField& patchLe = Le.boundaryFieldRef()[patchi];
+        const vectorField& patchECentres = eCentres.boundaryField()[patchi];
 
-        forAll(patchLe, edgeI)
+        forAll(patchLe, bndEdgei)
         {
-            patchLe[edgeI] =
-                bndEdges[edgeI].vec(pPoints) ^ bndEdgeNormals[edgeI];
+            vector& leVector = patchLe[bndEdgei];
 
-            patchLe[edgeI] *=
-              - sign
-                (
-                    patchLe[edgeI]&
-                    (
-                        fCentresInternal[bndEdgeFaces[edgeI]]
-                      - patchECentres[edgeI]
-                    )
-                );
+            vector edgeVec = bndEdges[bndEdgei].vec(localPoints);
+            const scalar magEdge(mag(edgeVec));
 
-            patchLe[edgeI] *=
-                magLe().boundaryField()[patchI][edgeI]
-                /mag(patchLe[edgeI]);
+            if (magEdge < ROOTVSMALL)
+            {
+                // Too small
+                leVector = Zero;
+                continue;
+            }
+
+            const vector& edgeCtr = patchECentres[bndEdgei];
+            const vector& edgeNorm = bndEdgeNormals[bndEdgei];
+            const vector& ownCentre = fCentres[bndEdgeFaces[bndEdgei]];
+
+            leVector = magEdge*normalised(edgeVec ^ edgeNorm);
+            leVector *= sign(leVector & (edgeCtr - ownCentre));
         }
     }
 }
@@ -289,23 +521,30 @@ void Foam::faMesh::calcMagLe() const
 
     const pointField& localPoints = points();
 
-    label edgei = 0;
-    for (const edge& e : patch().internalEdges())
+    // Internal (edge length)
     {
-        magLe.ref()[edgei] = e.mag(localPoints);
-        ++edgei;
+        auto iter = magLe.primitiveFieldRef().begin();
+
+        for (const edge& e : internalEdges())
+        {
+            *iter = e.mag(localPoints);
+            ++iter;
+        }
     }
 
-
-    forAll(boundary(), patchI)
+    // Boundary (edge length)
     {
-        const edgeList::subList patchEdges =
-            boundary()[patchI].patchSlice(edges());
+        auto& bfld = magLe.boundaryFieldRef();
 
-        forAll(patchEdges, edgeI)
+        forAll(boundary(), patchi)
         {
-            magLe.boundaryFieldRef()[patchI][edgeI] =
-                patchEdges[edgeI].mag(localPoints);
+            auto iter = bfld[patchi].begin();
+
+            for (const edge& e : boundary()[patchi].patchSlice(edges_))
+            {
+                *iter = e.mag(localPoints);
+                ++iter;
+            }
         }
     }
 }
@@ -343,21 +582,25 @@ void Foam::faMesh::calcAreaCentres() const
     const faceList& localFaces = faces();
 
     // Internal (face centres)
-    forAll(localFaces, faceI)
+    // Could also obtain from volume mesh faceCentres()
+    forAll(localFaces, facei)
     {
-        centres.ref()[faceI] = localFaces[faceI].centre(localPoints);
+        centres.ref()[facei] = localFaces[facei].centre(localPoints);
     }
 
     // Boundary (edge centres)
-    forAll(boundary(), patchI)
     {
-        const edgeList::subList patchEdges =
-            boundary()[patchI].patchSlice(edges());
+        auto& bfld = centres.boundaryFieldRef();
 
-        forAll(patchEdges, edgeI)
+        forAll(boundary(), patchi)
         {
-            centres.boundaryFieldRef()[patchI][edgeI] =
-                patchEdges[edgeI].centre(localPoints);
+            auto iter = bfld[patchi].begin();
+
+            for (const edge& e : boundary()[patchi].patchSlice(edges_))
+            {
+                *iter = e.centre(localPoints);
+                ++iter;
+            }
         }
     }
 }
@@ -389,28 +632,34 @@ void Foam::faMesh::calcEdgeCentres() const
             dimLength
         );
 
-    edgeVectorField& edgeCentres = *edgeCentresPtr_;
+    edgeVectorField& centres = *edgeCentresPtr_;
 
     const pointField& localPoints = points();
 
-    // Internal edges
-    label edgei = 0;
-    for (const edge& e : patch().internalEdges())
+    // Internal (edge centres)
     {
-        edgeCentres.ref()[edgei] = e.centre(localPoints);
-        ++edgei;
+        auto iter = centres.primitiveFieldRef().begin();
+
+        for (const edge& e : internalEdges())
+        {
+            *iter = e.centre(localPoints);
+            ++iter;
+        }
     }
 
-    // Boundary edges
-    forAll(boundary(), patchI)
+    // Boundary (edge centres)
     {
-        const edgeList::subList patchEdges =
-            boundary()[patchI].patchSlice(edges());
+        auto& bfld = centres.boundaryFieldRef();
 
-        forAll(patchEdges, edgeI)
+        forAll(boundary(), patchi)
         {
-            edgeCentres.boundaryFieldRef()[patchI][edgeI] =
-                patchEdges[edgeI].centre(localPoints);
+            auto iter = bfld[patchi].begin();
+
+            for (const edge& e : boundary()[patchi].patchSlice(edges_))
+            {
+                *iter = e.centre(localPoints);
+                ++iter;
+            }
         }
     }
 }
@@ -446,9 +695,10 @@ void Foam::faMesh::calcS() const
     const pointField& localPoints = points();
     const faceList& localFaces = faces();
 
-    forAll(S, faceI)
+    // Could also obtain from volume mesh faceAreas()
+    forAll(S, facei)
     {
-        S[faceI] = localFaces[faceI].mag(localPoints);
+        S[facei] = localFaces[facei].mag(localPoints);
     }
 }
 
@@ -485,6 +735,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
     const faceList& localFaces = faces();
 
     // Internal (faces)
+    // Could also obtain from volume mesh Sf() + normalise
     vectorField& nInternal = faceNormals.ref();
     forAll(localFaces, faceI)
     {
@@ -493,8 +744,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
 
     // Boundary - copy from edges
 
-    const edgeVectorField::Boundary& edgeNormalsBoundary =
-        edgeAreaNormals().boundaryField();
+    const auto& edgeNormalsBoundary = edgeAreaNormals().boundaryField();
 
     forAll(boundary(), patchI)
     {
@@ -528,46 +778,50 @@ void Foam::faMesh::calcEdgeAreaNormals() const
             *this,
             dimless
         );
-
     edgeVectorField& edgeAreaNormals = *edgeAreaNormalsPtr_;
 
 
-    // Point area normals
+    // Starting from point area normals
     const vectorField& pointNormals = pointAreaNormals();
 
+
+    // Internal edges
     forAll(edgeAreaNormals.internalField(), edgei)
     {
-        const edge& e = edges()[edgei];
+        const edge& e = edges_[edgei];
         const vector edgeVec = e.unitVec(points());
 
-        vector& n = edgeAreaNormals.ref()[edgei];
+        vector& edgeNorm = edgeAreaNormals.ref()[edgei];
 
-        n = (pointNormals[e.first()] + pointNormals[e.second()]);
+        // Average of both ends
+        edgeNorm = (pointNormals[e.first()] + pointNormals[e.second()]);
 
-        n -= edgeVec*(edgeVec & n);
-
-        n.normalise();
+        edgeNorm.removeCollinear(edgeVec);
+        edgeNorm.normalise();
     }
 
+    // Boundary
+    auto& bfld = edgeAreaNormals.boundaryFieldRef();
+
     forAll(boundary(), patchi)
     {
+        auto& pfld = bfld[patchi];
+
         const edgeList::subList patchEdges =
-            boundary()[patchi].patchSlice(edges());
+            boundary()[patchi].patchSlice(edges_);
 
-        vectorField& edgeNorms = edgeAreaNormals.boundaryFieldRef()[patchi];
-
-        forAll(patchEdges, edgei)
+        forAll(patchEdges, bndEdgei)
         {
-            const edge& e = patchEdges[edgei];
+            const edge& e = patchEdges[bndEdgei];
             const vector edgeVec = e.unitVec(points());
 
-            vector& n = edgeNorms[edgei];
+            vector& edgeNorm = pfld[bndEdgei];
 
-            n = (pointNormals[e.first()] + pointNormals[e.second()]);
+            // Average of both ends
+            edgeNorm = (pointNormals[e.first()] + pointNormals[e.second()]);
 
-            n -= edgeVec*(edgeVec & n);
-
-            n.normalise();
+            edgeNorm.removeCollinear(edgeVec);
+            edgeNorm.normalise();
         }
     }
 }
@@ -624,9 +878,14 @@ void Foam::faMesh::calcEdgeTransformTensors() const
             << abort(FatalError);
     }
 
-    edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges());
-    FieldField<Field, tensor>& edgeTransformTensors =
-        *edgeTransformTensorsPtr_;
+    edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges_);
+    auto& edgeTransformTensors = *edgeTransformTensorsPtr_;
+
+    // Initialize all transformation tensors
+    for (label edgei = 0; edgei < nEdges_; ++edgei)
+    {
+        edgeTransformTensors.set(edgei, new Field<tensor>(3, tensor::I));
+    }
 
     const areaVectorField& Nf = faceAreaNormals();
     const areaVectorField& Cf = areaCentres();
@@ -635,193 +894,129 @@ void Foam::faMesh::calcEdgeTransformTensors() const
     const edgeVectorField& Ce = edgeCentres();
 
     // Internal edges transformation tensors
-    for (label edgeI=0; edgeI<nInternalEdges(); ++edgeI)
+    for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
     {
-        edgeTransformTensors.set(edgeI, new Field<tensor>(3, I));
+        const label ownFacei = owner()[edgei];
+        const label neiFacei = neighbour()[edgei];
+        auto& tensors = edgeTransformTensors[edgei];
 
-        vector E = Ce.internalField()[edgeI];
+        vector edgeCtr = Ce.internalField()[edgei];
 
         if (skew())
         {
-            E -= skewCorrectionVectors().internalField()[edgeI];
+            edgeCtr -= skewCorrectionVectors().internalField()[edgei];
         }
 
         // Edge transformation tensor
-        vector il = E - Cf.internalField()[owner()[edgeI]];
-
-        il -= Ne.internalField()[edgeI]
-            *(Ne.internalField()[edgeI]&il);
-
-        il /= mag(il);
-
-        vector kl = Ne.internalField()[edgeI];
-        vector jl = kl^il;
-
-        edgeTransformTensors[edgeI][0] =
-            tensor
+        tensors[0] =
+            rotation_e3e1
             (
-                il.x(), il.y(), il.z(),
-                jl.x(), jl.y(), jl.z(),
-                kl.x(), kl.y(), kl.z()
+                Ne.internalField()[edgei],
+                (edgeCtr - Cf.internalField()[ownFacei])
             );
 
         // Owner transformation tensor
-        il = E - Cf.internalField()[owner()[edgeI]];
-
-        il -= Nf.internalField()[owner()[edgeI]]
-            *(Nf.internalField()[owner()[edgeI]]&il);
-
-        il /= mag(il);
-
-        kl = Nf.internalField()[owner()[edgeI]];
-        jl = kl^il;
-
-        edgeTransformTensors[edgeI][1] =
-            tensor
+        tensors[1] =
+            rotation_e3e1
             (
-                il.x(), il.y(), il.z(),
-                jl.x(), jl.y(), jl.z(),
-                kl.x(), kl.y(), kl.z()
+                Nf.internalField()[ownFacei],
+                (edgeCtr - Cf.internalField()[ownFacei])
             );
 
         // Neighbour transformation tensor
-        il = Cf.internalField()[neighbour()[edgeI]] - E;
-
-        il -= Nf.internalField()[neighbour()[edgeI]]
-            *(Nf.internalField()[neighbour()[edgeI]]&il);
-
-        il /= mag(il);
-
-        kl = Nf.internalField()[neighbour()[edgeI]];
-        jl = kl^il;
-
-        edgeTransformTensors[edgeI][2] =
-            tensor
+        tensors[2] =
+            rotation_e3e1
             (
-                il.x(), il.y(), il.z(),
-                jl.x(), jl.y(), jl.z(),
-                kl.x(), kl.y(), kl.z()
+                Nf.internalField()[neiFacei],
+                (Cf.internalField()[neiFacei] - edgeCtr)
             );
     }
 
     // Boundary edges transformation tensors
-    forAll(boundary(), patchI)
+    forAll(boundary(), patchi)
     {
-        if (boundary()[patchI].coupled())
+        const labelUList& edgeFaces = boundary()[patchi].edgeFaces();
+
+        if (boundary()[patchi].coupled())
         {
-            const labelUList& edgeFaces =
-                boundary()[patchI].edgeFaces();
+            vectorField nbrCf(Cf.boundaryField()[patchi].patchNeighbourField());
+            vectorField nbrNf(Nf.boundaryField()[patchi].patchNeighbourField());
 
-            vectorField ngbCf(Cf.boundaryField()[patchI].patchNeighbourField());
-
-            vectorField ngbNf(Nf.boundaryField()[patchI].patchNeighbourField());
-
-            forAll(edgeFaces, edgeI)
+            forAll(edgeFaces, bndEdgei)
             {
-                edgeTransformTensors.set
-                (
-                    boundary()[patchI].start() + edgeI,
-                    new Field<tensor>(3, I)
-                );
+                const label ownFacei = edgeFaces[bndEdgei];
+                const label meshEdgei(boundary()[patchi].start() + bndEdgei);
 
-                vector E = Ce.boundaryField()[patchI][edgeI];
+                auto& tensors = edgeTransformTensors[meshEdgei];
+
+                vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
 
                 if (skew())
                 {
-                    E -= skewCorrectionVectors()
-                        .boundaryField()[patchI][edgeI];
+                    edgeCtr -= skewCorrectionVectors()
+                        .boundaryField()[patchi][bndEdgei];
                 }
 
                 // Edge transformation tensor
-                vector il = E - Cf.internalField()[edgeFaces[edgeI]];
-
-                il -= Ne.boundaryField()[patchI][edgeI]
-                   *(Ne.boundaryField()[patchI][edgeI]&il);
-
-                il /= mag(il);
-
-                vector kl = Ne.boundaryField()[patchI][edgeI];
-                vector jl = kl^il;
-
-                edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
-                    tensor(il, jl, kl);
+                tensors[0] =
+                    rotation_e3e1
+                    (
+                        Ne.boundaryField()[patchi][bndEdgei],
+                        (edgeCtr - Cf.internalField()[ownFacei])
+                    );
 
                 // Owner transformation tensor
-                il = E - Cf.internalField()[edgeFaces[edgeI]];
-
-                il -= Nf.internalField()[edgeFaces[edgeI]]
-                   *(Nf.internalField()[edgeFaces[edgeI]]&il);
-
-                il /= mag(il);
-
-                kl = Nf.internalField()[edgeFaces[edgeI]];
-                jl = kl^il;
-
-                edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
-                    tensor(il, jl, kl);
+                tensors[1] =
+                    rotation_e3e1
+                    (
+                        Nf.internalField()[ownFacei],
+                        (edgeCtr - Cf.internalField()[ownFacei])
+                    );
 
                 // Neighbour transformation tensor
-                il = ngbCf[edgeI] - E;
-
-                il -= ngbNf[edgeI]*(ngbNf[edgeI]&il);
-
-                il /= mag(il);
-
-                kl = ngbNf[edgeI];
-
-                jl = kl^il;
-
-                edgeTransformTensors[boundary()[patchI].start() + edgeI][2] =
-                    tensor(il, jl, kl);
+                tensors[2] =
+                    rotation_e3e1
+                    (
+                        nbrNf[bndEdgei],
+                        (nbrCf[bndEdgei] - edgeCtr)
+                    );
             }
         }
         else
         {
-            const labelUList& edgeFaces = boundary()[patchI].edgeFaces();
-
-            forAll(edgeFaces, edgeI)
+            forAll(edgeFaces, bndEdgei)
             {
-                edgeTransformTensors.set
-                (
-                    boundary()[patchI].start() + edgeI,
-                    new Field<tensor>(3, I)
-                );
+                const label ownFacei = edgeFaces[bndEdgei];
+                const label meshEdgei(boundary()[patchi].start() + bndEdgei);
 
-                vector E = Ce.boundaryField()[patchI][edgeI];
+                auto& tensors = edgeTransformTensors[meshEdgei];
+
+                vector edgeCtr = Ce.boundaryField()[patchi][bndEdgei];
 
                 if (skew())
                 {
-                    E -= skewCorrectionVectors()
-                        .boundaryField()[patchI][edgeI];
+                    edgeCtr -= skewCorrectionVectors()
+                        .boundaryField()[patchi][bndEdgei];
                 }
 
                 // Edge transformation tensor
-                vector il = E - Cf.internalField()[edgeFaces[edgeI]];
-
-                il -= Ne.boundaryField()[patchI][edgeI]
-                   *(Ne.boundaryField()[patchI][edgeI]&il);
-
-                il /= mag(il);
-
-                vector kl = Ne.boundaryField()[patchI][edgeI];
-                vector jl = kl^il;
-
-                edgeTransformTensors[boundary()[patchI].start() + edgeI][0] =
-                    tensor(il, jl, kl);
+                tensors[0] =
+                    rotation_e3e1
+                    (
+                        Ne.boundaryField()[patchi][bndEdgei],
+                        (edgeCtr - Cf.internalField()[ownFacei])
+                    );
 
                 // Owner transformation tensor
-                il = E - Cf.internalField()[edgeFaces[edgeI]];
+                tensors[1] =
+                    rotation_e3e1
+                    (
+                        Nf.internalField()[ownFacei],
+                        (edgeCtr - Cf.internalField()[ownFacei])
+                    );
 
-                il -= Nf.internalField()[edgeFaces[edgeI]]
-                   *(Nf.internalField()[edgeFaces[edgeI]]&il);
-
-                il /= mag(il);
-
-                kl = Nf.internalField()[edgeFaces[edgeI]];
-                jl = kl^il;
-
-                edgeTransformTensors[boundary()[patchI].start() + edgeI][1] =
-                    tensor(il, jl, kl);
+                // Neighbour transformation tensor
+                tensors[2] = tensor::I;
             }
         }
     }
@@ -882,7 +1077,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
     DebugInFunction
         << "Calculating pointAreaNormals : face dual method" << endl;
 
-    result.resize(nPoints());
+    result.resize_nocopy(nPoints());
     result = Zero;
 
     const pointField& points = patch().localPoints();
@@ -945,7 +1140,7 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
 
             for (const label pti : patchPoints)
             {
-                result[pti] -= N*(N&result[pti]);
+                result[pti].removeCollinear(N);
             }
 
             // Axis point correction
@@ -1124,17 +1319,13 @@ void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
         forAllConstIters(fpNormals, iter)
         {
             const label pointi = iter.key();
-            vector fpnorm = iter.val();
+            vector fpnorm = normalised(iter.val());
 
-            fpnorm.normalise();
-            result[pointi] -= fpnorm*(fpnorm & result[pointi]);
+            result[pointi].removeCollinear(fpnorm);
         }
     }
 
-    for (vector& n : result)
-    {
-        n.normalise();
-    }
+    result.normalise();
 }
 
 
@@ -1469,7 +1660,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
                 const vector& origin = points[curPoint];
                 const vector axis = normalised(result[curPoint]);
                 vector dir(allPoints[0] - points[curPoint]);
-                dir -= axis*(axis&dir);
+                dir.removeCollinear(axis);
                 dir.normalise();
 
                 const coordinateSystem cs("cs", origin, axis, dir);
@@ -1633,7 +1824,7 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
                 const vector& origin = points[curPoint];
                 const vector axis = normalised(result[curPoint]);
                 vector dir(allPoints[0] - points[curPoint]);
-                dir -= axis*(axis&dir);
+                dir.removeCollinear(axis);
                 dir.normalise();
 
                 coordinateSystem cs("cs", origin, axis, dir);
@@ -1710,22 +1901,19 @@ Foam::tmp<Foam::edgeScalarField> Foam::faMesh::edgeLengthCorrection() const
     DebugInFunction
         << "Calculating edge length correction" << endl;
 
-    tmp<edgeScalarField> tcorrection
+    auto tcorrection = tmp<edgeScalarField>::New
     (
-        new edgeScalarField
+        IOobject
         (
-            IOobject
-            (
-                "edgeLengthCorrection",
-                mesh().pointsInstance(),
-                meshSubDir,
-                mesh()
-            ),
-            *this,
-            dimless
-        )
+            "edgeLengthCorrection",
+            mesh().pointsInstance(),
+            meshSubDir,
+            mesh()
+        ),
+        *this,
+        dimless
     );
-    edgeScalarField& correction = tcorrection.ref();
+    auto& correction = tcorrection.ref();
 
     const vectorField& pointNormals = pointAreaNormals();
 
diff --git a/src/finiteArea/faMesh/faMeshI.H b/src/finiteArea/faMesh/faMeshI.H
index dc1b238c3d..dfdbc00f79 100644
--- a/src/finiteArea/faMesh/faMeshI.H
+++ b/src/finiteArea/faMesh/faMeshI.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,7 +73,7 @@ inline Foam::label Foam::faMesh::nInternalEdges() const noexcept
 
 inline Foam::label Foam::faMesh::nBoundaryEdges() const noexcept
 {
-    return nEdges_ - nInternalEdges_;
+    return (nEdges_ - nInternalEdges_);
 }
 
 
@@ -95,6 +95,12 @@ inline const Foam::edgeList& Foam::faMesh::edges() const noexcept
 }
 
 
+inline const Foam::edgeList::subList Foam::faMesh::internalEdges() const
+{
+    return edgeList::subList(edges_, nInternalEdges_);
+}
+
+
 inline const Foam::faceList& Foam::faMesh::faces() const
 {
     return patch().localFaces();
@@ -139,6 +145,12 @@ inline Foam::uindirectPrimitivePatch& Foam::faMesh::patch()
 }
 
 
+inline bool Foam::faMesh::hasInternalEdgeLabels() const noexcept
+{
+    return false;
+}
+
+
 inline const Foam::List<Foam::labelPair>&
 Foam::faMesh::boundaryConnections() const
 {
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.C b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.C
new file mode 100644
index 0000000000..5a088e3a23
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.C
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faMeshSubset.H"
+#include "boolList.H"
+#include "BitOps.H"
+#include "Pstream.H"
+#include "emptyFaPatch.H"
+#include "cyclicFaPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::word Foam::faMeshSubset::exposedPatchName("oldInternalEdges");
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::faMeshSubset::checkHasSubMesh() const
+{
+    if (!subMeshPtr_)
+    {
+        FatalErrorInFunction
+            << "Mesh is not subsetted!" << nl
+            << abort(FatalError);
+
+        return false;
+    }
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faMeshSubset::faMeshSubset(const faMesh& baseMesh)
+:
+    baseMesh_(baseMesh),
+    subMeshPtr_(nullptr),
+    edgeFlipMapPtr_(nullptr),
+    pointMap_(),
+    faceMap_(),
+    cellMap_(),
+    patchMap_()
+{}
+
+
+Foam::faMeshSubset::faMeshSubset(const faMesh& baseMesh, const Foam::zero)
+:
+    faMeshSubset(baseMesh)
+{
+    reset(Foam::zero{});
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faMeshSubset::clear()
+{
+    subMeshPtr_.reset(nullptr);
+    edgeFlipMapPtr_.reset(nullptr);
+
+    pointMap_.clear();
+    faceMap_.clear();
+    cellMap_.clear();
+    patchMap_.clear();
+}
+
+
+void Foam::faMeshSubset::reset()
+{
+    clear();
+}
+
+
+void Foam::faMeshSubset::reset(const Foam::zero)
+{
+    clear();
+
+    // Create zero-sized subMesh
+    subMeshPtr_.reset
+    (
+        new faMesh
+        (
+            baseMesh_.mesh(),   // The polyMesh
+            // IOobject
+            // (
+            //     baseMesh_.name(),
+            //     baseMesh_.time().timeName(),
+            //     baseMesh_.time(),
+            //     IOobject::READ_IF_PRESENT,  // Read fa* if present
+            //     IOobject::NO_WRITE
+            // ),
+            Foam::zero{}                    // zero-sized
+        )
+    );
+    auto& newSubMesh = subMeshPtr_();
+
+
+    // Clone non-processor patches
+    {
+        const faBoundaryMesh& oldBoundary = baseMesh_.boundary();
+        const faBoundaryMesh& newBoundary = newSubMesh.boundary();
+
+        faPatchList newPatches(oldBoundary.nNonProcessor());
+
+        patchMap_ = identity(newPatches.size());
+
+        forAll(newPatches, patchi)
+        {
+            newPatches.set
+            (
+                patchi,
+                oldBoundary[patchi].clone
+                (
+                    newBoundary,
+                    labelList(),   // edgeLabels
+                    patchi,
+                    oldBoundary[patchi].ngbPolyPatchIndex()
+                )
+            );
+        }
+
+        newSubMesh.addFaPatches(newPatches);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H
new file mode 100644
index 0000000000..c02ed434f1
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubset.H
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::faMeshSubset
+
+Description
+    Holds a reference to the original mesh (the baseMesh)
+    and optionally to a subset of that mesh (the subMesh)
+    with mapping lists for points, faces, and cells.
+
+Caution
+    Currently not really functional for subsetting beyond handling
+    a simple zero-sized mesh.
+
+SourceFiles
+    faMeshSubset.C
+    faMeshSubsetI.H
+    faMeshSubsetTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_faMeshSubset_H
+#define Foam_faMeshSubset_H
+
+#include "areaFaMesh.H"
+#include "edgeFaMesh.H"
+#include "GeometricField.H"
+#include "bitSet.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                        Class faMeshSubset Declaration
+\*---------------------------------------------------------------------------*/
+
+class faMeshSubset
+{
+    // Private Data
+
+        //- The base mesh to subset from
+        const faMesh& baseMesh_;
+
+        //- Demand-driven subset mesh (pointer)
+        autoPtr<faMesh> subMeshPtr_;
+
+        //- Optional edge mapping array with flip encoded (-1/+1)
+        mutable autoPtr<labelList> edgeFlipMapPtr_;
+
+        //- Point mapping array
+        labelList pointMap_;
+
+        //- Face mapping array
+        labelList faceMap_;
+
+        //- Cell mapping array
+        labelList cellMap_;
+
+        //- Patch mapping array
+        labelList patchMap_;
+
+
+    // Private Member Functions
+
+        //- Calculate face flip map
+        void calcEdgeFlipMap() const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- FatalError if subset has not been performed
+        bool checkHasSubMesh() const;
+
+        //- No copy construct
+        faMeshSubset(const faMeshSubset&) = delete;
+
+        //- No copy assignment
+        void operator=(const faMeshSubset&) = delete;
+
+
+public:
+
+    // Static Data Members
+
+        //- Name for exposed internal edges (default: oldInternalEdges)
+        static word exposedPatchName;
+
+
+    // Constructors
+
+        //- Construct using the entire mesh (no subset)
+        explicit faMeshSubset(const faMesh& baseMesh);
+
+        //- Construct a zero-sized subset mesh, non-processor patches only
+        faMeshSubset(const faMesh& baseMesh, const Foam::zero);
+
+
+    // Member Functions
+
+    // Access
+
+        //- Original mesh
+        inline const faMesh& baseMesh() const noexcept;
+
+        //- Return baseMesh or subMesh, depending on the current state.
+        inline const faMesh& mesh() const noexcept;
+
+        //- Have subMesh?
+        inline bool hasSubMesh() const noexcept;
+
+        //- Return reference to subset mesh
+        inline const faMesh& subMesh() const;
+
+        //- Return reference to subset mesh
+        inline faMesh& subMesh();
+
+        //- Return point map
+        inline const labelList& pointMap() const;
+
+        //- Return face map
+        inline const labelList& faceMap() const;
+
+        //- Return edge map with sign to encode flipped edges
+        inline const labelList& edgeFlipMap() const;
+
+        //- Return cell map
+        inline const labelList& cellMap() const;
+
+        //- Return patch map
+        inline const labelList& patchMap() const;
+
+
+    // Edit
+
+        //- Reset subMesh and all maps
+        void clear();
+
+        //- Reset subMesh and all maps. Same as clear()
+        void reset();
+
+        //- Reset to a zero-sized subset mesh, non-processor patches only
+        void reset(const Foam::zero);
+
+
+    // Field Mapping (static functions)
+
+        //- Map area field.
+        //  Optionally allow unmapped faces not to produce a warning
+        template<class Type>
+        static tmp<GeometricField<Type, faPatchField, areaMesh>>
+        interpolate
+        (
+            const GeometricField<Type, faPatchField, areaMesh>&,
+            const faMesh& sMesh,
+            const bool allowUnmapped = false
+        );
+
+
+    // Field Mapping
+
+        //- Map area field.
+        //  Optionally allow unmapped faces not to produce a warning
+        template<class Type>
+        tmp<GeometricField<Type, faPatchField, areaMesh>>
+        interpolate
+        (
+            const GeometricField<Type, faPatchField, areaMesh>&,
+            const bool allowUnmapped = false
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "faMeshSubsetI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "faMeshSubsetTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetI.H b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetI.H
new file mode 100644
index 0000000000..8a1d67cd39
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetI.H
@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::faMesh& Foam::faMeshSubset::baseMesh() const noexcept
+{
+    return baseMesh_;
+}
+
+
+inline const Foam::faMesh& Foam::faMeshSubset::mesh() const noexcept
+{
+    return subMeshPtr_ ? *subMeshPtr_ : baseMesh_;
+}
+
+
+inline bool Foam::faMeshSubset::hasSubMesh() const noexcept
+{
+    return bool(subMeshPtr_);
+}
+
+
+inline const Foam::faMesh& Foam::faMeshSubset::subMesh() const
+{
+    checkHasSubMesh();
+
+    return *subMeshPtr_;
+}
+
+
+inline Foam::faMesh& Foam::faMeshSubset::subMesh()
+{
+    checkHasSubMesh();
+
+    return *subMeshPtr_;
+}
+
+
+inline const Foam::labelList& Foam::faMeshSubset::pointMap() const
+{
+    checkHasSubMesh();
+
+    return pointMap_;
+}
+
+
+inline const Foam::labelList& Foam::faMeshSubset::faceMap() const
+{
+    checkHasSubMesh();
+
+    return faceMap_;
+}
+
+
+inline const Foam::labelList& Foam::faMeshSubset::edgeFlipMap() const
+{
+    if (!edgeFlipMapPtr_)
+    {
+        // calcEdgeFlipMap();
+    }
+
+    return *edgeFlipMapPtr_;
+}
+
+
+inline const Foam::labelList& Foam::faMeshSubset::cellMap() const
+{
+    checkHasSubMesh();
+
+    return cellMap_;
+}
+
+
+inline const Foam::labelList& Foam::faMeshSubset::patchMap() const
+{
+    checkHasSubMesh();
+
+    return patchMap_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C
new file mode 100644
index 0000000000..c245505ca7
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshSubset/faMeshSubsetTemplates.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faMeshSubset.H"
+#include "areaFaMesh.H"
+#include "edgeFaMesh.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "emptyFaPatchFields.H"
+#include "directFaPatchFieldMapper.H"
+#include "flipOp.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>
+>
+Foam::faMeshSubset::interpolate
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const faMesh& sMesh,
+    const bool allowUnmapped
+)
+{
+    // 1. Create the complete field with dummy patch fields
+    PtrList<faPatchField<Type>> patchFields(sMesh.boundary().size());
+
+    forAll(patchFields, patchi)
+    {
+        patchFields.set
+        (
+            patchi,
+            faPatchField<Type>::New
+            (
+                calculatedFaPatchField<Type>::typeName,
+                sMesh.boundary()[patchi],
+                DimensionedField<Type, areaMesh>::null()
+            )
+        );
+    }
+
+    auto tresult = tmp<GeometricField<Type, faPatchField, areaMesh>>::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(), cellMap),
+        patchFields
+    );
+    auto& result = tresult.ref();
+    result.oriented() = vf.oriented();
+
+
+    // 2. Change the faPatchFields 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);
+
+        // allowUnmapped : special mode for if we do not want to be
+        // warned for unmapped faces (e.g. from faMeshDistribute).
+        const bool hasUnmapped = mapper.hasUnmapped();
+        if (allowUnmapped)
+        {
+            mapper.hasUnmapped() = false;
+        }
+
+        bf.set
+        (
+            patchi,
+            faPatchField<Type>::New
+            (
+                vf.boundaryField()[patchi],
+                subPatch,
+                result(),
+                mapper
+            )
+        );
+
+        if (allowUnmapped && hasUnmapped)
+        {
+            // Set unmapped values to zeroGradient. This is the default
+            // action for unmapped faPatchFields. Note that this bypasses
+            // any special logic for handling unmapped faPatchFields but
+            // since this is only used inside faMeshDistribute ...
+
+            tmp<Field<Type>> tfld(bf[patchi].patchInternalField());
+            const Field<Type>& fld = tfld();
+
+            Field<Type> value(bf[patchi]);
+            forAll(directAddressing, i)
+            {
+                if (directAddressing[i] == -1)
+                {
+                    value[i] = fld[i];
+                }
+            }
+            bf[patchi].faPatchField<Type>::operator=(value);
+        }
+    }
+
+    return tresult;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>
+>
+Foam::faMeshSubset::interpolate
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const bool allowUnmapped
+) const
+{
+    if (subMeshPtr_)
+    {
+        return interpolate(vf, *subMeshPtr_);
+    }
+
+    return vf;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshTools.C b/src/finiteArea/faMesh/faMeshTools/faMeshTools.C
new file mode 100644
index 0000000000..24b46edf5e
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshTools.C
@@ -0,0 +1,470 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2012-2016 OpenFOAM Foundation
+    Copyright (C) 2015-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 "faMeshTools.H"
+#include "faBoundaryMeshEntries.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "polyMesh.H"
+#include "processorFaPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::faMeshTools::unregisterMesh(const faMesh& mesh)
+{
+    auto& obr = const_cast<objectRegistry&>(mesh.thisDb());
+
+    // Checkout by name (casting ambiguity)
+    obr.checkOut(faMesh::typeName);
+    obr.checkOut("faBoundaryMesh");
+    obr.checkOut("faSchemes");
+    obr.checkOut("faSolution");
+}
+
+
+void Foam::faMeshTools::forceDemandDriven(faMesh& mesh)
+{
+    (void)mesh.globalData();
+
+    (void)mesh.Le();
+    (void)mesh.magLe();
+    (void)mesh.areaCentres();
+    (void)mesh.edgeCentres();
+
+    (void)mesh.faceAreaNormals();
+    (void)mesh.edgeAreaNormals();
+    (void)mesh.pointAreaNormals();
+    (void)mesh.faceCurvatures();
+    (void)mesh.edgeTransformTensors();
+}
+
+
+Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
+(
+    const IOobject& io,
+    const polyMesh& pMesh,
+    const bool masterOnlyReading,
+    const bool verbose
+)
+{
+    // Region name
+    // ~~~~~~~~~~~
+
+    const fileName meshSubDir
+    (
+        (pMesh.name() == polyMesh::defaultRegion ? word::null : pMesh.name())
+      / faMesh::meshSubDir
+    );
+
+
+    fileName facesInstance;
+
+    // Patch types
+    // ~~~~~~~~~~~
+    // Read and scatter master patches (without reading master mesh!)
+
+    PtrList<entry> patchEntries;
+    if (Pstream::master())
+    {
+        const bool oldParRun = Pstream::parRun(false);
+
+        facesInstance = io.time().findInstance
+        (
+            meshSubDir,
+            "faceLabels",
+            IOobject::MUST_READ
+        );
+
+        patchEntries = faBoundaryMeshEntries
+        (
+            IOobject
+            (
+                "faBoundary",
+                facesInstance,
+                meshSubDir,
+                io.db(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        Pstream::parRun(oldParRun);
+    }
+
+    // Broadcast information to all
+    Pstream::broadcasts
+    (
+        UPstream::worldComm,
+        patchEntries,
+        facesInstance
+    );
+
+
+    // Dummy meshes
+    // ~~~~~~~~~~~~
+
+    // Set up to read-if-present. Note: does not search for mesh so set
+    // instance explicitly
+
+    IOobject meshIO(io);
+    meshIO.instance() = facesInstance;
+    meshIO.readOpt(IOobject::READ_IF_PRESENT);
+
+    // For mesh components (faceLabels, ...)
+    IOobject cmptIO(meshIO, "faceLabels", meshSubDir);
+    cmptIO.readOpt(IOobject::MUST_READ);
+    cmptIO.writeOpt(IOobject::NO_WRITE);
+    cmptIO.registerObject(false);
+
+
+    // Check who has a mesh
+
+    const fileName meshDir = io.time().path()/facesInstance/meshSubDir;
+    bool haveMesh = isDir(meshDir);
+    if (masterOnlyReading && !Pstream::master())
+    {
+        haveMesh = false;
+        meshIO.readOpt(IOobject::NO_READ);
+    }
+
+    if (!haveMesh)
+    {
+        cmptIO.readOpt(IOobject::NO_READ);
+    }
+
+
+    // Read mesh
+    // ~~~~~~~~~
+    // Now all processors use supplied points,faces etc
+    // Note: solution, schemes are also using the supplied IOobject so
+    //       on slave will be NO_READ, on master READ_IF_PRESENT. This will
+    //       conflict with e.g. timeStampMaster reading so switch off.
+
+    const auto oldCheckType = IOobject::fileModificationChecking;
+    IOobject::fileModificationChecking = IOobject::timeStamp;
+
+
+    // faceLabels
+    cmptIO.rename("faceLabels");
+    labelIOList faceLabels(cmptIO);
+
+
+    auto meshPtr = autoPtr<faMesh>::New
+    (
+        pMesh,
+        std::move(faceLabels),
+        meshIO
+    );
+    auto& mesh = *meshPtr;
+
+    IOobject::fileModificationChecking = oldCheckType;
+
+
+    // Some processors without patches? - add patches
+
+    if (returnReduce(mesh.boundary().empty(), orOp<bool>()))
+    {
+        // Use patchEntries, which were read on master and broadcast
+
+        faPatchList patches(patchEntries.size());
+        label nPatches = 0;
+
+        const bool isEmptyMesh = (mesh.faceLabels().empty());
+
+        forAll(patchEntries, patchi)
+        {
+            const entry& e = patchEntries[patchi];
+            const word type(e.dict().get<word>("type"));
+            const word& name = e.keyword();
+
+            if
+            (
+                type == processorFaPatch::typeName
+            )
+            {
+                // Stop at the first processor patch.
+                // - logic will not work with inter-mixed proc-patches anyhow
+                break;
+            }
+            else
+            {
+                dictionary patchDict(e.dict());
+
+                if (isEmptyMesh)
+                {
+                    patchDict.set("edgeLabels", labelList());
+                }
+
+                patches.set
+                (
+                    patchi,
+                    faPatch::New
+                    (
+                        name,
+                        patchDict,
+                        nPatches++,
+                        mesh.boundary()
+                    )
+                );
+            }
+        }
+
+        patches.resize(nPatches);
+        mesh.addFaPatches(patches, false);  // No parallel comms
+    }
+
+    // Recreate basic geometry, globalMeshData etc.
+    mesh.init(false);
+    (void)mesh.globalData();
+
+    return meshPtr;
+}
+
+
+Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::loadOrCreateMesh
+(
+    const IOobject& io,
+    const polyMesh& pMesh,
+    const bool decompose,
+    const bool verbose
+)
+{
+    // Region name
+    // ~~~~~~~~~~~
+
+    const fileName meshSubDir
+    (
+        (pMesh.name() == polyMesh::defaultRegion ? word::null : pMesh.name())
+      / faMesh::meshSubDir
+    );
+
+
+    // Patch types
+    // ~~~~~~~~~~~
+    // Read and scatter master patches (without reading master mesh!)
+
+    PtrList<entry> patchEntries;
+    if (Pstream::master())
+    {
+        const bool oldParRun = Pstream::parRun(false);
+
+        patchEntries = faBoundaryMeshEntries
+        (
+            IOobject
+            (
+                "faBoundary",
+                io.instance(),
+                meshSubDir,
+                io.db(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        Pstream::parRun(oldParRun);
+    }
+
+    // Broadcast: send patches to all
+    Pstream::broadcast(patchEntries);  // == 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 haveMesh =
+    (
+        decompose
+      ? Pstream::master()
+      : fileHandler().isFile
+        (
+            fileHandler().filePath
+            (
+                io.time().path()/io.instance()/meshSubDir/"faceLabels"
+            )
+        )
+    );
+
+
+    if (!haveMesh)
+    {
+        const bool oldParRun = Pstream::parRun(false);
+
+        // Create dummy mesh - on procs that don't already have a mesh
+        faMesh dummyMesh
+        (
+            pMesh,
+            labelList(),
+            IOobject(io, IOobject::NO_READ, IOobject::AUTO_WRITE)
+        );
+
+        // Add patches
+        faPatchList patches(patchEntries.size());
+        label nPatches = 0;
+
+        forAll(patchEntries, patchi)
+        {
+            const entry& e = patchEntries[patchi];
+            const word type(e.dict().get<word>("type"));
+            const word& name = e.keyword();
+
+            if
+            (
+                type == processorFaPatch::typeName
+            )
+            {
+                // Stop at the first processor patch.
+                // - logic will not work with inter-mixed proc-patches anyhow
+                break;
+            }
+            else
+            {
+                dictionary patchDict(e.dict());
+                patchDict.set("edgeLabels", labelList());
+
+                patches.set
+                (
+                    patchi,
+                    faPatch::New
+                    (
+                        name,
+                        patchDict,
+                        nPatches++,
+                        dummyMesh.boundary()
+                    )
+                );
+            }
+        }
+
+        patches.resize(nPatches);
+        dummyMesh.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());
+
+        //Pout<< "Writing dummy mesh to " << dummyMesh.boundary().path()
+        //    << endl;
+        dummyMesh.write();
+
+        Pstream::parRun(oldParRun);  // Restore parallel state
+    }
+
+    // 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;
+
+    // Sync patches
+    // ~~~~~~~~~~~~
+
+    if (!Pstream::master() && haveMesh)
+    {
+        // Check master names against mine
+
+        const faBoundaryMesh& patches = mesh.boundary();
+
+        forAll(patchEntries, patchi)
+        {
+            const entry& e = patchEntries[patchi];
+            const word type(e.dict().get<word>("type"));
+            const word& name = e.keyword();
+
+            if
+            (
+                type == processorFaPatch::typeName
+            )
+            {
+                break;
+            }
+
+            if (patchi >= patches.size())
+            {
+                FatalErrorInFunction
+                    << "Non-processor patches not synchronised."
+                    << endl
+                    << "Processor " << Pstream::myProcNo()
+                    << " has only " << patches.size()
+                    << " patches, master has "
+                    << patchi
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                type != patches[patchi].type()
+             || name != patches[patchi].name()
+            )
+            {
+                FatalErrorInFunction
+                    << "Non-processor patches not synchronised."
+                    << endl
+                    << "Master patch " << patchi
+                    << " name:" << type
+                    << " type:" << type << endl
+                    << "Processor " << Pstream::myProcNo()
+                    << " patch " << patchi
+                    << " has name:" << patches[patchi].name()
+                    << " type:" << patches[patchi].type()
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    // Recreate basic geometry, globalMeshData etc.
+    mesh.init(false);
+    (void)mesh.globalData();
+
+    /// #if 0
+    /// faMeshTools::forceDemandDriven(mesh);
+    /// faMeshTools::unregisterMesh(mesh);
+    /// #endif
+
+    // Do some checks.
+
+    // Check if the boundary definition is unique
+    // and processor patches are correct
+    mesh.boundary().checkDefinition(verbose);
+    mesh.boundary().checkParallelSync(verbose);
+
+    return meshPtr;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshTools.H b/src/finiteArea/faMesh/faMeshTools/faMeshTools.H
new file mode 100644
index 0000000000..82a6513a78
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshTools.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::faMeshTools
+
+Description
+    A collection of tools for operating on an faMesh.
+
+SourceFiles
+    faMeshTools.C
+    faMeshToolsProcAddr.C
+    faMeshToolsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_faMeshTools_H
+#define Foam_faMeshTools_H
+
+#include "faMesh.H"
+#include "areaFieldsFwd.H"
+#include "edgeFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class mapDistributePolyMesh;
+class writeHandler;
+class polyMesh;
+class IOobject;
+
+/*---------------------------------------------------------------------------*\
+                         Class faMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class faMeshTools
+{
+public:
+
+    //- Unregister the faMesh from its associated polyMesh
+    //- to prevent triggering on polyMesh changes etc.
+    static void unregisterMesh(const faMesh& mesh);
+
+    //- Force creation of everything that might vaguely be used by patches.
+    //  This is fairly horrible.
+    static void forceDemandDriven(faMesh& mesh);
+
+
+    //- Read mesh or create dummy mesh (0 faces, >0 patches).
+    //  Works in two modes according to masterOnlyReading:
+    //  true : create a dummy mesh for all procs
+    //  false: checks locally for mesh directories and only creates dummy mesh
+    //         if not present
+    static autoPtr<faMesh> newMesh
+    (
+        const IOobject& io,
+        const polyMesh& pMesh,
+        const bool masterOnlyReading,
+        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,
+        const bool decompose,
+        const bool verbose = false
+    );
+
+
+    //- Read decompose/reconstruct addressing
+    static mapDistributePolyMesh readProcAddressing
+    (
+        const faMesh& mesh,
+        const autoPtr<faMesh>& baseMeshPtr
+    );
+
+    //- Write decompose/reconstruct addressing
+    //
+    //  \param mesh
+    //  \param faDistMap
+    //  \param decompose  running in decompose vs reconstruct mode
+    //  \param writeHandler file handler
+    //  \param procMesh   (optional) processor mesh in reconstruct mode
+    //
+    //  \note Since the faMesh holds a reference to a polyMesh,
+    //  in reconstruct mode it will refer to the base mesh, but
+    //  we need a means to proc addressing into the processor locations.
+    //  This is the purpose of the additional procMesh reference
+    static void writeProcAddressing
+    (
+        const faMesh& mesh,
+        const mapDistributePolyMesh& faDistMap,
+        const bool decompose,
+        autoPtr<fileOperation>&& writeHandler,
+        const faMesh* procMesh = nullptr
+    );
+
+
+    //- Flatten an edge field into linear addressing
+    //  Optionally use primitive patch edge ordering
+    template<class Type>
+    static tmp<Field<Type>> flattenEdgeField
+    (
+        const GeometricField<Type, faePatchField, edgeMesh>& fld,
+        const bool primitiveOrdering = false
+    );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "faMeshToolsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshToolsProcAddr.C b/src/finiteArea/faMesh/faMeshTools/faMeshToolsProcAddr.C
new file mode 100644
index 0000000000..b8a42738d0
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshToolsProcAddr.C
@@ -0,0 +1,468 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faMeshTools.H"
+#include "BitOps.H"
+#include "fileOperation.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "IOmapDistributePolyMesh.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Create a reconstruct map.
+
+static mapDistributePolyMesh createReconstructMap
+(
+    const faMesh& mesh,
+    const autoPtr<faMesh>& baseMeshPtr,
+    const labelUList& faceProcAddr,
+    const labelUList& edgeProcAddr,
+    const labelUList& pointProcAddr,
+    const labelUList& boundaryProcAddr
+)
+{
+    const label nOldPoints = mesh.nPoints();
+    const label nOldFaces = mesh.nFaces();
+    const label nOldEdges = mesh.nEdges();
+
+    ///Pout<< "old sizes"
+    ///    << " points:" << nOldPoints
+    ///    << " faces:" << nOldFaces
+    ///    << " edges:" << nOldEdges << nl;
+
+    const faBoundaryMesh& oldBndMesh = mesh.boundary();
+    labelList oldPatchStarts(oldBndMesh.patchStarts());
+
+    // Patches: purge -1 entries
+    labelList patchProcAddr
+    (
+        IndirectList<label>::subset_if
+        (
+            boundaryProcAddr,
+            labelRange::ge0()
+        )
+    );
+
+
+    labelListList faceSubMap(Pstream::nProcs());
+    faceSubMap[Pstream::masterNo()] = identity(nOldFaces);
+
+    labelListList edgeSubMap(Pstream::nProcs());
+    edgeSubMap[Pstream::masterNo()] = identity(nOldEdges);
+
+    labelListList pointSubMap(Pstream::nProcs());
+    pointSubMap[Pstream::masterNo()] = identity(nOldPoints);
+
+    labelListList patchSubMap(Pstream::nProcs());
+    patchSubMap[Pstream::masterNo()] = patchProcAddr;
+
+
+    // Gather addressing on the master
+    labelListList faceAddressing(Pstream::nProcs());
+    faceAddressing[Pstream::myProcNo()] = faceProcAddr;
+    Pstream::gatherList(faceAddressing);
+
+    labelListList edgeAddressing(Pstream::nProcs());
+    edgeAddressing[Pstream::myProcNo()] = edgeProcAddr;
+    Pstream::gatherList(edgeAddressing);
+
+    labelListList pointAddressing(Pstream::nProcs());
+    pointAddressing[Pstream::myProcNo()] = pointProcAddr;
+    Pstream::gatherList(pointAddressing);
+
+    labelListList patchAddressing(Pstream::nProcs());
+    patchAddressing[Pstream::myProcNo()] = patchProcAddr;
+    Pstream::gatherList(patchAddressing);
+
+
+    // NB: can only have a reconstruct on master!
+    if (Pstream::master() && baseMeshPtr && baseMeshPtr->nFaces())
+    {
+        const faMesh& baseMesh = *baseMeshPtr;
+
+        const label nNewPoints = baseMesh.nPoints();
+        const label nNewFaces = baseMesh.nFaces();
+        const label nNewEdges = baseMesh.nEdges();
+        const label nNewPatches = baseMesh.boundary().size();
+
+        /// Pout<< "new sizes"
+        ///     << " points:" << nNewPoints
+        ///     << " faces:" << nNewFaces
+        ///     << " edges:" << nNewEdges << nl;
+
+        mapDistribute faFaceMap
+        (
+            nNewFaces,
+            std::move(faceSubMap),
+            std::move(faceAddressing),
+            false,  // subHasFlip
+            false   // constructHasFlip
+        );
+
+        mapDistribute faEdgeMap
+        (
+            nNewEdges,
+            std::move(edgeSubMap),
+            std::move(edgeAddressing),
+            false,  // subHasFlip
+            false   // constructHasFlip
+        );
+
+        mapDistribute faPointMap
+        (
+            nNewPoints,
+            std::move(pointSubMap),
+            std::move(pointAddressing)
+        );
+
+        mapDistribute faPatchMap
+        (
+            nNewPatches,
+            std::move(patchSubMap),
+            std::move(patchAddressing)
+        );
+
+        return mapDistributePolyMesh
+        (
+            // Mesh before changes
+            nOldPoints,
+            nOldEdges,          // area: nOldEdges (volume: nOldFaces)
+            nOldFaces,          // area: nOldFaces (volume: nOldCells)
+
+            std::move(oldPatchStarts),
+            labelList(),        // oldPatchNMeshPoints [unused]
+
+            mapDistribute(std::move(faPointMap)),
+            mapDistribute(std::move(faEdgeMap)), // edgeMap (volume: faceMap)
+            mapDistribute(std::move(faFaceMap)), // faceMap (volume: cellMap)
+            mapDistribute(std::move(faPatchMap))
+        );
+    }
+    else
+    {
+        mapDistribute faFaceMap
+        (
+            0,  // nNewFaces
+            std::move(faceSubMap),
+            labelListList(Pstream::nProcs()),   // constructMap
+            false,  // subHasFlip
+            false   // constructHasFlip
+        );
+
+        mapDistribute faEdgeMap
+        (
+            0,  // nNewEdges
+            std::move(edgeSubMap),
+            labelListList(Pstream::nProcs()),   // constructMap
+            false,  // subHasFlip
+            false   // constructHasFlip
+        );
+
+        mapDistribute faPointMap
+        (
+            0,  // nNewPoints
+            std::move(pointSubMap),
+            labelListList(Pstream::nProcs())    // constructMap
+        );
+
+        mapDistribute faPatchMap
+        (
+            0,  // nNewPatches
+            std::move(patchSubMap),
+            labelListList(Pstream::nProcs())    // constructMap
+        );
+
+        return mapDistributePolyMesh
+        (
+            // Mesh before changes
+            nOldPoints,
+            nOldEdges,          // area: nOldEdges (volume: nOldFaces)
+            nOldFaces,          // area: nOldFaces (volume: nOldCells)
+
+            std::move(oldPatchStarts),
+            labelList(),        // oldPatchNMeshPoints [unused]
+
+            mapDistribute(std::move(faPointMap)),
+            mapDistribute(std::move(faEdgeMap)), // edgeMap (volume: faceMap)
+            mapDistribute(std::move(faFaceMap)), // faceMap (volume: cellMap)
+            mapDistribute(std::move(faPatchMap))
+        );
+    }
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::mapDistributePolyMesh
+Foam::faMeshTools::readProcAddressing
+(
+    const faMesh& mesh,
+    const autoPtr<faMesh>& baseMeshPtr
+)
+{
+    IOobject ioAddr
+    (
+        "procAddressing",
+        mesh.facesInstance(),
+        faMesh::meshSubDir,
+        mesh.thisDb(),
+        IOobject::READ_IF_PRESENT,
+        IOobject::NO_WRITE,
+        false  // no register
+    );
+
+    //if (ioAddr.typeHeaderOk<labelIOList>(true))
+    //{
+    //    Pout<< "Reading addressing from " << io.name() << " at "
+    //        << mesh.facesInstance() << nl << endl;
+    //    distMap.reset(new IOmapDistributePolyMesh(io));
+    //}
+    //else
+
+    {
+        Info<< "Reading (face|edge|face|point|boundary)ProcAddressing from "
+            << mesh.facesInstance().c_str() << '/'
+            << faMesh::meshSubDir << nl << endl;
+
+        ioAddr.rename("faceProcAddressing");
+        labelIOList faceProcAddressing(ioAddr, Zero);
+
+        ioAddr.rename("edgeProcAddressing");
+        labelIOList edgeProcAddressing(ioAddr, Zero);
+
+        ioAddr.rename("pointProcAddressing");
+        labelIOList pointProcAddressing(ioAddr, Zero);
+
+        ioAddr.rename("boundaryProcAddressing");
+        labelIOList boundaryProcAddressing(ioAddr, Zero);
+
+        if
+        (
+            mesh.nFaces() != faceProcAddressing.size()
+         || mesh.nEdges() != edgeProcAddressing.size()
+         || mesh.nPoints() != pointProcAddressing.size()
+         || mesh.boundary().size() != boundaryProcAddressing.size()
+        )
+        {
+            FatalErrorInFunction
+                << "Read addressing inconsistent with mesh sizes" << nl
+                << "faces:" << mesh.nFaces()
+                << " addressing:" << faceProcAddressing.objectRelPath()
+                << " size:" << faceProcAddressing.size() << nl
+                << "edges:" << mesh.nEdges()
+                << " addressing:" << edgeProcAddressing.objectRelPath()
+                << " size:" << edgeProcAddressing.size() << nl
+                << "points:" << mesh.nPoints()
+                << " addressing:" << pointProcAddressing.objectRelPath()
+                << " size:" << pointProcAddressing.size()
+                << "patches:" << mesh.boundary().size()
+                << " addressing:" << boundaryProcAddressing.objectRelPath()
+                << " size:" << boundaryProcAddressing.size()
+                << exit(FatalError);
+        }
+
+        return createReconstructMap
+        (
+            mesh,
+            baseMeshPtr,
+            faceProcAddressing,
+            edgeProcAddressing,
+            pointProcAddressing,
+            boundaryProcAddressing
+        );
+    }
+}
+
+
+void Foam::faMeshTools::writeProcAddressing
+(
+    const faMesh& mesh,
+    const mapDistributePolyMesh& map,
+    const bool decompose,
+    autoPtr<fileOperation>&& writeHandler,
+    const faMesh* procMesh
+)
+{
+    Info<< "Writing ("
+        << (decompose ? "decompose" : "reconstruct")
+        << ") procAddressing files to "
+        << mesh.facesInstance().c_str() << '/'
+        << faMesh::meshSubDir << endl;
+
+    IOobject ioAddr
+    (
+        "procAddressing",
+        mesh.facesInstance(),
+        faMesh::meshSubDir,
+        (procMesh && !decompose ? procMesh->thisDb() : mesh.thisDb()),
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false  // no register
+    );
+
+
+    // faceProcAddressing (faMesh)
+    ioAddr.rename("faceProcAddressing");
+    labelIOList faceMap(ioAddr, Zero);
+
+    // edgeProcAddressing (faMesh)
+    ioAddr.rename("edgeProcAddressing");
+    labelIOList edgeMap(ioAddr, Zero);
+
+    // pointProcAddressing (faMesh)
+    ioAddr.rename("pointProcAddressing");
+    labelIOList pointMap(ioAddr, Zero);
+
+    // boundaryProcAddressing (faMesh)
+    ioAddr.rename("boundaryProcAddressing");
+    labelIOList patchMap(ioAddr, Zero);
+
+    if (decompose)
+    {
+        // Decompose
+        // - forward map:  [undecomposed] -> [decomposed]
+
+        // area:faces (volume:cells)
+        faceMap = identity(map.nOldCells());
+        map.cellMap().distribute(faceMap);
+
+        // area:edges (volume:faces)
+        edgeMap = identity(map.nOldFaces());
+        map.faceMap().distribute(edgeMap);
+
+        pointMap = identity(map.nOldPoints());
+        map.distributePointData(pointMap);
+
+        patchMap = identity(map.patchMap().constructSize());
+        map.patchMap().mapDistributeBase::distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            label(-1),  // nullValue for new patches...
+            patchMap,
+            flipOp()    // negate op
+        );
+    }
+    else    // reconstruct
+    {
+        // Reconstruct
+        // - reverse map:  [undecomposed] <- [decomposed]
+
+        // area:faces (volume:cells)
+        faceMap = identity(mesh.nFaces());
+        map.cellMap().reverseDistribute(map.nOldCells(), faceMap);
+
+        // area:edges (volume:faces)
+        edgeMap = identity(mesh.patch().nEdges());
+        map.faceMap().reverseDistribute(map.nOldFaces(), edgeMap);
+
+        pointMap = identity(mesh.nPoints());
+        map.pointMap().reverseDistribute(map.nOldPoints(), pointMap);
+
+        patchMap = identity(mesh.boundary().size());
+        map.patchMap().mapDistributeBase::reverseDistribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            map.oldPatchSizes().size(),
+            label(-1),  // nullValue for unmapped patches...
+            patchMap
+        );
+    }
+
+    autoPtr<fileOperation> defaultHandler;
+    if (writeHandler)
+    {
+        defaultHandler = fileHandler(std::move(writeHandler));
+    }
+
+
+    // If we want procAddressing, need to manually write it ourselves
+    // since it was not registered anywhere
+
+    IOmapDistributePolyMeshRef procAddrMap
+    (
+        IOobject
+        (
+            "procAddressing",
+            mesh.facesInstance(),
+            faMesh::meshSubDir,
+            mesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false  // no register
+        ),
+        map
+    );
+
+
+    if (decompose)
+    {
+        // Write into proc directories
+        procAddrMap.write();
+    }
+    else
+    {
+        // Reconstruct: "procAddressing" only meaningful for rank 0
+        // and written into base (serial) location (if at all).
+
+        if (Pstream::master())
+        {
+            const bool oldParRun = Pstream::parRun(false);
+            procAddrMap.write();
+            Pstream::parRun(oldParRun);
+        }
+    }
+
+
+    const bool faceOk = faceMap.write();
+    const bool edgeOk = edgeMap.write();
+    const bool pointOk = pointMap.write();
+    const bool patchOk = patchMap.write();
+
+    if (defaultHandler)
+    {
+        writeHandler = fileHandler(std::move(defaultHandler));
+    }
+
+    if (!edgeOk || !faceOk || !pointOk || !patchOk)
+    {
+        WarningInFunction
+            << "Failed to write some of "
+            << faceMap.objectRelPath() << ", "
+            << edgeMap.objectRelPath() << ", "
+            << pointMap.objectRelPath() << ", "
+            << patchMap.objectRelPath() << endl;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshTools/faMeshToolsTemplates.C b/src/finiteArea/faMesh/faMeshTools/faMeshToolsTemplates.C
new file mode 100644
index 0000000000..74d482a1d9
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshTools/faMeshToolsTemplates.C
@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "edgeFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::faMeshTools::flattenEdgeField
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& fld,
+    const bool primitiveOrdering
+)
+{
+    const faMesh& mesh = fld.mesh();
+
+    auto tresult = tmp<Field<Type>>::New(mesh.nEdges(), Zero);
+    auto& result = tresult.ref();
+
+    // Internal field
+    result.slice(0, fld.size()) = fld;
+
+    if (primitiveOrdering)
+    {
+        // Boundary field in primitive patch order
+
+        forAll(fld.boundaryField(), patchi)
+        {
+            UIndirectList<Type>
+            (
+                result,
+                mesh.boundary()[patchi].edgeLabels()
+            ) = fld.boundaryField()[patchi];
+        }
+    }
+    else
+    {
+        // Boundary field in sub-list (slice) order
+
+        label start = fld.size();
+
+        forAll(fld.boundaryField(), patchi)
+        {
+            const label len = mesh.boundary()[patchi].size();
+
+            result.slice(start, len) = fld.boundaryField()[patchi];
+
+            start += len;
+        }
+    }
+
+    return tresult;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshUpdate.C b/src/finiteArea/faMesh/faMeshUpdate.C
index 1e8593b261..b1136484bf 100644
--- a/src/finiteArea/faMesh/faMeshUpdate.C
+++ b/src/finiteArea/faMesh/faMeshUpdate.C
@@ -37,8 +37,6 @@ License
 
 void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
 {
-    DebugInFunction << "Updating mesh" << endl;
-
     // if (!mpm.morphing())
     // {
     //     // No topo change
@@ -48,18 +46,14 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     // Create fa mesh mapper, using the old mesh
     const faMeshMapper mapper(*this, mpm);
 
-
     // Rebuild mesh
-
-    // Cast away const for interface reasons.  HJ, 12/Aug/2011
-    faMesh& m = const_cast<faMesh&>(*this);
-
+    // ~~~~~~~~~~~~
 
     // Clear existing mesh data
     clearOut();
 
     // Set new labels
-    m.faceLabels_ = mapper.areaMap().newFaceLabels();
+    faceLabels_ = mapper.areaMap().newFaceLabels();
 
     const uindirectPrimitivePatch& bp = patch();
 
@@ -121,12 +115,12 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     }
 
     // Set new edges for all patches
-    forAll(m.boundary_, patchI)
+    forAll(boundary_, patchI)
     {
-        m.boundary_[patchI].resetEdges(patchEdges[patchI]);
+        boundary_[patchI].resetEdges(patchEdges[patchI]);
     }
 
-    m.setPrimitiveMeshData();
+    setPrimitiveMeshData();
 
     // Create global mesh data
     if (Pstream::parRun())
@@ -135,10 +129,10 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     }
 
     // Calculate topology for the patches (processor-processor comms etc.)
-    m.boundary_.updateMesh();
+    boundary_.updateMesh();
 
     // Calculate the geometry for the patches (transformation tensors etc.)
-    m.boundary_.calcGeometry();
+    boundary_.calcGeometry();
 
 
     // Map fields
@@ -156,7 +150,7 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
 
 void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
 {
-    // Map all the areaFields in the objectRegistry
+    // Map areaFields in the objectRegistry
     MapGeometricFields<scalar, faPatchField, faMeshMapper, areaMesh>(mapper);
     MapGeometricFields<vector, faPatchField, faMeshMapper, areaMesh>(mapper);
     MapGeometricFields<sphericalTensor, faPatchField, faMeshMapper, areaMesh>
@@ -165,7 +159,7 @@ void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
         (mapper);
     MapGeometricFields<tensor, faPatchField, faMeshMapper, areaMesh>(mapper);
 
-    // Map all the edgeFields in the objectRegistry
+    // Map edgeFields in the objectRegistry
     MapGeometricFields<scalar, faePatchField, faMeshMapper, edgeMesh>(mapper);
     MapGeometricFields<vector, faePatchField, faMeshMapper, edgeMesh>(mapper);
     MapGeometricFields<sphericalTensor, faePatchField, faMeshMapper, edgeMesh>
@@ -227,7 +221,6 @@ void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
             }
         }
     }
-
 }