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/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
+    );
 }