diff --git a/applications/test/foamCellZoneToVTK/Make/files b/applications/test/foamCellZoneToVTK/Make/files
new file mode 100644
index 0000000000..a4be1f0e87
--- /dev/null
+++ b/applications/test/foamCellZoneToVTK/Make/files
@@ -0,0 +1,3 @@
+foamCellZoneToVTK.C
+
+EXE = $(FOAM_USER_APPBIN)/foamCellZoneToVTK
diff --git a/applications/test/foamCellZoneToVTK/Make/options b/applications/test/foamCellZoneToVTK/Make/options
new file mode 100644
index 0000000000..d820876f44
--- /dev/null
+++ b/applications/test/foamCellZoneToVTK/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lmeshTools
diff --git a/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C b/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C
new file mode 100644
index 0000000000..a2f265db03
--- /dev/null
+++ b/applications/test/foamCellZoneToVTK/foamCellZoneToVTK.C
@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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/>.
+
+Application
+    foamCellZoneToVTK.C
+
+Description
+    Write tet-decomposed OpenFOAM mesh in VTK format.
+    For diagnostic purposes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "timeSelector.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "foamVtkInternalMeshWriter.H"
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Write OpenFOAM cellZone mesh to VTK"
+    );
+    argList::addOption
+    (
+        "cellZone",
+        "name",
+        "Convert mesh subset corresponding to specified cellZone"
+    );
+    argList::addBoolOption
+    (
+        "list",
+        "List names of cellZones and exit"
+    );
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+
+    word cellZoneName;
+    args.readIfPresent("cellZone", cellZoneName);
+
+    const bool optList = args.found("list");
+
+    if (optList)
+    {
+        if (!cellZoneName.empty())
+        {
+            Info<< "specify -list or -cellZone, but not both!" << nl;
+            return 1;
+        }
+    }
+    else if (cellZoneName.empty())
+    {
+        Info<< "Did not specify a cellZone!!" << nl;
+        return 1;
+    }
+
+    #include "createTime.H"
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+    fileName exportName("zonemesh-" + cellZoneName);
+
+    #include "createPolyMesh.H"
+
+    forAll(timeDirs, timei)
+    {
+        runTime.setTime(timeDirs[timei], timei);
+
+        polyMesh::readUpdateState state = mesh.readUpdate();
+
+        if (!timei || state != polyMesh::UNCHANGED)
+        {
+            fileName meshName(exportName);
+            if (state != polyMesh::UNCHANGED)
+            {
+                meshName += '_' + runTime.timeName();
+            }
+
+            if (optList)
+            {
+                Info<< "cellZones:" << nl;
+                for (const word& name : mesh.cellZones().names())
+                {
+                    Info<< "    " << name << nl;
+                }
+            }
+            else
+            {
+                const cellZone* zonePtr =
+                    mesh.cellZones().cfindZone(cellZoneName);
+
+                Info<< "cellZone " << cellZoneName;
+                if (!zonePtr)
+                {
+                    Info<< " ... not found" << nl;
+                    continue;
+                }
+                Info<< nl;
+
+                const cellZone& zn = *zonePtr;
+
+                // Define a subset
+                vtk::vtuCells vtuCells;
+                vtuCells.reset(mesh, zn);
+
+                vtk::internalMeshWriter writer
+                (
+                    mesh,
+                    vtuCells,
+                    fileName
+                    (
+                        mesh.time().globalPath() / meshName
+                    )
+                );
+
+                writer.writeGeometry();
+
+                writer.beginCellData();
+                writer.writeProcIDs();
+
+                Info<< "Wrote " << writer.output().name() << nl;
+            }
+        }
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/foamMeshToTet-vtk/Make/files b/applications/test/foamMeshToTet-vtk/Make/files
new file mode 100644
index 0000000000..576d08e8c1
--- /dev/null
+++ b/applications/test/foamMeshToTet-vtk/Make/files
@@ -0,0 +1,4 @@
+foamMeshToTet-vtk.C
+writeVTKtetMesh.C
+
+EXE = $(FOAM_USER_APPBIN)/foamMeshToTet-vtk
diff --git a/applications/test/foamMeshToTet-vtk/Make/options b/applications/test/foamMeshToTet-vtk/Make/options
new file mode 100644
index 0000000000..d820876f44
--- /dev/null
+++ b/applications/test/foamMeshToTet-vtk/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lmeshTools
diff --git a/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C b/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C
new file mode 100644
index 0000000000..d1d078949e
--- /dev/null
+++ b/applications/test/foamMeshToTet-vtk/foamMeshToTet-vtk.C
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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/>.
+
+Application
+    foamMeshToTet-vtk
+
+Description
+    Write tet-decomposed OpenFOAM mesh in VTK format.
+    For diagnostic purposes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "timeSelector.H"
+#include "Time.H"
+#include "polyMesh.H"
+
+namespace Foam
+{
+    void writeVTKtetMesh(const fileName& output, const polyMesh& mesh);
+}
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Write tet-decomposed OpenFOAM mesh in VTK"
+    );
+    argList::noParallel();
+    timeSelector::addOptions();
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+    fileName exportName = "tetmesh";
+    if (args.found("case"))
+    {
+        exportName += '-' + args.globalCaseName();
+    }
+
+    #include "createPolyMesh.H"
+
+    forAll(timeDirs, timei)
+    {
+        runTime.setTime(timeDirs[timei], timei);
+
+        polyMesh::readUpdateState state = mesh.readUpdate();
+
+        if (!timei || state != polyMesh::UNCHANGED)
+        {
+            fileName meshName(exportName);
+            if (state != polyMesh::UNCHANGED)
+            {
+                meshName += '_' + runTime.timeName();
+            }
+
+            writeVTKtetMesh(meshName, mesh);
+        }
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C b/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C
new file mode 100644
index 0000000000..ea7157fac0
--- /dev/null
+++ b/applications/test/foamMeshToTet-vtk/writeVTKtetMesh.C
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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 "polyMesh.H"
+#include "Fstream.H"
+#include "tetMatcher.H"
+#include "foamVtkInternalMeshWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+void writeVTKtetMesh(const fileName& output, const polyMesh& mesh_)
+{
+    #if WM_LABEL_SIZE == 64
+    # define FOAM_VTK_LABEL_NAME "vtktypeint64"
+    #else
+    # define FOAM_VTK_LABEL_NAME "int"
+    #endif
+
+    const faceList& faces = mesh_.faces();
+    const labelList& faceOwner = mesh_.faceOwner();
+
+    label nTets = 0;
+
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
+    {
+        nTets += 2 * faces[facei].nTriangles();
+    }
+    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
+    {
+        nTets += faces[facei].nTriangles();
+    }
+    const label nPoints = (mesh_.nPoints() + mesh_.nCells());
+
+    OFstream os(output + ".vtk");
+
+    Info<< "Write: " << os.name() << nl;
+
+    os  << "# vtk DataFile Version 5.1" << nl
+        << "tet-mesh" << nl
+        << "ASCII" << nl
+        << "DATASET UNSTRUCTURED_GRID" << nl
+        << nl;
+
+    os  << "POINTS " << nPoints << " float" << nl;
+    for (const point& p : mesh_.points())
+    {
+        os  << float(p.x()) << ' '
+            << float(p.y()) << ' '
+            << float(p.z()) << nl;
+    }
+    for (const point& p : mesh_.cellCentres())
+    {
+        os  << float(p.x()) << ' '
+            << float(p.y()) << ' '
+            << float(p.z()) << nl;
+    }
+    os  << nl;
+
+    os  << "CELLS " << (1 + nTets) << ' ' << (4 * nTets) << nl;
+
+    os  << "OFFSETS " << FOAM_VTK_LABEL_NAME << nl
+        << 0;  // begin offset = 0
+    {
+        label offset = 0;
+        for (label teti = 0; teti < nTets; ++teti)
+        {
+            offset += 4;
+            os << ' ' << offset;
+        }
+        os << nl << nl;
+    }
+
+    labelList nLocalTets(mesh_.nCells(), Zero);
+
+    os  << nl
+        << "CONNECTIVITY " << FOAM_VTK_LABEL_NAME << nl;
+
+    for (label celli = 0; celli < mesh_.nCells(); ++celli)
+    {
+        const cell& cFaces = mesh_.cells()[celli];
+
+        if (tetMatcher::test(mesh_, celli))
+        {
+            // Tet: no cell-centre decomposition
+
+            const label facei = cFaces[0];
+            const face& f0 = faces[facei];
+
+            // Get the other point from f1. Tbd: check if not duplicate face
+            // (ACMI / ignoreBoundaryFaces_).
+            const face& f1 = faces[cFaces[1]];
+            label apexi = -1;
+            forAll(f1, fp)
+            {
+                apexi = f1[fp];
+                if (!f0.found(apexi))
+                {
+                    break;
+                }
+            }
+
+            const label p0 = f0[0];
+            label p1 = f0[1];
+            label p2 = f0[2];
+
+            if (faceOwner[facei] == celli)
+            {
+                std::swap(p1, p2);
+            }
+
+            ++nLocalTets[celli];
+            os << p0 << ' ' << p1 << ' ' << p2 << ' ' << apexi << nl;
+        }
+        else
+        {
+            for (const label facei : cFaces)
+            {
+                const face& f = faces[facei];
+
+                label fp0 = mesh_.tetBasePtIs()[facei];
+
+                // Fallback
+                if (fp0 < 0)
+                {
+                    fp0 = 0;
+                }
+
+                const label p0 = f[fp0];
+                label fp = f.fcIndex(fp0);
+                for (label i = 2; i < f.size(); ++i)
+                {
+                    label p1 = f[fp];
+                    fp = f.fcIndex(fp);
+                    label p2 = f[fp];
+
+                    if (faceOwner[facei] == celli)
+                    {
+                        std::swap(p1, p2);
+                    }
+
+                    ++nLocalTets[celli];
+                    os  << p0 << ' ' << p1 << ' ' << p2 << ' '
+                        << (mesh_.nPoints() + celli) << nl;
+                }
+            }
+        }
+    }
+
+    os  << nl
+        << "CELL_TYPES " << nTets << nl;
+
+    for (label teti = 0; teti < nTets; ++teti)
+    {
+        if (teti) os  << ' ';
+        os  << vtk::cellType::VTK_TETRA;
+    }
+    os << nl;
+
+    os << nl << "CELL_DATA " << nTets << nl
+       << "FIELD FieldData " << 1 << nl;
+
+    os << "cellID " << 1 << ' ' << nTets << " int" << nl;
+    forAll(nLocalTets, celli)
+    {
+        label n = nLocalTets[celli];
+        while (n--)
+        {
+            os << ' ' << celli;
+        }
+        os << nl;
+    }
+}
+
+} // End namespace Foam
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
index c98ba86456..add182f7ec 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
@@ -264,7 +264,7 @@ Description
         // Begin CellData
         if (internalWriter)
         {
-            // Optionally with cellID and procID fields
+            // Optionally with (cellID, procID) fields
             internalWriter->beginCellData
             (
                 (withMeshIds ? 1 + (internalWriter->parallel() ? 1 : 0) : 0)
@@ -278,20 +278,21 @@ Description
             }
         }
 
-        if (nVolFields)
+        if (nVolFields || withMeshIds)
         {
             for (vtk::patchWriter& writer : patchWriters)
             {
-                // Optionally with patchID field
+                // Optionally with (patchID, procID) fields
                 writer.beginCellData
                 (
-                    (withMeshIds ? 1 : 0)
+                    (withMeshIds ? 2 : 0)
                   + nVolFields
                 );
 
                 if (withMeshIds)
                 {
                     writer.writePatchIDs();
+                    writer.writeProcIDs();
                 }
             }
         }
diff --git a/src/conversion/vtk/output/foamVtkPatchWriterTemplates.C b/src/conversion/vtk/output/foamVtkPatchWriterTemplates.C
index 9dc41948e4..0219a3b835 100644
--- a/src/conversion/vtk/output/foamVtkPatchWriterTemplates.C
+++ b/src/conversion/vtk/output/foamVtkPatchWriterTemplates.C
@@ -149,7 +149,7 @@ void Foam::vtk::patchWriter::write
             << exit(FatalError);
     }
 
-    label nFaces = nLocalFaces_;
+    label nFaces = nLocalPolys_;
 
     if (parallel_)
     {
diff --git a/src/fileFormats/Make/files b/src/fileFormats/Make/files
index 2686cefc43..f7fd57b853 100644
--- a/src/fileFormats/Make/files
+++ b/src/fileFormats/Make/files
@@ -59,6 +59,7 @@ vtk/part/foamVtuCells.C
 vtk/part/foamVtuSizing.C
 
 vtk/read/vtkUnstructuredReader.C
+vtk/write/foamVtkLineWriter.C
 vtk/write/foamVtkPolyWriter.C
 vtk/write/foamVtkSurfaceWriter.C
 
diff --git a/src/fileFormats/vtk/base/foamVtkCore.C b/src/fileFormats/vtk/base/foamVtkCore.C
index 14f832fec6..b9488e8824 100644
--- a/src/fileFormats/vtk/base/foamVtkCore.C
+++ b/src/fileFormats/vtk/base/foamVtkCore.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,6 +113,8 @@ Foam::vtk::dataArrayAttrNames
 });
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 // Legacy
 
 const Foam::word Foam::vtk::legacy::fileExtension("vtk");
@@ -132,11 +134,27 @@ const Foam::Enum
 <
     Foam::vtk::fileTag
 >
-Foam::vtk::legacy::dataTypeNames
+Foam::vtk::legacy::fileTagNames
 ({
-    { vtk::fileTag::CELL_DATA,  "CELL_DATA" },
+    { vtk::fileTag::POINTS, "POINTS" },
+    { vtk::fileTag::CELLS, "CELLS" },
+    { vtk::fileTag::POLYS, "POLYGONS" },
+    { vtk::fileTag::VERTS, "VERTICES" },
+    { vtk::fileTag::LINES, "LINES" },
+    { vtk::fileTag::CELL_DATA, "CELL_DATA" },
     { vtk::fileTag::POINT_DATA, "POINT_DATA" },
 });
 
 
+const Foam::Enum
+<
+    Foam::vtk::dataArrayAttr
+>
+Foam::vtk::legacy::dataArrayAttrNames
+({
+    { vtk::dataArrayAttr::OFFSETS, "OFFSETS" },
+    { vtk::dataArrayAttr::CONNECTIVITY, "CONNECTIVITY" },
+});
+
+
 // ************************************************************************* //
diff --git a/src/fileFormats/vtk/base/foamVtkCore.H b/src/fileFormats/vtk/base/foamVtkCore.H
index 2524964a0a..08615850ec 100644
--- a/src/fileFormats/vtk/base/foamVtkCore.H
+++ b/src/fileFormats/vtk/base/foamVtkCore.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -183,8 +183,12 @@ namespace legacy
     //- Legacy content names (POLYDATA, UNSTRUCTURED_GRID)
     extern const Foam::Enum<vtk::fileTag> contentNames;
 
-    //- Legacy data type names (CELL_DATA, POINT_DATA)
-    extern const Foam::Enum<vtk::fileTag> dataTypeNames;
+    //- Legacy file tags (eg, LINES, CELL_DATA, POINT_DATA, ...)
+    extern const Foam::Enum<vtk::fileTag> fileTagNames;
+
+    //- Legacy attributes (eg, OFFSETS)
+    extern const Foam::Enum<dataArrayAttr> dataArrayAttrNames;
+
 
 } // End namespace legacy
 
diff --git a/src/fileFormats/vtk/file/foamVtkFileWriter.C b/src/fileFormats/vtk/file/foamVtkFileWriter.C
index aff45b9017..3cc8b7e72a 100644
--- a/src/fileFormats/vtk/file/foamVtkFileWriter.C
+++ b/src/fileFormats/vtk/file/foamVtkFileWriter.C
@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "foamVtkFileWriter.H"
+#include "globalIndex.H"
 #include "OSspecific.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -342,7 +343,7 @@ bool Foam::vtk::fileWriter::open(const fileName& file, bool parallel)
 
     // Open a file and attach a formatter
     // - on master (always)
-    // - on slave if not parallel
+    // - on subproc (if not parallel)
     //
     // This means we can always check if format_ is defined to know if output
     // is desired on any particular process.
@@ -515,7 +516,7 @@ void Foam::vtk::fileWriter::writeTimeValue(scalar timeValue)
             << exit(FatalError);
     }
 
-    // No collectives - can skip on slave processors
+    // No collectives - can skip on sub-procs
     if (!format_) return;
 
     if (legacy())
@@ -529,4 +530,66 @@ void Foam::vtk::fileWriter::writeTimeValue(scalar timeValue)
 }
 
 
+bool Foam::vtk::fileWriter::writeProcIDs(const label nValues)
+{
+    // Write procIDs whenever running in parallel
+
+    if (!Pstream::parRun())
+    {
+        return false;  // Non-parallel: skip
+    }
+
+    if (isState(outputState::CELL_DATA))
+    {
+        ++nCellData_;
+    }
+    else
+    {
+        reportBadState(FatalErrorInFunction, outputState::CELL_DATA)
+            << " for procID field" << nl << endl
+            << exit(FatalError);
+
+        return false;
+    }
+
+
+    const globalIndex procSizes
+    (
+        parallel_
+      ? globalIndex(nValues)
+      : globalIndex()
+    );
+
+    const label totalCount = (parallel_ ? procSizes.size() : nValues);
+
+    this->beginDataArray<label>("procID", totalCount);
+
+    bool good = false;
+
+    if (parallel_)
+    {
+        if (Pstream::master())
+        {
+            // Per-processor ids
+            for (const int proci : Pstream::allProcs())
+            {
+                vtk::write(format(), label(proci), procSizes.localSize(proci));
+            }
+            good = true;
+        }
+    }
+    else
+    {
+        vtk::write(format(), label(Pstream::myProcNo()), totalCount);
+        good = true;
+    }
+
+
+    this->endDataArray();
+
+    // MPI barrier
+    return parallel_ ? returnReduce(good, orOp<bool>()) : good;
+}
+
+
 // ************************************************************************* //
diff --git a/src/fileFormats/vtk/file/foamVtkFileWriter.H b/src/fileFormats/vtk/file/foamVtkFileWriter.H
index 19a1d4e150..ea25b30323 100644
--- a/src/fileFormats/vtk/file/foamVtkFileWriter.H
+++ b/src/fileFormats/vtk/file/foamVtkFileWriter.H
@@ -107,10 +107,10 @@ protected:
         //- The output file name
         fileName outputFile_;
 
-        //- The VTK formatter in use (master process)
+        //- The VTK formatter in use (only valid on master process)
         autoPtr<vtk::formatter> format_;
 
-        //- The backend ostream in use (master process)
+        //- The backend ostream in use (only opened on master process)
         std::ofstream os_;
 
 
@@ -123,16 +123,16 @@ protected:
         Ostream& reportBadState(Ostream&, outputState, outputState) const;
 
         //- The backend ostream in use
-        inline std::ofstream& os();
+        inline std::ofstream& os() noexcept;
 
         //- The VTK formatter in use
         inline vtk::formatter& format();
 
-        //- True if the output state corresponds to the test state.
-        inline bool isState(outputState test) const;
+        //- True if output state corresponds to the test state.
+        inline bool isState(outputState test) const noexcept;
 
-        //- True if the output state does not correspond to the test state.
-        inline bool notState(outputState test) const;
+        //- True if output state does not correspond to the test state.
+        inline bool notState(outputState test) const noexcept;
 
         //- Start of a field or DataArray output (legacy or non-legacy).
         template<class Type>
@@ -151,26 +151,6 @@ protected:
         //- End of a POINTS DataArray
         void endPoints();
 
-
-        //- Write uniform field content.
-        //  No context checking (eg, file-open, CellData, PointData, etc)
-        template<class Type>
-        void writeUniform
-        (
-            const word& fieldName,
-            const Type& val,
-            const label nValues
-        );
-
-        //- Write basic (primitive) field content
-        //  No context checking (eg, file-open, CellData, PointData, etc)
-        template<class Type>
-        void writeBasicField
-        (
-            const word& fieldName,
-            const UList<Type>& field
-        );
-
         //- Trigger change state to Piece. Resets nCellData_, nPointData_.
         bool enter_Piece();
 
@@ -194,6 +174,34 @@ protected:
         bool exit_File();
 
 
+    // Field writing
+
+        //- Write uniform field content.
+        //  No context checking (eg, file-open, CellData, PointData, etc)
+        //  The value and count can be different on each processor
+        template<class Type>
+        void writeUniform
+        (
+            const word& fieldName,
+            const Type& val,
+            const label nValues
+        );
+
+        //- Write basic (primitive) field content
+        //  No context checking (eg, file-open, CellData, PointData, etc)
+        template<class Type>
+        void writeBasicField
+        (
+            const word& fieldName,
+            const UList<Type>& field
+        );
+
+        //- Write nValues of processor ids as CellData (no-op in serial)
+        bool writeProcIDs(const label nValues);
+
+
+    // Other
+
         //- No copy construct
         fileWriter(const fileWriter&) = delete;
 
diff --git a/src/fileFormats/vtk/file/foamVtkFileWriterI.H b/src/fileFormats/vtk/file/foamVtkFileWriterI.H
index 62c7d86600..192eac45a7 100644
--- a/src/fileFormats/vtk/file/foamVtkFileWriterI.H
+++ b/src/fileFormats/vtk/file/foamVtkFileWriterI.H
@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-inline std::ofstream& Foam::vtk::fileWriter::os()
+inline std::ofstream& Foam::vtk::fileWriter::os() noexcept
 {
     return os_;
 }
@@ -39,13 +39,13 @@ inline Foam::vtk::formatter& Foam::vtk::fileWriter::format()
 }
 
 
-inline bool Foam::vtk::fileWriter::isState(outputState test) const
+inline bool Foam::vtk::fileWriter::isState(outputState test) const noexcept
 {
     return (test == state_);
 }
 
 
-inline bool Foam::vtk::fileWriter::notState(outputState test) const
+inline bool Foam::vtk::fileWriter::notState(outputState test) const noexcept
 {
     return (test != state_);
 }
diff --git a/src/fileFormats/vtk/file/foamVtkFileWriterTemplates.C b/src/fileFormats/vtk/file/foamVtkFileWriterTemplates.C
index 79a774fce9..21ec488823 100644
--- a/src/fileFormats/vtk/file/foamVtkFileWriterTemplates.C
+++ b/src/fileFormats/vtk/file/foamVtkFileWriterTemplates.C
@@ -89,14 +89,25 @@ void Foam::vtk::fileWriter::writeUniform
 (
     const word& fieldName,
     const Type& val,
-    const label nValues
+    const label nLocalValues
 )
 {
-    this->beginDataArray<Type>(fieldName, nValues);
+    label nTotal = nLocalValues;
 
-    if (format_)
+    if (parallel_)
     {
-        vtk::write(format(), val, nValues);
+        reduce(nTotal, sumOp<label>());
+    }
+
+    this->beginDataArray<Type>(fieldName, nTotal);
+
+    if (parallel_)
+    {
+        vtk::writeValueParallel(format_.ref(), val, nLocalValues);
+    }
+    else
+    {
+        vtk::write(format(), val, nLocalValues);
     }
 
     this->endDataArray();
diff --git a/src/fileFormats/vtk/output/foamVtkOutput.C b/src/fileFormats/vtk/output/foamVtkOutput.C
index 7e32b87dcf..08a00cf6e4 100644
--- a/src/fileFormats/vtk/output/foamVtkOutput.C
+++ b/src/fileFormats/vtk/output/foamVtkOutput.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,7 +39,7 @@ License
 #include "globalIndex.H"
 #include "instant.H"
 #include "Fstream.H"
-#include "Pstream.H"
+#include "PstreamBuffers.H"
 #include "OSspecific.H"
 
 // * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * * //
@@ -99,7 +99,7 @@ void Foam::vtk::writeIdentity
     label start
 )
 {
-    // No nComponents for label, so use fmt.write() directly
+    // No nComponents for label, can use fmt.write() directly
     for (label i=0; i < len; ++i)
     {
         fmt.write(start);
@@ -114,7 +114,7 @@ void Foam::vtk::writeList
     const UList<uint8_t>& values
 )
 {
-    // No nComponents for char, so use fmt.write() directly
+    // No nComponents for char, can use fmt.write() directly
     for (const uint8_t val : values)
     {
         fmt.write(val);
@@ -122,42 +122,6 @@ void Foam::vtk::writeList
 }
 
 
-void Foam::vtk::writeListParallel
-(
-    vtk::formatter& fmt,
-    const UList<uint8_t>& values
-)
-{
-    if (Pstream::master())
-    {
-        vtk::writeList(fmt, values);
-
-        List<uint8_t> recv;
-
-        // Receive and write
-        for (const int slave : Pstream::subProcs())
-        {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
-
-            fromSlave >> recv;
-
-            vtk::writeList(fmt, recv);
-        }
-    }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
-
-        toMaster << values;
-    }
-}
-
-
 void Foam::vtk::writeListParallel
 (
     vtk::formatter& fmt,
@@ -165,45 +129,55 @@ void Foam::vtk::writeListParallel
     const globalIndex& procOffset
 )
 {
+    // List sizes
+    const globalIndex sizes(values.size());
+
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
+    {
+        UOPstream os(Pstream::masterNo(), pBufs);
+        os.write
+        (
+            reinterpret_cast<const char*>(values.cdata()),
+            values.size_bytes()
+        );
+    }
+
+    pBufs.finishedSends();
+
     if (Pstream::master())
     {
+        // Master data
+
         // Write with offset
         const label offsetId = procOffset.offset(0);
-
         for (const label val : values)
         {
             vtk::write(fmt, val + offsetId);
         }
 
-        labelList recv;
-
         // Receive and write
-        for (const int slave : Pstream::subProcs())
+        for (const int proci : Pstream::subProcs())
         {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+            List<label> recv(sizes.localSize(proci));
 
-            fromSlave >> recv;
-
-            const label offsetId = procOffset.offset(slave);
+            UIPstream is(proci, pBufs);
+            is.read
+            (
+                reinterpret_cast<char*>(recv.data()),
+                recv.size_bytes()
+            );
 
             // Write with offset
+            const label offsetId = procOffset.offset(proci);
             for (const label val : recv)
             {
                 vtk::write(fmt, val + offsetId);
             }
         }
     }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
-
-        toMaster << values;
-    }
 }
 
 
