diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.C b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.C
index bade666a7a..c2b9793dc1 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.C
@@ -38,6 +38,7 @@ License
 #include "vtkTextActor.h"
 #include "vtkTextProperty.h"
 #include "vtkSmartPointer.h"
+#include "vtkInformation.h"
 
 // Templates (only needed here)
 #include "vtkPVFoamUpdateTemplates.C"
@@ -117,6 +118,70 @@ void Foam::vtkPVFoam::resetCounters()
 }
 
 
+template<class Container>
+bool Foam::vtkPVFoam::addOutputBlock
+(
+    vtkMultiBlockDataSet* output,
+    const HashTable<Container, string>& cache,
+    const arrayRange& selector,
+    const bool singleDataset
+) const
+{
+    const auto blockNo = output->GetNumberOfBlocks();
+    vtkSmartPointer<vtkMultiBlockDataSet> block;
+    int datasetNo = 0;
+
+    for (auto partId : selector)
+    {
+        if (selectedPartIds_.found(partId))
+        {
+            const auto& longName = selectedPartIds_[partId];
+            const word shortName = getPartName(partId);
+
+            auto iter = cache.find(longName);
+            if (iter.found() && (*iter).vtkmesh)
+            {
+                auto dataset = (*iter).vtkmesh;
+
+                if (singleDataset)
+                {
+                    output->SetBlock(blockNo, dataset);
+                    output->GetMetaData(blockNo)->Set
+                    (
+                        vtkCompositeDataSet::NAME(),
+                        shortName.c_str()
+                    );
+
+                    ++datasetNo;
+                    break;
+                }
+                else if (datasetNo == 0)
+                {
+                    block = vtkSmartPointer<vtkMultiBlockDataSet>::New();
+                    output->SetBlock(blockNo, block);
+                    output->GetMetaData(blockNo)->Set
+                    (
+                        vtkCompositeDataSet::NAME(),
+                        selector.name()
+                    );
+                }
+
+                block->SetBlock(datasetNo, dataset);
+                block->GetMetaData(datasetNo)->Set
+                (
+                    vtkCompositeDataSet::NAME(),
+                    shortName.c_str()
+                );
+
+                ++datasetNo;
+            }
+        }
+    }
+
+    return datasetNo;
+}
+
+
 int Foam::vtkPVFoam::setTime(int nRequest, const double requestTimes[])
 {
     Time& runTime = dbPtr_();
@@ -417,11 +482,6 @@ void Foam::vtkPVFoam::Update
 
     // Set up mesh parts selection(s)
     {
-        if (debug)
-        {
-            Info<< "<beg> updateMeshPartsStatus" << endl;
-        }
-
         vtkDataArraySelection* selection = reader_->GetPartSelection();
         const int n = selection->GetNumberOfArrays();
 
@@ -464,11 +524,6 @@ void Foam::vtkPVFoam::Update
                     << " : " << str << nl;
             }
         }
-
-        if (debug)
-        {
-            Info<< "<end> updateMeshPartsStatus" << endl;
-        }
     }
 
     reader_->UpdateProgress(0.15);
@@ -538,49 +593,51 @@ void Foam::vtkPVFoam::Update
     cachedVtp_.retain(nowActive);
     cachedVtu_.retain(nowActive);
 
-    // Reset (expire) dataset ids
-    partDataset_.clear();
-
-    // Convert meshes - start port0 at block=0
-    int blockNo = 0;
-
-    convertMeshVolume(output, blockNo);
-    convertMeshPatches(output, blockNo);
+    convertMeshVolume();
+    convertMeshPatches();
     reader_->UpdateProgress(0.6);
 
     if (reader_->GetIncludeZones())
     {
-        convertMeshCellZones(output, blockNo);
-        convertMeshFaceZones(output, blockNo);
-        convertMeshPointZones(output, blockNo);
+        convertMeshCellZones();
+        convertMeshFaceZones();
+        convertMeshPointZones();
         reader_->UpdateProgress(0.65);
     }
 
     if (reader_->GetIncludeSets())
     {
-        convertMeshCellSets(output, blockNo);
-        convertMeshFaceSets(output, blockNo);
-        convertMeshPointSets(output, blockNo);
+        convertMeshCellSets();
+        convertMeshFaceSets();
+        convertMeshPointSets();
         reader_->UpdateProgress(0.7);
     }
 
-    if (outputLagrangian)
-    {
-        // Lagrangian dual port - restart port1 at block=0
-        blockNo = 0;
-        convertMeshLagrangian(outputLagrangian, blockNo);
-    }
-    else
-    {
-        convertMeshLagrangian(output, blockNo);
-    }
+    convertMeshLagrangian();
 
     reader_->UpdateProgress(0.8);
 
     // Update fields
-    convertVolFields(output);
-    convertPointFields(output);
-    convertLagrangianFields(outputLagrangian ? outputLagrangian : output);
+    convertVolFields();
+    convertPointFields();
+    convertLagrangianFields();
+
+
+    // Assemble multiblock output
+    addOutputBlock(output, cachedVtu_, rangeVolume_, true); // One dataset
+    addOutputBlock(output, cachedVtp_, rangePatches_);
+    addOutputBlock(output, cachedVtu_, rangeCellZones_);
+    addOutputBlock(output, cachedVtp_, rangeFaceZones_);
+    addOutputBlock(output, cachedVtp_, rangePointZones_);
+    addOutputBlock(output, cachedVtu_, rangeCellSets_);
+    addOutputBlock(output, cachedVtp_, rangeFaceSets_);
+    addOutputBlock(output, cachedVtp_, rangePointSets_);
+    addOutputBlock
+    (
+        (outputLagrangian ? outputLagrangian : output),
+        cachedVtp_,
+        rangeLagrangian_
+    );
 
     if (debug)
     {
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.H b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.H
index 4a320affd2..0c520515ba 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.H
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoam.H
@@ -55,6 +55,10 @@ SourceFiles
 #include "volPointInterpolation.H"
 #include "foamPvCore.H"
 
+#include "vtkSmartPointer.h"
+#include "vtkPolyData.h"
+#include "vtkUnstructuredGrid.h"
+
 // * * * * * * * * * * * * * Forward Declarations  * * * * * * * * * * * * * //
 
 class vtkDataArraySelection;
@@ -65,11 +69,8 @@ class vtkPVFoamReader;
 class vtkRenderer;
 class vtkTextActor;
 class vtkMultiBlockDataSet;
-class vtkPolyData;
-class vtkUnstructuredGrid;
 class vtkIndent;
 
-template<class T> class vtkSmartPointer;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,6 +83,7 @@ class Time;
 class fvMesh;
 class IOobjectList;
 class polyPatch;
+class fvMeshSubset;
 
 template<class Type> class IOField;
 template<class Type> class Field;
@@ -160,16 +162,19 @@ class vtkPVFoam
             }
         };
 
