diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
index 3198116b9f0e2988f550eb19ff84025b639111b2..9d907d286ab19faa72292ad9cd567df1b972e134 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
@@ -118,7 +118,7 @@ inline static IOobject selectReadIO(const word& name, const Time& runTime)
 
 void Foam::sampledMeshedSurface::setZoneMap()
 {
-    // Ensure zoneIds_ are correctly populated
+    // Populate zoneIds_ based on surfZone information
 
     const meshedSurface& s = static_cast<const meshedSurface&>(*this);
 
@@ -279,19 +279,10 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
     s = surface_.subsetMesh(facesToSubset, pointMap, faceMap);
 
+
     // Ensure zoneIds_ are indeed correct
     setZoneMap();
 
-    // This is currently only partially useful
-    if (keepIds_)
-    {
-        originalIds_ = faceMap;
-    }
-    else
-    {
-        originalIds_.clear();
-    }
-
 
     // Subset cellOrFaceLabels (for compact faces)
     cellOrFaceLabels = labelUIndList(cellOrFaceLabels, faceMap)();
@@ -492,8 +483,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     ),
     sampleSource_(sampleSource),
     needsUpdate_(true),
-    keepIds_(false),
-    originalIds_(),
+    keepIds_(true),
     zoneIds_(),
     sampleElements_(),
     samplePoints_()
@@ -524,8 +514,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface
     ),
     sampleSource_(samplingSourceNames_.get("source", dict)),
     needsUpdate_(true),
-    keepIds_(dict.lookupOrDefault("keepIds", false)),
-    originalIds_(),
+    keepIds_(dict.getOrDefault("keepIds", true)),
     zoneIds_(),
     sampleElements_(),
     samplePoints_()
@@ -610,7 +599,6 @@ bool Foam::sampledMeshedSurface::expire()
     MeshStorage::clear();
     zoneIds_.clear();
 
-    originalIds_.clear();
     sampleElements_.clear();
     samplePoints_.clear();
 
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
index db0c18b3168f58937eecb79cedead0b0c88d7c3a..bc4bdfb094d630e60d2d0834b8dae221d26ce282 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
@@ -86,7 +86,7 @@ Usage
         surface  | surface name in triSurface/              | yes |
         patches  | Limit to named surface regions (wordRes) | no  |
         source   | cells/insideCells/boundaryFaces          | yes |
-        keepIds  | pass through id numbering                | no  | false
+        keepIds  | pass through id numbering                | no  | true
         file     | Alternative file name                    | no  |
         fileType | The surface format                       | no  | (extension)
         scale    | Surface scaling factor                   | no  | 0
@@ -156,10 +156,6 @@ private:
         //- Retain element ids/order of original surface
         bool keepIds_;
 
-        //- List of element ids/order of the original surface,
-        //- when keepIds is active.
-        labelList originalIds_;
-
         //- For compatibility with the meshSurf interface
         labelList zoneIds_;
 
@@ -258,6 +254,12 @@ public:
             return zoneIds_;
         }
 
+        //- Per-face identifier (eg, element Id)
+        virtual const labelList& faceIds() const
+        {
+            return MeshStorage::faceIds();
+        }
+
         //- Face area vectors
         virtual const vectorField& Sf() const
         {
@@ -279,14 +281,7 @@ public:
         //- If element ids/order of the original surface are kept
         virtual bool hasFaceIds() const
         {
-            return keepIds_;
-        }
-
-        //- List of element ids/order of the original surface,
-        //  when keepIds is active.
-        virtual const labelList& originalIds() const
-        {
-            return originalIds_;
+            return keepIds_ && !MeshStorage::faceIds().empty();
         }
 
         //- Sampling boundary values instead of cell values
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index 6dc8d28e69a16ce827dd484e4d779e208a637bd4..b21df8349e1c51d603f394750f3b7b861441f565 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -340,13 +340,6 @@ public:
             return false;
         }
 
-        //- List of element ids/order of the original surface,
-        //- when hasFaceIds is true.
-        virtual const labelList& originalIds() const
-        {
-            return labelList::null();
-        }
-
 
     // General registry storage (optional)
 
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index 0f133b402c0a732d2dc56be2c8861d0e3970fac4..a87a22d4801ddb50eecea74cd416dce5e6f200a6 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -558,10 +558,15 @@ bool Foam::sampledSurfaces::performAction(unsigned request)
             // Write original ids
             if (s.hasFaceIds() && !s.interpolate())
             {
+                // This is still quite horrible.
+
+                Field<label> ids(s.faceIds());
+                ids += label(1);  // From 0-based to 1-based
+
                 writeSurface
                 (
                     outWriter,
-                    Field<label>(s.originalIds()),
+                    ids,
                     "Ids"
                 );
             }
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index a5dd2fe7a5fd55cf6f760acbd91a6c3b8987f1aa..641a7bf2bc36d101b7cd5e86fd2a0e4aa90c4ee2 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -69,7 +69,7 @@ Description
         }
 
         surfaces
