diff --git a/applications/utilities/mesh/manipulation/setSet/Make/files b/applications/utilities/mesh/manipulation/setSet/Make/files
index a14c37f5ee..75c578a1e5 100644
--- a/applications/utilities/mesh/manipulation/setSet/Make/files
+++ b/applications/utilities/mesh/manipulation/setSet/Make/files
@@ -1,5 +1,3 @@
-writePointSet.C
-writePatch.C
 setSet.C
 
 EXE = $(FOAM_APPBIN)/setSet
diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C
index e63a8b03e3..861e754e7b 100644
--- a/applications/utilities/mesh/manipulation/setSet/setSet.C
+++ b/applications/utilities/mesh/manipulation/setSet/setSet.C
@@ -41,8 +41,7 @@ Description
 #include "OFstream.H"
 #include "IFstream.H"
 #include "demandDrivenData.H"
-#include "writePatch.H"
-#include "writePointSet.H"
+#include "vtkWritePolyData.H"
 #include "IOobjectList.H"
 #include "cellZoneSet.H"
 #include "faceZoneSet.H"
@@ -80,7 +79,7 @@ void writeVTK
         // Faces of set with OpenFOAM faceID as value
 
         faceList setFaces(currentSet.size());
-        labelList faceValues(currentSet.size());
+        labelField faceValues(currentSet.size());
         label setFacei = 0;
 
         forAllConstIter(topoSet, currentSet, iter)
@@ -92,14 +91,18 @@ void writeVTK
 
         primitiveFacePatch fp(setFaces, mesh.points());
 
-        vtkWriteOps::writePatch
+        vtkWritePolyData::write
         (
-            true,
+            mesh.time().path()/vtkName,
             currentSet.name(),
-            fp,
+            true,
+            fp.localPoints(),
+            labelList(),
+            edgeList(),
+            fp.localFaces(),
             "faceID",
-            faceValues,
-            mesh.time().path()/vtkName
+            false,
+            faceValues
         );
     }
     else if (isA<cellSet>(currentSet))
@@ -140,7 +143,7 @@ void writeVTK
         }
 
         faceList setFaces(cellFaces.size());
-        labelList faceValues(cellFaces.size());
+        labelField faceValues(cellFaces.size());
         label setFacei = 0;
 
         forAllConstIter(Map<label>, cellFaces, iter)
@@ -152,24 +155,32 @@ void writeVTK
 
         primitiveFacePatch fp(setFaces, mesh.points());
 
-        vtkWriteOps::writePatch
+        vtkWritePolyData::write
         (
-            true,
+            mesh.time().path()/vtkName,
             currentSet.name(),
-            fp,
+            true,
+            fp.localPoints(),
+            labelList(),
+            edgeList(),
+            fp.localFaces(),
             "cellID",
-            faceValues,
-            mesh.time().path()/vtkName
+            false,
+            faceValues
         );
     }
     else if (isA<pointSet>(currentSet))
     {
-        vtkWriteOps::writePointSet
+        std::ofstream os(mesh.time().path()/vtkName);
+        vtkWritePolyData::write
         (
+            mesh.time().path()/vtkName,
+            currentSet.name(),
             true,
-            mesh,
-            currentSet,
-            mesh.time().path()/vtkName
+            pointField(mesh.points(), currentSet.toc()),
+            identity(currentSet.size()),
+            edgeList(),
+            faceList()
         );
     }
     else