@@ -219,6 +193,9 @@ void Foam::vtk::legacy::fileHeader
     // Line 1:
     os  << "# vtk DataFile Version 2.0" << nl;
 
+    // OR
+    // os  << "# vtk DataFile Version 5.1" << nl;
+
     // Line 2: title
 
     const auto truncate = title.find('\n');
diff --git a/src/fileFormats/vtk/output/foamVtkOutput.H b/src/fileFormats/vtk/output/foamVtkOutput.H
index ac8c3b6852..ad7d6bf774 100644
--- a/src/fileFormats/vtk/output/foamVtkOutput.H
+++ b/src/fileFormats/vtk/output/foamVtkOutput.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,7 +59,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class instant;
 class globalIndex;
 
@@ -94,14 +94,20 @@ namespace vtk
     //  The output does not include the payload size.
     void writeList(vtk::formatter& fmt, const UList<uint8_t>& values);
 
-    //- Write a list of uint8_t values.
-    //  The output does not include the payload size.
-    void writeListParallel(vtk::formatter& fmt, const UList<uint8_t>& values);
-
     //- Component-wise write of a value (N times)
     template<class Type>
     inline void write(vtk::formatter& fmt, const Type& val, const label n=1);
 
+    //- Component-wise write of a value (N times) in parallel
+    //  The value and count may differ on each processor
+    template<class Type>
+    inline void writeValueParallel
+    (
+        vtk::formatter& fmt,
+        const Type& val,
+        const label count=1
+    );
+
 
     //- Write a list of values.
     //  The output does not include the payload size.
@@ -215,7 +221,7 @@ namespace legacy
 
 // Functions
 
-    //- Emit header for legacy file.
+    //- Emit header for legacy file (vtk DataFile Version 2.0)
     //  Writes "ASCII" or "BINARY" depending on specified type.
     void fileHeader(std::ostream& os, const std::string& title, bool binary);
 
@@ -248,6 +254,12 @@ namespace legacy
     //- Emit header for POINTS (with trailing newline).
     inline void beginPoints(std::ostream& os, label nPoints);
 
+    //- Emit header for LINES (with trailing newline).
+    //  The nConnectivity is the sum of all connectivity points used,
+    //  but \b without additional space for the size prefixes.
+    //  The additional prefix sizes are added internally.
+    inline void beginLines(std::ostream& os, label nLines, label nConnectivity);
+
     //- Emit header for POLYGONS (with trailing newline).
     //  The nConnectivity is the sum of all connectivity points used,
     //  but \b without additional space for the size prefixes.
diff --git a/src/fileFormats/vtk/output/foamVtkOutputI.H b/src/fileFormats/vtk/output/foamVtkOutputI.H
index 036432e3c0..9789eecc41 100644
--- a/src/fileFormats/vtk/output/foamVtkOutputI.H
+++ b/src/fileFormats/vtk/output/foamVtkOutputI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,7 +113,23 @@ inline void Foam::vtk::legacy::fileHeader
 inline void Foam::vtk::legacy::beginPoints(std::ostream& os, label nPoints)
 {
     os  << nl
-        << "POINTS " << nPoints << " float" << nl;
+        << legacy::fileTagNames[vtk::fileTag::POINTS]
+        << ' ' << nPoints
+        << " float" << nl;
+}
+
+
+inline void Foam::vtk::legacy::beginLines
+(
+    std::ostream& os,
+    label nLines,
+    label nConnectivity
+)
+{
+    os  << nl
+        << legacy::fileTagNames[vtk::fileTag::LINES]
+        << ' ' << nLines
+        << ' ' << (nLines + nConnectivity) << nl;
 }
 
 
@@ -125,7 +141,9 @@ inline void Foam::vtk::legacy::beginPolys
 )
 {
     os  << nl
-        << "POLYGONS " << nPolys << ' ' << (nPolys + nConnectivity) << nl;
+        << legacy::fileTagNames[vtk::fileTag::POLYS]
+        << ' ' << nPolys
+        << ' ' << (nPolys + nConnectivity) << nl;
 }
 
 
@@ -159,7 +177,7 @@ inline void Foam::vtk::legacy::beginCellData
 {
     fmt.os()
         << nl
-        << legacy::dataTypeNames[vtk::fileTag::CELL_DATA]
+        << legacy::fileTagNames[vtk::fileTag::CELL_DATA]
         << ' ' << nCells << nl;
     legacy::fieldData(fmt, nFields);
 }
@@ -174,7 +192,7 @@ inline void Foam::vtk::legacy::beginPointData
 {
     fmt.os()
         << nl
-        << legacy::dataTypeNames[vtk::fileTag::POINT_DATA]
+        << legacy::fileTagNames[vtk::fileTag::POINT_DATA]
         << ' ' << nPoints << nl;
     legacy::fieldData(fmt, nFields);
 }
diff --git a/src/fileFormats/vtk/output/foamVtkOutputTemplates.C b/src/fileFormats/vtk/output/foamVtkOutputTemplates.C
index e296a2b2b5..59bdb0e88f 100644
--- a/src/fileFormats/vtk/output/foamVtkOutputTemplates.C
+++ b/src/fileFormats/vtk/output/foamVtkOutputTemplates.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,7 +25,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Pstream.H"
+#include "globalIndex.H"
+#include "PstreamBuffers.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -122,6 +123,44 @@ void Foam::vtk::writeLists
 }
 
 
+template<class Type>
+void Foam::vtk::writeValueParallel
+(
+    vtk::formatter& fmt,
+    const Type& val,
+    const label count
+)
+{
+    if (Pstream::master())
+    {
+        vtk::write(fmt, val, count);
+
+        label subCount;
+        Type subValue;
+
+        // Receive each [size, value] tuple
+        for (const int proci : Pstream::subProcs())
+        {
+            IPstream is(Pstream::commsTypes::blocking, proci);
+            is >> subCount >> subValue;
+
+            vtk::write(fmt, subValue, subCount);
+        }
+    }
+    else
+    {
+        OPstream os
+        (
+            Pstream::commsTypes::blocking,
+            Pstream::masterNo()
+        );
+
+        // Send [size, value] tuple
+        os << count << val;
+    }
+}
+
+
 template<class Type>
 void Foam::vtk::writeListParallel
 (
@@ -129,32 +168,59 @@ void Foam::vtk::writeListParallel
     const UList<Type>& values
 )
 {
-    if (Pstream::master())
+    // List sizes
+    const globalIndex sizes(values.size());
+
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
     {
-        vtk::writeList(fmt, values);
-
-        List<Type> recv;
-
-        // Receive and write
-        for (const int slave : Pstream::subProcs())
+        UOPstream os(Pstream::masterNo(), pBufs);
+        if (is_contiguous<Type>::value)
         {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
-
-            fromSlave >> recv;
-
-            vtk::writeList(fmt, recv);
+            os.write
+            (
+                reinterpret_cast<const char*>(values.cdata()),
+                values.size_bytes()
+            );
+        }
+        else
+        {
+            os << values;
         }
     }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
 
-        toMaster << values;
+    pBufs.finishedSends();
+
+    if (Pstream::master())
+    {
+        // Write master data
+        vtk::writeList(fmt, values);
+
+        // Receive and write
+        for (const int proci : Pstream::subProcs())
+        {
+            UIPstream is(proci, pBufs);
+
+            {
+                List<Type> recv(sizes.localSize(proci));
+
+                if (is_contiguous<Type>::value)
+                {
+                    is.read
+                    (
+                        reinterpret_cast<char*>(recv.data()),
+                        recv.size_bytes()
+                    );
+                }
+                else
+                {
+                    is >> recv;
+                }
+                vtk::writeList(fmt, recv);
+            }
+        }
     }
 }
 
@@ -167,33 +233,39 @@ void Foam::vtk::writeListParallel
     const labelUList& addressing
 )
 {
+    UIndirectList<Type> send(values, addressing);
+
+    // List sizes
+    const globalIndex sizes(send.size());
+
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
+    {
+        UOPstream os(Pstream::masterNo(), pBufs);
+        os << send;
+    }
+
+    pBufs.finishedSends();
+
     if (Pstream::master())
     {
+        // Write master data
         vtk::writeList(fmt, values, addressing);
 
-        List<Type> recv;
-
         // Receive and write
-        for (const int slave : Pstream::subProcs())
+        for (const int proci : Pstream::subProcs())
         {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+            UIPstream is(proci, pBufs);
 
-            fromSlave >> recv;
-
-            vtk::writeList(fmt, recv);
+            {
+                List<Type> recv;
+                is >> recv;
+                vtk::writeList(fmt, recv);
+            }
         }
     }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
-
-        toMaster << List<Type>(values, addressing);
-    }
 }
 
 
@@ -205,32 +277,66 @@ void Foam::vtk::writeListParallel
     const bitSet& selected
 )
 {
-    if (Pstream::master())
+    List<Type> send;
+    if (!Pstream::master())
     {
-        vtk::writeList(fmt, values, selected);
+        send = subset(selected, values);
+    }
 
-        List<Type> recv;
+    // List sizes.
+    // NOTE okay to skip proc0 since we only need sizes (not offsets)
+    const globalIndex sizes(send.size());
 
-        // Receive and write
-        for (const int slave : Pstream::subProcs())
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
+    {
+        UOPstream os(Pstream::masterNo(), pBufs);
+        if (is_contiguous<Type>::value)
         {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
-
-            fromSlave >> recv;
-
-            vtk::writeList(fmt, recv);
+            os.write
+            (
+                reinterpret_cast<const char*>(send.cdata()),
+                send.size_bytes()
+            );
+        }
+        else
+        {
+            os << send;
         }
     }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
 
-        toMaster << subset(selected, values);
+    pBufs.finishedSends();
+
+    if (Pstream::master())
+    {
+        // Write master data
+        vtk::writeList(fmt, values, selected);
+
+        // Receive and write
+        for (const int proci : Pstream::subProcs())
+        {
+            UIPstream is(proci, pBufs);
+
+            {
+                List<Type> recv(sizes.localSize(proci));
+
+                if (is_contiguous<Type>::value)
+                {
+                    is.read
+                    (
+                        reinterpret_cast<char*>(recv.data()),
+                        recv.size_bytes()
+                    );
+                }
+                else
+                {
+                    is >> recv;
+                }
+                vtk::writeList(fmt, recv);
+            }
+        }
     }
 }
 
@@ -243,34 +349,90 @@ void Foam::vtk::writeListsParallel
     const UList<Type>& values2
 )
 {
+    // List sizes
+    const globalIndex sizes1(values1.size());
+    const globalIndex sizes2(values2.size());
+
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
+    {
+        UOPstream os(Pstream::masterNo(), pBufs);
+        if (is_contiguous<Type>::value)
+        {
+            os.write
+            (
+                reinterpret_cast<const char*>(values1.cdata()),
+                values1.size_bytes()
+            );
+            os.write
+            (
+                reinterpret_cast<const char*>(values2.cdata()),
+                values2.size_bytes()
+            );
+        }
+        else
+        {
+            os << values1 << values2;
+        }
+    }
+
+    pBufs.finishedSends();
+
     if (Pstream::master())
     {
+        // Write master data
         vtk::writeList(fmt, values1);
         vtk::writeList(fmt, values2);
 
-        List<Type> recv1, recv2;
-
-        // Receive and write
-        for (const int slave : Pstream::subProcs())
-        {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
-
-            fromSlave >> recv1 >> recv2;
-
-            vtk::writeList(fmt, recv1);
-            vtk::writeList(fmt, recv2);
-        }
-    }
-    else
-    {
-        // Send to master
-        OPstream toMaster
+        // Reserve max receive size
+        DynamicList<Type> recv
         (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
+            max(sizes1.maxNonLocalSize(), sizes2.maxNonLocalSize())
         );
 
-        toMaster << values1 << values2;
+        // Receive and write
+        for (const int proci : Pstream::subProcs())
+        {
+            UIPstream is(proci, pBufs);
+
+            // values1
+            {
+                List<Type> recv(sizes1.localSize(proci));
+                if (is_contiguous<Type>::value)
+                {
+                    is.read
+                    (
+                        reinterpret_cast<char*>(recv.data()),
+                        recv.size_bytes()
+                    );
+                }
+                else
+                {
+                    is >> recv;
+                }
+                vtk::writeList(fmt, recv);
+            }
+
+            // values2
+            {
+                List<Type> recv(sizes2.localSize(proci));
+                if (is_contiguous<Type>::value)
+                {
+                    is.read
+                    (
+                        reinterpret_cast<char*>(recv.data()),
+                        recv.size_bytes()
+                    );
+                }
+                else
+                {
+                    is >> recv;
+                }
+                vtk::writeList(fmt, recv);
+            }
+        }
     }
 }
 
@@ -284,35 +446,45 @@ void Foam::vtk::writeListsParallel
     const labelUList& addressing
 )
 {
+    UIndirectList<Type> send2(values2, addressing);
+
+    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+    // Send to master
+    if (!Pstream::master())
+    {
+        UOPstream os(Pstream::masterNo(), pBufs);
+        os << values1 << send2;
+    }
+
+    pBufs.finishedSends();
+
     if (Pstream::master())
     {
+        // Write master data
         vtk::writeList(fmt, values1);
         vtk::writeList(fmt, values2, addressing);
 
-        List<Type> recv1, recv2;
-
         // Receive and write
-        for (const int slave : Pstream::subProcs())
+        for (const int proci : Pstream::subProcs())
         {
-            IPstream fromSlave(Pstream::commsTypes::blocking, slave);
+            UIPstream is(proci, pBufs);
 
-            fromSlave >> recv1 >> recv2;
+            // values1
+            {
+                List<Type> recv;
+                is >> recv;
+                vtk::writeList(fmt, recv);
+            }
 
-            vtk::writeList(fmt, recv1);
-            vtk::writeList(fmt, recv2);
+            // values2 (send2)
+            {
+                List<Type> recv;
+                is >> recv;
+                vtk::writeList(fmt, recv);
+            }
         }
     }
-    else
-    {
-        // Send to master
-        OPstream toMaster
-        (
-            Pstream::commsTypes::blocking,
-            Pstream::masterNo()
-        );
-
-        toMaster << values1 << List<Type>(values2, addressing);
-    }
 }
 
 
diff --git a/src/fileFormats/vtk/part/foamVtkMeshMaps.H b/src/fileFormats/vtk/part/foamVtkMeshMaps.H
index c2c54cc049..474c2d158b 100644
--- a/src/fileFormats/vtk/part/foamVtkMeshMaps.H
+++ b/src/fileFormats/vtk/part/foamVtkMeshMaps.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -75,11 +75,8 @@ public:
 
     // Constructors
 
-        //- Default construct: zero-sized, no reserved size
-        inline foamVtkMeshMaps();
-
-        //- Construct with reserved size
-        inline explicit foamVtkMeshMaps(const label size);
+        //- Default construct: zero-sized, no reserved sizes
+        foamVtkMeshMaps() = default;
 
 
     // Member Functions
@@ -87,26 +84,26 @@ public:
     // Access
 
         //- Original cell ids for all cells (regular and decomposed).
-        //  A regular mesh comprising only primitive cell types, this will just
+        //  For regular mesh comprising only primitive cell types, this will
         //  be an identity list. However, for subsetted meshes and decomposed
         //  cells this becomes a useful means of mapping from the original mesh.
-        inline const labelList& cellMap() const;
+        inline const labelList& cellMap() const noexcept;
 
         //- Write access to original cell ids
-        inline DynamicList<label>& cellMap();
+        inline DynamicList<label>& cellMap() noexcept;
 
         //- Point labels for subsetted meshes
-        inline const labelList& pointMap() const;
+        inline const labelList& pointMap() const noexcept;
 
         //- Write access to point labels for subsetted meshes
-        inline DynamicList<label>& pointMap();
+        inline DynamicList<label>& pointMap() noexcept;
 
         //- Any additional (user) labels.
         //  Eg, cell-centre labels for additional points of decomposed cells
-        inline const labelList& additionalIds() const;
+        inline const labelList& additionalIds() const noexcept;
 
         //- Write access to additional (user) labels.
-        inline DynamicList<label>& additionalIds();
+        inline DynamicList<label>& additionalIds() noexcept;
 
 
     // Edit
@@ -115,7 +112,7 @@ public:
         inline void clear();
 
         //- Renumber cell ids (cellMap and additionalIds) to account for
-        //  subset meshes
+        //- subset meshes
         void renumberCells(const labelUList& mapping);
 
         //- Renumber point ids (pointMap) to account for subset meshes
diff --git a/src/fileFormats/vtk/part/foamVtkMeshMapsI.H b/src/fileFormats/vtk/part/foamVtkMeshMapsI.H
index 3934a7f946..2e39034ab1 100644
--- a/src/fileFormats/vtk/part/foamVtkMeshMapsI.H
+++ b/src/fileFormats/vtk/part/foamVtkMeshMapsI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,24 +27,6 @@ License
 
 #include "foamVtkMeshMaps.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-inline Foam::foamVtkMeshMaps::foamVtkMeshMaps()
-:
-    cellMap_(0),
-    pointMap_(0),
-    additionalIds_(0)
-{}
-
-
-inline Foam::foamVtkMeshMaps::foamVtkMeshMaps(const label size)
-:
-    cellMap_(size),
-    pointMap_(size),
-    additionalIds_(size)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline void Foam::foamVtkMeshMaps::clear()
@@ -56,42 +38,42 @@ inline void Foam::foamVtkMeshMaps::clear()
 
 
 inline const Foam::labelList&
-Foam::foamVtkMeshMaps::cellMap() const
+Foam::foamVtkMeshMaps::cellMap() const noexcept
 {
     return cellMap_;
 }
 
 
 inline Foam::DynamicList<Foam::label>&
-Foam::foamVtkMeshMaps::cellMap()
+Foam::foamVtkMeshMaps::cellMap() noexcept
 {
     return cellMap_;
 }
 
 
 inline const Foam::labelList&
-Foam::foamVtkMeshMaps::pointMap() const
+Foam::foamVtkMeshMaps::pointMap() const noexcept
 {
     return pointMap_;
 }
 
 
 inline Foam::DynamicList<Foam::label>&
-Foam::foamVtkMeshMaps::pointMap()
+Foam::foamVtkMeshMaps::pointMap() noexcept
 {
     return pointMap_;
 }
 
 
 inline const Foam::labelList&
-Foam::foamVtkMeshMaps::additionalIds() const
+Foam::foamVtkMeshMaps::additionalIds() const noexcept
 {
     return additionalIds_;
 }
 
 
 inline Foam::DynamicList<Foam::label>&
-Foam::foamVtkMeshMaps::additionalIds()
+Foam::foamVtkMeshMaps::additionalIds() noexcept
 {
     return additionalIds_;
 }