-        (
+        {
             f0surf
             {
                 type        meshedSurface;
@@ -88,7 +88,7 @@ Description
                 // Optional: registry storage
                 store       true
             }
-        );
+        }
     }
     \endverbatim
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 11194b8144cae02d65498c186a62cf3e8aa38493..991664f2b5c86451ad5763f068f1c2eb64fa6d88 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -183,6 +183,7 @@ template<class Face>
 Foam::MeshedSurface<Face>::MeshedSurface()
 :
     MeshReference(List<Face>(), pointField()),
+    faceIds_(),
     zones_()
 {}
 
@@ -194,6 +195,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     MeshReference(surf.surfFaces(), surf.points()),
+    faceIds_(surf.faceIds_),
     zones_(surf.zones_)
 {}
 
@@ -204,20 +206,40 @@ Foam::MeshedSurface<Face>::MeshedSurface
     const UnsortedMeshedSurface<Face>& surf
 )
 :
-    MeshReference(List<Face>(), surf.points())  // Copy points only
+    MeshReference(List<Face>(), surf.points()),  // Copy points only
+    faceIds_(),
+    zones_()
 {
     labelList faceMap;
     this->storedZones() = surf.sortedZones(faceMap);
 
+    // Faces, in the sorted order
     const List<Face>& origFaces = surf;
-    List<Face> newFaces(origFaces.size());
-
-    forAll(newFaces, facei)
     {
-        newFaces[faceMap[facei]] = origFaces[facei];
+        List<Face> newFaces(origFaces.size());
+
+        forAll(newFaces, facei)
+        {
+            newFaces[faceMap[facei]] = origFaces[facei];
+        }
+
+        this->storedFaces().transfer(newFaces);
     }
 
-    this->storedFaces().transfer(newFaces);
+    // FaceIds, in the sorted order
+    const labelList& origIds = surf.faceIds();
+
+    if (origIds.size() == origFaces.size())
+    {
+        labelList newFaceIds(origIds.size());
+
+        forAll(newFaceIds, facei)
+        {
+            newFaceIds[faceMap[facei]] = origIds[facei];
+        }
+
+        this->storedFaceIds().transfer(newFaceIds);
+    }
 }
 
 
@@ -254,6 +276,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     MeshReference(faceLst, pointLst), // Copy construct
+    faceIds_(),
     zones_(zoneLst)
 {}
 
@@ -267,6 +290,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     MeshReference(faceLst, pointLst, true), // Move construct
+    faceIds_(),
     zones_(zoneLst)
 {}
 
@@ -281,6 +305,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     MeshReference(faceLst, pointLst), // Copy construct
+    faceIds_(),
     zones_()
 {
     if (zoneSizes.size())
@@ -307,6 +332,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     MeshReference(faceLst, pointLst, true), // Move construct
+    faceIds_(),
     zones_()
 {
     if (zoneSizes.size())
@@ -571,6 +597,7 @@ void Foam::MeshedSurface<Face>::clear()
 
     storedPoints().clear();
     storedFaces().clear();
+    storedFaceIds().clear();
     storedZones().clear();
 }
 
@@ -690,6 +717,17 @@ bool Foam::MeshedSurface<Face>::stitchFaces
         faceMap.resize(newFacei);
         faceLst.resize(newFacei);
 
+        // The faceMap is a newToOld mapping and only removes elements
+        if (faceIds_.size())
+        {
+            forAll(faceMap, facei)
+            {
+                faceIds_[facei] = faceIds_[faceMap[facei]];
+            }
+
+            faceIds_.resize(newFacei);
+        }
+
         remapFaces(faceMap);
     }
     faceMap.clear();
@@ -849,6 +887,17 @@ bool Foam::MeshedSurface<Face>::checkFaces
         faceMap.resize(newFacei);
         faceLst.resize(newFacei);
 
+        // The faceMap is a newToOld mapping and only removes elements
+        if (faceIds_.size())
+        {
+            forAll(faceMap, facei)
+            {
+                faceIds_[facei] = faceIds_[faceMap[facei]];
+            }
+
+            faceIds_.resize(newFacei);
+        }
+
         remapFaces(faceMap);
     }
     faceMap.clear();
@@ -979,6 +1028,8 @@ Foam::label Foam::MeshedSurface<Face>::triangulate
         return 0;
     }
 
+    this->storedFaceIds().clear();  // Invalid or misleading
+
     List<Face>  newFaces(nTri);
     List<label> faceMap;
 
@@ -1118,11 +1169,20 @@ Foam::MeshedSurface<Face>::subsetMeshImpl
         zone.size() = newFacei - zone.start();
     }
 
