diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 3442889069920a5cbc37654598d20e839a85450d..079294190d6d87f8844848f099a7f87193d3df0a 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -9,6 +9,7 @@ sampledSet/midPoint/midPointSet.C
 sampledSet/midPointAndFace/midPointAndFaceSet.C
 sampledSet/sampledSet/sampledSet.C
 sampledSet/sampledSets/sampledSets.C
+sampledSet/sampledSets/sampledSetsGrouping.C
 sampledSet/sampledSetsFunctionObject/sampledSetsFunctionObject.C
 sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
 sampledSet/uniform/uniformSet.C
@@ -34,6 +35,7 @@ sampledSurface/distanceSurface/distanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C
+sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
 sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C
 sampledSurface/thresholdCellFaces/thresholdCellFaces.C
 sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index a964cab15293f9e5ccbf107391496604eccee903..600f89bb9b0389c65514532e26598c2b30941c99 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -34,103 +34,12 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-namespace Foam
-{
-    defineTypeNameAndDebug(sampledSets, 0);
-}
-
+defineTypeNameAndDebug(Foam::sampledSets, 0);
 bool Foam::sampledSets::verbose_ = false;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::sampledSets::checkFieldTypes()
-{
-    wordList fieldTypes(fieldNames_.size());
-
-    // check files for a particular time
-    if (loadFromFiles_)
-    {
-        forAll(fieldNames_, fieldi)
-        {
-            IOobject io
-            (
-                fieldNames_[fieldi],
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            );
-
-            if (io.headerOk())
-            {
-                fieldTypes[fieldi] = io.headerClassName();
-            }
-            else
-            {
-                fieldTypes[fieldi] = "(notFound)";
-            }
-        }
-    }
-    else
-    {
-        // check objectRegistry
-        forAll(fieldNames_, fieldi)
-        {
-            objectRegistry::const_iterator iter =
-                mesh_.find(fieldNames_[fieldi]);
-
-            if (iter != mesh_.objectRegistry::end())
-            {
-                fieldTypes[fieldi] = iter()->type();
-            }
-            else
-            {
-                fieldTypes[fieldi] = "(notFound)";
-            }
-        }
-    }
-
-
-    label nFields = 0;
-
-    // classify fieldTypes
-    nFields += grep(scalarFields_, fieldTypes);
-    nFields += grep(vectorFields_, fieldTypes);
-    nFields += grep(sphericalTensorFields_, fieldTypes);
-    nFields += grep(symmTensorFields_, fieldTypes);
-    nFields += grep(tensorFields_, fieldTypes);
-
-    if (Pstream::master())
-    {
-        if (debug)
-        {
-            Pout<< "timeName = " << mesh_.time().timeName() << nl
-                << "scalarFields    " << scalarFields_ << nl
-                << "vectorFields    " << vectorFields_ << nl
-                << "sphTensorFields " << sphericalTensorFields_ << nl
-                << "symTensorFields " << symmTensorFields_ <<nl
-                << "tensorFields    " << tensorFields_ <<nl;
-        }
-
-        if (nFields > 0)
-        {
-            if (debug)
-            {
-                Pout<< "Creating directory "
-                    << outputPath_/mesh_.time().timeName()
-                    << nl << endl;
-            }
-
-            mkDir(outputPath_/mesh_.time().timeName());
-        }
-    }
-
-    return nFields > 0;
-}
-
-
 void Foam::sampledSets::combineSampledSets
 (
     PtrList<coordSet>& masterSampledSets,
@@ -147,9 +56,9 @@ void Foam::sampledSets::combineSampledSets
 
     const PtrList<sampledSet>& sampledSets = *this;
 
-    forAll(sampledSets, seti)
+    forAll(sampledSets, setI)
     {
-        const sampledSet& samplePts = sampledSets[seti];
+        const sampledSet& samplePts = sampledSets[setI];
 
         // Collect data from all processors
         List<List<point> > gatheredPts(Pstream::nProcs());
@@ -190,7 +99,7 @@ void Foam::sampledSets::combineSampledSets
 
         // Sort curveDist and use to fill masterSamplePts
         SortableList<scalar> sortedDist(allCurveDist);
-        indexSets[seti] = sortedDist.indices();
+        indexSets[setI] = sortedDist.indices();
 
         // Get reference point (note: only master has all points)
         point refPt;
@@ -207,12 +116,12 @@ void Foam::sampledSets::combineSampledSets
 
         masterSampledSets.set
         (
-            seti,
+            setI,
             new coordSet
             (
                 samplePts.name(),
                 samplePts.axis(),
-                List<point>(UIndirectList<point>(allPts, indexSets[seti])),
+                List<point>(UIndirectList<point>(allPts, indexSets[setI])),
                 refPt
             )
         );
@@ -236,7 +145,7 @@ Foam::sampledSets::sampledSets
     loadFromFiles_(loadFromFiles),
     outputPath_(fileName::null),
     searchEngine_(mesh_, true),
-    fieldNames_(),
+    fieldSelection_(),
     interpolationScheme_(word::null),
     writeFormat_(word::null)
 {
@@ -285,13 +194,43 @@ void Foam::sampledSets::end()
 
 void Foam::sampledSets::write()
 {
-    if (size() && checkFieldTypes())
+    if (size())
     {
-        sampleAndWrite(scalarFields_);
-        sampleAndWrite(vectorFields_);
-        sampleAndWrite(sphericalTensorFields_);
-        sampleAndWrite(symmTensorFields_);
-        sampleAndWrite(tensorFields_);
+        const label nFields = classifyFields();
+
+        if (Pstream::master())
+        {
+            if (debug)
+            {
+                Pout<< "timeName = " << mesh_.time().timeName() << nl
+                    << "scalarFields    " << scalarFields_ << nl
+                    << "vectorFields    " << vectorFields_ << nl
+                    << "sphTensorFields " << sphericalTensorFields_ << nl
+                    << "symTensorFields " << symmTensorFields_ <<nl
+                    << "tensorFields    " << tensorFields_ <<nl;
+            }
+
+            if (nFields)
+            {
+                if (debug)
+                {
+                    Pout<< "Creating directory "
+                        << outputPath_/mesh_.time().timeName()
+                            << nl << endl;
+                }
+
+                mkDir(outputPath_/mesh_.time().timeName());
+            }
+        }
+
+        if (nFields)
+        {
+            sampleAndWrite(scalarFields_);
+            sampleAndWrite(vectorFields_);
+            sampleAndWrite(sphericalTensorFields_);
+            sampleAndWrite(symmTensorFields_);
+            sampleAndWrite(tensorFields_);
+        }
     }
 }
 
@@ -299,20 +238,20 @@ void Foam::sampledSets::write()
 void Foam::sampledSets::read(const dictionary& dict)
 {
     dict_ = dict;
+    dict_.lookup("fields") >> fieldSelection_;
+    clearFieldGroups();
 
-    fieldNames_ = wordList(dict_.lookup("fields"));
-
-    interpolationScheme_ = "cell";
-    dict_.readIfPresent("interpolationScheme", interpolationScheme_);
-
-    writeFormat_ = "null";
-    dict_.readIfPresent("setFormat", writeFormat_);
+    interpolationScheme_ = dict.lookupOrDefault<word>
+    (
+        "interpolationScheme",
+        "cell"
+    );
+    writeFormat_ = dict.lookupOrDefault<word>
+    (
+        "setFormat",
+        "null"
+    );
 
-    scalarFields_.clear();
-    vectorFields_.clear();
-    sphericalTensorFields_.clear();
-    symmTensorFields_.clear();
-    tensorFields_.clear();
 
     PtrList<sampledSet> newList
     (
@@ -324,7 +263,7 @@ void Foam::sampledSets::read(const dictionary& dict)
 
     if (Pstream::master() && debug)
     {
-        Pout<< "sample fields:" << fieldNames_ << nl
+        Pout<< "sample fields:" << fieldSelection_ << nl
             << "sample sets:" << nl << "(" << nl;
 
         forAll(*this, si)
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.H b/src/sampling/sampledSet/sampledSets/sampledSets.H
index da38f9e1c50ef5d43adf9426fb0517dbadd94a9a..7302d26deb9f6119a5072dbada454c7e2c1dd3d4 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.H
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "interpolation.H"
 #include "coordSet.H"
 #include "writer.H"
+#include "wordReList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,7 +55,7 @@ class dictionary;
 class fvMesh;
 
 /*---------------------------------------------------------------------------*\
-                      Class sampledSets Declaration
+                         Class sampledSets Declaration
 \*---------------------------------------------------------------------------*/
 
 class sampledSets
@@ -67,25 +68,40 @@ class sampledSets
         template<class Type>
         class fieldGroup
         :
-            public wordList
+            public DynamicList<word>
         {
         public:
 
-            //- Set formatter
-            autoPtr<writer<Type> > formatter;
+            //- The set formatter
+            autoPtr< writer<Type> > formatter;
 
             //- Construct null
             fieldGroup()
             :
-                wordList(0),
+                DynamicList<word>(0),
                 formatter(NULL)
             {}
 
+            //- Construct for a particular format
+            fieldGroup(const word& writeFormat)
+            :
+                DynamicList<word>(0),
+                formatter(writer<Type>::New(writeFormat))
+            {}
+
+            //- Reset format and field list
             void clear()
             {
-                wordList::clear();
+                DynamicList<word>::clear();
                 formatter.clear();
             }
+
+            //- Assign a new formatter
+            void operator=(const word& writeFormat)
+            {
+                formatter = writer<Type>::New(writeFormat);
+            }
+
         };
 
 
@@ -161,7 +177,7 @@ class sampledSets
         // Read from dictonary
 
             //- Names of fields to sample
-            wordList fieldNames_;
+            wordReList fieldSelection_;
 
             //- Interpolation scheme to use
             word interpolationScheme_;
@@ -187,16 +203,14 @@ class sampledSets
 
     // Private Member Functions
 
-        //- Classify field types, return true if nFields > 0
-        bool checkFieldTypes();
+        //- Clear old field groups
+        void clearFieldGroups();
 
-        //- Find the fields in the list of the given type, return count
-        template<class Type>
-        label grep
-        (
-            fieldGroup<Type>& fieldList,
-            const wordList& fieldTypes
-        ) const;
+        //- Append fieldName to the appropriate group
+        label appendFieldGroup(const word& fieldName, const word& fieldType);
+
+        //- Classify field types, returns the number of fields
+        label classifyFields();
 
         //- Combine points from all processors. Sort by curveDist and produce
         //  index list. Valid result only on master processor.
diff --git a/src/sampling/sampledSet/sampledSets/sampledSetsGrouping.C b/src/sampling/sampledSet/sampledSets/sampledSetsGrouping.C
new file mode 100644
index 0000000000000000000000000000000000000000..6e9fb68f0ec1a1c938a0fafe5f04838689a6430e
--- /dev/null
+++ b/src/sampling/sampledSet/sampledSets/sampledSetsGrouping.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledSets.H"
+#include "volFields.H"
+#include "IOobjectList.H"
+#include "stringListOps.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledSets::clearFieldGroups()
+{
+    scalarFields_.clear();
+    vectorFields_.clear();
+    sphericalTensorFields_.clear();
+    symmTensorFields_.clear();
+    tensorFields_.clear();
+}
+
+
+Foam::label Foam::sampledSets::appendFieldGroup
+(
+    const word& fieldName,
+    const word& fieldType
+)
+{
+    if (fieldType == volScalarField::typeName)
+    {
+        scalarFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volVectorField::typeName)
+    {
+        vectorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volSphericalTensorField::typeName)
+    {
+        sphericalTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volSymmTensorField::typeName)
+    {
+        symmTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volTensorField::typeName)
+    {
+        tensorFields_.append(fieldName);
+        return 1;
+    }
+
+    return 0;
+}
+
+
+Foam::label Foam::sampledSets::classifyFields()
+{
+    label nFields = 0;
+    clearFieldGroups();
+
+    if (loadFromFiles_)
+    {
+        // check files for a particular time
+        IOobjectList objects(mesh_, mesh_.time().timeName());
+        wordList allFields = objects.sortedNames();
+
+        labelList indices = findStrings(fieldSelection_, allFields);
+
+        forAll(indices, fieldI)
+        {
+            const word& fieldName = allFields[indices[fieldI]];
+
+            nFields += appendFieldGroup
+            (
+                fieldName,
+                objects.find(fieldName)()->headerClassName()
+            );
+        }
+    }
+    else
+    {
+        // check currently available fields
+        wordList allFields = mesh_.sortedNames();
+        labelList indices = findStrings(fieldSelection_, allFields);
+
+        forAll(indices, fieldI)
+        {
+            const word& fieldName = allFields[indices[fieldI]];
+
+            nFields += appendFieldGroup
+            (
+                fieldName,
+                mesh_.find(fieldName)()->type()
+            );
+        }
+    }
+
+    return nFields;
+}
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSet/sampledSets/sampledSetsTemplates.C b/src/sampling/sampledSet/sampledSets/sampledSetsTemplates.C
index c5cfecb8a8fb6c850e638676ca6f6fd255a81cb7..1450c8c7bf3fedfcff568906ccc8bebf4046c84f 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSetsTemplates.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSetsTemplates.C
@@ -46,10 +46,10 @@ Foam::sampledSets::volFieldSampler<Type>::volFieldSampler
         interpolation<Type>::New(interpolationScheme, field)
     );
 
-    forAll(samplers, seti)
+    forAll(samplers, setI)
     {
-        Field<Type>& values = this->operator[](seti);
-        const sampledSet& samples = samplers[seti];
+        Field<Type>& values = this->operator[](setI);
+        const sampledSet& samples = samplers[setI];
 
         values.setSize(samples.size());
         forAll(samples, samplei)
@@ -79,10 +79,10 @@ Foam::sampledSets::volFieldSampler<Type>::volFieldSampler
     List<Field<Type> >(samplers.size()),
     name_(field.name())
 {
-    forAll(samplers, seti)
+    forAll(samplers, setI)
     {
-        Field<Type>& values = this->operator[](seti);
-        const sampledSet& samples = samplers[seti];
+        Field<Type>& values = this->operator[](setI);
+        const sampledSet& samples = samplers[setI];
 
         values.setSize(samples.size());
         forAll(samples, samplei)
@@ -105,41 +105,12 @@ Foam::sampledSets::volFieldSampler<Type>::volFieldSampler
 {}
 
 
-template<class Type>
-Foam::label Foam::sampledSets::grep
-(
-    fieldGroup<Type>& fieldList,
-    const wordList& fieldTypes
-) const
-{
-    fieldList.setSize(fieldNames_.size());
-    label nFields = 0;
-
-    forAll(fieldNames_, fieldi)
-    {
-        if
-        (
-            fieldTypes[fieldi]
-         == GeometricField<Type, fvPatchField, volMesh>::typeName
-        )
-        {
-            fieldList[nFields] = fieldNames_[fieldi];
-            nFields++;
-        }
-    }
-
-    fieldList.setSize(nFields);
-
-    return nFields;
-}
-
-
 template<class Type>
 void Foam::sampledSets::writeSampleFile
 (
     const coordSet& masterSampleSet,
     const PtrList<volFieldSampler<Type> >& masterFields,
-    const label seti,
+    const label setI,
     const fileName& timeDir,
     const writer<Type>& formatter
 )
@@ -150,7 +121,7 @@ void Foam::sampledSets::writeSampleFile
     forAll(masterFields, fieldi)
     {
         valueSetNames[fieldi] = masterFields[fieldi].name();
-        valueSets[fieldi] = &masterFields[fieldi][seti];
+        valueSets[fieldi] = &masterFields[fieldi][setI];
     }
 
     fileName fName
@@ -180,11 +151,11 @@ void Foam::sampledSets::combineSampledValues
     {
         List<Field<T> > masterValues(indexSets.size());
 
-        forAll(indexSets, seti)
+        forAll(indexSets, setI)
         {
             // Collect data from all processors
             List<Field<T> > gatheredData(Pstream::nProcs());
-            gatheredData[Pstream::myProcNo()] = sampledFields[fieldi][seti];
+            gatheredData[Pstream::myProcNo()] = sampledFields[fieldi][setI];
             Pstream::gatherList(gatheredData);
 
             if (Pstream::master())
@@ -198,10 +169,10 @@ void Foam::sampledSets::combineSampledValues
                     )
                 );
 
-                masterValues[seti] = UIndirectList<T>
+                masterValues[setI] = UIndirectList<T>
                 (
                     allData,
-                    indexSets[seti]
+                    indexSets[setI]
                 )();
             }
         }
@@ -232,7 +203,7 @@ void Foam::sampledSets::sampleAndWrite
         // Create or use existing writer
         if (fields.formatter.empty())
         {
-            fields.formatter = writer<Type>::New(writeFormat_);
+            fields = writeFormat_;
         }
 
         // Storage for interpolated values
@@ -326,13 +297,13 @@ void Foam::sampledSets::sampleAndWrite
 
         if (Pstream::master())
         {
-            forAll(masterSampledSets_, seti)
+            forAll(masterSampledSets_, setI)
             {
                 writeSampleFile
                 (
-                    masterSampledSets_[seti],
+                    masterSampledSets_[setI],
                     masterFields,
-                    seti,
+                    setI,
                     outputPath_/mesh_.time().timeName(),
                     fields.formatter()
                 );
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index 7969b010cdb1f66c4d0fdbbc20618db93a60322c..a5155a9f49f24eef79aafc8e2d988adf5aa47405 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -33,11 +33,12 @@ License
 #include "mergePoints.H"
 #include "volPointInterpolation.H"
 
-#include "IOobjectList.H"
-#include "stringListOps.H"
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+defineTypeNameAndDebug(Foam::sampledSurfaces, 0);
+bool Foam::sampledSurfaces::verbose_ = false;
+Foam::scalar Foam::sampledSurfaces::mergeTol_ = 1e-10;
+
 namespace Foam
 {
     //- Used to offset faces in Pstream::combineOffset
@@ -63,102 +64,11 @@ namespace Foam
         }
     };
 
-
-    defineTypeNameAndDebug(sampledSurfaces, 0);
 }
 
 
-bool Foam::sampledSurfaces::verbose_(false);
-Foam::scalar Foam::sampledSurfaces::mergeTol_(1e-10);
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::sampledSurfaces::appendFieldType
-(
-    const word& fieldName,
-    const word& fieldType
-)
-{
-    if (fieldType == volScalarField::typeName)
-    {
-        scalarFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volVectorField::typeName)
-    {
-        vectorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volSphericalTensorField::typeName)
-    {
-        sphericalTensorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volSymmTensorField::typeName)
-    {
-        symmTensorFields_.append(fieldName);
-        return 1;
-    }
-    else if (fieldType == volTensorField::typeName)
-    {
-        tensorFields_.append(fieldName);
-        return 1;
-    }
-
-    return 0;
-}
-
-
-Foam::label Foam::sampledSurfaces::classifyFieldTypes()
-{
-    label nFields = 0;
-
-    scalarFields_.clear();
-    vectorFields_.clear();
-    sphericalTensorFields_.clear();
-    symmTensorFields_.clear();
-    tensorFields_.clear();
-
-    // check files for a particular time
-    if (loadFromFiles_)
-    {
-        IOobjectList objects(mesh_, mesh_.time().timeName());
-        wordList allFields = objects.sortedNames();
-
-        labelList indices = findStrings(fieldSelection_, allFields);
-
-        forAll(indices, fieldI)
-        {
-            const word& fieldName = allFields[indices[fieldI]];
-
-            nFields += appendFieldType
-            (
-                fieldName,
-                objects.find(fieldName)()->headerClassName()
-            );
-        }
-    }
-    else
-    {
-        wordList allFields = mesh_.sortedNames();
-        labelList indices = findStrings(fieldSelection_, allFields);
-
-        forAll(indices, fieldI)
-        {
-            const word& fieldName = allFields[indices[fieldI]];
-
-            nFields += appendFieldType
-            (
-                fieldName,
-                mesh_.find(fieldName)()->type()
-            );
-        }
-    }
-
-    return nFields;
-}
-
-
 void Foam::sampledSurfaces::writeGeometry() const
 {
     // Write to time directory under outputPath_
@@ -269,7 +179,7 @@ void Foam::sampledSurfaces::write()
         // finalize surfaces, merge points etc.
         update();
 
-        const label nFields = classifyFieldTypes();
+        const label nFields = classifyFields();
 
         if (Pstream::master())
         {
@@ -290,8 +200,8 @@ void Foam::sampledSurfaces::write()
             mkDir(outputPath_/mesh_.time().timeName());
         }
 
-        // write geometry first if required, or when no fields would otherwise
-        // be written
+        // write geometry first if required,
+        // or when no fields would otherwise be written
         if (nFields == 0 || genericFormatter_->separateFiles())
         {
             writeGeometry();
@@ -309,15 +219,7 @@ void Foam::sampledSurfaces::write()
 void Foam::sampledSurfaces::read(const dictionary& dict)
 {
     dict.lookup("fields") >> fieldSelection_;
-
-    // might be okay for a size estimate, but we don't really know
-    const label nFields = fieldSelection_.size();
-
-    scalarFields_.reset(nFields);
-    vectorFields_.reset(nFields);
-    sphericalTensorFields_.reset(nFields);
-    symmTensorFields_.reset(nFields);
-    tensorFields_.reset(nFields);
+    clearFieldGroups();
 
     interpolationScheme_ = dict.lookupOrDefault<word>
     (
@@ -340,7 +242,6 @@ void Foam::sampledSurfaces::read(const dictionary& dict)
         dict.lookup("surfaces"),
         sampledSurface::iNew(mesh_)
     );
-
     transfer(newList);
 
     if (Pstream::parRun())
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index cd28e1f2ca18f63ef3c546ecfff1cdebb52f5285..66b6472ff923f2ceee9f255c7555133bd005e7c1 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -69,7 +69,7 @@ class sampledSurfaces
         {
         public:
 
-            //- Surface formatter
+            //- The surface formatter
             autoPtr< surfaceWriter<Type> > formatter;
 
             //- Construct null
@@ -86,34 +86,19 @@ class sampledSurfaces
                 formatter(surfaceWriter<Type>::New(writeFormat))
             {}
 
-            //- Construct for a particular surface format and a list of field
-            //  names
-            fieldGroup
-            (
-                const word& writeFormat,
-                const wordList& fieldNames
-            )
-            :
-                DynamicList<word>(fieldNames),
-                formatter(surfaceWriter<Type>::New(writeFormat))
-            {}
-
-            void reset(const label nElem)
+            //- Reset format and field list
+            void clear()
             {
-                formatter.clear();
-                DynamicList<word>::reserve(nElem);
                 DynamicList<word>::clear();
+                formatter.clear();
             }
 
+            //- Assign a new formatter
             void operator=(const word& writeFormat)
             {
                 formatter = surfaceWriter<Type>::New(writeFormat);
             }
 
-            void operator=(const wordList& fieldNames)
-            {
-                DynamicList<word>::operator=(fieldNames);
-            }
         };
 
 
@@ -192,11 +177,14 @@ class sampledSurfaces
 
     // Private Member Functions
 
+        //- Clear old field groups
+        void clearFieldGroups();
+
         //- Append fieldName to the appropriate group
-        label appendFieldType(const word& fieldName, const word& fieldType);
+        label appendFieldGroup(const word& fieldName, const word& fieldType);
 
         //- Classify field types, returns the number of fields
-        label classifyFieldTypes();
+        label classifyFields();
 
         //- Write geometry only
         void writeGeometry() const;
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
new file mode 100644
index 0000000000000000000000000000000000000000..497c36f8423b386bdbbc601a8fdd4f5973cec45e
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesGrouping.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledSurfaces.H"
+#include "volFields.H"
+#include "IOobjectList.H"
+#include "stringListOps.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledSurfaces::clearFieldGroups()
+{
+    scalarFields_.clear();
+    vectorFields_.clear();
+    sphericalTensorFields_.clear();
+    symmTensorFields_.clear();
+    tensorFields_.clear();
+}
+
+
+Foam::label Foam::sampledSurfaces::appendFieldGroup
+(
+    const word& fieldName,
+    const word& fieldType
+)
+{
+    if (fieldType == volScalarField::typeName)
+    {
+        scalarFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volVectorField::typeName)
+    {
+        vectorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volSphericalTensorField::typeName)
+    {
+        sphericalTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volSymmTensorField::typeName)
+    {
+        symmTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == volTensorField::typeName)
+    {
+        tensorFields_.append(fieldName);
+        return 1;
+    }
+
+    return 0;
+}
+
+
+Foam::label Foam::sampledSurfaces::classifyFields()
+{
+    label nFields = 0;
+    clearFieldGroups();
+
+    // check files for a particular time
+    if (loadFromFiles_)
+    {
+        IOobjectList objects(mesh_, mesh_.time().timeName());
+        wordList allFields = objects.sortedNames();
+
+        labelList indices = findStrings(fieldSelection_, allFields);
+
+        forAll(indices, fieldI)
+        {
+            const word& fieldName = allFields[indices[fieldI]];
+
+            nFields += appendFieldGroup
+            (
+                fieldName,
+                objects.find(fieldName)()->headerClassName()
+            );
+        }
+    }
+    else
+    {
+        wordList allFields = mesh_.sortedNames();
+        labelList indices = findStrings(fieldSelection_, allFields);
+
+        forAll(indices, fieldI)
+        {
+            const word& fieldName = allFields[indices[fieldI]];
+
+            nFields += appendFieldGroup
+            (
+                fieldName,
+                mesh_.find(fieldName)()->type()
+            );
+        }
+    }
+
+    return nFields;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
index 6da55cfe8736220134b921afdef9403ea9f3fb6e..b5dc2d418ce45f63e8ac49e5470eb0d8dca65dd2 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
@@ -178,7 +178,7 @@ void Foam::sampledSurfaces::sampleAndWrite
 
                 if
                 (
-                    iter != mesh_.objectRegistry::end()
+                    iter != objectRegistry::end()
                  && iter()->type()
                  == GeometricField<Type, fvPatchField, volMesh>::typeName
                 )