+
         //- Bookkeeping for vtkPolyData
         class foamVtpData
         {
+        public:
+            vtkSmartPointer<vtkPolyData> vtkmesh;
         };
 
         //- Bookkeeping for vtkUnstructuredGrid
-        class foamVtuData
-        :
-            public polyDecomp
+        class foamVtuData : public polyDecomp
         {
+        public:
+            vtkSmartPointer<vtkUnstructuredGrid> vtkmesh;
         };
 
 
@@ -202,9 +207,6 @@ class vtkPVFoam
         //- The index of selected parts mapped to their names
         Map<string> selectedPartIds_;
 
-        //- Datasets corresponding to selected geometrical pieces
-        Map<label> partDataset_;
-
         //- Any information for 2D (VTP) geometries
         HashTable<foamVtpData, string> cachedVtp_;
 
@@ -213,7 +215,7 @@ class vtkPVFoam
 
         //- First instance and size of various mesh parts
         //  used to index into selectedPartIds and thus indirectly into
-        //  cachedVtp, cachedVtu and partDataset
+        //  cachedVtp, cachedVtu
         arrayRange rangeVolume_;
         arrayRange rangePatches_;
         arrayRange rangeLagrangian_;
@@ -230,6 +232,16 @@ class vtkPVFoam
 
     // Private Member Functions
 
+        template<class Container>
+        bool addOutputBlock
+        (
+            vtkMultiBlockDataSet* output,
+            const HashTable<Container, string>& cache,
+            const arrayRange& selector,
+            const bool singleDataset = false
+        ) const;
+
+
         //- Reset data counters
         void resetCounters();
 
@@ -279,31 +291,38 @@ class vtkPVFoam
       // Mesh conversion functions
 
         //- Convert volume mesh
-        void convertMeshVolume(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshVolume();
 
         //- Convert Lagrangian points
-        void convertMeshLagrangian(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshLagrangian();
 
         //- Convert mesh patches
-        void convertMeshPatches(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshPatches();
+
+        //- Convert subsetted mesh
+        void convertMeshSubset
+        (
+            const fvMeshSubset& subsetter,
+            const string& longName
+        );
 
         //- Convert cell zones
-        void convertMeshCellZones(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshCellZones();
 
         //- Convert face zones
-        void convertMeshFaceZones(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshFaceZones();
 
         //- Convert point zones
-        void convertMeshPointZones(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshPointZones();
 
         //- Convert cell sets
-        void convertMeshCellSets(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshCellSets();
 
         //- Convert face sets
-        void convertMeshFaceSets(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshFaceSets();
 
         //- Convert point sets
-        void convertMeshPointSets(vtkMultiBlockDataSet* output, int& blockNo);
+        void convertMeshPointSets();
 
 
       // Add mesh functions
@@ -359,13 +378,13 @@ class vtkPVFoam
 
 
         //- Convert volume fields
-        void convertVolFields(vtkMultiBlockDataSet* output);
+        void convertVolFields();
 
         //- Convert point fields
-        void convertPointFields(vtkMultiBlockDataSet* output);
+        void convertPointFields();
 
         //- Convert Lagrangian fields
-        void convertLagrangianFields(vtkMultiBlockDataSet* output);
+        void convertLagrangianFields();
 
 
       // Convert OpenFOAM fields
@@ -375,8 +394,7 @@ class vtkPVFoam
         void convertVolField
         (
             const PtrList<patchInterpolator>& patchInterpList,
-            const GeometricField<Type, fvPatchField, volMesh>& fld,
-            vtkMultiBlockDataSet* output
+            const GeometricField<Type, fvPatchField, volMesh>& fld
         );
 
         //- Volume fields - all types
@@ -385,8 +403,7 @@ class vtkPVFoam
         (
             const fvMesh& mesh,
             const PtrList<patchInterpolator>& patchInterpList,
-            const IOobjectList& objects,
-            vtkMultiBlockDataSet* output
+            const IOobjectList& objects
         );
 
         //- Volume internal fields (DimensionedField)- all types
@@ -395,8 +412,7 @@ class vtkPVFoam
         (
             const fvMesh& mesh,
             const PtrList<patchInterpolator>& patchInterpList,
-            const IOobjectList& objects,
-            vtkMultiBlockDataSet* output
+            const IOobjectList& objects
         );
 
         //- Volume field - all selected parts
@@ -405,7 +421,6 @@ class vtkPVFoam
         (
             const GeometricField<Type, fvPatchField, volMesh>& fld,
             autoPtr<GeometricField<Type, pointPatchField, pointMesh>>& ptfPtr,
-            vtkMultiBlockDataSet* output,
             const arrayRange& range
         );
 
@@ -414,8 +429,7 @@ class vtkPVFoam
         void convertLagrangianFields
         (
             const IOobjectList& objects,
-            vtkMultiBlockDataSet* output,
-            const label datasetNo
+            vtkPolyData* vtkmesh
         );
 
         //- Point fields - all types
@@ -423,8 +437,7 @@ class vtkPVFoam
         void convertPointFields
         (
             const pointMesh& pMesh,
-            const IOobjectList& objects,
-            vtkMultiBlockDataSet* output
+            const IOobjectList& objectst
         );
 
         //- Point field - all selected parts
@@ -432,7 +445,6 @@ class vtkPVFoam
         void convertPointFieldBlock
         (
             const GeometricField<Type, pointPatchField, pointMesh>& pfld,
-            vtkMultiBlockDataSet* output,
             const arrayRange& range
         );
 
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFieldTemplates.C b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFieldTemplates.C
index 70d0bdabe4..3df2eb9925 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFieldTemplates.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFieldTemplates.C
@@ -53,8 +53,7 @@ template<class Type>
 void Foam::vtkPVFoam::convertVolField
 (
     const PtrList<patchInterpolator>& patchInterpList,
-    const GeometricField<Type, fvPatchField, volMesh>& fld,
-    vtkMultiBlockDataSet* output
+    const GeometricField<Type, fvPatchField, volMesh>& fld
 )
 {
     const fvMesh& mesh = fld.mesh();
@@ -78,38 +77,12 @@ void Foam::vtkPVFoam::convertVolField
         );
     }
 
-    // Convert activated internalMesh regions
-    convertVolFieldBlock
-    (
-        fld,
-        ptfPtr,
-        output,
-        rangeVolume_
-    );
-
-    // Convert activated cellZones
-    convertVolFieldBlock
-    (
-        fld,
-        ptfPtr,
-        output,
-        rangeCellZones_
-    );
-
-    // Convert activated cellSets
-    convertVolFieldBlock
-    (
-        fld,
-        ptfPtr,
-        output,
-        rangeCellSets_
-    );
+    convertVolFieldBlock(fld, ptfPtr, rangeVolume_);    // internalMesh
+    convertVolFieldBlock(fld, ptfPtr, rangeCellZones_); // cellZones
+    convertVolFieldBlock(fld, ptfPtr, rangeCellSets_);  // cellSets
 
 
-    //
-    // Convert patches - if activated
-    // - skip field conversion for groups
-    //
+    // Patches - skip field conversion for groups
     for (auto partId : rangePatches_)
     {
         if
@@ -122,7 +95,6 @@ void Foam::vtkPVFoam::convertVolField
         }
 
         const word patchName = getPartName(partId);
-        const label datasetNo = partDataset_.lookup(partId, -1);
         const label patchId = patches.findPatchID(patchName);
 
         if (patchId < 0)
@@ -130,15 +102,12 @@ void Foam::vtkPVFoam::convertVolField
             continue;
         }
 
-        vtkPolyData* vtkmesh = getDataFromBlock<vtkPolyData>
-        (
-            output, rangePatches_, datasetNo
-        );
-
-        if (!vtkmesh)
+        auto iter = cachedVtp_.find(selectedPartIds_[partId]);
+        if (!iter.found() || !(*iter).vtkmesh)
         {
             continue;
         }
+        auto vtkmesh = (*iter).vtkmesh;
 
         const fvPatchField<Type>& ptf = fld.boundaryField()[patchId];
 
@@ -201,9 +170,8 @@ void Foam::vtkPVFoam::convertVolField
         }
     }
 
-    //
-    // Convert face zones - if activated
-    //
+
+    // Face Zones
     for (auto partId : rangeFaceZones_)
     {
         if (!selectedPartIds_.found(partId))
@@ -211,14 +179,14 @@ void Foam::vtkPVFoam::convertVolField
             continue;
         }
 
-        const word zoneName = getPartName(partId);
-        const label datasetNo = partDataset_.lookup(partId, -1);
-
-        if (datasetNo < 0)
+        auto iter = cachedVtu_.find(selectedPartIds_[partId]);
+        if (!iter.found() || !(*iter).vtkmesh)
         {
             continue;
         }
+        auto vtkmesh = (*iter).vtkmesh;
 
+        const word zoneName = getPartName(partId);
         const faceZoneMesh& zMesh = mesh.faceZones();
         const label zoneId = zMesh.findZoneID(zoneName);
 
@@ -227,28 +195,19 @@ void Foam::vtkPVFoam::convertVolField
             continue;
         }
 
-        vtkPolyData* vtkmesh = getDataFromBlock<vtkPolyData>
+        vtkSmartPointer<vtkFloatArray> cdata = convertFaceFieldToVTK
         (
-            output, rangeFaceZones_, datasetNo
+            fld,
+            zMesh[zoneId]
         );
 
-        if (vtkmesh)
-        {
-            vtkSmartPointer<vtkFloatArray> cdata = convertFaceFieldToVTK
-            (
-                fld,
-                zMesh[zoneId]
-            );
-
-            vtkmesh->GetCellData()->AddArray(cdata);
-        }
+        vtkmesh->GetCellData()->AddArray(cdata);
 
         // TODO: points
     }
 
-    //
-    // Convert face sets - if activated
-    //
+
+    // Face Sets
     for (auto partId : rangeFaceSets_)
     {
         if (!selectedPartIds_.found(partId))
@@ -256,19 +215,14 @@ void Foam::vtkPVFoam::convertVolField
             continue;
         }
 
-        const word selectName = getPartName(partId);
-        const label datasetNo = partDataset_.lookup(partId, -1);
-
-        vtkPolyData* vtkmesh = getDataFromBlock<vtkPolyData>
-        (
-            output, rangeFaceSets_, datasetNo
-        );
-
-        if (!vtkmesh)
+        auto iter = cachedVtu_.find(selectedPartIds_[partId]);
+        if (!iter.found() || !(*iter).vtkmesh)
         {
             continue;
         }
+        auto vtkmesh = (*iter).vtkmesh;
 
+        const word selectName = getPartName(partId);
         const faceSet fSet(mesh, selectName);
 
         vtkSmartPointer<vtkFloatArray> cdata = convertFaceFieldToVTK
@@ -289,13 +243,12 @@ void Foam::vtkPVFoam::convertVolFields
 (
     const fvMesh& mesh,
     const PtrList<patchInterpolator>& patchInterpList,
-    const IOobjectList& objects,
-    vtkMultiBlockDataSet* output
+    const IOobjectList& objects
 )
 {
     forAllConstIters(objects, iter)
     {
-        // restrict to GeometricField<Type, ...>
+        // Restrict to GeometricField<Type, ...>
         if
         (
             iter()->headerClassName()
@@ -313,7 +266,7 @@ void Foam::vtkPVFoam::convertVolFields
         );
 
         // Convert
-        convertVolField(patchInterpList, fld, output);
+        convertVolField(patchInterpList, fld);
     }
 }
 
@@ -323,15 +276,14 @@ void Foam::vtkPVFoam::convertDimFields
 (
     const fvMesh& mesh,
     const PtrList<patchInterpolator>& patchInterpList,
-    const IOobjectList& objects,
-    vtkMultiBlockDataSet* output
+    const IOobjectList& objects
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
 
     forAllConstIters(objects, iter)
     {
-        // restrict to DimensionedField<Type, ...>
+        // Restrict to DimensionedField<Type, ...>
         if
         (
             iter()->headerClassName()
@@ -374,7 +326,7 @@ void Foam::vtkPVFoam::convertDimFields
         );
         volFld.correctBoundaryConditions();
 
-        convertVolField(patchInterpList, volFld, output);
+        convertVolField(patchInterpList, volFld);
     }
 }
 
@@ -384,7 +336,6 @@ void Foam::vtkPVFoam::convertVolFieldBlock
 (
     const GeometricField<Type, fvPatchField, volMesh>& fld,
     autoPtr<GeometricField<Type, pointPatchField, pointMesh>>& ptfPtr,
-    vtkMultiBlockDataSet* output,
     const arrayRange& range
 )
 {
@@ -395,47 +346,39 @@ void Foam::vtkPVFoam::convertVolFieldBlock
             continue;
         }
 
-        const auto& longName = selectedPartIds_[partId];
-        const label datasetNo = partDataset_.lookup(partId, -1);
-
-        vtkUnstructuredGrid* vtkmesh =
-            getDataFromBlock<vtkUnstructuredGrid>(output, range, datasetNo);
-
-        if (!vtkmesh || !cachedVtu_.found(longName))
+        auto iter = cachedVtu_.find(selectedPartIds_[partId]);
+        if (!iter.found() || !(*iter).vtkmesh)
         {
             continue;
         }
+        auto vtuData = (*iter);
+        auto vtkmesh = (*iter).vtkmesh;
 
         vtkSmartPointer<vtkFloatArray> cdata = convertVolFieldToVTK
         (
             fld,
-            cachedVtu_[longName]
+            vtuData
         );
         vtkmesh->GetCellData()->AddArray(cdata);
 
         if (ptfPtr.valid())
         {
-            convertPointField
-            (
-                vtkmesh, ptfPtr(), fld,
-                cachedVtu_[longName]
-            );
+            convertPointField(vtkmesh, ptfPtr(), fld, vtuData);
         }
     }
 }
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 //
 // point-fields
 //