diff --git a/src/fileFormats/vtk/part/foamVtuCells.C b/src/fileFormats/vtk/part/foamVtuCells.C
index 1bbf7c8d52..8248de41c6 100644
--- a/src/fileFormats/vtk/part/foamVtuCells.C
+++ b/src/fileFormats/vtk/part/foamVtuCells.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,30 +89,23 @@ Foam::vtk::vtuCells::vtuCells
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::vtk::vtuCells::clear()
+void Foam::vtk::vtuCells::resize_all()
 {
-    vtuSizing::clear();
-    cellTypes_.clear();
-    vertLabels_.clear();
-    vertOffset_.clear();
-    faceLabels_.clear();
-    faceOffset_.clear();
-
-    maps_.clear();
-}
-
-
-void Foam::vtk::vtuCells::repopulate(const polyMesh& mesh)
-{
-    // vtuSizing::reset() called prior to this method
-
     cellTypes_.resize(nFieldCells());
     vertLabels_.resize(sizeOf(output_, slotType::CELLS));
     vertOffset_.resize(sizeOf(output_, slotType::CELLS_OFFSETS));
     faceLabels_.resize(sizeOf(output_, slotType::FACES));
     faceOffset_.resize(sizeOf(output_, slotType::FACES_OFFSETS));
+}
+
+
+void Foam::vtk::vtuCells::populateOutput(const polyMesh& mesh)
+{
+    // Already called
+    // - vtuSizing::reset
+    // - resize_all();
 
     switch (output_)
     {
@@ -163,10 +156,103 @@ void Foam::vtk::vtuCells::repopulate(const polyMesh& mesh)
 }
 
 
+void Foam::vtk::vtuCells::populateOutput(const UList<cellShape>& shapes)
+{
+    if (output_ != contentType::LEGACY && output_ != contentType::XML)
+    {
+        WarningInFunction
+            << "Internal formats not supported for shape cells - using XML"
+            << nl << nl;
+
+        output_ = contentType::XML;
+    }
+
+    vtuSizing::resetShapes(shapes);
+
+    maps_.clear();
+    resize_all();
+    // Done in populate routine:
+    /// maps_.cellMap() = identity(vtuSizing::nCells());
+
+    switch (output_)
+    {
+        case contentType::LEGACY:
+        {
+            populateShapesLegacy
+            (
+                shapes,
+                cellTypes_,
+                vertLabels_,
+                maps_
+            );
+            break;
+        }
+
+        case contentType::XML:
+        {
+            populateShapesXml
+            (
+                shapes,
+                cellTypes_,
+                vertLabels_,
+                vertOffset_,
+                faceLabels_,
+                faceOffset_,
+                maps_
+            );
+            break;
+        }
+
+        default:
+        {
+            FatalErrorInFunction
+                << "Unhandled VTK format " << int(output_) << nl
+                << exit(FatalError);
+            break;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::vtk::vtuCells::clear()
+{
+    vtuSizing::clear();
+    cellTypes_.clear();
+    vertLabels_.clear();
+    vertOffset_.clear();
+    faceLabels_.clear();
+    faceOffset_.clear();
+
+    maps_.clear();
+}
+
+
 void Foam::vtk::vtuCells::reset(const polyMesh& mesh)
 {
     vtuSizing::reset(mesh, decomposeRequest_);
-    repopulate(mesh);
+    resize_all();
+
+    populateOutput(mesh);
+}
+
+
+void Foam::vtk::vtuCells::reset
+(
+    const polyMesh& mesh,
+    const labelUList& subsetCellsIds
+)
+{
+    vtuSizing::reset(mesh, subsetCellsIds, decomposeRequest_);
+    resize_all();
+
+    if (selectionMode() == selectionModeType::SUBSET_MESH)
+    {
+        maps_.cellMap() = subsetCellsIds;
+    }
+
+    populateOutput(mesh);
 }
 
 
@@ -184,6 +270,75 @@ void Foam::vtk::vtuCells::reset
 }
 
 
+void Foam::vtk::vtuCells::resetShapes
+(
+    const UList<cellShape>& shapes
+)
+{
+    if (output_ != contentType::LEGACY && output_ != contentType::XML)
+    {
+        WarningInFunction
+            << "VTK internal format is not supported for shape cells"
+            << " switching to xml" << nl << nl;
+
+        output_ = contentType::XML;
+    }
+
+    decomposeRequest_ = false;
+
+    vtuSizing::resetShapes(shapes);
+
+    maps_.clear();
+    resize_all();
+    maps_.cellMap() = identity(vtuSizing::nCells());
+
+    switch (output_)
+    {
+        case contentType::LEGACY:
+        {
+            populateShapesLegacy
+            (
+                shapes,
+                cellTypes_,
+                vertLabels_,
+                maps_
+            );
+            break;
+        }
+
+        case contentType::XML:
+        {
+            populateShapesXml
+            (
+                shapes,
+                cellTypes_,
+                vertLabels_,
+                vertOffset_,
+                faceLabels_,
+                faceOffset_,
+                maps_
+            );
+            break;
+        }
+
+        default:
+        {
+            FatalErrorInFunction
+                << "Unhandled VTK format " << int(output_) << nl
+                << exit(FatalError);
+            break;
+        }
+    }
+}
+
+
+void Foam::vtk::vtuCells::addPointCellLabels(const labelUList& cellIds)
+{
+    maps_.additionalIds() = cellIds;
+    setNumAddPoints(maps_.additionalIds().size());
+}
+
+
 void Foam::vtk::vtuCells::renumberCells(const labelUList& mapping)
 {
     maps_.renumberCells(mapping);
diff --git a/src/fileFormats/vtk/part/foamVtuCells.H b/src/fileFormats/vtk/part/foamVtuCells.H
index 8718caa796..197190c6d5 100644
--- a/src/fileFormats/vtk/part/foamVtuCells.H
+++ b/src/fileFormats/vtk/part/foamVtuCells.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2014-2020 OpenCFD Ltd.
+    Copyright (C) 2014-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,9 +59,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
-class polyMesh;
-
 namespace vtk
 {
 // Forward Declarations
@@ -110,9 +107,14 @@ class vtuCells
 
     // Private Member Functions
 
-        //- Create the geometry using the previously requested output and
-        //- decomposition types.
-        void repopulate(const polyMesh& mesh);
+        //- Resize storage
+        void resize_all();
+
+        //- Create the geometry entries
+        void populateOutput(const polyMesh& mesh);
+
+        //- Create the geometry entries from shapes
+        void populateOutput(const UList<cellShape>& shapes);
 
         //- No copy construct
         vtuCells(const vtuCells&) = delete;
@@ -162,10 +164,10 @@ public:
     // Access
 
         //- The output content type
-        inline enum contentType content() const;
+        inline enum contentType content() const noexcept;
 
         //- Query the polyhedral decompose requested flag.
-        inline bool decomposeRequested() const;
+        inline bool decomposeRequested() const noexcept;
 
         //- True if no cellTypes are populated.
         inline bool empty() const noexcept;
@@ -183,6 +185,14 @@ public:
         //- decomposition types.
         void reset(const polyMesh& mesh);
 
+        //- Create the geometry for a mesh subset,
+        //- using previously requested output and decomposition types.
+        void reset
+        (
+            const polyMesh& mesh,
+            const labelUList& subsetCellsIds
+        );
+
         //- Respecify requested output and decomposition type prior to
         //- creating the geometry
         void reset
@@ -192,35 +202,42 @@ public:
             const bool decompose
         );
 
+        //- Reset sizing using primitive shapes only (ADVANCED USAGE)
+        //  Effectively removes any polyhedrals!
+        void resetShapes(const UList<cellShape>& shapes);
+
         //- Renumber cell ids to account for subset meshes
         void renumberCells(const labelUList& mapping);
 
         //- Renumber point ids to account for subset meshes
         void renumberPoints(const labelUList& mapping);
 
+        //- Define which additional cell-centres are to be used (ADVANCED!)
+        void addPointCellLabels(const labelUList& cellIds);
+
 
     // Storage Access
 
         //- Values for "types" (XML) and "CELL_TYPES" (legacy)
-        inline const List<uint8_t>& cellTypes() const;
+        inline const List<uint8_t>& cellTypes() const noexcept;
 
         //- Values for "connectivity" (XML) or "CELLS" (legacy)
-        inline const labelList& vertLabels() const;
+        inline const labelList& vertLabels() const noexcept;
 
         //- Values for "offsets" (XML only)
-        inline const labelList& vertOffsets() const;
+        inline const labelList& vertOffsets() const noexcept;
 
         //- Values for "faces" (XML only)
-        inline const labelList& faceLabels() const;
+        inline const labelList& faceLabels() const noexcept;
 
         //- Values for "faceoffset" (XML only)
-        inline const labelList& faceOffsets() const;
+        inline const labelList& faceOffsets() const noexcept;
 
         //- Additional point addressing (from added point to original cell)
-        inline const labelList& addPointCellLabels() const;
+        inline const labelList& addPointCellLabels() const noexcept;
 
         //- Original cell ids for all cells (regular and decomposed).
-        inline const labelList& cellMap() const;
+        inline const labelList& cellMap() const noexcept;
 };
 
 
diff --git a/src/fileFormats/vtk/part/foamVtuCellsI.H b/src/fileFormats/vtk/part/foamVtuCellsI.H
index e6c4d91436..2734351e54 100644
--- a/src/fileFormats/vtk/part/foamVtuCellsI.H
+++ b/src/fileFormats/vtk/part/foamVtuCellsI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,13 +30,13 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline enum Foam::vtk::vtuCells::contentType
-Foam::vtk::vtuCells::content() const
+Foam::vtk::vtuCells::content() const noexcept
 {
     return output_;
 }
 
 
-inline bool Foam::vtk::vtuCells::decomposeRequested() const
+inline bool Foam::vtk::vtuCells::decomposeRequested() const noexcept
 {
     return decomposeRequest_;
 }
@@ -55,49 +55,49 @@ inline Foam::label Foam::vtk::vtuCells::size() const noexcept
 
 
 inline const Foam::List<uint8_t>&
-Foam::vtk::vtuCells::cellTypes() const
+Foam::vtk::vtuCells::cellTypes() const noexcept
 {
     return cellTypes_;
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::vertLabels() const
+Foam::vtk::vtuCells::vertLabels() const noexcept
 {
     return vertLabels_;
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::vertOffsets() const
+Foam::vtk::vtuCells::vertOffsets() const noexcept
 {
     return vertOffset_;
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::faceLabels() const
+Foam::vtk::vtuCells::faceLabels() const noexcept
 {
     return faceLabels_;
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::faceOffsets() const
+Foam::vtk::vtuCells::faceOffsets() const noexcept
 {
     return faceOffset_;
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::addPointCellLabels() const
+Foam::vtk::vtuCells::addPointCellLabels() const noexcept
 {
     return maps_.additionalIds();
 }
 
 
 inline const Foam::labelList&
-Foam::vtk::vtuCells::cellMap() const
+Foam::vtk::vtuCells::cellMap() const noexcept
 {
     return maps_.cellMap();
 }
diff --git a/src/fileFormats/vtk/part/foamVtuSizing.C b/src/fileFormats/vtk/part/foamVtuSizing.C
index 369b358b1e..c1a1688e7a 100644
--- a/src/fileFormats/vtk/part/foamVtuSizing.C
+++ b/src/fileFormats/vtk/part/foamVtuSizing.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +31,7 @@ License
 #include "cellShape.H"
 
 // Only used in this file
-#include "foamVtuSizingTemplates.C"
+#include "foamVtuSizingImpl.C"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -42,6 +42,162 @@ void Foam::vtk::vtuSizing::presizeMaps(foamVtkMeshMaps& maps) const
 }
 
 
+void Foam::vtk::vtuSizing::checkSizes
+(
+    const vtk::vtuSizing& sizing,
+
+    const label cellTypes_size,
+    const label vertLabels_size,
+    const label vertOffset_size,
+    const label faceLabels_size,
+    const label faceOffset_size,
+
+    const enum contentType output,
+    const label cellMap_size,
+    const label addPointsIds_size
+)
+{
+    label nErrors = 0;
+
+    #undef  CHECK_SIZING
+    #define CHECK_SIZING(what, sizeInput, sizeExpected)        \
+    if (sizeInput != sizeExpected)                             \
+    {                                                          \
+        if (!nErrors++)                                        \
+        {                                                      \
+            FatalErrorInFunction << "VTK sizing error" << nl;  \
+        }                                                      \
+        FatalError                                             \
+            << "    " << what << " size=" << sizeInput         \
+            << " expected " << sizeExpected << nl;             \
+    }
+
+
+    CHECK_SIZING("cellTypes", cellTypes_size, sizing.nFieldCells());
+    CHECK_SIZING("cellMap", cellMap_size, sizing.nFieldCells());
+    CHECK_SIZING("addPointsIds", addPointsIds_size, sizing.nAddPoints());
+
+    switch (output)
+    {
+        case contentType::LEGACY:
+        {
+            CHECK_SIZING("legacy", vertLabels_size, sizing.sizeLegacy());
+            break;
+        }
+
+        case contentType::XML:
+        {
+            // XML uses connectivity/offset pair.
+            CHECK_SIZING
+            (
+                "connectivity",
+                vertLabels_size,
+                sizing.sizeXml(slotType::CELLS)
+            );
+            CHECK_SIZING
+            (
+                "offsets",
+                vertOffset_size,
+                sizing.sizeXml(slotType::CELLS_OFFSETS)
+            );
+            if (sizing.nFaceLabels())
+            {
+                CHECK_SIZING
+                (
+                    "faces",
+                    faceLabels_size,
+                    sizing.sizeXml(slotType::FACES)
+                );
+
+                CHECK_SIZING
+                (
+                    "faceOffsets",
+                    faceOffset_size,
+                    sizing.sizeXml(slotType::FACES_OFFSETS)
+                );
+            }
+            break;
+        }
+
+        case contentType::INTERNAL1:
+        {
+            // VTK-internal1 connectivity/offset pair.
+            CHECK_SIZING
+            (
+                "connectivity",
+                vertLabels_size,
+                sizing.sizeInternal1(slotType::CELLS)
+            );
+            CHECK_SIZING
+            (
+                "offsets",
+                vertOffset_size,
+                sizing.sizeInternal1(slotType::CELLS_OFFSETS)
+            );
+            if (sizing.nFaceLabels())
+            {
+                CHECK_SIZING
+                (
+                    "faces",
+                    faceLabels_size,
+                    sizing.sizeInternal1(slotType::FACES)
+                );
+                CHECK_SIZING
+                (
+                    "faceOffsets",
+                    faceOffset_size,
+                    sizing.sizeInternal1(slotType::FACES_OFFSETS)
+                );
+            }
+            break;
+        }
+
+        case contentType::INTERNAL2:
+        {
+            // VTK-internal2 connectivity/offset pair.
+            CHECK_SIZING
+            (
+                "connectivity",
+                vertLabels_size,
+                sizing.sizeInternal2(slotType::CELLS)
+            );
+            CHECK_SIZING
+            (
+                "offsets",
+                vertOffset_size,
+                sizing.sizeInternal2(slotType::CELLS_OFFSETS)
+            );
+            if (sizing.nFaceLabels())
+            {
+                CHECK_SIZING
+                (
+                    "faces",
+                    faceLabels_size,
+                    sizing.sizeInternal2(slotType::FACES)
+                );
+                CHECK_SIZING
+                (
+                    "faceOffsets",
+                    faceOffset_size,
+                    sizing.sizeInternal2(slotType::FACES_OFFSETS)
+                );
+            }
+            break;
+        }
+    }
+
+    if (nErrors)
+    {
+        FatalError
+            << nl
+            << "Total of " << nErrors << " sizing errors encountered!"
+            << exit(FatalError);
+    }
+
+    #undef CHECK_SIZING
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::vtk::vtuSizing::vtuSizing() noexcept
@@ -66,6 +222,7 @@ Foam::vtk::vtuSizing::vtuSizing
 void Foam::vtk::vtuSizing::clear() noexcept
 {
     decompose_   = false;
+    selectionMode_ = FULL_MESH;
     nCells_      = 0;
     nPoints_     = 0;
     nVertLabels_ = 0;
@@ -86,20 +243,56 @@ void Foam::vtk::vtuSizing::reset
     const bool decompose
 )
 {
-    const cellModel& tet      = cellModel::ref(cellModel::TET);
-    const cellModel& pyr      = cellModel::ref(cellModel::PYR);
-    const cellModel& prism    = cellModel::ref(cellModel::PRISM);
+    reset(mesh, labelUList::null(), decompose);
+}
+
+
+void Foam::vtk::vtuSizing::reset
+(
+    const polyMesh& mesh,
+    const labelUList& subsetCellsIds,
+    const bool decompose
+)
+{
+    // References to cell shape models
+    const cellModel& tet   = cellModel::ref(cellModel::TET);
+    const cellModel& pyr   = cellModel::ref(cellModel::PYR);
+    const cellModel& prism = cellModel::ref(cellModel::PRISM);
+    const cellModel& hex   = cellModel::ref(cellModel::HEX);
     const cellModel& wedge    = cellModel::ref(cellModel::WEDGE);
     const cellModel& tetWedge = cellModel::ref(cellModel::TETWEDGE);
-    const cellModel& hex      = cellModel::ref(cellModel::HEX);
 
     const cellShapeList& shapes = mesh.cellShapes();
 
     // Unique vertex labels per polyhedral
     labelHashSet hashUniqId(2*256);
 
-    decompose_ = decompose;
-    nCells_    = mesh.nCells();
+
+    // Special treatment for mesh subsets.
+    const bool isSubsetMesh
+    (
+        notNull(subsetCellsIds)
+    );
+
+    if (isSubsetMesh)
+    {
+        decompose_  = false;  // Disallow decomposition for subset mode
+        selectionMode_ = selectionModeType::SUBSET_MESH;
+    }
+    else
+    {
+        decompose_  = decompose;  // Disallow decomposition
+        selectionMode_ = selectionModeType::FULL_MESH;
+    }
+
+    const label nInputCells =
+    (
+        isSubsetMesh
+      ? subsetCellsIds.size()
+      : shapes.size()
+    );
+
+    nCells_    = nInputCells;
     nPoints_   = mesh.nPoints();
     nAddCells_ = 0;
     nAddVerts_ = 0;
@@ -109,10 +302,10 @@ void Foam::vtk::vtuSizing::reset
     nFaceLabels_ = 0;
     nVertPoly_   = 0;
 
-    const label len = mesh.nCells();
-
-    for (label celli=0; celli < len; ++celli)
+    for (label inputi = 0; inputi < nInputCells; ++inputi)
     {
+        const label celli(isSubsetMesh ? subsetCellsIds[inputi] : inputi);
+
         const cellShape& shape = shapes[celli];
         const cellModel& model = shape.model();
 
@@ -128,15 +321,15 @@ void Foam::vtk::vtuSizing::reset
             --nCellsPoly_;
             nVertLabels_ += shape.size();
         }
-        else if (model == tetWedge && decompose)
+        else if (model == tetWedge && decompose_)
         {
             nVertLabels_ += 6;  // Treat as squeezed prism (VTK_WEDGE)
         }
-        else if (model == wedge && decompose)
+        else if (model == wedge && decompose_)
         {
             nVertLabels_ += 8;  // Treat as squeezed hex
         }
-        else if (decompose)
+        else if (decompose_)
         {
             // Polyhedral: Decompose into tets + pyramids.
             ++nAddPoints_;
@@ -199,7 +392,78 @@ void Foam::vtk::vtuSizing::reset
     }
 
     // Requested and actually required
-    decompose_ = (decompose && nCellsPoly_);
+    decompose_ = (decompose_ && nCellsPoly_);
+}
+
+
+// Synchronize changes here with the following:
+// - vtuSizing::resetShapes
+// - vtuSizing::populateArrays
+//
+void Foam::vtk::vtuSizing::resetShapes
+(
+    const UList<cellShape>& shapes
+)
+{
+    const cellModel& tet      = cellModel::ref(cellModel::TET);
+    const cellModel& pyr      = cellModel::ref(cellModel::PYR);
+    const cellModel& prism    = cellModel::ref(cellModel::PRISM);
+    const cellModel& hex      = cellModel::ref(cellModel::HEX);
+
+    decompose_ = false;  // Disallow decomposition
+    selectionMode_ = SHAPE_MESH;
+
+    const label nInputCells = shapes.size();
+
+    nCells_    = nInputCells;
+    nPoints_   = 0;
+    nAddCells_ = 0;
+    nAddVerts_ = 0;
+
+    nCellsPoly_  = 0;
+    nVertLabels_ = 0;
+    nFaceLabels_ = 0;
+    nVertPoly_   = 0;
+
+    label nIgnored = 0;
+
+    for (label inputi = 0; inputi < nInputCells; ++inputi)
+    {
+        const cellShape& shape = shapes[inputi];
+        const cellModel& model = shape.model();
+
+        if
+        (
+            model == tet
+         || model == pyr
+         || model == prism
+         || model == hex
+        )
+        {
+            nVertLabels_ += shape.size();
+
+            // Guess for number of addressed points
+            nPoints_ = max(nPoints_, max(shape));
+        }
+        else
+        {
+            --nCells_;
+            ++nIgnored;
+        }
+    }
+
+    if (nIgnored)
+    {
+        FatalErrorInFunction
+            << "Encountered " << nIgnored << " unsupported cell shapes"
+            << " ... this is likely not good" << nl
+            << exit(FatalError);
+    }
+
+    if (nCells_)
+    {
+        ++nPoints_;
+    }
 }
 
 
@@ -339,6 +603,35 @@ void Foam::vtk::vtuSizing::populateLegacy
 }
 
 
+void Foam::vtk::vtuSizing::populateShapesLegacy
+(
+    const UList<cellShape>& shapes,
+    UList<uint8_t>& cellTypes,
+    labelUList& vertLabels,
+    foamVtkMeshMaps& maps
+) const
+{
+    // Leave as zero-sized so that populateArrays doesn't fill it.
+    List<label> unused;
+
+    presizeMaps(maps);
+
+    populateArrays
+    (
+        shapes,
+        *this,
+        cellTypes,
+        vertLabels,
+        unused, // offsets
+        unused, // faces
+        unused, // facesOffsets
+        contentType::LEGACY,
+        maps.cellMap(),
+        maps.additionalIds()
+    );
+}
+
+
 void Foam::vtk::vtuSizing::populateXml
 (
     const polyMesh& mesh,
@@ -368,6 +661,38 @@ void Foam::vtk::vtuSizing::populateXml
 }
 
 
+void Foam::vtk::vtuSizing::populateShapesXml
+(
+    const UList<cellShape>& shapes,
+    UList<uint8_t>& cellTypes,
+    labelUList& connectivity,
+    labelUList& offsets,
+    labelUList& faces,
+    labelUList& facesOffsets,
+    foamVtkMeshMaps& maps
+) const
+{
+    // Leave as zero-sized so that populateArrays doesn't fill it.
+    List<label> unused;
+
+    presizeMaps(maps);
+
+    populateArrays
+    (
+        shapes,
+        *this,
+        cellTypes,
+        connectivity,
+        offsets,
+        unused, // faces
+        unused, // facesOffsets
+        contentType::XML,
+        maps.cellMap(),
+        maps.additionalIds()
+    );
+}
+
+
 #undef  definePopulateInternalMethod
 #define definePopulateInternalMethod(Type)                                   \
                                                                              \
diff --git a/src/fileFormats/vtk/part/foamVtuSizing.H b/src/fileFormats/vtk/part/foamVtuSizing.H
index 826d7f7a02..7b2c5a6081 100644
--- a/src/fileFormats/vtk/part/foamVtuSizing.H
+++ b/src/fileFormats/vtk/part/foamVtuSizing.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,7 @@ Description
 
     These are the storage entries normally associate with each output-type:
     \table
-        legacy output
+        legacy output (vtk DataFile Version 2.0)
         \c types    | vtk cell type (1-255)
         \c cells    | nLabels and unique vertex labels used by the cell, or
                     | [nLabels nFaces, nFace0Pts, id1, id2, ..., nFace1Pts, id1, id2, ...]
@@ -118,6 +118,7 @@ namespace Foam
 {
 
 // Forward Declarations
+class cellShape;
 class polyMesh;
 
 namespace vtk
@@ -134,7 +135,7 @@ public:
     // Public Data
 
     //- Types of content that the storage may represent
-    enum contentType
+    enum contentType : char
     {
         LEGACY,     //!< Legacy VTK content
         XML,        //!< XML (VTU) content
@@ -143,7 +144,7 @@ public:
     };
 
     //- The possible storage 'slots' that can be used
-    enum slotType
+    enum slotType : char
     {
         CELLS,         //!< Cell connectivity (ALL)
         CELLS_OFFSETS, //!< Cell end-offsets (XML), locations (INTERNAL1)
@@ -152,6 +153,14 @@ public:
         FACES_OFFSETS  //!< Faces end-offsets (XML) or locations (INTERNAL1)
     };
 
+    //- How the mesh cells have been selected or defined
+    enum selectionModeType : char
+    {
+        FULL_MESH,
+        SUBSET_MESH,
+        SHAPE_MESH
+    };
+
 
 private:
 
@@ -160,6 +169,9 @@ private:
         //- Polyhedral decomposition requested
         bool decompose_;
 
+        //- How the mesh cells have been selected or defined
+        selectionModeType selectionMode_;
+
         //- Number of cells in the mesh
         label nCells_;
 
@@ -199,8 +211,35 @@ private:
         //- set-size for cellMap and additionalIds
         void presizeMaps(foamVtkMeshMaps& maps) const;
 
+        //- Verify list sizes
+        static void checkSizes
+        (
+            const vtk::vtuSizing& sizing,
+
+            const label cellTypes_size,
+            const label vertLabels_size,
+            const label vertOffset_size,
+            const label faceLabels_size,
+            const label faceOffset_size,
+
+            const enum contentType output,
+
+            const label cellMap_size,
+            const label addPointsIds_size
+        );
+
+        //- Adjust cell/face offsets
+        template<class LabelType>
+        static void adjustOffsets
+        (
+            UList<LabelType>& vertOffset,
+            UList<LabelType>& faceOffset,
+            const enum contentType output,
+            const bool hasFaceStream
+        );
+
         //- Populate lists. For (legacy | xml | internal) VTK representations
-        template<class LabelType, class LabelType2>
+        template<class LabelType>
         static void populateArrays
         (
             const polyMesh& mesh,
@@ -211,10 +250,52 @@ private:
             UList<LabelType>& faceLabels,
             UList<LabelType>& faceOffset,
             const enum contentType output,
-            UList<LabelType2>& cellMap,
-            UList<LabelType2>& addPointsIds
+            labelUList& cellMap,
+            labelUList& addPointsIds
         );
 
+        //- Populate lists - Primitive mesh shapes only!!
+        template<class LabelType>
+        static void populateArrays
+        (
+            const UList<cellShape>& shapes,
+            const vtk::vtuSizing& sizing,
+            UList<uint8_t>& cellTypes,
+            UList<LabelType>& vertLabels,
+            UList<LabelType>& vertOffset,
+            UList<LabelType>& faceLabels,    // unused
+            UList<LabelType>& faceOffset,    // unused
+            const enum contentType output,
+            labelUList& cellMap,
+            labelUList& addPointsIds  // unused
+        );
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Reset list for primitive shapes only (ADVANCED USAGE)
+        void populateShapesLegacy
+        (
+            const UList<cellShape>& shapes,
+            UList<uint8_t>& cellTypes,
+            labelUList& connectivity,
+            foamVtkMeshMaps& maps
+        ) const;
+
+        //- Reset list for primitive shapes only (ADVANCED USAGE)
+        void populateShapesXml
+        (
+            const UList<cellShape>& shapes,
+            UList<uint8_t>& cellTypes,
+            labelUList& connectivity,
+            labelUList& offsets,
+            labelUList& faces,
+            labelUList& facesOffsets,
+            foamVtkMeshMaps& maps
+        ) const;
+
 
 public:
 
@@ -240,6 +321,20 @@ public:
         //  Optionally with polyhedral decomposition.
         void reset(const polyMesh& mesh, const bool decompose=false);
 
+        //- Reset sizing by analyzing a subset of the mesh.
+        //  A labelUList::null() selection corresponds to the entire mesh.
+        //  Polyhedral decomposition is disabled for subsets.
+        void reset
+        (
+            const polyMesh& mesh,
+            const labelUList& subsetCellsIds,
+            const bool decompose = false
+        );
+
+        //- Reset sizing using primitive shapes only (ADVANCED USAGE)
+        //  Effectively removes any polyhedrals!
+        void resetShapes(const UList<cellShape>& shapes);
+
         //- Reset all sizes to zero.
         void clear() noexcept;
 
@@ -247,44 +342,56 @@ public:
     // Access
 
         //- Query the decompose flag (normally off)
-        inline bool decompose() const;
+        inline bool decompose() const noexcept;
+
+        //- Query how the mesh cells have been selected or defined
+        inline selectionModeType selectionMode() const noexcept;
 
         //- Number of cells for the mesh
-        inline label nCells() const;
+        inline label nCells() const noexcept;
 
         //- Number of points for the mesh
-        inline label nPoints() const;
+        inline label nPoints() const noexcept;
 
         //- Number of vertex labels for the mesh
-        inline label nVertLabels() const;
+        inline label nVertLabels() const noexcept;
 
         //- Number of polyhedral face labels for the mesh
-        inline label nFaceLabels() const;
+        inline label nFaceLabels() const noexcept;
 
         //- Number of polyhedral cells for the mesh
-        inline label nCellsPoly() const;
+        inline label nCellsPoly() const noexcept;
 
         //- Number of vertex labels for polyhedral cells of the mesh
-        inline label nVertPoly() const;
+        inline label nVertPoly() const noexcept;
 
         //- Number of additional (decomposed) cells for the mesh
-        inline label nAddCells() const;
+        inline label nAddCells() const noexcept;
 
         //- Number of additional (decomposed) points for the mesh
-        inline label nAddPoints() const;
+        inline label nAddPoints() const noexcept;
 
         //- Number of additional (decomposed) vertices for the mesh
-        inline label nAddVerts() const;
+        inline label nAddVerts() const noexcept;
 
 
         //- Number of field cells = nCells + nAddCells
-        inline label nFieldCells() const;
+        inline label nFieldCells() const noexcept;
 
         //- Number of field points = nPoints + nAddPoints
-        inline label nFieldPoints() const;
+        inline label nFieldPoints() const noexcept;
 
 
-    // Derived sizes
+    // Edit
+
+        //- Alter number of mesh points (ADVANCED USAGE)
+        void setNumPoints(label n) noexcept { nPoints_ = n; }
+
+        //- Alter number of additional (cell-centre) points (ADVANCED USAGE)
+        void setNumAddPoints(label n) noexcept { nAddPoints_ = n; }
+
+
+    // Derived Sizes
 
         //- Return the required size for the storage slot
         label sizeOf
diff --git a/src/fileFormats/vtk/part/foamVtuSizingI.H b/src/fileFormats/vtk/part/foamVtuSizingI.H
index e3c2cc55e8..6c20571487 100644
--- a/src/fileFormats/vtk/part/foamVtuSizingI.H
+++ b/src/fileFormats/vtk/part/foamVtuSizingI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,75 +29,82 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline bool Foam::vtk::vtuSizing::decompose() const
+inline bool Foam::vtk::vtuSizing::decompose() const noexcept
 {
     return decompose_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nCells() const
+inline Foam::vtk::vtuSizing::selectionModeType
+Foam::vtk::vtuSizing::selectionMode() const noexcept
+{
+    return selectionMode_;
+}
+
+
+inline Foam::label Foam::vtk::vtuSizing::nCells() const noexcept
 {
     return nCells_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nPoints() const
+inline Foam::label Foam::vtk::vtuSizing::nPoints() const noexcept
 {
     return nPoints_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nVertLabels() const
+inline Foam::label Foam::vtk::vtuSizing::nVertLabels() const noexcept
 {
     return nVertLabels_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nFaceLabels() const
+inline Foam::label Foam::vtk::vtuSizing::nFaceLabels() const noexcept
 {
     return nFaceLabels_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nCellsPoly() const
+inline Foam::label Foam::vtk::vtuSizing::nCellsPoly() const noexcept
 {
     return nCellsPoly_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nVertPoly() const
+inline Foam::label Foam::vtk::vtuSizing::nVertPoly() const noexcept
 {
     return nVertPoly_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nAddCells() const
+inline Foam::label Foam::vtk::vtuSizing::nAddCells() const noexcept
 {
     return nAddCells_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nAddPoints() const
+inline Foam::label Foam::vtk::vtuSizing::nAddPoints() const noexcept
 {
     return nAddPoints_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nAddVerts() const
+inline Foam::label Foam::vtk::vtuSizing::nAddVerts() const noexcept
 {
     return nAddVerts_;
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nFieldCells() const
+inline Foam::label Foam::vtk::vtuSizing::nFieldCells() const noexcept
 {
-    return nCells_ + nAddCells_;
+    return (nCells_ + nAddCells_);
 }
 
 
-inline Foam::label Foam::vtk::vtuSizing::nFieldPoints() const
+inline Foam::label Foam::vtk::vtuSizing::nFieldPoints() const noexcept
 {
-    return nPoints_ + nAddPoints_;
+    return (nPoints_ + nAddPoints_);
 }
 
 
diff --git a/src/fileFormats/vtk/part/foamVtuSizingTemplates.C b/src/fileFormats/vtk/part/foamVtuSizingImpl.C
similarity index 68%
rename from src/fileFormats/vtk/part/foamVtuSizingTemplates.C
rename to src/fileFormats/vtk/part/foamVtuSizingImpl.C
index 7770f2a257..69abf42df5 100644
--- a/src/fileFormats/vtk/part/foamVtuSizingTemplates.C
+++ b/src/fileFormats/vtk/part/foamVtuSizingImpl.C
@@ -32,21 +32,170 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class LabelType, class LabelType2>
+template<class LabelType>
+void Foam::vtk::vtuSizing::adjustOffsets
+(
+    UList<LabelType>& vertOffset,
+    UList<LabelType>& faceOffset,
+    const enum contentType output,
+    const bool hasFaceStream
+)
+{
+    // ===========================================
+    // Adjust vertOffset for all cells
+    // A second pass is needed for several reasons.
+    // - Additional (decomposed) cells are placed out of sequence
+    // - INTERNAL1 connectivity has size prefixed
+    //
+    // Cell offsets:
+    // - XML format expects end-offsets,
+    // - INTERNAL1 expects begin-offsets
+    // - INTERNAL2 expects begin/end-offsets
+
+    switch (output)
+    {
+        case contentType::LEGACY: // Nothing to do
+            break;
+
+        case contentType::XML:
+        {
+            // Transform cell sizes (vertOffset) into begin offsets
+
+            // vertOffset[0] already contains its size, leave untouched
+            for (label i = 1; i < vertOffset.size(); ++i)
+            {
+                vertOffset[i] += vertOffset[i-1];
+            }
+
+            // The end face offsets, leaving -1 untouched
+            if (hasFaceStream)
+            {
+                LabelType prev(0);
+
+                for (LabelType& off : faceOffset)
+                {
+                    const LabelType sz(off);
+                    if (sz > 0)
+                    {
+                        prev += sz;
+                        off = prev;
+                    }
+                }
+            }
+            break;
+        }
+
+        case contentType::INTERNAL1:
+        {
+            // Transform cell sizes (vertOffset) into begin offsets
+            {
+                LabelType beg(0);
+
+                for (LabelType& off : vertOffset)
+                {
+                    const LabelType sz(off);
+                    off = beg;
+                    beg += 1 + sz;  // Additional 1 to skip embedded prefix
+                }
+            }
+
+            // The begin face offsets, leaving -1 untouched
+            if (hasFaceStream)
+            {
+                LabelType beg(0);
+
+                for (LabelType& off : faceOffset)
+                {
+                    const LabelType sz(off);
+                    if (sz > 0)
+                    {
+                        off = beg;
+                        beg += sz;
+                    }
+                }
+            }
+            break;
+        }
+
+        case contentType::INTERNAL2:
+        {
+            // Transform cell sizes (vertOffset) into begin/end offsets
+            // input    [n1, n2, n3, ..., 0]
+            // becomes  [0, n1, n1+n2, n1+n2+n3, ..., nTotal]
+
+            // The last entry of vertOffset was initialized as zero and
+            // never revisited, so the following loop is OK
+            {
+                LabelType total(0);
+
+                for (LabelType& off : vertOffset)
+                {
+                    const LabelType sz(off);
+                    off = total;
+                    total += sz;
+                }
+            }
+
+            // The begin face offsets, leaving -1 untouched
+            if (hasFaceStream)
+            {
+                LabelType beg(0);
+
+                for (LabelType& off : faceOffset)
+                {
+                    const LabelType sz(off);
+                    if (sz > 0)
+                    {
+                        off = beg;
+                        beg += sz;
+                    }
+                }
+            }
+            break;
+        }
+    }
+}
+
+
+template<class LabelType>
 void Foam::vtk::vtuSizing::populateArrays
 (
     const polyMesh& mesh,
     const vtk::vtuSizing& sizing,
+
     UList<uint8_t>& cellTypes,
     UList<LabelType>& vertLabels,
     UList<LabelType>& vertOffset,
     UList<LabelType>& faceLabels,
     UList<LabelType>& faceOffset,
     const enum contentType output,
-    UList<LabelType2>& cellMap,
-    UList<LabelType2>& addPointsIds
+    labelUList& cellMap,
+    labelUList& addPointsIds
 )
 {
+    if (sizing.selectionMode() == selectionModeType::SHAPE_MESH)
+    {
+        FatalErrorInFunction
+            << "Programming error ... attempting to populate a VTU mesh"
+            << " but it was originally sized using independent cell shapes"
+            << exit(FatalError);
+    }
+
+    // Verify storage sizes
+    checkSizes
+    (
+        sizing,
+
+        cellTypes.size(),
+        vertLabels.size(), vertOffset.size(),
+        faceLabels.size(), faceOffset.size(),
+
+        output,
+
+        cellMap.size(),
+        addPointsIds.size()
+    );
+
     // Characteristics
 
     // Are vertLabels prefixed with the size?
@@ -58,228 +207,6 @@ void Foam::vtk::vtuSizing::populateArrays
     ) ? 1 : 0;
 
 
-    // STAGE 1: Verify storage sizes
-
-    if (cellTypes.size() != sizing.nFieldCells())
-    {
-        FatalErrorInFunction
-            << " cellTypes size=" << cellTypes.size()
-            << " expected " << sizing.nFieldCells()
-            << exit(FatalError);
-    }
-
-    if (cellMap.size() != sizing.nFieldCells())
-    {
-        FatalErrorInFunction
-            << " cellMap size=" << cellMap.size()
-            << " expected " << sizing.nFieldCells()
-            << exit(FatalError);
-    }
-
-    if (addPointsIds.size() != sizing.nAddPoints())
-    {
-        FatalErrorInFunction
-            << " addPointsIds size=" << addPointsIds.size()
-            << " expected " << sizing.nAddPoints()
-            << exit(FatalError);
-    }
-
-    switch (output)
-    {
-        case contentType::LEGACY:
-        {
-            if (vertLabels.size() != sizing.sizeLegacy())
-            {
-                FatalErrorInFunction
-                    << " legacy size=" << vertLabels.size()
-                    << " expected " << sizing.sizeLegacy()
-                    << exit(FatalError);
-            }
-            break;
-        }
-
-        case contentType::XML:
-        {
-            // XML uses connectivity/offset pair.
-            if
-            (
-                vertLabels.size()
-             != sizing.sizeXml(slotType::CELLS)
-            )
-            {
-                FatalErrorInFunction
-                    << " connectivity size=" << vertLabels.size()
-                    << " expected "
-                    << sizing.sizeXml(slotType::CELLS)
-                    << exit(FatalError);
-            }
-
-            if
-            (
-                vertOffset.size()
-             != sizing.sizeXml(slotType::CELLS_OFFSETS)
-            )
-            {
-                FatalErrorInFunction
-                    << " offsets size=" << vertOffset.size()
-                    << " expected "
-                    << sizing.sizeXml(slotType::CELLS_OFFSETS)
-                    << exit(FatalError);
-            }
-
-            if (sizing.nFaceLabels())
-            {
-                if
-                (
-                    faceLabels.size()
-                 != sizing.sizeXml(slotType::FACES)
-                )
-                {
-                    FatalErrorInFunction
-                        << " faces size=" << faceLabels.size()
-                        << " expected "
-                        << sizing.sizeXml(slotType::FACES)
-                        << exit(FatalError);
-                }
-
-                if
-                (
-                    faceOffset.size()
-                 != sizing.sizeXml(slotType::FACES_OFFSETS)
-                )
-                {
-                    FatalErrorInFunction
-                        << " facesOffsets size=" << faceOffset.size()
-                        << " expected "
-                        << sizing.sizeXml(slotType::FACES_OFFSETS)
-                        << exit(FatalError);
-                }
-            }
-            break;
-        }
-
-        case contentType::INTERNAL1:
-        {
-            // VTK-internal1 connectivity/offset pair.
-            if
-            (
-                vertLabels.size()
-             != sizing.sizeInternal1(slotType::CELLS)
-            )
-            {
-                FatalErrorInFunction
-                    << " connectivity size=" << vertLabels.size()
-                    << " expected "
-                    << sizing.sizeInternal1(slotType::CELLS)
-                    << exit(FatalError);
-            }
-
-            if
-            (
-                vertOffset.size()
-             != sizing.sizeInternal1(slotType::CELLS_OFFSETS)
-            )
-            {
-                FatalErrorInFunction
-                    << " offsets size=" << vertOffset.size()
-                    << " expected "
-                    << sizing.sizeInternal1(slotType::CELLS_OFFSETS)
-                    << exit(FatalError);
-            }
-
-            if (sizing.nFaceLabels())
-            {
-                if
-                (
-                    faceLabels.size()
-                 != sizing.sizeInternal1(slotType::FACES)
-                )
-                {
-                    FatalErrorInFunction
-                        << " faces size=" << faceLabels.size()
-                        << " expected "
-                        << sizing.sizeInternal1(slotType::FACES)
-                        << exit(FatalError);
-                }
-
-                if
-                (
-                    faceOffset.size()
-                 != sizing.sizeInternal1(slotType::FACES_OFFSETS)
-                )
-                {
-                    FatalErrorInFunction
-                        << " facesOffsets size=" << faceOffset.size()
-                        << " expected "
-                        << sizing.sizeInternal1(slotType::FACES_OFFSETS)
-                        << exit(FatalError);
-                }
-            }
-            break;
-        }
-
-        case contentType::INTERNAL2:
-        {
-            // VTK-internal2 connectivity/offset pair.
-            if
-            (
-                vertLabels.size()
-             != sizing.sizeInternal2(slotType::CELLS)
-            )
-            {
-                FatalErrorInFunction
-                    << " connectivity size=" << vertLabels.size()
-                    << " expected "
-                    << sizing.sizeInternal2(slotType::CELLS)
-                    << exit(FatalError);
-            }
-
-            if
-            (
-                vertOffset.size()
-             != sizing.sizeInternal2(slotType::CELLS_OFFSETS)
-            )
-            {
-                FatalErrorInFunction
-                    << " offsets size=" << vertOffset.size()
-                    << " expected "
-                    << sizing.sizeInternal2(slotType::CELLS_OFFSETS)
-                    << exit(FatalError);
-            }
-
-            if (sizing.nFaceLabels())
-            {
-                if
-                (
-                    faceLabels.size()
-                 != sizing.sizeInternal2(slotType::FACES)
-                )
-                {
-                    FatalErrorInFunction
-                        << " faces size=" << faceLabels.size()
-                        << " expected "
-                        << sizing.sizeInternal2(slotType::FACES)
-                        << exit(FatalError);
-                }
-
-                if
-                (
-                    faceOffset.size()
-                 != sizing.sizeInternal2(slotType::FACES_OFFSETS)
-                )
-                {
-                    FatalErrorInFunction
-                        << " facesOffsets size=" << faceOffset.size()
-                        << " expected "
-                        << sizing.sizeInternal2(slotType::FACES_OFFSETS)
-                        << exit(FatalError);
-                }
-            }
-            break;
-        }
-    }
-
-
     // Initialization
 
     faceOffset = -1;
@@ -297,9 +224,9 @@ void Foam::vtk::vtuSizing::populateArrays
     const cellModel& tet      = cellModel::ref(cellModel::TET);
     const cellModel& pyr      = cellModel::ref(cellModel::PYR);
     const cellModel& prism    = cellModel::ref(cellModel::PRISM);
+    const cellModel& hex      = cellModel::ref(cellModel::HEX);
     const cellModel& wedge    = cellModel::ref(cellModel::WEDGE);
     const cellModel& tetWedge = cellModel::ref(cellModel::TETWEDGE);
-    const cellModel& hex      = cellModel::ref(cellModel::HEX);
 
     const cellShapeList& shapes = mesh.cellShapes();
 
@@ -316,11 +243,11 @@ void Foam::vtk::vtuSizing::populateArrays
     // Index into vertLabels for decomposed polys
     label nVertDecomp = sizing.nVertLabels() + prefix*sizing.nCells();
 
-    // Placement of decomposed cells
+    // Placement of additional decomposed cells
     label nCellDecomp = mesh.nCells();
 
     // Placement of additional point labels
-    label nPointDecomp = 0;
+    label nPointDecomp = mesh.nPoints();
 
     // Non-decomposed polyhedral are represented as a face-stream.
     // For legacy format, this stream replaces the normal connectivity
@@ -346,23 +273,48 @@ void Foam::vtk::vtuSizing::populateArrays
     // the per-cell vertLabels entries, and the faceOffset contains the *size*
     // associated with the per-cell faceLabels.
 
-    const label len = shapes.size();
 
-    for (label celli=0; celli < len; ++celli)
+    // Special treatment for mesh subsets
+    // Here the cellMap is the list of input cells!
+
+    const bool isSubsetMesh
+    (
+        sizing.selectionMode() == selectionModeType::SUBSET_MESH
+    );
+
+    const label nInputCells =
+    (
+        isSubsetMesh
+      ? cellMap.size()
+      : shapes.size()
+    );
+
+
+    for
+    (
+        label inputi = 0, cellIndex = 0; // cellIndex: the ouput location
+        inputi < nInputCells;
+        ++inputi, ++cellIndex
+    )
     {
+        const label celli(isSubsetMesh ? cellMap[inputi] : inputi);
+
         const cellShape& shape = shapes[celli];
         const cellModel& model = shape.model();
 
-        cellMap[celli] = celli;
+        if (!isSubsetMesh)
+        {
+            cellMap[cellIndex] = celli;
+        }
 
         if (model == tet)
         {
-            cellTypes[celli] = vtk::cellType::VTK_TETRA;
+            cellTypes[cellIndex] = vtk::cellType::VTK_TETRA;
             constexpr label nShapePoints = 4; // OR shape.size();
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -376,12 +328,12 @@ void Foam::vtk::vtuSizing::populateArrays
         }
         else if (model == pyr)
         {
-            cellTypes[celli] = vtk::cellType::VTK_PYRAMID;
+            cellTypes[cellIndex] = vtk::cellType::VTK_PYRAMID;
             constexpr label nShapePoints = 5; // OR shape.size();
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -395,12 +347,12 @@ void Foam::vtk::vtuSizing::populateArrays
         }
         else if (model == hex)
         {
-            cellTypes[celli] = vtk::cellType::VTK_HEXAHEDRON;
+            cellTypes[cellIndex] = vtk::cellType::VTK_HEXAHEDRON;
             constexpr label nShapePoints = 8; // OR shape.size();
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -414,12 +366,12 @@ void Foam::vtk::vtuSizing::populateArrays
         }
         else if (model == prism)
         {
-            cellTypes[celli] = vtk::cellType::VTK_WEDGE;
+            cellTypes[cellIndex] = vtk::cellType::VTK_WEDGE;
             constexpr label nShapePoints = 6; // OR shape.size();
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -437,12 +389,12 @@ void Foam::vtk::vtuSizing::populateArrays
         else if (model == tetWedge && sizing.decompose())
         {
             // Treat as squeezed prism
-            cellTypes[celli] = vtk::cellType::VTK_WEDGE;
+            cellTypes[cellIndex] = vtk::cellType::VTK_WEDGE;
             constexpr label nShapePoints = 6;
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -459,12 +411,12 @@ void Foam::vtk::vtuSizing::populateArrays
         else if (model == wedge && sizing.decompose())
         {
             // Treat as squeezed hex
-            cellTypes[celli] = vtk::cellType::VTK_HEXAHEDRON;
+            cellTypes[cellIndex] = vtk::cellType::VTK_HEXAHEDRON;
             constexpr label nShapePoints = 8;
 
             if (vertOffset.size())
             {
-                vertOffset[celli] = nShapePoints;
+                vertOffset[cellIndex] = nShapePoints;
             }
             if (prefix)
             {
@@ -492,7 +444,7 @@ void Foam::vtk::vtuSizing::populateArrays
 
             // Mapping from additional point to cell, and the new vertex from
             // the cell-centre
-            const label newVertexLabel = mesh.nPoints() + nPointDecomp;
+            const label newVertexLabel = nPointDecomp;
 
             addPointsIds[nPointDecomp++] = celli;
 
@@ -526,7 +478,7 @@ void Foam::vtk::vtuSizing::populateArrays
                     if (firstCell)
                     {
                         firstCell = false;
-                        celLoc = celli;
+                        celLoc = cellIndex;
                         vrtLoc = nVertLabels;
                         nVertLabels += prefix + nShapePoints;
                     }
@@ -578,7 +530,7 @@ void Foam::vtk::vtuSizing::populateArrays
                     if (firstCell)
                     {
                         firstCell = false;
-                        celLoc = celli;
+                        celLoc = cellIndex;
                         vrtLoc = nVertLabels;
                         nVertLabels += prefix + nShapePoints;
                     }
@@ -626,7 +578,7 @@ void Foam::vtk::vtuSizing::populateArrays
 
             // face-stream
             //   [nFaces, nFace0Pts, id1, id2, ..., nFace1Pts, id1, id2, ...]
-            cellTypes[celli] = vtk::cellType::VTK_POLYHEDRON;
+            cellTypes[cellIndex] = vtk::cellType::VTK_POLYHEDRON;
             const labelList& cFaces = mesh.cells()[celli];
 
             const label startLabel = faceIndexer;
@@ -676,9 +628,9 @@ void Foam::vtk::vtuSizing::populateArrays
             else
             {
                 // Size for face stream
-                faceOffset[celli] = (faceIndexer - startLabel);
+                faceOffset[cellIndex] = (faceIndexer - startLabel);
 
-                vertOffset[celli] = hashUniqId.size();
+                vertOffset[cellIndex] = hashUniqId.size();
                 if (prefix)
                 {
                     vertLabels[nVertLabels++] = hashUniqId.size();
@@ -703,108 +655,224 @@ void Foam::vtk::vtuSizing::populateArrays
     // - INTERNAL1 expects begin-offsets
     // - INTERNAL2 expects begin/end-offsets
 
-    switch (output)
+    adjustOffsets<LabelType>
+    (
+        vertOffset,
+        faceOffset,
+        output,
+        sizing.nFaceLabels()  // hasFaceStream
+    );
+}
+
+
+
+// Synchronize changes here with the following:
+// - vtuSizing::resetShapes
+// - vtuSizing::populateArrays
+
+template<class LabelType>
+void Foam::vtk::vtuSizing::populateArrays
+(
+    const UList<cellShape>& shapes,
+    const vtk::vtuSizing& sizing,
+
+    UList<uint8_t>& cellTypes,
+    UList<LabelType>& vertLabels,
+    UList<LabelType>& vertOffset,
+    UList<LabelType>& faceLabels,
+    UList<LabelType>& faceOffset,
+    const enum contentType output,
+    labelUList& cellMap,
+    labelUList& addPointsIds
+)
+{
+    if (sizing.selectionMode() != selectionModeType::SHAPE_MESH)
     {
-        case contentType::LEGACY: // nothing to do
-            break;
+        FatalErrorInFunction
+            << "Programming error ... attempting to populate a VTU mesh"
+            << " from cell shapes, but sizing originated from a different"
+            << " representation" << nl
+            << exit(FatalError);
+    }
 
-        case contentType::XML:
+    // Verify storage sizes
+    checkSizes
+    (
+        sizing,
+
+        cellTypes.size(),
+        vertLabels.size(), vertOffset.size(),
+        faceLabels.size(), faceOffset.size(),
+
+        output,
+
+        cellMap.size(),
+        addPointsIds.size()
+    );
+
+    // Characteristics
+
+    // Are vertLabels prefixed with the size?
+    // Also use as the size of the prefixed information
+    const int prefix =
+    (
+        output == contentType::LEGACY
+     || output == contentType::INTERNAL1
+    ) ? 1 : 0;
+
+
+    // Initialization
+
+    faceOffset = -1;
+
+    // For INTERNAL2, the vertOffset is (nFieldCells+1), which means that
+    // the last entry is never visited. Set as zero now.
+
+    if (vertOffset.size())
+    {
+        vertOffset.first() = 0;
+        vertOffset.last() = 0;
+    }
+
+
+    const cellModel& tet      = cellModel::ref(cellModel::TET);
+    const cellModel& pyr      = cellModel::ref(cellModel::PYR);
+    const cellModel& prism    = cellModel::ref(cellModel::PRISM);
+    const cellModel& hex      = cellModel::ref(cellModel::HEX);
+
+    // Index into vertLabels for normal cells
+    label nVertLabels = 0;
+
+    // ===========================================
+    // STAGE 2: Rewrite in VTK form
+    // During this stage, the vertOffset contains the *size* associated with
+    // the per-cell vertLabels entries, and the faceOffset contains the *size*
+    // associated with the per-cell faceLabels.
+
+    const label nInputCells = shapes.size();
+
+    label nIgnored = 0;
+
+    for
+    (
+        label inputi = 0, cellIndex = 0; // cellIndex: the ouput location
+        inputi < nInputCells;
+        ++inputi, ++cellIndex
+    )
+    {
+        const cellShape& shape = shapes[inputi];
+        const cellModel& model = shape.model();
+
+        if (model == tet)
         {
-            // Transform cell sizes (vertOffset) into begin offsets
+            cellTypes[cellIndex] = vtk::cellType::VTK_TETRA;
+            constexpr label nShapePoints = 4; // OR shape.size();
 
-            // vertOffset[0] already contains its size, leave untouched
-            for (label i = 1; i < vertOffset.size(); ++i)
+            if (vertOffset.size())
             {
-                vertOffset[i] += vertOffset[i-1];
+                vertOffset[cellIndex] = nShapePoints;
+            }
+            if (prefix)
+            {
+                vertLabels[nVertLabels++] = nShapePoints;
             }
 
-            // The end face offsets, leaving -1 untouched
-            if (sizing.nFaceLabels())
+            for (const label cpi : shape)
             {
-                LabelType prev(0);
-
-                for (LabelType& off : faceOffset)
-                {
-                    const LabelType sz(off);
-                    if (sz > 0)
-                    {
-                        prev += sz;
-                        off = prev;
-                    }
-                }
+                vertLabels[nVertLabels++] = cpi;
             }
-            break;
         }
-
-        case contentType::INTERNAL1:
+        else if (model == pyr)
         {
-            // Transform cell sizes (vertOffset) into begin offsets
-            {
-                LabelType beg(0);
+            cellTypes[cellIndex] = vtk::cellType::VTK_PYRAMID;
+            constexpr label nShapePoints = 5; // OR shape.size();
 
-                for (LabelType& off : vertOffset)
-                {
-                    const LabelType sz(off);
-                    off = beg;
-                    beg += 1 + sz;  // Additional 1 to skip embedded prefix
-                }
+            if (vertOffset.size())
+            {
+                vertOffset[cellIndex] = nShapePoints;
+            }
+            if (prefix)
+            {
+                vertLabels[nVertLabels++] = nShapePoints;
             }
 
-            // The begin face offsets, leaving -1 untouched
-            if (sizing.nFaceLabels())
+            for (const label cpi : shape)
             {
-                LabelType beg(0);
-
-                for (LabelType& off : faceOffset)
-                {
-                    const LabelType sz(off);
-                    if (sz > 0)
-                    {
-                        off = beg;
-                        beg += sz;
-                    }
-                }
+                vertLabels[nVertLabels++] = cpi;
             }
-            break;
         }
-
-        case contentType::INTERNAL2:
+        else if (model == hex)
         {
-            // Transform cell sizes (vertOffset) into begin/end offsets
-            // input    [n1, n2, n3, ..., 0]
-            // becomes  [0, n1, n1+n2, n1+n2+n3, ..., nTotal]
+            cellTypes[cellIndex] = vtk::cellType::VTK_HEXAHEDRON;
+            constexpr label nShapePoints = 8; // OR shape.size();
 
-            // The last entry of vertOffset was initialized as zero and
-            // never revisited, so the following loop is OK
+            if (vertOffset.size())
             {
-                LabelType total(0);
-
-                for (LabelType& off : vertOffset)
-                {
-                    const LabelType sz(off);
-                    off = total;
-                    total += sz;
-                }
+                vertOffset[cellIndex] = nShapePoints;
+            }
+            if (prefix)
+            {
+                vertLabels[nVertLabels++] = nShapePoints;
             }
 
-            // The begin face offsets, leaving -1 untouched
-            if (sizing.nFaceLabels())
+            for (const label cpi : shape)
             {
-                LabelType beg(0);
-
-                for (LabelType& off : faceOffset)
-                {
-                    const LabelType sz(off);
-                    if (sz > 0)
-                    {
-                        off = beg;
-                        beg += sz;
-                    }
-                }
+                vertLabels[nVertLabels++] = cpi;
             }
-            break;
+        }
+        else if (model == prism)
+        {
+            cellTypes[cellIndex] = vtk::cellType::VTK_WEDGE;
+            constexpr label nShapePoints = 6; // OR shape.size();
+
+            if (vertOffset.size())
+            {
+                vertOffset[cellIndex] = nShapePoints;
+            }
+            if (prefix)
+            {
+                vertLabels[nVertLabels++] = nShapePoints;
+            }
+
+            // VTK_WEDGE triangles point outwards (swap 1<->2, 4<->5)
+            vertLabels[nVertLabels++] = shape[0];
+            vertLabels[nVertLabels++] = shape[2];
+            vertLabels[nVertLabels++] = shape[1];
+            vertLabels[nVertLabels++] = shape[3];
+            vertLabels[nVertLabels++] = shape[5];
+            vertLabels[nVertLabels++] = shape[4];
+        }
+        else
+        {
+            // Silent here.
+            // - already complained (and skipped) during initial sizing
+            --cellIndex;
+            ++nIgnored;
         }
     }
+
+    // May have been done by caller,
+    // but for additional safety set an identity mapping
+    ListOps::identity(cellMap);
+
+    // ===========================================
+    // Adjust vertOffset for all cells
+    // A second pass is needed for several reasons.
+    // - Additional (decomposed) cells are placed out of sequence
+    // - INTERNAL1 connectivity has size prefixed
+    //
+    // Cell offsets:
+    // - XML format expects end-offsets,
+    // - INTERNAL1 expects begin-offsets
+    // - INTERNAL2 expects begin/end-offsets
+
+    adjustOffsets<LabelType>
+    (
+        vertOffset,
+        faceOffset,
+        output,
+        sizing.nFaceLabels()  // hasFaceStream
+    );
 }
 
 
diff --git a/src/fileFormats/vtk/write/foamVtkLineWriter.C b/src/fileFormats/vtk/write/foamVtkLineWriter.C
new file mode 100644
index 0000000000..a5f42b5ad9
--- /dev/null
+++ b/src/fileFormats/vtk/write/foamVtkLineWriter.C
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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 "foamVtkLineWriter.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::vtk::lineWriter::lineWriter
+(
+    const pointField& points,
+    const edgeList& edges,
+    const vtk::outputOptions opts
+)
+:
+    vtk::polyWriter(opts),
+
+    points_(std::cref<pointField>(points)),
+    edges_(std::cref<edgeList>(edges)),
+    instant_()
+{}
+
+
+Foam::vtk::lineWriter::lineWriter
+(
+    const pointField& points,
+    const edgeList& edges,
+    const fileName& file,
+    bool parallel
+)
+:
+    lineWriter(points, edges)
+{
+    open(file, parallel);
+}
+
+
+Foam::vtk::lineWriter::lineWriter
+(
+    const pointField& points,
+    const edgeList& edges,
+    const vtk::outputOptions opts,
+    const fileName& file,
+    bool parallel
+)
+:
+    lineWriter(points, edges, opts)
+{
+    open(file, parallel);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::vtk::lineWriter::setTime(const instant& inst)
+{
+    instant_ = inst;
+}
+
+
+bool Foam::vtk::lineWriter::beginFile(std::string title)
+{
+    if (title.size())
+    {
+        return vtk::fileWriter::beginFile(title);
+    }
+
+    if (!instant_.name().empty())
+    {
+        return vtk::fileWriter::beginFile
+        (
+            "time='" + instant_.name() + "'"
+        );
+    }
+
+    // Provide default title
+    return vtk::fileWriter::beginFile("edges");
+}
+
+
+bool Foam::vtk::lineWriter::writeGeometry()
+{
+    return writeLineGeometry(points_.get(), edges_.get());
+}
+
+
+void Foam::vtk::lineWriter::writeTimeValue()
+{
+    if (!instant_.name().empty())
+    {
+        vtk::fileWriter::writeTimeValue(instant_.value());
+    }
+}
+
+
+void Foam::vtk::lineWriter::piece
+(
+    const pointField& points,
+    const edgeList& edges
+)
+{
+    endPiece();
+
+    points_ = std::cref<pointField>(points);
+    edges_ = std::cref<edgeList>(edges);
+}
+
+
+bool Foam::vtk::lineWriter::writeProcIDs()
+{
+    return vtk::fileWriter::writeProcIDs(nLocalLines_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/fileFormats/vtk/write/foamVtkLineWriter.H b/src/fileFormats/vtk/write/foamVtkLineWriter.H
new file mode 100644
index 0000000000..36245558f2
--- /dev/null
+++ b/src/fileFormats/vtk/write/foamVtkLineWriter.H
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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::vtk::lineWriter
+
+Description
+    Write edge/points (optionally with fields)
+    as a vtp file or a legacy vtk file.
+
+    The file output states are managed by the Foam::vtk::fileWriter class.
+    FieldData (eg, TimeValue) must appear before any geometry pieces.
+
+Note
+    Parallel output is combined into a single Piece without point merging,
+    which is similar to using multi-piece data sets, but allows more
+    convenient creation as a streaming process.
+    In the future, the duplicate points at processor connections
+    may be addressed using ghost points.
+
+SourceFiles
+    foamVtkLineWriter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vtk_lineWriter_H
+#define Foam_vtk_lineWriter_H
+
+#include "foamVtkPolyWriter.H"
+#include "instant.H"
+#include <functional>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace vtk
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class vtk::lineWriter Declaration
+\*---------------------------------------------------------------------------*/
+
+class lineWriter
+:
+    public vtk::polyWriter
+{
+    // Private Member Data
+
+        //- Reference to the points
+        std::reference_wrapper<const pointField> points_;
+
+        //- Reference to the edges
+        std::reference_wrapper<const edgeList> edges_;
+
+        //- Time name/value
+        instant instant_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        lineWriter(const lineWriter&) = delete;
+
+        //- No copy assignment
+        void operator=(const lineWriter&) = delete;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components (default format INLINE_BASE64)
+        lineWriter
+        (
+            const pointField& pts,
+            const edgeList& edges,
+            const vtk::outputOptions opts = vtk::formatType::INLINE_BASE64
+        );
+
+        //- Construct from components (default format INLINE_BASE64),
+        //- and open the file for writing.
+        //  The file name is with/without an extension.
+        lineWriter
+        (
+            const pointField& pts,
+            const edgeList& edges,
+            const fileName& file,
+            bool parallel = Pstream::parRun()
+        );
+
+        //- Construct from components and open the file for writing.
+        //  The file name is with/without an extension.
+        lineWriter
+        (
+            const pointField& pts,
+            const edgeList& edges,
+            const vtk::outputOptions opts,
+            const fileName& file,
+            bool parallel = Pstream::parRun()
+        );
+
+
+    //- Destructor
+    virtual ~lineWriter() = default;
+
+
+    // Member Functions
+
+        //- Define a time name/value for the output
+        virtual void setTime(const instant& inst);
+
+        //- Write file header (non-collective)
+        //  \note Expected calling states: (OPENED).
+        virtual bool beginFile(std::string title = "");
+
+        //- Write patch topology
+        //  Also writes the file header if not previously written.
+        //  \note Must be called prior to writing CellData or PointData
+        virtual bool writeGeometry();
+
+        //- Write "TimeValue" FieldData (name as per Catalyst output)
+        //  Must be called within the FIELD_DATA state.
+        //  \note As a convenience this can also be called from
+        //      (OPENED | DECLARED) states, in which case it invokes
+        //      beginFieldData(1) internally.
+        using vtk::fileWriter::writeTimeValue;
+
+        //- Write the currently set time as "TimeValue" FieldData
+        void writeTimeValue();
+
+        //- Reset point/edge references to begin a new piece
+        void piece(const pointField& points, const edgeList& edges);
+
+
+        //- Write processor ids for each line as CellData
+        //- (no-op in serial)
+        bool writeProcIDs();
+
+        //- Write a uniform field of Cell (Line) or Point values
+        template<class Type>
+        void writeUniform(const word& fieldName, const Type& val)
+        {
+            polyWriter::writeUniformValue<Type>(nLocalLines_, fieldName, val);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace vtk
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fileFormats/vtk/write/foamVtkPolyWriter.C b/src/fileFormats/vtk/write/foamVtkPolyWriter.C
index f11894a554..28de456d9c 100644
--- a/src/fileFormats/vtk/write/foamVtkPolyWriter.C
+++ b/src/fileFormats/vtk/write/foamVtkPolyWriter.C
@@ -29,8 +29,73 @@ License
 #include "foamVtkOutput.H"
 #include "globalIndex.H"
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// The connectivity count for a list of edges
+static inline label countConnectivity(const edgeList& edges)
+{
+    return 2 * edges.size();  // An edge always has two ends
+}
+
+
+// The connectivity count for a list of faces
+static label countConnectivity(const faceList& faces)
+{
+    label nConnectivity = 0;
+
+    for (const face& f : faces)
+    {
+        nConnectivity += f.size();
+    }
+
+    return nConnectivity;
+}
+
+} // End namespace Foam
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::vtk::polyWriter::beginPiece
+(
+    const pointField& points,
+    const edgeList& edges
+)
+{
+    // Basic sizes
+    nLocalPoints_ = points.size();
+    nLocalLines_  = edges.size();
+    nLocalPolys_  = 0;
+
+    numberOfPoints_ = nLocalPoints_;
+    numberOfCells_  = nLocalLines_;
+
+    if (parallel_)
+    {
+        reduce(numberOfPoints_, sumOp<label>());
+        reduce(numberOfCells_,  sumOp<label>());
+    }
+
+
+    // Nothing else to do for legacy
+    if (legacy()) return;
+
+    if (format_)
+    {
+        format().tag
+        (
+            vtk::fileTag::PIECE,
+            vtk::fileAttr::NUMBER_OF_POINTS, numberOfPoints_,
+            vtk::fileAttr::NUMBER_OF_LINES,  numberOfCells_
+            // AND: vtk::fileAttr::NUMBER_OF_POLYS,  0
+        );
+    }
+}
+
+
 void Foam::vtk::polyWriter::beginPiece
 (
     const pointField& points,
@@ -39,16 +104,11 @@ void Foam::vtk::polyWriter::beginPiece
 {
     // Basic sizes
     nLocalPoints_ = points.size();
-    nLocalFaces_  = faces.size();
-    nLocalVerts_  = 0;
-
-    for (const face& f : faces)
-    {
-        nLocalVerts_ += f.size();
-    }
+    nLocalLines_  = 0;
+    nLocalPolys_  = faces.size();
 
     numberOfPoints_ = nLocalPoints_;
-    numberOfCells_  = nLocalFaces_;
+    numberOfCells_  = nLocalPolys_;
 
     if (parallel_)
     {
@@ -67,6 +127,7 @@ void Foam::vtk::polyWriter::beginPiece
             vtk::fileTag::PIECE,
             vtk::fileAttr::NUMBER_OF_POINTS, numberOfPoints_,
             vtk::fileAttr::NUMBER_OF_POLYS,  numberOfCells_
+            // AND: vtk::fileAttr::NUMBER_OF_LINES,  0
         );
     }
 }
@@ -79,48 +140,219 @@ void Foam::vtk::polyWriter::writePoints
 {
     this->beginPoints(numberOfPoints_);
 
-    if (parallel_ ? Pstream::master() : true)
+    if (parallel_)
     {
-        {
-            vtk::writeList(format(), points);
-        }
+        vtk::writeListParallel(format_.ref(), points);
     }
+    else
+    {
+        vtk::writeList(format(), points);
+
+    }
+
+    this->endPoints();
+}
+
+
+void Foam::vtk::polyWriter::writeLinesLegacy
+(
+    const edgeList& edges,
+    const label pointOffset
+)
+{
+    // Connectivity count without additional storage (done internally)
+    const label nLocalConns = countConnectivity(edges);
+
+    label nLines = nLocalLines_;
+    label nConns = nLocalConns;
 
     if (parallel_)
     {
-        if (Pstream::master())
+        reduce(nLines, sumOp<label>());
+        reduce(nConns, sumOp<label>());
+    }
+
+    if (nLines != numberOfCells_)
+    {
+        FatalErrorInFunction
+            << "Expecting " << numberOfCells_
+            << " edges, but found " << nLines
+            << exit(FatalError);
+    }
+
+    legacy::beginLines(os_, nLines, nConns);
+
+    labelList vertLabels(nLocalLines_ + nLocalConns);
+
+    {
+        // Legacy: size + connectivity together
+        // [nPts, id1, id2, ..., nPts, id1, id2, ...]
+
+        auto iter = vertLabels.begin();
+
+        const label off = pointOffset;
+
+        for (const edge& e : edges)
         {
-            pointField recv;
+            *iter = e.size();   // The size prefix (always 2 for an edge)
+            ++iter;
 
-            // Receive each point field and write
-            for (const int subproci : Pstream::subProcs())
-            {
-                IPstream fromProc(Pstream::commsTypes::blocking, subproci);
+            *iter = off + e.first();    // Vertex labels
+            ++iter;
 
-                {
-                    fromProc >> recv;
-
-                    vtk::writeList(format(), recv);
-                }
-            }
-        }
-        else
-        {
-            // Send
-            OPstream toProc
-            (
-                Pstream::commsTypes::blocking,
-                Pstream::masterNo()
-            );
-
-            {
-                toProc << points;
-            }
+            *iter = off + e.second();
+            ++iter;
         }
     }
 
 
-    this->endPoints();
+    if (parallel_)
+    {
+        vtk::writeListParallel(format_.ref(), vertLabels);
+    }
+    else
+    {
+        vtk::writeList(format(), vertLabels);
+    }
+
+    if (format_)
+    {
+        format().flush();
+    }
+}
+
+
+void Foam::vtk::polyWriter::writeLines
+(
+    const edgeList& edges,
+    const label pointOffset
+)
+{
+    // Connectivity count without additional storage (done internally)
+    const label nLocalConns = countConnectivity(edges);
+
+    if (format_)
+    {
+        format().tag(vtk::fileTag::LINES);
+    }
+
+    //
+    // 'connectivity'
+    //
+    {
+        labelList vertLabels(nLocalConns);
+
+        label nConns = nLocalConns;
+
+        if (parallel_)
+        {
+            reduce(nConns, sumOp<label>());
+        }
+
+        if (format_)
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nConns);
+
+            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
+            format().writeSize(payLoad * sizeof(label));
+        }
+
+        {
+            // XML: connectivity only
+            // [id1, id2, ..., id1, id2, ...]
+
+            auto iter = vertLabels.begin();
+
+            const label off = pointOffset;
+
+            for (const edge& e : edges)
+            {
+                // Edge vertex labels
+                *iter = off + e.first();
+                ++iter;
+
+                *iter = off + e.second();
+                ++iter;
+            }
+        }
+
+
+        if (parallel_)
+        {
+            vtk::writeListParallel(format_.ref(), vertLabels);
+        }
+        else
+        {
+            vtk::writeList(format(), vertLabels);
+        }
+
+        if (format_)
+        {
+            format().flush();
+            format().endDataArray();
+        }
+    }
+
+
+    //
+    // 'offsets'  (connectivity offsets)
+    //
+    {
+        labelList vertOffsets(nLocalLines_);
+        label nOffs = vertOffsets.size();
+
+        if (parallel_)
+        {
+            reduce(nOffs, sumOp<label>());
+        }
+
+        if (format_)
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nOffs);
+
+            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
+            format().writeSize(payLoad);
+        }
+
+
+        // processor-local connectivity offsets
+        label off =
+        (
+            parallel_ ? globalIndex(nLocalConns).localStart() : 0
+        );
+
+
+        auto iter = vertOffsets.begin();
+
+        for (const edge& e : edges)
+        {
+            off += e.size();   // End offset
+            *iter = off;
+            ++iter;
+        }
+
+
+        if (parallel_)
+        {
+            vtk::writeListParallel(format_.ref(), vertOffsets);
+        }
+        else
+        {
+            vtk::writeList(format_.ref(), vertOffsets);
+        }
+
+
+        if (format_)
+        {
+            format().flush();
+            format().endDataArray();
+        }
+    }
+
+    if (format_)
+    {
+        format().endTag(vtk::fileTag::LINES);
+    }
 }
 
 
@@ -131,27 +363,28 @@ void Foam::vtk::polyWriter::writePolysLegacy
 )
 {
     // Connectivity count without additional storage (done internally)
+    const label nLocalConns = countConnectivity(faces);
 
-    label nFaces = nLocalFaces_;
-    label nVerts = nLocalVerts_;
+    label nPolys = nLocalPolys_;
+    label nConns = nLocalConns;
 
     if (parallel_)
     {
-        reduce(nFaces, sumOp<label>());
-        reduce(nVerts, sumOp<label>());
+        reduce(nPolys, sumOp<label>());
+        reduce(nConns, sumOp<label>());
     }
 
-    if (nFaces != numberOfCells_)
+    if (nPolys != numberOfCells_)
     {
         FatalErrorInFunction
             << "Expecting " << numberOfCells_
-            << " faces, but found " << nFaces
+            << " faces, but found " << nPolys
             << exit(FatalError);
     }
 
-    legacy::beginPolys(os_, nFaces, nVerts);
+    legacy::beginPolys(os_, nPolys, nConns);
 
-    labelList vertLabels(nLocalFaces_ + nLocalVerts_);
+    labelList vertLabels(nLocalPolys_ + nLocalConns);
 
     {
         // Legacy: size + connectivity together
@@ -159,21 +392,18 @@ void Foam::vtk::polyWriter::writePolysLegacy
 
         auto iter = vertLabels.begin();
 
-        label off = pointOffset;
+        const label off = pointOffset;
 
+        for (const face& f : faces)
         {
-            for (const face& f : faces)
-            {
-                *iter = f.size();       // The size prefix
-                ++iter;
+            *iter = f.size();       // The size prefix
+            ++iter;
 
-                for (const label pfi : f)
-                {
-                    *iter = pfi + off;  // Face vertex label
-                    ++iter;
-                }
+            for (const label id : f)
+            {
+                *iter = id + off;   // Vertex label
+                ++iter;
             }
-            // off += points.size();
         }
     }
 
@@ -200,6 +430,9 @@ void Foam::vtk::polyWriter::writePolys
     const label pointOffset
 )
 {
+    // Connectivity count without additional storage (done internally)
+    const label nLocalConns = countConnectivity(faces);
+
     if (format_)
     {
         format().tag(vtk::fileTag::POLYS);
@@ -209,18 +442,18 @@ void Foam::vtk::polyWriter::writePolys
     // 'connectivity'
     //
     {
-        labelList vertLabels(nLocalVerts_);
+        labelList vertLabels(nLocalConns);
 
-        label nVerts = nLocalVerts_;
+        label nConns = nLocalConns;
 
         if (parallel_)
         {
-            reduce(nVerts, sumOp<label>());
+            reduce(nConns, sumOp<label>());
         }
 
         if (format_)
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
+            const uint64_t payLoad = vtk::sizeofData<label>(nConns);
 
             format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
             format().writeSize(payLoad * sizeof(label));
@@ -234,16 +467,13 @@ void Foam::vtk::polyWriter::writePolys
 
             label off = pointOffset;
 
+            for (const face& f : faces)
             {
-                for (const face& f : faces)
+                for (const label id : f)
                 {
-                    for (const label pfi : f)
-                    {
-                        *iter = pfi + off;  // Face vertex label
-                        ++iter;
-                    }
+                    *iter = id + off;  // Face vertex label
+                    ++iter;
                 }
-                // off += points.size();
             }
         }
 
@@ -269,7 +499,7 @@ void Foam::vtk::polyWriter::writePolys
     // 'offsets'  (connectivity offsets)
     //
     {
-        labelList vertOffsets(nLocalFaces_);
+        labelList vertOffsets(nLocalPolys_);
         label nOffs = vertOffsets.size();
 
         if (parallel_)
@@ -289,19 +519,17 @@ void Foam::vtk::polyWriter::writePolys
         // processor-local connectivity offsets
         label off =
         (
-            parallel_ ? globalIndex(nLocalVerts_).localStart() : 0
+            parallel_ ? globalIndex(nLocalConns).localStart() : 0
         );
 
 
         auto iter = vertOffsets.begin();
 
+        for (const face& f : faces)
         {
-            for (const face& f : faces)
-            {
-                off += f.size();   // End offset
-                *iter = off;
-                ++iter;
-            }
+            off += f.size();  // End offset
+            *iter = off;
+            ++iter;
         }
 
 
@@ -340,8 +568,8 @@ Foam::vtk::polyWriter::polyWriter
     numberOfPoints_(0),
     numberOfCells_(0),
     nLocalPoints_(0),
-    nLocalFaces_(0),
-    nLocalVerts_(0)
+    nLocalLines_(0),
+    nLocalPolys_(0)
 {
     // We do not currently support append mode
     opts_.append(false);
@@ -387,6 +615,36 @@ bool Foam::vtk::polyWriter::writeGeometry()
 }
 
 
+bool Foam::vtk::polyWriter::writeLineGeometry
+(
+    const pointField& points,
+    const edgeList& edges
+)
+{
+    enter_Piece();
+
+    beginPiece(points, edges);
+
+    writePoints(points);
+
+    const label pointOffset =
+    (
+        parallel_ ? globalIndex(nLocalPoints_).localStart() : 0
+    );
+
+    if (legacy())
+    {
+        writeLinesLegacy(edges, pointOffset);
+    }
+    else
+    {
+        writeLines(edges, pointOffset);
+    }
+
+    return true;
+}
+
+
 bool Foam::vtk::polyWriter::writePolyGeometry
 (
     const pointField& points,
diff --git a/src/fileFormats/vtk/write/foamVtkPolyWriter.H b/src/fileFormats/vtk/write/foamVtkPolyWriter.H
index 8dadb17f75..4c10b0f2bf 100644
--- a/src/fileFormats/vtk/write/foamVtkPolyWriter.H
+++ b/src/fileFormats/vtk/write/foamVtkPolyWriter.H
@@ -51,6 +51,7 @@ SourceFiles
 
 #include "foamVtkFileWriter.H"
 #include "pointField.H"
+#include "edgeList.H"
 #include "faceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -75,30 +76,53 @@ protected:
         //- The number of field points for the current Piece
         label numberOfPoints_;
 
-        //- The number of field cells (faces) for the current Piece
+        //- The number of field cells (edges or faces) for the current Piece
         label numberOfCells_;
 
         //- Local number of points
         label nLocalPoints_;
 
-        //- Local number of faces
-        label nLocalFaces_;
+        //- Local number of lines (edges)
+        label nLocalLines_;
 
-        //- Local face vertices (connectivity) count. Sum of face sizes.
-        label nLocalVerts_;
+        //- Local number of polys (faces)
+        label nLocalPolys_;
 
 
+    // Protected Member Functions
+
+        //- Write a uniform field of Cell (Poly or Line) or Point values
+        template<class Type>
+        void writeUniformValue
+        (
+            const label nCellValues,  // Could be Poly or Line!
+            const word& fieldName,
+            const Type& val
+        );
+
 private:
 
     // Private Member Functions
 
-        //- Determine sizes (nLocalPoints_, nLocalFaces_, nLocalVerts_),
+        //- Determine sizes (nLocalPoints_, nLocalLines_),
+        //- and begin piece
+        void beginPiece(const pointField& points, const edgeList& edges);
+
+        //- Determine sizes (nLocalPoints_, nLocalPolys_),
         //- and begin piece
         void beginPiece(const pointField& points, const faceList& faces);
 
         //- Write points
         void writePoints(const pointField& points);
 
+        //- Write lines, legacy format
+        //  \param pointOffset processor-local point offset
+        void writeLinesLegacy(const edgeList& edges, const label pointOffset);
+
+        //- Write lines
+        //  \param pointOffset processor-local point offset
+        void writeLines(const edgeList& edges, const label pointOffset);
+
         //- Write faces, legacy format
         //  \param pointOffset processor-local point offset
         void writePolysLegacy(const faceList& faces, const label pointOffset);
@@ -163,6 +187,16 @@ public:
         //  their own geomety, but use writePolyGeometry() directly
         virtual bool writeGeometry();
 
+        //- Low-level write edge/point topology.
+        //- Normally used by writeGeometry() in a derived class
+        //  Also writes the file header if not previously written.
+        //  \note Must be called prior to writing CellData or PointData
+        bool writeLineGeometry
+        (
+            const pointField& points,
+            const edgeList& edges
+        );
+
         //- Low-level write face/point topology.
         //- Normally used by writeGeometry() in a derived class
         //  Also writes the file header if not previously written.
@@ -194,11 +228,7 @@ public:
 
     // Write
 
-        //- Write a uniform field of Cell (Face) or Point values
-        template<class Type>
-        void writeUniform(const word& fieldName, const Type& val);
-
-        //- Write a list of Cell (Face) or Point values
+        //- Write a list of Cell (Poly or Line) or Point values
         template<class Type>
         void write(const word& fieldName, const UList<Type>& field);
 };
diff --git a/src/fileFormats/vtk/write/foamVtkPolyWriterTemplates.C b/src/fileFormats/vtk/write/foamVtkPolyWriterTemplates.C
index b0db6d1cfe..319836b463 100644
--- a/src/fileFormats/vtk/write/foamVtkPolyWriterTemplates.C
+++ b/src/fileFormats/vtk/write/foamVtkPolyWriterTemplates.C
@@ -25,24 +25,27 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-void Foam::vtk::polyWriter::writeUniform
+void Foam::vtk::polyWriter::writeUniformValue
 (
+    const label nCellValues,
     const word& fieldName,
     const Type& val
 )
 {
+    label nValues(0);
+
     if (isState(outputState::CELL_DATA))
     {
         ++nCellData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfCells_);
+        nValues = nCellValues;
     }
     else if (isState(outputState::POINT_DATA))
     {
         ++nPointData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfPoints_);
+        nValues = nLocalPoints_;
     }
     else
     {
@@ -54,10 +57,16 @@ void Foam::vtk::polyWriter::writeUniform
         )
             << " for uniform field " << fieldName << nl << endl
             << exit(FatalError);
+
+        return;
     }
+
+    vtk::fileWriter::writeUniform<Type>(fieldName, val, nValues);
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 template<class Type>
 void Foam::vtk::polyWriter::write
 (
@@ -66,18 +75,18 @@ void Foam::vtk::polyWriter::write
 )
 {
     // Could check sizes:
-    //     nValues == nLocalFaces (CELL_DATA)
-    //     nValues == nLocalPoints (POINT_DATA)
+    // CELL_DATA:   nValues == (nLocalPolys | nLocalLines)
+    // POINT_DATA:  nValues == nLocalPoints
+
+    // const label nValues = field.size();
 
     if (isState(outputState::CELL_DATA))
     {
         ++nCellData_;
-        vtk::fileWriter::writeBasicField<Type>(fieldName, field);
     }
     else if (isState(outputState::POINT_DATA))
     {
         ++nPointData_;
-        vtk::fileWriter::writeBasicField<Type>(fieldName, field);
     }
     else
     {
@@ -89,7 +98,10 @@ void Foam::vtk::polyWriter::write
         )
             << " for field " << fieldName << nl << endl
             << exit(FatalError);
+        return;
     }
+
+    vtk::fileWriter::writeBasicField<Type>(fieldName, field);
 }
 
 
diff --git a/src/fileFormats/vtk/write/foamVtkSurfaceWriter.C b/src/fileFormats/vtk/write/foamVtkSurfaceWriter.C
index 9ee9709d86..7b07cace36 100644
--- a/src/fileFormats/vtk/write/foamVtkSurfaceWriter.C
+++ b/src/fileFormats/vtk/write/foamVtkSurfaceWriter.C
@@ -88,7 +88,7 @@ bool Foam::vtk::surfaceWriter::beginFile(std::string title)
         return vtk::fileWriter::beginFile(title);
     }
 
-    if (instant_.name().size())
+    if (!instant_.name().empty())
     {
         return vtk::fileWriter::beginFile
         (
@@ -129,4 +129,10 @@ void Foam::vtk::surfaceWriter::piece
 }
 
 
+bool Foam::vtk::surfaceWriter::writeProcIDs()
+{
+    return vtk::fileWriter::writeProcIDs(nLocalPolys_);
+}
+
+
 // ************************************************************************* //
diff --git a/src/fileFormats/vtk/write/foamVtkSurfaceWriter.H b/src/fileFormats/vtk/write/foamVtkSurfaceWriter.H
index 830ee6132d..0a9d2ccdf3 100644
--- a/src/fileFormats/vtk/write/foamVtkSurfaceWriter.H
+++ b/src/fileFormats/vtk/write/foamVtkSurfaceWriter.H
@@ -42,7 +42,6 @@ Note
 
 SourceFiles
     foamVtkSurfaceWriter.C
-    foamVtkSurfaceWriterTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -152,8 +151,20 @@ public:
         //- Write the currently set time as "TimeValue" FieldData
         void writeTimeValue();
 
-        //- Reset point, face references to begin a new piece
+        //- Reset point/face references to begin a new piece
         void piece(const pointField& points, const faceList& faces);
+
+
+        //- Write processor ids for each poly as CellData
+        //- (no-op in serial)
+        bool writeProcIDs();
+
+        //- Write a uniform field of Cell (Poly) or Point values
+        template<class Type>
+        void writeUniform(const word& fieldName, const Type& val)
+        {
+            polyWriter::writeUniformValue<Type>(nLocalPolys_, fieldName, val);
+        }
 };
 
 
diff --git a/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.C b/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.C
index 7f1951a971..782ce50b0e 100644
--- a/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.C
+++ b/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,50 +28,9 @@ License
 
 #include "VTKedgeFormat.H"
 #include "Fstream.H"
-#include "clock.H"
-#include "vtkUnstructuredReader.H"
 #include "Time.H"
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-void Foam::fileFormats::VTKedgeFormat::writeHeader
-(
-    Ostream& os,
-    const pointField& pointLst
-)
-{
-    // Write header
-    os  << "# vtk DataFile Version 2.0" << nl
-        << "featureEdgeMesh written " << clock::dateTime().c_str() << nl
-        << "ASCII" << nl
-        << nl
-        << "DATASET POLYDATA" << nl;
-
-    // Write vertex coords
-    os  << "POINTS " << pointLst.size() << " double" << nl;
-    for (const point& pt : pointLst)
-    {
-        os  << float(pt.x()) << ' '
-            << float(pt.y()) << ' '
-            << float(pt.z()) << nl;
-    }
-}
-
-
-void Foam::fileFormats::VTKedgeFormat::writeEdges
-(
-    Ostream& os,
-    const UList<edge>& edgeLst
-)
-{
-    os  << "LINES " << edgeLst.size() << ' ' << 3*edgeLst.size() << nl;
-
-    for (const edge& e : edgeLst)
-    {
-        os  << "2 " << e[0] << ' ' << e[1] << nl;
-    }
-}
-
+#include "foamVtkLineWriter.H"
+#include "vtkUnstructuredReader.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -122,19 +81,21 @@ bool Foam::fileFormats::VTKedgeFormat::read
     storedPoints().transfer(reader.points());
 
     label nEdges = 0;
-    forAll(reader.lines(), lineI)
+    for (const auto& lineVerts : reader.lines())
     {
-        nEdges += reader.lines()[lineI].size()-1;
+        if (lineVerts.size() > 1)
+        {
+            nEdges += (lineVerts.size()-1);
+        }
     }
-    storedEdges().setSize(nEdges);
+    storedEdges().resize(nEdges);
 
     nEdges = 0;
-    forAll(reader.lines(), lineI)
+    for (const auto& lineVerts : reader.lines())
     {
-        const labelList& verts = reader.lines()[lineI];
-        for (label i = 1; i < verts.size(); i++)
+        for (label i = 1; i < lineVerts.size(); ++i)
         {
-            storedEdges()[nEdges++] = edge(verts[i-1], verts[i]);
+            storedEdges()[nEdges++] = edge(lineVerts[i-1], lineVerts[i]);
         }
     }
 
@@ -148,16 +109,20 @@ void Foam::fileFormats::VTKedgeFormat::write
     const edgeMesh& eMesh
 )
 {
-    OFstream os(filename);
-    if (!os.good())
-    {
-        FatalErrorInFunction
-            << "Cannot open file for writing " << filename
-            << exit(FatalError);
-    }
+    // NB: restrict output to legacy ascii so that we are still able
+    // to read it with vtkUnstructuredReader
 
-    writeHeader(os, eMesh.points());
-    writeEdges(os, eMesh.edges());
+    vtk::lineWriter writer
+    (
+        eMesh.points(),
+        eMesh.edges(),
+        vtk::formatType::LEGACY_ASCII,
+        filename,
+        false  // non-parallel write (edgeMesh already serialized)
+    );
+
+    writer.beginFile("OpenFOAM edgeMesh");
+    writer.writeGeometry();
 }
 
 
diff --git a/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.H b/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.H
index ac88e0bc8d..084828d639 100644
--- a/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.H
+++ b/src/meshTools/edgeMesh/edgeMeshFormats/vtk/VTKedgeFormat.H
@@ -54,37 +54,12 @@ class VTKedgeFormat
 :
     public edgeMesh
 {
-    // Private Member Functions
-
-        //- No copy construct
-        VTKedgeFormat(const VTKedgeFormat&) = delete;
-
-        //- No copy assignment
-        void operator=(const VTKedgeFormat&) = delete;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Write header information with points
-        static void writeHeader
-        (
-            Ostream&,
-            const pointField&
-        );
-
-        //- Write edges
-        static void writeEdges(Ostream&, const UList<edge>&);
-
-
 public:
 
-
     // Constructors
 
-        //- Construct from file name
-        VTKedgeFormat(const fileName&);
+        //- Read construct from file name
+        explicit VTKedgeFormat(const fileName& filename);
 
 
     // Selectors
@@ -103,12 +78,12 @@ public:
     // Member Functions
 
         //- Write surface mesh components by proxy
-        static void write(const fileName&, const edgeMesh&);
+        static void write(const fileName& name, const edgeMesh&);
 
         //- Read from file
-        virtual bool read(const fileName&);
+        virtual bool read(const fileName& name);
 
-        //- Write object file
+        //- Write object to file
         virtual void write(const fileName& name) const
         {
             write(name, *this);
diff --git a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.C b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.C
index be2bf99fa4..a5ee5abd52 100644
--- a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.C
+++ b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.C
@@ -45,16 +45,24 @@ void Foam::vtk::internalMeshWriter::beginPiece()
 
     if (parallel_)
     {
+        if (debug > 1)
+        {
+            PoutInFunction
+                << ": nPoints=" << numberOfPoints_
+                << " nCells=" << numberOfCells_ << nl;
+        }
+
         reduce(numberOfPoints_, sumOp<label>());
         reduce(numberOfCells_,  sumOp<label>());
     }
 
+    DebugInFunction
+        << "nPoints=" << numberOfPoints_
+        << " nCells=" << numberOfCells_ << nl;
+
     // Nothing else to do for legacy
     if (legacy()) return;
 
-    DebugInFunction
-        << "nPoints=" << numberOfPoints_ << " nCells=" << numberOfCells_ << nl;
-
     if (format_)
     {
         format()
@@ -178,7 +186,10 @@ void Foam::vtk::internalMeshWriter::writeCellsLegacy(const label pointOffset)
 }
 
 
-void Foam::vtk::internalMeshWriter::writeCellsConnectivity(const label pointOffset)
+void Foam::vtk::internalMeshWriter::writeCellsConnectivity
+(
+    const label pointOffset
+)
 {
     //
     // 'connectivity'
@@ -316,7 +327,10 @@ void Foam::vtk::internalMeshWriter::writeCellsConnectivity(const label pointOffs
 }
 
 
-void Foam::vtk::internalMeshWriter::writeCellsFaces(const label pointOffset)
+void Foam::vtk::internalMeshWriter::writeCellsFaces
+(
+    const label pointOffset
+)
 {
     label nFaceLabels = vtuCells_.faceLabels().size();
 
@@ -628,39 +642,7 @@ bool Foam::vtk::internalMeshWriter::writeProcIDs()
         return false;
     }
 
-    if (isState(outputState::CELL_DATA))
-    {
-        ++nCellData_;
-    }
-    else
-    {
-        reportBadState(FatalErrorInFunction, outputState::CELL_DATA)
-            << " for procID field" << nl << endl
-            << exit(FatalError);
-    }
-
-    const globalIndex procMaps(vtuCells_.nFieldCells());
-
-    this->beginDataArray<label>("procID", procMaps.size());
-
-    bool good = false;
-
-    if (Pstream::master())
-    {
-        // Per-processor ids
-        for (const int proci : Pstream::allProcs())
-        {
-            vtk::write(format(), label(proci), procMaps.localSize(proci));
-        }
-
-        good = true;
-    }
-
-    this->endDataArray();
-
-
-    // MPI barrier
-    return returnReduce(good, orOp<bool>());
+    return vtk::fileWriter::writeProcIDs(vtuCells_.nFieldCells());
 }
 
 
diff --git a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.H b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.H
index 17c6989303..00da9739a7 100644
--- a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.H
+++ b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriter.H
@@ -225,6 +225,10 @@ public:
         //- Write primitive field of CellData
         template<class Type>
         void writeCellData(const word& fieldName, const UList<Type>& field);
+
+        //- Write primitive field of PointData
+        template<class Type>
+        void writePointData(const word& fieldName, const UList<Type>& field);
 };
 
 
diff --git a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriterTemplates.C b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriterTemplates.C
index 866117286a..ecf419c903 100644
--- a/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriterTemplates.C
+++ b/src/meshTools/output/vtk/mesh/foamVtkInternalMeshWriterTemplates.C
@@ -37,15 +37,17 @@ void Foam::vtk::internalMeshWriter::writeUniform
     const Type& val
 )
 {
+    label nValues(0);
+
     if (isState(outputState::CELL_DATA))
     {
         ++nCellData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfCells_);
+        nValues = vtuCells_.nFieldCells();
     }
     else if (isState(outputState::POINT_DATA))
     {
         ++nPointData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfPoints_);
+        nValues = vtuCells_.nFieldPoints();
     }
     else
     {
@@ -54,9 +56,14 @@ void Foam::vtk::internalMeshWriter::writeUniform
             FatalErrorInFunction,
             outputState::CELL_DATA,
             outputState::POINT_DATA
-        )   << " for field " << fieldName << nl << endl
+        )
+            << " for uniform field " << fieldName << nl << endl
             << exit(FatalError);
+
+        return;
     }
+
+    vtk::fileWriter::writeUniform<Type>(fieldName, val, nValues);
 }
 
 
@@ -95,4 +102,37 @@ void Foam::vtk::internalMeshWriter::writeCellData
 }
 
 
+template<class Type>
+void Foam::vtk::internalMeshWriter::writePointData
+(
+    const word& fieldName,
+    const UList<Type>& field
+)
+{
+    if (isState(outputState::POINT_DATA))
+    {
+        ++nPointData_;
+    }
+    else
+    {
+        reportBadState(FatalErrorInFunction, outputState::POINT_DATA)
+            << " for field " << fieldName << nl << endl
+            << exit(FatalError);
+    }
+
+    this->beginDataArray<Type>(fieldName, numberOfPoints_);
+
+    if (parallel_)
+    {
+        vtk::writeListParallel(format_.ref(), field);
+    }
+    else
+    {
+        vtk::writeList(format(), field);
+    }
+
+    this->endDataArray();
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.C b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.C
index acec853d6b..63f85246a6 100644
--- a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.C
+++ b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.C
@@ -38,14 +38,15 @@ void Foam::vtk::patchMeshWriter::beginPiece()
     // Basic sizes
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    nLocalPoints_ = nLocalFaces_ = nLocalVerts_ = 0;
+    nLocalPoints_ = nLocalPolys_ = 0;
+    nLocalVerts_ = 0;
 
     for (const label patchId : patchIDs_)
     {
         const polyPatch& pp = patches[patchId];
 
         nLocalPoints_ += pp.nPoints();
-        nLocalFaces_  += pp.size();
+        nLocalPolys_  += pp.size();
 
         for (const face& f : pp)
         {
@@ -54,7 +55,7 @@ void Foam::vtk::patchMeshWriter::beginPiece()
     }
 
     numberOfPoints_ = nLocalPoints_;
-    numberOfCells_ = nLocalFaces_;
+    numberOfCells_ = nLocalPolys_;
 
     if (parallel_)
     {
@@ -148,26 +149,26 @@ void Foam::vtk::patchMeshWriter::writePolysLegacy(const label pointOffset)
 
     // Connectivity count without additional storage (done internally)
 
-    label nFaces = nLocalFaces_;
+    label nPolys = nLocalPolys_;
     label nVerts = nLocalVerts_;
 
     if (parallel_)
     {
-        reduce(nFaces, sumOp<label>());
+        reduce(nPolys, sumOp<label>());
         reduce(nVerts, sumOp<label>());
     }
 
-    if (nFaces != numberOfCells_)
+    if (nPolys != numberOfCells_)
     {
         FatalErrorInFunction
             << "Expecting " << numberOfCells_
-            << " faces, but found " << nFaces
+            << " faces, but found " << nPolys
             << exit(FatalError);
     }
 
-    legacy::beginPolys(os_, nFaces, nVerts);
+    legacy::beginPolys(os_, nPolys, nVerts);
 
-    labelList vertLabels(nLocalFaces_ + nLocalVerts_);
+    labelList vertLabels(nLocalPolys_ + nLocalVerts_);
 
     {
         // Legacy: size + connectivity together
@@ -186,9 +187,9 @@ void Foam::vtk::patchMeshWriter::writePolysLegacy(const label pointOffset)
                 *iter = f.size();       // The size prefix
                 ++iter;
 
-                for (const label pfi : f)
+                for (const label id : f)
                 {
-                    *iter = pfi + off;  // Face vertex label
+                    *iter = id + off;   // Face vertex label
                     ++iter;
                 }
             }
@@ -258,9 +259,9 @@ void Foam::vtk::patchMeshWriter::writePolys(const label pointOffset)
 
                 for (const face& f : pp.localFaces())
                 {
-                    for (const label pfi : f)
+                    for (const label id : f)
                     {
-                        *iter = pfi + off;  // Face vertex label
+                        *iter = id + off;   // Face vertex label
                         ++iter;
                     }
                 }
@@ -290,7 +291,7 @@ void Foam::vtk::patchMeshWriter::writePolys(const label pointOffset)
     // 'offsets'  (connectivity offsets)
     //
     {
-        labelList vertOffsets(nLocalFaces_);
+        labelList vertOffsets(nLocalPolys_);
         label nOffs = vertOffsets.size();
 
         if (parallel_)
@@ -367,7 +368,7 @@ Foam::vtk::patchMeshWriter::patchMeshWriter
     numberOfPoints_(0),
     numberOfCells_(0),
     nLocalPoints_(0),
-    nLocalFaces_(0),
+    nLocalPolys_(0),
     nLocalVerts_(0),
 
     mesh_(mesh),
@@ -512,15 +513,15 @@ void Foam::vtk::patchMeshWriter::writePatchIDs()
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    label nFaces = nLocalFaces_;
+    label nPolys = nLocalPolys_;
 
     if (parallel_)
     {
-        reduce(nFaces, sumOp<label>());
+        reduce(nPolys, sumOp<label>());
     }
 
 
-    this->beginDataArray<label>("patchID", nFaces);
+    this->beginDataArray<label>("patchID", nPolys);
 
     if (parallel_ ? Pstream::master() : true)
     {
@@ -584,66 +585,7 @@ void Foam::vtk::patchMeshWriter::writePatchIDs()
 
 bool Foam::vtk::patchMeshWriter::writeProcIDs()
 {
-    // This is different than for internalWriter.
-    // Here we allow procIDs whenever running in parallel, even if the
-    // output is serial. This allow diagnosis of processor patches.
-
-    if (!Pstream::parRun())
-    {
-        // Skip in non-parallel
-        return false;
-    }
-
-    if (isState(outputState::CELL_DATA))
-    {
-        ++nCellData_;
-    }
-    else
-    {
-        reportBadState(FatalErrorInFunction, outputState::CELL_DATA)
-            << " for patchID field" << nl << endl
-            << exit(FatalError);
-    }
-
-    label nFaces = nLocalFaces_;
-
-    if (parallel_)
-    {
-        reduce(nFaces, sumOp<label>());
-    }
-
-
-    this->beginDataArray<label>("procID", nFaces);
-
-    bool good = false;
-
-    if (parallel_)
-    {
-        globalIndex procSizes(nLocalFaces_);
-
-        if (Pstream::master())
-        {
-            // Per-processor ids
-            for (const int proci : Pstream::allProcs())
-            {
-                vtk::write(format(), label(proci), procSizes.localSize(proci));
-            }
-
-            good = true;
-        }
-    }
-    else
-    {
-        vtk::write(format(), label(Pstream::myProcNo()), nLocalFaces_);
-
-        good = true;
-    }
-
-
-    this->endDataArray();
-
-    // MPI barrier
-    return parallel_ ? returnReduce(good, orOp<bool>()) : good;
+    return vtk::fileWriter::writeProcIDs(nLocalPolys_);
 }
 
 
@@ -668,15 +610,15 @@ bool Foam::vtk::patchMeshWriter::writeNeighIDs()
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    label nFaces = nLocalFaces_;
+    label nPolys = nLocalPolys_;
 
     if (parallel_)
     {
-        reduce(nFaces, sumOp<label>());
+        reduce(nPolys, sumOp<label>());
     }
 
 
-    this->beginDataArray<label>("neighID", nFaces);
+    this->beginDataArray<label>("neighID", nPolys);
 
     bool good = false;
 
diff --git a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.H b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.H
index 2a26482be1..0622d66d80 100644
--- a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.H
+++ b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,8 +82,8 @@ protected:
         //- Local number of points
         label nLocalPoints_;
 
-        //- Local number of faces
-        label nLocalFaces_;
+        //- Local number of polys (faces)
+        label nLocalPolys_;
 
         //- Local face vertices (connectivity) count. Sum of face sizes.
         label nLocalVerts_;
@@ -97,7 +97,7 @@ protected:
 
     // Private Member Functions
 
-        //- Determine sizes (nLocalPoints_, nLocalFaces_, nLocalVerts_),
+        //- Determine sizes (nLocalPoints_, nLocalPolys_),
         //- and begin piece.
         void beginPiece();
 
@@ -172,7 +172,7 @@ public:
         }
 
         //- The patch IDs
-        inline const labelList& patchIDs() const
+        const labelList& patchIDs() const noexcept
         {
             return patchIDs_;
         }
diff --git a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriterTemplates.C b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriterTemplates.C
index ef81dcb782..349a6b28fc 100644
--- a/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriterTemplates.C
+++ b/src/meshTools/output/vtk/patch/foamVtkPatchMeshWriterTemplates.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,23 +37,33 @@ void Foam::vtk::patchMeshWriter::writeUniform
     const Type& val
 )
 {
+    label nValues(0);
+
     if (isState(outputState::CELL_DATA))
     {
         ++nCellData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfCells_);
+        nValues = nLocalPolys_;
     }
     else if (isState(outputState::POINT_DATA))
     {
         ++nPointData_;
-        vtk::fileWriter::writeUniform<Type>(fieldName, val, numberOfPoints_);
+        nValues = nLocalPoints_;
     }
     else
     {
-        WarningInFunction
-            << "Ignore bad writer state (" << stateNames[state_]
-            << ") for field " << fieldName << nl << endl
+        reportBadState
+        (
+            FatalErrorInFunction,
+            outputState::CELL_DATA,
+            outputState::POINT_DATA
+        )
+            << " for uniform field " << fieldName << nl << endl
             << exit(FatalError);
+
+        return;
     }
+
+    vtk::fileWriter::writeUniform<Type>(fieldName, val, nValues);
 }
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBase.H b/src/sampling/surface/isoSurface/isoSurfaceBase.H
index 3e29657bd1..a856f06f26 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceBase.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceBase.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -127,26 +127,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Check for tet values above/below given (iso) value
-        //  Result encoded as a single integer
-        inline static constexpr int getTetCutIndex
-        (
-            const scalar a,
-            const scalar b,
-            const scalar c,
-            const scalar d,
-            const scalar isoval
-        ) noexcept
-        {
-            return
-            (
-                (a < isoval ? 1 : 0)
-              | (b < isoval ? 2 : 0)
-              | (c < isoval ? 4 : 0)
-              | (d < isoval ? 8 : 0)
-            );
-        }
-
         //- Count the number of cuts matching the mask type
         //  Checks as bitmask or as zero.
         static label countCutType
diff --git a/src/sampling/surface/isoSurface/isoSurfaceParams.C b/src/sampling/surface/isoSurface/isoSurfaceParams.C
index cf9494ed5d..60b003dcb4 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceParams.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceParams.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,10 +52,12 @@ const Foam::Enum
 Foam::isoSurfaceParams::filterNames
 ({
     { filterType::NONE, "none" },
-    { filterType::CELL, "cell" },
-    { filterType::DIAGCELL, "diagcell" },
     { filterType::PARTIAL, "partial" },
     { filterType::FULL, "full" },
+    { filterType::CLEAN, "clean" },
+
+    { filterType::CELL, "cell" },
+    { filterType::DIAGCELL, "diagcell" },
 });
 
 
@@ -137,6 +139,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
 :
     algo_(algo),
     filter_(filter),
+    snap_(true),
     mergeTol_(1e-6),
     clipBounds_(boundBox::invertedBox)
 {}
@@ -152,6 +155,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
 {
     algo_ = getAlgorithmType(dict, algo_);
     filter_ = getFilterType(dict, filter_);
+    snap_ = dict.getOrDefault("snap", true);
     dict.readIfPresent("mergeTol", mergeTol_);
     dict.readIfPresent("bounds", clipBounds_);
 }
diff --git a/src/sampling/surface/isoSurface/isoSurfaceParams.H b/src/sampling/surface/isoSurface/isoSurfaceParams.H
index 1b794f7059..d73da00c27 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceParams.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceParams.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,12 +35,22 @@ Description
         isoMethod   | Algorithm (cell/topo/point/default)   | no  | default
         regularise  | Face simplification (enum or bool)    | no  | true
         mergeTol    | Point merge tolerance (cell/point)    | no  | 1e-6
+        snap        | Point snapping (topo)                 | no  | true
         bounds      | Optional clip bounds                  | no  | inverted
     \endtable
 
     The default algorithm denotes the use of the current \em standard
     algorithm.
 
+Filtering types (for topological iso-surface)
+    - \c none : leave tet cuts untouched
+    - \c partial , \c cell : Combine intra-cell faces
+    - \c full , \c diagcell : Perform \c partial and remove face-diagonal
+      points
+    - \c clean : Perform \c full and eliminate open edges as well.
+      (<b>May cause excessive erosion!</b>)
+    .
+
 SourceFiles
     isoSurfaceParams.C
 
@@ -85,8 +95,12 @@ public:
             NONE = 0,   //!< No filtering
             CELL,       //!< Remove pyramid edge points
             DIAGCELL,   //!< Remove pyramid edge points, face-diagonals
+            NONMANIFOLD, //!< Remove pyramid edge points, face-diagonals
+                         //!< and non-manifold faces
+
             PARTIAL = CELL,     //!< Same as CELL
             FULL = DIAGCELL,    //!< Same as DIAGCELL
+            CLEAN = NONMANIFOLD //!< Same as NONMANIFOLD
         };
 
 
@@ -100,6 +114,9 @@ private:
         //- Filtering for iso-surface faces/points
         filterType filter_;
 
+        //- Point snapping enabled
+        bool snap_;
+
         //- Merge tolerance for cell/point (default: 1e-6)
         scalar mergeTol_;
 
@@ -186,6 +203,18 @@ public:
             filter_ = fltr;
         }
 
+        //- Get point snapping flag
+        bool snap() const noexcept
+        {
+            return snap_;
+        }
+
+        //- Set point snapping flag
+        void snap(bool on) noexcept
+        {
+            snap_ = on;
+        }
+
         //- Get current merge tolerance
         scalar mergeTol() const noexcept
         {
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index e8ff86266f..c0fe1a8077 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -29,13 +29,16 @@ License
 #include "isoSurfaceTopo.H"
 #include "polyMesh.H"
 #include "volFields.H"
+#include "edgeHashes.H"
 #include "tetCell.H"
-#include "tetMatcher.H"
 #include "tetPointRef.H"
 #include "DynamicField.H"
 #include "syncTools.H"
 #include "uindirectPrimitivePatch.H"
 #include "polyMeshTetDecomposition.H"
+#include "foamVtkInternalMeshWriter.H"
+#include "foamVtkLineWriter.H"
+#include "foamVtkSurfaceWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,452 +54,624 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// Get/set snapIndex (0, 1 or 2) at given position
+// 0 = no snap
+// 1 = snap to first edge end
+// 2 = snap to second edge end
+// NB: 4 lower bits left free for regular tet-cut information
 
-Foam::label Foam::isoSurfaceTopo::generatePoint
+#undef  SNAP_END_VALUE
+#undef  SNAP_END_ENCODE
+#define SNAP_END_ENCODE(pos, val)   (((val) << (4 + 2 * pos)))
+#define SNAP_END_VALUE(pos, val)    (((val) >> (4 + 2 * pos)) & 0x3)
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Check for tet values above/below given (iso) value
+// Result encoded as an integer, with possible snapping information too
+inline static int getTetCutIndex
 (
-    const label facei,
-    const bool edgeIsDiag,
-    const edge& vertices,
+    scalar p0,
+    scalar p1,
+    scalar p2,
+    scalar p3,
+    const scalar val,
+    const bool doSnap
+) noexcept
+{
+    int cutIndex
+    (
+        (p0 < val ? 1 : 0)  // point 0
+      | (p1 < val ? 2 : 0)  // point 1
+      | (p2 < val ? 4 : 0)  // point 2
+      | (p3 < val ? 8 : 0)  // point 3
+    );
 
-    DynamicList<edge>& pointToVerts,
-    DynamicList<label>& pointToFace,
-    DynamicList<bool>& pointFromDiag,
-    EdgeMap<label>& vertsToPoint
-) const
+    if (doSnap && cutIndex && cutIndex != 0xF)
+    {
+        // Not all below or all
+
+        // Calculate distances (for snapping)
+        p0 -= val; if (cutIndex & 1) p0 *= -1;
+        p1 -= val; if (cutIndex & 2) p1 *= -1;
+        p2 -= val; if (cutIndex & 4) p2 *= -1;
+        p3 -= val; if (cutIndex & 8) p3 *= -1;
+
+        // Add snap index into regular edge cut index
+        // Snap to end if less than approx 1% of the distance.
+        // - only valid if there is also a corresponding sign change
+        #undef  ADD_SNAP_INDEX
+        #define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2)                \
+        switch (cutIndex & (idx1 | idx2))                              \
+        {                                                              \
+            case idx1 : /* first below, second above */                \
+            case idx2 : /* first above, second below */                \
+                cutIndex |= SNAP_END_ENCODE                            \
+                (                                                      \
+                    pos,                                               \
+                    ((d1 * 100 < d2) ? 1 : (d2 * 100 < d1) ? 2 : 0)    \
+                );                                                     \
+                break;                                                 \
+        }
+
+        ADD_SNAP_INDEX(0, p0, p1, 1, 2);    // Edge 0: 0 -> 1
+        ADD_SNAP_INDEX(1, p0, p2, 1, 4);    // Edge 1: 0 -> 2
+        ADD_SNAP_INDEX(2, p0, p3, 1, 8);    // Edge 2: 0 -> 3
+        ADD_SNAP_INDEX(3, p3, p1, 8, 2);    // Edge 3: 3 -> 1
+        ADD_SNAP_INDEX(4, p1, p2, 2, 4);    // Edge 4: 1 -> 2
+        ADD_SNAP_INDEX(5, p3, p2, 8, 4);    // Edge 5: 3 -> 2
+        #undef ADD_SNAP_INDEX
+    }
+
+    return cutIndex;
+}
+
+
+// Append three labels to list.
+// Filter out degenerate (eg, snapped) tris. Flip face as requested
+inline static void appendTriLabels
+(
+    DynamicList<label>& verts,
+    const label a,
+    const label b,
+    const label c,
+    const bool flip  // Flip normals
+)
+{
+    if (a != b && b != c && c != a)
+    {
+        verts.append(a);
+        if (flip)
+        {
+            verts.append(c);
+            verts.append(b);
+        }
+        else
+        {
+            verts.append(b);
+            verts.append(c);
+        }
+    }
+}
+
+
+// Return point reference to mesh points or cell-centres
+inline static const point& getMeshPointRef
+(
+    const polyMesh& mesh,
+    const label pointi
+)
+{
+    return
+    (
+        pointi < mesh.nPoints()
+      ? mesh.points()[pointi]
+      : mesh.cellCentres()[pointi - mesh.nPoints()]
+    );
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoSurfaceTopo::tetCutAddressing::tetCutAddressing
+(
+    const label nCutCells,
+    const bool useSnap,
+    const bool useDebugCuts
+)
+:
+    vertsToPointLookup_(12*nCutCells),
+    snapVertsLookup_(0),
+
+    pointToFace_(10*nCutCells),
+    pointFromDiag_(10*nCutCells),
+
+    pointToVerts_(10*nCutCells),
+    cutPoints_(12*nCutCells),
+
+    debugCutTets_(),
+    debugCutTetsOn_(useDebugCuts)
+{
+    // Per cell: 5 pyramids cut, each generating 2 triangles
+
+    // Per cell: number of intersected edges:
+    //          - four faces cut so 4 mesh edges + 4 face-diagonal edges
+    //          - 4 of the pyramid edges
+
+    if (useSnap)
+    {
+        // Some, but not all, cells may have point snapping
+        snapVertsLookup_.resize(4*nCutCells);
+    }
+    if (debugCutTetsOn_)
+    {
+        debugCutTets_.reserve(6*nCutCells);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::isoSurfaceTopo::tetCutAddressing::clearDebug()
+{
+    debugCutTets_.clearStorage();
+}
+
+
+void Foam::isoSurfaceTopo::tetCutAddressing::clearDiagonal()
+{
+    pointToFace_.clearStorage();
+    pointFromDiag_.clearStorage();
+}
+
+
+void Foam::isoSurfaceTopo::tetCutAddressing::clearHashes()
+{
+    vertsToPointLookup_.clear();
+    snapVertsLookup_.clear();
+}
+
+
+Foam::label Foam::isoSurfaceTopo::tetCutAddressing::generatePoint
+(
+    label facei,
+    bool edgeIsDiagonal,
+    const int snapEnd,
+    const edge& vertices
+)
 {
     // Generate new point, unless it already exists for edge
+    // or corresponds to a snapped point (from another edge)
 
-    label pointi = vertsToPoint.lookup(vertices, -1);
+    label pointi = vertsToPointLookup_.lookup(vertices, -1);
     if (pointi == -1)
     {
-        pointi = pointToVerts.size();
+        bool addNewPoint(true);
 
-        pointToVerts.append(vertices);
-        pointToFace.append(facei);
-        pointFromDiag.append(edgeIsDiag);
-        vertsToPoint.insert(vertices, pointi);
+        const label snapPointi =
+        (
+            (snapEnd == 1) ? vertices.first()
+          : (snapEnd == 2) ? vertices.second()
+          : -1
+        );
+
+        if (snapPointi == -1)
+        {
+            // No snapped point
+            pointi = pointToVerts_.size();
+            pointToVerts_.append(vertices);
+        }
+        else
+        {
+            // Snapped point. No corresponding face or diagonal
+            facei = -1;
+            edgeIsDiagonal = false;
+
+            pointi = snapVertsLookup_.lookup(snapPointi, -1);
+            addNewPoint = (pointi == -1);
+            if (addNewPoint)
+            {
+                pointi = pointToVerts_.size();
+                snapVertsLookup_.insert(snapPointi, pointi);
+                pointToVerts_.append(edge(snapPointi, snapPointi));
+            }
+        }
+
+        if (addNewPoint)
+        {
+            pointToFace_.append(facei);
+            pointFromDiag_.append(edgeIsDiagonal);
+        }
+
+        vertsToPointLookup_.insert(vertices, pointi);
     }
 
     return pointi;
 }
 
 
-void Foam::isoSurfaceTopo::generateTriPoints
+bool Foam::isoSurfaceTopo::tetCutAddressing::generatePoints
 (
     const label facei,
     const int tetCutIndex,
     const tetCell& tetLabels,
-    const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
 
-    DynamicList<edge>& pointToVerts,
-    DynamicList<label>& pointToFace,
-    DynamicList<bool>& pointFromDiag,
-
-    EdgeMap<label>& vertsToPoint,
-    DynamicList<label>& verts       // Every three verts is new triangle
-) const
+    // Per tet edge whether is face diag etc
+    const FixedList<bool, 6>& edgeIsDiagonal
+)
 {
+    bool flip(false);
+    const label nCutPointsOld(cutPoints_.size());
+
     // Form the vertices of the triangles for each case
-    switch (tetCutIndex)
+    switch (tetCutIndex & 0x0F)
     {
         case 0x00:
         case 0x0F:
         break;
 
-        case 0x01:
-        case 0x0E:
+        // Cut point 0
+        case 0x0E: flip = true; [[fallthrough]];    // Point 0 above cut
+        case 0x01:                                  // Point 0 below cut
         {
-            verts.append
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[0],
-                    tetLabels.edge(0),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[0],
+                    SNAP_END_VALUE(0, tetCutIndex),
+                    tetLabels.edge(0)  // 0 -> 1
                 )
             );
-            verts.append
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[1],
-                    tetLabels.edge(1),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[1],
+                    SNAP_END_VALUE(1, tetCutIndex),
+                    tetLabels.edge(1)  // 0 -> 2
                 )
             );
-            verts.append
+            const label cutC
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[2],
-                    tetLabels.edge(2),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[2],
+                    SNAP_END_VALUE(2, tetCutIndex),
+                    tetLabels.edge(2)  // 0 -> 3
                 )
             );
 
-            if (tetCutIndex == 0x0E)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
         }
         break;
 
-        case 0x02:
-        case 0x0D:
+        // Cut point 1
+        case 0x0D: flip = true; [[fallthrough]];    // Point 1 above cut
+        case 0x02:                                  // Point 1 below cut
         {
-            verts.append
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[0],
-                    tetLabels.reverseEdge(0),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[0],
+                    SNAP_END_VALUE(0, tetCutIndex),
+                    tetLabels.edge(0)  // 0 -> 1
                 )
             );
-            verts.append
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[3],
-                    tetLabels.reverseEdge(3),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[3],
+                    SNAP_END_VALUE(3, tetCutIndex),
+                    tetLabels.edge(3)  // 3 -> 1
                 )
             );
-            verts.append
+            const label cutC
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[4],
-                    tetLabels.edge(4),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[4],
+                    SNAP_END_VALUE(4, tetCutIndex),
+                    tetLabels.edge(4)  // 1 -> 2
                 )
             );
 
-            if (tetCutIndex == 0x0D)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
         }
         break;
 
-        case 0x03:
-        case 0x0C:
+        // Cut point 0/1 | 2/3
+        case 0x0C: flip = true; [[fallthrough]];    // Point 0/1 above cut
+        case 0x03:                                  // Point 0/1 below cut
         {
-            const label p0p2
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[1],
-                    tetLabels.edge(1),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[1],
+                    SNAP_END_VALUE(1, tetCutIndex),
+                    tetLabels.edge(1)  // 0 -> 2
                 )
             );
-            const label p1p3
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[3],
-                    tetLabels.reverseEdge(3),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[2],
+                    SNAP_END_VALUE(2, tetCutIndex),
+                    tetLabels.edge(2)  // 0 -> 3
+                )
+            );
+            const label cutC
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[3],
+                    SNAP_END_VALUE(3, tetCutIndex),
+                    tetLabels.edge(3)  // 3 -> 1
+                )
+            );
+            const label cutD
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[4],
+                    SNAP_END_VALUE(4, tetCutIndex),
+                    tetLabels.edge(4)  // 1 -> 2
                 )
             );
 
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[2],
-                    tetLabels.edge(2),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-            verts.append(p1p3);
-            verts.append(p0p2);
-            verts.append(p1p3);
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[4],
-                    tetLabels.edge(4),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-            verts.append(p0p2);
-
-            if (tetCutIndex == 0x0C)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-5], verts[sz-4]);
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
+            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
         }
         break;
 