+
+    // Subset of faceIds. Can be empty.
+    labelList newFaceIds;
+    if (faceIds_.size())
+    {
+        newFaceIds = labelUIndList(faceIds_, faceMap);
+    }
+
     // Construct the sub-surface
     MeshedSurface<Face> newSurf;
     newSurf.storedFaces().transfer(newFaces);
     newSurf.storedPoints().transfer(newPoints);
     newSurf.storedZones().transfer(newZones);
+    newSurf.storedFaceIds().transfer(newFaceIds);
 
     return newSurf;
 }
@@ -1224,6 +1284,7 @@ void Foam::MeshedSurface<Face>::swap
     this->storedPoints().swap(surf.storedPoints());
     this->storedFaces().swap(surf.storedFaces());
     this->storedZones().swap(surf.storedZones());
+    this->storedFaceIds().swap(surf.storedFaceIds());
 }
 
 
@@ -1239,6 +1300,7 @@ void Foam::MeshedSurface<Face>::transfer
     this->storedPoints().transfer(pointLst);
     this->storedFaces().transfer(faceLst);
     this->storedZones().clear();
+    this->storedFaceIds().clear();  // Likely to be invalid
 }
 
 
@@ -1258,6 +1320,7 @@ void Foam::MeshedSurface<Face>::transfer
     this->storedPoints().transfer(surf.storedPoints());
     this->storedFaces().transfer(surf.storedFaces());
     this->storedZones().transfer(surf.storedZones());
+    this->storedFaceIds().transfer(surf.storedFaceIds());
 
     surf.clear();
 }
@@ -1316,6 +1379,8 @@ void Foam::MeshedSurface<Face>::swapFaces(List<Face>& faces)
 {
     MeshReference::clearOut();  // Topology changes
 
+    this->storedFaceIds().clear();  // Likely to be invalid
+
     this->storedFaces().swap(faces);
 }
 
@@ -1401,6 +1466,7 @@ void Foam::MeshedSurface<Face>::operator=(const MeshedSurface<Face>& surf)
 
     this->storedPoints() = surf.points();
     this->storedFaces()  = surf.surfFaces();
+    this->storedFaceIds() = surf.faceIds();
     this->storedZones()  = surf.surfZones();
 }
 
@@ -1419,7 +1485,9 @@ Foam::MeshedSurface<Face>::operator Foam::MeshedSurfaceProxy<Face>() const
     (
         this->points(),
         this->surfFaces(),
-        this->surfZones()
+        this->surfZones(),
+        labelUList::null(), // faceMap = none
+        this->faceIds()
     );
 }
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index 68e2d0a9c9a903d3aefaee2f62277b583ee5e01e..0fbeb652e3d201d98d11eb1870805ac829ea588b 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -115,6 +115,11 @@ private:
 
     // Private Data
 
+        //- Face ids.
+        //  If these exist, they are typically arise from reading a mesh
+        //  format from another CAE software (eg, NASTRAN, STARCD, ...)
+        labelList faceIds_;
+
         //- Zone information
         // (face ordering nFaces/startFace only used during reading/writing)
         surfZoneList zones_;
@@ -164,6 +169,12 @@ protected:
             return static_cast<List<Face>&>(*this);
         }
 
+        //- Non-const access to face ids
+        labelList& storedFaceIds()
+        {
+            return faceIds_;
+        }
+
         //- Non-const access to the zones
         surfZoneList& storedZones()
         {
@@ -175,6 +186,7 @@ protected:
         (
             DynamicList<Face>& unsortedFaces,
             DynamicList<label>& zoneIds,
+            DynamicList<label>& elemIds,
             bool sorted
         );
 
@@ -394,6 +406,14 @@ public:
             return static_cast<const List<Face>&>(*this);
         }
 
+        //- Return const access to faces ids
+        //  If these exist, they are typically arise from reading a mesh
+        //  format from another CAE software (eg, NASTRAN, STARCD, ...)
+        const labelList& faceIds() const
+        {
+            return faceIds_;
+        }
+
         //- Const access to the surface zones.
         //  If zones are defined, they must be contiguous and cover the
         //  entire surface
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C b/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
index 20961c4f9af32437e1d41d801c5ca83737548703..aed06ff9355d41a3704209e7da0d0a685b1bc0e0 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,9 +52,11 @@ namespace Foam
         MeshedSurface<face>& surf
     )
     {
-        // First triangulate
-        // - slightly wasteful for space, but adjusts the zones too!
+        // First triangulate.
+        // Potentially wasteful of space, but adjusts zones and
+        // invalidates the faceIds
         surf.triangulate();
+
         this->storedPoints().transfer(surf.storedPoints());
         this->storedZones().transfer(surf.storedZones());
 
@@ -83,9 +85,11 @@ namespace Foam
         MeshedSurface<face>& surf
     )
     {
-        // First triangulate
-        // - slightly wasteful for space, but adjusts the zones too!
+        // First triangulate.
+        // Potentially wasteful of space, but adjusts zones and
+        // invalidates the faceIds
         surf.triangulate();
+
         this->storedPoints().transfer(surf.storedPoints());
         this->storedZones().transfer(surf.storedZones());
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
index 223c28afcd6975549013b4c03bd357f8efb70048..051a413cbec0ec093f8f696d31012d4bdf115f1b 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
@@ -71,6 +71,7 @@ void Foam::MeshedSurface<Face>::sortFacesAndStore
 (
     DynamicList<Face>& unsortedFaces,
     DynamicList<label>& zoneIds,
+    DynamicList<label>& elemIds,
     bool sorted
 )
 {
@@ -84,10 +85,16 @@ void Foam::MeshedSurface<Face>::sortFacesAndStore
         sorted = true;
     }
 
+    if (elemIds.size() != nInputFaces)
+    {
+        elemIds.clear();
+    }
+
     if (sorted)
     {
         // No additional sorting required
         this->storedFaces().transfer(unsortedFaces);
+        this->storedFaceIds().transfer(elemIds);
         return;
     }
 
@@ -106,6 +113,15 @@ void Foam::MeshedSurface<Face>::sortFacesAndStore
         // Can use transfer, faceMap is unique
         newFaces[facei].transfer(unsortedFaces[faceMap[facei]]);
     }
