diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files
index aaec9f26ab05eec617f82b9c505be4af6616b175..f10482758a05b8e8333e0cea08e1c2c52d71edf4 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/files
@@ -1,4 +1,3 @@
-foamVtkLagrangianWriter.C
 foamToVTK.C
 
 EXE = $(FOAM_APPBIN)/foamToVTK
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options
index 27111f09d4689f20ba067aaad33eab7a6073299b..bae30bd5d502855b3c827547bf69a72f87f9ec0c 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options
@@ -2,14 +2,17 @@ EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/conversion/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
 
 EXE_LIBS = \
     -lconversion \
     -ldynamicMesh \
     -lfiniteArea \
-    -llagrangian \
+    -llagrangianIntermediate \
+    -lregionModels \
     -lgenericPatchFields
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertAreaFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertAreaFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..6f2a17942adbdb661c42b3f1d10f6b62dd8cef35
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertAreaFields.H
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Description
+    Code chunk for converting finite-area - included by foamToVTK.
+
+\*---------------------------------------------------------------------------*/
+
+//
+// Finite-area mesh and fields - need not exist
+//
+
+// No subsetting!
+if (doFiniteArea)
+{
+    using reportFields = foamToVtkReportFields;
+
+    autoPtr<faMesh> faMeshPtr;
+
+    const label nAreaFields =
+        objects.count(stringListOps::foundOp<word>(fieldTypes::area));
+
+    if (nAreaFields)
+    {
+        const bool throwing = FatalError.throwExceptions();
+        try
+        {
+            faMeshPtr.reset(new faMesh(meshProxy.baseMesh()));
+        }
+        catch (Foam::error& err)
+        {
+            faMeshPtr.clear();
+        }
+        FatalError.throwExceptions(throwing);
+    }
+
+    if (faMeshPtr.valid() && nAreaFields)
+    {
+        reportFields::area(Info, objects);
+
+        const auto& pp = faMeshPtr->patch();
+
+        vtk::surfaceMeshWriter writer
+        (
+            pp,
+            writeOpts,
+            (
+                outputDir/regionPrefix/"finite-area"
+              / "finiteArea" + timeDesc
+            ),
+            Pstream::parRun()
+        );
+
+        writer.beginFile(faMeshPtr->name());
+
+        writer.writeTimeValue(timeValue);
+        writer.writeGeometry();
+
+        writer.beginCellData(nAreaFields);
+
+        writeAllAreaFields
+        (
+            writer,
+            *faMeshPtr,
+            objects,
+            true // syncPar
+        );
+
+        fileName outputName(writer.output());
+
+        writer.close();
+
+        if (Pstream::master())
+        {
+            // Add to file-series and emit as JSON
+
+            fileName seriesName(vtk::seriesWriter::base(outputName));
+
+            vtk::seriesWriter& series = vtkSeries(seriesName);
+
+            // First time?
+            // Load from file, verify against filesystem,
+            // prune time >= currentTime
+            if (series.empty())
+            {
+                series.load(seriesName, true, timeValue);
+            }
+
+            series.append(timeValue, outputName);
+            series.write(seriesName);
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertLagrangian.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertLagrangian.H
new file mode 100644
index 0000000000000000000000000000000000000000..9a2596b826182c16ed77432c14d1736e49091455
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertLagrangian.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Description
+    Code chunk for post-processing conversion of cloud(s) to VTK PolyData
+    (.vtp extension).
+
+\*---------------------------------------------------------------------------*/
+
+if (doLagrangian)
+{
+    const fileName cloudPrefix = (regionPrefix/cloud::prefix);
+
+    wordList cloudNames = ListOps::create<word>
+    (
+        readDir
+        (
+            mesh.time().path()/mesh.time().timeName()/cloudPrefix,
+            fileName::DIRECTORY
+        ),
+        nameOp<fileName>()
+    );
+
+    // Synchronise cloud names
+    Pstream::combineGather(cloudNames, ListOps::uniqueEqOp<word>());
+    Pstream::combineScatter(cloudNames);
+
+    // Consistent order
+    Foam::sort(cloudNames);
+
+
+    for (const word& cloudName : cloudNames)
+    {
+        IOobjectList cloudObjs(mesh, runTime.timeName(), cloudPrefix/cloudName);
+
+        bool isCloud = false;
+        if (cloudObjs.erase("coordinates"))
+        {
+            isCloud = true;
+        }
+        if (cloudObjs.erase("positions"))
+        {
+            isCloud = true;
+        }
+
+        if (!returnReduce(isCloud, orOp<bool>()))
+        {
+            continue;
+        }
+
+        // Limited to basic IOField types
+        cloudObjs.filterClasses
+        (
+            stringListOps::foundOp<word>(fieldTypes::basic)
+        );
+
+        // Are there cloud fields (globally)?
+        if (returnReduce(cloudObjs.empty(), andOp<bool>()))
+        {
+            continue;
+        }
+
+        vtk::lagrangianWriter writer
+        (
+            meshProxy.baseMesh(),
+            cloudName,
+            writeOpts,
+            // Output name for the cloud
+            (
+                outputDir/regionPrefix/cloud::prefix
+              / cloudName/cloudName + timeDesc
+            ),
+            Pstream::parRun()
+        );
+
+        Info<< "    Lagrangian: "
+            << writer.output().relative(runTime.globalPath()) << nl;
+
+        writer.writeTimeValue(mesh.time().value());
+        writer.writeGeometry();
+
+        // Begin CellData/PointData
+        writer.beginParcelData();
+
+        writer.writeFields<label>(cloudObjs);
+        writer.writeFields<scalar>(cloudObjs);
+        writer.writeFields<vector>(cloudObjs);
+        writer.writeFields<sphericalTensor>(cloudObjs);
+        writer.writeFields<symmTensor>(cloudObjs);
+        writer.writeFields<tensor>(cloudObjs);
+
+        fileName outputName(writer.output());
+
+        writer.close();
+
+        if (Pstream::master())
+        {
+            // Add to file-series and emit as JSON
+
+            fileName seriesName(vtk::seriesWriter::base(outputName));
+
+            vtk::seriesWriter& series = vtkSeries(seriesName);
+
+            // First time?
+            // Load from file, verify against filesystem,
+            // prune time >= currentTime
+            if (series.empty())
+            {
+                series.load(seriesName, true, timeValue);
+            }
+
+            series.append(timeValue, outputName);
+            series.write(seriesName);
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertSurfaceFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertSurfaceFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..cc5e02a793579e5142f897245b1e2fc675cca0db
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertSurfaceFields.H
@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Description
+    Code chunk for post-processing surface fields to VTK PolyData
+
+\*---------------------------------------------------------------------------*/
+
+{
+    using reportFields = foamToVtkReportFields;
+
+    // Load only once when possible
+    label nSurfaceScalarField = -1;
+    label nSurfaceVectorField = -1;
+
+    PtrList<const surfaceScalarField> sScalars;
+    PtrList<const surfaceVectorField> sVectors;
+
+    // Surface Fields
+    if (doSurfaceFields)
+    {
+        if (nSurfaceScalarField == -1)
+        {
+            sScalars = readFields<surfaceScalarField>
+            (
+                meshProxy,
+                objects,
+                selectedFields
+            );
+
+            reportFields::print("    surfScalar   :", Info, sScalars);
+            nSurfaceScalarField = sScalars.size();
+        }
+        else
+        {
+            sScalars.resize(nSurfaceScalarField); // Consistent sizing
+        }
+
+        if (nSurfaceVectorField == -1)
+        {
+            sVectors = readFields<surfaceVectorField>
+            (
+                meshProxy,
+                objects,
+                selectedFields
+            );
+
+            reportFields::print("    surfVector   :", Info, sVectors);
+            nSurfaceVectorField = sVectors.size();
+        }
+        else
+        {
+            sVectors.resize(nSurfaceVectorField); // Consistent sizing
+        }
+
+        if (sScalars.size())
+        {
+            // Change scalar fields into vector fields, but leave
+            // the count of vector fields unchanged. This allows us to
+            // easily delete these synthetic fields later.
+
+            surfaceVectorField unitNorm(mesh.Sf()/mesh.magSf());
+
+            sVectors.resize(nSurfaceVectorField + nSurfaceScalarField);
+
+            label nExtra = 0;
+            for (const auto& ssf : sScalars)
+            {
+                surfaceVectorField* tsvfPtr = (ssf * unitNorm).ptr();
+                tsvfPtr->rename(ssf.name());
+                sVectors.set(nSurfaceVectorField + nExtra, tsvfPtr);
+                ++nExtra;
+            }
+        }
+
+        if (sVectors.size())
+        {
+            vtk::surfaceFieldWriter writer
+            (
+                meshProxy.mesh(),
+                writeOpts,
+                (
+                    outputDir/regionPrefix
+                  / "surface-fields"/"surfaceFields" + timeDesc
+                ),
+                Pstream::parRun()
+            );
+
+            Info<< "    Surface   : "
+                << writer.output().relative(runTime.globalPath()) << nl;
+
+
+            writer.writeTimeValue(timeValue);
+            writer.writeGeometry();
+
+            writer.beginPointData(sVectors.size());
+
+            for (const auto& fld : sVectors)
+            {
+                writer.write(fld);
+            }
+
+            fileName outputName(writer.output());
+
+            writer.close();
+
+            if (Pstream::master())
+            {
+                // Add to file-series and emit as JSON
+
+                fileName seriesName(vtk::seriesWriter::base(outputName));
+
+                vtk::seriesWriter& series = vtkSeries(seriesName);
+
+                // First time?
+                // Load from file, verify against filesystem,
+                // prune time >= currentTime
+                if (series.empty())
+                {
+                    series.load(seriesName, true, timeValue);
+                }
+
+                series.append(timeValue, outputName);
+                series.write(seriesName);
+            }
+        }
+    }
+
+
+    // Write faceZones (POLYDATA file, one for each zone)
+
+    if (doFaceZones && !mesh.faceZones().empty())
+    {
+        if (nSurfaceScalarField == -1)
+        {
+            sScalars = readFields<surfaceScalarField>
+            (
+                meshProxy,
+                objects,
+                selectedFields
+            );
+            nSurfaceScalarField = sScalars.size();
+
+            reportFields::print("    surfScalar   :", Info, sScalars);
+        }
+        else
+        {
+            sScalars.resize(nSurfaceScalarField); // Consistent sizing
+        }
+
+        if (nSurfaceVectorField == -1)
+        {
+            sVectors = readFields<surfaceVectorField>
+            (
+                meshProxy,
+                objects,
+                selectedFields
+            );
+            nSurfaceVectorField = sVectors.size();
+
+            reportFields::print("    surfVector   :", Info, sVectors);
+        }
+        else
+        {
+            sVectors.resize(nSurfaceVectorField); // Consistent sizing
+        }
+
+        for (const faceZone& fz : mesh.faceZones())
+        {
+            indirectPrimitivePatch pp
+            (
+                IndirectList<face>(mesh.faces(), fz),
+                mesh.points()
+            );
+
+            vtk::surfaceMeshWriter writer
+            (
+                pp,
+                writeOpts,
+                (
+                    outputDir/regionPrefix/fz.name()
+                  / (meshProxy.useSubMesh() ? meshProxy.name() : fz.name())
+                  + timeDesc
+                ),
+                Pstream::parRun()
+            );
+
+            Info<< "    FaceZone  : "
+                << writer.output().relative(runTime.globalPath()) << nl;
+
+
+            writer.beginFile(fz.name());
+            writer.writeTimeValue(timeValue);
+            writer.writeGeometry();
+
+            writer.beginCellData(sScalars.size() + sVectors.size());
+
+            for (const auto& fld : sScalars)
+            {
+                writer.write(fld);
+            }
+            for (const auto& fld : sVectors)
+            {
+                writer.write(fld);
+            }
+
+            fileName outputName(writer.output());
+
+            writer.close();
+
+            if (Pstream::master())
+            {
+                // Add to file-series and emit as JSON
+
+                fileName seriesName(vtk::seriesWriter::base(outputName));
+
+                vtk::seriesWriter& series = vtkSeries(seriesName);
+
+                // First time?
+                // Load from file, verify against filesystem,
+                // prune time >= currentTime
+                if (series.empty())
+                {
+                    series.load(seriesName, true, timeValue);
+                }
+
+                series.append(timeValue, outputName);
+                series.write(seriesName);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertTopoSet.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertTopoSet.H
new file mode 100644
index 0000000000000000000000000000000000000000..49111ea890cbe4ba603e9086527a92ce7900eda0
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertTopoSet.H
@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Description
+    Code chunk for converting face and point sets - included by foamToVTK.
+
+\*---------------------------------------------------------------------------*/
+
+// flag to top-level code to signal early stop.
+bool wroteTopoSet = false;
+
+// Write faceSet (as PolyData)
+if (faceSetName.size())
+{
+    // Load
+    faceSet set(mesh, faceSetName);
+
+    fileName outputName
+    (
+        outputDir/regionPrefix/"face-set"
+      / set.name()/set.name() + timeDesc
+    );
+
+    Info<< "    faceSet   : "
+        << outputName.relative(runTime.globalPath()) << nl;
+
+    vtk::writeFaceSet
+    (
+        mesh,
+        set,
+        writeOpts,
+        outputName,
+        Pstream::parRun()
+    );
+
+    wroteTopoSet = true;
+}
+
+
+// Write pointSet (as PolyData)
+if (pointSetName.size())
+{
+    // Load
+    pointSet set(mesh, pointSetName);
+
+    fileName outputName
+    (
+        outputDir/regionPrefix/"point-set"
+      / set.name()/set.name() + timeDesc
+    );
+
+    Info<< "    pointSet  : "
+        << outputName.relative(runTime.globalPath()) << nl;
+
+    vtk::writePointSet
+    (
+        mesh,
+        set,
+        writeOpts,
+        outputName,
+        Pstream::parRun()
+    );
+
+    wroteTopoSet = true;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..ad8802f6a6f88b6a83e70a4924a02f757fee5e95
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/convertVolumeFields.H
@@ -0,0 +1,462 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Description
+    Code chunk for converting volume and dimensioned fields
+    included by foamToVTK.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+{
+    using reportFields = foamToVtkReportFields;
+
+    const label nVolFields =
+    (
+        (doInternal || doBoundary)
+      ? objects.count(stringListOps::foundOp<word>(fieldTypes::volume))
+      : 0
+    );
+
+    const label nDimFields =
+    (
+        (doInternal || doBoundary)
+      ? objects.count(stringListOps::foundOp<word>(fieldTypes::internal))
+      : 0
+    );
+
+    const label nPointFields =
+    (
+        noPointValues
+      ? 0
+      : objects.count(stringListOps::foundOp<word>(fieldTypes::point))
+    );
+
+
+    reportFields::volume(Info, objects);
+    reportFields::internal(Info, objects);
+
+
+    // Setup for the vtm writer.
+    // For legacy format, the information added is simply ignored.
+
+    fileName vtmOutputBase
+    (
+        outputDir/regionPrefix/vtkName + timeDesc
+    );
+
+    // Combined internal + boundary in a vtm file
+    vtk::vtmWriter vtmWriter;
+
+    // Collect individual boundaries into a vtm file
+    vtk::vtmWriter vtmBoundaries;
+
+    // Setup the internal writer
+    autoPtr<vtk::internalWriter> internalWriter;
+
+    // Interpolator for volume and dimensioned fields
+    autoPtr<volPointInterpolation> pInterp;
+
+    if (doInternal)
+    {
+        if (!noPointValues)
+        {
+            pInterp.reset(new volPointInterpolation(mesh));
+        }
+
+        if (vtuMeshCells.empty())
+        {
+            // Use the appropriate mesh (baseMesh or subMesh)
+            vtuMeshCells.reset(meshProxy.mesh());
+        }
+
+        internalWriter = autoPtr<vtk::internalWriter>::New
+        (
+            meshProxy.mesh(),
+            vtuMeshCells,
+            writeOpts,
+            // The output base name for internal
+            (
+                writeOpts.legacy()
+              ? vtmOutputBase
+              : (vtmOutputBase / "internal")
+            ),
+            Pstream::parRun()
+        );
+
+        // No sub-block for internal
+        vtmWriter.append_vtu
+        (
+            "internal",
+            vtmOutputBase.name()/"internal"
+        );
+
+        Info<< "    Internal  : "
+            << internalWriter->output().relative(runTime.globalPath()) << nl;
+
+        internalWriter->writeTimeValue(mesh.time().value());
+        internalWriter->writeGeometry();
+    }
+
+
+    // Setup the patch writers
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    PtrList<vtk::patchWriter> patchWriters;
+    PtrList<PrimitivePatchInterpolation<primitivePatch>> patchInterps;
+
+    labelList patchIds;
+    if (doBoundary)
+    {
+        patchIds = getSelectedPatches(patches, includePatches, excludePatches);
+    }
+
+    if (oneBoundary && patchIds.size())
+    {
+        auto writer = autoPtr<vtk::patchWriter>::New
+        (
+            meshProxy.mesh(),
+            patchIds,
+            writeOpts,
+            nearCellValue,
+            // Output one patch: "boundary"
+            (
+                writeOpts.legacy()
+              ?
+                (
+                    outputDir/regionPrefix/"boundary"
+                  / (meshProxy.useSubMesh() ? meshProxy.name() : "boundary")
+                  + timeDesc
+                )
+              : (vtmOutputBase / "boundary")
+            ),
+            Pstream::parRun()
+        );
+
+        // No sub-block for one-patch
+        vtmWriter.append_vtp
+        (
+            "boundary",
+            vtmOutputBase.name()/"boundary"
+        );
+
+        Info<< "    Boundaries: "
+            << writer->output().relative(runTime.globalPath()) << nl;
+
+        writer->writeTimeValue(timeValue);
+        writer->writeGeometry();
+
+        // Transfer writer to list for later use
+        patchWriters.resize(1);
+        patchWriters.set(0, writer);
+
+        // Avoid patchInterpolation for each sub-patch
+        patchInterps.resize(1); // == nullptr
+    }
+    else if (patchIds.size())
+    {
+        patchWriters.resize(patchIds.size());
+        if (!noPointValues)
+        {
+            patchInterps.resize(patchIds.size());
+        }
+
+        label nPatchWriters = 0;
+        label nPatchInterps = 0;
+
+        for (const label patchId : patchIds)
+        {
+            const polyPatch& pp = patches[patchId];
+
+            auto writer = autoPtr<vtk::patchWriter>::New
+            (
+                meshProxy.mesh(),
+                labelList(one(), pp.index()),
+                writeOpts,
+                nearCellValue,
+                // Output patch: "boundary"/name
+                (
+                    writeOpts.legacy()
+                  ?
+                    (
+                        outputDir/regionPrefix/pp.name()
+                      / (meshProxy.useSubMesh() ? meshProxy.name() : pp.name())
+                      + timeDesc
+                    )
+                  : (vtmOutputBase / "boundary" / pp.name())
+                ),
+                Pstream::parRun()
+            );
+
+            if (!nPatchWriters)
+            {
+                vtmWriter.beginBlock("boundary");
+                vtmBoundaries.beginBlock("boundary");
+            }
+
+            vtmWriter.append_vtp
+            (
+                pp.name(),
+                vtmOutputBase.name()/"boundary"/pp.name()
+            );
+
+            vtmBoundaries.append_vtp
+            (
+                pp.name(),
+                "boundary"/pp.name()
+            );
+
+            Info<< "    Boundary  : "
+                << writer->output().relative(runTime.globalPath()) << nl;
+
+            writer->writeTimeValue(timeValue);
+            writer->writeGeometry();
+
+            // Transfer writer to list for later use
+            patchWriters.set(nPatchWriters++, writer);
+
+            if (patchInterps.size())
+            {
+                patchInterps.set
+                (
+                    nPatchInterps++,
+                    new PrimitivePatchInterpolation<primitivePatch>(pp)
+                );
+            }
+        }
+
+        if (nPatchWriters)
+        {
+            vtmWriter.endBlock("boundary");
+        }
+
+        patchWriters.resize(nPatchWriters);
+        patchInterps.resize(nPatchInterps);
+    }
+
+
+    // CellData
+    {
+        // Begin CellData
+        if (internalWriter.valid())
+        {
+            // cellIds + procIds (parallel)
+            internalWriter->beginCellData
+            (
+                1 + (internalWriter->parallel() ? 1 : 0)
+              + nVolFields + nDimFields
+            );
+            internalWriter->writeCellIDs();
+            internalWriter->writeProcIDs(); // parallel only
+        }
+
+        if (nVolFields)
+        {
+            for (vtk::patchWriter& writer : patchWriters)
+            {
+                writer.beginCellData(1 + nVolFields);
+                writer.writePatchIDs();
+            }
+        }
+
+        writeAllVolFields
+        (
+            internalWriter,
+            patchWriters,
+            meshProxy,
+            objects,
+            true  // syncPar
+        );
+
+        writeAllDimFields
+        (
+            internalWriter,
+            meshProxy,
+            objects,
+            true  // syncPar
+        );
+
+        // End CellData is implicit
+    }
+
+
+    // PointData
+    // - only construct pointMesh on request since it constructs
+    //   edge addressing
+    if (!noPointValues)
+    {
+        // Begin PointData
+        if (internalWriter.valid())
+        {
+            internalWriter->beginPointData
+            (
+                nVolFields + nDimFields + nPointFields
+            );
+        }
+
+        forAll(patchWriters, writeri)
+        {
+            const label nPatchFields =
+            (
+                (
+                    writeri < patchInterps.size() && patchInterps.set(writeri)
+                  ? nVolFields
+                  : 0
+                )
+              + nPointFields
+            );
+
+            if (nPatchFields)
+            {
+                patchWriters[writeri].beginPointData(nPatchFields);
+            }
+        }
+
+        writeAllVolFields
+        (
+            internalWriter, pInterp,
+            patchWriters,   patchInterps,
+            meshProxy,
+            objects,
+            true  // syncPar
+        );
+
+        writeAllDimFields
+        (
+            internalWriter, pInterp,
+            meshProxy,
+            objects,
+            true  // syncPar
+        );
+
+        writeAllPointFields
+        (
+            internalWriter,
+            patchWriters,
+            meshProxy,
+            objects,
+            true  // syncPar
+        );
+
+        // End PointData is implicit
+    }
+
+
+    // Finish writers
+    if (internalWriter.valid())
+    {
+        internalWriter->close();
+    }
+
+    for (vtk::patchWriter& writer : patchWriters)
+    {
+        writer.close();
+    }
+
+    pInterp.clear();
+    patchWriters.clear();
+    patchInterps.clear();
+
+
+    // Collective output
+
+    if (Pstream::master())
+    {
+        // Naming for vtm, file series etc.
+        fileName outputName(vtmOutputBase);
+
+        if (writeOpts.legacy())
+        {
+            if (doInternal)
+            {
+                // Add to file-series and emit as JSON
+
+                outputName.ext(vtk::legacy::fileExtension);
+
+                fileName seriesName(vtk::seriesWriter::base(outputName));
+
+                vtk::seriesWriter& series = vtkSeries(seriesName);
+
+                // First time?
+                // Load from file, verify against filesystem,
+                // prune time >= currentTime
+                if (series.empty())
+                {
+                    series.load(seriesName, true, timeValue);
+                }
+
+                series.append(timeValue, timeDesc);
+                series.write(seriesName);
+            }
+        }
+        else
+        {
+            if (vtmWriter.size())
+            {
+                // Emit ".vtm"
+
+                outputName.ext(vtmWriter.ext());
+
+                fileName seriesName(vtk::seriesWriter::base(outputName));
+
+                vtmWriter.setTime(timeValue);
+                vtmWriter.write(outputName);
+
+                // Add to file-series and emit as JSON
+
+                vtk::seriesWriter& series = vtkSeries(seriesName);
+
+                // First time?
+                // Load from file, verify against filesystem,
+                // prune time >= currentTime
+                if (series.empty())
+                {
+                    series.load(seriesName, true, timeValue);
+                }
+
+                series.append(timeValue, outputName);
+                series.write(seriesName);
+
+                // Add to multi-region vtm
+                vtmMultiRegion.add(regionName, regionPrefix, vtmWriter);
+            }
+
+            if (vtmBoundaries.size())
+            {
+                // Emit "boundary.vtm" with collection of boundaries
+
+                // Naming for vtm
+                fileName outputName(vtmOutputBase / "boundary");
+                outputName.ext(vtmBoundaries.ext());
+
+                vtmBoundaries.setTime(timeValue);
+                vtmBoundaries.write(outputName);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/createMeshes.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/createMeshes.H
new file mode 100644
index 0000000000000000000000000000000000000000..65897c0885dd0401f57764b3f684a6474d778c4d
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/createMeshes.H
@@ -0,0 +1,58 @@
+PtrList<fvMesh> meshes(regionNames.size());
+PtrList<fvMeshSubsetProxy> meshProxies(regionNames.size());
+PtrList<vtk::vtuCells> vtuMappings(regionNames.size());
+
+forAll(regionNames, regioni)
+{
+    const word& regionName = regionNames[regioni];
+
+    Info<< "Create mesh";
+
+    if (regionName != fvMesh::defaultRegion)
+    {
+        Info<< ' ' << regionName;
+    }
+    Info<< " for time = " << runTime.timeName() << nl << endl;
+
+    meshes.set
+    (
+        regioni,
+        new fvMesh
+        (
+            IOobject
+            (
+                regionName,
+                runTime.timeName(),
+                runTime,
+                IOobject::MUST_READ
+            )
+        )
+    );
+
+    // Mesh subsetting, or pass through
+    meshProxies.set
+    (
+        regioni,
+        new fvMeshSubsetProxy
+        (
+            meshes[regioni],
+            cellSubsetType,
+            cellSelectionName
+        )
+    );
+
+    // VTU sizing and decomposition information
+    vtuMappings.set
+    (
+        regioni,
+        new vtk::vtuCells(writeOpts, decomposePoly)
+    );
+}
+
+
+Info<< "VTK mesh topology: "
+    << timer.cpuTimeIncrement() << " s, "
+    << mem.update().size() << " kB" << endl;
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H
deleted file mode 100644
index 2b293f1711f286ae595e67a7daef3a13b60aff97..0000000000000000000000000000000000000000
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H
+++ /dev/null
@@ -1,77 +0,0 @@
-// check all time directories for the following:
-
-// Any cloud names:
-wordHashSet allCloudDirs;
-
-if (timeDirs.size() && !noLagrangian)
-{
-    const fileName& baseDir = mesh.time().path();
-    const fileName cloudPrefix(regionPrefix/cloud::prefix);
-
-    Info<< "Searching for lagrangian ... " << flush;
-
-    forAll(timeDirs, timeI)
-    {
-        const word& timeName = timeDirs[timeI].name();
-
-        // DO NOT USE -->>  runTime.setTime(timeDirs[timeI], timeI);  <<--
-        // It incurs a large overhead when done so frequently.
-
-        fileNameList cloudDirs = readDir
-        (
-            baseDir/timeName/cloudPrefix,
-            fileName::DIRECTORY
-        );
-
-        for (const word& cloudName : cloudDirs)
-        {
-            IOobjectList cloudObjs
-            (
-                mesh,
-                timeName,
-                cloudPrefix/cloudName
-            );
-
-            // Clouds require "coordinates".
-            // The "positions" are for v1706 and lower.
-            if (cloudObjs.found("coordinates") || cloudObjs.found("positions"))
-            {
-                if (allCloudDirs.insert(cloudName))
-                {
-                    Info<< nl << "    At time: " << timeName
-                        << " detected cloud directory : " << cloudName
-                        << flush;
-                }
-            }
-        }
-    }
-
-    if (allCloudDirs.empty())
-    {
-        Info<< "none detected." << endl;
-    }
-
-    if (Pstream::parRun())
-    {
-        Pstream::combineGather(allCloudDirs, HashSetOps::plusEqOp<word>());
-        Pstream::combineScatter(allCloudDirs);
-    }
-}
-
-
-// Sorted list of cloud names
-const wordList cloudNames(allCloudDirs.sortedToc());
-
-if (cloudNames.size())
-{
-    // Complete the echo information
-    Info<< "(";
-    for (const word& cloudName : cloudNames)
-    {
-        Info<< ' ' << cloudName;
-    }
-    Info<< " ) " << endl;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
index 3fdb93d05c1ac9f63e8150032be040173373a136..a587530b1b5e459c4a3ec04820474a3d85f6efbc 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
@@ -28,16 +28,16 @@ Group
     grpPostProcessingUtilities
 
 Description
-    VTK file format writer.
+    General OpenFOAM to VTK file writer.
 
+    Other bits
     - Handles volFields, pointFields, surfaceScalarField, surfaceVectorField
       fields.
     - Mesh topo changes.
     - Both ascii and binary.
     - Single time step writing.
     - Write subset only.
-    - Automatic decomposition of cells; polygons on boundary undecomposed since
-      handled by vtk.
+    - Optional decomposition of cells.
 
 Usage
     \b foamToVTK [OPTION]
@@ -46,11 +46,8 @@ Usage
       - \par -ascii
         Write VTK data in ASCII format instead of binary.
 
-      - \par -xml
-        Write VTK data in XML format instead of legacy format
-
-      - \par -mesh \<name\>
-        Use a different mesh name (instead of -region)
+      - \par -legacy
+        Write VTK data in legacy format instead of XML format
 
       - \par -fields \<fields\>
         Convert selected fields only. For example,
@@ -74,11 +71,14 @@ Usage
       - \par -nearCellValue
         Output cell value on patches instead of patch value itself
 
+      - \par -noBoundary
+        Suppress output for all boundary patches
+
       - \par -noInternal
-        Do not generate file for mesh, only for patches
+        Suppress output for internal (volume) mesh
 
       - \par -noLagrangian
-        Suppress writing lagrangian positions and fields.
+        Suppress writing Lagrangian positions and fields.
 
       - \par -noPointValues
         No pointFields
@@ -86,26 +86,21 @@ Usage
       - \par -noFaceZones
         No faceZones
 
-      - \par -noLinks
-        (in parallel) do not link processor files to master
-
-      - \par poly
-        write polyhedral cells without tet/pyramid decomposition
+      - \par -poly-decomp
+        Decompose polyhedral cells into tets/pyramids
 
       - \par -allPatches
-        Combine all patches into a single file
+        Combine all patches into a single boundary file
+
+      - \par -patch \<patchNames\>
+        Specify which patch or patches (name or regex) to convert.
 
       - \par -excludePatches \<patchNames\>
-        Specify patches (wildcards) to exclude. For example,
+        Specify which patch or patches (name or regex) not to convert.
+        For example,
         \verbatim
           -excludePatches '( inlet_1 inlet_2 "proc.*")'
         \endverbatim
-        The quoting is required to avoid shell expansions and to pass the
-        information as a single argument. The double quotes denote a regular
-        expression.
-
-      - \par -useTimeName
-        use the time index in the VTK file name instead of the time index
 
 Note
     The mesh subset is handled by fvMeshSubsetProxy. Slight inconsistency in
@@ -116,184 +111,102 @@ Note
     fvMeshSubset.interpolate function to directly interpolate the
     whole-mesh values onto the subset patch.
 
-Note
-    \par new file format:
-    no automatic timestep recognition.
-    However can have .pvd file format which refers to time simulation
-    if XML *.vtu files are available:
-
-    \verbatim
-      <?xml version="1.0"?>
-      <VTKFile type="Collection" version="0.1" byte_order="LittleEndian">
-        <Collection>
-          <DataSet timestep="50" file="pitzDaily_2.vtu"/>
-          <DataSet timestep="100" file="pitzDaily_3.vtu"/>
-          <DataSet timestep="150" file="pitzDaily_4.vtu"/>
-          <DataSet timestep="200" file="pitzDaily_5.vtu"/>
-          <DataSet timestep="250" file="pitzDaily_6.vtu"/>
-          <DataSet timestep="300" file="pitzDaily_7.vtu"/>
-          <DataSet timestep="350" file="pitzDaily_8.vtu"/>
-          <DataSet timestep="400" file="pitzDaily_9.vtu"/>
-          <DataSet timestep="450" file="pitzDaily_10.vtu"/>
-          <DataSet timestep="500" file="pitzDaily_11.vtu"/>
-        </Collection>
-      </VTKFile>
-    \endverbatim
-
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "pointMesh.H"
-#include "volPointInterpolation.H"
 #include "emptyPolyPatch.H"
-#include "PstreamCombineReduceOps.H"
-#include "HashOps.H"
-#include "labelIOField.H"
-#include "scalarIOField.H"
-#include "sphericalTensorIOField.H"
-#include "symmTensorIOField.H"
-#include "tensorIOField.H"
+#include "volPointInterpolation.H"
 #include "faceZoneMesh.H"
-#include "Cloud.H"
-#include "passiveParticle.H"
-#include "wordRes.H"
 #include "areaFields.H"
 #include "fvMeshSubsetProxy.H"
-#include "readFields.H"
 #include "faceSet.H"
 #include "pointSet.H"
+#include "HashOps.H"
+#include "regionProperties.H"
+#include "stringListOps.H"
 
-#include "foamVtkOutputOptions.H"
+#include "Cloud.H"
+#include "readFields.H"
+#include "reportFields.H"
+
+#include "foamVtmWriter.H"
 #include "foamVtkInternalWriter.H"
 #include "foamVtkPatchWriter.H"
 #include "foamVtkSurfaceMeshWriter.H"
 #include "foamVtkLagrangianWriter.H"
+#include "foamVtkSurfaceFieldWriter.H"
 #include "foamVtkWriteTopoSet.H"
-#include "foamVtkWriteSurfFields.H"
+#include "foamVtkSeriesWriter.H"
+
+#include "writeAreaFields.H"
+#include "writeDimFields.H"
+#include "writeVolFields.H"
+#include "writePointFields.H"
 
 #include "memInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class GeoField>
-void print(const char* msg, Ostream& os, const UPtrList<const GeoField>& flds)
-{
-    if (flds.size())
-    {
-        os  << msg;
-        forAll(flds, i)
-        {
-            os  << ' ' << flds[i].name();
-        }
-        os  << endl;
-    }
-}
-
-
-void print(const char* msg, Ostream& os, const wordList& flds)
-{
-    if (flds.size())
-    {
-        os  << msg;
-        forAll(flds, i)
-        {
-            os  << ' ' << flds[i];
-        }
-        os  << endl;
-    }
-}
-
-
 labelList getSelectedPatches
 (
     const polyBoundaryMesh& patches,
-    const wordRes& excludePatches
+    const wordRes& whitelist,
+    const wordRes& blacklist
 )
 {
     DynamicList<label> patchIDs(patches.size());
 
-    Info<< "Combining patches:" << endl;
-
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const polyPatch& pp = patches[patchi];
-
-        if
-        (
-            isType<emptyPolyPatch>(pp)
-         || (Pstream::parRun() && isType<processorPolyPatch>(pp))
-        )
-        {
-            Info<< "    discarding empty/processor patch " << patchi
-                << " " << pp.name() << endl;
-        }
-        else if (excludePatches.match(pp.name()))
+        if (isType<emptyPolyPatch>(pp))
         {
-            Info<< "    excluding patch " << patchi
-                << " " << pp.name() << endl;
+            continue;
         }
-        else
+        else if (Pstream::parRun() && isType<processorPolyPatch>(pp))
         {
-            patchIDs.append(patchi);
-            Info<< "    patch " << patchi << " " << pp.name() << endl;
+            break; // No processor patches for parallel output
         }
-    }
 
-    return patchIDs.shrink();
-}
-
-
-HashTable<wordHashSet> candidateObjects
-(
-    const IOobjectList& objects,
-    const wordHashSet& supportedTypes,
-    const bool specifiedFields,
-    const wordHashSet& selectedFields
-)
-{
-    // Special case = no fields
-    if (specifiedFields && selectedFields.empty())
-    {
-        return HashTable<wordHashSet>();
-    }
+        const word& patchName = pp.name();
 
-    HashTable<wordHashSet> usable = objects.classes();
+        bool accept = false;
 
-    // Limited to types that we explicitly handle
-    usable.retain(supportedTypes);
+        if (whitelist.size())
+        {
+            const auto matched = whitelist.matched(patchName);
 
-    // If specified, further limit to selected fields
-    if (specifiedFields)
-    {
-        forAllIters(usable, iter)
+            accept =
+            (
+                matched == wordRe::LITERAL
+              ? true
+              : (matched == wordRe::REGEX && !blacklist.match(patchName))
+            );
+        }
+        else
         {
-            iter.object().retain(selectedFields);
+            accept = !blacklist.match(patchName);
         }
 
-        // Prune entries without any fields
-        usable.filterValues
-        (
-            [](const wordHashSet& vals){ return !vals.empty(); }
-        );
+        if (accept)
+        {
+            patchIDs.append(pp.index());
+        }
     }
 
-    return usable;
+    return patchIDs.shrink();
 }
 
 
 //
-// Process args for output options
-// Default from foamVtkOutputOptions is inline ASCII xml
+// Process args for output options (-ascii, -legacy)
 //
 vtk::outputOptions getOutputOptions(const argList& args)
 {
+    // Default is inline ASCII xml
     vtk::outputOptions opts;
 
-    if (args.found("xml"))
-    {
-        opts.ascii(args.found("ascii"));
-    }
-    else
+    if (args.found("legacy"))
     {
         opts.legacy(true);
 
@@ -314,6 +227,10 @@ vtk::outputOptions getOutputOptions(const argList& args)
             }
         }
     }
+    else
+    {
+        opts.ascii(args.found("ascii"));
+    }
 
     return opts;
 }
@@ -323,147 +240,202 @@ vtk::outputOptions getOutputOptions(const argList& args)
 
 int main(int argc, char *argv[])
 {
-    argList::addNote("VTK file format writer");
+    argList::addNote("OpenFOAM to VTK writer");
     timeSelector::addOptions();
 
-    #include "addRegionOption.H"
+    // Infrequently needed, mark as advanced.
+    argList::setAdvanced("noFunctionObjects");
 
-    argList::addOption
+    argList::addBoolOption
     (
-        "fields",
-        "wordList",
-        "Only convert the specified fields - eg '(p T U)'"
+        "ascii",
+        "Write in ASCII format instead of binary"
+    );
+    argList::addBoolOption
+    (
+        "legacy",
+        "Write legacy format instead of xml"
+    );
+    argList::addBoolOption
+    (
+        "poly-decomp",
+        "Decompose polyhedral cells into tets/pyramids",
+        true  // mark as an advanced option
+    );
+    argList::ignoreOptionCompat
+    (
+        {"xml", 1806},  // xml is now default, use -legacy to toggle
+        false           // bool option, no argument
+    );
+    argList::ignoreOptionCompat
+    (
+        {"poly", 1806}, // poly is now default, use -poly-decomp to toggle
+        false           // bool option, no argument
     );
+
     argList::addOption
     (
         "cellSet",
         "name",
-        "Convert a mesh subset corresponding to the specified cellSet"
+        "Convert mesh subset corresponding to specified cellSet",
+        true  // mark as an advanced option
     );
     argList::addOption
     (
         "cellZone",
         "name",
-        "Convert a mesh subset corresponding to the specified cellZone"
+        "Convert mesh subset corresponding to specified cellZone",
+        true  // mark as an advanced option
     );
     argList::addOption
     (
         "faceSet",
         "name",
-        "Restrict conversion to the specified faceSet"
+        "Convert specified faceSet only",
+        true  // mark as an advanced option
     );
     argList::addOption
     (
         "pointSet",
         "name",
-        "Restrict conversion to the specified pointSet"
+        "Convert specified pointSet only",
+        true  // mark as an advanced option
     );
-    argList::addBoolOption
+
+    argList::addOption
     (
-        "ascii",
-        "Write in ASCII format instead of binary"
+        "fields",
+        "wordRes",
+        "Only convert specified fields\n"
+        "Eg, '(p T U \"alpha.*\")'"
     );
+
     argList::addBoolOption
     (
-        "xml",
-        "Write VTK xml instead of legacy format"
+        "surfaceFields",
+        "Write surfaceScalarFields (eg, phi)",
+        true  // mark as an advanced option
     );
     argList::addBoolOption
     (
-        "poly",
-        "Write polyhedral cells without tet/pyramid decomposition"
+        "finiteAreaFields",
+        "Write finite area fields",
+        true  // mark as an advanced option
     );
     argList::addBoolOption
     (
-        "surfaceFields",
-        "Write surfaceScalarFields (e.g., phi)"
+        "nearCellValue",
+        "Use cell value on patches instead of patch value itself",
+        true  // mark as an advanced option
     );
     argList::addBoolOption
     (
-        "finiteAreaFields",
-        "Write finite area fields"
+        "noBoundary",  // no-boundary
+        "Suppress output for boundary patches"
     );
     argList::addBoolOption
     (
-        "nearCellValue",
-        "Use cell value on patches instead of patch value itself"
+        "noInternal",  // no-internal
+        "Suppress output for internal volume mesh"
     );
     argList::addBoolOption
     (
-        "noInternal",
-        "Do not generate file for mesh, only for patches"
+        "noLagrangian",  // no-lagrangian
+        "Suppress writing lagrangian positions and fields"
     );
     argList::addBoolOption
     (
-        "noLagrangian",
-        "Suppress writing lagrangian positions and fields"
+        "noPointValues",  // no-point-data
+        "No pointFields and no interpolated PointData"
     );
     argList::addBoolOption
     (
-        "noPointValues",
-        "No pointFields"
+        "allPatches",  // one-boundary
+        "Combine all patches into a single file",
+        true  // mark as an advanced option
+    );
+
+    #include "addRegionOption.H"
+
+    argList::addOption
+    (
+        "regions",
+        "wordRes",
+        "Operate on selected regions from regionProperties.\n"
+        "Eg, '( gas \"solid.*\" )'"
     );
     argList::addBoolOption
     (
-        "allPatches",
-        "Combine all patches into a single file"
+        "allRegions",
+        "Operate on all regions in regionProperties"
+    );
+
+    argList::addOption
+    (
+        "patches",
+        "wordRes",
+        "A list of patches to include.\n"
+        "Eg, '( front \".*back\" )'"
     );
     argList::addOption
     (
         "excludePatches",
         "wordRes",
-        "A list of patches to exclude - eg '( inlet \".*Wall\" )' "
+        "A list of patches to exclude\n"
+        "Eg, '( inlet \".*Wall\" )'"
     );
+
     argList::addBoolOption
     (
         "noFaceZones",
-        "No faceZones"
+        "No faceZones",
+        true  // mark as an advanced option
     );
-    argList::addBoolOption
+    argList::ignoreOptionCompat
+    (
+        {"noLinks", 1806},      // ignore never make any links
+        false                   // bool option, no argument
+    );
+    argList::ignoreOptionCompat
     (
-        "noLinks",
-        "Do not link processor VTK files - parallel only"
+        {"useTimeName", 1806},  // ignore - works poorly with VTM formats
+        false                   // bool option, no argument
     );
     argList::addBoolOption
     (
-        "useTimeName",
-        "Use time name instead of the time index when naming files"
+        "overwrite",
+        "Remove any existing VTK output directory"
     );
     argList::addOption
     (
         "name",
         "subdir",
-        "Sub-directory name for VTK output (default: 'VTK')"
+        "Directory name for VTK output (default: 'VTK')"
     );
 
     #include "setRootCase.H"
 
-    cpuTime timer;
-    memInfo mem;
-    Info<< "Initial memory "
-        << mem.update().size() << " kB" << endl;
+    const bool decomposePoly = args.found("poly-decomp");
+    const bool doBoundary    = !args.found("noBoundary");
+    const bool doInternal    = !args.found("noInternal");
+    const bool doLagrangian  = !args.found("noLagrangian");
+    const bool doFiniteArea  = args.found("finiteAreaFields");
+    const bool doSurfaceFields = args.found("surfaceFields");
 
-    #include "createTime.H"
+    const bool doFaceZones   = !args.found("noFaceZones") && doInternal;
+    const bool oneBoundary   = args.found("allPatches") && doBoundary;
+    const bool nearCellValue = args.found("nearCellValue") && doBoundary;
+    const bool allRegions    = args.found("allRegions");
 
-    const bool decomposePoly   = !args.found("poly");
-    const bool doWriteInternal = !args.found("noInternal");
-    const bool doFaceZones     = !args.found("noFaceZones");
-    const bool doLinks         = !args.found("noLinks");
-    const bool useTimeName     = args.found("useTimeName");
-    const bool noLagrangian    = args.found("noLagrangian");
-    const bool nearCellValue   = args.found("nearCellValue");
-
-    const vtk::outputOptions fmtType = getOutputOptions(args);
+    const vtk::outputOptions writeOpts = getOutputOptions(args);
 
     if (nearCellValue)
     {
-        WarningInFunction
-            << "Using neighbouring cell value instead of patch value"
+        Info<< "Using neighbouring cell value instead of patch value"
             << nl << endl;
     }
 
     const bool noPointValues = args.found("noPointValues");
-
     if (noPointValues)
     {
         Info<< "Outputting cell values only."
@@ -471,1114 +443,312 @@ int main(int argc, char *argv[])
             << nl;
     }
 
-    const bool allPatches = args.found("allPatches");
-
-    wordRes excludePatches;
-    if (args.readListIfPresent<wordRe>("excludePatches", excludePatches))
+    wordRes includePatches, excludePatches;
+    if (doBoundary)
     {
-        Info<< "Not including patches " << excludePatches << nl << endl;
-    }
-
-    string vtkName = runTime.caseName();
-
-    fvMeshSubsetProxy::subsetType cellSubsetType = fvMeshSubsetProxy::NONE;
-    word cellSubsetName;
-    if (args.readIfPresent("cellSet", cellSubsetName))
-    {
-        vtkName = cellSubsetName;
-        cellSubsetType = fvMeshSubsetProxy::SET;
-    }
-    else if (args.readIfPresent("cellZone", cellSubsetName))
-    {
-        vtkName = cellSubsetName;
-        cellSubsetType = fvMeshSubsetProxy::ZONE;
-    }
-    else if (Pstream::parRun())
-    {
-        // Strip off leading casename, leaving just processor_DDD ending.
-        const auto i = vtkName.rfind("processor");
-
-        if (i != string::npos)
+        if (args.readListIfPresent<wordRe>("patches", includePatches))
         {
-            vtkName = vtkName.substr(i);
+            Info<< "Including patches " << flatOutput(includePatches)
+                << nl << endl;
         }
-    }
-
-    word faceSetName;
-    args.readIfPresent("faceSet", faceSetName);
-
-    word pointSetName;
-    args.readIfPresent("pointSet", pointSetName);
-
-    // Define sub-directory name to use for VTK data.
-    const word vtkDirName = args.lookupOrDefault<word>("name", "VTK");
-
-    #include "createNamedMesh.H"
-
-    // VTK/ directory in the case
-    fileName fvPath(runTime.path()/vtkDirName);
-
-    // Directory of mesh (region0 gets filtered out)
-    fileName regionPrefix;
-    if (regionName != polyMesh::defaultRegion)
-    {
-        fvPath = fvPath/regionName;
-        regionPrefix = regionName;
-    }
-
-    if (isDir(fvPath))
-    {
-        if
-        (
-            args.found("time")
-         || args.found("latestTime")
-         || cellSubsetName.size()
-         || faceSetName.size()
-         || pointSetName.size()
-         || regionName != polyMesh::defaultRegion
-        )
+        if (args.readListIfPresent<wordRe>("excludePatches", excludePatches))
         {
-            Info<< "Keeping old VTK files in " << fvPath << nl << endl;
-        }
-        else
-        {
-            Info<< "Deleting old VTK files in " << fvPath << nl << endl;
-
-            rmDir(fvPath);
+            Info<< "Excluding patches " << flatOutput(excludePatches)
+                << nl << endl;
         }
     }
 
-    mkDir(fvPath);
-
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
-    // Mesh wrapper: does subsetting and decomposition
-    fvMeshSubsetProxy meshProxy(mesh, cellSubsetType, cellSubsetName);
-
-    // Collect decomposition information etc.
-    vtk::vtuCells vtuMeshCells(fmtType, decomposePoly);
+    wordRes selectedFields;
+    const bool useFieldFilter =
+        args.readListIfPresent<wordRe>("fields", selectedFields);
 
-    Info<< "VTK mesh topology: "
-        << timer.cpuTimeIncrement() << " s, "
-        << mem.update().size() << " kB" << endl;
+    #include "createTime.H"
 
-    #include "findClouds.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
 
-    // Supported volume field types
-    const wordHashSet vFieldTypes
-    {
-        volScalarField::typeName,
-        volVectorField::typeName,
-        volSphericalTensorField::typeName,
-        volSymmTensorField::typeName,
-        volTensorField::typeName
-    };
-
-    // Supported dimensioned field types
-    const wordHashSet dFieldTypes
-    {
-        volScalarField::Internal::typeName,
-        volVectorField::Internal::typeName,
-        volSphericalTensorField::Internal::typeName,
-        volSymmTensorField::Internal::typeName,
-        volTensorField::Internal::typeName
-    };
-
-    // Supported point field types
-    const wordHashSet pFieldTypes
-    {
-        pointScalarField::typeName,
-        pointVectorField::typeName,
-        pointSphericalTensorField::typeName,
-        pointSymmTensorField::typeName,
-        pointTensorField::typeName
-    };
-
-    // Supported cloud (lagrangian) field types
-    const wordHashSet cFieldTypes
-    {
-        labelIOField::typeName,
-        scalarIOField::typeName,
-        vectorIOField::typeName,
-        symmTensorIOField::typeName,
-        tensorIOField::typeName
-    };
+    // Information for file series
+    HashTable<vtk::seriesWriter, fileName> vtkSeries;
 
-    forAll(timeDirs, timei)
+    wordList regionNames;
+    wordRes selectRegions;
+    if (allRegions)
     {
-        runTime.setTime(timeDirs[timei], timei);
-
-        Info<< "Time: " << runTime.timeName() << endl;
-
-        const word timeDesc =
-            useTimeName ? runTime.timeName() : Foam::name(runTime.timeIndex());
-
-        // Check for new polyMesh/ and update mesh, fvMeshSubset and cell
-        // decomposition.
-        polyMesh::readUpdateState meshState = meshProxy.readUpdate();
+        regionNames =
+            regionProperties(runTime, IOobject::READ_IF_PRESENT).names();
 
-        const fvMesh& mesh = meshProxy.mesh();
-        if
-        (
-            meshState == polyMesh::TOPO_CHANGE
-         || meshState == polyMesh::TOPO_PATCH_CHANGE
-        )
-        {
-            // Trigger change for vtk cells too
-            vtuMeshCells.clear();
-        }
-
-
-        // Attempt topoSets first
-        bool shortCircuit = false;
-
-        // If faceSet: write faceSet only (as polydata)
-        if (faceSetName.size())
+        if (regionNames.empty())
         {
-            // Load
-            faceSet set(mesh, faceSetName);
-
-            // Filename as if patch with same name.
-            const fileName outputName
-            (
-                fvPath/set.name()
-              / set.name() + "_" + timeDesc
-            );
-
-            mkDir(outputName.path());
-
-            Info<< "    faceSet   : "
-                << outputName.relative(runTime.path()) << nl;
-
-            vtk::writeFaceSet
-            (
-                meshProxy.mesh(),
-                set,
-                fmtType,
-                outputName,
-                false // In parallel (later)
-            );
+            Info<<"Warning: "
+                << "No regionProperties - assuming default region"
+                << nl << endl;
 
-            shortCircuit = true;
+            regionNames.resize(1);
+            regionNames.first() = fvMesh::defaultRegion;
         }
-
-        // If pointSet: write pointSet only (as polydata)
-        if (pointSetName.size())
+        else
         {
-            // Load
-            pointSet set(mesh, pointSetName);
-
-            // Filename as if patch with same name.
-            const fileName outputName
-            (
-                fvPath/set.name()
-              / set.name() + "_" + timeDesc
-            );
-
-            mkDir(outputName.path());
-
-            Info<< "    pointSet  : "
-                << outputName.relative(runTime.path()) << nl;
-
-            vtk::writePointSet
-            (
-                meshProxy.mesh(),
-                set,
-                fmtType,
-                outputName,
-                false // In parallel (later)
-            );
-
-            shortCircuit = true;
+            Info<< "Using all regions in regionProperties" << nl
+                << "    "<< flatOutput(regionNames) << nl;
         }
-
-        if (shortCircuit)
+    }
+    else if (args.readListIfPresent<wordRe>("regions", selectRegions))
+    {
+        if (selectRegions.empty())
         {
-            continue;
+            regionNames.resize(1);
+            regionNames.first() = fvMesh::defaultRegion;
         }
-
-
-        // Search for list of objects for this time
-        IOobjectList objects(mesh, runTime.timeName());
-
-        wordHashSet selectedFields;
-        const bool specifiedFields = args.readIfPresent
-        (
-            "fields",
-            selectedFields
-        );
-
-        // Construct the vol fields
-        // References the original mesh, but uses subsetted portion only.
-
-        PtrList<const volScalarField> vScalarFld;
-        PtrList<const volVectorField> vVectorFld;
-        PtrList<const volSphericalTensorField> vSphTensorf;
-        PtrList<const volSymmTensorField> vSymTensorFld;
-        PtrList<const volTensorField> vTensorFld;
-
-        if
+        else if
         (
-            candidateObjects
-            (
-                objects,
-                vFieldTypes,
-                specifiedFields,
-                selectedFields
-            ).size()
+            selectRegions.size() == 1 && !selectRegions.first().isPattern()
         )
         {
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                vScalarFld
-            );
-            print("    volScalar        :", Info, vScalarFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                vVectorFld
-            );
-            print("    volVector        :", Info, vVectorFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                vSphTensorf
-            );
-            print("    volSphTensor     :", Info, vSphTensorf);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                vSymTensorFld
-            );
-            print("    volSymmTensor    :", Info, vSymTensorFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                vTensorFld
-            );
-            print("    volTensor        :", Info, vTensorFld);
+            regionNames.resize(1);
+            regionNames.first() = selectRegions.first();
         }
-
-        const label nVolFields =
-        (
-            vScalarFld.size()
-          + vVectorFld.size()
-          + vSphTensorf.size()
-          + vSymTensorFld.size()
-          + vTensorFld.size()
-        );
-
-
-        // Construct dimensioned fields
-        PtrList<const volScalarField::Internal> dScalarFld;
-        PtrList<const volVectorField::Internal> dVectorFld;
-        PtrList<const volSphericalTensorField::Internal> dSphTensorFld;
-        PtrList<const volSymmTensorField::Internal> dSymTensorFld;
-        PtrList<const volTensorField::Internal> dTensorFld;
-
-        if
-        (
-            candidateObjects
-            (
-                objects,
-                dFieldTypes,
-                specifiedFields,
-                selectedFields
-            ).size()
-        )
+        else
         {
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                dScalarFld
-            );
-            print("    volScalar::Internal      :", Info, dScalarFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                dVectorFld
-            );
-            print("    volVector::Internal      :", Info, dVectorFld);
+            regionNames =
+                regionProperties(runTime, IOobject::READ_IF_PRESENT).names();
 
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                dSphTensorFld
-            );
-            print("    volSphTensor::Internal   :", Info, dSphTensorFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                dSymTensorFld
-            );
-            print("    volSymmTensor::Internal  :", Info, dSymTensorFld);
-
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                dTensorFld
-            );
-            print("    volTensor::Internal      :", Info, dTensorFld);
-        }
-
-        const label nDimFields =
-        (
-            dScalarFld.size()
-          + dVectorFld.size()
-          + dSphTensorFld.size()
-          + dSymTensorFld.size()
-          + dTensorFld.size()
-        );
-
-
-        // Finite-area mesh and fields - need not exist
-
-        if (args.found("finiteAreaFields"))
-        {
-            autoPtr<faMesh> aMeshPtr;
+            if (regionNames.empty())
             {
-                const bool throwing = FatalError.throwExceptions();
-                try
-                {
-                    aMeshPtr.reset(new faMesh(meshProxy.baseMesh()));
-                }
-                catch (Foam::error& err)
-                {
-                    aMeshPtr.clear();
-                }
-                FatalError.throwExceptions(throwing);
-            }
+                Info<<"Warning: "
+                    << "No regionProperties - assuming default region"
+                    << nl << endl;
 
-            if (aMeshPtr.valid())
+                regionNames.resize(1);
+                regionNames.first() = fvMesh::defaultRegion;
+            }
+            else
             {
-                // Construct the area fields
+                inplaceSubsetStrings(selectRegions, regionNames);
 
-                PtrList<const areaScalarField> aScalarFld;
-                PtrList<const areaVectorField> aVectorFld;
-                PtrList<const areaSphericalTensorField> aSphTensorf;
-                PtrList<const areaSymmTensorField> aSymTensorFld;
-                PtrList<const areaTensorField> aTensorFld;
-
-                const faMesh& aMesh = aMeshPtr();
-
-                if (!specifiedFields || selectedFields.size())
+                if (regionNames.empty())
                 {
-                    readFields
-                    (
-                        aMesh,
-                        objects,
-                        selectedFields,
-                        aScalarFld
-                    );
-                    print("    areaScalar           :", Info, aScalarFld);
-
-                    readFields
-                    (
-                        aMesh,
-                        objects,
-                        selectedFields,
-                        aVectorFld
-                    );
-                    print("    areaVector           :", Info, aVectorFld);
-
-                    readFields
-                    (
-                        aMesh,
-                        objects,
-                        selectedFields,
-                        aSphTensorf
-                    );
-                    print("    areaSphericalTensor   :", Info, aSphTensorf);
-
-                    readFields
-                    (
-                        aMesh,
-                        objects,
-                        selectedFields,
-                        aSymTensorFld
-                    );
-                    print("    areaSymmTensor        :", Info, aSymTensorFld);
-
-                    readFields
-                    (
-                        aMesh,
-                        objects,
-                        selectedFields,
-                        aTensorFld
-                    );
-                    print("    areaTensor            :", Info, aTensorFld);
+                    Info<< "No matching regions ... stopping" << nl << endl;
+                    return 1;
                 }
 
-                const label nAreaFields =
-                (
-                    aScalarFld.size()
-                  + aVectorFld.size()
-                  + aSphTensorf.size()
-                  + aSymTensorFld.size()
-                  + aTensorFld.size()
-                );
-
-                fileName outputName(fvPath/"finiteArea");
-
-                mkDir(outputName);
+                Info<< "Using matching regions: "
+                    << flatOutput(regionNames) << nl;
+            }
+        }
+    }
+    else
+    {
+        regionNames.resize(1);
+        regionNames.first() =
+            args.lookupOrDefault<word>("region", fvMesh::defaultRegion);
+    }
 
-                const auto& pp = aMesh.patch();
 
-                vtk::surfaceMeshWriter writer
-                (
-                    pp,
-                    aMesh.name(),
-                    outputName/"finiteArea" + "_" + timeDesc,
-                    fmtType
-                );
+    // Names for sets and zones
+    word cellSelectionName;
+    word faceSetName;
+    word pointSetName;
 
-                // Number of fields
-                writer.beginCellData(nAreaFields);
+    fvMeshSubsetProxy::subsetType cellSubsetType = fvMeshSubsetProxy::NONE;
 
-                writer.write(aScalarFld);
-                writer.write(aVectorFld);
-                writer.write(aSphTensorf);
-                writer.write(aSymTensorFld);
-                writer.write(aTensorFld);
+    string vtkName = runTime.globalCaseName();
 
-                writer.endCellData();
+    if (regionNames.size() == 1)
+    {
+        if (args.readIfPresent("cellSet", cellSelectionName))
+        {
+            vtkName = cellSelectionName;
+            cellSubsetType = fvMeshSubsetProxy::SET;
 
-                writer.writeFooter();
-            }
+            Info<< "Converting cellSet " << cellSelectionName
+                << " only. New outside faces as \"oldInternalFaces\"."
+                << nl;
         }
+        else if (args.readIfPresent("cellZone", cellSelectionName))
+        {
+            vtkName = cellSelectionName;
+            cellSubsetType = fvMeshSubsetProxy::ZONE;
 
-        PtrList<const pointScalarField> pScalarFld;
-        PtrList<const pointVectorField> pVectorFld;
-        PtrList<const pointSphericalTensorField> pSphTensorFld;
-        PtrList<const pointSymmTensorField> pSymTensorFld;
-        PtrList<const pointTensorField> pTensorFld;
+            Info<< "Converting cellZone " << cellSelectionName
+                << " only. New outside faces as \"oldInternalFaces\"."
+                << nl;
+        }
 
-        // Construct pointMesh only if necessary since it constructs edge
-        // addressing (expensive on polyhedral meshes)
-        if
+        args.readIfPresent("faceSet", faceSetName);
+        args.readIfPresent("pointSet", pointSetName);
+    }
+    else
+    {
+        for
         (
-            !noPointValues
-         && candidateObjects
-            (
-                objects,
-                pFieldTypes,
-                specifiedFields,
-                selectedFields
-            ).size()
+            const word& opt
+          : { "cellSet", "cellZone", "faceSet", "pointSet" }
         )
         {
-            const pointMesh& ptMesh = pointMesh::New(meshProxy.baseMesh());
+            if (args.found(opt))
+            {
+                Info<< "Ignoring -" << opt << " for multi-regions" << nl;
+            }
+        }
+    }
 
-            readFields
-            (
-                meshProxy,
-                ptMesh,
-                objects,
-                selectedFields,
-                pScalarFld
-            );
-            print("    pointScalar      :", Info, pScalarFld);
 
-            readFields
-            (
-                meshProxy,
-                ptMesh,
-                objects,
-                selectedFields,
-                pVectorFld
-            );
-            print("    pointVector      :", Info, pVectorFld);
+    cpuTime timer;
+    memInfo mem;
+    Info<< "Initial memory " << mem.update().size() << " kB" << endl;
 
-            readFields
-            (
-                meshProxy,
-                ptMesh,
-                objects,
-                selectedFields,
-                pSphTensorFld
-            );
-            print("    pointSphTensor   : ", Info, pSphTensorFld);
+    #include "createMeshes.H"
 
-            readFields
-            (
-                meshProxy,
-                ptMesh,
-                objects,
-                selectedFields,
-                pSymTensorFld
-            );
-            print("    pointSymmTensor  :", Info, pSymTensorFld);
+    // Directory management
 
-            readFields
-            (
-                meshProxy,
-                ptMesh,
-                objects,
-                selectedFields,
-                pTensorFld
-            );
-            print("    pointTensor      :", Info, pTensorFld);
-        }
+    // Sub-directory for output
+    const word vtkDirName = args.lookupOrDefault<word>("name", "VTK");
 
-        const label nPointFields =
-            pScalarFld.size()
-          + pVectorFld.size()
-          + pSphTensorFld.size()
-          + pSymTensorFld.size()
-          + pTensorFld.size();
+    const fileName outputDir(runTime.globalPath()/vtkDirName);
 
-        if (doWriteInternal)
+    if (Pstream::master())
+    {
+        for (const word& regionName : regionNames)
         {
-            if (vtuMeshCells.empty())
-            {
-                // subMesh or baseMesh
-                vtuMeshCells.reset(meshProxy.mesh());
-            }
-
-            // Create file and write header
-            fileName outputName
-            (
-                fvPath/vtkName
-              + "_"
-              + timeDesc
-            );
-            Info<< "    Internal  : "
-                << outputName.relative(runTime.path()) << nl;
-
-            // Write mesh
-            vtk::internalWriter writer
-            (
-                meshProxy.mesh(),
-                vtuMeshCells,
-                outputName,
-                fmtType
-            );
+            // VTK/regionName  directory in the case
 
-            // CellData
+            fileName regionDir;
+            if (regionName != polyMesh::defaultRegion)
             {
-                writer.beginCellData(1 + nVolFields + nDimFields);
-
-                // Write cellID field
-                writer.writeCellIDs();
-
-                // Write volFields
-                writer.write(vScalarFld);
-                writer.write(vVectorFld);
-                writer.write(vSphTensorf);
-                writer.write(vSymTensorFld);
-                writer.write(vTensorFld);
-
-                // Write dimensionedFields
-                writer.write(dScalarFld);
-                writer.write(dVectorFld);
-                writer.write(dSphTensorFld);
-                writer.write(dSymTensorFld);
-                writer.write(dTensorFld);
-
-                writer.endCellData();
+                regionDir = outputDir / regionName;
             }
 
-            // PointData
-            if (!noPointValues)
+            if (args.found("overwrite") && isDir(regionDir))
             {
-                writer.beginPointData(nVolFields + nDimFields + nPointFields);
-
-                // pointFields
-                writer.write(pScalarFld);
-                writer.write(pVectorFld);
-                writer.write(pSphTensorFld);
-                writer.write(pSymTensorFld);
-                writer.write(pTensorFld);
-
-                // Interpolated volFields
-                volPointInterpolation pInterp(mesh);
-
-                writer.write(pInterp, vScalarFld);
-                writer.write(pInterp, vVectorFld);
-                writer.write(pInterp, vSphTensorf);
-                writer.write(pInterp, vSymTensorFld);
-                writer.write(pInterp, vTensorFld);
-
-                writer.write(pInterp, dScalarFld);
-                writer.write(pInterp, dVectorFld);
-                writer.write(pInterp, dSphTensorFld);
-                writer.write(pInterp, dSymTensorFld);
-                writer.write(pInterp, dTensorFld);
-
-                writer.endPointData();
+                Info<< "Deleting old VTK files in "
+                    << regionDir.relative(runTime.globalPath())
+                    << nl << endl;
+                rmDir(regionDir);
             }
-
-            writer.writeFooter();
+            mkDir(regionDir);
         }
+    }
 
-        //---------------------------------------------------------------------
-        //
-        // Write surface fields
-        //
-        //---------------------------------------------------------------------
 
-        if (args.found("surfaceFields"))
-        {
-            PtrList<const surfaceScalarField> sScalarFld;
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                sScalarFld
-            );
-            print("    surfScalar   :", Info, sScalarFld);
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-            PtrList<const surfaceVectorField> sVectorFld;
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                sVectorFld
-            );
-            print("    surfVector   :", Info, sVectorFld);
+    forAll(timeDirs, timei)
+    {
+        runTime.setTime(timeDirs[timei], timei);
 
-            if (sScalarFld.size())
-            {
-                // Rework the scalar fields into vector fields.
-                const label sz = sVectorFld.size();
+        const word timeDesc = "_" + Foam::name(runTime.timeIndex());
+        const scalar timeValue = runTime.value();
 
-                sVectorFld.setSize(sz + sScalarFld.size());
+        Info<< "Time: " << runTime.timeName() << endl;
 
-                surfaceVectorField n(mesh.Sf()/mesh.magSf());
 
-                forAll(sScalarFld, i)
-                {
-                    surfaceVectorField* ssfPtr = (sScalarFld[i]*n).ptr();
-                    ssfPtr->rename(sScalarFld[i].name());
-                    sVectorFld.set(sz+i, ssfPtr);
-                }
-                sScalarFld.clear();
-            }
+        // Accumulate information for multi-region VTM
+        vtk::vtmWriter vtmMultiRegion;
 
-            if (sVectorFld.size())
-            {
-                mkDir(fvPath / "surfaceFields");
+        // vtmMultiRegion.set(vtkDir/vtkName + timeDesc)
 
-                fileName outputName
-                (
-                    fvPath
-                  / "surfaceFields"
-                  / "surfaceFields"
-                  + "_"
-                  + timeDesc
-                );
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
 
-                vtk::writeSurfFields
-                (
-                    meshProxy.mesh(),
-                    outputName,
-                    fmtType,
-                    sVectorFld
-                );
+            fileName regionPrefix;
+            if (regionName != polyMesh::defaultRegion)
+            {
+                regionPrefix = regionName;
             }
-        }
 
+            auto& meshProxy = meshProxies[regioni];
+            auto& vtuMeshCells = vtuMappings[regioni];
 
-        //---------------------------------------------------------------------
-        //
-        // Write patches (POLYDATA file, one for each patch)
-        //
-        //---------------------------------------------------------------------
+            // polyMesh::readUpdateState meshState = mesh.readUpdate();
 
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+            // Check for new polyMesh/ and update mesh, fvMeshSubset
+            // and cell decomposition.
+            polyMesh::readUpdateState meshState =
+                meshProxy.readUpdate();
 
-        if (allPatches)
-        {
-            mkDir(fvPath/"allPatches");
+            const fvMesh& mesh = meshProxy.mesh();
 
-            fileName outputName
+            if
             (
-                fvPath/"allPatches"
-              / (meshProxy.useSubMesh() ? cellSubsetName : "allPatches")
-              + "_"
-              + timeDesc
-            );
-            Info<< "    Combined patches    : "
-                << outputName.relative(runTime.path()) << nl;
-
-            vtk::patchWriter writer
-            (
-                meshProxy.mesh(),
-                outputName,
-                fmtType,
-                nearCellValue,
-                getSelectedPatches(patches, excludePatches)
-            );
-
-            // CellData
+                meshState == polyMesh::TOPO_CHANGE
+             || meshState == polyMesh::TOPO_PATCH_CHANGE
+            )
             {
-                writer.beginCellData(1 + nVolFields);
-
-                // Write patchID field
-                writer.writePatchIDs();
-
-                // Write volFields
-                writer.write(vScalarFld);
-                writer.write(vVectorFld);
-                writer.write(vSphTensorf);
-                writer.write(vSymTensorFld);
-                writer.write(vTensorFld);
-
-                writer.endCellData();
+                // Trigger change for vtk cells too
+                vtuMeshCells.clear();
             }
 
-            // PointData
-            if (!noPointValues)
+            // Write topoSets before attempting anything else
             {
-                writer.beginPointData(nPointFields);
-
-                // Write pointFields
-                writer.write(pScalarFld);
-                writer.write(pVectorFld);
-                writer.write(pSphTensorFld);
-                writer.write(pSymTensorFld);
-                writer.write(pTensorFld);
-
-                // no interpolated volFields to avoid creating
-                // patchInterpolation for all subpatches.
-
-                writer.endPointData();
-            }
-
-            writer.writeFooter();
-        }
-        else
-        {
-            forAll(patches, patchi)
-            {
-                const polyPatch& pp = patches[patchi];
-
-                if (excludePatches.match(pp.name()))
+                #include "convertTopoSet.H"
+                if (wroteTopoSet)
                 {
-                    // Skip excluded patch
                     continue;
                 }
-
-                mkDir(fvPath/pp.name());
-
-                fileName outputName
-                (
-                    fvPath/pp.name()
-                  / (meshProxy.useSubMesh() ? cellSubsetName : pp.name())
-                  + "_"
-                  + timeDesc
-                );
-                Info<< "    Patch     : "
-                    << outputName.relative(runTime.path()) << nl;
-
-                vtk::patchWriter writer
-                (
-                    meshProxy.mesh(),
-                    outputName,
-                    fmtType,
-                    nearCellValue,
-                    labelList{patchi}
-                );
-
-                if (!isA<emptyPolyPatch>(pp))
-                {
-                    // VolFields + patchID
-                    writer.beginCellData(1 + nVolFields);
-
-                    // Write patchID field
-                    writer.writePatchIDs();
-
-                    // Write volFields
-                    writer.write(vScalarFld);
-                    writer.write(vVectorFld);
-                    writer.write(vSphTensorf);
-                    writer.write(vSymTensorFld);
-                    writer.write(vTensorFld);
-
-                    writer.endCellData();
-
-                    if (!noPointValues)
-                    {
-                        writer.beginPointData(nVolFields + nPointFields);
-
-                        // Write pointFields
-                        writer.write(pScalarFld);
-                        writer.write(pVectorFld);
-                        writer.write(pSphTensorFld);
-                        writer.write(pSymTensorFld);
-                        writer.write(pTensorFld);
-
-                        PrimitivePatchInterpolation<primitivePatch> pInter
-                        (
-                            pp
-                        );
-
-                        // Write interpolated volFields
-                        writer.write(pInter, vScalarFld);
-                        writer.write(pInter, vVectorFld);
-                        writer.write(pInter, vSphTensorf);
-                        writer.write(pInter, vSymTensorFld);
-                        writer.write(pInter, vTensorFld);
-
-                        writer.endPointData();
-                    }
-                }
-
-                writer.writeFooter();
             }
-        }
-
-        //---------------------------------------------------------------------
-        //
-        // Write faceZones (POLYDATA file, one for each zone)
-        //
-        //---------------------------------------------------------------------
 
-        if (doFaceZones && !mesh.faceZones().empty())
-        {
-            PtrList<const surfaceScalarField> sScalarFld;
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                sScalarFld
-            );
-            print("    surfScalar   :", Info, sScalarFld);
-
-            PtrList<const surfaceVectorField> sVectorFld;
-            readFields
-            (
-                meshProxy,
-                meshProxy.baseMesh(),
-                objects,
-                selectedFields,
-                sVectorFld
-            );
-            print("    surfVector   :", Info, sVectorFld);
+            // Search for list of objects for this time
+            IOobjectList objects(meshProxy.baseMesh(), runTime.timeName());
 
-            for (const faceZone& fz : mesh.faceZones())
+            if (useFieldFilter)
             {
-                mkDir(fvPath/fz.name());
-
-                fileName outputName =
-                (
-                    fvPath/fz.name()
-                  / (meshProxy.useSubMesh() ? cellSubsetName : fz.name())
-                  + "_"
-                  + timeDesc
-                );
-                Info<< "    FaceZone  : "
-                    << outputName.relative(runTime.path()) << nl;
+                objects.filterObjects(selectedFields);
+            }
 
-                indirectPrimitivePatch pp
-                (
-                    IndirectList<face>(mesh.faces(), fz),
-                    mesh.points()
-                );
+            // Prune restart fields
+            objects.prune_0();
 
-                vtk::surfaceMeshWriter writer
+            if (noPointValues)
+            {
+                // Prune point fields unless specifically requested
+                objects.filterClasses
                 (
-                    pp,
-                    fz.name(),
-                    outputName,
-                    fmtType
+                    [](const word& clsName)
+                    {
+                        return fieldTypes::point.found(clsName);
+                    },
+                    true // prune
                 );
+            }
 
-                // Number of fields
-                writer.beginCellData(sScalarFld.size() + sVectorFld.size());
+            // Volume, internal, point fields
+            #include "convertVolumeFields.H"
 
-                writer.write(sScalarFld);
-                writer.write(sVectorFld);
+            // Surface fields
+            #include "convertSurfaceFields.H"
 
-                writer.endCellData();
+            // Finite-area mesh and fields - need not exist
+            #include "convertAreaFields.H"
 
-                writer.writeFooter();
-            }
+            // Write lagrangian data
+            #include "convertLagrangian.H"
         }
 
-
-        //---------------------------------------------------------------------
-        //
-        // Write lagrangian data
-        //
-        //---------------------------------------------------------------------
-
-        for (const word& cloudName : cloudNames)
+        // Emit multi-region vtm
+        if (Pstream::master() && regionNames.size() > 1)
         {
-            // Always create the cloud directory.
-            mkDir(fvPath/cloud::prefix/cloudName);
-
             fileName outputName
             (
-                fvPath/cloud::prefix/cloudName/cloudName
-              + "_" + timeDesc
+                outputDir/vtkName + "-regions" + timeDesc + ".vtm"
             );
-            Info<< "    Lagrangian: "
-                << outputName.relative(runTime.path()) << nl;
-
-            IOobjectList cloudObjs
-            (
-                mesh,
-                runTime.timeName(),
-                cloud::prefix/cloudName
-            );
-
-            // Clouds require "coordinates".
-            // The "positions" are for v1706 and lower.
-            bool cloudExists =
-            (
-                cloudObjs.found("coordinates")
-             || cloudObjs.found("positions")
-            );
-            reduce(cloudExists, orOp<bool>());
-
-            if (cloudExists)
-            {
-                // Limited to types that we explicitly handle
-                HashTable<wordHashSet> cloudFields = cloudObjs.classes();
-                cloudFields.retain(cFieldTypes);
-
-                // The number of cloud fields (locally)
-                label nCloudFields = 0;
-                forAllConstIters(cloudFields, citer)
-                {
-                    nCloudFields += citer.object().size();
-                }
-
-                // Ensure all processes have identical information
-                if (Pstream::parRun())
-                {
-                    Pstream::mapCombineGather
-                    (
-                        cloudFields,
-                        HashSetOps::plusEqOp<word>()
-                    );
-                    Pstream::mapCombineScatter(cloudFields);
-                }
-
-
-                // Build lists of field names and echo some information
-
-                const wordList labelNames
-                (
-                    cloudFields(labelIOField::typeName).sortedToc()
-                );
-                print("        labels      :", Info, labelNames);
-
-                const wordList scalarNames
-                (
-                    cloudFields(scalarIOField::typeName).sortedToc()
-                );
-                print("        scalars     :", Info, scalarNames);
 
-                const wordList vectorNames
-                (
-                    cloudFields(vectorIOField::typeName).sortedToc()
-                );
-                print("        vectors     :", Info, vectorNames);
-
-                const wordList sphNames
-                (
-                    cloudFields(sphericalTensorIOField::typeName).sortedToc()
-                );
-                print("        sphTensors  :", Info, sphNames);
-
-                const wordList symmNames
-                (
-                    cloudFields(symmTensorIOField::typeName).sortedToc()
-                );
-                print("        symmTensors :", Info, symmNames);
-
-                const wordList tensorNames
-                (
-                    cloudFields(tensorIOField::typeName).sortedToc()
-                );
-                print("        tensors     :", Info, tensorNames);
-
-                vtk::lagrangianWriter writer
-                (
-                    meshProxy.mesh(),
-                    cloudName,
-                    outputName,
-                    fmtType
-                );
+            vtmMultiRegion.setTime(timeValue);
+            vtmMultiRegion.write(outputName);
 
-                // Write number of fields (on this processor)
-                writer.beginParcelData(nCloudFields);
+            fileName seriesName(vtk::seriesWriter::base(outputName));
 
-                // Fields
-                writer.writeIOField<label>(labelNames);
-                writer.writeIOField<scalar>(scalarNames);
-                writer.writeIOField<vector>(vectorNames);
-                writer.writeIOField<sphericalTensor>(sphNames);
-                writer.writeIOField<symmTensor>(symmNames);
-                writer.writeIOField<tensor>(tensorNames);
+            vtk::seriesWriter& series = vtkSeries(seriesName);
 
-                writer.endParcelData();
-
-                writer.writeFooter();
-            }
-            else
+            // First time?
+            // Load from file, verify against filesystem,
+            // prune time >= currentTime
+            if (series.empty())
             {
-                vtk::lagrangianWriter writer
-                (
-                    meshProxy.mesh(),
-                    cloudName,
-                    outputName,
-                    fmtType,
-                    true
-                );
-
-                // Write number of fields
-                writer.beginParcelData(0);
-
-                writer.endParcelData();
-
-                writer.writeFooter();
+                series.load(seriesName, true, timeValue);
             }
+
+            series.append(timeValue, outputName);
+            series.write(seriesName);
         }
 
         Info<< "Wrote in "
@@ -1587,68 +757,6 @@ int main(int argc, char *argv[])
     }
 
 
-    //---------------------------------------------------------------------
-    //
-    // Link parallel outputs back to undecomposed case for ease of loading
-    //
-    //---------------------------------------------------------------------
-
-    if (Pstream::parRun() && doLinks)
-    {
-        mkDir(runTime.path()/".."/vtkDirName);
-        chDir(runTime.path()/".."/vtkDirName);
-
-        Info<< "Linking all processor files to "
-            << runTime.path()/".."/vtkDirName
-            << endl;
-
-        // Get list of vtk files
-        fileName procVTK
-        (
-            fileName("..")
-          / "processor" + Foam::name(Pstream::myProcNo())
-          / vtkDirName
-        );
-
-        fileNameList dirs(readDir(procVTK, fileName::DIRECTORY));
-        label sz = dirs.size();
-        dirs.setSize(sz+1);
-        dirs[sz] = ".";
-
-        for (const fileName& subDir : dirs)
-        {
-            fileNameList subFiles(readDir(procVTK/subDir, fileName::FILE));
-
-            for (const fileName& subFile : subFiles)
-            {
-                fileName procFile(procVTK/subDir/subFile);
-
-                if (exists(procFile))
-                {
-                    // Could likely also use Foam::ln() directly
-                    List<string> cmd
-                    {
-                        "ln",
-                        "-s",
-                        procFile,
-                        (
-                            "processor"
-                          + Foam::name(Pstream::myProcNo())
-                          + "_"
-                          + procFile.name()
-                        )
-                    };
-
-                    if (Foam::system(cmd) == -1)
-                    {
-                        WarningInFunction
-                            << "Could not execute command " << cmd << endl;
-                    }
-                }
-            }
-        }
-    }
-
     Info<< "\nEnd: "
         << timer.elapsedCpuTime() << " s, "
         << mem.update().peak() << " kB (peak)\n" << endl;
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
index d67d0b1602ae3dc5a2e60a22544db1bd19bdf813..2b5b93e86214779f2e832e327858aad3d3c177af 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,91 +23,147 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "readFields.H"
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class GeoField>
+Foam::tmp<GeoField> Foam::getField
+(
+    const IOobject* io,
+    const typename GeoField::Mesh& mesh,
+    const bool syncPar
+)
+{
+    if (io)
+    {
+        return tmp<GeoField>::New(*io, mesh);
+    }
+
+    return nullptr;
+}
 
-namespace Foam
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::getField
+(
+    const IOobject* io,
+    const fvMeshSubsetProxy& proxy,
+    const bool syncPar
+)
 {
+    return
+        proxy.interpolate
+        (
+            getField<GeoField>(io, proxy.baseMesh(), syncPar)
+        );
+}
+
+
+template<class GeoField>
+Foam::tmp<GeoField> Foam::getField
+(
+    const typename GeoField::Mesh& mesh,
+    const IOobjectList& objects,
+    const word& fieldName,
+    const bool syncPar
+)
+{
+    // Can do something with syncPar on failure ...
+
+    return getField<GeoField>(objects.findObject(fieldName), mesh, syncPar);
+}
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 template<class GeoField>
-label readFields
+Foam::tmp<GeoField> Foam::getField
 (
     const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const word& fieldName,
+    const bool syncPar
+)
+{
+    // Can do something with syncPar on failure ...
+
+    return getField<GeoField>(objects.findObject(fieldName), proxy, syncPar);
+}
+
+
+template<class GeoField>
+Foam::PtrList<const GeoField> Foam::readFields
+(
     const typename GeoField::Mesh& mesh,
     const IOobjectList& objects,
-    const wordHashSet& selectedFields,
-    PtrList<const GeoField>& fields
+    const wordRes& selection
 )
 {
+    const bool syncPar = true;
+
     // Available fields of type GeoField, sorted order
     const wordList fieldNames =
     (
-        selectedFields.empty()
-      ? objects.sortedNames(GeoField::typeName)
-      : objects.sortedNames(GeoField::typeName, selectedFields)
+        selection.empty()
+      ? objects.sortedNames<GeoField>()
+      : objects.sortedNames<GeoField>(selection)
     );
 
     // Construct the fields
-    fields.resize(fieldNames.size());
+    PtrList<const GeoField> fields(fieldNames.size());
 
     label nFields = 0;
 
     for (const word& fieldName : fieldNames)
     {
-        fields.set
-        (
-            nFields++,
-            proxy.interpolate
-            (
-                GeoField(*(objects[fieldName]), mesh)
-            ).ptr()
-        );
+        auto tfield =
+            getField<GeoField>(mesh, objects, fieldName, syncPar);
+
+        if (tfield.valid())
+        {
+            fields.set(nFields++, tfield.ptr());
+        }
     }
 
-    return nFields;
+    fields.resize(nFields);
+    return fields;
 }
 
 
 template<class GeoField>
-label readFields
+Foam::PtrList<const GeoField> Foam::readFields
 (
-    const typename GeoField::Mesh& mesh,
+    const fvMeshSubsetProxy& proxy,
     const IOobjectList& objects,
-    const wordHashSet& selectedFields,
-    PtrList<const GeoField>& fields
+    const wordRes& selection
 )
 {
+    const bool syncPar = true;
+
     // Available fields of type GeoField, sorted order
     const wordList fieldNames =
     (
-        selectedFields.empty()
-      ? objects.sortedNames(GeoField::typeName)
-      : objects.sortedNames(GeoField::typeName, selectedFields)
+        selection.empty()
+      ? objects.sortedNames<GeoField>()
+      : objects.sortedNames<GeoField>(selection)
     );
 
     // Construct the fields
-    fields.resize(fieldNames.size());
+    PtrList<const GeoField> fields(fieldNames.size());
 
     label nFields = 0;
 
     for (const word& fieldName : fieldNames)
     {
-        fields.set
-        (
-            nFields++,
-            new GeoField(*(objects[fieldName]), mesh)
-        );
+        auto tfield =
+            getField<GeoField>(proxy, objects, fieldName, syncPar);
+
+        if (tfield.valid())
+        {
+            fields.set(nFields++, tfield.ptr());
+        }
     }
 
-    return nFields;
+    fields.resize(nFields);
+    return fields;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
index 672994d8ccfd961f35414ba54c5ef3bf279c6e7b..88be7e30bf000a3f41d5a1c121a62592b8807c02 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/readFields.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,8 @@ InClass
     Foam::readFields
 
 Description
+    Helper routines for reading a field or fields, optionally with
+    a mesh subset (using fvMeshSubsetProxy).
 
 SourceFiles
     readFields.C
@@ -35,35 +37,73 @@ SourceFiles
 #define readFields_H
 
 #include "fvMeshSubsetProxy.H"
-#include "HashSet.H"
-#include "PtrList.H"
 #include "IOobjectList.H"
+#include "PtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-//- Read the fields, optionally subset, and place on the pointer list
+//- Get the field or return nullptr
+template<class GeoField>
+tmp<GeoField> getField
+(
+    const IOobject* io,
+    const typename GeoField::Mesh& mesh,
+    const bool syncPar
+);
+
+
+//- Get the field and subset it, or return nullptr
 template<class GeoField>
-label readFields
+tmp<GeoField> getField
 (
+    const IOobject* io,
     const fvMeshSubsetProxy& proxy,
+    const bool syncPar
+);
+
+
+//- Get the named field from the objects, or return nullptr.
+template<class GeoField>
+tmp<GeoField> getField
+(
     const typename GeoField::Mesh& mesh,
     const IOobjectList& objects,
-    const wordHashSet& selectedFields,
-    PtrList<const GeoField>& fields
+    const word& fieldName,
+    const bool syncPar
 );
 
 
-//- Read the fields, and place on the pointer list
+//- Get the named field from the objects and subset it, or return nullptr
 template<class GeoField>
-label readFields
+tmp<GeoField> getField
+(
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const word& fieldName,
+    const bool syncPar
+);
+
+
+//- Read the fields, and return as a pointer list
+template<class GeoField>
+PtrList<const GeoField> readFields
 (
     const typename GeoField::Mesh& mesh,
     const IOobjectList& objects,
-    const wordHashSet& selectedFields,
-    PtrList<const GeoField>& fields
+    const wordRes& selection
+);
+
+
+//- Read the fields, optionally subset, and return as a pointer list
+template<class GeoField>
+PtrList<const GeoField> readFields
+(
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const wordRes& selection
 );
 
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/reportFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/reportFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..badf38232daa700a5fa0a9e8d29cd79db55b4da8
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/reportFields.H
@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::foamToVtkReportFields
+
+Description
+    Collection of simple static methods for reporting field names
+    by category, which is used by foamToVTK.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef foamToVtkReportFields_H
+#define foamToVtkReportFields_H
+
+#include "UPtrList.H"
+#include "Ostream.H"
+#include "areaFields.H"
+#include "volFields.H"
+#include "pointFields.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class foamToVtkReportFields Declaration
+\*---------------------------------------------------------------------------*/
+
+struct foamToVtkReportFields
+{
+    template<class GeoField>
+    static void print
+    (
+        const char* msg,
+        Ostream& os,
+        const UPtrList<const GeoField>& flds
+    )
+    {
+        if (flds.size())
+        {
+            os  << msg;
+            for (const GeoField& fld : flds)
+            {
+                os  << ' ' << fld.name();
+            }
+            os  << endl;
+        }
+    }
+
+
+    static void print
+    (
+        const char* msg,
+        Ostream& os,
+        const wordList& fieldNames
+    )
+    {
+        if (fieldNames.size())
+        {
+            os  << msg;
+            for (const word& fieldName : fieldNames)
+            {
+                os  << ' ' << fieldName;
+            }
+            os  << endl;
+        }
+    }
+
+
+    template<class FieldType>
+    static void print
+    (
+        const char* msg,
+        Ostream& os,
+        const IOobjectList& objects
+    )
+    {
+        print(msg, os, objects.sortedNames<FieldType>());
+    }
+
+
+    //- Supported volume field types
+    static void volume(Ostream& os, const IOobjectList& objects)
+    {
+        print<volScalarField>
+        (
+            "    volScalar    :", os, objects
+        );
+        print<volVectorField>
+        (
+            "    volVector    :", os, objects
+        );
+        print<volSphericalTensorField>
+        (
+            "    volSphTensor :", os, objects
+        );
+        print<volSymmTensorField>
+        (
+            "    volSymTensor :", os, objects
+        );
+        print<volTensorField>
+        (
+            "    volTensor    :", os, objects
+        );
+    }
+
+
+    //- Supported dimensioned field types
+    static void internal(Ostream& os, const IOobjectList& objects)
+    {
+        print<volScalarField::Internal>
+        (
+            "    volScalar:Internal    :", os, objects
+        );
+        print<volVectorField::Internal>
+        (
+            "    volVector:Internal    :", os, objects
+        );
+        print<volSphericalTensorField::Internal>
+        (
+            "    volSphTensor:Internal :", os, objects
+        );
+        print<volSymmTensorField::Internal>
+        (
+            "    volSymTensor:Internal :", os, objects
+        );
+        print<volTensorField::Internal>
+        (
+            "    volTensor:Internal    :", os, objects
+        );
+    }
+
+
+    //- Supported point field types
+    static void point(Ostream& os, const IOobjectList& objects)
+    {
+    }
+
+
+    //- Supported area field types
+    static void area(Ostream& os, const IOobjectList& objects)
+    {
+        print<areaScalarField>
+        (
+            "    areaScalar    :", os, objects
+        );
+        print<areaVectorField>
+        (
+            "    areaVector    :", os, objects
+        );
+        print<areaSphericalTensorField>
+        (
+            "    areaSphTensor :", os, objects
+        );
+        print<areaSymmTensorField>
+        (
+            "    areaSymTensor :", os, objects
+        );
+        print<areaTensorField>
+        (
+            "    areaTensor    :", os, objects
+        );
+    }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..469865f44310d6825800df7e47b2619a4e699610
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+InNamespace
+    Foam
+
+Description
+    Read finite-area fields from disk
+    and write with vtk::surfaceMeshWriter
+
+SourceFiles
+    writeAreaFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef writeAreaFields_H
+#define writeAreaFields_H
+
+#include "readFields.H"
+#include "foamVtkSurfaceMeshWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class GeoField>
+bool writeAreaField
+(
+    vtk::surfaceMeshWriter& writer,
+    const tmp<GeoField>& tfield
+)
+{
+    if (tfield.valid())
+    {
+        writer.write(tfield());
+        tfield.clear();
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class GeoField>
+label writeAreaFields
+(
+    vtk::surfaceMeshWriter& writer,
+    const typename GeoField::Mesh& mesh,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writeAreaField<GeoField>
+            (
+                writer,
+                getField<GeoField>(mesh, objects, fieldName, syncPar)
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+label writeAllAreaFields
+(
+    vtk::surfaceMeshWriter& writer,
+    const faMesh& mesh,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writeAreaFields<FieldType>              \
+        (                                       \
+            writer,                             \
+            mesh,                               \
+            objects,                            \
+            syncPar                             \
+        )
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(areaScalarField);
+    count += foamToVtk_WRITE_FIELD(areaVectorField);
+    count += foamToVtk_WRITE_FIELD(areaSphericalTensorField);
+    count += foamToVtk_WRITE_FIELD(areaSymmTensorField);
+    count += foamToVtk_WRITE_FIELD(areaTensorField);
+
+    #undef foamToVtk_WRITE_FIELD
+    return count;
+}
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeDimFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeDimFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..8a5e1720cdede7e2bbe3afe9e502465be1c3abc1
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeDimFields.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+InNamespace
+    Foam
+
+Description
+    Read dimensioned fields from disk
+    and write with vtk::internalWriter
+
+SourceFiles
+    writeDimFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef writeDimFields_H
+#define writeDimFields_H
+
+#include "readFields.H"
+#include "foamVtkInternalWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class GeoField>
+bool writeDimField
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const tmp<GeoField>& tfield
+)
+{
+    if (!tfield.valid())
+    {
+        return false;
+    }
+
+    const auto& field = tfield();
+
+    if (internalWriter.valid())
+    {
+        internalWriter->write(field);
+    }
+
+    tfield.clear();
+    return true;
+}
+
+
+template<class GeoField>
+bool writeDimField
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+    const tmp<GeoField>& tfield
+)
+{
+    if (!tfield.valid())
+    {
+        return false;
+    }
+
+    const auto& field = tfield();
+
+    if (internalWriter.valid() && pInterp.valid())
+    {
+        internalWriter->write(field, *pInterp);
+    }
+
+    tfield.clear();
+    return true;
+}
+
+
+template<class GeoField>
+label writeDimFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writeDimField<GeoField>
+            (
+                internalWriter,
+                getField<GeoField>(proxy, objects, fieldName, syncPar)
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+template<class GeoField>
+label writeDimFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writeDimField<GeoField>
+            (
+                internalWriter, pInterp,
+                getField<GeoField>(proxy, objects, fieldName, syncPar)
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+label writeAllDimFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writeDimFields<FieldType>               \
+        (                                       \
+            internalWriter,                     \
+            proxy,                              \
+            objects,                            \
+            syncPar                             \
+        )
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(volScalarField::Internal);
+    count += foamToVtk_WRITE_FIELD(volVectorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volSphericalTensorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volSymmTensorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volTensorField::Internal);
+
+    #undef foamToVTK_WRITE_FIELD
+    return count;
+}
+
+
+label writeAllDimFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writeDimFields<FieldType>               \
+        (                                       \
+            internalWriter, pInterp,            \
+            proxy,                              \
+            objects,                            \
+            syncPar                             \
+        )
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(volScalarField::Internal);
+    count += foamToVtk_WRITE_FIELD(volVectorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volSphericalTensorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volSymmTensorField::Internal);
+    count += foamToVtk_WRITE_FIELD(volTensorField::Internal);
+
+    #undef foamToVTK_WRITE_FIELD
+    return count;
+}
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..6d053ec78df81b9f376ec64e639188c003488c52
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writePointFields.H
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+InNamespace
+    Foam
+
+Description
+    Read point fields from disk
+    and write with vtk::internalWriter and vtk::patchWriter
+
+SourceFiles
+    writePointFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef writePointFields_H
+#define writePointFields_H
+
+#include "readFields.H"
+#include "foamVtkInternalWriter.H"
+#include "foamVtkPatchWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class GeoField>
+bool writePointField
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+
+    const tmp<GeoField>& tfield,
+    const fvMeshSubsetProxy& proxy
+)
+{
+    if (!tfield.valid())
+    {
+        return false;
+    }
+
+    tmp<GeoField> tproxied;
+    if (proxy.useSubMesh())
+    {
+        tproxied = proxy.interpolate(tfield());
+        tfield.clear();
+    }
+    else
+    {
+        tproxied = tfield;
+    }
+
+    if (!tproxied.valid())
+    {
+        // Or Error?
+        return false;
+    }
+
+
+    const auto& field = tproxied();
+
+    // Internal
+    if (internalWriter.valid())
+    {
+        internalWriter->write(field);
+    }
+
+    // Boundary
+    for (vtk::patchWriter& writer : patchWriters)
+    {
+        writer.write(field);
+    }
+
+
+    tproxied.clear();
+
+    return true;
+}
+
+
+template<class GeoField>
+label writePointFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+
+    const fvMeshSubsetProxy& proxy,
+    const typename GeoField::Mesh& ptMesh,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writePointField<GeoField>
+            (
+                internalWriter,
+                patchWriters,
+                getField<GeoField>(ptMesh, objects, fieldName, syncPar),
+                proxy
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+label writeAllPointFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    const pointMesh& ptMesh = pointMesh::New(proxy.baseMesh());
+
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writePointFields<FieldType>             \
+        (                                       \
+            internalWriter,                     \
+            patchWriters,                       \
+            proxy, ptMesh,                      \
+            objects,                            \
+            syncPar                             \
+        )
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(pointScalarField);
+    count += foamToVtk_WRITE_FIELD(pointVectorField);
+    count += foamToVtk_WRITE_FIELD(pointSphericalTensorField);
+    count += foamToVtk_WRITE_FIELD(pointSymmTensorField);
+    count += foamToVtk_WRITE_FIELD(pointTensorField);
+
+    #undef foamToVTK_WRITE_FIELD
+    return count;
+}
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeVolFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeVolFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..7dc84b88cf11f9c49945475152c02d890dca992e
--- /dev/null
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeVolFields.H
@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+InNamespace
+    Foam
+
+Description
+    Read volume fields from disk
+    and write with vtk::internalWriter and vtk::patchWriter
+
+SourceFiles
+    writeVolFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef writeVolFields_H
+#define writeVolFields_H
+
+#include "readFields.H"
+#include "foamVtkInternalWriter.H"
+#include "foamVtkPatchWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class GeoField>
+bool writeVolField
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+    const tmp<GeoField>& tfield
+)
+{
+    if (!tfield.valid())
+    {
+        return false;
+    }
+
+    const auto& field = tfield();
+
+    // Internal
+    if (internalWriter.valid())
+    {
+        internalWriter->write(field);
+    }
+
+    // Boundary
+    for (vtk::patchWriter& writer : patchWriters)
+    {
+        writer.write(field);
+    }
+
+    tfield.clear();
+    return true;
+}
+
+
+template<class GeoField>
+bool writeVolField
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+
+    UPtrList<vtk::patchWriter>& patchWriters,
+    const UPtrList<PrimitivePatchInterpolation<primitivePatch>>& patchInterps,
+
+    const tmp<GeoField>& tfield
+)
+{
+    if (!tfield.valid())
+    {
+        return false;
+    }
+
+    const auto& field = tfield();
+
+    // Internal
+    if (internalWriter.valid() && pInterp.valid())
+    {
+        internalWriter->write(field, *pInterp);
+    }
+
+    // Boundary
+    label writeri = 0;
+    for (vtk::patchWriter& writer : patchWriters)
+    {
+        if (writeri < patchInterps.size() && patchInterps.set(writeri))
+        {
+            writer.write(field, patchInterps[writeri]);
+        }
+        ++writeri;
+    }
+
+    tfield.clear();
+    return true;
+}
+
+
+template<class GeoField>
+label writeVolFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writeVolField<GeoField>
+            (
+                internalWriter,
+                patchWriters,
+                getField<GeoField>(proxy, objects, fieldName, syncPar)
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+template<class GeoField>
+label writeVolFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+
+    UPtrList<vtk::patchWriter>& patchWriters,
+    const UPtrList<PrimitivePatchInterpolation<primitivePatch>>& patchInterps,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    label count = 0;
+
+    for (const word& fieldName : objects.sortedNames<GeoField>())
+    {
+        if
+        (
+            writeVolField<GeoField>
+            (
+                internalWriter, pInterp,
+                patchWriters,   patchInterps,
+                getField<GeoField>(proxy, objects, fieldName, syncPar)
+            )
+        )
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+label writeAllVolFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    UPtrList<vtk::patchWriter>& patchWriters,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writeVolFields<FieldType>               \
+        (                                       \
+            internalWriter,                     \
+            patchWriters,                       \
+            proxy,                              \
+            objects,                            \
+            syncPar                             \
+        )
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(volScalarField);
+    count += foamToVtk_WRITE_FIELD(volVectorField);
+    count += foamToVtk_WRITE_FIELD(volSphericalTensorField);
+    count += foamToVtk_WRITE_FIELD(volSymmTensorField);
+    count += foamToVtk_WRITE_FIELD(volTensorField);
+
+    #undef foamToVTK_WRITE_FIELD
+    return count;
+}
+
+
+label writeAllVolFields
+(
+    autoPtr<vtk::internalWriter>& internalWriter,
+    const autoPtr<volPointInterpolation>& pInterp,
+
+    UPtrList<vtk::patchWriter>& patchWriters,
+    const UPtrList<PrimitivePatchInterpolation<primitivePatch>>& patchInterps,
+
+    const fvMeshSubsetProxy& proxy,
+    const IOobjectList& objects,
+    const bool syncPar
+)
+{
+    #undef  foamToVtk_WRITE_FIELD
+    #define foamToVtk_WRITE_FIELD(FieldType)    \
+        writeVolFields<FieldType>               \
+        (                                       \
+            internalWriter, pInterp,            \
+            patchWriters,   patchInterps,       \
+            proxy,                              \
+            objects,                            \
+            syncPar                             \
+        )
+
+
+    label count = 0;
+    count += foamToVtk_WRITE_FIELD(volScalarField);
+    count += foamToVtk_WRITE_FIELD(volVectorField);
+    count += foamToVtk_WRITE_FIELD(volSphericalTensorField);
+    count += foamToVtk_WRITE_FIELD(volSymmTensorField);
+    count += foamToVtk_WRITE_FIELD(volTensorField);
+
+    #undef foamToVtk_WRITE_FIELD
+    return count;
+}
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //