diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 7db0fb783e..3b225adc6a 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -78,6 +78,32 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::sampledTriSurfaceMesh::setZoneMap
+(
+    const surfZoneList& zoneLst,
+    labelList& zoneIds
+)
+{
+    label sz = 0;
+    forAll(zoneLst, zonei)
+    {
+        const surfZone& zn = zoneLst[zonei];
+        sz += zn.size();
+    }
+
+    zoneIds.setSize(sz);
+    forAll(zoneLst, zonei)
+    {
+        const surfZone& zn = zoneLst[zonei];
+
+        // Assign sub-zone Ids
+        SubList<label>(zoneIds, zn.size(), zn.start()) = zonei;
+    }
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 const Foam::indexedOctree<Foam::treeDataFace>&
@@ -146,16 +172,11 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
 
     // Global numbering for cells/faces - only used to uniquely identify local
     // elements
-    globalIndex globalCells
-    (
-        (sampleSource_ == cells || sampleSource_ == insideCells)
-      ? mesh().nCells()
-      : mesh().nFaces()
-    );
+    globalIndex globalCells(onBoundary() ? mesh().nFaces() : mesh().nCells());
 
     forAll(nearest, i)
     {
-        nearest[i].first() = GREAT;
+        nearest[i].first()  = GREAT;
         nearest[i].second() = labelMax;
     }
 
@@ -174,7 +195,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
             );
             if (nearInfo.hit())
             {
-                nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
+                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
                 nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
             }
         }
@@ -192,7 +213,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
                 const label index = cellTree.findInside(fc[triI]);
                 if (index != -1)
                 {
-                    nearest[triI].first() = 0.0;
+                    nearest[triI].first()  = 0.0;
                     nearest[triI].second() = globalCells.toGlobal(index);
                 }
             }