-
 template<class Type>
 void Foam::vtkPVFoam::convertPointFields
 (
     const pointMesh& pMesh,
-    const IOobjectList& objects,
-    vtkMultiBlockDataSet* output
+    const IOobjectList& objects
 )
 {
     const polyMesh& mesh = pMesh.mesh();
@@ -444,7 +387,7 @@ void Foam::vtkPVFoam::convertPointFields
     forAllConstIters(objects, iter)
     {
         const word& fieldName = iter()->name();
-        // restrict to this GeometricField<Type, ...>
+        // Restrict to this GeometricField<Type, ...>
         if
         (
             iter()->headerClassName()
@@ -462,34 +405,12 @@ void Foam::vtkPVFoam::convertPointFields
         GeometricField<Type, pointPatchField, pointMesh> pfld(*iter(), pMesh);
 
 
-        // Convert internalMesh (if selected)
-        convertPointFieldBlock
-        (
-            pfld,
-            output,
-            rangeVolume_
-        );
-
-        // Convert (selected) cellZones
-        convertPointFieldBlock
-        (
-            pfld,
-            output,
-            rangeCellZones_
-        );
-
-        // Convert (selected) cellSets
-        convertPointFieldBlock
-        (
-            pfld,
-            output,
-            rangeCellSets_
-        );
+        convertPointFieldBlock(pfld, rangeVolume_);     // internalMesh
+        convertPointFieldBlock(pfld, rangeCellZones_);  // cellZones
+        convertPointFieldBlock(pfld, rangeCellSets_);   // cellSets
 
 
-        //
-        // Convert patches (if selected)
-        //
+        // Patches)
         for (auto partId : rangePatches_)
         {
             if (!selectedPartIds_.found(partId))
@@ -505,7 +426,6 @@ void Foam::vtkPVFoam::convertPointFields
             }
 
             const word patchName = getPartName(partId);
-            const label datasetNo = partDataset_.lookup(partId, -1);
             const label patchId = patches.findPatchID(patchName);
 
             if (patchId < 0)
@@ -513,21 +433,20 @@ void Foam::vtkPVFoam::convertPointFields
                 continue;
             }
 
-            vtkPolyData* vtkmesh = getDataFromBlock<vtkPolyData>
+            auto iter = cachedVtp_.find(selectedPartIds_[partId]);
+            if (!iter.found() || !(*iter).vtkmesh)
+            {
+                continue;
+            }
+            auto vtkmesh = (*iter).vtkmesh;
+
+            vtkSmartPointer<vtkFloatArray> pdata = convertFieldToVTK
             (
-                output, rangePatches_, datasetNo
+                fieldName,
+                pfld.boundaryField()[patchId].patchInternalField()()
             );
 
-            if (vtkmesh)
-            {
-                vtkSmartPointer<vtkFloatArray> pdata = convertFieldToVTK
-                (
-                    fieldName,
-                    pfld.boundaryField()[patchId].patchInternalField()()
-                );
-
-                vtkmesh->GetPointData()->AddArray(pdata);
-            }
+            vtkmesh->GetPointData()->AddArray(pdata);
         }
 
         //
@@ -540,8 +459,14 @@ void Foam::vtkPVFoam::convertPointFields
                 continue;
             }
 