diff --git a/applications/utilities/mesh/manipulation/setSet/writePatch.C b/applications/utilities/mesh/manipulation/setSet/writePatch.C
deleted file mode 100644
index 0aea556367..0000000000
--- a/applications/utilities/mesh/manipulation/setSet/writePatch.C
+++ /dev/null
@@ -1,119 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-     \\/     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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "writePatch.H"
-#include "OFstream.H"
-#include "vtkWriteOps.H"
-#include "primitiveFacePatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::vtkWriteOps::writePatch
-(
-    const bool binary,
-    const word& setName,
-    const primitiveFacePatch& fp,
-    const word& fieldName,
-    labelList& fieldValues,
-    const fileName& fileName
-)
-{
-    std::ofstream pStream(fileName.c_str());
-
-    pStream
-        << "# vtk DataFile Version 2.0" << std::endl
-        << setName << std::endl;
-    if (binary)
-    {
-        pStream << "BINARY" << std::endl;
-    }
-    else
-    {
-        pStream << "ASCII" << std::endl;
-    }
-    pStream << "DATASET POLYDATA" << std::endl;
-
-
-    //------------------------------------------------------------------
-    //
-    // Write topology
-    //
-    //------------------------------------------------------------------
-
-
-    // Write points and faces as polygons
-
-    pStream << "POINTS " << fp.nPoints() << " float" << std::endl;
-
-    DynamicList<floatScalar> ptField(3*fp.nPoints());
-
-    vtkWriteOps::insert(fp.localPoints(), ptField);
-
-    vtkWriteOps::write(pStream, binary, ptField);
-
-
-    label nFaceVerts = 0;
-
-    forAll(fp.localFaces(), facei)
-    {
-        nFaceVerts += fp.localFaces()[facei].size() + 1;
-    }
-    pStream << "POLYGONS " << fp.size() << ' ' << nFaceVerts
-        << std::endl;
-
-
-    DynamicList<label> vertLabels(nFaceVerts);
-
-    forAll(fp.localFaces(), facei)
-    {
-        const face& f = fp.localFaces()[facei];
-
-        vertLabels.append(f.size());
-
-        vtkWriteOps::insert(f, vertLabels);
-    }
-    vtkWriteOps::write(pStream, binary, vertLabels);
-
-
-    //-----------------------------------------------------------------
-    //
-    // Write data
-    //
-    //-----------------------------------------------------------------
-
-    // Write faceID
-
-    pStream
-        << "CELL_DATA " << fp.size() << std::endl
-        << "FIELD attributes 1" << std::endl;
-
-    // Cell ids first
-    pStream << fieldName << " 1 " << fp.size() << " int" << std::endl;
-
-    vtkWriteOps::write(pStream, binary, fieldValues);
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/setSet/writePatch.H b/applications/utilities/mesh/manipulation/setSet/writePatch.H
deleted file mode 100644
index 9bf616e3fc..0000000000
--- a/applications/utilities/mesh/manipulation/setSet/writePatch.H
+++ /dev/null
@@ -1,68 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-     \\/     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/>.
-
-InClass
-    Foam::writePatch
-
-Description
-    Write faceSet to vtk polydata file. Only one data which is original
-    faceID.
-
-SourceFiles
-    writePatch.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef writePatch_H
-#define writePatch_H
-
-#include "primitiveMesh.H"
-#include "primitiveFacePatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace vtkWriteOps
-{
-
-//- Write patch
-void writePatch
-(
-    const bool binary,
-    const word& setName,
-    const primitiveFacePatch& fp,
-    const word& fieldName,
-    labelList& fieldValues,
-    const fileName& fileName
-);
-
-} // End namespace vtkWriteOps
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/setSet/writePointSet.C b/applications/utilities/mesh/manipulation/setSet/writePointSet.C
deleted file mode 100644
index 0b3e52b08f..0000000000
--- a/applications/utilities/mesh/manipulation/setSet/writePointSet.C
+++ /dev/null
@@ -1,109 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-     \\/     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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "writePointSet.H"
-#include "OFstream.H"
-#include "vtkWriteOps.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::vtkWriteOps::writePointSet
-(
-    const bool binary,
-    const primitiveMesh& mesh,
-    const topoSet& set,
-    const fileName& fileName
-)
-{
-    std::ofstream pStream(fileName.c_str());
-
-    pStream
-        << "# vtk DataFile Version 2.0" << std::endl
-        << set.name() << std::endl;
-    if (binary)
-    {
-        pStream << "BINARY" << std::endl;
-    }
-    else
-    {
-        pStream << "ASCII" << std::endl;
-    }
-    pStream << "DATASET POLYDATA" << std::endl;
-
-
-    //------------------------------------------------------------------
-    //
-    // Write topology
-    //
-    //------------------------------------------------------------------
-
-
-    labelList pointLabels(set.toc());
-
-    pointField setPoints(mesh.points(), pointLabels);
-
-    // Write points
-
-    pStream << "POINTS " << pointLabels.size() << " float" << std::endl;
-
-    DynamicList<floatScalar> ptField(3*pointLabels.size());
-
-    vtkWriteOps::insert(setPoints, ptField);
-
-    vtkWriteOps::write(pStream, binary, ptField);
-
-    // Write vertices
-
-    pStream << "VERTICES " << pointLabels.size() << ' '
-        << 2*pointLabels.size() << std::endl;
-
-    DynamicList<label> vertexPoints(2*pointLabels.size());
-    forAll(pointLabels, pointi)
-    {
-        vertexPoints.append(1);
-        vertexPoints.append(pointi);
-    }
-    vtkWriteOps::write(pStream, binary, vertexPoints);
-
-    //-----------------------------------------------------------------
-    //
-    // Write data
-    //
-    //-----------------------------------------------------------------
-
-    // Write pointID
-
-    pStream
-        << "POINT_DATA " << pointLabels.size() << std::endl
-        << "FIELD attributes 1" << std::endl;
-
-    // Cell ids first
-    pStream << "pointID 1 " << pointLabels.size() << " int" << std::endl;
-
-    vtkWriteOps::write(pStream, binary, pointLabels);
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/setSet/writePointSet.H b/applications/utilities/mesh/manipulation/setSet/writePointSet.H
deleted file mode 100644
index 6309461971..0000000000
--- a/applications/utilities/mesh/manipulation/setSet/writePointSet.H
+++ /dev/null
@@ -1,67 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-     \\/     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
-    Write pointSet to vtk polydata file. Only one data which is original
-    pointID.
-
-SourceFiles
-    writePointSet.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef writePointSet_H
-#define writePointSet_H
-
-#include "primitiveMesh.H"
-#include "pointSet.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace vtkWriteOps
-{
-
-//- Write pointSet to vtk polydata file.
-//  Only one data which is original pointID.
-void writePointSet
-(
-    const bool binary,
-    const primitiveMesh& mesh,
-    const topoSet& set,
-    const fileName& fileName
-);
-
-} // End namespace vtkWriteOps
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/fileFormats/vtk/vtkWritePolyData.H b/src/fileFormats/vtk/vtkWritePolyData.H
new file mode 100644
index 0000000000..685f163428
--- /dev/null
+++ b/src/fileFormats/vtk/vtkWritePolyData.H
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     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::vtkWritePolyData
+
+Description
+    General write functions for vtk polygonal data files
+
+SourceFiles
+    vtkWritePolyData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vtkWritePolyData_H
+#define vtkWritePolyData_H
+
+#include "vtkWriteOps.H"
+#include "fileName.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace vtkWritePolyData Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace vtkWritePolyData
+{
+
+//- Helper for templated write
+template<class Types, class Type>
+inline void setFieldTypeValue
+(
+    UPtrList<const Field<Types>>& fieldTypeValues,
+    const label fieldi,
+    const Field<Type>& fieldTypeValue
+)
+{}
+
+
+//- Helper for templated write
+template<class Type>
+inline void setFieldTypeValue
+(
+    UPtrList<const Field<Type>>& fieldTypeValues,
+    const label fieldi,
+    const Field<Type>& fieldTypeValue
+)
+{
+    fieldTypeValues.set(fieldi, &fieldTypeValue);
+}
+
+
+//- Helper for templated write
+inline void unpackFieldTypeValues
+(
+    wordList& fieldNames,
+    boolList& fieldIsPointValues,
+    UPtrList<const Field<label>>& fieldLabelValues
+    #define FieldTypeValuesNonConstArg(Type, nullArg) \
+    , UPtrList<const Field<Type>>& field##Type##Values
+    FOR_ALL_FIELD_TYPES(FieldTypeValuesNonConstArg)
+    #undef FieldTypeValuesNonConstArg
+)
+{}
+
+
+//- Helper for templated write
+template<class Type, class ... Args>
+inline void unpackFieldTypeValues
+(
+    wordList& fieldNames,
+    boolList& fieldIsPointValues,
+    UPtrList<const Field<label>>& fieldLabelValues
+    #define FieldTypeValuesNonConstArg(Type, nullArg) \
+    , UPtrList<const Field<Type>>& field##Type##Values
+    FOR_ALL_FIELD_TYPES(FieldTypeValuesNonConstArg),
+    #undef FieldTypeValuesNonConstArg
+    const word& fieldName,
+    const bool fieldIsPointValue,
+    const Field<Type>& fieldTypeValue,
+    Args& ... args
+)
+{
+    const label fieldi = fieldNames.size() - 1 - sizeof...(Args)/3;
+
+    fieldNames[fieldi] = fieldName;
+    fieldIsPointValues[fieldi] = fieldIsPointValue;
+    setFieldTypeValue(fieldLabelValues, fieldi, fieldTypeValue);
+    #define SetFieldTypeValue(Type, nullArg) \
+        setFieldTypeValue(field##Type##Values, fieldi, fieldTypeValue);
+    FOR_ALL_FIELD_TYPES(SetFieldTypeValue);
+    #undef SetFieldTypeValue
+
+    unpackFieldTypeValues
+    (
+        fieldNames,
+        fieldIsPointValues,
+        fieldLabelValues
+        #define FieldTypeValuesParameter(Type, nullArg) , field##Type##Values
+        FOR_ALL_FIELD_TYPES(FieldTypeValuesParameter),
+        #undef FieldTypeValuesParameter
+        args ...
+    );
+}
+
+
+//- Write the field values out for a type
+template<class Type, class DataType>
+void writeFieldTypeValues
+(
+    std::ostream& os,
+    const bool binary,
+    const wordList& fieldNames,
+    const boolList& fieldIsPointValues,
+    const UPtrList<const Field<Type>>& fieldTypeValues,
+    const bool writePointValues
+);
+
+
+//- Write VTK polygonal data to a file. Takes a PtrList of fields of labels and
+//  of every primitive type. Each PtrList should be the length of the total
+//  number of fields and only one PtrList should be non-null for each field
+//  index.
+template<class PointField, class VertexList, class LineList, class FaceList>
+void write
+(
+    const fileName& file,
+    const word& title,
+    const bool binary,
+    const PointField& points,
+    const VertexList& vertices,
+    const LineList& lines,
+    const FaceList& faces,
+    const wordList& fieldNames,
+    const boolList& fieldIsPointValues,
+    const UPtrList<const Field<label>>& fieldLabelValues
+    #define FieldTypeValuesConstArg(Type, nullArg) \
+        , const UPtrList<const Field<Type>>& field##Type##Values
+    FOR_ALL_FIELD_TYPES(FieldTypeValuesConstArg)
+    #undef FieldTypeValuesConstArg
+);
+
+
+//- Write VTK polygonal data to a file. Takes any number of name,
+//  isPointValues, values arguments at the end. E.g.;
+//
+//  write
+//  (
+//      // Output options
+//      "myPolyData.vtk", "myPolyData", false,
+//
+//      // Geometry
+//      pp.localPoints(), labelList(), labelListList(), pp.localFaces(),
+//
+//      // Fields
+//      "faceIDs", true, pp.addressing(),
+//      "facePressures", false, Field<scalar>(pp.size(), ...),
+//      "pointVelocities", true, Field<vector>(pp.nPoints(), ...)
+//  );
+//
+template
+<
+    class PointField,
+    class VertexList,
+    class LineList,
+    class FaceList,
+    class ... Args
+>
+inline void write
+(
+    const fileName& file,
+    const word& title,
+    const bool binary,
+    const PointField& points,
+    const VertexList& vertices,
+    const LineList& lines,
+    const FaceList& faces,
+    const Args& ... args
+);
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "vtkWritePolyDataTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fileFormats/vtk/vtkWritePolyDataTemplates.C b/src/fileFormats/vtk/vtkWritePolyDataTemplates.C
new file mode 100644
index 0000000000..94596eebeb
--- /dev/null
+++ b/src/fileFormats/vtk/vtkWritePolyDataTemplates.C
@@ -0,0 +1,297 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vtkWritePolyData.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DataType>
+void Foam::vtkWritePolyData::writeFieldTypeValues
+(
+    std::ostream& os,
+    const bool binary,
+    const wordList& fieldNames,
+    const boolList& fieldIsPointValues,
+    const UPtrList<const Field<Type>>& fieldTypeValues,
+    const bool writePointValues
+)
+{
+    forAll(fieldNames, fieldi)
+    {
+        if
+        (
+            fieldIsPointValues[fieldi] == writePointValues
+         && fieldTypeValues.set(fieldi)
+        )
+        {
+            const label nCmpt = pTraits<Type>::nComponents;
+
+            os  << fieldNames[fieldi] << ' ' << pTraits<Type>::nComponents
+                << ' ' << fieldTypeValues[fieldi].size() << ' '
+                << (std::is_integral<DataType>::value ? "int" : "float") << nl;
+
+            List<DataType> data(nCmpt*fieldTypeValues[fieldi].size());
+            label i = 0;
+            forAll(fieldTypeValues[fieldi], fieldValuei)
+            {
+                for (direction cmpt = 0; cmpt < nCmpt; ++ cmpt)
+                {
+                    data[i ++] =
+                        component
+                        (
+                            fieldTypeValues[fieldi][fieldValuei],
+                            cmpt
+                        );
+                }
+            }
+
+            vtkWriteOps::write(os, binary, data);
+        }
+    }
+}
+
+
+template<class PointField, class VertexList, class LineList, class FaceList>
+void Foam::vtkWritePolyData::write
+(
+    const fileName& file,
+    const word& title,
+    const bool binary,
+    const PointField& points,
+    const VertexList& vertices,
+    const LineList& lines,
+    const FaceList& faces,
+    const wordList& fieldNames,
+    const boolList& fieldIsPointValues,
+    const UPtrList<const Field<label>>& fieldLabelValues
+    #define FieldTypeValuesConstArg(Type, nullArg) \
+        , const UPtrList<const Field<Type>>& field##Type##Values
+    FOR_ALL_FIELD_TYPES(FieldTypeValuesConstArg)
+    #undef FieldTypeValuesConstArg
+)
+{
+    // Open the file
+    std::ofstream os(file, std::ios::binary);
+
+    // Write the header
+    vtkWriteOps::writeHeader(os, binary, title);
+    os << "DATASET POLYDATA" << nl;
+
+    // Write the points
+    {
+        os  << "POINTS " << points.size() << " float" << nl;
+        List<floatScalar> coordinates(points.size()*3);
+        forAll(points, pointi)
+        {
+            const point& p = points[pointi];
+            forAll(p, i)
+            {
+                coordinates[3*pointi + i] = float(p[i]);
+            }
+        }
+        vtkWriteOps::write(os, binary, coordinates);
+    }
+
+    // Write the vertices
+    if (vertices.size())
+    {
+        os  << "VERTICES " << vertices.size() << ' '
+            << 2*vertices.size() << nl;
+        labelList data(2*vertices.size());
+        forAll(vertices, vertexi)
+        {
+            data[2*vertexi] = 1;
+            data[2*vertexi + 1] = vertices[vertexi];
+        }
+        vtkWriteOps::write(os, binary, data);
+    }
+
+    // Write the lines
+    if (lines.size())
+    {
+        label nLineNodes = 0;
+        forAll(lines, facei)
+        {
+            nLineNodes += lines[facei].size();
+        }
+        os  << "LINES " << lines.size() << ' '
+            << lines.size() + nLineNodes << nl;
+        labelList data(lines.size() + nLineNodes);
+        label i = 0;
+        forAll(lines, linei)
+        {
+            data[i ++] = lines[linei].size();
+            forAll(lines[linei], linePointi)
+            {
+                data[i ++] = lines[linei][linePointi];
+            }
+        }
+        vtkWriteOps::write(os, binary, data);
+    }
+
+    // Write the faces
+    if (faces.size())
+    {
+        label nFaceNodes = 0;
+        forAll(faces, facei)
+        {
+            nFaceNodes += faces[facei].size();
+        }
+        os  << "POLYGONS " << faces.size() << ' '
+            << faces.size() + nFaceNodes << nl;
+        labelList data(faces.size() + nFaceNodes);
+        label i = 0;
+        forAll(faces, facei)
+        {
+            data[i ++] = faces[facei].size();
+            forAll(faces[facei], facePointi)
+            {
+                data[i ++] = faces[facei][facePointi];
+            }
+        }
+        vtkWriteOps::write(os, binary, data);
+    }
+
+    // Write the fields
+    const label nPointFields = count(fieldIsPointValues, true);
+    const label nFaceFields = count(fieldIsPointValues, false);
+    if (nPointFields > 0)
+    {
+        os  << "POINT_DATA " << points.size() << nl
+            << "FIELD attributes " << nPointFields << nl;
+        writeFieldTypeValues<label, label>
+        (
+            os,
+            binary,
+            fieldNames,
+            fieldIsPointValues,
+            fieldLabelValues,
+            true
+        );
+        #define WriteFieldTypeValues(Type, nullArg)                            \
+            writeFieldTypeValues<Type, floatScalar>                            \
+            (                                                                  \
+                os,                                                            \
+                binary,                                                        \
+                fieldNames,                                                    \
+                fieldIsPointValues,                                            \
+                field##Type##Values,                                           \
+                true                                                           \
+            );
+        FOR_ALL_FIELD_TYPES(WriteFieldTypeValues)
+        #undef WriteFieldTypeValues
+    }
+    if (nFaceFields > 0)
+    {
+        os  << "CELL_DATA "
+            << vertices.size() + lines.size() + faces.size() << nl
+            << "FIELD attributes " << nFaceFields << nl;
+        writeFieldTypeValues<label, label>
+        (
+            os,
+            binary,
+            fieldNames,
+            fieldIsPointValues,
+            fieldLabelValues,
+            false
+        );
+        #define WriteFieldTypeValues(Type, nullArg)                            \
+            writeFieldTypeValues<Type, floatScalar>                            \
+            (                                                                  \
+                os,                                                            \
+                binary,                                                        \
+                fieldNames,                                                    \
+                fieldIsPointValues,                                            \
+                field##Type##Values,                                           \
+                false                                                          \
+            );
+        FOR_ALL_FIELD_TYPES(WriteFieldTypeValues)
+        #undef WriteFieldTypeValues
+    }
+}
+
+
+template
+<
+    class PointField,
+    class VertexList,
+    class LineList,
+    class FaceList,
+    class ... Args
+>
+inline void Foam::vtkWritePolyData::write
+(
+    const fileName& file,
+    const word& title,
+    const bool binary,
+    const PointField& points,
+    const VertexList& vertices,
+    const LineList& lines,
+    const FaceList& faces,
+    const Args& ... args
+)
+{
+    const label nFields = sizeof...(Args)/3;
+
+    wordList fieldNames(nFields);
+    boolList fieldIsPointValues(nFields);
+    UPtrList<const Field<label>> fieldLabelValues(nFields);
+    #define DeclareFieldTypeValues(Type, nullArg) \
+        UPtrList<const Field<Type>> field##Type##Values(nFields);
+    FOR_ALL_FIELD_TYPES(DeclareFieldTypeValues);
+    #undef DeclareFieldTypeValues
+
+    unpackFieldTypeValues
+    (
+        fieldNames,
+        fieldIsPointValues,
+        fieldLabelValues
+        #define FieldTypeValuesParameter(Type, nullArg) , field##Type##Values
+        FOR_ALL_FIELD_TYPES(FieldTypeValuesParameter),
+        #undef FieldTypeValuesParameter
+        args ...
+    );
+
+    write
+    (
+        file,
+        title,
+        binary,
+        points,
+        vertices,
+        lines,
+        faces,
+        fieldNames,
+        fieldIsPointValues,
+        fieldLabelValues
+        #define FieldTypeValuesParameter(Type, nullArg) , field##Type##Values
+        FOR_ALL_FIELD_TYPES(FieldTypeValuesParameter)
+        #undef FieldTypeValuesParameter
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
index 8f02d97f8b..a1d7d08c8e 100644
--- a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,7 +27,7 @@ License
 #include "OFstream.H"
 #include "OSspecific.H"
 #include "makeSurfaceWriterMethods.H"
-#include "vtkWriteOps.H"
+#include "vtkWritePolyData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -39,59 +39,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::vtkSurfaceWriter::writeGeometry
-(
-    std::ostream& os,
-    const pointField& points,
-    const faceList& faces
-) const
-{
-    const bool binary = (writeFormat_ == IOstream::BINARY);
-
-    // VTK header
-    vtkWriteOps::writeHeader(os, binary, "sampleSurface");
-    os << "DATASET POLYDATA" << nl;
-
-    // Write vertex coords
-    os  << "POINTS " << points.size() << " float" << nl;
-
-    List<floatScalar> po(points.size()*3);
-    label ind = 0;
-    forAll(points, pointi)
-    {
-        const point& pt = points[pointi];
-        forAll(pt, cmpt)
-        {
-            po[ind++] = float(pt[cmpt]);
-        }
-    }
-    vtkWriteOps::write(os, binary, po);
-
-    // Write faces
-    label nNodes = 0;
-    forAll(faces, facei)
-    {
-        nNodes += faces[facei].size();
-    }
-
-    os  << "POLYGONS " << faces.size() << ' '
-        << faces.size() + nNodes << nl;
-
-    labelList polygons(faces.size() + nNodes);
-    ind = 0;
-    forAll(faces, facei)
-    {
-        const face& f = faces[facei];
-        polygons[ind++] = f.size();
-        forAll(f, fp)
-        {
-            polygons[ind++] = f[fp];
-        }
-    }
-    vtkWriteOps::write(os, binary, polygons);
-}
-
-
 template<class Type>
 void Foam::vtkSurfaceWriter::Write
 (
@@ -104,53 +51,24 @@ void Foam::vtkSurfaceWriter::Write
     const bool isNodeValues
 ) const
 {
-    const bool binary = (writeFormat_ == IOstream::BINARY);
-
     if (!isDir(outputDir))
     {
         mkDir(outputDir);
     }
 
-    const word filePath = outputDir/fieldName + '_' + surfaceName + ".vtk";
-
-    ofstream os(filePath, std::ios::binary);
-
-    if (debug)
-    {
-        Info<< "Writing field " << fieldName << " to " << filePath << endl;
-    }
-
-    writeGeometry(os, points, faces);
-
-    // Write data
-    if (isNodeValues)
-    {
-        os  << "POINT_DATA ";
-    }
-    else
-    {
-        os  << "CELL_DATA ";
-    }
-
-    os  << values.size() << nl
-        << "FIELD attributes 1" << nl
-        << fieldName << " ";
-
-    const label nComp = pTraits<Type>::nComponents;
-
-    os  << nComp << " " << values.size() << " float" << nl;
-
-    List<floatScalar> vals(values.size()*nComp);
-    label ind = 0;
-    forAll(values, elemI)
-    {
-        for (direction cmpt=0; cmpt < nComp; ++cmpt)
-        {
-            vals[ind++] = component(values[elemI], cmpt);
-        }
-    }
-
-    vtkWriteOps::write(os, binary, vals);
+    vtkWritePolyData::write
+    (
+        outputDir/fieldName + '_' + surfaceName + ".vtk",
+        "sampleSurface",
+        writeFormat_ == IOstream::BINARY,
+        points,
+        labelList(),
+        edgeList(),
+        faces,
+        fieldName,
+        isNodeValues,
+        values
+    );
 }
 
 
@@ -186,15 +104,16 @@ void Foam::vtkSurfaceWriter::write
         mkDir(outputDir);
     }
 
-    word filePath =  outputDir/surfaceName + ".vtk";
-    ofstream os(filePath, std::ios::binary);
-
-    if (debug)
-    {
-        Info<< "Writing geometry to " << filePath << endl;
-    }
-
-    writeGeometry(os, points, faces);
+    vtkWritePolyData::write
+    (
+        outputDir/surfaceName + ".vtk",
+        "sampleSurface",
+        writeFormat_ == IOstream::BINARY,
+        points,
+        labelList(),
+        edgeList(),
+        faces
+    );
 }
 
 
diff --git a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
index a44fed6d49..514ca31075 100644
--- a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
+++ b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,13 +53,6 @@ class vtkSurfaceWriter
 {
     // Private Member Functions
 
-        void writeGeometry
-        (
-            std::ostream&,
-            const pointField&,
-            const faceList&
-        ) const;
-
         //- Templated write operation
         template<class Type>
         void Write
@@ -101,7 +94,6 @@ public:
             const faceList& faces
         ) const;
 
-
         //- Write scalarField for a single surface to file.
         //  One value per face or vertex (isNodeValues = true)
         virtual void write