-        case 0x04:
-        case 0x0B:
+        // Cut point 2
+        case 0x0B: flip = true; [[fallthrough]];    // Point 2 above cut
+        case 0x04:                                  // Point 2 below cut
         {
-            verts.append
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[1],
-                    tetLabels.reverseEdge(1),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[1],
+                    SNAP_END_VALUE(1, tetCutIndex),
+                    tetLabels.edge(1)  // 0 -> 2
                 )
             );
-            verts.append
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[4],
-                    tetLabels.reverseEdge(4),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[4],
+                    SNAP_END_VALUE(4, tetCutIndex),
+                    tetLabels.edge(4)  // 1 -> 2
                 )
             );
-            verts.append
+            const label cutC
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[5],
-                    tetLabels.reverseEdge(5),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[5],
+                    SNAP_END_VALUE(5, tetCutIndex),
+                    tetLabels.edge(5)  // 3 -> 2
                 )
             );
 
-            if (tetCutIndex == 0x0B)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
         }
         break;
 
-        case 0x05:
-        case 0x0A:
+        // Cut point 0/2 | 1/3
+        case 0x0A: flip = true; [[fallthrough]];    // Point 0/2 above cut
+        case 0x05:                                  // Point 0/2 below cut
         {
-            const label p0p1
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[0],
-                    tetLabels.edge(0),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[0],
+                    SNAP_END_VALUE(0, tetCutIndex),
+                    tetLabels.edge(0)  // 0 -> 1
                 )
             );