@@ -216,7 +237,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
             );
             if (nearInfo.hit())
             {
-                nearest[triI].first() = magSqr(nearInfo.hitPoint()-fc[triI]);
+                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
                 nearest[triI].second() = globalCells.toGlobal
                 (
                     bTree.shapes().faceLabels()[nearInfo.index()]
@@ -270,6 +291,37 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
     // From original point to compact points
     labelList reversePointMap(s.points().size(), -1);
 
+    // Handle region-wise sorting (makes things slightly more complicated)
+    zoneIds_.setSize(s.size(), -1);
+
+    // Better not to use triSurface::sortedZones here,
+    // since we'll sort ourselves
+
+    // Get zone/region sizes used, store under the original region Id
+    Map<label> zoneSizes;
+
+    // Recover region names from the input surface
+    Map<word> zoneNames;
+    {
+        const geometricSurfacePatchList& patches = s.patches();
+
+        forAll(patches, patchi)
+        {
+            zoneNames.set
+            (
+                patchi,
+                (
+                    patches[patchi].name().empty()
+                  ? Foam::name("patch%d", patchi)
+                  : patches[patchi].name()
+                )
+            );
+
+            zoneSizes.set(patchi, 0);
+        }
+    }
+
+
     {
         label newPointi = 0;
         label newFacei = 0;
@@ -278,57 +330,139 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
         {
             if (cellOrFaceLabels[facei] != -1)
             {
-                faceMap[newFacei++] = facei;
-
                 const triSurface::FaceType& f = s[facei];
+                const label regionid = f.region();
+
+                Map<label>::iterator fnd = zoneSizes.find(regionid);
+                if (fnd != zoneSizes.end())
+                {
+                    fnd()++;
+                }
+                else
+                {
+                    // This shouldn't happen
+                    zoneSizes.insert(regionid, 1);
+                    zoneNames.set
+                    (
+                        regionid,
+                        Foam::name("patch%d", regionid)
+                    );
+                }
+
+                // Store new faces compact
+                faceMap[newFacei] = facei;
+                zoneIds_[newFacei] = regionid;
+                ++newFacei;
+
+                // Renumber face points
                 forAll(f, fp)
                 {
-                    if (reversePointMap[f[fp]] == -1)
+                    const label labI = f[fp];
+
+                    if (reversePointMap[labI] == -1)
                     {
-                        pointMap[newPointi] = f[fp];
-                        reversePointMap[f[fp]] = newPointi++;
+                        pointMap[newPointi] = labI;
+                        reversePointMap[labI] = newPointi++;
                     }
                 }
             }
         }
+
+        // Trim
         faceMap.setSize(newFacei);
+        zoneIds_.setSize(newFacei);
         pointMap.setSize(newPointi);
     }
 
 
+    // Assign start/size (and name) to the newZones
+    // re-use the lookup to map (zoneId => zoneI)
+    surfZoneList zoneLst(zoneSizes.size());
+    label start = 0;
+    label zoneI = 0;
+    forAllIter(Map<label>, zoneSizes, iter)
+    {
+        // No negative regionids, so Map<label> sorts properly
+        const label regionid = iter.key();
+
+        word name;
+        Map<word>::const_iterator fnd = zoneNames.find(regionid);
+        if (fnd != zoneNames.end())
+        {
+            name = fnd();
+        }
+        if (name.empty())
+        {
+            name = ::Foam::name("patch%d", regionid);
+        }
+
+        zoneLst[zoneI] = surfZone
+        (
+            name,
+            0,           // initialize with zero size
+            start,
+            zoneI
+        );
+
+        // Adjust start for the next zone and save (zoneId => zoneI) mapping
+        start += iter();
+        iter() = zoneI++;
+    }
+
+
+    // At this stage:
+    // - faceMap to map the (unsorted) compact to original triangle
+    // - zoneIds for the next sorting
+    // - zoneSizes contains region -> count information
+
+    // Rebuild the faceMap for the sorted order
+    labelList sortedFaceMap(faceMap.size());
+
+    forAll(zoneIds_, facei)
+    {
+        label zonei = zoneIds_[facei];
+        label sortedFacei = zoneLst[zonei].start() + zoneLst[zonei].size()++;
+        sortedFaceMap[sortedFacei] = faceMap[facei];
+    }
+
+    // zoneIds are now simply flat values
+    setZoneMap(zoneLst, zoneIds_);
+
+    // Replace the faceMap with the properly sorted face map
+    faceMap.transfer(sortedFaceMap);
+
     if (keepIds_)
     {
         originalIds_ = faceMap;
     }
 
-    // Subset cellOrFaceLabels
+    // Subset cellOrFaceLabels (for compact faces)
     cellOrFaceLabels = UIndirectList<label>(cellOrFaceLabels, faceMap)();
 
     // Store any face per point (without using pointFaces())
     labelList pointToFace(pointMap.size());
 
     // Create faces and points for subsetted surface
-    faceList& faces = this->storedFaces();
-    faces.setSize(faceMap.size());
+    faceList& surfFaces = this->storedFaces();
+    surfFaces.setSize(faceMap.size());
 
-    labelList& zoneIds = this->storedZoneIds();
-    zoneIds.setSize(faceMap.size());
+    this->storedZones().transfer(zoneLst);
 
-    forAll(faceMap, i)
+    forAll(faceMap, facei)
     {
-        const labelledTri& f = s[faceMap[i]];
-        triFace newF
-        (
-            reversePointMap[f[0]],
-            reversePointMap[f[1]],
-            reversePointMap[f[2]]
-        );
-        faces[i]   = newF.triFaceFace();
-        zoneIds[i] = f.region(); // preserve zone information
+        const labelledTri& origF = s[faceMap[facei]];
+        face& f = surfFaces[facei];
 
-        forAll(newF, fp)
+        f = triFace
+        (
+            reversePointMap[origF[0]],
+            reversePointMap[origF[1]],
+            reversePointMap[origF[2]]
+        );
+
+        forAll(f, fp)
         {
-            pointToFace[newF[fp]] = i;
+            pointToFace[f[fp]] = facei;
         }
     }
 
@@ -418,7 +552,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
                 const point& pt = points()[pointi];
                 label facei = cellOrFaceLabels[pointToFace[pointi]];
                 sampleElements_[pointi] = facei;
-                samplePoints_[pointi] =  mesh().faces()[facei].nearestPoint
+                samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
                 (
                     pt,
                     mesh().points()
@@ -429,16 +563,16 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
     else
     {
         // if sampleSource_ == cells:
-        //      samplePoints_   : n/a
         //      sampleElements_ : per surface triangle the cell
+        //      samplePoints_   : n/a
         // if sampleSource_ == insideCells:
-        //      samplePoints_   : n/a
         //      sampleElements_ : -1 or per surface triangle the cell
-        // else:
         //      samplePoints_   : n/a
+        // else:
         //      sampleElements_ : per surface triangle the boundary face
-        samplePoints_.clear();
+        //      samplePoints_   : n/a
         sampleElements_.transfer(cellOrFaceLabels);
+        samplePoints_.clear();
     }
 
 
@@ -448,77 +582,47 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
         OFstream str(mesh().time().path()/"surfaceToMesh.obj");
         Info<< "Dumping correspondence from local surface (points or faces)"
             << " to mesh (cells or faces) to " << str.name() << endl;
-        label vertI = 0;
+
+        const vectorField& centres =
+        (
+            onBoundary()
+          ? mesh().faceCentres()
+          : mesh().cellCentres()
+        );
 
         if (sampledSurface::interpolate())
         {
-            if (sampleSource_ == cells || sampleSource_ == insideCells)
+            label vertI = 0;
+            forAll(samplePoints_, pointi)
             {
-                forAll(samplePoints_, pointi)
-                {
-                    meshTools::writeOBJ(str, points()[pointi]);
-                    vertI++;
+                meshTools::writeOBJ(str, points()[pointi]);
+                vertI++;
 
-                    meshTools::writeOBJ(str, samplePoints_[pointi]);
-                    vertI++;
+                meshTools::writeOBJ(str, samplePoints_[pointi]);
+                vertI++;
 
-                    label celli = sampleElements_[pointi];
-                    meshTools::writeOBJ(str, mesh().cellCentres()[celli]);
-                    vertI++;
-                    str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
-                        << nl;
-                }
-            }
-            else
-            {
-                forAll(samplePoints_, pointi)
-                {
-                    meshTools::writeOBJ(str, points()[pointi]);
-                    vertI++;
-
-                    meshTools::writeOBJ(str, samplePoints_[pointi]);
-                    vertI++;
-
-                    label facei = sampleElements_[pointi];
-                    meshTools::writeOBJ(str, mesh().faceCentres()[facei]);
-                    vertI++;
-                    str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
-                        << nl;
-                }
+                label elemi = sampleElements_[pointi];
+                meshTools::writeOBJ(str, centres[elemi]);
+                vertI++;
+                str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
             }
         }
         else
         {
-            if (sampleSource_ == cells || sampleSource_ == insideCells)
+            label vertI = 0;
+            forAll(sampleElements_, triI)
             {
-                forAll(sampleElements_, triI)
-                {
-                    meshTools::writeOBJ(str, faceCentres()[triI]);
-                    vertI++;
+                meshTools::writeOBJ(str, faceCentres()[triI]);
+                vertI++;
 
-                    label celli = sampleElements_[triI];
-                    meshTools::writeOBJ(str, mesh().cellCentres()[celli]);
-                    vertI++;
-                    str << "l " << vertI-1 << ' ' << vertI << nl;
-                }
-            }
-            else
-            {
-                forAll(sampleElements_, triI)
-                {
-                    meshTools::writeOBJ(str, faceCentres()[triI]);
-                    vertI++;
-
-                    label facei = sampleElements_[triI];
-                    meshTools::writeOBJ(str, mesh().faceCentres()[facei]);
-                    vertI++;
-                    str << "l " << vertI-1 << ' ' << vertI << nl;
-                }
+                label elemi = sampleElements_[triI];
+                meshTools::writeOBJ(str, centres[elemi]);
+                vertI++;
+                str << "l " << vertI-1 << ' ' << vertI << nl;
             }
         }
     }
 
-
     needsUpdate_ = false;
     return true;
 }
@@ -553,6 +657,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     needsUpdate_(true),
     keepIds_(false),
     originalIds_(),
+    zoneIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -584,6 +689,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     needsUpdate_(true),
     keepIds_(dict.lookupOrDefault<Switch>("keepIds", false)),
     originalIds_(),
+    zoneIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -617,6 +723,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     needsUpdate_(true),
     keepIds_(false),
     originalIds_(),
+    zoneIds_(),
     sampleElements_(0),
     samplePoints_(0)
 {}
@@ -646,6 +753,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
 
     sampledSurface::clearGeom();
     MeshStorage::clear();
+    zoneIds_.clear();
 
     originalIds_.clear();
     boundaryTreePtr_.clear();
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
index 6cd96a4770..b46d7649f6 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
@@ -62,6 +62,7 @@ Description
 
 SourceFiles
     sampledTriSurfaceMesh.C
+    sampledTriSurfaceMeshTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -71,6 +72,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "triSurfaceMesh.H"
 #include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,7 +89,7 @@ class meshSearch;
 class sampledTriSurfaceMesh
 :
     public sampledSurface,
-    public UnsortedMeshedSurface<face>
+    public meshedSurface
 {
 public:
         //- Types of communications
@@ -95,13 +97,13 @@ public:
         {
             cells,
             insideCells,
-            boundaryFaces,
+            boundaryFaces
         };
 
 private:
 
     //- Private typedefs for convenience
-        typedef UnsortedMeshedSurface<face> MeshStorage;
+        typedef meshedSurface MeshStorage;
 
 
     // Private data
@@ -127,6 +129,9 @@ private:
         //- Search tree for all non-coupled boundary faces
         mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
 
+        //- For compatibility with the meshSurf interface
+        labelList zoneIds_;
+
         //- From local surface triangle to mesh cell/face.
         labelList sampleElements_;
 
@@ -194,6 +199,10 @@ public:
 
     // Member Functions
 
+        //- Set new zoneIds list based on the surfZoneList information
+        static void setZoneMap(const surfZoneList&, labelList& zoneIds);
+
+
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -226,7 +235,7 @@ public:
         //- Const access to per-face zone/region information
         virtual const labelList& zoneIds() const
         {
-            return MeshStorage::zoneIds();
+            return zoneIds_;
         }
 
         //- Face area vectors
@@ -260,6 +269,12 @@ public:
             return originalIds_;
         }
 
+        //- Sampling boundary values instead of cell values
+        bool onBoundary() const
+        {
+            return sampleSource_ == boundaryFaces;
+        }
+
 
         //- Sample field on surface
         virtual tmp<scalarField> sample
@@ -323,7 +338,7 @@ public:
             const interpolation<tensor>&
         ) const;
 
-        //- Write
+        //- Write information
         virtual void print(Ostream&) const;
 };
 
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 880b7fd5b7..9cfba3874c 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,7 +38,7 @@ Foam::sampledTriSurfaceMesh::sampleField
     tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
     Field<Type>& values = tvalues.ref();
 
-    if (sampleSource_ == cells || sampleSource_ == insideCells)
+    if (!onBoundary())
     {
         // Sample cells
 
@@ -52,7 +52,7 @@ Foam::sampledTriSurfaceMesh::sampleField
         // Sample boundary faces
 
         const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-        label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+        const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
 
         // Create flat boundary field
 
@@ -94,7 +94,7 @@ Foam::sampledTriSurfaceMesh::interpolateField
     tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
     Field<Type>& values = tvalues.ref();
 
-    if (sampleSource_ == cells || sampleSource_ == insideCells)
+    if (!onBoundary())
     {
         // Sample cells.