+
+    auto& newFaceIds = this->storedFaceIds();
+    newFaceIds.resize(elemIds.size());
+
+    // Element ids in sorted order
+    forAll(newFaceIds, facei)
+    {
+        newFaceIds[facei] = elemIds[faceMap[facei]];
+    }
 }
 
 
diff --git a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
index 9ec84f1520d335e1b1e111c5e7e407810a83986e..4f560f67d2d867e2fed2028c0f8235f3ab565658 100644
--- a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
+++ b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
@@ -238,13 +238,15 @@ Foam::MeshedSurfaceProxy<Face>::MeshedSurfaceProxy
     const pointField& pointLst,
     const UList<Face>& faceLst,
     const UList<surfZone>& zoneLst,
-    const labelUList& faceMap
+    const labelUList& faceMap,
+    const labelUList& faceIdsLst
 )
 :
     points_(pointLst),
     faces_(faceLst),
     zones_(zoneLst),
-    faceMap_(faceMap)
+    faceMap_(faceMap),
+    faceIds_(faceIdsLst)
 {}
 
 
diff --git a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
index 660d94756186fbdf89444a4db1158248e9b559a2..1da08c85da58992097ed8b21f9d66a1b6cb17a0d 100644
--- a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
+++ b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
@@ -31,6 +31,9 @@ Description
     A proxy for writing MeshedSurface, UnsortedMeshedSurface and surfMesh
     to various file formats.
 
+    The constructor interface is fat and ugly, but is largely encapsulated
+    by conversion operators in other classes.
+
 SourceFiles
     MeshedSurfaceProxy.C
     MeshedSurfaceProxys.C
@@ -73,6 +76,9 @@ class MeshedSurfaceProxy
 
         const UList<label>& faceMap_;
 
+        const UList<label>& faceIds_;
+
+
 public:
 
     // Public Typedefs
@@ -104,8 +110,9 @@ public:
         (
             const pointField& pointLst,
             const UList<Face>& faceLst,
-            const UList<surfZone>& zoneLst = List<surfZone>(),
-            const labelUList& faceMap = labelUList::null()
+            const UList<surfZone>& zoneLst = UList<surfZone>::null(),
+            const labelUList& faceMap = labelUList::null(),
+            const labelUList& faceIdLst = labelUList::null()
         );
 
 
@@ -187,12 +194,24 @@ public:
                 return faceMap_;
             }
 
+            //- Const access to the faceIds, zero-sized when unused
+            const labelUList& faceIds() const
+            {
+                return faceIds_;
+            }
+
             //- Can/should use faceMap?
             bool useFaceMap() const
             {
                 return faceMap_.size() == faces_.size();
             }
 
+            //- Possible to use faceIds?
+            bool useFaceIds() const
+            {
+                return faceIds_.size() == faces_.size();
+            }
+
             //- Count number of triangles.
             inline label nTriangles() const;
 
diff --git a/src/surfMesh/mergedSurf/mergedSurf.C b/src/surfMesh/mergedSurf/mergedSurf.C
index fc0f2ed500c0c1f17257f367dbadaba64ee2098f..2db9ada9b222476ae4248028eab3ac2d53427b98 100644
--- a/src/surfMesh/mergedSurf/mergedSurf.C
+++ b/src/surfMesh/mergedSurf/mergedSurf.C
@@ -61,6 +61,7 @@ Foam::mergedSurf::mergedSurf
     const pointField& unmergedPoints,
     const faceList& unmergedFaces,
     const labelList& origZoneIds,
+    const labelList& origFaceIds,
     const scalar mergeDim
 )
 :