-            const label p2p3
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[5],
-                    tetLabels.reverseEdge(5),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[4],
+                    SNAP_END_VALUE(4, tetCutIndex),
+                    tetLabels.edge(4)  // 1 -> 2
+                )
+            );
+            const label cutC
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[5],
+                    SNAP_END_VALUE(5, tetCutIndex),
+                    tetLabels.edge(5)  // 3 -> 2
+                )
+            );
+            const label cutD
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[2],
+                    SNAP_END_VALUE(2, tetCutIndex),
+                    tetLabels.edge(2)  // 0 -> 3
                 )
             );
 
-            verts.append(p0p1);
-            verts.append(p2p3);
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[2],
-                    tetLabels.edge(2),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-            verts.append(p0p1);
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[4],
-                    tetLabels.edge(4),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-            verts.append(p2p3);
-
-            if (tetCutIndex == 0x0A)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-5], verts[sz-4]);
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
+            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
         }
         break;
 
-        case 0x06:
-        case 0x09:
+        // Cut point 1/2 | 0/3
+        case 0x09: flip = true; [[fallthrough]];    // Point 1/2 above cut
+        case 0x06:                                  // Point 1/2 below cut
         {
-            const label p0p1
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[0],
-                    tetLabels.edge(0),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[0],
+                    SNAP_END_VALUE(0, tetCutIndex),
+                    tetLabels.edge(0)  // 0 -> 1
                 )
             );
