diff --git a/applications/utilities/finiteArea/checkFaMesh/Make/options b/applications/utilities/finiteArea/checkFaMesh/Make/options
index b69f32d758..2a0d2e96ce 100644
--- a/applications/utilities/finiteArea/checkFaMesh/Make/options
+++ b/applications/utilities/finiteArea/checkFaMesh/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
-    -lfiniteVolume \
+    -lfileFormats \
+    -lsurfMesh \
     -lmeshTools
diff --git a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
index d83aaeb7f2..e15b830083 100644
--- a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
+++ b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,14 +30,22 @@ Application
 Description
     Check a finiteArea mesh
 
-Author
+Original Authors
     Zeljko Tukovic, FAMENA
     Hrvoje Jasak, Wikki Ltd.
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
-#include "faCFD.H"
+#include "Time.H"
+#include "argList.H"
+#include "faMesh.H"
+#include "polyMesh.H"
+#include "areaFaMesh.H"
+#include "edgeFaMesh.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "processorFaPatch.H"
+#include "foamVtkUIndPatchWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,37 +60,32 @@ int main(int argc, char *argv[])
         "Check a finiteArea mesh"
     );
 
-    #include "addRegionOption.H"
+    argList::addBoolOption
+    (
+        "write-vtk",
+        "Write mesh as a vtp (vtk) file for display or debugging"
+    );
 
+    #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createNamedMesh.H"
-    #include "createFaMesh.H"
+    #include "createNamedPolyMesh.H"
+
+
+    // Create
+    faMesh aMesh(mesh);
 
     Info<< "Time = " << runTime.timeName() << nl << endl;
 
-    // General mesh statistics
-    Info<< "Number of points: " << aMesh.nPoints() << nl
-        << "Number of internal edges: " << aMesh.nInternalEdges() << nl
-        << "Number of edges: " << aMesh.nEdges() << nl
-        << "Number of faces: " << aMesh.nFaces() << nl
-        << endl;
+    #include "printMeshSummary.H"
 
-    // Check geometry
-    Info<< "Face area:  min = " << min(aMesh.S().field())
-        << " max = "  << max(aMesh.S().field()) << nl
-        << "Internal edge length: min = "
-        << min(aMesh.magLe().internalField()) << nl
-        << " max = "  << max(aMesh.magLe().internalField()) << nl
-        << "Edge length: min = "
-        << min(aMesh.magLe()).value() << nl
-        << " max = "  << max(aMesh.magLe()).value() << nl
-        << "Face area normals:  min = " << min(aMesh.faceAreaNormals().field())
-        << " max = "  << max(aMesh.faceAreaNormals().field()) << nl
-        << endl;
+    if (args.found("write-vtk"))
+    {
+        #include "faMeshWriteVTK.H"
+    }
 
+    Info<< "\nEnd\n" << endl;
 
-    Info << "\nEnd" << endl;
     return 0;
 }
 