@@ -71,6 +72,7 @@ Foam::mergedSurf::mergedSurf
         unmergedPoints,
         unmergedFaces,
         origZoneIds,
+        origFaceIds,
         mergeDim
     );
 }
@@ -91,6 +93,7 @@ void Foam::mergedSurf::clear()
     pointsMap_.clear();
 
     zoneIds_.clear();
+    faceIds_.clear();
 }
 
 
@@ -106,6 +109,7 @@ bool Foam::mergedSurf::merge
             unmergedSurface.points(),
             unmergedSurface.faces(),
             unmergedSurface.zoneIds(),
+            unmergedSurface.faceIds(),
             mergeDim
         );
 }
@@ -124,6 +128,7 @@ bool Foam::mergedSurf::merge
             unmergedPoints,
             unmergedFaces,
             labelList(),
+            labelList(),
             mergeDim
         );
 }
@@ -134,6 +139,7 @@ bool Foam::mergedSurf::merge
     const pointField& unmergedPoints,
     const faceList& unmergedFaces,
     const labelList& origZoneIds,
+    const labelList& origFaceIds,
     const scalar mergeDim
 )
 {
@@ -160,6 +166,7 @@ bool Foam::mergedSurf::merge
     // Now handle per-face information
 
     globalIndex::gatherOp(origZoneIds, zoneIds_);
+    globalIndex::gatherOp(origFaceIds, faceIds_);
 
     return true;
 }
diff --git a/src/surfMesh/mergedSurf/mergedSurf.H b/src/surfMesh/mergedSurf/mergedSurf.H
index f2adba45295930a3a06daa3628328db41b51a278..1b03e8e874308f1423ecc23168a19331d925f04e 100644
--- a/src/surfMesh/mergedSurf/mergedSurf.H
+++ b/src/surfMesh/mergedSurf/mergedSurf.H
@@ -96,6 +96,7 @@ public:
             const pointField& unmergedPoints,
             const faceList& unmergedFaces,
             const labelList& origZoneIds,
+            const labelList& origFaceIds,
             const scalar mergeDim
         );
 
@@ -133,6 +134,12 @@ public:
             return zoneIds_;
         }
 
+        //- Per-face identifier (eg, element Id)
+        virtual const labelList& faceIds() const
+        {
+            return faceIds_;
+        }
+
         //- Map for reordered points (old-to-new)
         const labelList& pointsMap() const
         {
@@ -166,6 +173,7 @@ public:
             const pointField& unmergedPoints,
             const faceList& unmergedFaces,
             const labelList& origZoneIds,
+            const labelList& origFaceIds,
             const scalar mergeDim
         );
 
diff --git a/src/surfMesh/meshedSurf/meshedSurf.H b/src/surfMesh/meshedSurf/meshedSurf.H
index 36420372e36008b31bcecaab9ea6bb4f771f5d26..6933b5eed21d4af3fc266493c71ce25d29a352fd 100644
--- a/src/surfMesh/meshedSurf/meshedSurf.H
+++ b/src/surfMesh/meshedSurf/meshedSurf.H
@@ -79,6 +79,12 @@ public:
         {
             return labelList::null();
         }
+
+        //- Per-face identifier (eg, element Id)
+        virtual const labelList& faceIds() const
+        {
+            return labelList::null();
+        }
 };
 
 
diff --git a/src/surfMesh/meshedSurf/meshedSurfRef.H b/src/surfMesh/meshedSurf/meshedSurfRef.H
index 7ec13c203b87dd223befd9fbea960127f776dae4..c52d35363b35592b7ef00330493cd71478a87846 100644
--- a/src/surfMesh/meshedSurf/meshedSurfRef.H
+++ b/src/surfMesh/meshedSurf/meshedSurfRef.H
@@ -53,6 +53,7 @@ class meshedSurfRef
     std::reference_wrapper<const pointField> points_;
     std::reference_wrapper<const faceList> faces_;
     std::reference_wrapper<const labelList> zoneIds_;
+    std::reference_wrapper<const labelList> faceIds_;
 
 
 public:
@@ -64,7 +65,8 @@ public:
         :
             points_(std::cref<pointField>(pointField::null())),
             faces_(std::cref<faceList>(faceList::null())),
-            zoneIds_(std::cref<labelList>(labelList::null()))
+            zoneIds_(std::cref<labelList>(labelList::null())),
+            faceIds_(std::cref<labelList>(labelList::null()))
         {}
 
 
@@ -73,12 +75,14 @@ public:
         (
             const pointField& pointLst,
             const faceList& faceLst,
-            const labelList& zoneIdLst = labelList::null()
+            const labelList& zoneIdLst = labelList::null(),
+            const labelList& faceIdLst = labelList::null()
         )
         :
             points_(std::cref<pointField>(pointLst)),
             faces_(std::cref<faceList>(faceLst)),