-            const label p2p3
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[5],
-                    tetLabels.reverseEdge(5),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[3],
+                    SNAP_END_VALUE(3, tetCutIndex),
+                    tetLabels.edge(3)  // 3 -> 1
+                )
+            );
+            const label cutC
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[5],
+                    SNAP_END_VALUE(5, tetCutIndex),
+                    tetLabels.edge(5)  // 3 -> 2
+                )
+            );
+            const label cutD
+            (
+                generatePoint
+                (
+                    facei,
+                    edgeIsDiagonal[1],
+                    SNAP_END_VALUE(1, tetCutIndex),
+                    tetLabels.edge(1)  // 0 -> 2
                 )
             );
 
-            verts.append(p0p1);
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[3],
-                    tetLabels.reverseEdge(3),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-            verts.append(p2p3);
-            verts.append(p0p1);
-            verts.append(p2p3);
-            verts.append
-            (
-                generatePoint
-                (
-                    facei,
-                    edgeIsDiag[1],
-                    tetLabels.edge(1),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
-                )
-            );
-
-            if (tetCutIndex == 0x09)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-5], verts[sz-4]);
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
+            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
         }
         break;
 
-        case 0x08:
-        case 0x07:
+        // Cut point 3
+        case 0x07: flip = true; [[fallthrough]];    // Point 3 above cut
+        case 0x08:                                  // Point 3 below cut
         {
-            verts.append
+            const label cutA
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[2],
-                    tetLabels.reverseEdge(2),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[2],
+                    SNAP_END_VALUE(2, tetCutIndex),
+                    tetLabels.edge(2)  // 0 -> 3
                 )
             );