diff --git a/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H b/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
new file mode 100644
index 0000000000..1df8b4d312
--- /dev/null
+++ b/applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    VTK output of faMesh with some geometric or debug fields
+
+\*---------------------------------------------------------------------------*/
+
+{
+    // finiteArea
+    vtk::uindirectPatchWriter writer
+    (
+        aMesh.patch(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // CellData
+    writer.beginCellData(3);
+    writer.writeProcIDs();
+    writer.write("area", aMesh.S().field());
+    writer.write("normal", aMesh.faceAreaNormals());
+
+    // PointData
+    writer.beginPointData(1);
+    writer.write("normal", aMesh.pointAreaNormals());
+
+    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
new file mode 100644
index 0000000000..256d2b1371
--- /dev/null
+++ b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Summary of faMesh information
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const faBoundaryMesh& patches = aMesh.boundary();
+    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;
+
+    if (Pstream::parRun())
+    {
+        for (const faPatch& fap : patches)
+        {
+            const auto* cpp = isA<processorFaPatch>(fap);
+
+            if (cpp)
+            {
+                nProcEdges += fap.nEdges();
+            }
+        }
+    }
+
+    const label nBoundEdges = aMesh.nBoundaryEdges() - nProcEdges;
+
+    Info<< "  Number of boundary edges: "
+        << returnReduce(nBoundEdges - nProcEdges, sumOp<label>()) << nl;
+
+    if (Pstream::parRun())
+    {
+        Info<< "  Number of processor edges: "
+            << returnReduce(nProcEdges, sumOp<label>()) << nl;
+    }
+
+
+    Info<< "----------------" << nl
+        << "Patches" << nl
+        << "----------------" << nl;
+
+    for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+    {
+        const faPatch& p = patches[patchi];
+
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
+    }
+
+
+    // Geometry information
+    Info<< nl;
+    {
+        scalarMinMax limit(gMinMax(aMesh.S().field()));
+        Info<< "Face area:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    {
+        scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
+
+        // Include processor boundaries into 'internal' edges
+        if (Pstream::parRun())
+        {
+            for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
+            {
+                limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+            }
+
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length (internal):" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+
+
+        // Include (non-processor) boundaries
+        for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+        {
+            limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    // Not particularly meaningful
+    #if 0
+    {
+        MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
+
+        Info<< "Face area normals:" << nl
+            << "    min = " << limit.min() << " max = " << limit.max() << nl;
+    }
+    #endif
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/finiteArea/makeFaMesh/Make/options b/applications/utilities/finiteArea/makeFaMesh/Make/options
index 7a28b2c4b6..3f949d5344 100644
--- a/applications/utilities/finiteArea/makeFaMesh/Make/options
+++ b/applications/utilities/finiteArea/makeFaMesh/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/faDecompose/lnInclude \
     -I$(LIB_SRC)/parallel/reconstruct/faReconstruct/lnInclude
@@ -8,6 +9,7 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteArea \
     -lfileFormats \
+    -lsurfMesh \
     -lmeshTools \
     -lfaDecompose \
     -lfaReconstruct
diff --git a/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H b/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
index 5c38c9b4f1..3b91b260ca 100644
--- a/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
+++ b/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
@@ -17,9 +17,11 @@ Description
 
 if (Pstream::parRun())
 {
-    faMeshReconstructor reconstructor(areaMesh);
+    faMeshReconstructor reconstructor(aMesh);
     reconstructor.writeAddressing();
 
+    Info<< "Wrote proc-addressing" << nl << endl;
+
     // Handle area fields
     // ------------------
 
@@ -72,22 +74,53 @@ if (Pstream::parRun())
 
     if (nAreaFields)
     {
-        Info<< "Decomposing " << nAreaFields << " area fields" << nl << endl;
+        Info<< "Decomposing " << nAreaFields << " area fields" << nl;
 
         faFieldDecomposer fieldDecomposer
         (
             fullMesh,
-            areaMesh,
+            aMesh,
             reconstructor.edgeProcAddressing(),
             reconstructor.faceProcAddressing(),
             reconstructor.boundaryProcAddressing()
         );
 
-        fieldDecomposer.decomposeFields(areaScalarFields);
-        fieldDecomposer.decomposeFields(areaVectorFields);
-        fieldDecomposer.decomposeFields(areaSphTensorFields);
-        fieldDecomposer.decomposeFields(areaSymmTensorFields);
-        fieldDecomposer.decomposeFields(areaTensorFields);
+        if (areaScalarFields.size())
+        {
+            Info<< "  scalars: "
+                << flatOutput(PtrListOps::names(areaScalarFields)) << nl;
+            fieldDecomposer.decomposeFields(areaScalarFields);
+        }
+
+        if (areaVectorFields.size())
+        {
+            Info<< "  vectors: "
+                << flatOutput(PtrListOps::names(areaVectorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaVectorFields);
+        }
+
+        if (areaSphTensorFields.size())
+        {
+            Info<< "  sphTensors: "
+                << flatOutput(PtrListOps::names(areaSphTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaSphTensorFields);
+        }
+
+        if (areaSymmTensorFields.size())
+        {
+            Info<< "  symmTensors: "
+                << flatOutput(PtrListOps::names(areaSymmTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaSymmTensorFields);
+        }
+
+        if (areaTensorFields.size())
+        {
+            Info<< "  tensors: "
+                << flatOutput(PtrListOps::names(areaTensorFields)) << nl;
+            fieldDecomposer.decomposeFields(areaTensorFields);
+        }
+
+        Info<< endl;
     }
 }
 
diff --git a/applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H b/applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H
index a50d18df3c..a5f4192d39 100644
--- a/applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H
+++ b/applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H
@@ -16,15 +16,18 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    Info<< "Writing edges in obj format" << endl;
+    Info<< nl
+        << "Writing edges in obj format" << endl;
 
-    word outputName("faMesh-edges.obj");
+    word outputName("finiteArea-edges.obj");
 
     if (Pstream::parRun())
     {
         outputName = word
         (
-            "faMesh-edges-" + Foam::name(Pstream::myProcNo()) + ".obj"
+            "finiteArea-edges-proc"
+          + Foam::name(Pstream::myProcNo())
+          + ".obj"
         );
     }
 
@@ -36,7 +39,7 @@ Description
         false
     );
 
-    os.write(areaMesh.patch().edges(), areaMesh.patch().localPoints());
+    os.write(aMesh.patch().edges(), aMesh.patch().localPoints());
 }
 
 
diff --git a/applications/utilities/finiteArea/makeFaMesh/faMeshWriteVTK.H b/applications/utilities/finiteArea/makeFaMesh/faMeshWriteVTK.H
new file mode 100644
index 0000000000..1df8b4d312
--- /dev/null
+++ b/applications/utilities/finiteArea/makeFaMesh/faMeshWriteVTK.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    VTK output of faMesh with some geometric or debug fields
+
+\*---------------------------------------------------------------------------*/
+
+{
+    // finiteArea
+    vtk::uindirectPatchWriter writer
+    (
+        aMesh.patch(),
+        // vtk::formatType::INLINE_ASCII,
+        fileName
+        (
+            aMesh.mesh().time().globalPath() / "finiteArea"
+        )
+    );
+
+    writer.writeGeometry();
+
+    // CellData
+    writer.beginCellData(3);
+    writer.writeProcIDs();
+    writer.write("area", aMesh.S().field());
+    writer.write("normal", aMesh.faceAreaNormals());
+
+    // PointData
+    writer.beginPointData(1);
+    writer.write("normal", aMesh.pointAreaNormals());
+
+    Info<< nl
+        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
index e80ed303d9..f93685ea48 100644
--- a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
+++ b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
@@ -32,7 +32,7 @@ Description
     When called in parallel, it will also try to act like decomposePar,
     create procAddressing and decompose serial finite-area fields.
 
-Author
+Original Authors
     Zeljko Tukovic, FAMENA
     Hrvoje Jasak, Wikki Ltd.
 
@@ -48,6 +48,8 @@ Author
 #include "areaFields.H"
 #include "faFieldDecomposer.H"
 #include "faMeshReconstructor.H"
+#include "PtrListOps.H"
+#include "foamVtkUIndPatchWriter.H"
 #include "OBJstream.H"
 
 using namespace Foam;
@@ -69,11 +71,21 @@ int main(int argc, char *argv[])
     );
     argList::addOption("dict", "file", "Alternative faMeshDefinition");
 
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Create but do not write"
+    );
+    argList::addBoolOption
+    (
+        "write-vtk",
+        "Write mesh as a vtp (vtk) file for display or debugging"
+    );
     argList::addBoolOption
     (
         "write-edges-obj",
-        "Write mesh edges as obj files and exit",
-        false  // could make an advanced option
+        "Write mesh edges as obj files (one per processor)",
+        true  // advanced option (debugging only)
     );
 
     #include "addRegionOption.H"
@@ -81,6 +93,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createNamedPolyMesh.H"
 
+    const bool dryrun = args.found("dry-run");
+
     // Reading faMeshDefinition dictionary
     #include "findMeshDefinitionDict.H"
 
@@ -91,33 +105,40 @@ int main(int argc, char *argv[])
         meshDefDict.add("emptyPatch", patchName, true);
     }
 
-    // Create
-    faMesh areaMesh(mesh, meshDefDict);
 
-    bool quickExit = false;
+    // Create
+    faMesh aMesh(mesh, meshDefDict);
+
+    // Mesh information
+    #include "printMeshSummary.H"
 
     if (args.found("write-edges-obj"))
     {
-        quickExit = true;
         #include "faMeshWriteEdgesOBJ.H"
     }
 
-    if (quickExit)
+    if (args.found("write-vtk"))
     {
-        Info<< "\nEnd\n" << endl;
-        return 0;
+        #include "faMeshWriteVTK.H"
     }
 
-    // Set the precision of the points data to 10
-    IOstream::defaultPrecision(10);
+    if (dryrun)
+    {
+        Info<< "\ndry-run: not writing mesh or decomposing fields\n" << nl;
+    }
+    else
+    {
+        // Set the precision of the points data to 10
+        IOstream::defaultPrecision(10);
 
-    Info<< nl << "Write finite area mesh." << nl;
-    areaMesh.write();
-    Info<< endl;
+        Info<< nl << "Write finite area mesh." << nl;
+        aMesh.write();
 
-    #include "decomposeFaFields.H"
+        Info<< endl;
+        #include "decomposeFaFields.H"
+    }
 
-    Info << "\nEnd\n" << endl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H b/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
new file mode 100644
index 0000000000..5dd75c8804
--- /dev/null
+++ b/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Summary of faMesh information
+
+\*---------------------------------------------------------------------------*/
+
+{
+    const faBoundaryMesh& patches = aMesh.boundary();
+    const label nNonProcessor = patches.nNonProcessor();
+    const label nPatches = patches.size();
+
+    Info<< "----------------" << nl
+        << "Mesh Information" << nl
+        << "----------------" << nl
+        << "  " << "boundingBox: " << boundBox(aMesh.points()) << nl
+        << "  " << "nFaces: " << returnReduce(aMesh.nFaces(), sumOp<label>())
+        << nl;
+
+
+    Info<< "----------------" << nl
+        << "Patches" << nl
+        << "----------------" << nl;
+
+    for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+    {
+        const faPatch& p = patches[patchi];
+
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
+    }
+
+    // Geometry information
+    Info<< nl;
+    {
+        scalarMinMax limit(gMinMax(aMesh.S().field()));
+        Info<< "Face area:" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+    }
+
+    {
+        scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
+
+        // Include processor boundaries into 'internal' edges
+        if (Pstream::parRun())
+        {
+            for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
+            {
+                limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+            }
+
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length (internal):" << nl
+            << "    min = " << limit.min() << " max = "  << limit.max() << nl;
+
+
+        // Include (non-processor) boundaries
+        for (label patchi = 0; patchi < nNonProcessor; ++patchi)
+        {
+            limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
+        }
+
+        if (Pstream::parRun())
+        {
+            reduce(limit, minMaxOp<scalar>());
+        }
+
+        Info<< "Edge length:" << nl
+            << "    min = " << limit.min()
+            << " max = "  << limit.max() << nl;
+    }
+
+    // Not particularly meaningful
+    #if 0
+    {
+        MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
+
+        Info<< "Face area normals:" << nl
+            << "    min = " << limit.min() << " max = " << limit.max() << nl;
+    }
+    #endif
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
index c40ff59d4e..cd4eee3537 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@@ -75,8 +75,15 @@ Foam::faBoundaryMesh::faBoundaryMesh
     regIOobject(io),
     mesh_(mesh)
 {
-    if (readOpt() == IOobject::MUST_READ)
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+    )
     {
+        // Warn for MUST_READ_IF_MODIFIED
+        warnNoRereading<faBoundaryMesh>();
+
         faPatchList& patches = *this;
 
         // Read faPatch list
diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C
index fea1ea36e8..a98c49060d 100644
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@@ -343,7 +343,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(pMesh.time().timePath()/mesh().dbDir()/"S0"))
+    if (fileHandler().isFile(pMesh.time().timePath()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (
@@ -508,7 +508,7 @@ Foam::faMesh::faMesh
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(mesh().time().timePath()/"S0"))
+    if (fileHandler().isFile(pMesh.time().timePath()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (
diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
index 0bee133c30..32f3e4ca5c 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,34 +47,33 @@ const Foam::scalar Foam::cyclicFaPatch::matchTol_ = 1e-3;
 
 void Foam::cyclicFaPatch::calcTransforms()
 {
+    const label sizeby2 = this->size()/2;
     if (size() > 0)
     {
-        // const label sizeby2 = this->size()/2;
-        pointField half0Ctrs(size()/2);
-        pointField half1Ctrs(size()/2);
-        for (label i=0; i<size()/2; ++i)
+        pointField half0Ctrs(sizeby2);
+        pointField half1Ctrs(sizeby2);
+        for (label edgei=0; edgei < sizeby2; ++edgei)
         {
-            half0Ctrs[i] = this->edgeCentres()[i];
-            half1Ctrs[i] = this->edgeCentres()[i+size()/2];
+            half0Ctrs[edgei] = this->edgeCentres()[edgei];
+            half1Ctrs[edgei] = this->edgeCentres()[edgei + sizeby2];
         }
 
-        vectorField half0Normals(size()/2);
-        vectorField half1Normals(size()/2);
+        vectorField half0Normals(sizeby2);
+        vectorField half1Normals(sizeby2);
 
         const vectorField eN(edgeNormals()*magEdgeLengths());
 
         scalar maxMatchError = 0;
         label errorEdge = -1;
 
-        for (label edgei = 0; edgei < size()/2; ++edgei)
+        for (label edgei = 0; edgei < sizeby2; ++edgei)
         {
             half0Normals[edgei] = eN[edgei];
-            label nbrEdgei = edgei + size()/2;
-            half1Normals[edgei] = eN[nbrEdgei];
+            half1Normals[edgei] = eN[edgei + sizeby2];
 
             scalar magLe = mag(half0Normals[edgei]);
             scalar nbrMagLe = mag(half1Normals[edgei]);
-            scalar avLe = (magLe + nbrMagLe)/2.0;
+            scalar avLe = 0.5*(magLe + nbrMagLe);
 
             if (magLe < ROOTVSMALL && nbrMagLe < ROOTVSMALL)
             {
@@ -101,10 +100,10 @@ void Foam::cyclicFaPatch::calcTransforms()
         // Check for error in edge matching
         if (maxMatchError > matchTol_)
         {
-            label nbrEdgei = errorEdge + size()/2;
+            label nbrEdgei = errorEdge + sizeby2;
             scalar magLe = mag(half0Normals[errorEdge]);
             scalar nbrMagLe = mag(half1Normals[errorEdge]);
-            scalar avLe = (magLe + nbrMagLe)/2.0;
+            scalar avLe = 0.5*(magLe + nbrMagLe);
 
             FatalErrorInFunction
                 << "edge " << errorEdge
@@ -162,7 +161,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
 
     for (label edgei = 0; edgei < sizeby2; ++edgei)
     {
-        scalar avL = (magL[edgei] + magL[edgei + sizeby2])/2.0;
+        scalar avL = 0.5*(magL[edgei] + magL[edgei + sizeby2]);
 
         if
         (
@@ -191,7 +190,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
     // Check for error in matching
     if (maxMatchError > matchTol_)
     {
-        scalar avL = (magL[errorEdge] + magL[errorEdge + sizeby2])/2.0;
+        scalar avL = 0.5*(magL[errorEdge] + magL[errorEdge + sizeby2]);
 
         FatalErrorInFunction
             << "edge " << errorEdge << " and " << errorEdge + sizeby2
@@ -313,7 +312,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::transfer
 
     const label sizeby2 = this->size()/2;
 
-    for (label edgei=0; edgei<sizeby2; ++edgei)
+    for (label edgei=0; edgei < sizeby2; ++edgei)
     {
         pnf[edgei] = interfaceData[edgei + sizeby2];
         pnf[edgei + sizeby2] = interfaceData[edgei];
@@ -345,7 +344,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::internalFieldTransfer
 
     const label sizeby2 = this->size()/2;
 
-    for (label edgei=0; edgei<sizeby2; ++edgei)
+    for (label edgei=0; edgei < sizeby2; ++edgei)
     {
         pnf[edgei] = iF[edgeCells[edgei + sizeby2]];
         pnf[edgei + sizeby2] = iF[edgeCells[edgei]];
diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
index d257e4a199..4ac1c79f22 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
@@ -80,12 +80,7 @@ void Foam::processorFaPatch::makeNonGlobalPatchPoints() const
      || !boundaryMesh().mesh()().globalData().nGlobalPoints()
     )
     {
-        nonGlobalPatchPointsPtr_ = new labelList(nPoints());
-        labelList& ngpp = *nonGlobalPatchPointsPtr_;
-        forAll(ngpp, i)
-        {
-            ngpp[i] = i;
-        }
+        nonGlobalPatchPointsPtr_ = new labelList(identity(nPoints()));
     }
     else
     {
diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
index 5924eaae01..b21be3af93 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -177,19 +177,19 @@ public:
     // Member Functions
 
         //- Return interface size
-        virtual label interfaceSize() const
+        virtual label interfaceSize() const noexcept
         {
             return size();
         }
 
         //- Return processor number
-        int myProcNo() const
+        int myProcNo() const noexcept
         {
             return myProcNo_;
         }
 
         //- Return neighbour processor number
-        int neighbProcNo() const
+        int neighbProcNo() const noexcept
         {
             return neighbProcNo_;
         }
@@ -197,18 +197,11 @@ public:
         //- Return true if running parallel
         virtual bool coupled() const
         {
-            if (Pstream::parRun())
-            {
-                return true;
-            }
-            else
-            {
-                return false;
-            }
+            return Pstream::parRun();
         }
 
         //- Is this the master side?
-        virtual bool master() const
+        virtual bool master() const noexcept
         {
             return (myProcNo_ < neighbProcNo_);
         }
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index 81ddba7ee6..f6e8b98aae 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -210,88 +210,6 @@ const Foam::labelList& Foam::faPatch::pointLabels() const
 }
 
 
-void Foam::faPatch::calcPointLabels() const
-{
-    SLList<label> labels;
-
-    const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
-
-    forAll(edges, edgeI)
-    {
-        bool existStart = false;
-        bool existEnd = false;
-
-        forAllIters(labels, iter)
-        {
-            if (*iter == edges[edgeI].start())
-            {
-                existStart = true;
-            }
-
-            if (*iter == edges[edgeI].end())
-            {
-                existEnd = true;
-            }
-        }
-
-        if (!existStart)
-        {
-            labels.append(edges[edgeI].start());
-        }
-
-        if (!existEnd)
-        {
-            labels.append(edges[edgeI].end());
-        }
-    }
-
-    pointLabelsPtr_ = new labelList(labels);
-}
-
-
-void Foam::faPatch::calcPointEdges() const
-{
-    const labelList& points = pointLabels();
-
-    const edgeList::subList e = patchSlice(boundaryMesh().mesh().edges());
-
-    // set up storage for pointEdges
-    List<SLList<label>> pointEdgs(points.size());
-
-    forAll(e, edgeI)
-    {
-        const edge& curPoints = e[edgeI];
-
-        forAll(curPoints, pointI)
-        {
-            const label localPointIndex = points.find(curPoints[pointI]);
-
-            pointEdgs[localPointIndex].append(edgeI);
-        }
-    }
-
-    // sort out the list
-    pointEdgesPtr_ = new labelListList(pointEdgs.size());
-    labelListList& pEdges = *pointEdgesPtr_;
-
-    forAll(pointEdgs, pointI)
-    {
-        pEdges[pointI].setSize(pointEdgs[pointI].size());
-
-        label i = 0;
-        for
-        (
-            SLList<label>::iterator curEdgesIter = pointEdgs[pointI].begin();
-            curEdgesIter != pointEdgs[pointI].end();
-            ++curEdgesIter, ++i
-        )
-        {
-            pEdges[pointI][i] = curEdgesIter();
-        }
-    }
-}
-
-
 const Foam::labelListList& Foam::faPatch::pointEdges() const
 {
     if (!pointEdgesPtr_)
@@ -303,6 +221,94 @@ const Foam::labelListList& Foam::faPatch::pointEdges() const
 }
 
 
+void Foam::faPatch::calcPointLabels() const
+{
+    const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
+
+    // Walk boundary edges.
+    // The edge orientation corresponds to the face orientation
+    // (outwards normal).
+
+    // Note: could combine this with calcPointEdges for more efficiency
+
+    // Map<label> markedPoints(4*edges.size());
+    labelHashSet markedPoints(4*edges.size());
+    DynamicList<label> dynEdgePoints(2*edges.size());
+
+    for (const edge& e : edges)
+    {
+        // if (markedPoints.insert(e.first(), markedPoints.size()))
+        if (markedPoints.insert(e.first()))
+        {
+            dynEdgePoints.append(e.first());
+        }
+        // if (markedPoints.insert(e.second(), markedPoints.size()))
+        if (markedPoints.insert(e.second()))
+        {
+            dynEdgePoints.append(e.second());
+        }
+    }
+
+    // Transfer to plain list (reuse storage)
+    pointLabelsPtr_ = new labelList(std::move(dynEdgePoints));
+    /// const labelList& edgePoints = *pointLabelsPtr_;
+    ///
+    /// // Cannot use invertManyToMany - we have non-local edge numbering
+    ///
+    /// // Intermediate storage for pointEdges.
+    /// // Points on the boundary will normally connect 1 or 2 edges only.
+    /// List<DynamicList<label,2>> dynPointEdges(edgePoints.size());
+    ///
+    /// forAll(edges, edgei)
+    /// {
+    ///     const edge& e = edges[edgei];
+    ///
+    ///     dynPointEdges[markedPoints[e.first()]].append(edgei);
+    ///     dynPointEdges[markedPoints[e.second()]].append(edgei);
+    /// }
+    ///
+    /// // Flatten to regular list
+    /// pointEdgesPtr_ = new labelListList(edgePoints.size());
+    /// auto& pEdges = *pointEdgesPtr_;
+    ///
+    /// forAll(pEdges, pointi)
+    /// {
+    ///     pEdges[pointi] = std::move(dynPointEdges[pointi]);
+    /// }
+}
+
+
+void Foam::faPatch::calcPointEdges() const
+{
+    const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
+
+    const labelList& edgePoints = pointLabels();
+
+    // Cannot use invertManyToMany - we have non-local edge numbering
+
+    // Intermediate storage for pointEdges.
+    // Points on the boundary will normally connect 1 or 2 edges only.
+    List<DynamicList<label,2>> dynPointEdges(edgePoints.size());
+
+    forAll(edges, edgei)
+    {
+        const edge& e = edges[edgei];
+
+        dynPointEdges[edgePoints.find(e.first())].append(edgei);
+        dynPointEdges[edgePoints.find(e.second())].append(edgei);
+    }
+
+    // Flatten to regular list
+    pointEdgesPtr_ = new labelListList(edgePoints.size());
+    auto& pEdges = *pointEdgesPtr_;
+
+    forAll(pEdges, pointi)
+    {
+        pEdges[pointi] = std::move(dynPointEdges[pointi]);
+    }
+}
+
+
 Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 {
     if (nbrPolyPatchId_ < 0)
@@ -383,7 +389,7 @@ const Foam::scalarField& Foam::faPatch::magEdgeLengths() const
 
 Foam::tmp<Foam::vectorField> Foam::faPatch::edgeNormals() const
 {
-    tmp<vectorField> tedgeNorm(edgeLengths());
+    auto tedgeNorm = tmp<vectorField>::New(edgeLengths());
 
     for (vector& n : tedgeNorm.ref())
     {
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
index e252942118..71df3aa500 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
@@ -80,17 +80,14 @@ class faPatch
         //- Reference to boundary mesh
         const faBoundaryMesh& boundaryMesh_;
 
+        //- Demand-driven: edge-face addressing
+        mutable labelList::subList* edgeFacesPtr_;
 
-        // Demand-driven private data
+        //- Demand-driven: local points labels
+        mutable labelList* pointLabelsPtr_;
 
-            //- Edge-face addressing
-            mutable labelList::subList* edgeFacesPtr_;
-
-            //- Local points labels
-            mutable labelList* pointLabelsPtr_;
-
-            //- Point-edge addressing
-            mutable labelListList* pointEdgesPtr_;
+        //- Demand-driven: point-edge addressing
+        mutable labelListList* pointEdgesPtr_;
 
 
     // Private Member Functions