-            zoneIds_(std::cref<labelList>(zoneIdLst))
+            zoneIds_(std::cref<labelList>(zoneIdLst)),
+            faceIds_(std::cref<labelList>(faceIdLst))
         {}
 
 
@@ -106,12 +110,19 @@ public:
             return zoneIds_.get();
         }
 
+        //- Per-face identifier (eg, element Id)
+        virtual const labelList& faceIds() const
+        {
+            return faceIds_.get();
+        }
+
         //- Remove all references by redirecting to null objects
         void clear()
         {
             points_ = std::cref<pointField>(pointField::null());
             faces_ = std::cref<faceList>(faceList::null());
             zoneIds_ = std::cref<labelList>(labelList::null());
+            faceIds_ = std::cref<labelList>(labelList::null());
         }
 
         //- Reset components
@@ -119,12 +130,14 @@ public:
         (
             const pointField& pointLst,
             const faceList& faceLst,
-            const labelList& zoneIdLst = labelList::null()
+            const labelList& zoneIdLst = labelList::null(),
+            const labelList& faceIdLst = labelList::null()
         )
         {
             points_ = std::cref<pointField>(pointLst);
             faces_ = std::cref<faceList>(faceLst);
             zoneIds_ = std::cref<labelList>(zoneIdLst);
+            faceIds_ = std::cref<labelList>(faceIdLst);
         }
 };
 
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
index 34632ffa60c0289b93f556608065ac524f49fbe1..b76b3562e7a0d5a39f4bb78ee63119e76aadcab3 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
@@ -117,6 +117,8 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
 
     DynamicList<label>  pointId;    // Nastran point id (1-based)
     DynamicList<point>  dynPoints;
+
+    DynamicList<label>  dynElemId;  // Nastran element id (1-based)
     DynamicList<Face>   dynFaces;
     DynamicList<label>  dynZones;
     DynamicList<label>  dynSizes;
@@ -127,6 +129,9 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     label zoneId = 0;
     bool sorted = true;
 
+    // Element id gets trashed with decompose into a triangle!
+    bool ignoreElemId = false;
+
     // Name for face group
     Map<word> nameLookup;
 
@@ -222,7 +227,7 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
 
         if (cmd == "CTRIA3")
         {
-            (void) nextNasField(line, linei, 8); // 8-16
+            label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
             label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
             const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
             const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
@@ -247,13 +252,15 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
                 // Info<< "zone" << zoneId << " => group " << groupId <<nl;
             }
 
+            --elemId;   // Convert 1-based -> 0-based
+            dynElemId.append(elemId);
             dynFaces.append(Face{a, b, c});
             dynZones.append(zoneId);
             dynSizes[zoneId]++;
         }
         else if (cmd == "CQUAD4")
         {
-            (void) nextNasField(line, linei, 8); // 8-16
+            label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
             label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
             const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
             const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
@@ -281,6 +288,9 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
 
             if (faceTraits<Face>::isTri())
             {
+                ignoreElemId = true;
+                dynElemId.clear();
+
                 dynFaces.append(Face{a, b, c});
                 dynFaces.append(Face{c, d, a});
                 dynZones.append(zoneId);
@@ -289,6 +299,9 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
             }
             else
             {
+                --elemId;   // Convert 1-based -> 0-based
+
+                dynElemId.append(elemId);
                 dynFaces.append(Face{a,b,c,d});
                 dynZones.append(zoneId);
                 dynSizes[zoneId]++;
@@ -358,6 +371,10 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     //        << " points:" << dynPoints.size()
     //        << endl;
 
+    if (ignoreElemId)
+    {
+        dynElemId.clear();
+    }
 
     // Transfer to normal lists
     this->storedPoints().transfer(dynPoints);
@@ -403,7 +420,7 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
         }
     }
 
-    this->sortFacesAndStore(dynFaces, dynZones, sorted);
+    this->sortFacesAndStore(dynFaces, dynZones, dynElemId, sorted);
 
     // Add zones (retaining empty ones)
     this->addZones(dynSizes, names);
@@ -428,6 +445,7 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
     const UList<point>& pointLst = surf.points();
     const UList<Face>&  faceLst  = surf.surfFaces();
     const UList<label>& faceMap  = surf.faceMap();
+    const UList<label>& elemIds  = surf.faceIds();
 
     // for no zones, suppress the group name
     const surfZoneList zones =
@@ -439,6 +457,25 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
 
     const bool useFaceMap = (surf.useFaceMap() && zones.size() > 1);
 
+    // Possible to use faceIds?
+    bool useOrigFaceIds =
+        (!useFaceMap && elemIds.size() == faceLst.size());
+
+    if (useOrigFaceIds)
+    {
+        // Not possible with on-the-fly face decomposition
+
+        for (const auto& f : faceLst)
+        {
+            if (f.size() > 4)
+            {
+                useOrigFaceIds = false;
+                break;
+            }
+        }
+    }
+
+
     OFstream os(filename, streamOpt);
     if (!os.good())
     {
@@ -489,6 +526,11 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
 
             const Face& f = faceLst[facei];
 
+            if (useOrigFaceIds)
+            {
+                elemId = elemIds[facei];
+            }
+
             elemId = writeShell(os, f, elemId, zoneIndex);
         }
 