-            verts.append
+            const label cutB
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[5],
-                    tetLabels.reverseEdge(5),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[5],
+                    SNAP_END_VALUE(5, tetCutIndex),
+                    tetLabels.edge(5)  // 3 -> 2
                 )
             );
-            verts.append
+            const label cutC
             (
                 generatePoint
                 (
                     facei,
-                    edgeIsDiag[3],
-                    tetLabels.edge(3),
-                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+                    edgeIsDiagonal[3],
+                    SNAP_END_VALUE(3, tetCutIndex),
+                    tetLabels.edge(3)  // 3 -> 1
                 )
             );
-            if (tetCutIndex == 0x07)
-            {
-                // Flip normals
-                const label sz = verts.size();
-                std::swap(verts[sz-2], verts[sz-1]);
-            }
+
+            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
         }
         break;
     }
+
+    const bool added(nCutPointsOld != cutPoints_.size());
+
+    if (added && debugCutTetsOn_)
+    {
+        debugCutTets_.append(tetLabels.shape());
+    }
+
+    return added;
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Requires mesh_, tetBasePtIs
 void Foam::isoSurfaceTopo::generateTriPoints
 (
     const label celli,
     const bool isTet,
-
-    DynamicList<edge>& pointToVerts,
-    DynamicList<label>& pointToFace,
-    DynamicList<bool>& pointFromDiag,
-
-    EdgeMap<label>& vertsToPoint,
-    DynamicList<label>& verts,
-    DynamicList<label>& faceLabels
+    const labelList& tetBasePtIs,
+    tetCutAddressing& tetCutAddr
 ) const
 {
     const faceList& faces = mesh_.faces();
     const labelList& faceOwner = mesh_.faceOwner();
     const cell& cFaces = mesh_.cells()[celli];
+    const bool doSnap = this->snap();
 
     if (isTet)
     {
         // For tets don't do cell-centre decomposition, just use the
         // tet points and values
 
-        const label startTrii = verts.size();
-
         const label facei = cFaces[0];
         const face& f0 = faces[facei];
 
         // Get the other point from f1. Tbd: check if not duplicate face
         // (ACMI / ignoreBoundaryFaces_).
         const face& f1 = faces[cFaces[1]];
-        label oppositeI = -1;
+        label apexi = -1;
         forAll(f1, fp)
         {
-            oppositeI = f1[fp];
-            if (!f0.found(oppositeI))
+            apexi = f1[fp];
+            if (!f0.found(apexi))
             {
                 break;
             }
         }
 
-        label p0 = f0[0];
+        const label p0 = f0[0];
         label p1 = f0[1];
         label p2 = f0[2];
 
@@ -505,31 +680,27 @@ void Foam::isoSurfaceTopo::generateTriPoints
             std::swap(p1, p2);
         }
 
-        generateTriPoints
+        const tetCell tetLabels(p0, p1, p2, apexi);
+        const int tetCutIndex
         (
-            facei,
             getTetCutIndex
             (
                 pVals_[p0],
                 pVals_[p1],
                 pVals_[p2],
-                pVals_[oppositeI],
-                iso_
-            ),
-            tetCell(p0, p1, p2, oppositeI),
-            FixedList<bool, 6>(false),  // edgeIsDiag = false
-
-            pointToVerts,
-            pointToFace,
-            pointFromDiag,
-            vertsToPoint,
-            verts       // Every three verts is new triangle
+                pVals_[apexi],
+                iso_,
+                doSnap
+            )
         );
 
-        for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
-        {
-            faceLabels.append(facei);
-        }
+        tetCutAddr.generatePoints
+        (
+            facei,
+            tetCutIndex,
+            tetLabels,
+            FixedList<bool, 6>(false)  // Not face diagonal
+        );
     }
     else
     {
@@ -544,11 +715,9 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 continue;
             }
 
-            const label startTrii = verts.size();
-
             const face& f = faces[facei];
 
-            label fp0 = tetBasePtIs_[facei];
+            label fp0 = tetBasePtIs[facei];
 
             // Fallback
             if (fp0 < 0)
@@ -556,55 +725,48 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 fp0 = 0;
             }
 
+            const label p0 = f[fp0];
             label fp = f.fcIndex(fp0);
             for (label i = 2; i < f.size(); ++i)
             {
-                const label nextFp = f.fcIndex(fp);
-
-                FixedList<bool, 6> edgeIsDiag(false);
-
-                label p0 = f[fp0];
                 label p1 = f[fp];
-                label p2 = f[nextFp];
+                fp = f.fcIndex(fp);
+                label p2 = f[fp];
+
+                FixedList<bool, 6> edgeIsDiagonal(false);
                 if (faceOwner[facei] == celli)
                 {
                     std::swap(p1, p2);
-                    if (i != 2) edgeIsDiag[1] = true;
-                    if (i != f.size()-1) edgeIsDiag[0] = true;
+                    if (i != 2) edgeIsDiagonal[1] = true;
+                    if (i != f.size()-1) edgeIsDiagonal[0] = true;
                 }
                 else
                 {
-                    if (i != 2) edgeIsDiag[0] = true;
-                    if (i != f.size()-1) edgeIsDiag[1] = true;
+                    if (i != 2) edgeIsDiagonal[0] = true;
+                    if (i != f.size()-1) edgeIsDiagonal[1] = true;
                 }
 
-                generateTriPoints
+                const tetCell tetLabels(p0, p1, p2, mesh_.nPoints()+celli);
+                const int tetCutIndex
                 (
-                    facei,
                     getTetCutIndex
                     (
                         pVals_[p0],
                         pVals_[p1],
                         pVals_[p2],
                         cVals_[celli],
-                        iso_
-                    ),
-                    tetCell(p0, p1, p2, mesh_.nPoints()+celli),
-                    edgeIsDiag,
-
-                    pointToVerts,
-                    pointToFace,
-                    pointFromDiag,
-                    vertsToPoint,
-                    verts       // Every three verts is new triangle
+                        iso_,
+                        doSnap
+                    )
                 );
 
-                fp = nextFp;
-            }
-
-            for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
-            {
-                faceLabels.append(facei);
+                tetCutAddr.generatePoints
+                (
+                    facei,
+                    tetCutIndex,
+                    tetLabels,
+                    edgeIsDiagonal
+                );
             }
         }
     }
@@ -708,11 +870,11 @@ void Foam::isoSurfaceTopo::removeInsidePoints
 
     const pointField& points = s.points();
 
-    DynamicList<face> compactFaces(s.size()/8);
+    DynamicList<face> compactFaces(s.size()/4);
 
     for (label celli = 0; celli < start.size()-1; ++celli)
     {
-        // All triangles for the current cell
+        // Triangles for the current cell
 
         const label nTris = start[celli+1]-start[celli];
 
@@ -730,7 +892,6 @@ void Foam::isoSurfaceTopo::removeInsidePoints
                 pp,
                 pointFromDiag,
                 pointToFace,
-                //protectedFace,
                 celli,
 
                 compactFaces,
@@ -791,9 +952,11 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     const bitSet& ignoreCells
 )
 :