+            auto iter = cachedVtp_.find(selectedPartIds_[partId]);
+            if (!iter.found() || !(*iter).vtkmesh)
+            {
+                continue;
+            }
+            auto vtkmesh = (*iter).vtkmesh;
+
             const word zoneName = getPartName(partId);
-            const label datasetNo = partDataset_.lookup(partId, -1);
             const label zoneId = mesh.faceZones().findZoneID(zoneName);
 
             if (zoneId < 0)
@@ -549,29 +474,21 @@ void Foam::vtkPVFoam::convertPointFields
                 continue;
             }
 
-            vtkPolyData* vtkmesh = getDataFromBlock<vtkPolyData>
+            // Extract the field on the zone
+            Field<Type> znfld
             (
-                output, rangeFaceZones_, datasetNo
+                pfld.primitiveField(),
+                mesh.faceZones()[zoneId]().meshPoints()
             );
 
-            if (vtkmesh)
-            {
-                // Extract the field on the zone
-                Field<Type> znfld
+            vtkSmartPointer<vtkFloatArray> pdata =
+                convertFieldToVTK
                 (
-                    pfld.primitiveField(),
-                    mesh.faceZones()[zoneId]().meshPoints()
+                    fieldName,
+                    znfld
                 );
 
-                vtkSmartPointer<vtkFloatArray> pdata =
-                    convertFieldToVTK
-                    (
-                        fieldName,
-                        znfld
-                    );
-
-                vtkmesh->GetPointData()->AddArray(pdata);
-            }
+            vtkmesh->GetPointData()->AddArray(pdata);
         }
     }
 }
@@ -581,7 +498,6 @@ template<class Type>
 void Foam::vtkPVFoam::convertPointFieldBlock
 (
     const GeometricField<Type, pointPatchField, pointMesh>& pfld,
-    vtkMultiBlockDataSet* output,
     const arrayRange& range
 )
 {
@@ -592,25 +508,21 @@ void Foam::vtkPVFoam::convertPointFieldBlock
             continue;
         }
 
-        const auto& longName = selectedPartIds_[partId];
-        const label datasetNo = partDataset_.lookup(partId, -1);
-
-        vtkUnstructuredGrid* vtkmesh = getDataFromBlock<vtkUnstructuredGrid>
-        (
-            output, range, datasetNo
-        );
-
-        if (!vtkmesh || !cachedVtu_.found(longName))
+        auto iter = cachedVtu_.find(selectedPartIds_[partId]);
+        if (!iter.found() || ! (*iter).vtkmesh)
         {
-            convertPointField
-            (
-                vtkmesh,
-                pfld,
-                GeometricField<Type, fvPatchField, volMesh>::null(),
-                cachedVtu_[longName]
-            );
+            continue;
         }
+        auto vtuData = (*iter);
+        auto vtkmesh = (*iter).vtkmesh;
 
+        convertPointField
+        (
+            vtkmesh,
+            pfld,
+            GeometricField<Type, fvPatchField, volMesh>::null(),
+            vtuData
+        );
     }
 }
 
@@ -739,32 +651,24 @@ template<class Type>
 void Foam::vtkPVFoam::convertLagrangianFields
 (
     const IOobjectList& objects,
-    vtkMultiBlockDataSet* output,
-    const label datasetNo
+    vtkPolyData* vtkmesh
 )
 {
-    const arrayRange& range = rangeLagrangian_;
-
     forAllConstIters(objects, iter)
     {
         // restrict to this IOField<Type>
         if (iter()->headerClassName() == IOField<Type>::typeName)
         {
-            vtkPolyData* vtkmesh =
-                getDataFromBlock<vtkPolyData>(output, range, datasetNo);
+            IOField<Type> fld(*iter());
 
-            if (vtkmesh)
-            {
-                IOField<Type> fld(*iter());
+            vtkSmartPointer<vtkFloatArray> fldData =
+                convertFieldToVTK
+                (
+                    fld.name(),
+                    fld
+                );
 
-                vtkSmartPointer<vtkFloatArray> fldData =
-                    convertFieldToVTK
-                    (
-                        fld.name(),
-                        fld
-                    );
-                vtkmesh->GetPointData()->AddArray(fldData);
-            }
+            vtkmesh->GetPointData()->AddArray(fldData);
         }
     }
 }
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFields.C b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFields.C
index 697737d74d..8f1ca82460 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFields.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamFields.C
@@ -62,10 +62,7 @@ void Foam::vtkPVFoam::pruneObjectList
 }
 
 
-void Foam::vtkPVFoam::convertVolFields
-(
-    vtkMultiBlockDataSet* output
-)
+void Foam::vtkPVFoam::convertVolFields()
 {
     const fvMesh& mesh = *meshPtr_;
 
@@ -121,17 +118,17 @@ void Foam::vtkPVFoam::convertVolFields
         }
     }
 
-    convertVolFields<scalar>(mesh, interpLst, objects, output);
-    convertVolFields<vector>(mesh, interpLst, objects, output);
-    convertVolFields<sphericalTensor>(mesh, interpLst, objects, output);
-    convertVolFields<symmTensor>(mesh, interpLst, objects, output);
-    convertVolFields<tensor>(mesh, interpLst, objects, output);
+    convertVolFields<scalar>(mesh, interpLst, objects);
+    convertVolFields<vector>(mesh, interpLst, objects);
+    convertVolFields<sphericalTensor>(mesh, interpLst, objects);
+    convertVolFields<symmTensor>(mesh, interpLst, objects);
+    convertVolFields<tensor>(mesh, interpLst, objects);
 
