diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4c85cd819fbc1c8572e63a21ee340aa98d03511e..965197d4727d8d21420e06f809a4a9546cb0aa1e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -202,6 +202,7 @@ ENDIF ()
 INCLUDE (ThirdPartyTinyxml)
 INCLUDE (ThirdPartyMetis)
 INCLUDE (ThirdPartyHDF5)
+INCLUDE (ThirdPartySIONLIB)
 INCLUDE (ThirdPartyScotch)
 INCLUDE (ThirdPartyZlib)
 INCLUDE (ThirdPartyBoost)
diff --git a/cmake/FindSIONLIB.cmake b/cmake/FindSIONLIB.cmake
new file mode 100755
index 0000000000000000000000000000000000000000..3653df188e6d30eb09495b57941f64b0d20a8249
--- /dev/null
+++ b/cmake/FindSIONLIB.cmake
@@ -0,0 +1,150 @@
+# Module that checks whether SIONLIB is available.
+#
+# Variables used by this module which you may want to set:
+# SIONLIB_ROOTDIR   Path list to search for SIONLIB
+# SIONLIB_SUFFIX    suffix to the library name , e.g. gcc or something
+# SIONLIB_INCDIR    directory with SIONLIB headers inside
+# SIONLIB_LIBDIR    directory with SIONLIB libraries inside
+#
+# Sets the following variables
+#
+# SIONLIB_INCLUDE_DIR    Path to the SIONLIB include dir
+# SIONLIB_MPI_LIBRARY
+# SIONLIB_GEN_LIBRARY
+# SIONLIB_SER_LIBRARY
+# SIONLIB_COM_LIBRARY
+# SIONLIB_COMLOCKNONE_LIBRARY
+# SIONLIB_LIBRARY_DIR
+# SIONLIB_FOUND          True if SIONLIB was found and usable
+# HAVE_SIONLIB           True if SIONLIB was found and usable
+# SIONLIB_LIBRARIES      Names of the SIONLIB libraries
+#
+
+set(SIONLIB_ROOTDIR $ENV{SIONLIB_DIR} CACHE PATH "Path list to search for SIONLIB" FORCE)
+set(SIONLIB_SUFFIX "_lib64" CACHE STRING "suffix to the library name , e.g. gcc or something")
+set(SIONLIB_INCDIR $ENV{SIONLIB_INC} CACHE PATH "directory with SIONLIB headers inside" FORCE)
+set(SIONLIB_LIBDIR $ENV{SIONLIB_LIB} CACHE PATH "directory with SIONLIB libraries inside" FORCE)
+
+MESSAGE(STATUS "SIONLIB_ROOTDIR=${SIONLIB_ROOTDIR}")
+MESSAGE(STATUS "SIONLIB_INCDIR=${SIONLIB_INCDIR}")
+MESSAGE(STATUS "SIONLIB_LIBDIR=${SIONLIB_LIBDIR}")
+
+mark_as_advanced(SIONLIB_ROOTDIR SIONLIB_SUFFIX SIONLIB_INCDIR SIONLIB_LIBDIR)
+
+
+#look for header files at positions given by the user
+find_path(SIONLIB_INCLUDE_DIR
+  NAMES "sion.h"
+  PATHS ${SIONLIB_INCDIR}
+  PATH_SUFFIXES "include" 
+  NO_DEFAULT_PATH
+  NO_SYSTEM_ENVIRONMENT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_INCLUDE_DIR=${SIONLIB_INCLUDE_DIR}")
+
+# check header usability
+include(CMakePushCheckState)
+cmake_push_check_state()
+set(CMAKE_REQUIRED_DEFINITIONS "${CMAKE_REQUIRED_DEFINITIONS} ${MPI_COMPILE_FLAGS}")
+set(CMAKE_REQUIRED_INCLUDES ${CMAKE_REQUIRED_INCLUDES} ${MPI_INCLUDE_PATH} ${SIONLIB_INCLUDE_DIR})
+set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} ${MPI_LIBRARIES})
+check_include_files(sion.h SIONLIB_HEADER_USABLE)
+
+
+find_library(SIONLIB_MPI_LIBRARY
+  NAMES "sionmpi_64"
+  PATHS ${SIONLIB_LIBDIR}
+  PATH_SUFFIXES "lib"
+  NO_DEFAULT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_MPI_LIBRARY=${SIONLIB_MPI_LIBRARY}")
+
+find_library(SIONLIB_GEN_LIBRARY
+  NAMES "siongen_64"
+  PATHS ${SIONLIB_LIBDIR}
+  PATH_SUFFIXES "lib"
+  NO_DEFAULT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_GEN_LIBRARY=${SIONLIB_GEN_LIBRARY}")
+
+find_library(SIONLIB_SER_LIBRARY
+  NAMES "sionser_64"
+  PATHS ${SIONLIB_LIBDIR}
+  PATH_SUFFIXES "lib"
+  NO_DEFAULT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_SER_LIBRARY=${SIONLIB_SER_LIBRARY}")
+
+find_library(SIONLIB_COM_LIBRARY
+  NAMES "sioncom_64"
+  PATHS ${SIONLIB_LIBDIR}
+  PATH_SUFFIXES "lib"
+  NO_DEFAULT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_COM_LIBRARY=${SIONLIB_COM_LIBRARY}")
+
+find_library(SIONLIB_COMLOCKNONE_LIBRARY
+  NAMES "sioncom_64_lock_none"
+  PATHS ${SIONLIB_LIBDIR}
+  PATH_SUFFIXES "lib"
+  NO_DEFAULT_PATH
+)
+
+MESSAGE(STATUS "SIONLIB_COMLOCKNONE_LIBRARY=${SIONLIB_COMLOCKNONE_LIBRARY}")
+
+
+cmake_pop_check_state()
+
+# behave like a CMake module is supposed to behave
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(
+  "SIONlib"
+  DEFAULT_MSG
+  SIONLIB_INCLUDE_DIR
+  SIONLIB_MPI_LIBRARY
+  SIONLIB_GEN_LIBRARY
+  SIONLIB_SER_LIBRARY
+  SIONLIB_COM_LIBRARY
+  SIONLIB_COMLOCKNONE_LIBRARY
+)
+
+mark_as_advanced(SIONLIB_INCLUDE_DIR SIONLIB_MPI_LIBRARY SIONLIB_GEN_LIBRARY SIONLIB_SER_LIBRARY SIONLIB_COM_LIBRARY SIONLIB_COMLOCKNONE_LIBRARY)
+
+# if both headers and library are found, store results
+if(SIONLIB_FOUND)
+  set(SIONLIB_LIBRARY_DIR ${SIONLIB_LIBDIR})
+  set(SIONLIB_LIBRARIES ${SIONLIB_MPI_LIBRARY} ${SIONLIB_GEN_LIBRARY} ${SIONLIB_SER_LIBRARY} ${SIONLIB_COM_LIBRARY} ${SIONLIB_COMLOCKNONE_LIBRARY})
+  # log result
+  file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeOutput.log
+    "Determining location of SIONLIB succeded:\n"
+    "Include directory: ${SIONLIB_INCLUDE_DIR}\n"
+    "Library directory: ${SIONLIB_LIBRARY_DIR}\n\n")
+  set(SIONLIB_COMPILE_FLAGS "-I${SIONLIB_INCLUDE_DIR} -D_SION_XT -DSION_MPI"
+    CACHE STRING "Compile Flags used by Nektar++ when compiling with SIONLIB libraries")
+  set(SIONLIB_LIBRARIES ${SIONLIB_LIBRARIES} 
+    CACHE STRING "Libraries used by Nektar++ when linking with SIONLIB libraries")
+else(SIONLIB_FOUND)
+  # log errornous result
+  file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log
+    "Determing location of SIONLIB failed:\n"
+    "Include directory: ${SIONLIB_INCLUDE_DIR}\n"
+    "Library directory: ${SIONLIB_LIBRARY_DIR}\n\n")
+endif(SIONLIB_FOUND)
+
+MESSAGE(STATUS "SIONLIB_LIBRARIES=${SIONLIB_LIBRARIES}")
+
+#set HAVE_SIONLIB for config.h
+set(HAVE_SIONLIB ${SIONLIB_FOUND})
+
+#add all sionlib related flags to ALL_PKG_FLAGS, this must happen regardless of a target using add_sionlib_flags
+if(SIONLIB_FOUND)
+  set_property(GLOBAL APPEND PROPERTY ALL_PKG_FLAGS "${SIONLIB_COMPILE_FLAGS}")
+  foreach(dir "${SIONLIB_INCLUDE_DIR}")
+    set_property(GLOBAL APPEND PROPERTY ALL_PKG_FLAGS "-I${dir}")
+  endforeach()
+endif()
\ No newline at end of file
diff --git a/cmake/ThirdPartyBoost.cmake b/cmake/ThirdPartyBoost.cmake
index 914282e494921ea8dd185f950a13bca2d77cdddb..dbcc16d6513e561328703f8d14610bd4fcedbf9e 100644
--- a/cmake/ThirdPartyBoost.cmake
+++ b/cmake/ThirdPartyBoost.cmake
@@ -8,8 +8,10 @@
 
 #If the user has not set BOOST_ROOT, look in a couple common places first.
 MESSAGE(STATUS "Searching for Boost:")
-SET(MIN_VER "1.56.0")
-SET(NEEDED_BOOST_LIBS thread iostreams filesystem system program_options regex)
+SET(MIN_VER "1.60.0")
+SET(NEEDED_BOOST_LIBS thread iostreams date_time filesystem system
+    program_options regex)
+SET(Boost_DEBUG 0)
 SET(Boost_NO_BOOST_CMAKE ON)
 IF( BOOST_ROOT )
     SET(Boost_NO_SYSTEM_PATHS ON)
@@ -164,7 +166,13 @@ IF (THIRDPARTY_BUILD_BOOST)
         STRING(TOUPPER ${BOOSTLIB} BOOSTLIB_UPPER)
         THIRDPARTY_LIBRARY(Boost_${BOOSTLIB_UPPER}_LIBRARY
             SHARED boost_${BOOSTLIB} DESCRIPTION "Boost ${BOOSTLIB} library")
+        THIRDPARTY_LIBRARY(Boost_${BOOSTLIB_UPPER}_LIBRARY_DEBUG
+            SHARED boost_${BOOSTLIB} DESCRIPTION "Boost ${BOOSTLIB} library, debug")
+        THIRDPARTY_LIBRARY(Boost_${BOOSTLIB_UPPER}_LIBRARY_RELEASE
+            SHARED boost_${BOOSTLIB} DESCRIPTION "Boost ${BOOSTLIB} library, release")
         MARK_AS_ADVANCED(Boost_${BOOSTLIB_UPPER}_LIBRARY)
+        MARK_AS_ADVANCED(Boost_${BOOSTLIB_UPPER}_LIBRARY_DEBUG)
+        MARK_AS_ADVANCED(Boost_${BOOSTLIB_UPPER}_LIBRARY_RELEASE)
         LIST(APPEND Boost_LIBRARIES ${Boost_${BOOSTLIB_UPPER}_LIBRARY})
     ENDFOREACH()
 
diff --git a/cmake/ThirdPartySIONLIB.cmake b/cmake/ThirdPartySIONLIB.cmake
new file mode 100755
index 0000000000000000000000000000000000000000..b9226bdbc75ea00ffd72b5868de7a75044b483a7
--- /dev/null
+++ b/cmake/ThirdPartySIONLIB.cmake
@@ -0,0 +1,33 @@
+########################################################################
+#
+# ThirdParty configuration for Nektar++
+#
+# SIONlib
+#
+########################################################################
+
+OPTION(NEKTAR_USE_SIONLIB
+    "Enable SIONlib I/O support." ON)
+
+IF (NEKTAR_USE_SIONLIB)
+    INCLUDE (CheckIncludeFiles)
+
+    IF (NOT NEKTAR_USE_MPI)
+        MESSAGE(FATAL_ERROR "SIONlib requires Nektar++ to be configured with NEKTAR_USE_MPI for MPI support.")
+    ENDIF()
+
+    # Try to find SIONlib package.
+    FIND_PACKAGE(SIONLIB)
+
+    IF (NOT SIONLIB_FOUND)
+        MESSAGE(FATAL_ERROR "SIONlib not detected.")
+    ENDIF()
+    
+    SET(BUILD_SIONLIB OFF)
+
+    SET(SIONLIB_CONFIG_INCLUDE_DIR ${SIONLIB_INCLUDE_DIR})
+    
+    MARK_AS_ADVANCED(SIONLIB_LIBRARIES)
+    MARK_AS_ADVANCED(SIONLIB_INCLUDE_DIR)
+    INCLUDE_DIRECTORIES(SYSTEM ${SIONLIB_INCLUDE_DIR})
+ENDIF()
diff --git a/library/Demos/LibUtilities/FieldIOBenchmarker.cpp b/library/Demos/LibUtilities/FieldIOBenchmarker.cpp
index c30f37ca6f1fa3752ef5d184d13c7a36404060bf..4c4bc9ed37055eaa3c4c96d8b426af9e17a79ca9 100644
--- a/library/Demos/LibUtilities/FieldIOBenchmarker.cpp
+++ b/library/Demos/LibUtilities/FieldIOBenchmarker.cpp
@@ -10,7 +10,7 @@
 //  Department of Aeronautics, Imperial College London (UK), and Scientific
 //  Computing and Imaging Institute, University of Utah (USA).
 //
-//  License for the specific language governing rights and laimitations under
+//  License for the specific language governing rights and limitations under
 //  Permission is hereby granted, free of charge, to any person obtaining a
 //  copy of this software and associated documentation files (the "Software"),
 //  to deal in the Software without restriction, including without limitation
@@ -29,37 +29,56 @@
 //  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 //  DEALINGS IN THE SOFTWARE.
 //
-//  Description: Measure the performance of FieldIO XML and HDF5 classes.
+//  Description: Measure the performance of FieldIO XML, HDF5 and SIONlib classes.
 //
 ////////////////////////////////////////////////////////////////////////////////
 
+#include <string>
+#include <boost/algorithm/string.hpp>
+#include <boost/program_options.hpp>
+#include <boost/bimap.hpp>
+#include <LibUtilities/BasicUtils/H5.h>
+#include <LibUtilities/BasicUtils/SIONlib.h>
 #include <LibUtilities/BasicUtils/FieldIO.h>
 #include <LibUtilities/BasicUtils/FieldIOXml.h>
+#include <LibUtilities/BasicUtils/FieldIOHdf5.h>
+#include <LibUtilities/BasicUtils/FieldIOSIONlib.h>
 #include <LibUtilities/BasicUtils/FileSystem.h>
 #include <LibUtilities/Communication/CommMpi.h>
-#include <boost/algorithm/string.hpp>
-#include <boost/program_options.hpp>
-#include <unordered_set>
-#include <string>
+#include <LibUtilities/BasicUtils/SessionReader.h>
 
+// Below, we'd like to use an unordered set for its faster lookup performance
+// However this is only available if C++11 is.
+//
+#if __cplusplus >= 201103L
+#include <unordered_set>
 typedef std::unordered_set<int> IntSet;
+#else
+#include <set>
+typedef std::set<int> IntSet;
+#endif
 
 using namespace Nektar;
 using namespace LibUtilities;
 
-namespace po    = boost::program_options;
+namespace po = boost::program_options;
 namespace berrc = boost::system::errc;
 
 typedef std::vector<FieldDefinitionsSharedPtr> DefVec;
 typedef std::vector<std::vector<NekDouble> > DatVec;
 
-/// Struct that contains experimental setup.
+typedef boost::shared_ptr<FieldIOXml> FieldIOXmlSharedPtr;
+typedef boost::shared_ptr<FieldIOHdf5> FieldIOHdf5SharedPtr;
+typedef boost::shared_ptr<FieldIOSIONlib> FieldIOSIONlibSharedPtr;
+
 struct Experiment
 {
     /// Test read (false) or write (true)
     bool write;
-    /// Test the HDF5 (true) or XML (false) reader/writer.
+    /// Test the HDF5 (true) or XML/SIONlib (false) reader/writer.
     bool hdf;
+    /// Test the SIONlib (true) or XML/HDF5 (false) reader/writer.
+    bool sionlib;
     /// If true, print additional debugging information.
     bool verbose;
     /// Number of writes to perform.
@@ -72,25 +91,38 @@ struct Experiment
     CommSharedPtr comm;
 };
 
-typedef std::vector<NekDouble> Results;
+IOSettingsSharedPtr sionlib_settings(new IOSettings);
+
+typedef std::vector<double> Results;
+
+const double MB = 1000000.0;
+const int AORTIC_ARCH = 0;
+const int RACING_CAR = 1;
 
+std::string SIONlib_GetIOBlocksPerChunk(int nprocs, int testid);
 Results TestRead(Experiment &exp);
 Results TestWrite(Experiment &exp);
-void PrintResults(Experiment &exp, Results &results);
+
 
 int main(int argc, char *argv[])
 {
     Experiment exp;
     exp.write   = false;
     exp.hdf     = false;
+    exp.sionlib = false;
     exp.verbose = false;
     exp.n       = 3;
     exp.comm    = GetCommFactory().CreateInstance("ParallelMPI", argc, argv);
 
+    sionlib_settings->insert( IOSettings::value_type("IOBlockSize", "65536") );
+    sionlib_settings->insert( IOSettings::value_type("IOBlocksPerChunk", SIONlib_GetIOBlocksPerChunk(exp.comm->GetSize(), AORTIC_ARCH)) );
+    sionlib_settings->insert( IOSettings::value_type("IOReadMode", "br,collective,collsize=-1") );
+    sionlib_settings->insert( IOSettings::value_type("IOWriteMode", "bw,collective,collectivemerge,collsize=-1") );
+        
     po::options_description desc("Available options");
     desc.add_options()("help,h", "Produce this help message.")(
         "mode,m", po::value<char>(),
-        "Choose r[ead] (default), x[ml write] or h[df5 write]")(
+        "Choose r[ead] (default), x[ml write] or h[df5 write] or s[sionlib write]")(
         "number,n", po::value<unsigned>(),
         "Number of iterations to perform, default 3")("verbose,v",
                                                       "Enable verbose mode.");
@@ -157,15 +189,22 @@ int main(int argc, char *argv[])
         switch (mode)
         {
             case 'r':
-                exp.write = false;
+                exp.write   = false;
                 break;
             case 'x':
-                exp.write = true;
-                exp.hdf   = false;
+                exp.write   = true;
+                exp.hdf     = false;
+                exp.sionlib = false;
                 break;
             case 'h':
-                exp.write = true;
-                exp.hdf   = true;
+                exp.write   = true;
+                exp.hdf     = true;
+                exp.sionlib = false;
+                break;
+            case 's':
+                exp.write   = true;
+                exp.hdf     = false;
+                exp.sionlib = true;
                 break;
             default:
                 std::cout << "Unrecognised mode: " << mode << std::endl;
@@ -194,10 +233,62 @@ int main(int argc, char *argv[])
         res = TestRead(exp);
     }
 
-    PrintResults(exp, res);
     exp.comm->Finalise();
 }
 
+
+/**
+ * @brief Return the number of file system blocks per SIONlib chunk.
+ *
+ * This result depends on the number of processes and on the size of the checkpoint
+ * file associated with the test case.
+ *
+ * @param nprocs  The number of MPI processes.
+ * @param testid  The test case --- this determines the size of the checkpoint file.
+ *
+ * @return string.
+ */
+std::string SIONlib_GetIOBlocksPerChunk(int nprocs, int testid)
+{
+    std::string bpc = "1.0";
+  
+    if (AORTIC_ARCH == testid)
+    {
+        switch (nprocs)
+        {
+            case 48:
+                bpc = "8.0";
+                break;
+            case 96:
+                bpc = "4.0";
+                break;
+            case 192:
+                bpc = "2.0";
+                break;
+        }
+    }
+    else if (RACING_CAR == testid)
+    {
+        switch (nprocs)
+        {
+            case 768:
+                bpc = "110.0";
+                break;   
+            case 1536:
+                bpc = "55.0";
+                break;
+            case 3072:
+                bpc = "29.0";
+                break;
+            case 6144:
+	        bpc = "16.0";
+                break;
+        }
+    }
+
+    return bpc;
+}
+    
 /**
  * @brief Read elemental IDs from XML field file format for this rank.
  *
@@ -274,91 +365,13 @@ Array<OneD, int> ReadIDsForThisRank(Experiment &exp, FieldIOSharedPtr fio)
 }
 
 /**
- * @brief Construct data for this rank based on element IDs and input data.
- *
- * Extract from @p inFieldDefs and @p inFieldData those elements which are in @p
- * ElementIDs and return them in the parameters @p outFieldDefs and @p
- * outFieldData.
+ * @brief  Get the IDs of the elements managed by this rank.
  *
- * @param inFieldDefs    Field definitions of input file
- * @param inFieldData    Input field data corresponding to the definitions
- * @param ElementIDs     Element IDs for this rank
- * @param outFieldDefs   Output field definitions containing this rank's data
- * @param outFieldData   Output field data corresponding to the definitions
- */
-void FilterDataForThisRank(const DefVec &inFieldDefs,
-                           const DatVec &inFieldData,
-                           Array<OneD, int> ElementIDs,
-                           DefVec &outFieldDefs,
-                           DatVec &outFieldData)
-{
-    // Create a set with all the IDs
-    IntSet IDs(ElementIDs.begin(), ElementIDs.end());
-
-    // Clear the output vectors
-    outFieldDefs.clear();
-    outFieldData.clear();
-
-    // Loop through all the loaded elements and copy over if in the requested
-    // set
-    DefVec::const_iterator inDefIt = inFieldDefs.begin();
-    DatVec::const_iterator inDatIt = inFieldData.begin();
-    for (; inDefIt != inFieldDefs.end(); ++inDefIt, ++inDatIt)
-    {
-        FieldDefinitionsSharedPtr inDef = *inDefIt;
-        // Use list to avoid endless reallocation.
-        std::list<unsigned int> elOut;
-        std::list<NekDouble> datOut;
-
-        unsigned dat_per_el = inDatIt->size() / inDef->m_elementIDs.size();
-
-        std::vector<unsigned int>::const_iterator elIt =
-            inDef->m_elementIDs.begin();
-        std::vector<NekDouble>::const_iterator datIt = inDatIt->begin();
-        for (; elIt != inDef->m_elementIDs.end(); ++elIt, datIt += dat_per_el)
-        {
-            if (IDs.find(*elIt) != IDs.end())
-            {
-                // Copy across element id
-                elOut.push_back(*elIt);
-                // and data
-                datOut.insert(datOut.end(), datIt, datIt + dat_per_el);
-            }
-        }
-        if (elOut.size())
-        {
-            // create the outFieldDefs
-            FieldDefinitionsSharedPtr defOut =
-                FieldDefinitionsSharedPtr(new FieldDefinitions(
-                    inDef->m_shapeType,
-                    std::vector<unsigned int>(elOut.begin(), elOut.end()),
-                    inDef->m_basis, inDef->m_uniOrder, inDef->m_numModes,
-                    inDef->m_fields, inDef->m_numHomogeneousDir,
-                    inDef->m_homogeneousLengths, inDef->m_homoStrips,
-                    inDef->m_homogeneousSIDs, inDef->m_homogeneousZIDs,
-                    inDef->m_homogeneousYIDs, inDef->m_points,
-                    inDef->m_pointsDef, inDef->m_numPoints,
-                    inDef->m_numPointsDef));
-            // Add to return
-            outFieldDefs.push_back(defOut);
-            // create the out data vector from our list
-            outFieldData.push_back(
-                std::vector<NekDouble>(datOut.begin(), datOut.end()));
-        }
-    }
-}
-
-/**
- * @brief Read all data in those files that this rank wants (and any other data
- * in them too).
+ * @param exp  Experimental setup.
  *
- * @param Exp  Experiment to be run
- * @param outFieldDefs   Output field definitions
- * @param outFieldData   Output data corresponding to the definitions
+ * @return IDs of the elements managed by this rank.
  */
-void ReadWholeFilesForThisRank(Experiment &exp,
-                               DefVec &outFieldDefs,
-                               DatVec &outFieldData)
+Array<OneD, int> ReadIDsForThisRank(Experiment& exp)
 {
     std::string ft = FieldIO::GetFileType(exp.dataSource, exp.comm);
     FieldIOSharedPtr fio =
@@ -366,43 +379,16 @@ void ReadWholeFilesForThisRank(Experiment &exp,
 
     if (fs::is_directory(exp.dataSource))
     {
-        Array<OneD, int> ElementIDs = ReadIDsForThisRank(exp, fio);
-        FieldMetaDataMap fieldmetadatamap;
-
-        // Load all the data from files that contain any of the IDs we want.
-        fio->Import(exp.dataSource, outFieldDefs, outFieldData,
-                    fieldmetadatamap, ElementIDs);
+        return ReadIDsForThisRank(exp, fio);
     }
-    else
+    else if (exp.hdf)
     {
-        fio->Import(exp.dataSource, outFieldDefs, outFieldData);
+        // ensure auto_selective
+        Array<OneD, int> ElementIDs;
+        return ElementIDs;
     }
-}
-
-/**
- * @brief Read only the data that this rank wants.
- *
- * @param exp           Experiment setup details
- * @param outFieldDefs  Resulting output field definitions for this rank.
- * @param outFieldData  Output field data for this rank.
- */
-void ReadDecomposed(Experiment &exp, DefVec &outFieldDefs, DatVec &outFieldData)
-{
-    std::string ft = FieldIO::GetFileType(exp.dataSource, exp.comm);
-    FieldIOSharedPtr fio =
-        GetFieldIOFactory().CreateInstance(ft, exp.comm, true);
-
-    Array<OneD, int> ElementIDs = ReadIDsForThisRank(exp, fio);
-    DefVec fileFieldDefs;
-    DatVec fileFieldData;
-    FieldMetaDataMap fieldmetadatamap;
 
-    // Load all the data from files that contain any of the IDs we want.
-    fio->Import(exp.dataSource, fileFieldDefs, fileFieldData, fieldmetadatamap,
-                ElementIDs);
-    // Filter it
-    FilterDataForThisRank(fileFieldDefs, fileFieldData, ElementIDs,
-                          outFieldDefs, outFieldData);
+    return NullInt1DArray;
 }
 
 /**
@@ -417,45 +403,60 @@ void ReadDecomposed(Experiment &exp, DefVec &outFieldDefs, DatVec &outFieldData)
  */
 Results TestRead(Experiment &exp)
 {
+    const std::string ft = FieldIO::GetFileType(exp.dataSource, exp.comm);
     if (exp.verbose)
     {
-        std::cout << "Beginning read experiment with " << exp.n << " loops."
-                  << std::endl;
+        std::cout << "Beginning read experiment..."
+                  << exp.n << " iterations." << std::endl;
         std::cout << "Determining file type... ";
+        std::cout << ft << std::endl;
     }
 
-    const std::string ft = FieldIO::GetFileType(exp.dataSource, exp.comm);
-    if (exp.verbose)
+    if (0 == ft.compare("Hdf5"))
+    {
+        exp.hdf = true;
+    }
+    else if (0 == ft.compare("SIONlib"))
     {
-        std::cout << ft << endl;
+        exp.sionlib = true;
     }
 
+    std::vector<FieldDefinitionsSharedPtr> fielddefs;
+    std::vector < std::vector<NekDouble> > fielddata;
+    FieldMetaDataMap fieldmetadatamap;
+    FieldIOSharedPtr fio = GetFieldIOFactory().CreateInstance(
+            ft, exp.comm, true);
+
+    if (exp.sionlib)
+    {
+        fio->InitFromBenchmarker(sionlib_settings);
+    }
+    
+    Array<OneD, int> ElementIDs = ReadIDsForThisRank(exp);
+        
     Results res(exp.n, 0.0);
     for (unsigned i = 0; i < exp.n; ++i)
     {
         if (exp.verbose)
         {
-            std::cout << "Test " << i << " of " << exp.n;
+            std::cout << "Test " << i+1 << " of " << exp.n;
         }
+                                
+        double t0 = MPI_Wtime();
+        unsigned long nRead = fio->Import(exp.dataSource,
+            fielddefs, fielddata,
+            fieldmetadatamap, ElementIDs);        
+        double t1 = MPI_Wtime() - t0;
 
-        std::vector<FieldDefinitionsSharedPtr> fielddefs;
-        std::vector<std::vector<NekDouble> > fielddata;
-        // Synchronise
-        exp.comm->Block();
-
-        NekDouble t0 = MPI_Wtime();
-
-        ReadWholeFilesForThisRank(exp, fielddefs, fielddata);
-
-        NekDouble t1 = MPI_Wtime();
-        t1 -= t0;
+        exp.comm->AllReduce(nRead, LibUtilities::ReduceSum);
 
         if (exp.verbose)
         {
-            std::cout << ": t = " << t1 << " s" << std::endl;
+            std::cout << ": read " << nRead/MB << " MB in "
+                      << t1 << " s." << std::endl;
         }
-
-        res[i] = t1;
+        
+        res[i] = (nRead/MB) / t1;
     }
     return res;
 }
@@ -472,112 +473,76 @@ Results TestRead(Experiment &exp)
  */
 Results TestWrite(Experiment &exp)
 {
+    const std::string ft = FieldIO::GetFileType(exp.dataSource, exp.comm);
     if (exp.verbose)
-        std::cout << "Reading in input: " << exp.dataSource << std::endl;
-
-    std::vector<FieldDefinitionsSharedPtr> fielddefs;
-    std::vector<std::vector<NekDouble> > fielddata;
-    ReadDecomposed(exp, fielddefs, fielddata);
-
-    std::string outtype;
-    if (exp.hdf)
-    {
-        outtype = "Hdf5";
-    }
-    else
     {
-        outtype = "Xml";
+        std::cout << "Beginning write experiment..."
+                  << exp.n << " iterations." << std::endl;
+        std::cout << "Determining file type... ";
+        std::cout << ft << std::endl;
     }
+    
+    std::vector<FieldDefinitionsSharedPtr> fielddefs;
+    std::vector < std::vector<NekDouble> > fielddata;
+    FieldMetaDataMap fieldmetadatamap;
+    FieldIOSharedPtr fio = GetFieldIOFactory().CreateInstance(
+            ft, exp.comm, true);
 
-    if (exp.verbose)
+    if (exp.sionlib)
     {
-        std::cout << "Beginning write (" << outtype << ") experiment with "
-                  << exp.n << " loops." << std::endl;
-        std::cout << "Writing to temp file: " << exp.dataDest << std::endl;
+        fio->InitFromBenchmarker(sionlib_settings);
     }
+    
+    Array<OneD, int> ElementIDs = ReadIDsForThisRank(exp);
+    fio->Import(exp.dataSource, fielddefs, fielddata,
+        fieldmetadatamap, ElementIDs);
+
+    exp.comm->Block();
 
     Results res(exp.n, 0);
+    fs::path specPath(exp.dataDest);
+
     for (unsigned i = 0; i < exp.n; ++i)
     {
         if (exp.verbose)
         {
-            std::cout << "Test " << i << " of " << exp.n << std::endl;
+            std::cout << "Test " << i+1 << " of " << exp.n;
         }
+        
+        double t0 = MPI_Wtime();
+        unsigned long nWritten = fio->Write(exp.dataDest, fielddefs, fielddata);
+        double t1 = MPI_Wtime() - t0;
 
-        // Synchronise - have to do this before removing any old data in case
-        // any ranks haven't closed their file yet.
-        exp.comm->Block();
-        // Remove any existing files
-        fs::path specPath(exp.dataDest);
-        try
+        exp.comm->AllReduce(nWritten, LibUtilities::ReduceSum);
+         
+        if (exp.verbose)
         {
-            fs::remove_all(specPath);
+            std::cout << ": written " << nWritten/MB << " MB in "
+                      << t1 << " s."  << std::endl;
         }
-        catch (fs::filesystem_error &e)
+        
+        res[i] = (nWritten/MB) / t1;
+        /*
+        if (exp.comm->GetRank() == 0)
         {
-            ASSERTL0(e.code().value() == berrc::no_such_file_or_directory,
-                     "Filesystem error: " + string(e.what()));
+            try
+            {
+                // remove any files that have just been written
+                // as part of test
+                fs::remove_all(specPath);  
+            }
+            catch (fs::filesystem_error& e)
+            {
+                ASSERTL0(
+                    e.code().value()
+                    == berrc::no_such_file_or_directory,
+                    "Filesystem error: " + string(e.what()));
+            }
         }
-
-        // Synchronise to make sure we're all at the same point.
+        */
         exp.comm->Block();
-
-        NekDouble t0 = MPI_Wtime();
-
-        FieldIOSharedPtr fio =
-            GetFieldIOFactory().CreateInstance(outtype, exp.comm, true);
-
-        fio->Write(exp.dataDest, fielddefs, fielddata);
-
-        NekDouble t1 = MPI_Wtime();
-        t1 -= t0;
-
-        if (exp.verbose)
-        {
-            std::cout << ": t = " << t1 << " s" << std::endl;
-        }
-
-        res[i] = t1;
     }
     return res;
 }
 
-/**
- * @brief Print out the results of the timing.
- *
- * Given an experimental setup @p exp and timing in @p results, this routine
- * prints the mean time per read/write.
- *
- * @param exp      Experimental setup.
- * @param results  Results of timing for @p exp
- */
-void PrintResults(Experiment &exp, Results &results)
-{
-    NekDouble sum   = 0.0;
-    NekDouble sumSq = 0.0;
-
-    for (Results::const_iterator it = results.begin(); it != results.end();
-         ++it)
-    {
-        NekDouble x = *it;
-        sum += x;
-        sumSq += x * x;
-    }
 
-    NekDouble mean = sum / exp.n;
-    // NekDouble var = sumSq / exp.n - mean*mean;
-    // NekDouble std = std::sqrt(var);
-
-    if (exp.comm->GetSize() > 1)
-    {
-        // Use all version since reduce to specified rank isn't wrapped.
-        exp.comm->AllReduce(mean, ReduceSum);
-        mean /= exp.comm->GetSize();
-    }
-
-    if (exp.comm->GetRank() == 0)
-    {
-        std::cout << "Mean: " << mean << std::endl;
-        // std::cout << "Std: " << std << std::endl;
-    }
-}
diff --git a/library/LibUtilities/BasicConst/NektarUnivTypeDefs.hpp b/library/LibUtilities/BasicConst/NektarUnivTypeDefs.hpp
index edfcf7c33008e44b15e18d0529e0e4adb2cd4444..0cbad4d21ab4221a4283ff26d60931c1e477304a 100644
--- a/library/LibUtilities/BasicConst/NektarUnivTypeDefs.hpp
+++ b/library/LibUtilities/BasicConst/NektarUnivTypeDefs.hpp
@@ -41,6 +41,7 @@
 
 namespace Nektar
 {
+    typedef unsigned char NekByte;
     typedef double NekDouble;
 
     typedef std::int32_t  NekInt;
diff --git a/library/LibUtilities/BasicUtils/FieldIO.cpp b/library/LibUtilities/BasicUtils/FieldIO.cpp
index 1e0867b2864f499d1cf87f90d18d64cd4dade2c6..8780f924877ccbc2ad020ea9b803c034c4469fb6 100644
--- a/library/LibUtilities/BasicUtils/FieldIO.cpp
+++ b/library/LibUtilities/BasicUtils/FieldIO.cpp
@@ -77,7 +77,8 @@ FieldIOFactory &GetFieldIOFactory()
 /// Enumerator for auto-detection of FieldIO types.
 enum FieldIOType {
     eXML,
-    eHDF5
+    eHDF5,
+    eSIONLIB
 };
 
 
@@ -97,7 +98,7 @@ const std::string FieldIO::GetFileType(const std::string &filename,
                                        CommSharedPtr comm)
 {
     FieldIOType ioType = eXML;
-    int size  = comm->GetSize();
+    int size = comm->GetSize();
     bool root = comm->TreatAsRankZero();
 
     if (size == 1 || root)
@@ -119,22 +120,44 @@ const std::string FieldIO::GetFileType(const std::string &filename,
         // Read first 8 bytes. If they correspond with magic bytes below it's an
         // HDF5 file. XML is potentially a nightmare with all the different
         // encodings so we'll just assume it's OK if it's not HDF.
-        const unsigned char magic[8] = {
+        const unsigned char hdf5_magic[8] = {
             0x89, 0x48, 0x44, 0x46, 0x0d, 0x0a, 0x1a, 0x0a};
 
+        const unsigned char sionlib_magic[8] = {
+            0x73, 0x69, 0x6f, 0x6e, 0x00, 0x00, 0x00, 0x00};
+
         std::ifstream datafile(datafilename.c_str(), ios_base::binary);
         ASSERTL0(datafile.good(), "Unable to open file: " + filename);
 
-        ioType = eHDF5;
+        bool is_hdf5 = true;
+        bool is_sionlib = true;
         for (unsigned i = 0; i < 8 && datafile.good(); ++i)
         {
             unsigned char byte = datafile.get();
-            if (byte != magic[i])
+            if (is_hdf5 && byte != hdf5_magic[i])
+            {
+                is_hdf5 = false;
+            }
+            if (is_sionlib && byte != sionlib_magic[i])
+            {
+                is_sionlib = false;
+            }
+
+            if (!is_hdf5 && !is_sionlib)
             {
-                ioType = eXML;
                 break;
             }
         }
+
+        ioType = eXML;
+        if (is_hdf5)
+        {
+            ioType = eHDF5;
+        }
+        else if (is_sionlib)
+        {
+            ioType = eSIONLIB;
+        }
     }
 
     if (size > 1)
@@ -153,6 +176,10 @@ const std::string FieldIO::GetFileType(const std::string &filename,
     {
         iofmt = "Hdf5";
     }
+    else if (ioType == eSIONLIB)
+    {
+        iofmt = "SIONlib";
+    }
     else
     {
         // Error
@@ -187,10 +214,14 @@ FieldIOSharedPtr FieldIO::CreateDefault(
         iofmt = session->GetCmdLineArgument<std::string>("io-format");
     }
 
-    return GetFieldIOFactory().CreateInstance(
+    FieldIOSharedPtr fieldio = GetFieldIOFactory().CreateInstance(
         iofmt,
         session->GetComm(),
         session->GetSharedFilesystem());
+
+    fieldio->Init(session);
+
+    return fieldio;
 }
 
 /**
@@ -210,10 +241,15 @@ FieldIOSharedPtr FieldIO::CreateForFile(
 {
     const std::string iofmt =
         FieldIO::GetFileType(filename, session->GetComm());
-    return GetFieldIOFactory().CreateInstance(
+    
+    FieldIOSharedPtr fieldio = GetFieldIOFactory().CreateInstance(
         iofmt,
         session->GetComm(),
         session->GetSharedFilesystem());
+
+    fieldio->Init(session);
+
+    return fieldio;
 }
 
 /**
@@ -226,12 +262,13 @@ FieldIOSharedPtr FieldIO::CreateForFile(
  * @param fielddata     Binary field data that stores the output corresponding
  *                      to @p fielddefs.
  * @param fieldinfomap  Associated field metadata map.
+ * @return The number of bytes written.
  */
-void Write(const std::string &outFile,
-           std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-           std::vector<std::vector<NekDouble> > &fielddata,
-           const FieldMetaDataMap &fieldinfomap,
-           const bool backup)
+uint64_t Write(const std::string &outFile,
+               std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+               std::vector<std::vector<NekDouble> > &fielddata,
+               const FieldMetaDataMap &fieldinfomap,
+               const bool backup)
 {
 #ifdef NEKTAR_USE_MPI
     int size;
@@ -253,7 +290,7 @@ void Write(const std::string &outFile,
 #endif
     CommSharedPtr c    = GetCommFactory().CreateInstance("Serial", 0, 0);
     FieldIOSharedPtr f = GetFieldIOFactory().CreateInstance("Xml", c, false);
-    f->Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
+    return f->Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
 }
 
 /**
@@ -270,8 +307,9 @@ void Write(const std::string &outFile,
  * @param ElementIDs    Element IDs that lie on this processor, which can be
  *                      optionally supplied to avoid reading the entire file on
  *                      each processor.
+ * @return The number of bytes read.
  */
-LIB_UTILITIES_EXPORT void Import(
+LIB_UTILITIES_EXPORT uint64_t Import(
     const std::string &infilename,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     std::vector<std::vector<NekDouble> > &fielddata,
@@ -299,7 +337,7 @@ LIB_UTILITIES_EXPORT void Import(
     CommSharedPtr c    = GetCommFactory().CreateInstance("Serial", 0, 0);
     const std::string iofmt = FieldIO::GetFileType(infilename, c);
     FieldIOSharedPtr f = GetFieldIOFactory().CreateInstance(iofmt, c, false);
-    f->Import(infilename, fielddefs, fielddata, fieldinfomap, ElementIDs);
+    return f->Import(infilename, fielddefs, fielddata, fieldinfomap, ElementIDs);
 }
 
 /**
@@ -393,7 +431,7 @@ std::string FieldIO::SetUpOutput(const std::string outname, bool perRank, bool b
     ASSERTL0(!outname.empty(), "Empty path given to SetUpOutput()");
 
     int nprocs = m_comm->GetSize();
-    bool root  = m_comm->TreatAsRankZero();
+    bool root = m_comm->TreatAsRankZero();
 
     // Path to output: will be directory if parallel, normal file if
     // serial.
@@ -501,7 +539,7 @@ std::string FieldIO::SetUpOutput(const std::string outname, bool perRank, bool b
 
     if (root)
     {
-        cout << "Writing: " << specPath;
+        //cout << "Writing: " << specPath << std::endl;
     }
 
     // serial processing just add ending.
diff --git a/library/LibUtilities/BasicUtils/FieldIO.h b/library/LibUtilities/BasicUtils/FieldIO.h
index 446b1f966183662f798b7d94dc9a6a0df6add8ab..6aedf4616a3a9a5d2fff645a80aab1e77ee3fc43 100644
--- a/library/LibUtilities/BasicUtils/FieldIO.h
+++ b/library/LibUtilities/BasicUtils/FieldIO.h
@@ -53,6 +53,9 @@ namespace LibUtilities
 typedef std::map<std::string, std::string> FieldMetaDataMap;
 static FieldMetaDataMap NullFieldMetaDataMap;
 
+typedef std::map<std::string, std::string> IOSettings;
+typedef std::shared_ptr<IOSettings> IOSettingsSharedPtr;
+ 
 /**
  * @brief Base class for writing hierarchical data (XML or HDF5).
  */
@@ -179,13 +182,14 @@ struct FieldDefinitions
 
 typedef std::shared_ptr<FieldDefinitions> FieldDefinitionsSharedPtr;
 
-LIB_UTILITIES_EXPORT void Write(
+LIB_UTILITIES_EXPORT uint64_t Write(
     const std::string &outFile,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     std::vector<std::vector<NekDouble> > &fielddata,
     const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
     const bool backup = false);
-LIB_UTILITIES_EXPORT void Import(
+
+LIB_UTILITIES_EXPORT uint64_t Import(
     const std::string &infilename,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     std::vector<std::vector<NekDouble> > &fielddata = NullVectorNekDoubleVector,
@@ -230,14 +234,20 @@ public:
     {
     }
 
-    LIB_UTILITIES_EXPORT inline void Write(
+    LIB_UTILITIES_EXPORT inline void Init(
+        const LibUtilities::SessionReaderSharedPtr session);
+
+    LIB_UTILITIES_EXPORT inline void InitFromBenchmarker(
+        const LibUtilities::IOSettingsSharedPtr iosettings);
+
+    LIB_UTILITIES_EXPORT inline uint64_t Write(
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
         const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
         const bool backup = false);
 
-    LIB_UTILITIES_EXPORT inline void Import(
+    LIB_UTILITIES_EXPORT inline uint64_t Import(
         const std::string &infilename,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata =
@@ -249,6 +259,12 @@ public:
         const std::string &filename,
         FieldMetaDataMap  &fieldmetadatamap);
 
+    LIB_UTILITIES_EXPORT void ImportMultiFldFileIDs(
+        const std::string &inFile,
+        std::vector<std::string> &fileNames,
+        std::vector<std::vector<unsigned int> > &elementList,
+        FieldMetaDataMap &fieldmetadatamap);
+
     LIB_UTILITIES_EXPORT static const std::string GetFileType(
         const std::string &filename, CommSharedPtr comm);
     LIB_UTILITIES_EXPORT virtual const std::string &GetClassName() const = 0;
@@ -282,8 +298,16 @@ protected:
     LIB_UTILITIES_EXPORT std::string SetUpOutput(
         const std::string outname, bool perRank, bool backup = false);
 
+    /// @copydoc FieldIO::Init
+    LIB_UTILITIES_EXPORT virtual void v_Init(
+        const LibUtilities::SessionReaderSharedPtr session) {};
+
+    /// @copydoc FieldIO::InitFromBenchmarker
+    LIB_UTILITIES_EXPORT virtual void v_InitFromBenchmarker(
+        const LibUtilities::IOSettingsSharedPtr iosettings) {};
+    
     /// @copydoc FieldIO::Write
-    LIB_UTILITIES_EXPORT virtual void v_Write(
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Write(
         const std::string                      &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> >   &fielddata,
@@ -291,7 +315,7 @@ protected:
         const bool                              backup = false) = 0;
 
     /// @copydoc FieldIO::Import
-    LIB_UTILITIES_EXPORT virtual void v_Import(
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Import(
         const std::string &infilename,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> >
@@ -302,10 +326,37 @@ protected:
     /// @copydoc FieldIO::ImportFieldMetaData
     LIB_UTILITIES_EXPORT virtual DataSourceSharedPtr v_ImportFieldMetaData(
         const std::string &filename, FieldMetaDataMap &fieldmetadatamap) = 0;
+
+    LIB_UTILITIES_EXPORT virtual void v_ImportMultiFldFileIDs(
+        const std::string &inFile,
+        std::vector<std::string> &fileNames,
+        std::vector<std::vector<unsigned int> > &elementList,
+        FieldMetaDataMap &fieldmetadatamap) {};
 };
 
 typedef std::shared_ptr<FieldIO> FieldIOSharedPtr;
 
+/**
+ * @brief Allow the FieldIO class an opportunity for further initialisation.
+ *
+ * @param pointer to session configuration
+ */
+inline void FieldIO::Init(const LibUtilities::SessionReaderSharedPtr session)
+{
+    return v_Init(session);
+}
+
+/**
+ * @brief Allow the FieldIO class an opportunity for further initialisation
+ * from outside a usual Nektar++ session, e.g., a benchmarker.
+ *
+ * @param pointer to stl map containing extra io settings
+ */
+inline void FieldIO::InitFromBenchmarker(const LibUtilities::IOSettingsSharedPtr iosettings)
+{
+    return v_InitFromBenchmarker(iosettings);
+}
+
 /**
  * @brief Write out the field information to the file @p outFile.
  *
@@ -314,14 +365,15 @@ typedef std::shared_ptr<FieldIO> FieldIOSharedPtr;
  * @param fielddata     Binary field data that stores the output corresponding
  *                      to @p fielddefs.
  * @param fieldinfomap  Associated field metadata map.
+ * @return The number of bytes written.
  */
-inline void FieldIO::Write(const std::string                      &outFile,
-                           std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                           std::vector<std::vector<NekDouble> > &fielddata,
-                           const FieldMetaDataMap                 &fieldinfomap,
-                           const bool                              backup)
+inline uint64_t FieldIO::Write(const std::string             &outFile,
+                               std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                               std::vector<std::vector<NekDouble> >   &fielddata,
+                               const FieldMetaDataMap                 &fieldinfomap,
+                               const bool                              backup)
 {
-    v_Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
+    return v_Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
 }
 
 /**
@@ -336,14 +388,15 @@ inline void FieldIO::Write(const std::string                      &outFile,
  * @param ElementIDs    Element IDs that lie on this processor, which can be
  *                      optionally supplied to avoid reading the entire file on
  *                      each processor.
+ * @return The number of bytes read.
  */
-inline void FieldIO::Import(const std::string                      &infilename,
-                            std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                            std::vector<std::vector<NekDouble> >   &fielddata,
-                            FieldMetaDataMap                       &fieldinfo,
-                            const Array<OneD, int>                 &ElementIDs)
+inline uint64_t FieldIO::Import(const std::string             &infilename,
+                                std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                                std::vector<std::vector<NekDouble> >   &fielddata,
+                                FieldMetaDataMap                       &fieldinfo,
+                                const Array<OneD, int>                 &ElementIDs)
 {
-    v_Import(infilename, fielddefs, fielddata, fieldinfo, ElementIDs);
+    return v_Import(infilename, fielddefs, fielddata, fieldinfo, ElementIDs);
 }
 
 /**
@@ -359,6 +412,24 @@ inline DataSourceSharedPtr FieldIO::ImportFieldMetaData(
     return v_ImportFieldMetaData(filename, fieldmetadatamap);
 }
 
+
+/**
+ * @brief Read file containing element ID to partition mapping.
+ *
+ * @param inFile             Input multi-field file name.
+ * @param fileNames          List of partition filenames.
+ * @param elementList        Vector of element IDs that lie on each process.
+ * @param fieldmetadatamap   Field metadata map that is read from @p inFile.
+ */
+inline void FieldIO::ImportMultiFldFileIDs(
+        const std::string &inFile,
+        std::vector<std::string> &fileNames,
+        std::vector<std::vector<unsigned int> > &elementList,
+        FieldMetaDataMap &fieldmetadatamap)
+{
+    v_ImportMultiFldFileIDs(inFile, fileNames, elementList, fieldmetadatamap);
+}
+ 
 }
 }
 #endif
diff --git a/library/LibUtilities/BasicUtils/FieldIOHdf5.cpp b/library/LibUtilities/BasicUtils/FieldIOHdf5.cpp
index b35ce6a4325a4d7d56ee5449ffafcab617c46e3f..0a6b347f6cd48be2bacf5034f6b19d4a5812affa 100644
--- a/library/LibUtilities/BasicUtils/FieldIOHdf5.cpp
+++ b/library/LibUtilities/BasicUtils/FieldIOHdf5.cpp
@@ -172,22 +172,26 @@ FieldIOHdf5::FieldIOHdf5(LibUtilities::CommSharedPtr pComm,
  * @param fielddefs         Input field definitions.
  * @param fielddata         Input field data.
  * @param fieldmetadatamap  Field metadata.
+ * @return The number of bytes written.
  */
-void FieldIOHdf5::v_Write(const std::string &outFile,
-                          std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                          std::vector<std::vector<NekDouble> > &fielddata,
-                          const FieldMetaDataMap &fieldmetadatamap,
-                          const bool backup)
+uint64_t FieldIOHdf5::v_Write(const std::string &outFile,
+                              std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                              std::vector<std::vector<NekDouble> > &fielddata,
+                              const FieldMetaDataMap &fieldmetadatamap,
+                              const bool backup)
 {
     std::stringstream prfx;
     prfx << m_comm->GetRank() << ": FieldIOHdf5::v_Write(): ";
+
+    uint64_t nWritten = 0;
+
     double tm0 = 0.0, tm1 = 0.0;
 
     if (m_comm->GetRank() == 0)
     {
         tm0 = m_comm->Wtime();
     }
-
+    
     SetUpOutput(outFile, false, backup);
 
     // We make a number of assumptions in this code:
@@ -200,10 +204,13 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
     ASSERTL1(fielddefs.size() == fielddata.size(),
              prfx.str() + "fielddefs and fielddata have incompatible lengths.");
 
-    size_t nFields    = fielddefs.size();
-    size_t nMaxFields = nFields;
+    std::size_t nFields    = fielddefs.size();
+    std::size_t nMaxFields = nFields;
+    std::size_t nMinFields = nFields;
     m_comm->AllReduce(nMaxFields, LibUtilities::ReduceMax);
-
+    m_comm->AllReduce(nMinFields, LibUtilities::ReduceMin);
+    //cout << prfx.str() << "nFields " << nFields << endl;
+    
     int root_rank = -1;
     bool amRoot = false;
     LibUtilities::CommSharedPtr max_fields_comm;
@@ -247,7 +254,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
     int homDim = -1;
     int varOrder = 0;
 
-    for (int f = 0; f < nFields; ++f)
+    for (std::size_t f = 0; f < nFields; ++f)
     {
         if (!fielddefs[f]->m_uniOrder)
         {
@@ -261,7 +268,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
     // Calculate the total number of elements handled by this MPI process and
     // the total number of bytes required to store the elements. Base the name
     // of each field on the hash of the field definition.
-    for (int f = 0; f < nFields; ++f)
+    for (std::size_t f = 0; f < nFields; ++f)
     {
         ASSERTL1(fielddata[f].size() > 0,
                  prfx.str() +
@@ -283,13 +290,13 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
         // Hash the field specification
         std::stringstream hashStream;
         std::size_t nSubFields = fielddefs[f]->m_fields.size();
-        for (int sf = 0; sf < nSubFields; ++sf)
+        for (std::size_t sf = 0; sf < nSubFields; ++sf)
         {
             hashStream << fielddefs[f]->m_fields[sf];
         }
 
         nSubFields = fielddefs[f]->m_basis.size();
-        for (int sf = 0; sf < nSubFields; ++sf)
+        for (std::size_t sf = 0; sf < nSubFields; ++sf)
         {
             hashStream << fielddefs[f]->m_basis[sf];
         }
@@ -326,7 +333,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
         {
             nSubFields = fielddefs[f]->m_homogeneousLengths.size();
             homoLengths[f].resize(nSubFields);
-            for (int sf = 0; sf < nSubFields; ++sf)
+            for (std::size_t sf = 0; sf < nSubFields; ++sf)
             {
                 NekDouble len = fielddefs[f]->m_homogeneousLengths[sf];
                 hashStream << len;
@@ -339,7 +346,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
                 homoYIDs[f].resize(nSubFields);
                 decomps[f * MAX_DCMPS + HOMY_DCMP_IDX] = nSubFields;
                 cnts[HOMY_CNT_IDX] += nSubFields;
-                for (int sf = 0; sf < nSubFields; ++sf)
+                for (std::size_t sf = 0; sf < nSubFields; ++sf)
                 {
                     homoYIDs[f][sf] = fielddefs[f]->m_homogeneousYIDs[sf];
                 }
@@ -351,7 +358,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
                 homoZIDs[f].resize(nSubFields);
                 decomps[f * MAX_DCMPS + HOMZ_DCMP_IDX] = nSubFields;
                 cnts[HOMZ_CNT_IDX] += nSubFields;
-                for (int sf = 0; sf < nSubFields; ++sf)
+                for (std::size_t sf = 0; sf < nSubFields; ++sf)
                 {
                     homoZIDs[f][sf] = fielddefs[f]->m_homogeneousZIDs[sf];
                 }
@@ -363,7 +370,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
                 homoSIDs[f].resize(nSubFields);
                 decomps[f * MAX_DCMPS + HOMS_DCMP_IDX] = nSubFields;
                 cnts[HOMS_CNT_IDX] += nSubFields;
-                for (int sf = 0; sf < nSubFields; ++sf)
+                for (std::size_t sf = 0; sf < nSubFields; ++sf)
                 {
                     homoSIDs[f][sf] = fielddefs[f]->m_homogeneousSIDs[sf];
                 }
@@ -446,6 +453,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
 
         // Record file format version as attribute in main group.
         root->SetAttribute("FORMAT_VERSION", FORMAT_VERSION);
+        nWritten += sizeof(FORMAT_VERSION);
 
         // Calculate the indexes to be used by each MPI process when reading the
         // IDS and DATA datasets
@@ -571,7 +579,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
 
     for (int n = 0; n < m_comm->GetSize(); ++n)
     {
-        for (int i = 0; i < nMaxFields; ++i)
+        for (std::size_t i = 0; i < nMaxFields; ++i)
         {
             uint64_t hash = all_hashes[n*nMaxFields + i];
 
@@ -586,9 +594,10 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
     }
 
     // Having constructed the map, go ahead and write the attributes out.
-    for (auto &sIt : writingProcs)
+    map<int, std::vector<uint64_t> >::iterator sIt;
+    for (sIt = writingProcs.begin(); sIt != writingProcs.end(); sIt++)
     {
-        int rank = sIt.first;
+        int rank = sIt->first;
 
         // Write out this rank's groups.
         if (m_comm->GetRank() == rank)
@@ -605,32 +614,41 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
 
             // Write a HDF5 group for each field
             hashToProc.clear();
-            for (int i = 0; i < sIt.second.size(); ++i)
+            for (std::size_t i = 0; i < sIt->second.size(); ++i)
             {
-                for (int f = 0; f < nFields; ++f)
+                for (std::size_t f = 0; f < nFields; ++f)
                 {
-                    if (sIt.second[i] !=
+                    if (sIt->second[i] !=
                         all_hashes[m_comm->GetRank() * nMaxFields + f] ||
-                        hashToProc.find(sIt.second[i]) != hashToProc.end())
+                        hashToProc.find(sIt->second[i]) != hashToProc.end())
                     {
                         continue;
                     }
 
-                    hashToProc.insert(sIt.second[i]);
+                    hashToProc.insert(sIt->second[i]);
 
                     // Just in case we've already written this
                     H5::GroupSharedPtr field_group =
                         root->CreateGroup(fieldNames[f]);
                     ASSERTL1(field_group,
                              prfx.str() + "cannot create field group.");
+
                     field_group->SetAttribute("FIELDS", fielddefs[f]->m_fields);
                     field_group->SetAttribute("BASIS", fielddefs[f]->m_basis);
                     field_group->SetAttribute("SHAPE", shapeStrings[f]);
 
+                    for (std::size_t j = 0; j < fielddefs[f]->m_fields.size(); ++j)
+                    {
+                        nWritten += fielddefs[f]->m_fields[j].size();
+                    }
+                    nWritten += fielddefs[f]->m_basis.size()*sizeof(LibUtilities::BasisType);
+                    nWritten += shapeStrings[f].size();
+
                     if (homoLengths[f].size() > 0)
                     {
                         field_group->SetAttribute("HOMOGENEOUSLENGTHS",
                                                   homoLengths[f]);
+                        nWritten += homoLengths[f].size()*sizeof(NekDouble);
                     }
 
                     // If the field has only uniform order, we write the order
@@ -641,12 +659,14 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
                     {
                         field_group->SetAttribute("NUMMODESPERDIR",
                                                   numModesPerDirUni[f]);
+                        nWritten += numModesPerDirUni[f].size();
                     }
                     else
                     {
                         std::string numModesPerDir = "MIXORDER";
                         field_group->SetAttribute("NUMMODESPERDIR",
                                                   numModesPerDir);
+                        nWritten += numModesPerDir.size();
                     }
                 }
             }
@@ -680,6 +700,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
 
         decomps_fspace->SelectRange(0, all_decomps.size());
         decomps_dset->Write(all_decomps, decomps_fspace, writeSR);
+        nWritten += all_decomps.size()*sizeof(uint64_t);
     }
 
     // Initialise the dataset indexes for all MPI processes
@@ -693,15 +714,11 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
 
     // Set properties for parallel file access (if we're in parallel)
     H5::PListSharedPtr parallelProps = H5::PList::Default();
-    H5::PListSharedPtr writePL = H5::PList::Default();
     if (m_comm->GetSize() > 1)
     {
         // Use MPI/O to access the file
         parallelProps = H5::PList::FileAccess();
         parallelProps->SetMpio(m_comm);
-        // Use collective IO
-        writePL = H5::PList::DatasetXfer();
-        writePL->SetDxMpioCollective();
     }
 
     // Reopen the file
@@ -762,107 +779,140 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
         ASSERTL1(homs_fspace, prfx.str() + "cannot open HOMOGENEOUSSIDS filespace.");
     }
 
-    // Write the data
-    for (int f = 0; f < nFields; ++f)
+    
+    // we use independent (non-collective) writes for all datasets that are actually part of the field definition,
+    // i.e., {order_dset,homz_dset,homy_dset,homs_dset}; otherwise, we would have to assume that all fields
+    // were defined such that they all used the same selection of field definition datasets
+    nWritten += WriteFieldDataInd(nFields, order_fspace, order_dset, order_i, numModesPerDirVar);
+    nWritten += WriteFieldDataInd(nFields, homy_fspace, homy_dset, homy_i, homoYIDs);
+    nWritten += WriteFieldDataInd(nFields, homz_fspace, homz_dset, homz_i, homoZIDs);
+    nWritten += WriteFieldDataInd(nFields, homs_fspace, homs_dset, homs_i, homoSIDs);
+    
+    // write all the element IDs and element data collectively, taking into account the fact that
+    // different processes maybe handling different numbers of fields
+    std::vector<std::vector<unsigned int> > elemIDs(nFields);
+    for (std::size_t f = 0; f < nFields; ++f)
     {
-        // write the element ids
-        std::size_t nFieldElems = fielddefs[f]->m_elementIDs.size();
-        ids_fspace->SelectRange(ids_i, nFieldElems);
-        ids_dset->Write(fielddefs[f]->m_elementIDs, ids_fspace, writePL);
-        ids_i += nFieldElems;
-
-        // write the element values
-        std::size_t nFieldVals = fielddata[f].size();
-        data_fspace->SelectRange(data_i, nFieldVals);
-        data_dset->Write(fielddata[f], data_fspace, writePL);
-        data_i += nFieldVals;
+        elemIDs[f] = fielddefs[f]->m_elementIDs;
     }
+    nWritten += WriteFieldData(nMinFields, nFields, ids_fspace, ids_dset, ids_i, elemIDs);
+    nWritten += WriteFieldData(nMinFields, nFields, data_fspace, data_dset, data_i, fielddata);
+            
+    
+    m_comm->Block();
 
-    if (order_dset)
-    {
-        for (int f = 0; f < nFields; ++f)
-        {
-            std::size_t nOrders = numModesPerDirVar[f].size();
-            order_fspace->SelectRange(order_i, nOrders);
-            order_dset->Write(numModesPerDirVar[f], order_fspace, writePL);
-            order_i += nOrders;
-        }
-    }
+    // all data has been written
+    //if (m_comm->GetRank() == 0)
+    //{
+        //tm1 = m_comm->Wtime();
+        //cout << " (" << tm1 - tm0 << "s, HDF5)" << endl;
+    //}
 
-    if (homy_dset)
-    {
-        for (int f = 0; f < nFields; ++f)
-        {
-            std::size_t nYIDs = homoYIDs[f].size();
-            homy_fspace->SelectRange(homy_i, nYIDs);
-            homy_dset->Write(homoYIDs[f], homy_fspace, writePL);
-            homy_i += nYIDs;
-        }
-    }
+    return nWritten;
+}
 
-    if (homz_dset)
-    {
-        for (int f = 0; f < nFields; ++f)
-        {
-            std::size_t nZIDs = homoZIDs[f].size();
-            homz_fspace->SelectRange(homz_i, nZIDs);
-            homz_dset->Write(homoZIDs[f], homz_fspace, writePL);
-            homz_i += nZIDs;
-        }
-    }
+/**
+ * @brief Write the contributions from a number of fields to a specific dataset.
+ *
+ * The data are written independently.
+ *
+ * @param nFields  the number of fields handled by this process
+ * @param fspace   hdf5 file space
+ * @param dset     hdf5 data set
+ * @param data_i   starting offset into hdf5 data set for this process
+ * @param data     vector of data vectors (one for each field) to be written to data set
+ * @return The number of bytes written.
+ */
+template <class T>
+uint64_t FieldIOHdf5::WriteFieldDataInd(std::size_t nFields,
+                                        H5::DataSpaceSharedPtr &fspace, H5::DataSetSharedPtr &dset,
+                                        uint64_t data_i, std::vector<std::vector<T> > &data)
+{
+    if (!fspace || !dset) return 0;
 
-    if (homs_dset)
+    std::size_t nDataItems = 0;
+    uint64_t nWritten = 0;
+
+    H5::PListSharedPtr writePL = H5::PList::DatasetXfer();
+    writePL->SetDxMpioIndependent();
+
+    for (std::size_t f = 0; f < nFields; ++f)
     {
-        for (int f = 0; f < nFields; ++f)
+        nDataItems = data[f].size();
+        if (nDataItems > 0)
         {
-            std::size_t nSIDs = homoSIDs[f].size();
-            homs_fspace->SelectRange(homs_i, nSIDs);
-            homs_dset->Write(homoSIDs[f], homs_fspace, writePL);
-            homs_i += nSIDs;
+            fspace->SelectRange(data_i, nDataItems);
+            dset->Write(data[f], fspace, writePL);
+            data_i += nDataItems;
+            nWritten += nDataItems*sizeof(T);
         }
     }
 
-    for (int f = nFields; f < nMaxFields; ++f)
-    {
-        // this MPI process is handling fewer than nMaxFields fields
-        // so, since this is a collective operation
-        // just rewrite the element ids and values of the last field
-        ids_dset->Write(
-            fielddefs[nFields - 1]->m_elementIDs, ids_fspace, writePL);
-        data_dset->Write(fielddata[nFields - 1], data_fspace, writePL);
-
-        if (order_dset)
-        {
-            order_dset->Write(numModesPerDirVar[nFields - 1],
-                              order_fspace, writePL);
-        }
+    return nWritten;
+}
 
-        if (homy_dset)
-        {
-            homy_dset->Write(homoYIDs[nFields - 1], homy_fspace, writePL);
-        }
+/**
+ * @brief Write the contributions from a number of fields to a specific dataset.
+ *
+ * The data are written collectively; hence, this method allows for the fact that
+ * different processes might be handling different numbers of fields.
+ *
+ * @param nMinFields the lowest number of fields handled by any process
+ * @param nFields    the number of fields handled by this process
+ * @param fspace     hdf5 file space
+ * @param dset       hdf5 data set
+ * @param data_i     starting offset into hdf5 data set for this process
+ * @param data       vector of data vectors (one for each field) to be written to data set
+ * @return The number of bytes written.
+ */
+template <class T>
+uint64_t FieldIOHdf5::WriteFieldData(std::size_t nMinFields, std::size_t nFields,
+                                     H5::DataSpaceSharedPtr &fspace, H5::DataSetSharedPtr &dset,
+                                     uint64_t data_i, std::vector<std::vector<T> > &data)
+{
+    if (!fspace || !dset) return 0;
 
-        if (homz_dset)
-        {
-            homz_dset->Write(homoZIDs[nFields - 1], homz_fspace, writePL);
-        }
+    bool concatenate_last_fields = nMinFields < nFields;
+    std::size_t nFirstFields = concatenate_last_fields ? nMinFields-1 : nFields;
+    std::size_t nDataItems = 0;
+    std::size_t f = 0;
+    uint64_t nWritten = 0;
 
-        if (homs_dset)
-        {
-            homs_dset->Write(homoSIDs[nFields - 1], homs_fspace, writePL);
-        }
+    H5::PListSharedPtr writePL = H5::PList::DatasetXfer();
+    writePL->SetDxMpioCollective();
+
+    for (; f < nFirstFields; ++f)
+    {
+        nDataItems = data[f].size();
+        fspace->SelectRange(data_i, nDataItems);
+        dset->Write(data[f], fspace, writePL);
+        data_i += nDataItems;
+        nWritten += nDataItems*sizeof(T);
     }
 
-    m_comm->Block();
+    if (!concatenate_last_fields) return nWritten;
 
-    // all data has been written
-    if (m_comm->GetRank() == 0)
+    std::vector<T> last_data;
+
+    nDataItems = data[f].size();
+    fspace->SelectRange(H5S_SELECT_SET, data_i, nDataItems);
+    last_data.insert(last_data.end(), data[f].begin(), data[f].end());
+    data_i += nDataItems;
+    f++;
+        
+    for (; f < nFields; ++f)
     {
-        tm1 = m_comm->Wtime();
-        cout << " (" << tm1 - tm0 << "s, HDF5)" << endl;
+        nDataItems = data[f].size();
+        fspace->SelectRange(H5S_SELECT_OR, data_i, nDataItems);
+        last_data.insert(last_data.end(), data[f].begin(), data[f].end());
+        data_i += nDataItems;
     }
+
+    dset->Write(last_data, fspace, writePL);
+    return nWritten + last_data.size()*sizeof(T);
 }
 
+
 /**
  * @brief Import a HDF5 format file.
  *
@@ -874,21 +924,24 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
  *                          this rank. The resulting field definitions will only
  *                          contain data for the element IDs specified in this
  *                          array.
+ * @return The number of bytes read.
  */
-void FieldIOHdf5::v_Import(const std::string &infilename,
-                           std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                           std::vector<std::vector<NekDouble> > &fielddata,
-                           FieldMetaDataMap &fieldinfomap,
-                           const Array<OneD, int> &ElementIDs)
+uint64_t FieldIOHdf5::v_Import(const std::string &infilename,
+                               std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                               std::vector<std::vector<NekDouble> > &fielddata,
+                               FieldMetaDataMap &fieldinfomap,
+                               const Array<OneD, int> &ElementIDs)
 {
     std::stringstream prfx;
+    prfx << m_comm->GetRank() << ": FieldIOHdf5::v_Import(): ";
     int nRanks = m_comm->GetSize();
 
+    uint64_t nRead = 0;
+    
     // Set properties for parallel file access (if we're in parallel)
     H5::PListSharedPtr parallelProps = H5::PList::Default();
     H5::PListSharedPtr readPL = H5::PList::Default();
-    H5::PListSharedPtr readPLInd = H5::PList::Default();
-
+    
     if (nRanks > 1)
     {
         // Use MPI/O to access the file
@@ -897,8 +950,6 @@ void FieldIOHdf5::v_Import(const std::string &infilename,
         // Use collective IO
         readPL = H5::PList::DatasetXfer();
         readPL->SetDxMpioCollective();
-        readPLInd = H5::PList::DatasetXfer();
-        readPLInd->SetDxMpioIndependent();
     }
 
     DataSourceSharedPtr dataSource = H5DataSource::create(
@@ -927,143 +978,333 @@ void FieldIOHdf5::v_Import(const std::string &infilename,
              "Unable to determine Nektar++ HDF5 file version");
     root->GetAttribute("FORMAT_VERSION", formatVersion);
 
+    nRead += sizeof(formatVersion);
+    
     ASSERTL0(formatVersion <= FORMAT_VERSION,
              "File format if " + infilename + " is higher than supported in "
              "this version of Nektar++");
 
-    // Open the datasets
-    H5::DataSetSharedPtr decomps_dset = root->OpenDataSet("DECOMPOSITION");
-    ASSERTL1(decomps_dset, prfx.str() + "cannot open DECOMPOSITION dataset.");
-
-    H5::DataSetSharedPtr ids_dset = root->OpenDataSet("ELEMENTIDS");
-    ASSERTL1(ids_dset, prfx.str() + "cannot open ELEMENTIDS dataset.");
-
-    H5::DataSetSharedPtr data_dset = root->OpenDataSet("DATA");
-    ASSERTL1(data_dset, prfx.str() + "cannot open DATA dataset.");
-
-    // Open the dataset file spaces
-    H5::DataSpaceSharedPtr decomps_fspace = decomps_dset->GetSpace();
-    ASSERTL1(decomps_fspace,
-             prfx.str() + "cannot open DECOMPOSITION filespace.");
-
-    H5::DataSpaceSharedPtr ids_fspace = ids_dset->GetSpace();
-    ASSERTL1(ids_fspace, prfx.str() + "cannot open ELEMENTIDS filespace.");
-
-    H5::DataSpaceSharedPtr data_fspace = data_dset->GetSpace();
-    ASSERTL1(data_fspace, prfx.str() + "cannot open DATA filespace.");
-
-    // Read entire IDS data set to get list of global IDs.
-    std::vector<uint64_t> ids;
-
-    ids_dset->Read(ids, ids_fspace, readPL);
-
-    std::unordered_set<uint64_t> toread;
-    if (ElementIDs != NullInt1DArray)
+    std::unordered_set<uint64_t> selectedIDs;
+    bool selective = false;
+    bool auto_selective = false;
+    if (&ElementIDs != &NullInt1DArray)
     {
-        for (uint64_t i = 0; i < ElementIDs.num_elements(); ++i)
+        if (ElementIDs.num_elements() > 0)
+        {
+            selective = true;
+            for (uint64_t i = 0; i < ElementIDs.num_elements(); ++i)
+            {
+                selectedIDs.insert(ElementIDs[i]);
+            }
+        }
+        else
         {
-            toread.insert(ElementIDs[i]);
+            auto_selective = true;
         }
     }
 
+    // Open and read in the entire the decomps dataset.
+    H5::DataSetSharedPtr decomps_dset = root->OpenDataSet("DECOMPOSITION");
+    ASSERTL1(decomps_dset, "Cannot open DECOMPOSITION dataset.");
+    H5::DataSpaceSharedPtr decomps_fspace = decomps_dset->GetSpace();
+    ASSERTL1(decomps_fspace, "Cannot open DECOMPOSITION filespace.");
+    
     std::vector<uint64_t> decomps;
     decomps_dset->Read(decomps, decomps_fspace, readPL);
+    nRead += decomps.size()*sizeof(uint64_t);
+    
+    uint64_t nDecomps = decomps.size() / MAX_DCMPS;
+    uint64_t i_autosel_min = 0, i_autosel_lim = 0;
+    if (auto_selective)
+    {
+        uint64_t nDecompsPerRank = nDecomps / nRanks;
+        i_autosel_min = m_comm->GetRank()*nDecompsPerRank;
+        i_autosel_lim = i_autosel_min + nDecompsPerRank; 
+    }
 
-    size_t nDecomps = decomps.size() / MAX_DCMPS;
-    size_t cnt = 0, cnt2 = 0;
 
-    // Mapping from each decomposition to offsets in the data array.
-    vector<OffsetHelper> decompsToOffsets (nDecomps);
+    // Open the element ids dataset.
+    H5::DataSetSharedPtr ids_dset = root->OpenDataSet("ELEMENTIDS");
+    ASSERTL1(ids_dset, "Cannot open ELEMENTIDS dataset.");
+    H5::DataSpaceSharedPtr ids_fspace = ids_dset->GetSpace();
+    ASSERTL1(ids_fspace, "Cannot open ELEMENTIDS filespace.");
+
 
+    // Mapping from each decomposition to offsets in the data array.
+    vector<OffsetHelper> decompsToOffsets(nDecomps);
+    
     // Mapping from each group's hash to a vector of element IDs. Note this has
     // to be unsigned int, since that's what we use in FieldDefinitions.
-    map<uint64_t, vector<unsigned int> > groupsToElmts;
-
+    map<uint64_t, vector<unsigned int> > groupsToElems;
+    
     // Mapping from each group's hash to each of the decompositions.
     map<uint64_t, set<uint64_t> > groupsToDecomps;
-
-    // True if we are pulling element IDs from ElementIDs.
-    bool selective = toread.size() > 0;
-
+    
     // Counters for data offsets
     OffsetHelper running;
 
-    for (size_t i = 0; i < nDecomps; ++i, cnt += MAX_DCMPS)
-    {
-        uint64_t nElmt     = decomps[cnt + ELEM_DCMP_IDX];
-        uint64_t groupHash = decomps[cnt + HASH_DCMP_IDX];
+    H5S_seloper_t h5_select = H5S_SELECT_SET;
 
-        vector<uint64_t> tmp;
+    uint64_t nAutoSelectedElems = 0;
+    uint64_t cnt = 0;
+    for (uint64_t i = 0; i < nDecomps; ++i, cnt += MAX_DCMPS)
+    {
+        uint64_t nElems = decomps[cnt + ELEM_DCMP_IDX];
+	uint64_t groupHash = decomps[cnt + HASH_DCMP_IDX];
 
-        if (selective)
+	if (auto_selective)
         {
-            for (size_t j = 0; j < nElmt; ++j)
+            // auto_selective means the number of ranks is equal to
+            // that used to write the checkpoint data.
+            if (i < i_autosel_min)
             {
-                uint64_t elmtId = ids[cnt2 + j];
-                if (toread.find(elmtId) != toread.end())
-                {
-                    tmp.push_back(elmtId);
-                }
+                // Have not yet reached this rank's
+                // decomposition data.
+                nElems = 0;
+            }
+            else if (i >= i_autosel_lim)
+            {
+                // Have reached the end of this rank's
+                // decomposition data.
+                break;
             }
-        }
-        else
-        {
-            tmp.insert(
-                tmp.begin(), ids.begin() + cnt2, ids.begin() + cnt2 + nElmt);
         }
 
-        vector<unsigned int> tmp2(nElmt);
-        for (size_t j = 0; j < nElmt; ++j)
+        if (nElems > 0)
         {
-            tmp2[j] = ids[cnt2+j];
-        }
-
-        cnt2 += nElmt;
+	    if (auto_selective)
+	    {
+	        nAutoSelectedElems += nElems;
+	        ids_fspace->SelectRange(h5_select, running.ids, nElems);
+                h5_select = H5S_SELECT_OR;
+
+                groupsToDecomps[groupHash].insert(i);
+                decompsToOffsets[i] = running;
+	    }
+            else
+	    {
+	        vector<unsigned int> ids(nElems);
+		ids_fspace->SelectRange(running.ids, nElems);
+                ids_dset->Read(ids, ids_fspace, readPL);
+                nRead += ids.size()*sizeof(unsigned int);
+
+                bool found = !selective;
+                if (selective)
+                {
+		    for (auto &it : ids)
+                    {
+		        if (selectedIDs.find(it) != selectedIDs.end())
+                        {
+			    found = true; 
+			    break;
+                        }
+		    }    
+                }
 
-        if (tmp.size() > 0)
-        {
-            groupsToDecomps[groupHash].insert(i);
-        }
+                if (found)
+                {  
+                    groupsToDecomps[groupHash].insert(i);
+                    groupsToElems[i] = ids;
+                    decompsToOffsets[i] = running;
 
-        groupsToElmts[i] = tmp2;
-        decompsToOffsets[i] = running;
+                    if (selective && selectedIDs.size() == 0)
+                    {
+                        break;
+                    }
+                }
+	    } // end of <if (auto_selective)> else clause
+        } // end of <if (nElems > 0)> clause
 
+	running.ids   += decomps[cnt + ELEM_DCMP_IDX];
         running.data  += decomps[cnt + VAL_DCMP_IDX];
         running.order += decomps[cnt + ORDER_DCMP_IDX];
         running.homy  += decomps[cnt + HOMY_DCMP_IDX];
         running.homz  += decomps[cnt + HOMZ_DCMP_IDX];
         running.homs  += decomps[cnt + HOMS_DCMP_IDX];
-    }
+	
+    } // end of <for (uint64_t i = 0; i < nDecomps; ++i, cnt += MAX_DCMPS)> loop
+
+    if (auto_selective)
+    {
+        vector<unsigned int> ids(nAutoSelectedElems);
+        ids_dset->Read(ids, ids_fspace, readPL);
 
+	uint64_t elemOffset = 0;
+	for (auto &gIt : groupsToDecomps)
+        {
+            for (auto &di : gIt.second)
+            {
+                uint64_t nElems = decomps[di*MAX_DCMPS + ELEM_DCMP_IDX];
+                vector<unsigned int> sids(&ids[elemOffset],&ids[elemOffset+nElems]);
+                groupsToElems[di] = sids;
+		elemOffset += nElems;
+            }  
+        }
+    }
+    
+ 
+    H5::DataSetSharedPtr data_dset;
+    H5::DataSpaceSharedPtr data_fspace;
+    data_dset = root->OpenDataSet("DATA");
+    ASSERTL1(data_dset, prfx.str() + "cannot open DATA dataset.");
+    data_fspace = data_dset->GetSpace();
+    ASSERTL1(data_fspace, prfx.str() + "cannot open DATA filespace.");
+    
+
+    // chain together all the necessary HDF5 range selections
+    uint64_t nDataItems = 0;
+    uint64_t nRealDataItems = 0;
+    h5_select = H5S_SELECT_SET;
+    map<uint64_t, uint64_t> dataCountMap;
+    
     for (auto &gIt : groupsToDecomps)
     {
-        // Select region from dataset for this decomposition.
-        for (auto &sIt : gIt.second)
+        for (auto &di : gIt.second)
         {
-            std::stringstream fieldNameStream;
-            fieldNameStream << gIt.first;
+            OffsetHelper offset = decompsToOffsets[di];
+
+            nDataItems = decomps[di*MAX_DCMPS + VAL_DCMP_IDX];
+            data_fspace->SelectRange(h5_select, offset.data, nDataItems);
+            dataCountMap[offset.data] = nDataItems;
+            nRealDataItems += nDataItems;
+            h5_select = H5S_SELECT_OR;            
+        }  
+    }
+
+    // now read the actual HDF5 data
+    std::vector<NekDouble> decompRealData;
+    if (nRealDataItems > 0)
+    {
+        data_dset->Read(decompRealData, data_fspace, readPL);
+        ASSERTL0(decompRealData.size() == nRealDataItems,
+            prfx.str() + "unexpected number of data items.");
+        nRead += nRealDataItems;
+    }
+
 
+    // calculate the adjusted offsets for the recently-read data
+    map<uint64_t, uint64_t> dataOffsetMap;
+    uint64_t offVal;
+
+    offVal = 0;
+    for (auto &offIt : dataCountMap)
+    {
+        dataOffsetMap[offIt.first] = offVal;
+        offVal += offIt.second;
+    }
+
+    // transfer the recently read data to the fielddata data structures
+    // and read in the rest of the field definition data
+    for (auto &gIt : groupsToDecomps)
+    {
+        std::stringstream fieldNameStream;
+        fieldNameStream << gIt.first;
+
+        for (auto &di : gIt.second)
+        {
             FieldDefinitionsSharedPtr fielddef =
                 MemoryManager<FieldDefinitions>::AllocateSharedPtr();
-            ImportFieldDef(readPLInd, root, decomps, sIt, decompsToOffsets[sIt],
-                           fieldNameStream.str(), fielddef);
+            nRead += ImportFieldDef(root, fieldNameStream.str(), fielddef);
 
-            fielddef->m_elementIDs = groupsToElmts[sIt];
+            OffsetHelper offset = decompsToOffsets[di];
+
+            if (fielddef->m_numHomogeneousDir >= 1)
+            {
+                nRead += ImportFieldDataInd(root, "HOMOGENEOUSZIDS", "homogeneousZIDs",
+                             decomps[di*MAX_DCMPS + HOMZ_DCMP_IDX], offset.homz,
+                             fielddef->m_homogeneousZIDs);
+            }
+
+            if (fielddef->m_numHomogeneousDir >= 2)
+            {
+                nRead += ImportFieldDataInd(root, "HOMOGENEOUSYIDS", "homogeneousYIDs",
+                             decomps[di*MAX_DCMPS + HOMY_DCMP_IDX], offset.homy,
+                             fielddef->m_homogeneousYIDs);
+            }
+
+            if (fielddef->m_homoStrips)
+            {
+                nRead += ImportFieldDataInd(root, "HOMOGENEOUSSIDS", "homogeneousSIDs",
+                             decomps[di*MAX_DCMPS + HOMS_DCMP_IDX], offset.homs,
+                             fielddef->m_homogeneousSIDs);
+            }
+
+            if (!fielddef->m_uniOrder)
+            {
+                nRead += ImportFieldDataInd(root, "POLYORDERS", "numModes",
+                             decomps[di*MAX_DCMPS + ORDER_DCMP_IDX], offset.order,
+                             fielddef->m_numModes);
+            }
+            
+            fielddef->m_elementIDs = groupsToElems[di];
             fielddefs.push_back(fielddef);
 
-            if (fielddata != NullVectorNekDoubleVector)
+            if (fielddata != NullVectorNekDoubleVector && nRealDataItems > 0)
             {
-                std::vector<NekDouble> decompFieldData;
-                ImportFieldData(
-                    readPLInd, data_dset, data_fspace,
-                    decompsToOffsets[sIt].data, decomps, sIt, fielddef,
-                    decompFieldData);
-                fielddata.push_back(decompFieldData);
+                vector<NekDouble>::const_iterator dataSelStart = decompRealData.begin() + dataOffsetMap[offset.data]; 
+                vector<NekDouble>::const_iterator dataSelEnd = dataSelStart + dataCountMap[offset.data];
+                vector<NekDouble> dataSel(dataSelStart, dataSelEnd);
+
+                uint64_t datasize = CheckFieldDefinition(fielddef);
+                ASSERTL0(dataSel.size() == datasize * fielddef->m_fields.size(),
+                    prfx.str() + "input data is not the same length as header information.");
+                
+                fielddata.push_back(dataSel);
             }
-        }
-    }
+            
+        }  // end of <for (auto &di : gIt.second)> loop
+        
+    } // end of <for (auto &gIt : groupsToDecomps)> loop
 
     m_comm->Block();
+
+    return nRead;
+    
+}
+
+/**
+ * @brief Import a number of data items from a specified position within a hdf5 dataset.
+ *
+ * The data are imported independently.
+ *
+ * @param root        root group containing field definitions
+ * @param dsetName    the name of the hdf5 dataset
+ * @param dataTag     data item descriptive tag
+ * @param nDataItems  the number of items to be imported from the dataset
+ * @param offset      the position within dataset from where reading should commence
+ * @param data        the data vector for storing data read
+ * @return The number of bytes read.
+ */
+template <class T>
+uint64_t FieldIOHdf5::ImportFieldDataInd(
+    H5::GroupSharedPtr root,
+    std::string dsetName,
+    std::string dataTag,
+    uint64_t nDataItems,
+    uint64_t offset,
+    std::vector<T> &data)
+{
+    if (nDataItems == 0) return 0;
+
+    std::stringstream prfx;
+    prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldDataInd(): ";
+
+    H5::PListSharedPtr readPL = H5::PList::Default();
+    readPL = H5::PList::DatasetXfer();
+    readPL->SetDxMpioIndependent();
+
+    H5::DataSetSharedPtr dset = root->OpenDataSet(dsetName);
+    ASSERTL1(dset, prfx.str() + "cannot open " + dsetName + " dataset.");
+    H5::DataSpaceSharedPtr fspace = dset->GetSpace();
+    ASSERTL1(fspace, prfx.str() + "cannot open " + dsetName + " filespace.");
+
+    fspace->SelectRange(offset, nDataItems);
+
+    dset->Read(data, fspace, readPL);
+    ASSERTL0(data.size() == nDataItems,
+        prfx.str() + "unexpected number of " + dataTag + " items.");
+
+    return nDataItems*sizeof(T);
 }
 
 /**
@@ -1073,19 +1314,17 @@ void FieldIOHdf5::v_Import(const std::string &infilename,
  * @param root         Root group containing field definitions.
  * @param group        Group name to process.
  * @param def          On output contains field definitions.
+ * @return The number of bytes read.
  */
-void FieldIOHdf5::ImportFieldDef(
-    H5::PListSharedPtr        readPL,
+uint64_t FieldIOHdf5::ImportFieldDef(
     H5::GroupSharedPtr        root,
-    std::vector<uint64_t>    &decomps,
-    uint64_t                  decomp,
-    OffsetHelper              offset,
     std::string               group,
     FieldDefinitionsSharedPtr def)
 {
     std::stringstream prfx;
-    prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldDefsHdf5(): ";
-
+    prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldDef(): ";
+    
+    uint64_t nRead = 0;
     H5::GroupSharedPtr field = root->OpenGroup(group);
     ASSERTL1(field, prfx.str() + "cannot open field group, " + group + '.');
 
@@ -1099,15 +1338,20 @@ void FieldIOHdf5::ImportFieldDef(
         if (attrName == "FIELDS")
         {
             field->GetAttribute(attrName, def->m_fields);
+            for (int i = 0; i < def->m_fields.size(); i++)
+            {  
+                nRead += def->m_fields[i].size();
+            }
         }
         else if (attrName == "SHAPE")
         {
             std::string shapeString;
             field->GetAttribute(attrName, shapeString);
+            nRead += shapeString.size();
 
             // check to see if homogeneous expansion and if so
             // strip down the shapeString definition
-            size_t loc;
+            int loc;
             //---> this finds the first location of 'n'!
             if (shapeString.find("Strips") != string::npos)
             {
@@ -1147,10 +1391,14 @@ void FieldIOHdf5::ImportFieldDef(
         else if (attrName == "BASIS")
         {
             field->GetAttribute(attrName, def->m_basis);
+            nRead += def->m_basis.size()*sizeof(LibUtilities::BasisType);
+            
             // check the basis is in range
-            for (auto &bIt : def->m_basis)
+            std::vector<BasisType>::const_iterator bIt  = def->m_basis.begin();
+            std::vector<BasisType>::const_iterator bEnd = def->m_basis.end();
+            for (; bIt != bEnd; ++bIt)
             {
-                BasisType bt = bIt;
+                BasisType bt = *bIt;
                 ASSERTL0(bt >= 0 && bt < SIZE_BasisType,
                          prfx.str() +
                          "unable to correctly parse the basis types.");
@@ -1159,11 +1407,13 @@ void FieldIOHdf5::ImportFieldDef(
         else if (attrName == "HOMOGENEOUSLENGTHS")
         {
             field->GetAttribute(attrName, def->m_homogeneousLengths);
+            nRead += def->m_homogeneousLengths.size()*sizeof(NekDouble);
         }
         else if (attrName == "NUMMODESPERDIR")
         {
             std::string numModesPerDir;
             field->GetAttribute(attrName, numModesPerDir);
+            nRead += numModesPerDir.size();
 
             if (strstr(numModesPerDir.c_str(), "UNIORDER:"))
             {
@@ -1179,6 +1429,7 @@ void FieldIOHdf5::ImportFieldDef(
         {
             std::string pointsString;
             field->GetAttribute(attrName, pointsString);
+            nRead += pointsString.size();
             def->m_pointsDef = true;
 
             std::vector<std::string> pointsStrings;
@@ -1215,6 +1466,7 @@ void FieldIOHdf5::ImportFieldDef(
         {
             std::string numPointsPerDir;
             field->GetAttribute(attrName, numPointsPerDir);
+            nRead += numPointsPerDir.size();
             def->m_numPointsDef = true;
 
             bool valid = ParseUtils::GenerateVector(
@@ -1231,79 +1483,9 @@ void FieldIOHdf5::ImportFieldDef(
         }
     }
 
-    if (def->m_numHomogeneousDir >= 1)
-    {
-        H5::DataSetSharedPtr homz_dset = root->OpenDataSet("HOMOGENEOUSZIDS");
-        H5::DataSpaceSharedPtr homz_fspace = homz_dset->GetSpace();
-        uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMZ_DCMP_IDX];
-        homz_fspace->SelectRange(offset.homz, dsize);
-        homz_dset->Read(def->m_homogeneousZIDs, homz_fspace, readPL);
-    }
-
-    if (def->m_numHomogeneousDir >= 2)
-    {
-        H5::DataSetSharedPtr homy_dset = root->OpenDataSet("HOMOGENEOUSYIDS");
-        H5::DataSpaceSharedPtr homy_fspace = homy_dset->GetSpace();
-        uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMY_DCMP_IDX];
-        homy_fspace->SelectRange(offset.homy, dsize);
-        homy_dset->Read(def->m_homogeneousYIDs, homy_fspace, readPL);
-    }
-
-    if (def->m_homoStrips)
-    {
-        H5::DataSetSharedPtr homs_dset = root->OpenDataSet("HOMOGENEOUSSIDS");
-        H5::DataSpaceSharedPtr homs_fspace = homs_dset->GetSpace();
-        uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMS_DCMP_IDX];
-        homs_fspace->SelectRange(offset.homs, dsize);
-        homs_dset->Read(def->m_homogeneousSIDs, homs_fspace, readPL);
-    }
-
-    if (!def->m_uniOrder)
-    {
-        H5::DataSetSharedPtr order_dset = root->OpenDataSet("POLYORDERS");
-        H5::DataSpaceSharedPtr order_fspace = order_dset->GetSpace();
-        uint64_t dsize = decomps[decomp * MAX_DCMPS + ORDER_DCMP_IDX];
-        order_fspace->SelectRange(offset.order, dsize);
-        order_dset->Read(def->m_numModes, order_fspace, readPL);
-    }
+    return nRead;
 }
 
-/**
- * @brief Import the field data from the HDF5 document.
- *
- * @param readPL       Reading parameter list.
- * @param data_dset    Pointer to the `DATA` dataset.
- * @param data_fspace  Pointer to the `DATA` data space.
- * @param data_i       Index in the `DATA` dataset to start reading from.
- * @param decomps      Information from the `DECOMPOSITION` dataset.
- * @param decomp       Index of the decomposition.
- * @param fielddef     Field definitions for this file
- * @param fielddata    On return contains resulting field data.
- */
-void FieldIOHdf5::ImportFieldData(
-    H5::PListSharedPtr               readPL,
-    H5::DataSetSharedPtr             data_dset,
-    H5::DataSpaceSharedPtr           data_fspace,
-    uint64_t                         data_i,
-    std::vector<uint64_t>           &decomps,
-    uint64_t                         decomp,
-    const FieldDefinitionsSharedPtr  fielddef,
-    std::vector<NekDouble>          &fielddata)
-{
-    std::stringstream prfx;
-    prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldData(): ";
-
-    uint64_t nElemVals  = decomps[decomp * MAX_DCMPS + VAL_DCMP_IDX];
-    uint64_t nFieldVals = nElemVals;
-
-    data_fspace->SelectRange(data_i, nFieldVals);
-    data_dset->Read(fielddata, data_fspace, readPL);
-    int datasize = CheckFieldDefinition(fielddef);
-    ASSERTL0(
-        fielddata.size() == datasize * fielddef->m_fields.size(),
-        prfx.str() +
-        "input data is not the same length as header information.");
-}
 
 /**
  * @brief Import field metadata from @p filename and return the data source
@@ -1319,6 +1501,7 @@ DataSourceSharedPtr FieldIOHdf5::v_ImportFieldMetaData(
     H5::PListSharedPtr parallelProps = H5::PList::Default();
     DataSourceSharedPtr ans = H5DataSource::create(filename, parallelProps);
     ImportHDF5FieldMetaData(ans, fieldmetadatamap);
+
     return ans;
 }
 
diff --git a/library/LibUtilities/BasicUtils/FieldIOHdf5.h b/library/LibUtilities/BasicUtils/FieldIOHdf5.h
index 9882a220a9d7fcc804a8244d66cfa6c18212df99..39c0ae34f95c3dfd22d2375b4864c47c1ad6d079 100644
--- a/library/LibUtilities/BasicUtils/FieldIOHdf5.h
+++ b/library/LibUtilities/BasicUtils/FieldIOHdf5.h
@@ -227,55 +227,55 @@ public:
 
 private:
     struct OffsetHelper {
-        OffsetHelper() : data(0), order(0), homy(0), homz(0), homs(0) {}
+        OffsetHelper() : ids(0), data(0), order(0), homy(0), homz(0), homs(0) {}
         OffsetHelper(const OffsetHelper &in) :
-            data(in.data), order(in.order), homy(in.homy), homz(in.homz),
-            homs(in.homs)
+            ids(in.ids), data(in.data), order(in.order),
+                homy(in.homy), homz(in.homz), homs(in.homs)
         {
         }
 
-        uint64_t data, order, homy, homz, homs;
+      uint64_t ids, data, order, homy, homz, homs;
     };
-
-    LIB_UTILITIES_EXPORT virtual void v_Write(
+    
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Write(
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
         const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
         const bool backup = false);
 
-    LIB_UTILITIES_EXPORT virtual void v_Import(
-        const std::string &infilename,
+    template <class T>
+    LIB_UTILITIES_EXPORT uint64_t WriteFieldDataInd(std::size_t nFields,
+        H5::DataSpaceSharedPtr &fspace, H5::DataSetSharedPtr &dset,
+        uint64_t data_i, std::vector<std::vector<T> > &data);
+    
+    template <class T>
+    LIB_UTILITIES_EXPORT uint64_t WriteFieldData(std::size_t nMinFields, std::size_t nFields,
+        H5::DataSpaceSharedPtr &fspace, H5::DataSetSharedPtr &dset,
+        uint64_t data_i, std::vector<std::vector<T> > &data);
+
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Import(const std::string &infilename,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-        std::vector<std::vector<NekDouble> > &fielddata =
-                                                      NullVectorNekDoubleVector,
+        std::vector<std::vector<NekDouble> > &fielddata = NullVectorNekDoubleVector,
         FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
         const Array<OneD, int> &ElementIDs = NullInt1DArray);
 
+    template <class T>
+    LIB_UTILITIES_EXPORT uint64_t ImportFieldDataInd(H5::GroupSharedPtr root,
+        std::string dsetName, std::string dataTag, uint64_t nDataItems,
+        uint64_t offset, std::vector<T> &data);
+
+    LIB_UTILITIES_EXPORT uint64_t ImportFieldDef(H5::GroupSharedPtr root,
+        std::string               group,
+        FieldDefinitionsSharedPtr def);
+    
     LIB_UTILITIES_EXPORT virtual DataSourceSharedPtr v_ImportFieldMetaData(
         const std::string &filename, FieldMetaDataMap &fieldmetadatamap);
 
     LIB_UTILITIES_EXPORT void ImportHDF5FieldMetaData(
         DataSourceSharedPtr dataSource, FieldMetaDataMap &fieldmetadatamap);
 
-    LIB_UTILITIES_EXPORT void ImportFieldDef(
-        H5::PListSharedPtr        readPL,
-        H5::GroupSharedPtr        root,
-        std::vector<uint64_t>    &decomps,
-        uint64_t                  decomp,
-        OffsetHelper              offset,
-        std::string               group,
-        FieldDefinitionsSharedPtr def);
-
-    LIB_UTILITIES_EXPORT void ImportFieldData(
-        H5::PListSharedPtr               readPL,
-        H5::DataSetSharedPtr             data_dset,
-        H5::DataSpaceSharedPtr           data_fspace,
-        uint64_t                         data_i,
-        std::vector<uint64_t>           &decomps,
-        uint64_t                         decomp,
-        const FieldDefinitionsSharedPtr  fielddef,
-        std::vector<NekDouble>          &fielddata);
+    
 };
 }
 }
diff --git a/library/LibUtilities/BasicUtils/FieldIOSIONlib.cpp b/library/LibUtilities/BasicUtils/FieldIOSIONlib.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..99ff709b4d268b06c7063cdf35e3d2859248dd32
--- /dev/null
+++ b/library/LibUtilities/BasicUtils/FieldIOSIONlib.cpp
@@ -0,0 +1,772 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: FieldIOSIONlib.cpp
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: I/O routines relating to Fields into SIONlib files.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "sion.h"
+
+#include <LibUtilities/BasicUtils/SIONlib.h>
+#include <LibUtilities/BasicUtils/FieldIOSIONlib.h>
+#include <LibUtilities/BasicUtils/ParseUtils.h>
+#include <boost/unordered_set.hpp>
+#include <LibUtilities/Communication/CommMpi.h>
+
+
+namespace berrc = boost::system::errc;
+
+namespace Nektar
+{
+namespace LibUtilities
+{ 
+
+std::string FieldIOSIONlib::className =
+    GetFieldIOFactory().RegisterCreatorFunction(
+        "SIONlib", FieldIOSIONlib::create, "SIONlib-based output of field data.");
+
+  
+const unsigned int FieldIOSIONlib::FORMAT_VERSION = 1;
+  
+const std::string FieldIOSIONlib::ATTRNAME_FIELDS = "FIELDS";
+const std::string FieldIOSIONlib::ATTRNAME_BASIS = "BASIS";
+const std::string FieldIOSIONlib::ATTRNAME_SHAPE = "SHAPE";
+const std::string FieldIOSIONlib::ATTRNAME_HOMOLENS = "HOMOGENEOUSLENGTHS";
+const std::string FieldIOSIONlib::ATTRNAME_NUMMODES = "NUMMODESPERDIR";
+const std::string FieldIOSIONlib::ATTRNAME_POINTSTYPE = "POINTSTYPE";
+const std::string FieldIOSIONlib::ATTRNAME_NUMPOINTS = "NUMPOINTSPERDIR";
+  
+const std::string FieldIOSIONlib::ATTRVALUE_MIXORDER = "MIXORDER";
+const std::string FieldIOSIONlib::ATTRVALUE_UNIORDER = "UNIORDER";
+
+sion_int32 FieldIOSIONlib::block_size = -1;
+sion_int64 FieldIOSIONlib::chunk_size = -1;
+std::string FieldIOSIONlib::read_mode = "";
+std::string FieldIOSIONlib::write_mode = "";
+
+  
+FieldIOSIONlib::FieldIOSIONlib(LibUtilities::CommSharedPtr pComm,
+                         bool sharedFilesystem)
+    : FieldIO(pComm, sharedFilesystem)
+{
+}
+
+
+void FieldIOSIONlib::v_Init(const LibUtilities::SessionReaderSharedPtr session)
+{
+    block_size = -1;
+    if (session->DefinesSolverInfo("IOBlockSize"))
+    {
+        std::string int_str = session->GetSolverInfo("IOBlockSize");
+        std::istringstream(int_str) >> block_size;
+    }
+
+    NekDouble blocks_per_chunk = 1.0;
+    if (session->DefinesSolverInfo("IOBlocksPerChunk"))
+    {
+        std::string int_str = session->GetSolverInfo("IOBlocksPerChunk");
+        std::istringstream(int_str) >> blocks_per_chunk;
+    }
+    chunk_size = (sion_int64) (block_size*blocks_per_chunk);
+
+    read_mode = "br";
+    if (session->DefinesSolverInfo("IOReadMode"))
+    {
+        read_mode = session->GetSolverInfo("IOReadMode");
+    }
+
+    write_mode = "bw";
+    if (session->DefinesSolverInfo("IOWriteMode"))
+    {
+        write_mode = session->GetSolverInfo("IOWriteMode");
+    }
+}
+
+
+void FieldIOSIONlib::v_InitFromBenchmarker(const LibUtilities::IOSettingsSharedPtr iosettings)
+{
+    LibUtilities::IOSettings::iterator it;
+    
+    block_size = -1;
+    it = iosettings->find("IOBlockSize");
+    if (iosettings->end() != it)
+    {
+        std::istringstream(it->second) >> block_size;
+    }
+
+    NekDouble blocks_per_chunk = 1.0;
+    it = iosettings->find("IOBlocksPerChunk");
+    if (iosettings->end() != it)
+    {
+        std::istringstream(it->second) >> blocks_per_chunk;
+    }
+    chunk_size = (sion_int64) (block_size*blocks_per_chunk);
+
+    read_mode = "br";
+    it = iosettings->find("IOReadMode");
+    if (iosettings->end() != it)
+    {
+        read_mode = it->second;
+    }
+
+    write_mode = "bw";
+    it = iosettings->find("IOWriteMode");
+    if (iosettings->end() != it)
+    {
+        write_mode = it->second;
+    }
+}
+
+  
+/**
+ * Open a SIONlib file for writing by constructing and returning
+ * a pointer to a SIONlib object, and also outputting the format version.
+ *
+ * @param outFilename  name of output file.
+ * @return Pointer to SIONlib object.
+ */
+SIONlib::SIONFile *FieldIOSIONlib::OpenFileForWriting(const std::string &outFilename)
+{
+    CommMpi *commMpiPtr =  (CommMpi*) m_comm.get();
+    SIONlib::SIONFile *fp = new SIONlib::SIONFile(outFilename, write_mode,
+        1, chunk_size, block_size,
+        m_comm->GetRank(), commMpiPtr->GetComm(), commMpiPtr->GetComm());
+
+    ASSERTL0(NULL != fp,
+             "Unable to construct SIONFile object for " + outFilename);
+
+    SIONlib::SIONFile &f = *fp;
+    f.open();
+    
+    ASSERTL0(-1 != f.getReturnCode(),
+             "Unable to open SIONlib file " + outFilename);
+
+    f << FORMAT_VERSION;
+    
+    return fp;
+}
+
+
+  
+void FieldIOSIONlib::WriteRawDataToBinaryVector(std::vector<NekByte> &v, const NekByte* dptr, const size_t dsize)
+{
+    std::vector<NekByte> v2(dsize);
+    std::copy(dptr, dptr+dsize, v2.begin());
+
+    v.reserve(v.size() + v2.size());
+    v.insert(v.end(), v2.begin(), v2.end());
+}
+
+void FieldIOSIONlib::WriteStringToBinaryVector(std::vector<NekByte> &v, const std::string &s)
+{
+    size_t ssize = s.size();
+    WriteRawDataToBinaryVector(v, (const NekByte*) &ssize, sizeof(size_t));
+    WriteRawDataToBinaryVector(v, (const NekByte*) s.c_str(), ssize);  
+}
+
+template<class T>   
+void FieldIOSIONlib::WriteDatumToBinaryVector(std::vector<NekByte> &v, const T d)
+{
+    WriteRawDataToBinaryVector(v, (const NekByte*) &d, sizeof(T));
+}
+
+template<class T>   
+void FieldIOSIONlib::WriteDataToBinaryVector(std::vector<NekByte> &v, const std::vector<T> &d)
+{
+    size_t dsize = d.size();
+    WriteRawDataToBinaryVector(v, (const NekByte*) &dsize, sizeof(size_t));
+    for (auto it = d.begin(); it != d.end(); ++it)
+    {
+        WriteRawDataToBinaryVector(v, (const NekByte*) &(*it), sizeof(T));
+    }
+}
+
+  
+unsigned long FieldIOSIONlib::v_Write(const std::string &outFilePrefix,
+    std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+    std::vector<std::vector<NekDouble> > &fielddata,
+    const FieldMetaDataMap &fieldmetadatamap,
+    const bool backup)
+{
+    // We make a number of assumptions in this code:
+    //   1. All element ids have the same type: unsigned int
+    //   2. All elements within a given field have the same number of values
+    //   3. All element values have the same type, NekDouble
+
+    ASSERTL1(fielddefs.size() == fielddata.size(),
+             "fielddefs and fielddata have incompatible lengths.");
+
+    size_t nFields = fielddefs.size();
+    size_t nMaxFields = nFields;
+
+    int homDim = -1;
+    int varOrder = 0;
+
+    for (size_t i = 0; i < nFields; ++i)
+    {
+        if (!fielddefs[i]->m_uniOrder)
+        {
+            varOrder = 1;
+            break;
+        }
+    }
+
+    m_comm->AllReduce(varOrder, LibUtilities::ReduceMax);
+    m_comm->AllReduce(nMaxFields, LibUtilities::ReduceMax);
+
+    SetUpOutput(outFilePrefix, false, backup);
+
+    // Each MPI process iterates through its fields and outputs field
+    // definitions and field data to the SIONlib file.
+    SIONlib::SIONFile* fp = OpenFileForWriting(outFilePrefix);
+    ASSERTL0(NULL != fp, "Cannot open SIONlib file.");
+    SIONlib::SIONFile &f = *fp; 
+
+    if (f.isModeCollectiveMerge())
+    {
+        f << (size_t) m_comm->GetRank();
+    }
+    f << nFields;
+
+    std::vector<NekByte> attrVector;
+    
+    for (size_t i = 0; i < nFields; ++i)
+    {
+        ASSERTL1(fielddata[i].size() > 0,
+            "fielddata vector must contain at least one value.");
+        ASSERTL1(fielddata[i].size() ==
+            fielddefs[i]->m_fields.size()*CheckFieldDefinition(fielddefs[i]),
+            "fielddata vector has invalid size.");
+
+        FieldDefinitionsSharedPtr def = fielddefs[i];
+
+        size_t nAttrs = (def->m_numHomogeneousDir ? 5 : 4);
+        if (def->m_pointsDef) nAttrs++;
+        if (def->m_numPointsDef) nAttrs++;
+        
+        WriteDatumToBinaryVector(attrVector, nAttrs);
+	
+        // FIELDS
+        ////////////////////////////////////////////////////////////////
+        WriteStringToBinaryVector(attrVector, ATTRNAME_FIELDS);
+        
+        size_t nAttrFields = def->m_fields.size();
+        WriteDatumToBinaryVector(attrVector, nAttrFields);
+
+        if (nAttrFields > 0)
+        {
+            for (size_t j = 0; j < nAttrFields; ++j)
+            { 
+                WriteStringToBinaryVector(attrVector, def->m_fields[j]);
+            }
+        }
+        ////////////////////////////////////////////////////////////////
+                
+        // BASIS
+        ////////////////////////////////////////////////////////////////
+        WriteStringToBinaryVector(attrVector, ATTRNAME_BASIS);
+
+	WriteDataToBinaryVector(attrVector, def->m_basis);
+        ////////////////////////////////////////////////////////////////
+        
+        // SHAPE
+        ////////////////////////////////////////////////////////////////
+        WriteStringToBinaryVector(attrVector, ATTRNAME_SHAPE);
+
+        std::stringstream shapeStringStream;
+        shapeStringStream << ShapeTypeMap[def->m_shapeType];
+
+        if (def->m_numHomogeneousDir > 0)
+        {
+            if (homDim == -1)
+            {
+                homDim = def->m_numHomogeneousDir;
+            }
+
+            ASSERTL1(homDim == def->m_numHomogeneousDir,
+                "HDF5 does not support variable homogeneous directions in "
+                "the same file.");
+
+            shapeStringStream << "-HomogenousExp"
+                << def->m_numHomogeneousDir << "D";
+        }
+
+        if (def->m_homoStrips)
+        {
+            shapeStringStream << "-Strips";
+        }
+        
+        WriteStringToBinaryVector(attrVector, shapeStringStream.str());
+        ////////////////////////////////////////////////////////////////
+        
+        // Determine HOMOGENEOUS attributes
+        if (def->m_numHomogeneousDir)
+        {
+            // HOMOGENEOUSLENGTHS
+            ////////////////////////////////////////////////////////////////
+	    WriteStringToBinaryVector(attrVector, ATTRNAME_HOMOLENS);
+	    WriteDataToBinaryVector(attrVector, def->m_homogeneousLengths);
+            ////////////////////////////////////////////////////////////////
+
+            // homo IDs are streamed from separate arrays
+            ////////////////////////////////////////////////////////////////
+	    WriteDataToBinaryVector(attrVector, def->m_homogeneousYIDs);
+	    WriteDataToBinaryVector(attrVector, def->m_homogeneousZIDs);
+	    WriteDataToBinaryVector(attrVector, def->m_homogeneousSIDs);
+            ////////////////////////////////////////////////////////////////
+        }
+        
+        // NUMMODESPERDIR
+        ////////////////////////////////////////////////////////////////
+        WriteStringToBinaryVector(attrVector, ATTRNAME_NUMMODES);
+
+        std::vector<unsigned int>& numModes = def->m_numModes;
+        size_t nNumModes = def->m_basis.size();
+        ASSERTL1(nNumModes == def->m_numModes.size(),
+             "basis and numModes have incompatible lengths.");
+        
+        if (def->m_uniOrder && !varOrder)
+        {
+	    WriteStringToBinaryVector(attrVector, ATTRVALUE_UNIORDER);
+        }
+        else
+        {
+	    WriteStringToBinaryVector(attrVector, ATTRVALUE_MIXORDER);
+        }
+
+        if (nNumModes > 0)
+        {
+	    WriteDataToBinaryVector(attrVector, def->m_numModes);
+        }
+        ////////////////////////////////////////////////////////////////
+
+        // POINTSTYPE
+        ////////////////////////////////////////////////////////////////
+        if (def->m_pointsDef)
+        {
+	    WriteStringToBinaryVector(attrVector, ATTRNAME_POINTSTYPE);
+
+            stringstream pointsTypeStringStream;
+            std::vector<LibUtilities::PointsType>& points = def->m_points;
+            size_t nPoints = points.size();
+            
+            for (size_t j = 0; j < nPoints; ++j)
+            {
+                if (j > 0)
+                {
+                    pointsTypeStringStream << ",";
+                }
+                pointsTypeStringStream << kPointsTypeStr[points[j]];
+            }
+
+            WriteStringToBinaryVector(attrVector, pointsTypeStringStream.str());
+        }
+        ////////////////////////////////////////////////////////////////
+
+        // NUMPOINTSPERDIR
+        ////////////////////////////////////////////////////////////////
+        if (def->m_numPointsDef)
+        {
+	    WriteStringToBinaryVector(attrVector, ATTRNAME_NUMPOINTS);
+
+            stringstream numPointsStringStream;
+            std::vector<unsigned int>& numPoints = def->m_numPoints;
+            size_t nNumPoints = numPoints.size();
+            
+            for (size_t j = 0; j < nNumPoints; ++j)
+            {
+                if (j > 0)
+                {
+                    numPointsStringStream << ",";
+                }
+                numPointsStringStream << numPoints[j];
+            }
+
+            WriteStringToBinaryVector(attrVector, numPointsStringStream.str());
+        }
+        ////////////////////////////////////////////////////////////////
+
+	// field attributes
+	////////////////////////////////////////////////////////////////
+        f << attrVector;
+	attrVector.clear();
+	////////////////////////////////////////////////////////////////
+
+	// elementIDs
+        ////////////////////////////////////////////////////////////////
+        f << def->m_elementIDs;
+        ////////////////////////////////////////////////////////////////
+
+        // data
+        ////////////////////////////////////////////////////////////////
+        f << fielddata[i];
+        ////////////////////////////////////////////////////////////////
+    }
+    
+    for (size_t i = nFields; i < nMaxFields; ++i)
+    {
+        const std::vector<NekByte> dumAttributes = {0};
+	const std::vector<unsigned int> dumElementIds = {0};
+	const std::vector<NekDouble> dumData = {0.0};
+        f << dumAttributes;
+	f << dumElementIds;
+	f << dumData;
+    }
+    
+    unsigned long nWritten = f.getBytesWritten();
+    f.close();
+    
+    m_comm->Block();
+    return nWritten;
+}
+  
+
+/**
+ * Open a SIONlib file for reading by constructing and returning
+ * a pointer to a SIONlib object, and also read in the format version,
+ * throwing an error if it is greater than FieldIOSIONlib::FORMAT_VERSION.
+ *
+ * @param inFilename  name of input file.
+ * @return Pointer to SIONlib object.
+ */
+SIONlib::SIONFile *FieldIOSIONlib::OpenFileForReading(const std::string &inFilename)
+{
+    CommMpi *commMpiPtr =  (CommMpi*) m_comm.get();
+    SIONlib::SIONFile *fp = new SIONlib::SIONFile(inFilename, read_mode,
+        1, chunk_size, block_size,
+        m_comm->GetRank(), commMpiPtr->GetComm(), commMpiPtr->GetComm());
+
+    ASSERTL0(NULL != fp,
+             "Unable to construct SIONFile object for " + inFilename);
+
+    SIONlib::SIONFile &f = *fp;
+    f.open();
+    
+    ASSERTL0(-1 != f.getReturnCode(),
+             "Unable to open SIONlib file " + inFilename);
+
+    unsigned int formatVersion = 0;
+    f >> formatVersion;
+    
+    std::stringstream version;
+    version << formatVersion;
+
+    ASSERTL0(formatVersion <= FORMAT_VERSION,
+             "File " + inFilename + " is at version " + version.str() +
+             ", which is higher than supported in this version of Nektar++.");
+    
+    return fp;
+}
+
+
+unsigned long FieldIOSIONlib::v_Import(const std::string &infilename,
+    std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+    std::vector<std::vector<NekDouble> > &fielddata,
+    FieldMetaDataMap &fieldinfomap,
+    const Array<OneD, int> &ElementIDs)  
+{
+    SIONlib::SIONFile* fp = OpenFileForReading(infilename);
+    ASSERTL0(NULL != fp, "Cannot open SIONlib file.");
+    SIONlib::SIONFile &f = *fp;
+    
+    if (f.isModeCollectiveMerge())
+    {
+        unsigned long bytesWritten = f.getBytesWritten();
+        unsigned long bytesRead = 0;
+    
+        if (bytesWritten > 0)
+        {
+            // this process is one of the collectors that originally wrote the checkpoint file
+            do
+	    {
+	        size_t rk = 0;
+	        f >> rk;
+	        bytesRead += sizeof(rk);
+	    
+	        if (rk != m_comm->GetRank())
+	        {
+	            // the data stored hereafter originated from process rk
+		  
+		    // TODO: send the number of fields to rk
+	            // TODO: for each field do the following,
+	            //         pack the data from rk into three buffers
+	            //             (i) field attributes, vector<NekByte>
+	            //            (ii) element ids, vector<unsigned int>
+	            //           (iii) element data, vector<double>
+	            //         send each buffer to process rk
+
+	            // TODO: update bytesRead
+	        }
+	        else
+	        {
+		    bytesRead += DirectImport(f, fielddefs,fielddata,fieldinfomap,ElementIDs);
+	        }
+	    }
+	    while (bytesWritten > bytesRead);
+	}
+	else
+	{
+	    // this process is not a collector, i.e., its data was written to the area of the
+	    // checkpoint file reserved for a collector
+	  
+	    // TODO: wait to receive the field count from an unknown collector
+	    // TODO: for each field do the following,
+	    //           (i) receive the field attributes
+	    //          (ii) receive the element ids
+	    //         (iii) receive the element data
+	}
+    }
+    else
+    {
+        return DirectImport(f, fielddefs,fielddata,fieldinfomap,ElementIDs);
+    }
+}
+
+
+unsigned long FieldIOSIONlib::DirectImport(SIONlib::SIONFile &f,
+    std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+    std::vector<std::vector<NekDouble> > &fielddata,
+    FieldMetaDataMap &fieldinfomap,
+    const Array<OneD, int> &ElementIDs)  
+{
+    size_t nFields = 0;
+    f >> nFields;
+    
+    for (size_t i = 0; i < nFields; ++i)
+    {
+        FieldDefinitionsSharedPtr def =
+            MemoryManager<FieldDefinitions>::AllocateSharedPtr();
+        
+        std::string attrName, attrValue;
+        size_t nAttrs, nAttrFields;
+
+        f >> nAttrs;
+        
+        for (size_t j = 0; j < nAttrs; ++j)
+        {         
+            f >> attrName;
+
+            if (attrName == ATTRNAME_FIELDS)
+            {
+                f >> nAttrFields;
+                
+                def->m_fields.resize(nAttrFields);
+
+                for (size_t k = 0; k < nAttrFields; ++k)
+                {
+                    f >> def->m_fields[k];
+                }
+            }
+            else if (attrName == ATTRNAME_BASIS)
+            {
+                f >> def->m_basis;
+
+                // check the basis is in range
+                std::vector<BasisType>::const_iterator bIt  = def->m_basis.begin();
+                std::vector<BasisType>::const_iterator bEnd = def->m_basis.end();
+                for (; bIt != bEnd; ++bIt)
+                {
+                    BasisType bt = *bIt;
+                    ASSERTL0(bt >= 0 && bt < SIZE_BasisType,
+                        "Unable to correctly parse the basis types.");
+                }
+            }
+            else if (attrName == ATTRNAME_SHAPE)
+            {
+                f >> attrValue;        
+
+                // check to see if homogeneous expansion and if so
+                // strip down the shapeString definition
+                size_t loc;
+                //---> this finds the first location of 'n'!
+                if (attrValue.find("Strips") != string::npos)
+                {
+                    def->m_homoStrips = true;
+                }
+
+                if ((loc = attrValue.find_first_of("-")) != string::npos)
+                {
+                    if (attrValue.find("Exp1D") != string::npos)
+                    {
+                        def->m_numHomogeneousDir = 1;
+                    }
+                    else // HomogeneousExp1D
+                    {
+                        def->m_numHomogeneousDir = 2;
+                    }
+
+                    attrValue.erase(loc, attrValue.length());
+                }
+
+                // get the geometrical shape
+                bool valid = false;
+                for (unsigned int k = 0; k < SIZE_ShapeType; k++)
+                {
+                    if (ShapeTypeMap[k] == attrValue)
+                    {
+                        def->m_shapeType = (ShapeType) k;
+                        valid = true;
+                        break;
+                    }
+                }
+
+                ASSERTL0(valid,
+                    std::string("Unable to correctly parse the shape type: ")
+                        .append(attrValue)
+                        .c_str());
+            }
+            else if (attrName == ATTRNAME_HOMOLENS)
+            {
+                def->m_numHomogeneousDir = true;
+          
+                f >> def->m_homogeneousLengths;
+                
+                f >> def->m_homogeneousYIDs;
+                f >> def->m_homogeneousZIDs;
+                f >> def->m_homogeneousSIDs;
+            }
+            else if (attrName == ATTRNAME_NUMMODES)
+            {
+                f >> attrValue;
+
+                if (attrValue == ATTRVALUE_UNIORDER)
+                {
+                    def->m_uniOrder = true;
+                    f >> def->m_numModes;
+                }
+                else if (attrValue == ATTRVALUE_MIXORDER)
+                {
+                    f >> def->m_numModes;
+                }
+                else
+                {
+                    std::string errstr("Unknown " + attrName + " value: ");
+                    errstr += attrValue;
+                    ASSERTL1(false, errstr.c_str());
+                }
+            }
+            else if (attrName == ATTRNAME_POINTSTYPE)
+            {
+                def->m_pointsDef = true;
+            
+                f >> attrValue;
+            
+                std::vector<std::string> pointsStrings;
+                bool valid = ParseUtils::GenerateVector(
+                    attrValue.c_str(), pointsStrings);
+
+                ASSERTL0(valid,
+                    "Unable to correctly parse the points types.");
+                
+                for (std::vector<std::string>::size_type k = 0;
+                     k < pointsStrings.size();
+                     ++k)
+                {
+                    valid = false;
+                    for (unsigned int l = 0; l < SIZE_PointsType; ++l)
+                    {
+                        if (kPointsTypeStr[l] == pointsStrings[k])
+                        {
+                            def->m_points.push_back((PointsType)l);
+                            valid = true;
+                            break;
+                        }
+                    }
+
+                    ASSERTL0(valid,
+                        std::string(
+                            "Unable to correctly parse the points type: ")
+                            .append(pointsStrings[k])
+                            .c_str());
+                }
+            }
+            else if (attrName == ATTRNAME_NUMPOINTS)
+            {
+                def->m_numPointsDef = true;
+            
+                f >> attrValue;
+
+                bool valid = ParseUtils::GenerateVector(
+                    attrValue.c_str(), def->m_numPoints);
+                ASSERTL0(valid,
+                    "Unable to correctly parse the number of points.");
+            }
+            else
+            {
+                std::string errstr("Unknown field attribute: ");
+                errstr += attrName;
+                ASSERTL1(false, errstr.c_str());
+            }
+        }
+        
+        // element IDs
+        ////////////////////////////////////////////////
+        f >> def->m_elementIDs;
+        ////////////////////////////////////////////////
+        
+        // element data
+        ////////////////////////////////////////////////
+        ASSERTL0(fielddata != NullVectorNekDoubleVector,
+            "Null fielddata.");
+        
+        std::vector<NekDouble> data;
+        f >> data;
+        
+        int datasize = CheckFieldDefinition(def);
+        ASSERTL0(data.size() == datasize*def->m_fields.size(),
+            "Input data is not the same length as header information.");
+        ////////////////////////////////////////////////
+
+        fielddefs.push_back(def);
+        fielddata.push_back(data);
+    }
+
+    unsigned long nRead = f.getBytesRead();
+    f.close();
+
+    m_comm->Block();
+    
+    return nRead;
+}
+
+  
+DataSourceSharedPtr FieldIOSIONlib::v_ImportFieldMetaData(
+    const std::string &filename,
+    FieldMetaDataMap &fieldmetadatamap)
+{
+    DataSourceSharedPtr ans = SIONlibDataSource::create(filename);
+    
+    return ans;
+}
+
+
+}
+}
diff --git a/library/LibUtilities/BasicUtils/FieldIOSIONlib.h b/library/LibUtilities/BasicUtils/FieldIOSIONlib.h
new file mode 100755
index 0000000000000000000000000000000000000000..20636cd622054266d7c6244ebfd7a2b95dda37b7
--- /dev/null
+++ b/library/LibUtilities/BasicUtils/FieldIOSIONlib.h
@@ -0,0 +1,177 @@
+///////////////////////////////////////////////////////////////////////////////
+//
+// File FieldIOSIONlib.h
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: Field IO to/from SIONlib files.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_FIELDIOSIONLIB_H
+#define NEKTAR_LIB_UTILITIES_BASIC_UTILS_FIELDIOSIONLIB_H
+
+#include <LibUtilities/BasicUtils/SessionReader.h>
+#include <LibUtilities/BasicUtils/FieldIO.h>
+#include <LibUtilities/BasicUtils/FileSystem.h>
+
+namespace Nektar
+{
+namespace LibUtilities
+{
+
+/**
+ * @class Class encapsulating simple SIONlib data source.
+ */
+class SIONlibDataSource : public DataSource
+{
+public:
+    /// Constructor based on filename.
+    SIONlibDataSource(const std::string &fn)
+    {
+    }
+
+    std::string Get()
+    {
+        return fn;
+    }
+
+    const std::string Get() const
+    {
+        return fn;
+    }
+
+    /// Static constructor for this data source.
+    static DataSourceSharedPtr create(const std::string &fn)
+    {
+        return DataSourceSharedPtr(new SIONlibDataSource(fn));
+    }
+
+private:
+    /// SION document.
+    std::string fn;
+};
+typedef std::shared_ptr<SIONlibDataSource> SIONlibDataSourceSharedPtr;
+
+ 
+class FieldIOSIONlib : public FieldIO
+{
+public:
+    static const unsigned int FORMAT_VERSION;
+    
+    /// Creates an instance of this class
+    LIB_UTILITIES_EXPORT static FieldIOSharedPtr create(
+        LibUtilities::CommSharedPtr pComm, bool sharedFilesystem)
+    {
+        return MemoryManager<FieldIOSIONlib>::AllocateSharedPtr(pComm,
+                                                             sharedFilesystem);
+    }
+
+    /// Name of class
+    LIB_UTILITIES_EXPORT static std::string className;
+
+    LIB_UTILITIES_EXPORT FieldIOSIONlib(
+        LibUtilities::CommSharedPtr pComm,
+        bool sharedFilesystem);
+
+    LIB_UTILITIES_EXPORT virtual ~FieldIOSIONlib()
+    {
+    }
+    
+    /// Get class name
+    inline virtual const std::string &GetClassName() const
+    {
+        return className;
+    }
+
+private:
+    static const std::string ATTRNAME_FIELDS;
+    static const std::string ATTRNAME_BASIS;
+    static const std::string ATTRNAME_SHAPE;
+    static const std::string ATTRNAME_HOMOLENS;
+    static const std::string ATTRNAME_NUMMODES;
+    static const std::string ATTRNAME_POINTSTYPE;
+    static const std::string ATTRNAME_NUMPOINTS;
+    
+    static const std::string ATTRVALUE_MIXORDER;
+    static const std::string ATTRVALUE_UNIORDER;
+
+    static sion_int32 block_size;
+    static sion_int64 chunk_size;
+    static std::string read_mode;
+    static std::string write_mode;
+    
+    SIONlib::SIONFile *OpenFileForWriting(const std::string &outFilename);
+
+    LIB_UTILITIES_EXPORT virtual void v_Init(
+        const LibUtilities::SessionReaderSharedPtr session);
+
+    LIB_UTILITIES_EXPORT virtual void v_InitFromBenchmarker(
+        const LibUtilities::IOSettingsSharedPtr iosettings);
+
+
+    void WriteRawDataToBinaryVector(std::vector<NekByte> &v, const NekByte* dptr, const size_t dsize);
+    void WriteStringToBinaryVector(std::vector<NekByte> &v, const std::string &s);
+
+    template<class T>   
+    void WriteDatumToBinaryVector(std::vector<NekByte> &v, const T d);
+
+    template<class T>   
+    void WriteDataToBinaryVector(std::vector<NekByte> &v, const std::vector<T> &d);
+ 
+    LIB_UTILITIES_EXPORT virtual unsigned long v_Write(
+        const std::string &outFileSuffix,
+        std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+        std::vector<std::vector<NekDouble> > &fielddata,
+        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+        const bool backup = false);
+
+    
+    SIONlib::SIONFile *OpenFileForReading(const std::string &inFilename);
+
+    unsigned long DirectImport(SIONlib::SIONFile &f,
+        std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+        std::vector<std::vector<NekDouble> > &fielddata,
+        FieldMetaDataMap &fieldinfomap,
+        const Array<OneD, int> &ElementiDs);
+        
+    LIB_UTILITIES_EXPORT virtual unsigned long v_Import(const std::string &infilename,
+        std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+        std::vector<std::vector<NekDouble> > &fielddata = NullVectorNekDoubleVector,
+        FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+        const Array<OneD, int> &ElementiDs = NullInt1DArray);
+
+    
+    LIB_UTILITIES_EXPORT virtual DataSourceSharedPtr v_ImportFieldMetaData(
+        const std::string &filename, FieldMetaDataMap &fieldmetadatamap);
+};
+
+}
+}
+
+#endif
diff --git a/library/LibUtilities/BasicUtils/FieldIOXml.cpp b/library/LibUtilities/BasicUtils/FieldIOXml.cpp
index 7030b7556362a8ea4d5aa4856c472711398f89b6..5fe6ed3c3802e44ef4cd4b18a81f430030d010b6 100644
--- a/library/LibUtilities/BasicUtils/FieldIOXml.cpp
+++ b/library/LibUtilities/BasicUtils/FieldIOXml.cpp
@@ -84,22 +84,25 @@ FieldIOXml::FieldIOXml(LibUtilities::CommSharedPtr pComm, bool sharedFilesystem)
  * @param fielddefs         Input field definitions.
  * @param fielddata         Input field data.
  * @param fieldmetadatamap  Field metadata.
+ * @return The number of bytes written.
  */
-void FieldIOXml::v_Write(const std::string &outFile,
-                         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                         std::vector<std::vector<NekDouble> > &fielddata,
-                         const FieldMetaDataMap &fieldmetadatamap,
-                         const bool backup)
+uint64_t FieldIOXml::v_Write(const std::string &outFile,
+                             std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                             std::vector<std::vector<NekDouble> > &fielddata,
+                             const FieldMetaDataMap &fieldmetadatamap,
+                             const bool backup)
 {
+    uint64_t nWritten = 0;
+    
     double tm0 = 0.0, tm1 = 0.0;
     if (m_comm->TreatAsRankZero())
     {
         tm0 = m_comm->Wtime();
     }
-
+    
     // Check everything seems sensible
     ASSERTL1(fielddefs.size() == fielddata.size(),
-             "Length of fielddefs and fielddata incompatible");
+             "Length of fielddefs and fielddata incompatible.");
     for (int f = 0; f < fielddefs.size(); ++f)
     {
         ASSERTL1(fielddata[f].size() > 0,
@@ -152,6 +155,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
             }
             fieldsString = fieldsStringStream.str();
         }
+        nWritten += fieldsString.size();
         elemTag->SetAttribute("FIELDS", fieldsString);
 
         // Write SHAPE
@@ -175,6 +179,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
 
             shapeString = shapeStringStream.str();
         }
+        nWritten += shapeString.size();
         elemTag->SetAttribute("SHAPE", shapeString);
 
         // Write BASIS
@@ -195,6 +200,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
             }
             basisString = basisStringStream.str();
         }
+        nWritten += basisString.size();
         elemTag->SetAttribute("BASIS", basisString);
 
         // Write homogeneuous length details
@@ -214,6 +220,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
                 }
                 homoLenString = homoLenStringStream.str();
             }
+            nWritten += homoLenString.size();
             elemTag->SetAttribute("HOMOGENEOUSLENGTHS", homoLenString);
         }
 
@@ -239,6 +246,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
                     }
                     homoYIDsString = homoYIDsStringStream.str();
                 }
+                nWritten += homoYIDsString.size();
                 elemTag->SetAttribute("HOMOGENEOUSYIDS", homoYIDsString);
             }
 
@@ -261,6 +269,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
                     }
                     homoZIDsString = homoZIDsStringStream.str();
                 }
+                nWritten += homoZIDsString.size();
                 elemTag->SetAttribute("HOMOGENEOUSZIDS", homoZIDsString);
             }
 
@@ -283,6 +292,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
                     }
                     homoSIDsString = homoSIDsStringStream.str();
                 }
+                nWritten += homoSIDsString.size();
                 elemTag->SetAttribute("HOMOGENEOUSSIDS", homoSIDsString);
             }
         }
@@ -328,6 +338,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
 
             numModesString = numModesStringStream.str();
         }
+        nWritten += numModesString.size();
         elemTag->SetAttribute("NUMMODESPERDIR", numModesString);
 
         // Write ID
@@ -338,20 +349,26 @@ void FieldIOXml::v_Write(const std::string &outFile,
             std::stringstream idStringStream;
             idString = ParseUtils::GenerateSeqString(fielddefs[f]->m_elementIDs);
         }
+        nWritten += idString.size();
         elemTag->SetAttribute("ID", idString);
+
+        std::string compressedString = LibUtilities::CompressData::GetCompressString();
+        nWritten += compressedString.size();
         elemTag->SetAttribute("COMPRESSED",
-                              LibUtilities::CompressData::GetCompressString());
+                              compressedString);
 
         // Add this information for future compatibility
         // issues, for exmaple in case we end up using a 128
         // bit machine.
-        elemTag->SetAttribute("BITSIZE",
-                              LibUtilities::CompressData::GetBitSizeStr());
+        std::string bitSizeString = LibUtilities::CompressData::GetBitSizeStr();
+        nWritten += bitSizeString.size();
+        elemTag->SetAttribute("BITSIZE" ,bitSizeString);
+        
         std::string base64string;
         ASSERTL0(Z_OK == CompressData::ZlibEncodeToBase64Str(fielddata[f],
                                                              base64string),
                  "Failed to compress field data.");
-
+        nWritten += base64string.size();
         elemTag->LinkEndChild(new TiXmlText(base64string));
     }
     doc.SaveFile(filename);
@@ -359,11 +376,13 @@ void FieldIOXml::v_Write(const std::string &outFile,
     m_comm->Block();
 
     // all data has been written
-    if (m_comm->TreatAsRankZero())
-    {
-        tm1 = m_comm->Wtime();
-        cout << " (" << tm1 - tm0 << "s, XML)" << endl;
-    }
+    //if (m_comm->TreatAsRankZero())
+    //{
+    //    tm1 = m_comm->Wtime();
+    //    cout << " (" << tm1 - tm0 << "s, XML)" << endl;
+    //}
+    
+    return nWritten;
 }
 
 /**
@@ -389,7 +408,7 @@ void FieldIOXml::WriteMultiFldFileIDs(
     doc.LinkEndChild(decl);
 
     ASSERTL0(fileNames.size() == elementList.size(),
-             "Outfile names and list of elements ids does not match");
+             "Outfile names and list of elements ids does not match.");
 
     TiXmlElement *root = new TiXmlElement("NEKTAR");
     doc.LinkEndChild(root);
@@ -426,7 +445,7 @@ void FieldIOXml::WriteMultiFldFileIDs(
  * @param elementList        Vector of element IDs that lie on each process.
  * @param fieldmetadatamap   Field metadata map that is read from @p inFile.
  */
-void FieldIOXml::ImportMultiFldFileIDs(
+void FieldIOXml::v_ImportMultiFldFileIDs(
     const std::string                       &inFile,
     std::vector<std::string>                &fileNames,
     std::vector<std::vector<unsigned int> > &elementList,
@@ -480,7 +499,7 @@ void FieldIOXml::ImportMultiFldFileIDs(
         std::string elementIDsStr(elementIDs);
 
         std::vector<unsigned int> idvec;
-        ParseUtils::GenerateSeqVector(elementIDsStr, idvec);
+        ParseUtils::GenerateSeqVector(elementIDsStr.c_str(), idvec);
 
         elementList.push_back(idvec);
 
@@ -499,13 +518,16 @@ void FieldIOXml::ImportMultiFldFileIDs(
  *                          this rank. The resulting field definitions will only
  *                          contain data for the element IDs specified in this
  *                          array.
+ * @return The number of bytes read.
  */
-void FieldIOXml::v_Import(const std::string &infilename,
-                          std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-                          std::vector<std::vector<NekDouble> > &fielddata,
-                          FieldMetaDataMap &fieldinfomap,
-                          const Array<OneD, int> &ElementIDs)
+uint64_t FieldIOXml::v_Import(const std::string &infilename,
+                              std::vector<FieldDefinitionsSharedPtr> &fielddefs,
+                              std::vector<std::vector<NekDouble> > &fielddata,
+                              FieldMetaDataMap &fieldinfomap,
+                              const Array<OneD, int> &ElementIDs)
 {
+    uint64_t nRead = 0;
+    
     std::string infile = infilename;
 
     fs::path pinfilename(infilename);
@@ -535,10 +557,10 @@ void FieldIOXml::v_Import(const std::string &infilename,
                 fullpath                       = pinfilename / pfilename;
                 string fname                   = PortablePath(fullpath);
                 DataSourceSharedPtr dataSource = XmlDataSource::create(fname);
-                ImportFieldDefs(dataSource, fielddefs, false);
+                nRead += ImportFieldDefs(dataSource, fielddefs, false);
                 if (fielddata != NullVectorNekDoubleVector)
                 {
-                    ImportFieldData(dataSource, fielddefs, fielddata);
+                    nRead += ImportFieldData(dataSource, fielddefs, fielddata);
                 }
             }
         }
@@ -576,10 +598,10 @@ void FieldIOXml::v_Import(const std::string &infilename,
                 fullpath                       = pinfilename / pfilename;
                 string fname                   = PortablePath(fullpath);
                 DataSourceSharedPtr dataSource = XmlDataSource::create(fname);
-                ImportFieldDefs(dataSource, fielddefs, false);
+                nRead += ImportFieldDefs(dataSource, fielddefs, false);
                 if (fielddata != NullVectorNekDoubleVector)
                 {
-                    ImportFieldData(dataSource, fielddefs, fielddata);
+                    nRead += ImportFieldData(dataSource, fielddefs, fielddata);
                 }
             }
         }
@@ -588,12 +610,15 @@ void FieldIOXml::v_Import(const std::string &infilename,
     {
         // serial format case
         DataSourceSharedPtr doc = ImportFieldMetaData(infilename, fieldinfomap);
-        ImportFieldDefs(doc, fielddefs, false);
+        nRead += ImportFieldDefs(doc, fielddefs, false);
         if (fielddata != NullVectorNekDoubleVector)
         {
-            ImportFieldData(doc, fielddefs, fielddata);
+            nRead += ImportFieldData(doc, fielddefs, fielddata);
         }
     }
+
+    m_comm->Block();
+    return nRead;
 }
 
 /**
@@ -765,12 +790,14 @@ void FieldIOXml::SetUpFieldMetaData(
  * @param fielddefs   Output vector that will contain read field definitions.
  * @param expChild    Determines if the field definitions are defined by
  *                    `<EXPANSIONS>` or in `<NEKTAR>`.
+ * @return The number of bytes read.
  */
-void FieldIOXml::ImportFieldDefs(
+uint64_t FieldIOXml::ImportFieldDefs(
     DataSourceSharedPtr dataSource,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     bool expChild)
 {
+    uint64_t nRead = 0;
     XmlDataSourceSharedPtr xml =
         std::static_pointer_cast<XmlDataSource>(dataSource);
     TiXmlElement *master =
@@ -884,6 +911,8 @@ void FieldIOXml::ImportFieldDefs(
                     ASSERTL1(false, errstr.c_str());
                 }
 
+                nRead += attr->ValueStr().size();
+                
                 // Get the next attribute.
                 attr = attr->Next();
             }
@@ -1091,6 +1120,8 @@ void FieldIOXml::ImportFieldDefs(
         }
         loopXml = loopXml->NextSiblingElement(strLoop);
     }
+
+    return nRead;
 }
 
 /**
@@ -1099,12 +1130,15 @@ void FieldIOXml::ImportFieldDefs(
  * @param dataSource  Target XML file
  * @param fielddefs   Field definitions for file
  * @param fielddata   On return, contains field data for each field.
+ * @return The number of bytes read.
  */
-void FieldIOXml::ImportFieldData(
+uint64_t FieldIOXml::ImportFieldData(
     DataSourceSharedPtr dataSource,
     const std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     std::vector<std::vector<NekDouble> > &fielddata)
 {
+    uint64_t nRead = 0;
+    
     int cntdumps = 0;
     XmlDataSourceSharedPtr xml =
         std::static_pointer_cast<XmlDataSource>(dataSource);
@@ -1150,6 +1184,8 @@ void FieldIOXml::ImportFieldData(
                           " but got " + string(CompressStr));
             }
 
+            nRead += elementStr.size();
+            
             ASSERTL0(Z_OK == CompressData::ZlibDecodeFromBase64Str(
                                  elementStr, elementFieldData),
                      "Failed to decompress field data.");
@@ -1159,14 +1195,17 @@ void FieldIOXml::ImportFieldData(
             ASSERTL0(
                 fielddata[cntdumps].size() ==
                     datasize * fielddefs[cntdumps]->m_fields.size(),
-                "Input data is not the same length as header infoarmation");
+                "Input data is not the same length as header information");
 
             cntdumps++;
 
             element = element->NextSiblingElement("ELEMENTS");
         }
         master = master->NextSiblingElement("NEKTAR");
+
     }
+
+    return nRead;
 }
 }
 }
diff --git a/library/LibUtilities/BasicUtils/FieldIOXml.h b/library/LibUtilities/BasicUtils/FieldIOXml.h
index fe7360a8f7326743c0cb22379fcbbd2dd4758325..607b6034750ea81a68a588fc5eeb65fb967285b9 100644
--- a/library/LibUtilities/BasicUtils/FieldIOXml.h
+++ b/library/LibUtilities/BasicUtils/FieldIOXml.h
@@ -220,13 +220,13 @@ public:
     LIB_UTILITIES_EXPORT virtual ~FieldIOXml()
     {
     }
-
-    LIB_UTILITIES_EXPORT void ImportFieldDefs(
+    
+    LIB_UTILITIES_EXPORT uint64_t ImportFieldDefs(
         DataSourceSharedPtr dataSource,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         bool expChild);
 
-    LIB_UTILITIES_EXPORT void ImportFieldData(
+    LIB_UTILITIES_EXPORT uint64_t ImportFieldData(
         DataSourceSharedPtr dataSource,
         const std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata);
@@ -240,21 +240,21 @@ public:
     LIB_UTILITIES_EXPORT void SetUpFieldMetaData(
         const std::string &outname,
         const std::vector<FieldDefinitionsSharedPtr> &fielddefs,
-        const FieldMetaDataMap &fieldmetadatamap);
-
-    LIB_UTILITIES_EXPORT void ImportMultiFldFileIDs(
-        const std::string &inFile,
-        std::vector<std::string> &fileNames,
-        std::vector<std::vector<unsigned int> > &elementList,
-        FieldMetaDataMap &fieldmetadatamap);
-
-    LIB_UTILITIES_EXPORT void v_Import(
+        const FieldMetaDataMap &fieldmetadatamap);    
+        
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Import(
         const std::string &infilename,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata =
             NullVectorNekDoubleVector,
         FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
-        const Array<OneD, int> &ElementIDs = NullInt1DArray);
+        const Array<OneD, int> &ElementiDs = NullInt1DArray);
+        
+    LIB_UTILITIES_EXPORT virtual void v_ImportMultiFldFileIDs(
+        const std::string &inFile,
+        std::vector<std::string> &fileNames,
+        std::vector<std::vector<unsigned int> > &elementList,
+        FieldMetaDataMap &fieldmetadatamap);
 
     /// Returns the class name.
     inline virtual const std::string &GetClassName() const
@@ -262,8 +262,8 @@ public:
         return className;
     }
 
-private:
-    LIB_UTILITIES_EXPORT virtual void v_Write(
+private:  
+    LIB_UTILITIES_EXPORT virtual uint64_t v_Write(
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
diff --git a/library/LibUtilities/BasicUtils/H5.cpp b/library/LibUtilities/BasicUtils/H5.cpp
index b9e62170ad8745bfe0a0c495113c17c171cca8cf..bf486ed52852a2f7df9ed91f9c0201e118a5d56b 100644
--- a/library/LibUtilities/BasicUtils/H5.cpp
+++ b/library/LibUtilities/BasicUtils/H5.cpp
@@ -465,6 +465,12 @@ void DataSpace::SelectRange(const hsize_t start, const hsize_t count)
             (m_Id, H5S_SELECT_SET, &start, NULL, &count, NULL));
 }
 
+void DataSpace::SelectRange(const H5S_seloper_t op, const hsize_t start, const hsize_t count)
+{
+    H5_CALL(H5Sselect_hyperslab,
+            (m_Id, op, &start, NULL, &count, NULL));
+}
+
 hsize_t DataSpace::GetSize()
 {
     return H5Sget_simple_extent_npoints(m_Id);
diff --git a/library/LibUtilities/BasicUtils/H5.h b/library/LibUtilities/BasicUtils/H5.h
index 0cdccfc83a3a3a6055b6a9602c99c78c340794d4..342fccc0f713f208eeaefab5396ca76e02739061 100644
--- a/library/LibUtilities/BasicUtils/H5.h
+++ b/library/LibUtilities/BasicUtils/H5.h
@@ -326,6 +326,7 @@ public:
 
     void Close();
     void SelectRange(const hsize_t start, const hsize_t count);
+    void SelectRange(const H5S_seloper_t op, const hsize_t start, const hsize_t count);
 
     hsize_t GetSize();
     std::vector<hsize_t> GetDims();
diff --git a/library/LibUtilities/BasicUtils/SIONlib.cpp b/library/LibUtilities/BasicUtils/SIONlib.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..8995c29afcd174cf41a02595ccc4821b9712791e
--- /dev/null
+++ b/library/LibUtilities/BasicUtils/SIONlib.cpp
@@ -0,0 +1,624 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: SIONlib.cpp
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: Minimal SIONlib wrapper
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include <LibUtilities/BasicUtils/SIONlib.h>
+
+
+namespace Nektar
+{
+namespace LibUtilities
+{
+namespace SIONlib
+{
+
+const unsigned int SION_Base::LUSTRE_STRIPE_BASE_SIZE = 65536;  
+
+// Common functions
+char *SION_Base::getSionFileName() const
+{
+    return _sion_file_name;
+}
+
+void SION_Base::setDebug(bool debug)
+{
+    _debug = debug;
+}
+  
+void SION_Base::setMode(std::string mode)
+{
+    _mode = mode;
+}
+
+std::string SION_Base::getMode() const
+{
+    return _mode;
+}
+
+bool SION_Base::isModeCollective() const
+{
+    return _is_mode_collective;
+}
+
+bool SION_Base::isModeCollectiveMerge() const
+{
+    return _is_mode_collectivemerge;
+}
+
+
+void SION_Base::setNumberOfFiles(int num_files)
+{
+    _num_files = num_files;
+}
+
+int SION_Base::getNumberOfFiles() const
+{
+    return _num_files;
+}
+
+void SION_Base::setNumberOfTasks(int num_tasks)
+{
+    _num_tasks = num_tasks;
+}
+
+int SION_Base::getNumberOfTasks() const
+{
+    return _num_tasks;
+}
+
+void SION_Base::setRank(int rank)
+{
+    _rank = rank;
+}
+
+int SION_Base::getRank() const
+{
+    return _rank;
+}
+
+void SION_Base::setChunkSize(sion_int64 chunk_size)
+{
+    _chunk_size = chunk_size;
+}
+
+sion_int64 SION_Base::getChunkSize() const
+{
+    return _chunk_size;
+}
+
+void SION_Base::setChunkSizes(sion_int64 *chunk_sizes)
+{
+    _chunk_sizes = chunk_sizes;
+}
+
+sion_int64 *SION_Base::getChunkSizes() const
+{
+    return _chunk_sizes;
+}
+
+void SION_Base::setGlobalRanks(int *global_ranks)
+{
+    _global_ranks = global_ranks;
+}
+
+int *SION_Base::getGlobalRanks() const
+{
+    return _global_ranks;
+}
+
+void SION_Base::setFileSystemBlockSize(sion_int32 fs_blk_size)
+{
+    _fs_blk_size = fs_blk_size;
+}
+
+sion_int32 SION_Base::getFileSystemBlockSize() const
+{
+    return _fs_blk_size;
+}
+
+int SION_Base::getNumberOfSuccessfulReadElements() const
+{
+    return _number_of_elements_sucessfully_read;
+}
+
+int SION_Base::getReturnCode() const
+{
+    return _return_code;
+}
+
+int SION_Base::getSid() const
+{
+    return _sid;
+}
+
+
+/* get information (with sion datatypes) */
+int SION_Base::getFileEndianness() const
+{
+    return sion_get_file_endianness(_sid);
+}
+
+sion_int64 SION_Base::getBytesWritten() const
+{
+    return sion_get_bytes_written(_sid);
+}
+
+sion_int64 SION_Base::getBytesRead() const
+{
+    return sion_get_bytes_read(_sid);
+}
+
+sion_int64 SION_Base::getBytesAvailInBlock() const
+{
+    return sion_bytes_avail_in_block(_sid);
+}
+
+sion_int64 SION_Base::getBytesAvailInChunk() const
+{
+    return sion_bytes_avail_in_chunk(_sid);
+}
+
+sion_int64 SION_Base::getPosition() const
+{
+    return sion_get_position(_sid);
+}
+
+void SION_Base::seek()
+{
+    _return_code = sion_seek(_sid,
+        SION_CURRENT_RANK, SION_CURRENT_BLK, SION_CURRENT_POS);
+}
+
+  
+
+SIONFile::SIONFile(std::string sion_file_name, std::string mode,
+    int num_files, sion_int64 chunk_size, sion_int32 block_size,
+    int global_rank, MPI_Comm gComm, MPI_Comm lComm)
+{
+    size_t ncharacter = sion_file_name.length()+1;
+    const char *tmp_sion_file_name = sion_file_name.c_str();
+    _sion_file_name = new char[ncharacter];
+
+    strncpy(_sion_file_name, tmp_sion_file_name, ncharacter);
+
+    _debug = false;
+
+    _mode = mode;
+
+    _is_mode_collective = false;
+    _is_mode_collectivemerge = false;
+    std::istringstream mode_stream(mode);
+    std::string mode_comp;
+    while (getline(mode_stream, mode_comp, ','))
+    {
+        if (!_is_mode_collective) _is_mode_collective = (0 == mode_comp.compare("collective"));
+        if (!_is_mode_collectivemerge) _is_mode_collectivemerge = (0 == mode_comp.compare("collectivemerge"));
+    }
+    
+    _num_files = num_files;
+
+    _global_rank = global_rank;
+    
+    _g_comm = gComm;
+    _l_comm = lComm;
+
+    const sion_int32 base_size = LUSTRE_STRIPE_BASE_SIZE;
+        
+    if (block_size > 0)
+    {
+        sion_int32 base_cnt = block_size / base_size;
+
+        if (0 == base_cnt)
+        {
+            // block_size can never be smaller than base_size
+            block_size = base_size;
+        }
+        else
+        {
+            if (0 != block_size % base_size)
+            {
+                // block_size must be a multiple of base_size
+                base_cnt += 1;
+                block_size = base_size*base_cnt;  
+            }
+
+            if (base_cnt > 1 && 0 != base_cnt % 2)
+            {
+                // if block_size > base_size then
+                // block_size must be an even multiple of base_size
+                base_cnt += 1;
+                block_size = base_size*base_cnt;
+                // this ensures that Lustre file stripe size can be
+                // set equal to the SIONlib file system block size
+            }
+        }
+
+        _fs_blk_size = block_size;
+    }
+    else
+    {
+        block_size = base_size;
+        _fs_blk_size = -1;
+    }
+    
+
+    if (chunk_size > 0)
+    {
+        if (block_size > chunk_size)
+        {
+            // set chunk_size such that block_size is a multiple of chunk_size
+            if (0 != block_size % chunk_size)
+            {
+                sion_int32 chunk_cnt = block_size / chunk_size + 1;
+                chunk_size = block_size / chunk_cnt;
+            }
+        }
+        else if (chunk_size > block_size)
+        {
+            // force chunk_size to be a multiple of block_size
+            if (0 != chunk_size % block_size)
+            {
+                sion_int32 block_cnt = chunk_size / block_size + 1;
+                chunk_size = block_size*block_cnt;
+            }
+        }
+    }
+    else
+    {
+        chunk_size = base_size;
+    }
+    _chunk_size = chunk_size;
+
+    
+    _file_ptr = NULL;
+    _new_sion_file_name = new char[255];
+}
+
+SIONFile::~SIONFile()
+{
+    delete [] _sion_file_name;
+    _sion_file_name = NULL;
+    delete [] _new_sion_file_name;
+    _new_sion_file_name = NULL;
+}
+
+
+void SIONFile::setLocalCommunicator(MPI_Comm lComm)
+{
+    _l_comm = lComm;
+}
+
+MPI_Comm SIONFile::getLocalCommunicator() const
+{
+    return _l_comm;
+}
+
+void SIONFile::setGlobalCommunicator(MPI_Comm gComm)
+{
+    _g_comm = gComm;
+}
+
+MPI_Comm SIONFile::getGlobalCommunicator() const
+{
+    return _g_comm;
+}
+
+void SIONFile::setGlobalRank(int global_rank)
+{
+    _global_rank = global_rank;
+}
+
+int SIONFile::getGlobalRank() const
+{
+    return _global_rank;
+}
+
+char *SIONFile::getNewSionFileName() const
+{
+    return _new_sion_file_name;
+}
+
+void SIONFile::open()
+{
+    _sid = sion_paropen_mpi(_sion_file_name, _mode.c_str(), &_num_files,
+        _g_comm, &_l_comm, &_chunk_size, &_fs_blk_size,
+        &_global_rank, NULL, &_new_sion_file_name);
+    
+    _return_code = _sid;
+}
+
+void SIONFile::close()
+{
+    _return_code = sion_parclose_mpi(_sid);
+}
+
+    
+void SIONFile::ensureFreeSpace(long numbytes)
+{
+    _return_code = sion_ensure_free_space(_sid, numbytes);
+}
+
+void SIONFile::endOfFile()
+{
+    _return_code = sion_feof(_sid);
+}
+
+
+template<class T>
+void SIONFile::write(const T rhs)
+{
+    if (_is_mode_collective)
+    {
+        _return_code = sion_coll_fwrite(&rhs, sizeof(rhs), 1, _sid);
+    }
+    else
+    {
+        _return_code = sion_fwrite(&rhs, sizeof(rhs), 1, _sid);
+    }
+}
+
+SIONFile &operator<<(SIONFile &sf, const LibUtilities::BasisType &rhs)
+{
+    sf.write(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const unsigned int &rhs)
+{
+    sf.write(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const size_t &rhs)
+{
+    sf.write(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const double &rhs)
+{
+    sf.write(rhs);
+    return sf;
+}
+
+
+void SIONFile::write_string(const std::string &rhs)
+{
+    size_t rhs_size = rhs.size();
+    
+    if (_is_mode_collective)
+    {
+        _return_code = sion_coll_fwrite(&rhs_size, sizeof(rhs_size), 1, _sid);
+        _return_code = sion_coll_fwrite(rhs.c_str(), sizeof(char), rhs_size, _sid);
+    }
+    else
+    {
+        _return_code = sion_fwrite(&rhs_size, sizeof(rhs_size), 1, _sid);
+        _return_code = sion_fwrite(rhs.c_str(), sizeof(char), rhs_size, _sid);
+    }
+}
+  
+SIONFile &operator<<(SIONFile &sf, const std::string &rhs)
+{
+    sf.write_string(rhs);
+    return sf;
+}
+
+
+template<class T>
+void SIONFile::write_vector(const std::vector<T> &rhs)
+{
+    size_t rhs_size = rhs.size();
+    
+    if (_is_mode_collective)
+    {
+        _return_code = sion_coll_fwrite(&rhs_size, sizeof(rhs_size), 1, _sid);
+        _return_code = sion_coll_fwrite(rhs.data(), sizeof(T), rhs_size, _sid);
+    }
+    else
+    {
+        _return_code = sion_fwrite(&rhs_size, sizeof(rhs_size), 1, _sid);
+        _return_code = sion_fwrite(rhs.data(), sizeof(T), rhs_size, _sid);
+    }
+}
+
+SIONFile &operator<<(SIONFile &sf, const std::vector<LibUtilities::BasisType> &rhs)
+{
+    sf.write_vector(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const std::vector<unsigned int> &rhs)
+{
+    sf.write_vector(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const std::vector<double> &rhs)
+{
+    sf.write_vector(rhs);
+    return sf;
+}
+
+SIONFile &operator<<(SIONFile &sf, const std::vector<unsigned char> &rhs)
+{
+    sf.write_vector(rhs);
+    return sf;
+}
+
+template<class T>
+void SIONFile::read(T *rhs)
+{
+    if (NULL != rhs)
+    {
+        if (_is_mode_collective)
+        {
+            _return_code = sion_coll_fread(rhs, sizeof(T), 1, _sid);
+        }
+        else
+        {
+            _return_code = sion_fread(rhs, sizeof(T), 1, _sid);
+        }
+    }
+}
+
+SIONFile &operator>>(SIONFile &sf, LibUtilities::BasisType &rhs)
+{
+    sf.read(&rhs);
+    return sf;
+}
+
+SIONFile &operator>>(SIONFile &sf, unsigned int &rhs)
+{
+    sf.read(&rhs);
+    return sf;
+}
+
+SIONFile &operator>>(SIONFile &sf, size_t &rhs)
+{
+    sf.read(&rhs);
+    return sf;
+}
+
+SIONFile &operator>>(SIONFile &sf, double &rhs)
+{
+    sf.read(&rhs);
+    return sf;
+}
+
+
+void SIONFile::read_string(std::string &rhs)
+{
+    size_t rhs_size = 0;
+
+    if (_is_mode_collective)
+    {
+        _return_code = sion_coll_fread(&rhs_size, sizeof(rhs_size), 1, _sid);
+    }
+    else
+    {
+        _return_code = sion_fread(&rhs_size, sizeof(rhs_size), 1, _sid);
+    }
+        
+    if (1 == _return_code && rhs_size > 0)
+    {
+        char *rhs_buf = new char[rhs_size];
+        if (NULL != rhs_buf)
+        {
+            if (_is_mode_collective)
+            {
+                _return_code = sion_coll_fread(rhs_buf, sizeof(char), rhs_size, _sid);
+            }
+            else
+            {
+                _return_code = sion_fread(rhs_buf, sizeof(char), rhs_size, _sid);
+            }
+
+            if (rhs_size == _return_code)
+            {
+                rhs.resize(rhs_size);
+                rhs.assign(rhs_buf, rhs_size);
+            }
+            
+            delete [] rhs_buf;
+        }
+    }
+}
+
+SIONFile &operator>>(SIONFile &sf, std::string &rhs)
+{
+    sf.read_string(rhs);
+    return sf;
+}
+
+
+template<class T>
+void SIONFile::read_vector(std::vector<T> &rhs)
+{
+    size_t rhs_size = 0;
+    if (_is_mode_collective)
+    {
+        _return_code = sion_coll_fread(&rhs_size, sizeof(rhs_size), 1, _sid);
+    }
+    else
+    {
+        _return_code = sion_fread(&rhs_size, sizeof(rhs_size), 1, _sid);
+    }
+
+    if (1 == _return_code && rhs_size > 0)
+    {
+        T *rhs_buf = new T[rhs_size];
+        if (NULL != rhs_buf)
+        {
+            if (_is_mode_collective)
+            {
+                _return_code = sion_coll_fread(rhs_buf, sizeof(T), rhs_size, _sid);
+            }
+            else
+            {
+                _return_code = sion_fread(rhs_buf, sizeof(T), rhs_size, _sid);
+            }
+
+            if (rhs_size == _return_code)
+            {
+                rhs.resize(rhs_size);
+                rhs.assign(rhs_buf, rhs_buf+rhs_size);
+            }
+            
+            delete [] rhs_buf;
+        }
+    }
+}
+
+SIONFile &operator>>(SIONFile &sf, std::vector<LibUtilities::BasisType> &rhs)
+{
+    sf.read_vector(rhs);
+    return sf;
+}
+
+SIONFile &operator>>(SIONFile &sf, std::vector<unsigned int> &rhs)
+{
+    sf.read_vector(rhs);
+    return sf;
+}
+
+SIONFile &operator>>(SIONFile &sf, std::vector<double> &rhs)
+{
+    sf.read_vector(rhs);
+    return sf;
+}
+
+}
+}
+}
diff --git a/library/LibUtilities/BasicUtils/SIONlib.h b/library/LibUtilities/BasicUtils/SIONlib.h
new file mode 100755
index 0000000000000000000000000000000000000000..5500c82dec267b56d7e2469ccc770a4fec602617
--- /dev/null
+++ b/library/LibUtilities/BasicUtils/SIONlib.h
@@ -0,0 +1,226 @@
+///////////////////////////////////////////////////////////////////////////////
+//
+// File SIONlib.h
+//
+// For more information, please see: http://www.nektar.info
+//
+// The MIT License
+//
+// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+// Department of Aeronautics, Imperial College London (UK), and Scientific
+// Computing and Imaging Institute, University of Utah (USA).
+//
+// License for the specific language governing rights and limitations under
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the "Software"),
+// to deal in the Software without restriction, including without limitation
+// the rights to use, copy, modify, merge, publish, distribute, sublicense,
+// and/or sell copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+// DEALINGS IN THE SOFTWARE.
+//
+// Description: Simple OO wrapper around SIONlib
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_SIONLIB_H
+#define NEKTAR_LIB_UTILITIES_BASIC_UTILS_SIONLIB_H
+
+#include "mpi.h"
+#include "sion.h"
+
+#include <exception>
+#include <string>
+#include <string.h>
+#include <vector>
+#include <sstream>
+#include <iostream>
+#include <stdlib.h>
+#include <fstream>
+
+#include <LibUtilities/Foundations/BasisType.h>
+
+namespace Nektar
+{
+namespace LibUtilities
+{
+namespace SIONlib
+{
+  
+class SION_Base
+{
+
+public:
+    static const unsigned int LUSTRE_STRIPE_BASE_SIZE;
+    
+    SION_Base() : _sid(-9999) {}
+    ~SION_Base() {}
+
+    // Common functions
+    char *getSionFileName() const;
+
+    void setDebug(bool debug);
+    
+    void setMode(std::string mode);
+    std::string getMode() const;
+    bool isModeCollective() const;
+    bool isModeCollectiveMerge() const;
+
+    void setNumberOfFiles(int num_files);
+    int getNumberOfFiles() const;
+
+    void setNumberOfTasks(int num_tasks);
+    int getNumberOfTasks() const;
+
+    void setRank(int rank);
+    int getRank() const;
+
+    void setChunkSize(sion_int64 chunk_size);
+    sion_int64 getChunkSize() const;
+
+    void setChunkSizes(sion_int64 *chunk_sizes);
+    sion_int64 *getChunkSizes() const;
+
+    void setGlobalRanks(int *global_ranks);
+    int *getGlobalRanks() const;
+
+    void setFileSystemBlockSize(sion_int32 fs_blk_size);
+    sion_int32 getFileSystemBlockSize() const;
+
+    int getNumberOfSuccessfulReadElements() const;
+   
+    int getSid() const;
+
+    int getReturnCode() const;
+
+    void seek();
+    
+    /* get information (with sion datatypes) */
+    int getFileEndianness() const;
+    sion_int64 getBytesWritten() const;
+    sion_int64 getBytesRead() const;
+    sion_int64 getBytesAvailInBlock() const;
+    sion_int64 getBytesAvailInChunk() const;
+    sion_int64 getPosition() const;
+
+protected:
+    // Common attributes
+    char *_sion_file_name;
+    std::string _mode;
+    bool _is_mode_collective;
+    bool _is_mode_collectivemerge;
+    int _num_files;
+    int _num_tasks;
+    int _rank;
+    sion_int64 *_chunk_sizes;
+    sion_int64 _chunk_size;
+    sion_int32 _fs_blk_size;
+    int *_global_ranks;
+    FILE *_file_ptr;
+    int _number_of_elements_sucessfully_read;
+    int _return_code;
+    int _sid;
+    bool _debug;
+
+    //  get information (with sion datatypes) 
+    int _file_endianness;
+    sion_int64 _bytes_written;
+    sion_int64 _bytes_read;
+    sion_int64 _bytes_avail_in_block;
+    sion_int64 _bytes_avail_in_chunk;
+    sion_int64 _position;
+};
+
+
+class SIONFile : public SION_Base
+{
+
+public:
+      
+    SIONFile();
+    SIONFile(std::string sion_file_name, std::string mode = "bw",
+        int num_files = 1, sion_int64 chunk_size=LUSTRE_STRIPE_BASE_SIZE, sion_int32 block_size=-1,
+        int global_rank = 0, MPI_Comm gComm=MPI_COMM_WORLD, MPI_Comm lComm=MPI_COMM_WORLD);
+    virtual ~SIONFile();
+
+    char *getNewSionFileName() const;
+    void setLocalCommunicator(MPI_Comm lComm);
+    MPI_Comm getLocalCommunicator() const;
+
+    void setGlobalCommunicator(MPI_Comm gComm);
+    MPI_Comm getGlobalCommunicator() const;
+
+    void setGlobalRank(int global_rank);
+    int getGlobalRank() const;
+
+    void open();
+    void close();
+
+    void ensureFreeSpace(long numbytes);
+    void endOfFile();
+
+
+    template<class T>
+    void write(const T rhs);
+
+    friend SIONFile &operator<<(SIONFile &sf, const LibUtilities::BasisType &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const unsigned int &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const size_t &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const double &rhs);
+
+    void write_string(const std::string &rhs);
+    
+    friend SIONFile &operator<<(SIONFile &sf, const std::string &rhs);
+
+    template<class T>
+    void write_vector(const std::vector<T> &rhs);
+
+    friend SIONFile &operator<<(SIONFile &sf, const std::vector<LibUtilities::BasisType> &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const std::vector<unsigned int> &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const std::vector<double> &rhs);
+    friend SIONFile &operator<<(SIONFile &sf, const std::vector<unsigned char> &rhs);
+
+    
+    template<class T>
+    void read(T *rhs);
+
+    friend SIONFile &operator>>(SIONFile &sf, LibUtilities::BasisType &rhs);
+    friend SIONFile &operator>>(SIONFile &sf, unsigned int &rhs);
+    friend SIONFile &operator>>(SIONFile &sf, size_t &rhs);
+    friend SIONFile &operator>>(SIONFile &sf, double &rhs);
+
+    void read_string(std::string &rhs);
+    
+    friend SIONFile &operator>>(SIONFile &sf, std::string &rhs);
+
+    template<class T>
+    void read_vector(std::vector<T> &rhs);
+
+    friend SIONFile &operator>>(SIONFile &sf, std::vector<LibUtilities::BasisType> &rhs);
+    friend SIONFile &operator>>(SIONFile &sf, std::vector<unsigned int> &rhs);
+    friend SIONFile &operator>>(SIONFile &sf, std::vector<double> &rhs);
+
+    
+private:
+    char *_new_sion_file_name;
+    MPI_Comm _g_comm;
+    MPI_Comm _l_comm;
+    int _global_rank;
+};
+
+ 
+}
+}
+}
+
+#endif
diff --git a/library/LibUtilities/CMakeLists.txt b/library/LibUtilities/CMakeLists.txt
index c2812918867b397af969e6b7e6756475fb32545b..0dc8aff59250b807e3f090d91b28bdf246c65c54 100644
--- a/library/LibUtilities/CMakeLists.txt
+++ b/library/LibUtilities/CMakeLists.txt
@@ -64,6 +64,15 @@ IF( NEKTAR_USE_HDF5 )
         ./BasicUtils/FieldIOHdf5.cpp)
 ENDIF( NEKTAR_USE_HDF5 )
 
+# SIONlib
+IF( NEKTAR_USE_SIONLIB )
+    SET(BasicUtilsHeaders ${BasicUtilsHeaders}
+        ./BasicUtils/SIONlib.h
+        ./BasicUtils/FieldIOSIONlib.h)
+    SET(BasicUtilsSources ${BasicUtilsSources}
+        ./BasicUtils/SIONlib.cpp
+        ./BasicUtils/FieldIOSIONlib.cpp)
+ENDIF( NEKTAR_USE_SIONLIB )
 
 SET(CommunicationHeaders
     ./Communication/Comm.h
@@ -424,6 +433,13 @@ IF( NEKTAR_USE_HDF5 )
     ADD_DEPENDENCIES(LibUtilities hdf5-1.8.16)
 ENDIF( NEKTAR_USE_HDF5 )
 
+# SIONlib
+IF( NEKTAR_USE_SIONLIB )
+    SET_SOURCE_FILES_PROPERTIES(./BasicUtils/FieldIOSIONlib.cpp
+        PROPERTY COMPILE_FLAGS "${SIONLIB_COMPILE_FLAGS}")
+    TARGET_LINK_LIBRARIES(LibUtilities LINK_PRIVATE ${SIONLIB_LIBRARIES})
+ENDIF( NEKTAR_USE_SIONLIB )
+
 INSTALL(DIRECTORY ./
     DESTINATION ${NEKTAR_INCLUDE_DIR}/LibUtilities
     COMPONENT dev