diff --git a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
index 25928e9732ccdc3dcdade5dc9f0bab4fdf01ae6e..a2b43c28eadd13d6aa64073f6d37012ea2f19901 100644
--- a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
@@ -71,6 +71,8 @@ bool Foam::fileFormats::OBJsurfaceFormat<Face>::read
 
     DynamicList<point> dynPoints;
     DynamicList<label> dynVerts;
+
+    DynamicList<label> dynElemId; // unused
     DynamicList<Face>  dynFaces;
 
     DynamicList<word>  dynNames;
@@ -207,7 +209,7 @@ bool Foam::fileFormats::OBJsurfaceFormat<Face>::read
     // Transfer to normal lists
     this->storedPoints().transfer(dynPoints);
 
-    this->sortFacesAndStore(dynFaces, dynZones, sorted);
+    this->sortFacesAndStore(dynFaces, dynZones, dynElemId, sorted);
 
     // Add zones (retaining empty ones)
     this->addZones(dynSizes, dynNames);
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
index 797cee86520d898483e3bdf0ebf0963b2052bfa5..d1db77d26f7601b64d697dba228e6bf372241c8b 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
@@ -127,7 +127,9 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
 
     readHeader(is, STARCDCore::HEADER_CEL);
 
+    DynamicList<label> dynElemId;  // STARCD element id (1-based)
     DynamicList<Face>  dynFaces;
+
     DynamicList<label> dynZones;
     DynamicList<word>  dynNames;
     DynamicList<label> dynSizes;
@@ -137,6 +139,9 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
     bool sorted = true;
     label zoneId = 0;
 
+    // Element id gets trashed with decompose into a triangle!
+    bool ignoreElemId = false;
+
     label ignoredLabel, shapeId, nLabels, cellTableId, typeId;
     DynamicList<label> vertexLabels(64);
 
@@ -144,7 +149,9 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
 
     while (is.read(tok).good() && tok.isLabel())
     {
-        // const label starCellId = tok.labelToken();
+        // First token is the element id (1-based)
+        label elemId = tok.labelToken();
+
         is  >> shapeId
             >> nLabels
             >> cellTableId
@@ -203,6 +210,8 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
             if (faceTraits<Face>::isTri() && nLabels > 3)
             {
                 // The face needs triangulation
+                ignoreElemId = true;
+                dynElemId.clear();
 
                 face f(vertices);
 
@@ -220,6 +229,9 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
             }
             else if (nLabels >= 3)
             {
+                --elemId;   // Convert 1-based -> 0-based
+                dynElemId.append(elemId);
+
                 dynFaces.append(Face(vertices));
                 dynZones.append(zoneId);
                 dynSizes[zoneId]++;
@@ -228,7 +240,14 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
     }
     mapPointId.clear();
 
-    this->sortFacesAndStore(dynFaces, dynZones, sorted);
+
+    if (ignoreElemId)
+    {
+        dynElemId.clear();
+    }
+
+
+    this->sortFacesAndStore(dynFaces, dynZones, dynElemId, sorted);
 
     // Add zones (retaining empty ones)
     this->addZones(dynSizes, dynNames);
@@ -253,6 +272,7 @@ void Foam::fileFormats::STARCDsurfaceFormat<Face>::write
     const UList<point>& pointLst = surf.points();
     const UList<Face>&  faceLst  = surf.surfFaces();
     const UList<label>& faceMap  = surf.faceMap();
+    const UList<label>& elemIds  = surf.faceIds();
 
     const surfZoneList zones =
     (
@@ -263,6 +283,15 @@ void Foam::fileFormats::STARCDsurfaceFormat<Face>::write
 
     const bool useFaceMap = (surf.useFaceMap() && zones.size() > 1);
 
+    // Possible to use faceIds?
+    const bool useOrigFaceIds =
+    (
+        !useFaceMap
+     && surf.useFaceIds()
+     && elemIds.size() == faceLst.size()
+    );
+
+
     fileName baseName = filename.lessExt();
 
     // The .vrt file
@@ -287,6 +316,11 @@ void Foam::fileFormats::STARCDsurfaceFormat<Face>::write
 
             const Face& f = faceLst[facei];
 
+            if (useOrigFaceIds)
+            {
+                elemId = elemIds[facei];
+            }
+
             writeShell(os, f, elemId, zoneIndex);
             ++elemId;
         }
diff --git a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
index 068c6c77f231ac657a1f79971926b3a4a94b2185..38f19dc70e81b64ca0168a6df062a1ab99ff101c 100644
--- a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
@@ -175,6 +175,8 @@ bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
         }
     }
 
