diff --git a/applications/utilities/parallelProcessing/reconstructPar/Make/options b/applications/utilities/parallelProcessing/reconstructPar/Make/options
index e63099c8f5..ece7c8b030 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/Make/options
+++ b/applications/utilities/parallelProcessing/reconstructPar/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude
 
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index c4a15f97b6..7c82ef91a2 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -41,6 +41,10 @@ Description
 #include "pointFieldReconstructor.H"
 #include "reconstructLagrangian.H"
 
+#include "cellSet.H"
+#include "faceSet.H"
+#include "pointSet.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 bool haveAllTimes
@@ -99,6 +103,11 @@ int main(int argc, char *argv[])
         "skip reconstructing lagrangian positions and fields"
     );
     argList::addBoolOption
+    (
+        "sets",
+        "reconstruct cellSets, faceSets, pointSets"
+    );
+    argList::addBoolOption
     (
         "newTimes",
         "only reconstruct new times (i.e. that do not exist already)"
@@ -113,6 +122,9 @@ int main(int argc, char *argv[])
         args.optionLookup("fields")() >> selectedFields;
     }
 
+    const bool reconstructSets = args.optionFound("sets");
+
+
     const bool noLagrangian = args.optionFound("noLagrangian");
 
     HashSet<word> selectedLagrangianFields;
@@ -668,6 +680,148 @@ int main(int argc, char *argv[])
                     Info<< "No lagrangian fields" << nl << endl;
                 }
             }
+
+
+            if (reconstructSets)
+            {
+                // Scan to find all sets
+                HashTable<label> cSetNames;
+                HashTable<label> fSetNames;
+                HashTable<label> pSetNames;
+
+                forAll(procMeshes.meshes(), procI)
+                {
+                    const fvMesh& procMesh = procMeshes.meshes()[procI];
+
+                    IOobjectList objects
+                    (
+                        procMesh, procMesh.facesInstance(), "polyMesh/sets"
+                    );
+                    IOobjectList cSets(objects.lookupClass(cellSet::typeName));
+                    forAllConstIter(IOobjectList, cSets, iter)
+                    {
+                        cSetNames.insert(iter.key(), cSetNames.size());
+                    }
+
+                    IOobjectList fSets(objects.lookupClass(faceSet::typeName));
+                    forAllConstIter(IOobjectList, fSets, iter)
+                    {
+                        fSetNames.insert(iter.key(), fSetNames.size());
+                    }
+                    IOobjectList pSets(objects.lookupClass(pointSet::typeName));
+                    forAllConstIter(IOobjectList, pSets, iter)
+                    {
+                        pSetNames.insert(iter.key(), pSetNames.size());
+                    }
+                }
+
+                // Construct all sets
+                PtrList<cellSet> cellSets(cSetNames.size());
+                PtrList<faceSet> faceSets(fSetNames.size());
+                PtrList<pointSet> pointSets(pSetNames.size());
+
+                // Load sets
+                forAll(procMeshes.meshes(), procI)
+                {
+                    const fvMesh& procMesh = procMeshes.meshes()[procI];
+
+                    IOobjectList objects
+                    (
+                        procMesh, procMesh.facesInstance(), "polyMesh/sets"
+                    );
+
+                    // cellSets
+                    const labelList& cellMap =
+                        procMeshes.cellProcAddressing()[procI];
+
+                    IOobjectList cSets(objects.lookupClass(cellSet::typeName));
+                    forAllConstIter(IOobjectList, cSets, iter)
+                    {
+                        // Load cellSet
+                        const cellSet procSet(*iter());
+                        label setI = cSetNames[iter.key()];
+                        if (!cellSets.set(setI))
+                        {
+                            cellSets.set
+                            (
+                                setI,
+                                new cellSet(mesh, iter.key(), procSet.size())
+                            );
+                        }
+                        cellSet& cSet = cellSets[setI];
+
+                        forAllConstIter(cellSet, procSet, iter)
+                        {
+                            cSet.insert(cellMap[iter.key()]);
+                        }
+                    }
+
+                    // faceSets
+                    const labelList& faceMap =
+                        procMeshes.faceProcAddressing()[procI];
+
+                    IOobjectList fSets(objects.lookupClass(faceSet::typeName));
+                    forAllConstIter(IOobjectList, fSets, iter)
+                    {
+                        // Load faceSet
+                        const faceSet procSet(*iter());
+                        label setI = fSetNames[iter.key()];
+                        if (!faceSets.set(setI))
+                        {
+                            faceSets.set
+                            (
+                                setI,
+                                new faceSet(mesh, iter.key(), procSet.size())
+                            );
+                        }
+                        faceSet& fSet = faceSets[setI];
+
+                        forAllConstIter(faceSet, procSet, iter)
+                        {
+                            fSet.insert(mag(faceMap[iter.key()])-1);
+                        }
+                    }
+                    // pointSets
+                    const labelList& pointMap =
+                        procMeshes.pointProcAddressing()[procI];
+
+                    IOobjectList pSets(objects.lookupClass(pointSet::typeName));
+                    forAllConstIter(IOobjectList, pSets, iter)
+                    {
+                        // Load pointSet
+                        const pointSet propSet(*iter());
+                        label setI = pSetNames[iter.key()];
+                        if (!pointSets.set(setI))
+                        {
+                            pointSets.set
+                            (
+                                setI,
+                                new pointSet(mesh, iter.key(), propSet.size())
+                            );
+                        }
+                        pointSet& pSet = pointSets[setI];
+
+                        forAllConstIter(pointSet, propSet, iter)
+                        {
+                            pSet.insert(pointMap[iter.key()]);
+                        }
+                    }
+                }
+
+                // Write sets
+                forAll(cellSets, i)
+                {
+                    cellSets[i].write();
+                }
+                forAll(faceSets, i)
+                {
+                    faceSets[i].write();
+                }
+                forAll(pointSets, i)
+                {
+                    pointSets[i].write();
+                }
+            }
         }
     }