-    convertDimFields<scalar>(mesh, interpLst, objects, output);
-    convertDimFields<vector>(mesh, interpLst, objects, output);
-    convertDimFields<sphericalTensor>(mesh, interpLst, objects, output);
-    convertDimFields<symmTensor>(mesh, interpLst, objects, output);
-    convertDimFields<tensor>(mesh, interpLst, objects, output);
+    convertDimFields<scalar>(mesh, interpLst, objects);
+    convertDimFields<vector>(mesh, interpLst, objects);
+    convertDimFields<sphericalTensor>(mesh, interpLst, objects);
+    convertDimFields<symmTensor>(mesh, interpLst, objects);
+    convertDimFields<tensor>(mesh, interpLst, objects);
 
     if (debug)
     {
@@ -141,10 +138,7 @@ void Foam::vtkPVFoam::convertVolFields
 }
 
 
-void Foam::vtkPVFoam::convertPointFields
-(
-    vtkMultiBlockDataSet* output
-)
+void Foam::vtkPVFoam::convertPointFields()
 {
     const fvMesh& mesh = *meshPtr_;
 
@@ -186,11 +180,11 @@ void Foam::vtkPVFoam::convertPointFields
     // Construct interpolation on the raw mesh
     const pointMesh& pMesh = pointMesh::New(mesh);
 
-    convertPointFields<scalar>(pMesh, objects, output);
-    convertPointFields<vector>(pMesh, objects, output);
-    convertPointFields<sphericalTensor>(pMesh, objects, output);
-    convertPointFields<symmTensor>(pMesh, objects, output);
-    convertPointFields<tensor>(pMesh, objects, output);
+    convertPointFields<scalar>(pMesh, objects);
+    convertPointFields<vector>(pMesh, objects);
+    convertPointFields<sphericalTensor>(pMesh, objects);
+    convertPointFields<symmTensor>(pMesh, objects);
+    convertPointFields<tensor>(pMesh, objects);
 
     if (debug)
     {
@@ -200,10 +194,7 @@ void Foam::vtkPVFoam::convertPointFields
 }
 
 
-void Foam::vtkPVFoam::convertLagrangianFields
-(
-    vtkMultiBlockDataSet* output
-)
+void Foam::vtkPVFoam::convertLagrangianFields()
 {
     arrayRange& range = rangeLagrangian_;
     const fvMesh& mesh = *meshPtr_;
@@ -231,13 +222,14 @@ void Foam::vtkPVFoam::convertLagrangianFields
             continue;
         }
 
-        const word  cloudName = getPartName(partId);
-        const label datasetNo = partDataset_.lookup(partId, -1);
+        const word cloudName = getPartName(partId);
+        auto iter = cachedVtp_.find(selectedPartIds_[partId]);
 
-        if (datasetNo < 0)
+        if (!iter.found() || !(*iter).vtkmesh)
         {
             continue;
         }
+        auto vtkmesh = (*iter).vtkmesh;
 
         // Get the Lagrangian fields for this time and this cloud
         // but only keep selected fields
@@ -265,12 +257,12 @@ void Foam::vtkPVFoam::convertLagrangianFields
             }
         }
 
