diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index f15c61b36315f0d3bd89dcdb14334c7bbc6d7e70..39568568524f490ea0922ba49771d3358b106d4d 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -31,6 +31,7 @@ License
 #include "volPointInterpolation.H"
 #include "PatchTools.H"
 #include "mapPolyMesh.H"
+#include "sampledTriSurfaceMesh.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -84,6 +85,37 @@ void Foam::sampledSurfaces::writeGeometry() const
 }
 
 
+void Foam::sampledSurfaces::writeOriginalIds()
+{
+    const word fieldName = "Ids";
+    const fileName outputDir = outputPath_/time_.timeName();
+
+    forAll(*this, surfI)
+    {
+        const sampledSurface& s = operator[](surfI);
+
+        if (isA<sampledTriSurfaceMesh>(s))
+        {
+            const sampledTriSurfaceMesh& surf =
+                dynamicCast<const sampledTriSurfaceMesh&>(s);
+
+            if (surf.keepIds())
+            {
+                const labelList& idLst = surf.originalIds();
+
+                Field<scalar> ids(idLst.size());
+                forAll(idLst, i)
+                {
+                    ids[i] = idLst[i];
+                }
+
+                writeSurface(ids, surfI, fieldName, outputDir);
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledSurfaces::sampledSurfaces
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index 1d9c9e730b91a668999fa9b5292804f22c0d2128..8d610dbb5921504452b1c3076d480a98073a1fd6 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -116,6 +116,9 @@ class sampledSurfaces
         //- Write geometry only
         void writeGeometry() const;
 
+        //- Write scalar field with original ids
+        void writeOriginalIds();
+
         //- Write sampled fieldName on surface and on outputDir path
         template<class Type>
         void writeSurface
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
index 4e25d84af06d86b2f616e092d5afd9451760c021..143837a4eb3294936e916089e57103e155b95bd8 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
@@ -60,7 +60,6 @@ Foam::label Foam::sampledSurfaces::classifyFields()
     {
         // Check currently available fields
         wordList allFields = obr_.sortedNames();
-        labelList indices = findStrings(fieldSelection_, allFields);
 
         forAll(fieldSelection_, i)
         {
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
index ba39848aad2bd0d809c3e56b2d0a0e615f43b679..0dd8c248efdff0e0c493171f84d401da54d5e406 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
@@ -49,9 +49,7 @@ void Foam::sampledSurfaces::writeSurface
         gatheredValues[Pstream::myProcNo()] = values;
         Pstream::gatherList(gatheredValues);
 
-
         fileName sampleFile;
-
         if (Pstream::master())
         {
             // Combine values into single field
@@ -181,24 +179,21 @@ void Foam::sampledSurfaces::sampleAndWrite
 template<class GeoField>
 void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& objects)
 {
-    wordList names;
+    wordList fieldNames;
     if (loadFromFiles_)
     {
-        IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
-        names = fieldObjects.names();
+        fieldNames = objects.sortedNames(GeoField::typeName, fieldSelection_);
     }
     else
     {
-        names = mesh_.thisDb().names<GeoField>();
-    }
+        fieldNames = mesh_.thisDb().sortedNames<GeoField>(fieldSelection_);
 
-    labelList nameIDs(findStrings(fieldSelection_, names));
-
-    wordHashSet fieldNames(wordList(names, nameIDs));
+        writeOriginalIds();
+    }
 
-    forAllConstIter(wordHashSet, fieldNames, iter)
+    forAll(fieldNames, fieldi)
     {
-        const word& fieldName = iter.key();
+        const word& fieldName = fieldNames[fieldi];
 
         if ((Pstream::master()) && verbose_)
         {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index c1c08006d4542c0dc3caf1f5e41d78bb9d23f851..38fdc7c748f8d86602db4f299e1cb7c80fb1cc28 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -296,6 +296,11 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
     }
 
 
+    if (keepIds_)
+    {
+        originalIds_ = faceMap;
+    }
+
     // Subset cellOrFaceLabels
     cellOrFaceLabels = UIndirectList<label>(cellOrFaceLabels, faceMap)();
 
@@ -545,6 +550,8 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     ),
     sampleSource_(sampleSource),
     needsUpdate_(true),
+    keepIds_(false),
+    originalIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -573,6 +580,8 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     ),
     sampleSource_(samplingSourceNames_[dict.lookup("source")]),
     needsUpdate_(true),
+    keepIds_(dict.lookupOrDefault<Switch>("keepIds", false)),
+    originalIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -594,7 +603,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
             name,
             mesh.time().constant(), // instance
             "triSurface",           // local
-            mesh,                  // registry
+            mesh,                   // registry
             IOobject::NO_READ,
             IOobject::NO_WRITE,
             false
@@ -603,6 +612,8 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     ),
     sampleSource_(samplingSourceNames_[sampleSourceName]),
     needsUpdate_(true),
+    keepIds_(false),
+    originalIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -633,6 +644,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
     sampledSurface::clearGeom();
     MeshStorage::clear();
 
+    originalIds_.clear();
     boundaryTreePtr_.clear();
     sampleElements_.clear();
     samplePoints_.clear();
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index 7403bfc65e24a4888b7b1fb18282646e3d1d4b3d..6cd96a4770c0c49e8c19a7a871016fd3c6792e47 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -117,6 +117,13 @@ private:
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
+        //- Retain element ids/order of original surface
+        bool keepIds_;
+
+        //- List of element ids/order of the original surface,
+        //  when keepIds is active.
+        labelList originalIds_;
+
         //- Search tree for all non-coupled boundary faces
         mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
 
@@ -240,6 +247,19 @@ public:
             return MeshStorage::Cf();
         }
 
+        //- If element ids/order of the original surface are kept
+        bool keepIds() const
+        {
+            return keepIds_;
+        }
+
+        //- List of element ids/order of the original surface,
+        //  when keepIds is active.
+        const labelList& originalIds() const
+        {
+            return originalIds_;
+        }
+
 
         //- Sample field on surface
         virtual tmp<scalarField> sample