-    isoSurfaceBase(mesh, cellValues, pointValues, iso, params),
-    cellCutType_(mesh.nCells(), cutType::UNVISITED)
+    isoSurfaceBase(mesh, cellValues, pointValues, iso, params)
 {
+    // The cell cut type
+    List<cutType> cellCutType_(mesh.nCells(), cutType::UNVISITED);
+
     if (debug)
     {
         Pout<< "isoSurfaceTopo:" << nl
@@ -819,9 +982,11 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     nBlockedCells +=
         blockCells(cellCutType_, params.getClipBounds(), volumeType::OUTSIDE);
 
-
-    // Adjust tet base points to improve tet quality
-    tetBasePtIs_ = polyMeshTetDecomposition::adjustTetBasePtIs(mesh_, debug);
+    // Adjusted tet base points to improve tet quality
+    labelList tetBasePtIs
+    (
+        polyMeshTetDecomposition::adjustTetBasePtIs(mesh_, debug)
+    );
 
 
     // Determine cell cuts
@@ -851,7 +1016,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                 false
             ),
             fvmesh,
-            dimensionedScalar(dimless, Zero)
+            dimensionedScalar(dimless)
         );
 
         auto& debugFld = debugField.primitiveFieldRef();
@@ -861,68 +1026,133 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
             debugFld[celli] = cellCutType_[celli];
         }
 
-        Pout<< "Writing cut types:"
-            << debugField.objectPath() << endl;
-
+        Info<< "Writing cut types: " << debugField.objectRelPath() << nl;
         debugField.write();
     }
 
+    // Additional debugging
+    if (debug & 8)
+    {
+        const word timeDesc
+        (
+            word::printf("%08d", mesh_.time().timeIndex())
+        );
 
-    // Per cell: 5 pyramids cut, each generating 2 triangles
-    //  - pointToVerts : from generated iso point to originating mesh verts
-    DynamicList<edge> pointToVerts(10*nCutCells);
-    //  - pointToFace : from generated iso point to originating mesh face
-    DynamicList<label> pointToFace(10*nCutCells);
-    //  - pointFromDiag : if generated iso point is on face diagonal
-    DynamicList<bool> pointFromDiag(10*nCutCells);
+        // Write debug cuts cells in VTK format
+        {
+            constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
+            labelList debugCutCells(nCutCells, Zero);
 
-    // Per cell: number of intersected edges:
-    //          - four faces cut so 4 mesh edges + 4 face-diagonal edges
-    //          - 4 of the pyramid edges
-    EdgeMap<label> vertsToPoint(12*nCutCells);
-    DynamicList<label> verts(12*nCutCells);
-    // Per cell: 5 pyramids cut (since only one pyramid not cut)
-    DynamicList<label> faceLabels(5*nCutCells);
-    DynamicList<label> cellLabels(5*nCutCells);
+            label nout = 0;
+            forAll(cellCutType_, celli)
+            {
+                if ((cellCutType_[celli] & realCut) != 0)
+                {
+                    debugCutCells[nout] = celli;
+                    ++nout;
+                    if (nout >= nCutCells) break;
+                }
+            }
 
+            // The mesh subset cut
+            vtk::vtuCells vtuCells;
+            vtuCells.reset(mesh_, debugCutCells);
+
+            vtk::internalMeshWriter writer
+            (
+                mesh_,
+                vtuCells,
+                fileName
+                (
+                    mesh_.time().globalPath()
+                  / ("isoSurfaceTopo." + timeDesc + "-cutCells")
+                )
+            );
+
+            writer.writeGeometry();
+
+            // Cell qualities
+            writer.beginCellData();
+            writer.writeProcIDs();
+            writer.writeCellData("cutField", cVals_);
+
+            // Point qualities
+            writer.beginPointData();
+            writer.writePointData("cutField", pVals_);
+
+            Info<< "isoSurfaceTopo : (debug) wrote "
+                << returnReduce(nCutCells, sumOp<label>())
+                << " cut cells: "
+                << writer.output().name() << nl;
+        }
+    }
+
+
+    tetCutAddressing tetCutAddr
+    (
+        nCutCells,
+        this->snap(),
+        (debug & 8)  // Enable debug tets
+    );
 
     labelList startTri(mesh_.nCells()+1, Zero);
-
     for (label celli = 0; celli < mesh_.nCells(); ++celli)
     {
-        startTri[celli] = faceLabels.size();
+        startTri[celli] = tetCutAddr.nFaces();
         if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
         {
             generateTriPoints
             (
                 celli,
-
                 // Same as tetMatcher::test(mesh_, celli),
-                bool(cellCutType_[celli] & cutType::TETCUT),  // isTet
+                bool(cellCutType_[celli] & cutType::TETCUT),
 
-                pointToVerts,
-                pointToFace,
-                pointFromDiag,
-
-                vertsToPoint,
-                verts,
-                faceLabels
+                tetBasePtIs,
+                tetCutAddr
             );
-
-            for (label i = startTri[celli]; i < faceLabels.size(); ++i)
-            {
-                cellLabels.append(celli);
-            }
         }
     }
-    startTri[mesh_.nCells()] = faceLabels.size();
+    startTri.last() = tetCutAddr.nFaces();
+
+    // Information not needed anymore:
+    tetBasePtIs.clear();
+    tetCutAddr.clearHashes();
 
 
-    pointToVerts_.transfer(pointToVerts);
-    meshCells_.transfer(cellLabels);
-    pointToFace_.transfer(pointToFace);
+    // From list of vertices -> triangular faces
+    faceList allTriFaces(startTri.last());
+    {
+        auto& verts = tetCutAddr.cutPoints();
 
-    pointField allPoints
+        label verti = 0;
+        for (face& f : allTriFaces)
+        {
+            f.resize(3);
+            f[0] = verts[verti++];
+            f[1] = verts[verti++];
+            f[2] = verts[verti++];
+        }
+        verts.clearStorage();  // Not needed anymore
+    }
+
+
+    // The cells cut by the triangular faces
+    meshCells_.resize(startTri.last());
+    for (label celli = 0; celli < startTri.size()-1; ++celli)
+    {
+        // All triangles for the current cell
+        labelList::subList
+        (
+            meshCells_,
+            (startTri[celli+1] - startTri[celli]),
+            startTri[celli]
+        ) = celli;
+    }
+
+
+    pointToVerts_.transfer(tetCutAddr.pointToVerts());
+
+    pointField allTriPoints
     (
         this->interpolateTemplate
         (
@@ -933,37 +1163,16 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
 
 
     // Assign to MeshedSurface
-    faceList allTris(faceLabels.size());
-    label verti = 0;
-    for (face& allTri : allTris)
-    {
-        allTri.resize(3);
-        allTri[0] = verts[verti++];
-        allTri[1] = verts[verti++];
-        allTri[2] = verts[verti++];
-    }
-
-
-    surfZoneList allZones(one{}, surfZone("allFaces", allTris.size()));
 
     Mesh::clear();
     Mesh updated
     (
-        std::move(allPoints),
-        std::move(allTris),
-        std::move(allZones)
+        std::move(allTriPoints),
+        std::move(allTriFaces),
+        surfZoneList(one{}, surfZone("allFaces", allTriFaces.size()))
     );
     Mesh::transfer(updated);
 
-    // Now:
-    // - generated faces and points are assigned to *this
-    // - per point we know:
-    //  - pointOnDiag: whether it is on a face-diagonal edge
-    //  - pointToFace_: from what pyramid (cell+face) it was produced
-    //    (note that the pyramid faces are shared between multiple mesh faces)
-    //  - pointToVerts_ : originating mesh vertex or cell centre
-
-
     if (debug)
     {
         Pout<< "isoSurfaceTopo : generated "
@@ -972,6 +1181,16 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     }
 
 
+    // Now:
+    // - generated faces and points are assigned to *this
+    // - per point we know:
+    //  - pointOnDiag: whether it is on a face-diagonal edge
+    //  - pointToFace: from what pyramid (cell+face) it was produced
+    //    (note that the pyramid faces are shared between multiple mesh faces)
+    //  - pointToVerts_ : originating mesh vertex or cell centre
+
+    // Basic filtering
+
     if (params.filter() != filterType::NONE)
     {
         // Triangulate outside (filter edges to cell centres and optionally
@@ -982,22 +1201,21 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
         removeInsidePoints
         (
             *this,
-            (params.filter() == filterType::DIAGCELL ? true : false),
-
-            pointFromDiag,
-            pointToFace_,
+            // Filter face diagonal
+            (
+                params.filter() == filterType::DIAGCELL
+             || params.filter() == filterType::NONMANIFOLD
+            ),
+            tetCutAddr.pointFromDiag(),
+            tetCutAddr.pointToFace(),
             startTri,
             pointCompactMap,
             compactCellIDs
         );
 
         pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
-        pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
         meshCells_.transfer(compactCellIDs);
 
-        pointCompactMap.clearStorage();
-        compactCellIDs.clearStorage();
-
         if (debug)
         {
             Pout<< "isoSurfaceTopo :"
@@ -1008,33 +1226,35 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
         }
     }
 
-    // Not required after this stage
-    pointFromDiag.clearStorage();
+    // Diagonal filter information not needed anymore
+    tetCutAddr.clearDiagonal();
 
 
-    if (params.filter() == filterType::DIAGCELL)
+    // For more advanced filtering (eg, removal of open edges)
+    // need the boundary and other 'protected' points
+
+    bitSet isProtectedPoint;
+    if
+    (
+        (params.filter() == filterType::NONMANIFOLD)
+     || tetCutAddr.debugCutTetsOn()
+    )
     {
-        // We remove verts on face diagonals. This is in fact just
-        // straightening the edges of the face through the cell. This can
-        // close off 'pockets' of triangles and create open or
-        // multiply-connected triangles
+        // Mark points on mesh outside as 'protected'
+        // - never erode these edges
 
-        // Solved by eroding open-edges
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        isProtectedPoint.resize(mesh_.nPoints());
 
-        // Mark points on mesh outside.
-        bitSet isBoundaryPoint(mesh.nPoints());
         for
         (
-            label facei = mesh.nInternalFaces();
-            facei < mesh.nFaces();
+            label facei = mesh_.nInternalFaces();
+            facei < mesh_.nFaces();
             ++facei
         )
         {
-            isBoundaryPoint.set(mesh.faces()[facei]);
+            isProtectedPoint.set(mesh_.faces()[facei]);
         }
 
-
         // Include faces that would be exposed from mesh subset
         if (nBlockedCells)
         {
@@ -1051,11 +1271,226 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                  != (cellCutType_[faceNei[facei]] == cutType::BLOCKED)
                 )
                 {
-                    isBoundaryPoint.set(mesh.faces()[facei]);
+                    isProtectedPoint.set(mesh_.faces()[facei]);
                 }
             }
         }
+    }
 
+    // Initial cell cut information not needed anymore
+    cellCutType_.clear();
+
+
+    // Additional debugging
+    if (tetCutAddr.debugCutTetsOn())
+    {
+        const word timeDesc
+        (
+            word::printf("%08d", mesh_.time().timeIndex())
+        );
+
+        // Write debug cut tets in VTK format
+        {
+            const auto& debugCuts = tetCutAddr.debugCutTets();
+
+            // The TET shapes, using the mesh_ points information
+            vtk::vtuCells vtuCells;
+            vtuCells.resetShapes(debugCuts);
+
+            // Use all points and all cell-centres
+            vtuCells.setNumPoints(mesh_.nPoints());
+            vtuCells.addPointCellLabels(identity(mesh_.nCells()));
+
+            vtk::internalMeshWriter writer
+            (
+                mesh_,
+                vtuCells,
+                fileName
+                (
+                    mesh_.time().globalPath()
+                  / ("isoSurfaceTopo." + timeDesc + "-cutTets")
+                )
+            );
+
+            writer.writeGeometry();
+
+            writer.beginCellData();
+            writer.writeProcIDs();
+
+            // Quality of the cut tets
+            {
+                Field<scalar> cutTetQuality(debugCuts.size());
+                forAll(cutTetQuality, teti)
+                {
+                    cutTetQuality[teti] = tetPointRef
+                    (
+                        getMeshPointRef(mesh_, debugCuts[teti][0]),
+                        getMeshPointRef(mesh_, debugCuts[teti][1]),
+                        getMeshPointRef(mesh_, debugCuts[teti][2]),
+                        getMeshPointRef(mesh_, debugCuts[teti][3])
+                    ).quality();
+                }
+                writer.writeCellData("tetQuality", cutTetQuality);
+            }
+
+            if (this->snap())
+            {
+                labelList snapped(vtuCells.nFieldPoints(), Zero);
+
+                for (const edge& verts : pointToVerts_)
+                {
+                    if (verts.first() == verts.second())
+                    {
+                        // Duplicate index (ie, snapped)
+                        snapped[verts.first()] = 1;
+                    }
+                }
+
+                writer.beginPointData();
+                writer.writePointData("snapped", snapped);
+            }
+
+            Info<< "isoSurfaceTopo : (debug) wrote "
+                << returnReduce(debugCuts.size(), sumOp<label>())
+                << " cut tets: "
+                << writer.output().name() << nl;
+        }
+
+        // Determining open edges. Same logic as used later...
+
+        labelHashSet openEdgeIds(0);
+
+        {
+            const Mesh& s = *this;
+
+            const labelList& mp = s.meshPoints();
+            const edgeList& surfEdges = s.edges();
+            const labelListList& edgeFaces = s.edgeFaces();
+            openEdgeIds.resize(2*s.size());
+
+            forAll(edgeFaces, edgei)
+            {
+                const labelList& eFaces = edgeFaces[edgei];
+                if (eFaces.size() == 1)
+                {
+                    // Open edge (not originating from a boundary face)
+
+                    const edge& e = surfEdges[edgei];
+                    const edge& verts0 = pointToVerts_[mp[e.first()]];
+                    const edge& verts1 = pointToVerts_[mp[e.second()]];
+
+                    if
+                    (
+                        isProtectedPoint.test(verts0.first())
+                     && isProtectedPoint.test(verts0.second())
+                     && isProtectedPoint.test(verts1.first())
+                     && isProtectedPoint.test(verts1.second())
+                    )
+                    {
+                        // Open edge on boundary face. Keep
+                    }
+                    else
+                    {
+                        // Open edge
+                        openEdgeIds.insert(edgei);
+                    }
+                }
+            }
+
+            const label nOpenEdges
+            (
+                returnReduce(openEdgeIds.size(), sumOp<label>())
+            );
+
+            if (nOpenEdges)
+            {
+                const edgeList debugEdges
+                (
+                    surfEdges,
+                    openEdgeIds.sortedToc()
+                );
+
+                vtk::lineWriter writer
+                (
+                    s.points(),
+                    debugEdges,
+                    fileName
+                    (
+                        mesh_.time().globalPath()
+                      / ("isoSurfaceTopo." + timeDesc + "-openEdges")
+                    )
+                );
+
+                writer.writeGeometry();
+
+                writer.beginCellData();
+                writer.writeProcIDs();
+
+                Info<< "isoSurfaceTopo : (debug) wrote "
+                    << nOpenEdges << " open edges: "
+                    << writer.output().name() << nl;
+            }
+            else
+            {
+                Info<< "isoSurfaceTopo : no open edges" << nl;
+            }
+        }
+
+        // Write debug surface with snaps
+        if (this->snap())
+        {
+            const Mesh& s = *this;
+
+            vtk::surfaceWriter writer
+            (
+                s.points(),
+                s,
+                fileName
+                (
+                    mesh_.time().globalPath()
+                  / ("isoSurfaceTopo." + timeDesc + "-surface")
+                )
+            );
+
+            writer.writeGeometry();
+
+            writer.beginCellData();
+            writer.writeProcIDs();
+
+            {
+                labelList snapped(s.nPoints(), Zero);
+                forAll(pointToVerts_, i)
+                {
+                    const edge& verts =pointToVerts_[i];
+                    if (verts.first() == verts.second())
+                    {
+                        // Duplicate index (ie, snapped)
+                        snapped[i] = 1;
+                    }
+                }
+
+                writer.beginPointData();
+                writer.write("snapped", snapped);
+            }
+
+            Info<< "isoSurfaceTopo : (debug) wrote "
+                << returnReduce(s.size(), sumOp<label>())
+                << " faces : "
+                << writer.output().name() << nl;
+        }
+    }
+    tetCutAddr.clearDebug();
+
+
+    if (params.filter() == filterType::NONMANIFOLD)
+    {
+        // We remove verts on face diagonals. This is in fact just
+        // straightening the edges of the face through the cell. This can
+        // close off 'pockets' of triangles and create open or
+        // multiply-connected triangles
+
+        // Solved by eroding open-edges
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         // The list of surface faces that should be retained after erosion
         Mesh& surf = *this;
@@ -1086,19 +1521,18 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                 const labelList& eFaces = edgeFaces[edgei];
                 if (eFaces.size() == 1)
                 {
-                    // Open edge. Check that vertices do not originate
-                    // from a boundary face
-                    const edge& e = surfEdges[edgei];
+                    // Open edge (not originating from a boundary face)
 
+                    const edge& e = surfEdges[edgei];
                     const edge& verts0 = pointToVerts_[mp[e.first()]];
                     const edge& verts1 = pointToVerts_[mp[e.second()]];
 
                     if
                     (
-                        isBoundaryPoint.test(verts0.first())
-                     && isBoundaryPoint.test(verts0.second())
-                     && isBoundaryPoint.test(verts1.first())
-                     && isBoundaryPoint.test(verts1.second())
+                        isProtectedPoint.test(verts0.first())
+                     && isProtectedPoint.test(verts0.second())
+                     && isProtectedPoint.test(verts1.first())
+                     && isProtectedPoint.test(verts1.second())
                     )
                     {
                         // Open edge on boundary face. Keep
@@ -1160,7 +1594,6 @@ void Foam::isoSurfaceTopo::inplaceSubsetMesh(const bitSet& include)
     meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
 
     pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
-    pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
 }
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index 29fd36d8d4..cd5886d135 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -50,6 +50,7 @@ namespace Foam
 {
 
 // Forward Declarations
+class cellShape;
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -62,63 +63,103 @@ class isoSurfaceTopo
 {
     // Private Data
 
-        //- Corrected version of tetBasePtIs
-        labelList tetBasePtIs_;
-
-        //- Per point: originating mesh vertex/cc. See encoding above
+        //- Per point: originating mesh vertex/cell centre combination.
+        //  Vertices less than nPoints are mesh vertices,
+        //  duplicate vertices indicate a snapped point
         edgeList pointToVerts_;
 
-        //- For every point the originating face in mesh
-        labelList pointToFace_;
 
-        //- The cell cut type
-        List<cutType> cellCutType_;
+    // Private Classes
+
+        //- Handling, bookkeeping for tet cuts
+        class tetCutAddressing
+        {
+            // Bookkeeping hashes used during construction
+            EdgeMap<label> vertsToPointLookup_;
+            Map<label> snapVertsLookup_;
+
+            // Filter information for face diagonals
+            DynamicList<label> pointToFace_;
+            DynamicList<bool> pointFromDiag_;
+
+            // Final output
+            DynamicList<edge> pointToVerts_;
+            DynamicList<label> cutPoints_;
+
+            //- List of cut (decomposed) cell tets. Debug use only.
+            DynamicList<cellShape> debugCutTets_;
+
+            bool debugCutTetsOn_;
+
+
+        public:
+
+        // Constructors
+
+            //- Construct with reserved sizes
+            tetCutAddressing
+            (
+                const label nCutCells,
+                const bool useSnap,
+                const bool useDebugCuts = false
+            );
+
+
+        // Member Functions
+
+            //- Effective number of faces
+            label nFaces() const { return cutPoints_.size()/3; }
+
+            DynamicList<label>& pointToFace() { return pointToFace_; }
+            DynamicList<bool>& pointFromDiag() { return pointFromDiag_; }
+
+            DynamicList<edge>& pointToVerts() { return pointToVerts_; }
+            DynamicList<label>& cutPoints() { return cutPoints_; }
+            DynamicList<cellShape>& debugCutTets() { return debugCutTets_; }
+
+            //- Number of debug cut tets
+            label nDebugTets() const { return debugCutTets_.size(); }
+
+            //- Debugging cut tets active
+            bool debugCutTetsOn() const { return debugCutTetsOn_; }
+
+            void clearDebug();
+            void clearDiagonal();
+            void clearHashes();
+
+            //- Generate single point on edge
+            label generatePoint
+            (
+                label facei,          //!< Originating mesh face for cut-point
+                bool edgeIsDiagonal,  //!< Edge on face diagonal
+
+                // 0: no snap, 1: snap first, 2: snap second
+                const int snapEnd,
+
+                const edge& vertices
+            );
+
+            //- Generate triangles from tet cut
+            bool generatePoints
+            (
+                const label facei,
+                const int tetCutIndex,  //!< Encoded tet cut + tet snapping
+                const tetCell& tetLabels,
+                const FixedList<bool, 6>& edgeIsDiagonal
+            );
+        };
 
 
     // Private Member Functions
 
-        //- Generate single point on edge
-        label generatePoint
-        (
-            const label facei,
-            const bool edgeIsDiag,
-            const edge& vertices,
-
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
-            EdgeMap<label>& vertsToPoint
-        ) const;
-
-        //- Generate triangles from tet
-        void generateTriPoints
-        (
-            const label facei,
-            const int tetCutIndex,  //!< Encoded tet cuts. getTetCutIndex()
-            const tetCell& tetLabels,
-            const FixedList<bool, 6>& edgeIsDiag,
-
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
-
-            EdgeMap<label>& vertsToPoint,
-            DynamicList<label>& verts
-        ) const;
-
-        //- Generate triangles from cell
+        //- Generate triangle points from cell
         void generateTriPoints
         (
             const label celli,
             const bool isTet,
+            const labelList& tetBasePtIs,
 
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
-
-            EdgeMap<label>& vertsToPoint,
-            DynamicList<label>& verts,
-            DynamicList<label>& faceLabels
+            tetCutAddressing& tetCutAddr
         ) const;
 
 
@@ -165,12 +206,12 @@ protected:
 
     // Sampling
 
-        //- Interpolates cCoords,pCoords.
+        //- Interpolates cellData and pointData fields
         template<class Type>
         tmp<Field<Type>> interpolateTemplate
         (
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
+            const Field<Type>& cellData,
+            const Field<Type>& pointData
         ) const;
 
 public:
@@ -205,14 +246,10 @@ public:
 
     // Member Functions
 
-        //- For every point originating face (pyramid) in mesh
-        const labelList& pointToFace() const
-        {
-            return pointToFace_;
-        }
-
-        //- Per point: originating mesh vertex/cc. See encoding above
-        const edgeList& pointToVerts() const
+        //- Per point: originating mesh vertex/cell centre combination.
+        //  Vertices less than nPoints are mesh vertices (encoding above),
+        //  duplicate vertices indicate a snapped point
+        const edgeList& pointToVerts() const noexcept
         {
             return pointToVerts_;
         }
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
index 8ffec74e8d..5056fb38a5 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,8 +32,8 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::isoSurfaceTopo::interpolateTemplate
 (
-    const Field<Type>& cellCoords,
-    const Field<Type>& pointCoords
+    const Field<Type>& cellData,
+    const Field<Type>& pointData
 ) const
 {
     auto tfld = tmp<Field<Type>>::New(pointToVerts_.size());
@@ -41,41 +41,50 @@ Foam::isoSurfaceTopo::interpolateTemplate
 
     forAll(pointToVerts_, i)
     {
+        const edge& verts = pointToVerts_[i];
+        Type& val = fld[i];
+
         scalar s0;
-        Type p0;
+        Type v0;
         {
-            label idx = pointToVerts_[i].first();
+            label idx = verts.first();
             if (idx < mesh_.nPoints())
             {
                 // Point index
                 s0 = pVals_[idx];
-                p0 = pointCoords[idx];
+                v0 = pointData[idx];
             }
             else
             {
                 // Cell index
                 idx -= mesh_.nPoints();
                 s0 = cVals_[idx];
-                p0 = cellCoords[idx];
+                v0 = cellData[idx];
             }
         }
 
         scalar s1;
-        Type p1;
+        Type v1;
         {
-            label idx = pointToVerts_[i].second();
-            if (idx < mesh_.nPoints())
+            label idx = verts.second();
+            if (idx == verts.first())
+            {
+                // Duplicate index (ie, snapped)
+                val = v0;
+                continue;
+            }
+            else if (idx < mesh_.nPoints())
             {
                 // Point index
                 s1 = pVals_[idx];
-                p1 = pointCoords[idx];
+                v1 = pointData[idx];
             }
             else
             {
                 // Cell index
                 idx -= mesh_.nPoints();
                 s1 = cVals_[idx];
-                p1 = cellCoords[idx];
+                v1 = cellData[idx];
             }
         }
 
@@ -83,11 +92,11 @@ Foam::isoSurfaceTopo::interpolateTemplate
         if (mag(d) > VSMALL)
         {
             const scalar s = (iso_-s0)/d;
-            fld[i] = s*p1+(1.0-s)*p0;
+            val = s*v1+(1.0-s)*v0;
         }
         else
         {
-            fld[i] = 0.5*(p0+p1);
+            val = 0.5*(v0+v1);
         }
     }