diff --git a/applications/test/decomposePar/Make/files b/applications/test/decomposePar/Make/files
deleted file mode 100644
index 39e4256963..0000000000
--- a/applications/test/decomposePar/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-Test-decomposePar.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-decomposePar
diff --git a/applications/test/decomposePar/Make/options b/applications/test/decomposePar/Make/options
deleted file mode 100644
index 25f37cdfad..0000000000
--- a/applications/test/decomposePar/Make/options
+++ /dev/null
@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -ldecompose \
-    -ldecompositionMethods \
-    -L$(FOAM_LIBBIN)/dummy \
-    -lkahipDecomp -lmetisDecomp -lscotchDecomp
diff --git a/applications/test/decomposePar/Test-decomposePar.C b/applications/test/decomposePar/Test-decomposePar.C
deleted file mode 100644
index d94e29df39..0000000000
--- a/applications/test/decomposePar/Test-decomposePar.C
+++ /dev/null
@@ -1,308 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    Test-decomposePar
-
-Group
-    grpParallelUtilities
-
-Description
-    Like decomposePar -dry-run, but with additional options
-
-\*---------------------------------------------------------------------------*/
-
-#include "OSspecific.H"
-#include "fvCFD.H"
-#include "cpuTime.H"
-#include "volFields.H"
-#include "IOobjectList.H"
-#include "regionProperties.H"
-#include "decompositionInformation.H"
-#include "decompositionModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addNote
-    (
-        "Special-purpose version of decomposePar with additional"
-        " -domain and -method options."
-        " The '-dry-run' and '-cellDist' are implicit.\n"
-        "NB: The -domain/-method overrides may not work very well with regions"
-    );
-
-    argList::noParallel();
-    argList::addOption
-    (
-        "decomposeParDict",
-        "file",
-        "Use specified file for decomposePar dictionary"
-    );
-
-    #include "addAllRegionOptions.H"
-
-    argList::addBoolOption
-    (
-        "verbose",
-        "Additional verbosity"
-    );
-    argList::addOption
-    (
-        "domains",
-        "N",
-        "Override numberOfSubdomains"
-    );
-
-    argList::addOption
-    (
-        "method",
-        "name",
-        "Override decomposition method"
-    );
-
-
-    // These are implicit so just ignore them
-    argList::ignoreOptionCompat({"dry-run", 0}, false);
-    argList::addBoolOption
-    (
-        "cellDist",
-        "Write cell distribution as volScalarField for visualization"
-    );
-    argList::addBoolOption
-    (
-        "cellDist-internal",
-        "Write cell distribution (internal field) for visualization"
-    );
-
-    // Include explicit constant options, have zero from time range
-    timeSelector::addOptions(true, false);
-
-    #include "setRootCase.H"
-
-    const bool optRegion  = args.found("region");
-    const bool verbose    = args.found("verbose");
-
-    const label numSubdomains = args.getOrDefault<label>("domains", 0);
-    const word methodName = args.getOrDefault<word>("method", word::null);
-
-    // Set time from database
-    #include "createTime.H"
-    // Allow override of time
-    instantList times = timeSelector::selectIfPresent(runTime, args);
-
-    // Allow override of decomposeParDict location
-    const fileName decompDictFile =
-        args.getOrDefault<fileName>("decomposeParDict", "");
-
-    // Get region names
-    #include "getAllRegionOptions.H"
-
-    Info<< "Decomposing regions: "
-        << flatOutput(regionNames) << nl << endl;
-
-
-    forAll(regionNames, regioni)
-    {
-        const word& regionName = regionNames[regioni];
-        const word& regionDir =
-        (
-            regionName == polyMesh::defaultRegion ? word::null : regionName
-        );
-
-        Info<< "\n\nDecomposing mesh " << regionName << nl << endl;
-        Info<< "Create mesh..." << flush;
-
-        fvMesh mesh
-        (
-            IOobject
-            (
-                regionName,
-                runTime.timeName(),
-                runTime,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        );
-
-        Info<< " nCells = " << mesh.nCells() << endl;
-
-        Info<< "\nCalculating distribution of cells" << endl;
-        cpuTime decompositionTime;
-
-        const decompositionModel& model = decompositionModel::New
-        (
-            mesh,
-            decompDictFile
-        );
-
-        // Allow command-line override for quick testing
-
-        dictionary& modelDict = const_cast<decompositionModel&>(model);
-
-        if (numSubdomains)
-        {
-            modelDict.add
-            (
-                word("numberOfSubdomains"),
-                numSubdomains,
-                true
-            );
-        }
-
-        if (!methodName.empty())
-        {
-            modelDict.add
-            (
-                word("method"),
-                methodName,
-                true
-            );
-        }
-
-        scalarField cellWeights;
-        word weightName;
-        if (model.readIfPresent("weightField", weightName))
-        {
-            volScalarField weights
-            (
-                IOobject
-                (
-                    weightName,
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh
-            );
-            cellWeights = weights.primitiveField();
-        }
-
-        decompositionMethod& method = model.decomposer();
-
-        CompactListList<label> cellCells;
-        decompositionMethod::calcCellCells
-        (
-            mesh,
-            identity(mesh.nCells()),
-            mesh.nCells(),
-            false,
-            cellCells
-        );
-
-        labelList cellToProc = method.decompose(mesh, cellWeights);
-
-        Info<< "\nFinished decomposition into "
-            << method.nDomains() << " domains in "
-            << decompositionTime.elapsedCpuTime()
-            << " s" << nl << endl;
-
-        decompositionInformation info
-        (
-            cellCells,
-            cellToProc,
-            method.nDomains()
-        );
-
-        if (verbose)
-        {
-            info.printDetails(Info);
-            Info<< nl;
-        }
-        info.printSummary(Info);
-
-
-        if (args.found("cellDist-internal"))
-        {
-            volScalarField::Internal cellDist
-            (
-                IOobject
-                (
-                    "cellDist",
-                    runTime.timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                mesh,
-                dimensionedScalar("cellDist", dimless, -1)
-            );
-
-            forAll(cellToProc, celli)
-            {
-                cellDist[celli] = cellToProc[celli];
-            }
-
-            cellDist.write();
-
-            Info<< nl << "Wrote decomposition as volScalarField::Internal to "
-                << cellDist.name() << " for visualization."
-                << endl;
-
-            fileHandler().flush();
-        }
-        else if (args.found("cellDist"))
-        {
-            volScalarField cellDist
-            (
-                IOobject
-                (
-                    "cellDist",
-                    runTime.timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
-                ),
-                mesh,
-                dimensionedScalar("cellDist", dimless, -1),
-                zeroGradientFvPatchScalarField::typeName
-            );
-
-            forAll(cellToProc, celli)
-            {
-                cellDist[celli] = cellToProc[celli];
-            }
-
-            cellDist.correctBoundaryConditions();
-            cellDist.write();
-
-            Info<< nl << "Wrote decomposition as volScalarField to "
-                << cellDist.name() << " for visualization."
-                << endl;
-
-            fileHandler().flush();
-        }
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/files b/applications/utilities/parallelProcessing/decomposePar/Make/files
index 0aba0bb703..03abb7247f 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/files
+++ b/applications/utilities/parallelProcessing/decomposePar/Make/files
@@ -2,7 +2,11 @@ decomposePar.C
 domainDecomposition.C
 domainDecompositionMesh.C
 domainDecompositionDistribute.C
+domainDecompositionWrite.C
+
 domainDecompositionDryRun.C
+domainDecompositionDryRunWrite.C
+
 dimFieldDecomposer.C
 pointFieldDecomposer.C
 faMeshDecomposition.C
diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/options b/applications/utilities/parallelProcessing/decomposePar/Make/options
index efe085933a..104ca56d7c 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/options
+++ b/applications/utilities/parallelProcessing/decomposePar/Make/options
@@ -1,20 +1,22 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
     -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
     -llagrangian \
     -ldynamicMesh \
     -lgenericPatchFields \
-    -ldecompose \
     -ldecompositionMethods \
+    -ldecompose \
     -L$(FOAM_LIBBIN)/dummy \
     -lkahipDecomp -lmetisDecomp -lscotchDecomp
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index daf25f3894..7c9784ce8d 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -47,7 +47,7 @@ Usage
 
       - \par -cellDist
         Write the cell distribution as a labelList, for use with 'manual'
-        decomposition method and as a volScalarField for visualization.
+        decomposition method and as a VTK or volScalarField for visualization.
 
       - \par -constant
         Include the 'constant/' dir in the times list.
@@ -67,7 +67,7 @@ Usage
 
       - \par -dry-run
         Test without writing the decomposition. Changes -cellDist to
-        only write volScalarField.
+        only write VTK output.
 
       - \par -fields
         Use existing geometry decomposition and convert fields only.
@@ -301,13 +301,28 @@ int main(int argc, char *argv[])
     (
         "dry-run",
         "Test without writing the decomposition. "
-        "Changes -cellDist to only write volScalarField."
+        "Changes -cellDist to only write VTK output."
     );
     argList::addBoolOption
     (
         "verbose",
         "Additional verbosity"
     );
+    argList::addOption
+    (
+        "domains",
+        "N",
+        "Override numberOfSubdomains (-dry-run only)",
+        true  // Advanced option
+    );
+    argList::addOption
+    (
+        "method",
+        "name",
+        "Override decomposition method (-dry-run only)",
+        true  // Advanced option
+    );
+
     argList::addBoolOption
     (
         "cellDist",
@@ -421,7 +436,9 @@ int main(int argc, char *argv[])
                     IOobject::NO_WRITE,
                     false
                 ),
-                decompDictFile
+                decompDictFile,
+                args.getOrDefault<label>("domains", 0),
+                args.getOrDefault<word>("method", word::null)
             );
 
             decompTest.execute(writeCellDist, verbose);
@@ -580,33 +597,9 @@ int main(int argc, char *argv[])
             {
                 const labelList& procIds = mesh.cellToProc();
 
-                // Write decomposition as volScalarField for visualization
-                volScalarField cellDist
-                (
-                    IOobject
-                    (
-                        "cellDist",
-                        runTime.timeName(),
-                        mesh,
-                        IOobject::NO_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    mesh,
-                    dimensionedScalar("cellDist", dimless, -1),
-                    zeroGradientFvPatchScalarField::typeName
-                );
-
-                forAll(procIds, celli)
-                {
-                   cellDist[celli] = procIds[celli];
-                }
-
-                cellDist.correctBoundaryConditions();
-                cellDist.write();
-
-                Info<< nl << "Wrote decomposition as volScalarField to "
-                    << cellDist.name() << " for visualization."
-                    << endl;
+                // Write decomposition for visualization
+                mesh.writeVolField("cellDist");
+                //TBD: mesh.writeVTK("cellDist");
 
                 // Write decomposition as labelList for use with 'manual'
                 // decomposition method.
@@ -626,7 +619,7 @@ int main(int argc, char *argv[])
                 cellDecomposition.write();
 
                 Info<< nl << "Wrote decomposition to "
-                    << cellDecomposition.objectPath()
+                    << runTime.relativePath(cellDecomposition.objectPath())
                     << " for use in manual decomposition." << endl;
             }
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H
index 6a21b59dbb..af1f509fb6 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,8 @@ Description
 
 SourceFiles
     domainDecomposition.C
-    decomposeMesh.C
+    domainDecompositionDistribute.C
+    domainDecompositionMesh.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -40,16 +42,16 @@ SourceFiles
 
 #include "fvMesh.H"
 #include "labelList.H"
-#include "SLList.H"
-#include "PtrList.H"
 #include "point.H"
 #include "Time.H"
 #include "volFields.H"
 
-
 namespace Foam
 {
 
+// Forward Declarations
+class decompositionModel;
+
 /*---------------------------------------------------------------------------*\
                      Class domainDecomposition Declaration
 \*---------------------------------------------------------------------------*/
@@ -58,7 +60,7 @@ class domainDecomposition
 :
     public fvMesh
 {
-    // Private data
+    // Private Data
 
         //- Optional: points at the facesInstance
         autoPtr<pointIOField> facesInstancePointsPtr_;
@@ -117,6 +119,7 @@ class domainDecomposition
         //- Sub patch sizes for inter-processor patches
         List<labelListList> procProcessorPatchSubPatchStarts_;
 
+
     // Private Member Functions
 
         void distributeCells();
@@ -165,7 +168,7 @@ public:
         //  \param io the IOobject for mesh
         //  \param decompDictFile optional non-standard location for the
         //      decomposeParDict file
-        domainDecomposition
+        explicit domainDecomposition
         (
             const IOobject& io,
             const fileName& decompDictFile = ""
@@ -178,29 +181,68 @@ public:
 
     // Member Functions
 
+        //- The mesh
+        const fvMesh& mesh() const noexcept
+        {
+            return *this;
+        }
+
+
+    // Settings
+
         //- Number of processor in decomposition
-        label nProcs() const
+        label nProcs() const noexcept
         {
             return nProcs_;
         }
 
-        //- Is the decomposition data to be distributed for each processor
-        bool distributed() const
+        //- Is decomposition data to be distributed for each processor
+        bool distributed() const noexcept
         {
             return distributed_;
         }
 
+        //- Change distributed flag
+        bool distributed(const bool on) noexcept
+        {
+            bool old(distributed_);
+            distributed_ = on;
+            return old;
+        }
+
+        //- Return decomposition model used
+        const decompositionModel& model() const;
+
+        //- Update flags based on the decomposition model settings
+        //  Sets "distributed"
+        void updateParameters(const dictionary& params);
+
+
+    // Mappings
+
+        //- Cell-processor decomposition labels
+        const labelList& cellToProc() const noexcept
+        {
+            return cellToProc_;
+        }
+
+
+    // Decompose
+
         //- Decompose mesh
         void decomposeMesh();
 
         //- Write decomposition
         bool writeDecomposition(const bool decomposeSets);
 
-        //- Cell-processor decomposition labels
-        const labelList& cellToProc() const
-        {
-            return cellToProc_;
-        }
+
+    // Write
+
+        //- Write decomposition as volScalarField for visualization
+        void writeVolField(const word& timeName) const;
+
+        //- Write cell distribution as VTU file (binary)
+        void writeVTK(const fileName& file) const;
 };
 
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
index 00174b9879..b175c96e4f 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 #include "volFields.H"
 #include "decompositionModel.H"
 #include "decompositionInformation.H"
-#include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -137,32 +136,9 @@ void Foam::domainDecompositionDryRun::execute
 
     if (writeCellDist)
     {
-        // Write decomposition as volScalarField for visualization
-        volScalarField cellDist
-        (
-            IOobject
-            (
-                "cellDist",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar("cellDist", dimless, -1),
-            zeroGradientFvPatchScalarField::typeName
-        );
-
-        forAll(cellToProc, celli)
-        {
-            cellDist[celli] = cellToProc[celli];
-        }
-
-        cellDist.correctBoundaryConditions();
-        cellDist.write();
-
-        Info<< "Wrote decomposition as volScalarField for visualization:"
-            << nl << cellDist.objectPath() << nl;
+        // Write decomposition for visualization
+        // - write as VTU to avoid any impact
+        writeVTK("cellDist", cellToProc);
 
 // Less likely that this is actually required, but may be useful...
 //
@@ -183,9 +159,9 @@ void Foam::domainDecompositionDryRun::execute
 //        );
 //        cellDecomposition.write();
 //
-//        Info<< "Wrote decomposition for use in manual decomposition:"
-//            << nl << cellDecomposition.objectPath() << nl;
-// #endif
+//        Info<< nl << "Wrote decomposition to "
+//            << runTime.relativePath(cellDecomposition.objectPath())
+//            << " for use in manual decomposition." << endl;
 
         Info<< nl;
     }
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.H b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.H
index cedab0da4d..7764b233af 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.H
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRun.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -63,6 +63,23 @@ class domainDecompositionDryRun
         const word methodNameOverride_;
 
 
+    // Private Members
+
+        //- Write decomposition as volScalarField for visualization
+        void writeVolField
+        (
+            const word& timeName,
+            const labelUList& cellToProc
+        ) const;
+
+        //- Write cell distribution as VTU file (binary)
+        void writeVTK
+        (
+            const fileName& file,
+            const labelUList& cellToProc
+        ) const;
+
+
 public:
 
     // Constructors
@@ -88,6 +105,12 @@ public:
 
     // Member Functions
 
+        //- The mesh
+        const fvMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
+
         //- Perform dry-run
         void execute(const bool writeCellDist, const bool verbose = false);
 };
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C
new file mode 100644
index 0000000000..5e02b54f57
--- /dev/null
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "domainDecompositionDryRun.H"
+#include "foamVtkInternalMeshWriter.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::domainDecompositionDryRun::writeVolField
+(
+    const word& timeName,
+    const labelUList& procIds
+) const
+{
+    // Write decomposition as volScalarField for visualization
+    volScalarField cellDist
+    (
+        IOobject
+        (
+            "cellDist",
+            timeName,
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh(),
+        dimensionedScalar("cellDist", dimless, -1),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    forAll(procIds, celli)
+    {
+        cellDist[celli] = procIds[celli];
+    }
+
+    cellDist.correctBoundaryConditions();
+    cellDist.write();
+
+    Info<< nl << "Wrote decomposition to "
+        << this->mesh().time().relativePath(cellDist.objectPath())
+        << " (volScalarField) for visualization."
+        << endl;
+}
+
+
+void Foam::domainDecompositionDryRun::writeVTK
+(
+    const fileName& file,
+    const labelUList& cellToProc
+) const
+{
+    const vtk::vtuCells cells(this->mesh());
+
+    // not parallel
+    vtk::internalMeshWriter writer(this->mesh(), cells, file, false);
+
+    writer.writeGeometry();
+    writer.beginCellData();
+    writer.writeCellData("procID", cellToProc);
+
+    Info<< "Wrote decomposition to "
+        << this->mesh().time().relativePath(writer.output())
+        << " for visualization."
+        << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionWrite.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionWrite.C
new file mode 100644
index 0000000000..cc2ee88bcd
--- /dev/null
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionWrite.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "domainDecomposition.H"
+#include "foamVtkInternalMeshWriter.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::domainDecomposition::writeVolField
+(
+    const word& timeName
+) const
+{
+    const labelList& procIds = this->cellToProc();
+
+    // Write decomposition as volScalarField for visualization
+    volScalarField cellDist
+    (
+        IOobject
+        (
+            "cellDist",
+            timeName,
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh(),
+        dimensionedScalar("cellDist", dimless, -1),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    forAll(procIds, celli)
+    {
+        cellDist[celli] = procIds[celli];
+    }
+
+    cellDist.correctBoundaryConditions();
+    cellDist.write();
+
+    Info<< nl << "Wrote decomposition to "
+        << this->mesh().time().relativePath(cellDist.objectPath())
+        << " (volScalarField) for visualization."
+        << endl;
+}
+
+
+void Foam::domainDecomposition::writeVTK(const fileName& file) const
+{
+    const vtk::vtuCells cells(this->mesh());
+
+    // not parallel
+    vtk::internalMeshWriter writer(this->mesh(), cells, file, false);
+
+    writer.writeGeometry();
+    writer.beginCellData();
+    writer.writeCellData("procID", cellToProc());
+
+    Info<< "Wrote decomposition to "
+        << this->mesh().time().relativePath(writer.output())
+        << " for visualization."
+        << endl;
+}
+
+
+// ************************************************************************* //