+    DynamicList<label> dynElemId; // unused
+
     if (nTri > faces.size())
     {
         // We are here if the target surface needs triangles and
@@ -203,7 +205,7 @@ bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
             zoneSizes[zonei]++;
         }
 
-        this->sortFacesAndStore(dynFaces, dynZones, sorted);
+        this->sortFacesAndStore(dynFaces, dynZones, dynElemId, sorted);
 
         // Add zones (retaining empty ones)
         this->addZones(zoneSizes, zoneNames);
@@ -225,7 +227,7 @@ bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
             zoneSizes[zonei]++;
         }
 
-        this->sortFacesAndStore(dynFaces, dynZones, sorted);
+        this->sortFacesAndStore(dynFaces, dynZones, dynElemId, sorted);
 
         // Add zones (retaining empty ones)
         this->addZones(zoneSizes, zoneNames);
diff --git a/src/surfMesh/writers/nastran/nastranSurfaceWriter.C b/src/surfMesh/writers/nastran/nastranSurfaceWriter.C
index 9592e5bba920e4031ceb3293c463e1152f8e9efa..fcfae7d7607ce66c8ea266735d27ce2c03a667a8 100644
--- a/src/surfMesh/writers/nastran/nastranSurfaceWriter.C
+++ b/src/surfMesh/writers/nastran/nastranSurfaceWriter.C
@@ -211,6 +211,24 @@ void Foam::surfaceWriters::nastranWriter::writeGeometry
     const pointField& points = surf.points();
     const faceList&    faces = surf.faces();
     const labelList&   zones = surf.zoneIds();
+    const labelList& elemIds = surf.faceIds();
+
+    // Possible to use faceIds?
+    bool useOrigFaceIds = (elemIds.size() == faces.size());
+
+    if (useOrigFaceIds)
+    {
+        // Not possible with on-the-fly face decomposition
+        for (const auto& f : faces)
+        {
+            if (f.size() > 4)
+            {
+                useOrigFaceIds = false;
+                break;
+            }
+        }
+    }
+
 
     // Write points
 
@@ -238,6 +256,12 @@ void Foam::surfaceWriters::nastranWriter::writeGeometry
     {
         const face& f = faces[facei];
 
+        if (useOrigFaceIds)
+        {
+            // When available and not decomposed
+            elemId = elemIds[facei];
+        }
+
         // 1-offset for PID
         const label propId = 1 + (facei < zones.size() ? zones[facei] : 0);
 
diff --git a/src/surfMesh/writers/nastran/nastranSurfaceWriterImpl.C b/src/surfMesh/writers/nastran/nastranSurfaceWriterImpl.C
index 737240262c9a991ce2cfbe771dcb451ec2ca6000..56fee26250114c0cd0503aefd1376352f084d019 100644
--- a/src/surfMesh/writers/nastran/nastranSurfaceWriterImpl.C
+++ b/src/surfMesh/writers/nastran/nastranSurfaceWriterImpl.C
@@ -247,6 +247,16 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
 
         // Regular (undecomposed) faces
         const faceList& faces = surf.faces();
+        const labelList& elemIds = surf.faceIds();
+
+        // Possible to use faceIds?
+        // Not possible with on-the-fly face decomposition
+        const bool useOrigFaceIds =
+        (
+            elemIds.size() == faces.size()
+         && decompFaces.empty()
+        );
+
 
         label elemId = 0;
 
@@ -254,6 +264,12 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
         {
             forAll(faces, facei)
             {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
                 const label beginElemId = elemId;
 
                 // Any face decomposition
@@ -299,6 +315,12 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
 
             forAll(faces, facei)
             {
+                if (useOrigFaceIds)
+                {
+                    // When available and not decomposed
+                    elemId = elemIds[facei];
+                }
+
                 const Type v(varScale * *valIter);
                 ++valIter;
 
diff --git a/src/surfMesh/writers/starcd/starcdSurfaceWriter.C b/src/surfMesh/writers/starcd/starcdSurfaceWriter.C
index bad0418f448ec6843491299985b23fa8f6cbfff7..0aa30d2a81171078ceb3518f63f223fb466f1c28 100644
--- a/src/surfMesh/writers/starcd/starcdSurfaceWriter.C
+++ b/src/surfMesh/writers/starcd/starcdSurfaceWriter.C
@@ -145,7 +145,14 @@ Foam::fileName Foam::surfaceWriters::starcdWriter::write()
             mkDir(outputFile.path());
         }
 
-        MeshedSurfaceProxy<face>(surf.points(), surf.faces()).write
+        MeshedSurfaceProxy<face>
+        (
+            surf.points(),
+            surf.faces(),
+            UList<surfZone>::null(), // one zone
+            labelUList::null(),      // no faceMap
+            surf.faceIds()           // with face ids (if possible)
+        ).write
         (
             outputFile,
             "starcd",  // Canonical selection name