-        convertLagrangianFields<label>(objects, output, datasetNo);
-        convertLagrangianFields<scalar>(objects, output, datasetNo);
-        convertLagrangianFields<vector>(objects, output, datasetNo);
-        convertLagrangianFields<sphericalTensor>(objects, output, datasetNo);
-        convertLagrangianFields<symmTensor>(objects, output, datasetNo);
-        convertLagrangianFields<tensor>(objects, output, datasetNo);
+        convertLagrangianFields<label>(objects, vtkmesh);
+        convertLagrangianFields<scalar>(objects, vtkmesh);
+        convertLagrangianFields<vector>(objects, vtkmesh);
+        convertLagrangianFields<sphericalTensor>(objects, vtkmesh);
+        convertLagrangianFields<symmTensor>(objects, vtkmesh);
+        convertLagrangianFields<tensor>(objects, vtkmesh);
     }
 
     if (debug)
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMesh.C b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMesh.C
index 9e202a352e..e133417912 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMesh.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVFoamReader/vtkPVFoam/vtkPVFoamMesh.C
@@ -42,20 +42,14 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::vtkPVFoam::convertMeshVolume
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshVolume()
 {
     arrayRange& range = rangeVolume_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshVolume" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -63,106 +57,75 @@ void Foam::vtkPVFoam::convertMeshVolume
     // this looks like more than one part, but it isn't
     for (auto partId : range)
     {
-        if (!selectedPartIds_.found(partId))
+        if (selectedPartIds_.found(partId))
         {
-            continue;
+            const auto& longName = selectedPartIds_[partId];
+
+            vtkSmartPointer<vtkUnstructuredGrid> vtkmesh =
+                volumeVTKMesh
+                (
+                    mesh,
+                    cachedVtu_(longName)
+                );
+
+            cachedVtu_[longName].vtkmesh = vtkmesh;
         }
-
-        const auto& longName = selectedPartIds_[partId];
-        const word partName = getPartName(partId);
-
-        vtkSmartPointer<vtkUnstructuredGrid> vtkmesh = volumeVTKMesh
-        (
-            mesh,
-            cachedVtu_(longName)
-        );
-
-        if (vtkmesh)
-        {
-            addToBlock(output, vtkmesh, range, datasetNo, partName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshVolume" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshLagrangian
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshLagrangian()
 {
     arrayRange& range = rangeLagrangian_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshLagrangian" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
     for (auto partId : range)
     {
-        if (!selectedPartIds_.found(partId))
+        if (selectedPartIds_.found(partId))
         {
-            continue;
+            const auto& longName = selectedPartIds_[partId];
+            const word cloudName = getPartName(partId);
+
+            vtkSmartPointer<vtkPolyData> vtkmesh =
+                lagrangianVTKMesh
+                (
+                    mesh,
+                    cloudName
+                );
+
+            cachedVtp_(longName).vtkmesh = vtkmesh;
         }
-
-        const word cloudName = getPartName(partId);
-
-        vtkSmartPointer<vtkPolyData> vtkmesh =
-            lagrangianVTKMesh(mesh, cloudName);
-
-        if (vtkmesh)
-        {
-            addToBlock(output, vtkmesh, range, datasetNo, cloudName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshLagrangian" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshPatches
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshPatches()
 {
     arrayRange& range = rangePatches_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     if (debug)
     {
-        Info<< "<beg> convertMeshPatches" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -235,34 +198,51 @@ void Foam::vtkPVFoam::convertMeshPatches
 
         if (vtkmesh)
         {
-            addToBlock(output, vtkmesh, range, datasetNo, partName);
-            partDataset_.set(partId, datasetNo++);
+            cachedVtp_(longName).vtkmesh = vtkmesh;
         }
     }
 
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
-    }
-
     if (debug)
     {
-        Info<< "<end> convertMeshPatches" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshCellZones
+void Foam::vtkPVFoam::convertMeshSubset
 (
-    vtkMultiBlockDataSet* output,
-    int& blockNo
+    const fvMeshSubset& subsetter,
+    const string& longName
 )
 {
-    arrayRange& range = rangeCellZones_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
+    vtkSmartPointer<vtkUnstructuredGrid> vtkmesh = volumeVTKMesh
+    (
+        subsetter.subMesh(),
+        cachedVtu_(longName)
+    );
+
+    // Convert cellMap, addPointCellLabels to global cell ids
+    inplaceRenumber
+    (
+        subsetter.cellMap(),
+        cachedVtu_[longName].cellMap()
+    );
+    inplaceRenumber
+    (
+        subsetter.cellMap(),
+        cachedVtu_[longName].additionalIds()
+    );
+
+    // copy pointMap as well, otherwise pointFields fail
+    cachedVtu_[longName].pointMap() = subsetter.pointMap();
+    cachedVtu_[longName].vtkmesh = vtkmesh;
+}
+
+
+void Foam::vtkPVFoam::convertMeshCellZones()
+{
+    const arrayRange& range = rangeCellZones_;
     const fvMesh& mesh = *meshPtr_;
 
     if (range.empty())
@@ -272,7 +252,7 @@ void Foam::vtkPVFoam::convertMeshCellZones
 
     if (debug)
     {
-        Info<< "<beg> convertMeshCellZones" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -302,62 +282,25 @@ void Foam::vtkPVFoam::convertMeshCellZones
         fvMeshSubset subsetter(mesh);
         subsetter.setLargeCellSubset(zMesh[zoneId]);
 
-        vtkSmartPointer<vtkUnstructuredGrid> vtkmesh = volumeVTKMesh
-        (
-            subsetter.subMesh(),
-            cachedVtu_(longName)
-        );
-
-        if (vtkmesh)
-        {
-            // Convert cellMap, addPointCellLabels to global cell ids
-            inplaceRenumber
-            (
-                subsetter.cellMap(),
-                cachedVtu_[longName].cellMap()
-            );
-            inplaceRenumber
-            (
-                subsetter.cellMap(),
-                cachedVtu_[longName].additionalIds()
-            );
-
-            // copy pointMap as well, otherwise pointFields fail
-            cachedVtu_[longName].pointMap() = subsetter.pointMap();
-
-            addToBlock(output, vtkmesh, range, datasetNo, zoneName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
+        convertMeshSubset(subsetter, longName);
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshCellZones" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshCellSets
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshCellSets()
 {
-    arrayRange& range = rangeCellSets_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
+    const arrayRange& range = rangeCellSets_;
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshCellSets" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -380,57 +323,20 @@ void Foam::vtkPVFoam::convertMeshCellSets
         fvMeshSubset subsetter(mesh);
         subsetter.setLargeCellSubset(cSet);
 
-        vtkSmartPointer<vtkUnstructuredGrid> vtkmesh = volumeVTKMesh
-        (
-            subsetter.subMesh(),
-            cachedVtu_(longName)
-        );
-
-        if (vtkmesh)
-        {
-            // Convert cellMap, addPointCellLabels to global cell ids
-            inplaceRenumber
-            (
-                subsetter.cellMap(),
-                cachedVtu_[longName].cellMap()
-            );
-            inplaceRenumber
-            (
-                subsetter.cellMap(),
-                cachedVtu_[longName].additionalIds()
-            );
-
-            // copy pointMap as well, otherwise pointFields fail
-            cachedVtu_[longName].pointMap() = subsetter.pointMap();
-
-            addToBlock(output, vtkmesh, range, datasetNo, partName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
+        convertMeshSubset(subsetter, longName);
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshCellSets" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshFaceZones
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshFaceZones()
 {
     arrayRange& range = rangeFaceZones_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (range.empty())
@@ -440,7 +346,7 @@ void Foam::vtkPVFoam::convertMeshFaceZones
 
     if (debug)
     {
-        Info<< "<beg> convertMeshFaceZones" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -452,6 +358,7 @@ void Foam::vtkPVFoam::convertMeshFaceZones
             continue;
         }
 
+        const auto& longName = selectedPartIds_[partId];
         const word zoneName = getPartName(partId);
         const label  zoneId = zMesh.findZoneID(zoneName);
 
@@ -467,43 +374,31 @@ void Foam::vtkPVFoam::convertMeshFaceZones
         }
 
         vtkSmartPointer<vtkPolyData> vtkmesh =
-            patchVTKMesh(zoneName, zMesh[zoneId]());
+            patchVTKMesh
+            (
+                zoneName,
+                zMesh[zoneId]()
+            );
 
-        if (vtkmesh)
-        {
-            addToBlock(output, vtkmesh, range, datasetNo, zoneName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
+        cachedVtp_(longName).vtkmesh = vtkmesh;
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshFaceZones" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshFaceSets
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshFaceSets()
 {
     arrayRange& range = rangeFaceSets_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshFaceSets" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -514,6 +409,7 @@ void Foam::vtkPVFoam::convertMeshFaceSets
             continue;
         }
 
+        const auto& longName = selectedPartIds_[partId];
         const word partName = getPartName(partId);
 
         if (debug)
@@ -536,43 +432,31 @@ void Foam::vtkPVFoam::convertMeshFaceSets
         }
 
         vtkSmartPointer<vtkPolyData> vtkmesh =
-            patchVTKMesh("faceSet:" + partName, p);
+            patchVTKMesh
+            (
+                "faceSet:" + partName,
+                p
+            );
 
-        if (vtkmesh)
-        {
-            addToBlock(output, vtkmesh, range, datasetNo, partName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
+        cachedVtp_(longName).vtkmesh = vtkmesh;
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshFaceSets" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshPointZones
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshPointZones()
 {
     arrayRange& range = rangePointZones_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshPointZones" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -586,6 +470,7 @@ void Foam::vtkPVFoam::convertMeshPointZones
                 continue;
             }
 
+            const auto& longName = selectedPartIds_[partId];
             const word zoneName = getPartName(partId);
             const label zoneId = zMesh.findZoneID(zoneName);
 
@@ -616,42 +501,26 @@ void Foam::vtkPVFoam::convertMeshPointZones
 
             vtkmesh->SetPoints(vtkpoints);
 
-            if (vtkmesh)
-            {
-                addToBlock(output, vtkmesh, range, datasetNo, zoneName);
-                partDataset_.set(partId, datasetNo++);
-            }
+            cachedVtp_(longName).vtkmesh = vtkmesh;
         }
     }
 
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
-    }
-
     if (debug)
     {
-        Info<< "<end> convertMeshPointZones" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
 
 
-void Foam::vtkPVFoam::convertMeshPointSets
-(
-    vtkMultiBlockDataSet* output,
-    int& blockNo
-)
+void Foam::vtkPVFoam::convertMeshPointSets()
 {
     arrayRange& range = rangePointSets_;
-    range.block(blockNo);      // set output block
-    label datasetNo = 0;       // restart at dataset 0
     const fvMesh& mesh = *meshPtr_;
 
     if (debug)
     {
-        Info<< "<beg> convertMeshPointSets" << endl;
+        Info<< "<beg> " << FUNCTION_NAME << endl;
         printMemory();
     }
 
@@ -662,6 +531,7 @@ void Foam::vtkPVFoam::convertMeshPointSets
             continue;
         }
 
+        const auto& longName = selectedPartIds_[partId];
         const word partName = getPartName(partId);
 
         if (debug)
@@ -691,22 +561,12 @@ void Foam::vtkPVFoam::convertMeshPointSets
 
         vtkmesh->SetPoints(vtkpoints);
 
-        if (vtkmesh)
-        {
-            addToBlock(output, vtkmesh, range, datasetNo, partName);
-            partDataset_.set(partId, datasetNo++);
-        }
-    }
-
-    // anything added?
-    if (datasetNo)
-    {
-        ++blockNo;
+        cachedVtp_(longName).vtkmesh = vtkmesh;
     }
 
     if (debug)
     {
-        Info<< "<end> convertMeshPointSets" << endl;
+        Info<< "<end> " << FUNCTION_NAME << endl;
         printMemory();
     }
 }
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.C b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.C
index d9c0373460..7e67d10a11 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.C
@@ -103,28 +103,6 @@ void Foam::foamPvCore::addToBlock
 }
 
 
-int Foam::foamPvCore::getNumberOfDataSets
-(
-    vtkMultiBlockDataSet* output,
-    const arrayRange& selector
-)
-{
-    const int blockNo = selector.block();
-
-    vtkMultiBlockDataSet* block = vtkMultiBlockDataSet::SafeDownCast
-    (
-        output->GetBlock(blockNo)
-    );
-
-    if (block)
-    {
-        return block->GetNumberOfBlocks();
-    }
-
-    return 0;
-}
-
-
 int Foam::foamPvCore::getSelected
 (
     boolList& status,
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.H b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.H
index 6243d5343c..e4ab28316d 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.H
+++ b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCore.H
@@ -158,23 +158,6 @@ public:
         const std::string& datasetName
     );
 
-    //- Convenience method for the VTK multiblock API
-    //  Always returns a nullptr if datasetNo is negative
-    template<class Type=vtkDataSet>
-    static Type* getDataFromBlock
-    (
-        vtkMultiBlockDataSet* output,
-        const arrayRange& selector,
-        const label datasetNo
-    );
-
-    //- Convenience method for the VTK multiblock API
-    static int getNumberOfDataSets
-    (
-        vtkMultiBlockDataSet* output,
-        const arrayRange& selector
-    );
-
     //- Add objects of Type to array selection
     template<class Type>
     static label addToSelection
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCoreTemplates.C b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCoreTemplates.C
index 7908a7752d..46c4bf3eb7 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCoreTemplates.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/foamPv/foamPvCoreTemplates.C
@@ -25,36 +25,9 @@ License
 
 #include "IOobjectList.H"
 #include "vtkDataArraySelection.h"
-#include "vtkMultiBlockDataSet.h"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type>
-Type* Foam::foamPvCore::getDataFromBlock
-(
-    vtkMultiBlockDataSet* output,
-    const arrayRange& selector,
-    const label datasetNo
-)
-{
-    const int blockNo = selector.block();
-
-    vtkMultiBlockDataSet* block =
-    (
-        datasetNo < 0
-      ? nullptr
-      : vtkMultiBlockDataSet::SafeDownCast(output->GetBlock(blockNo))
-    );
-
-    if (block)
-    {
-        return Type::SafeDownCast(block->GetBlock(datasetNo));
-    }
-
-    return nullptr;
-}
-
-
 template<class Type>
 Foam::label Foam::foamPvCore::addToSelection
 (