diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C
index 1859654e..5d4e59d8 100644
--- a/Apps/FiniteDefElasticity/main_NonLinEl.C
+++ b/Apps/FiniteDefElasticity/main_NonLinEl.C
@@ -282,7 +282,6 @@ int main (int argc, char** argv)
 
   const double epsT = 1.0e-6;
   if (dtDump <= 0.0) dtDump = params.stopTime + 1.0;
-  if (format < 0)    dtSave = params.stopTime + 1.0;
   double nextDump = params.time.t + dtDump;
   double nextSave = params.time.t + dtSave;
 
@@ -302,7 +301,7 @@ int main (int argc, char** argv)
       std::cout <<"\nWriting HDF5 file "<< infile <<".hdf5"<< std::endl;
 
     writer = new DataExporter(true);
-    writer->registerField("u","solution",DataExporter::SIM, skip2nd ? 0 : -1);
+    writer->registerField("u","displacement",DataExporter::SIM, skip2nd ? 0 : -1);
     writer->setFieldValue("u",model,(void*)&simulator.getSolution());
     writer->registerWriter(new HDF5Writer(infile));
     writer->registerWriter(new XMLWriter(infile));
@@ -337,8 +336,9 @@ int main (int argc, char** argv)
     if (params.time.t + epsT*params.time.dt > nextSave)
     {
       // Save solution variables to VTF for visualization
-      if (!simulator.saveStep(++iStep,params.time.t,n,skip2nd))
-	return 6;
+      if (format >= 0)
+        if (!simulator.saveStep(++iStep,params.time.t,n,skip2nd))
+	  return 6;
 
       // Save solution variables to HDF5
       if (writer)
diff --git a/Apps/HDF5toVTF/HDF5toVTF.C b/Apps/HDF5toVTF/HDF5toVTF.C
new file mode 100644
index 00000000..b4793f3a
--- /dev/null
+++ b/Apps/HDF5toVTF/HDF5toVTF.C
@@ -0,0 +1,138 @@
+// $Id$
+//==============================================================================
+//!
+//! \file HDF5toVTF.C
+//!
+//! \date Apr 08 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Convert a HDF5 results database to VTF for visualisation.
+//!
+//==============================================================================
+
+#include "SIM1D.h"
+#include "SIM2D.h"
+#include "SIM3D.h"
+#include "HDF5Writer.h"
+#include "XMLWriter.h"
+#include "StringUtils.h"
+#include <sstream>
+#include <stdlib.h>
+#include <string.h>
+
+typedef std::map< std::string,std::vector<XMLWriter::Entry> > ProcessList;
+
+
+int main (int argc, char** argv)
+{
+  int format = 0;
+  int n[3] = { 5, 5, 5 };
+  int dims = 3;
+  char* infile = 0;
+  char* vtffile = 0;
+
+  for (int i = 1; i < argc; i++)
+    if (!strcmp(argv[i],"-vtf") && i < argc-1)
+      format = atoi(argv[++i]);
+    else if (!strcmp(argv[i],"-nviz") && i < argc-1)
+      n[0] = n[1] = n[2] = atoi(argv[++i]);
+    else if (!strcmp(argv[i],"-1D"))
+      dims = 1;
+    else if (!strcmp(argv[i],"-2D"))
+      dims = 2;
+    else if (!infile)
+      infile = argv[i];
+    else if (!vtffile)
+      vtffile = argv[i];
+    else
+      std::cerr <<"  ** Unknown option ignored: "<< argv[i] << std::endl;
+
+  if (!infile) {
+    std::cout <<"usage: "<< argv[0]
+              <<" <inputfile> [<vtffile>] [-nviz <nviz>] [-1D|-2D]"<< std::endl;
+    return 0;
+  }
+  else if (!vtffile)
+    vtffile = infile;
+
+  std::cout <<"\n >>> IFEM HDF5 to VTF converter <<<"
+            <<"\n ==================================\n"
+            <<"\nInput file: "<< infile
+            <<"\nOutput file: "<< vtffile
+            <<"\nNumber of visualization points: "
+            << n[0] <<" "<< n[1] << " " << n[2] << std::endl;
+
+  HDF5Writer hdf(strtok(infile,"."),true);
+  XMLWriter xml(infile);
+  xml.readInfo();
+
+  int levels = xml.getLastTimeLevel();
+  if (levels > 0) SIMinput::msgLevel = 1;
+  std::cout <<"Reading "<< infile <<": Time levels = "<< levels << std::endl;
+
+  const std::vector<XMLWriter::Entry>& entry = xml.getEntries();
+  std::vector<XMLWriter::Entry>::const_iterator it;
+
+  ProcessList processlist;
+  for (it = entry.begin(); it != entry.end(); ++it)
+    processlist[it->patchfile].push_back(*it);
+
+  ProcessList::const_iterator pit = processlist.begin();
+  for (int j = 1; pit != processlist.end(); ++pit, ++j)
+  {
+    SIMbase* sim;
+    if (dims == 1)
+      sim = new SIM1D;
+    else if (dims == 2)
+      sim = new SIM2D;
+    else
+      sim = new SIM3D;
+
+    std::stringstream dummy;
+    dummy << "PATCHFILE " << pit->first;
+    sim->parse(const_cast<char*>(dummy.str().c_str()),dummy);
+    std::stringstream dummy2;
+    std::string foo(pit->first);
+    dummy2 << "NODEFILE " << replaceAll(foo,".g2",".gno");
+    sim->parse(const_cast<char*>(dummy2.str().c_str()),dummy2);
+
+    if (!sim->preprocess(std::vector<int>(),false))
+      return 1;
+
+    bool ok = true;
+    if (processlist.size() > 1) {
+      std::string temp(strtok(vtffile,"."));
+      std::stringstream str;
+      str <<"-"<< j;
+      temp.append(str.str());
+      ok = sim->writeGlv(temp.c_str(),n,format);
+    }
+    else
+      ok = sim->writeGlv(vtffile,n,format);
+
+    int block = 0;
+    for (int i = 0; i <= levels && ok; ++i) {
+      int k = 20;
+      if (levels > 0) std::cout <<"\nTime level "<< i << std::endl;
+      for (it = pit->second.begin(); it != pit->second.end() && ok; ++it) {
+        Vector vec;
+        std::cout <<"Reading \""<< it->name <<"\""<< std::endl;
+        ok = hdf.readField(i,it->name,vec,sim,it->components);
+        if (it->description == "displacement")
+          ok &= sim->writeGlvS(vec,n,i+1,block,i,1,NULL,10,it->components);
+        else if (it->components < 2)
+          ok &= sim->writeGlvS(vec,n,i+1,block,i,2,it->name.c_str(),k++,1);
+        else
+          ok &= sim->writeGlvV(vec,it->name.c_str(),n,i+1,block);
+      }
+      if (ok)
+        ok = sim->writeGlvStep(i+1,i);
+    }
+    delete sim;
+
+    if (!ok) return 2;
+  }
+
+  return 0;
+}
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f5880ea0..4fa38c54 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -129,12 +129,11 @@ IF(NOT WIN32)
 ENDIF(NOT WIN32)
 
 # Make the IFEM library
-
 FILE(GLOB_RECURSE IFEM_SRCS ${PROJECT_SOURCE_DIR}/src/*.[Cf]
                             ${PROJECT_SOURCE_DIR}/3rdparty/*.C)
 ADD_LIBRARY(IFEM ${IFEM_SRCS})
 
-# Make the Apps
+# Make some Apps
 FILE(GLOB_RECURSE Poisson_SRCS ${PROJECT_SOURCE_DIR}/Apps/Poisson/*.C)
 ADD_EXECUTABLE(Poisson ${Poisson_SRCS})
 TARGET_LINK_LIBRARIES(Poisson IFEM ${DEPLIBS})
@@ -145,7 +144,13 @@ FILE(GLOB_RECURSE LinEl_SRCS
 ADD_EXECUTABLE(LinEl ${LinEl_SRCS})
 TARGET_LINK_LIBRARIES(LinEl IFEM ${DEPLIBS})
 
-# Tests
+IF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
+  FILE(GLOB_RECURSE HDF2VTF_SRCS ${PROJECT_SOURCE_DIR}/Apps/HDF5toVTF/*.C)
+  ADD_EXECUTABLE(HDF5toVTF ${HDF2VTF_SRCS})
+  TARGET_LINK_LIBRARIES(HDF5toVTF IFEM ${DEPLIBS})
+ENDIF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
+
+# Regression tests
 FILE(GLOB_RECURSE LINEL_TESTFILES "${PROJECT_SOURCE_DIR}/Apps/LinearElasticity/Test/*.reg")
 FOREACH(TESTFILE ${LINEL_TESTFILES})
   ADD_TEST(${TESTFILE} ${PROJECT_SOURCE_DIR}/test/regtest.sh "${CMAKE_BINARY_DIR}/${EXECUTABLE_OUTPUT_PATH}/LinEl" "${TESTFILE}")
diff --git a/HOWTO b/HOWTO
index 4e69d667..0932d926 100644
--- a/HOWTO
+++ b/HOWTO
@@ -1,17 +1,17 @@
 Dette må gjøres for å kompilere denne versjonen:
 
+sudo apt-get install libboost-dev
 sudo apt-get install libblas-dev
 sudo apt-get install liblapack-dev
 sudo apt-get install libarpack2-dev
 sudo apt-get install libsuperlu3-dev
-sudo apt-get install libboost-dev
 sudo apt-get install libpetsc3.0.0-dev
 
 Installer GoToolsCore og GoTrivariate fra GoTools-svn.
 Gjør så cmake ., make og sudo make install i disse katalogene
 
 CMake og out-of-tree builds:
-CMake støtter debug og release-builds samtidig vi det som kalles
+CMake støtter debug og release-builds samtidig via det som kalles
 out-of-tree builds. In-tree builds av App'ene forventer at har byggefilene
 i en underkatalog med samme navn som bygge-typen. Feks for å bygge Debug
 gjør vi:
@@ -20,26 +20,28 @@ mkdir Debug
 cd Debug
 cmake .. -DCMAKE_BUILD_TYPE:STRING=Debug
 
-På samme måte kan vi lage en release-katalog. 
+På samme måte kan vi lage en release-katalog.
 
 Merk: Hvis du har en CMakeCache.txt i root når du prøver dette, vil
 det ikke fungere.
 
 Flagg av interesse:
 Per default lenker vi mot et minimum av bibliotek. Det betyr ingen
-PETSc, ingen SuperLU, ingen VTFWriter og ingen SAMG. Disse kan slås på
-med opsjoner:
+PETSc og ingen SAMG. Disse kan slås på med opsjoner:
 
--DENABLE_SUPERLU:BOOL=1, -DENABLE_PETSC:BOOL=1, -DENABLE_VTFWRITER:BOOL=1
-og -DENABLE_SAMG:BOOL=1.
+-DENABLE_PETSC:BOOL=1 og -DENABLE_SAMG:BOOL=1.
 
-Vi bygger kun libIFEM og Apps/Poisson.
-Stokes og FiniteDefElasticity har egen CMakeLists.txt. Disse er
-satt opp til å bruke in-tree kopi av libIFEM per default, men sjekker
-system hvis den ikke finner in-tree. Du kan tvinge system med 
--DFORCE_SYSTEM_IFEM:BOOL=1. Merk at in-tree sjekkes både for <root>/<type> og
-i <root>.
+For å kompilere mot parallel Petsc bruker du -DENABLE_PARALLEL_PETSC:BOOL=1
+
+Ved å spesifisere -DDISABLE_SUPERLU:BOOL=1 kan du slå av SuperLU ligningsløseren.
+Ved å spesifisere -DENABLE_SUPERLU_MT:BOOL=1 kan du aktivere multi-threaded SuperLU
+istedet for den serielle versjonen.
 
 Ved å spesifisere -DDISABLE_HDF5:BOOL=1 kan du slå av HDF5-støtten.
 
-For å kompilere mot parallel Petsc bruker du -DENABLE_PARALLEL_PETSC:BOOL=1
+Vi bygger kun libIFEM, Apps/Poisson og Apps/LinearElasticity.
+Stokes og FiniteDefElasticity har egne CMakeLists.txt. Disse
+er satt opp til å bruke in-tree kopi av libIFEM per default,
+men sjekker system hvis den ikke finner in-tree.
+Du kan tvinge system med -DFORCE_SYSTEM_IFEM:BOOL=1.
+Merk at in-tree sjekkes både for <root>/<type> og i <root>.
diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index 24488517..c2b69672 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -484,10 +484,10 @@ bool ASMbase::injectNodeVec (const Vector& nodeVec, Vector& globRes,
 {
   if (nndof == 0) nndof = nf;
 
-  if (nodeVec.size() < MLGN.size()*nndof)
+  if (nodeVec.size() != MLGN.size()*nndof)
   {
     std::cerr <<" *** ASMbase::injectNodeVec:: Invalid patch vector, size = "
-	      << nodeVec.size() <<" < "<< MLGN.size()*nndof << std::endl;
+	      << nodeVec.size() <<" != "<< MLGN.size()*nndof << std::endl;
     return false;
   }
 
diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index 4cde90ce..71e9c0f1 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -1036,17 +1036,19 @@ bool SIMbase::writeGlvV (const Vector& vec, const char* fieldName,
 
 bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 			 int iStep, int& nBlock, double time,
-			 bool psolOnly, const char* vecName)
+			 char psolOnly, const char* pvecName,
+			 int idBlock, int psolComps)
 {
   if (psol.empty())
     return true;
   else if (!myVtf)
     return false;
 
-  if (!psolOnly)
+  if (!psolOnly && myProblem)
     myProblem->initResultPoints(time);
 
   Matrix field;
+  Vector lovec;
   size_t i, j;
   int geomID = 0;
   const size_t pMAX = 6;
@@ -1070,11 +1072,14 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 
     // 1. Evaluate primary solution variables
 
-    myModel[i]->extractNodeVec(psol,myProblem->getSolution());
-    if (!myModel[i]->evalSolution(field,myProblem->getSolution(),nViz))
+    Vector& locvec = myProblem ? myProblem->getSolution() : lovec;
+    myModel[i]->extractNodeVec(psol,locvec,psolComps);
+    if (!myModel[i]->evalSolution(field,locvec,nViz))
       return false;
 
-    if (!myVtf->writeVres(field,++nBlock,++geomID,nVcomp))
+    if (psolOnly > 1)
+      geomID++;
+    else if (!myVtf->writeVres(field,++nBlock,++geomID,nVcomp))
       return false;
     else
       vID[0].push_back(nBlock);
@@ -1085,7 +1090,7 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
       else
 	dID[j].push_back(nBlock);
 
-    if (psolOnly) continue; // skip secondary solution
+    if (psolOnly || !myProblem) continue; // skip secondary solution
 
     // 2. Direct evaluation of secondary solution variables
 
@@ -1154,22 +1159,32 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
   // Write result block identifications
 
   bool ok = true;
-  int idBlock = 10;
-  if (vecName)
-    ok = myVtf->writeVblk(vID[0],vecName,idBlock,iStep);
-  else
-    ok = myVtf->writeDblk(vID[0],"Solution",idBlock,iStep);
-  if (!ok) return false;
+  if (!vID[0].empty())
+    if (pvecName)
+      ok = myVtf->writeVblk(vID[0],pvecName,idBlock,iStep);
+    else
+      ok = myVtf->writeDblk(vID[0],"Solution",idBlock,iStep);
 
-  if (scalarEq && !psolOnly)
-    if (!myVtf->writeVblk(vID[1],"Flux",idBlock+1,iStep))
-      return false;
+  if (ok && !vID[1].empty())
+    ok = myVtf->writeVblk(vID[1],"Flux",idBlock+1,iStep);
 
-  for (j = 0; j < pMAX && !dID[j].empty(); j++)
-    if (!myVtf->writeSblk(dID[j],myProblem->getField1Name(j),++idBlock,iStep))
-      return false;
+  std::string pname;
+  if (!myProblem)
+  {
+    pname = pvecName ? pvecName : "Solution";
+    if (psolOnly < 2) pname += "_w";
+  }
 
-  if (psolOnly) return true;
+  for (j = 0; j < pMAX && !dID[j].empty() && ok; j++)
+  {
+    if (myProblem)
+      pname = myProblem->getField1Name(j);
+    else if (psolOnly < 2)
+      (*pname.rbegin()) ++;
+    ok = myVtf->writeSblk(dID[j],pname.c_str(),++idBlock,iStep);
+  }
+
+  if (psolOnly || !myProblem || !ok) return ok;
 
   size_t nf = myProblem->getNoFields(2);
   for (i = j = 0; i < nf && !sID[j].empty(); i++, j++)
diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h
index 896fdc8d..61f79d3f 100644
--- a/src/SIM/SIMbase.h
+++ b/src/SIM/SIMbase.h
@@ -296,15 +296,18 @@ public:
   //! \param nBlock Running result block counter
   //! \param[in] time Load/time step parameter
   //! \param[in] psolOnly If \e true, skip secondary solution field evaluation
-  //! \param[in] vecName Optional name of the primary vector field solution
+  //! \param[in] pvecName Optional name of the primary vector field solution
+  //! \param[in] idBlock Starting value of result block numbering
+  //! \param[in] psolComps Optional number of primary solution components
   //!
   //! \details The primary vector field solution is written as a deformation
-  //! plot if \a vecName is NULL. Otherwise, it is written as a named vector
+  //! plot if \a pvecName is NULL. Otherwise, it is written as a named vector
   //! field. If an analytical solution is provided, the exact secondary solution
   //! fields are written to the VTF-file as well.
   virtual bool writeGlvS(const Vector& psol, const int* nViz, int iStep,
-			 int& nBlock, double time = 0.0, bool psolOnly = false,
-			 const char* vecName = 0);
+			 int& nBlock, double time = 0.0, char psolOnly = 0,
+			 const char* pvecName = 0, int idBlock = 10,
+			 int psolComps = 0);
 
   //! \brief Writes a mode shape and associated eigenvalue to the VTF-file.
   //! \details The eigenvalue is used as a label on the step state info
diff --git a/src/Utility/DataExporter.C b/src/Utility/DataExporter.C
index 639e90a6..07472b56 100644
--- a/src/Utility/DataExporter.C
+++ b/src/Utility/DataExporter.C
@@ -3,15 +3,26 @@
 #include "DataExporter.h"
 #include <iostream>
 #include <algorithm>
+#ifdef PARALLEL_PETSC
+#include <mpi.h>
+#endif
 
 
-DataExporter::DataExporter(bool dynWrts) : deleteWriters(dynWrts), m_level(-1)
+DataWriter::DataWriter (const std::string& name) : m_name(name)
 {
+#ifdef PARALLEL_PETSC
+  MPI_Comm_size(MPI_COMM_WORLD,&m_size);
+  MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
+#else
+  m_size = 1;
+  m_rank = 0;
+#endif
 }
 
-DataExporter::~DataExporter()
+
+DataExporter::~DataExporter ()
 {
-  if (deleteWriters)
+  if (m_delete)
     for (size_t i = 0; i < m_writers.size(); i++)
       delete m_writers[i];
 }
@@ -40,6 +51,7 @@ bool DataExporter::registerWriter(DataWriter* writer)
   return true;
 }
 
+
 bool DataExporter::setFieldValue(const std::string& name,
 				 void* data, void* data2)
 {
@@ -52,10 +64,11 @@ bool DataExporter::setFieldValue(const std::string& name,
   return true;
 }
 
+
 bool DataExporter::dumpTimeLevel()
 {
   if (m_level == -1)
-    m_level = getWritersTimeLevel()+1;
+    m_level = this->getWritersTimeLevel()+1;
 
   std::map<std::string,FileEntry>::iterator it;
   std::vector<DataWriter*>::iterator it2;
@@ -72,7 +85,7 @@ bool DataExporter::dumpTimeLevel()
           (*it2)->writeSIM(m_level,*it);
           break;
         default:
-	  std::cout <<"DataExporter: Invalid field type registered, skipping"
+	  std::cerr <<"DataExporter: Invalid field type registered, skipping"
 		    << std::endl;
           break;
       }
@@ -84,52 +97,59 @@ bool DataExporter::dumpTimeLevel()
   return true;
 }
 
-bool DataExporter::loadTimeLevel(int level, DataWriter* input)
-{
-  if (!input && m_writers.empty())
-    return false;
-  if (!input)
-    input = m_writers.front();
-  if (level == -1)
-    level = m_level = input->getLastTimeLevel();
-  if (level == -1)
-    return false;
 
+bool DataExporter::loadTimeLevel (int level, DataWriter* input)
+{
+  if (!input)
+    if (m_writers.empty())
+      return false;
+    else
+      input = m_writers.front();
+
+  if (level == -1)
+    if ((m_level = input->getLastTimeLevel()) < 0)
+      return false;
+    else
+      level = m_level;
+
+  bool ok = true;
   input->openFile(level);
   std::map<std::string,FileEntry>::iterator it;
-  for (it = m_entry.begin(); it != m_entry.end(); ++it)
-  {
+  for (it = m_entry.begin(); it != m_entry.end() && ok; ++it)
     if (!it->second.data)
-      return false;
-    switch (it->second.field)
+      ok = false;
+    else switch (it->second.field)
     {
       case VECTOR:
-        input->readVector(level,*it);
+        ok = input->readVector(level,*it);
         break;
       case SIM:
-        input->readSIM(level,*it);
+        ok = input->readSIM(level,*it);
         break;
       default:
-	std::cout << "DataExporter::loadTimeLevel: Invalid field type "
-		  << "registered, skipping" << std::endl;
+        ok = false;
+	std::cerr <<" *** DataExporter: Invalid field type registered "
+		  << it->second.field << std::endl;
         break;
     }
-  }
+
   input->closeFile(level);
   m_level++;
 
-  return true;
+  return ok;
 }
 
-int DataExporter::getTimeLevel()
+
+int DataExporter::getTimeLevel ()
 {
   if (m_level == -1)
-    m_level = getWritersTimeLevel();
+    m_level = this->getWritersTimeLevel();
 
   return m_level;
 }
 
-int DataExporter::getWritersTimeLevel()
+
+int DataExporter::getWritersTimeLevel () const
 {
   std::vector<int> levels;
   std::vector<DataWriter*>::const_iterator it2;
diff --git a/src/Utility/DataExporter.h b/src/Utility/DataExporter.h
index c457d81f..c6b48131 100644
--- a/src/Utility/DataExporter.h
+++ b/src/Utility/DataExporter.h
@@ -9,51 +9,55 @@
 
 class DataWriter;
 
-class DataExporter {
-  public:
-    enum FieldType {
-      VECTOR,
-      SIM
-    };
+class DataExporter
+{
+ public:
+  enum FieldType {
+    VECTOR,
+    SIM
+  };
 
-    struct FileEntry {
-      std::string description;
-      FieldType field;
-      int size;
-      void* data;
-      void* data2;
-    };
+  struct FileEntry {
+    std::string description;
+    FieldType field;
+    int size;
+    void* data;
+    void* data2;
+  };
 
-    DataExporter(bool dynamicWriters = false);
-    virtual ~DataExporter();
+  //! \brief Default constructor.
+  //! \param[in] dynWriters If \e true, delete the writers on destruction.
+  DataExporter(bool dynWriters = false) : m_delete(dynWriters), m_level(-1) {}
+  //! \brief The desctructor deletes the writers if \a dynWriters was \e true.
+  ~DataExporter();
 
-    //! \brief Registers an entry for storage.
-    //! param[in] name Name of entry
-    //! param[in] description Description of entry
-    //! param[in] Type of entry
-    //! param[in] size  set to number of entries in an array,
-    //                 the time level to use for SIM
-    bool registerField(const std::string& name,
-                       const std::string& description,
-                       FieldType field, size_t size=0);
+  //! \brief Registers an entry for storage.
+  //! param[in] name Name of entry
+  //! param[in] description Description of entry
+  //! param[in] field Type of entry
+  //! param[in] size Number of entries in an array,
+  //! the time level to use for SIM
+  bool registerField(const std::string& name,
+		     const std::string& description,
+		     FieldType field, size_t size=0);
 
-    bool registerWriter(DataWriter* writer);
+  bool registerWriter(DataWriter* writer);
 
-    bool setFieldValue(const std::string& name, void* data, void* data2=NULL);
+  bool setFieldValue(const std::string& name, void* data, void* data2=NULL);
 
-    bool dumpTimeLevel();
+  bool dumpTimeLevel();
 
-    //! \brief Loads last time level with first registered writer by default.
-    bool loadTimeLevel(int level=-1, DataWriter* input=NULL);
-    int getTimeLevel();
+  //! \brief Loads last time level with first registered writer by default.
+  bool loadTimeLevel(int level=-1, DataWriter* input=NULL);
+  int getTimeLevel();
 
-  protected:
-    int getWritersTimeLevel();
+protected:
+  int getWritersTimeLevel() const;
 
-    std::map<std::string,FileEntry> m_entry;
-    std::vector<DataWriter*> m_writers;
-    bool deleteWriters;
-    int m_level;
+  std::map<std::string,FileEntry> m_entry;
+  std::vector<DataWriter*> m_writers;
+  bool m_delete;
+  int m_level;
 };
 
 
@@ -62,7 +66,7 @@ typedef std::pair<std::string,DataExporter::FileEntry> DataEntry;
 class DataWriter
 {
 protected:
-  DataWriter() {}
+  DataWriter(const std::string& name);
 
 public:
   virtual ~DataWriter() {}
@@ -73,7 +77,13 @@ public:
   virtual void closeFile(int level) = 0;
 
   virtual void writeVector(int level, const DataEntry& entry) = 0;
-  virtual void readVector(int level, const DataEntry& entry) = 0;
+  virtual bool readVector(int level, const DataEntry& entry) = 0;
   virtual void writeSIM(int level, const DataEntry& entry) = 0;
-  virtual void readSIM(int level, const DataEntry& entry) = 0;
+  virtual bool readSIM(int level, const DataEntry& entry) = 0;
+
+protected:
+  std::string m_name; //!< File name
+
+  int m_size; //!< Number of MPI nodes (processors)
+  int m_rank; //!< MPI rank (processor ID)
 };
diff --git a/src/Utility/HDF5Writer.C b/src/Utility/HDF5Writer.C
index 8bb936f1..9f99d218 100644
--- a/src/Utility/HDF5Writer.C
+++ b/src/Utility/HDF5Writer.C
@@ -4,7 +4,6 @@
 #include "SIMbase.h"
 #include "IntegrandBase.h"
 #include <sstream>
-#include <numeric>
 #include <sys/stat.h>
 
 #ifdef HAS_HDF5
@@ -16,30 +15,23 @@
 #endif
 
 
-HDF5Writer::HDF5Writer (const std::string& name) : m_hdf5(name+".hdf5")
+HDF5Writer::HDF5Writer (const std::string& name, bool append)
+  : DataWriter(name+".hdf5")
 {
 #ifdef HAS_HDF5
   struct stat temp;
   // file already exists - open and find the next group
-  if (stat(m_hdf5.c_str(),&temp) == 0)
+  if (append && stat(m_name.c_str(),&temp) == 0)
     m_flag = H5F_ACC_RDWR;
   else
     m_flag = H5F_ACC_TRUNC;
 #endif
-#ifdef PARALLEL_PETSC
-  MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
-  MPI_Comm_size(MPI_COMM_WORLD,&m_size);
-#else
-  m_rank = 0;
-#endif
 }
 
-HDF5Writer::~HDF5Writer()
-{
-}
 
 int HDF5Writer::getLastTimeLevel()
 {
+  int result = 0;
 #ifdef HAS_HDF5
   if (m_flag == H5F_ACC_TRUNC)
     return -1;
@@ -47,12 +39,11 @@ int HDF5Writer::getLastTimeLevel()
   hid_t acc_tpl = H5P_DEFAULT;
 #ifdef PARALLEL_PETSC
   MPI_Info info = MPI_INFO_NULL;
-  acc_tpl = H5Pcreate (H5P_FILE_ACCESS);
+  acc_tpl = H5Pcreate(H5P_FILE_ACCESS);
   H5Pset_fapl_mpio(acc_tpl, MPI_COMM_WORLD, info);
 #endif
 
-  m_file = H5Fopen(m_hdf5.c_str(),m_flag,acc_tpl);
-  int result=0;
+  m_file = H5Fopen(m_name.c_str(),m_flag,acc_tpl);
   while (1) {
     std::stringstream str;
     str << '/' << result;
@@ -61,11 +52,8 @@ int HDF5Writer::getLastTimeLevel()
     result++;
   }
   H5Fclose(m_file);
-
-  return result-1;
-#else
-  return -1;
 #endif
+  return result-1;
 }
 
 void HDF5Writer::openFile(int level)
@@ -78,9 +66,9 @@ void HDF5Writer::openFile(int level)
   H5Pset_fapl_mpio(acc_tpl, MPI_COMM_WORLD, info);
 #endif
   if (m_flag == H5F_ACC_TRUNC)
-    m_file = H5Fcreate(m_hdf5.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
+    m_file = H5Fcreate(m_name.c_str(),m_flag,H5P_DEFAULT,acc_tpl);
   else
-    m_file = H5Fopen(m_hdf5.c_str(),m_flag,acc_tpl);
+    m_file = H5Fopen(m_name.c_str(),m_flag,acc_tpl);
 
   std::stringstream str;
   str << '/' << level;
@@ -153,9 +141,10 @@ void HDF5Writer::writeArray(int group, const std::string& name,
 #endif
 }
 
-void HDF5Writer::readVector(int level, const DataEntry& entry)
+bool HDF5Writer::readVector(int level, const DataEntry& entry)
 {
 //  readArray(level,entry.first,entry.second.size,entry.second.data);
+  return true;
 }
 
 void HDF5Writer::writeVector(int level, const DataEntry& entry)
@@ -163,14 +152,15 @@ void HDF5Writer::writeVector(int level, const DataEntry& entry)
   writeArray(level,entry.first,entry.second.size,entry.second.data);
 }
 
-void HDF5Writer::readSIM(int level, const DataEntry& entry)
+bool HDF5Writer::readSIM(int level, const DataEntry& entry)
 {
+  bool ok = true;
 #ifdef HAS_HDF5
   SIMbase* sim = static_cast<SIMbase*>(entry.second.data);
   Vector* sol = static_cast<Vector*>(entry.second.data2);
-  if (!sol) return;
+  if (!sol) return false;
 
-  for (int i = 0; i < sim->getNoPatches(); ++i) {
+  for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
     std::stringstream str;
     str << level;
     str << '/';
@@ -180,12 +170,42 @@ void HDF5Writer::readSIM(int level, const DataEntry& entry)
     if (loc > 0) {
       double* tmp = NULL; int siz = 0;
       readArray(group2,entry.first,siz,tmp);
-      sim->injectPatchSolution(*sol,Vector(tmp,siz),loc-1,entry.second.size);
+      ok = sim->injectPatchSolution(*sol,Vector(tmp,siz),
+                                    loc-1,entry.second.size);
       delete[] tmp;
     }
     H5Gclose(group2);
   }
 #endif
+  return ok;
+}
+
+bool HDF5Writer::readField(int level, const std::string& name,
+                           Vector& vec, SIMbase* sim, int components)
+{
+  bool ok = true;
+  openFile(level);
+#ifdef HAS_HDF5
+  vec.resize(sim->getNoNodes()*components);
+  for (int i = 0; i < sim->getNoPatches() && ok; ++i) {
+    std::stringstream str;
+    str << level;
+    str << '/';
+    str << i+1;
+    hid_t group2 = H5Gopen2(m_file,str.str().c_str(),H5P_DEFAULT);
+    int loc = sim->getLocalPatchIndex(i+1);
+    if (loc > 0) {
+      double* tmp = NULL; int siz = 0;
+      readArray(group2,name,siz,tmp);
+      ok = sim->injectPatchSolution(vec,Vector(tmp,siz),
+                                    loc-1,components);
+      delete[] tmp;
+    }
+    H5Gclose(group2);
+  }
+#endif
+  closeFile(level);
+  return ok;
 }
 
 void HDF5Writer::writeSIM(int level, const DataEntry& entry)
@@ -208,17 +228,29 @@ void HDF5Writer::writeSIM(int level, const DataEntry& entry)
       group2 = H5Gcreate2(m_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT);
     int loc = sim->getLocalPatchIndex(i+1);
     if (loc > 0) // we own the patch
+    {
       sim->extractPatchSolution(*sol,loc-1);
-    if (entry.second.size == -1) {
-      Matrix field;
-      if (loc > 0)
-	sim->evalSecondarySolution(field,loc-1);
-      for (size_t j = 0; j < prob->getNoFields(2); ++j)
-        writeArray(group2,prob->getField2Name(j),field.cols(),
-		   field.getRow(j+1).ptr());
+      Vector& psol = const_cast<Integrand*>(prob)->getSolution();
+      writeArray(group2,entry.first,psol.size(),psol.ptr());
+      // TODO: For mixed methods we need to output each primary field
+      // TODO: The above is sufficient for simulation restarts.
+      // TODO: For mixed methods we need to output each primary field
+      // TODO: separately in addition, with reference to correct spline basis.
+      if (entry.second.size == -1) {
+        Matrix field;
+        sim->evalSecondarySolution(field,loc-1);
+        for (size_t j = 0; j < prob->getNoFields(2); ++j)
+          writeArray(group2,prob->getField2Name(j),field.cols(),
+                     field.getRow(j+1).ptr());
+      }
+    }
+    else // must write empty dummy records for the other patches
+    {
+      double dummy;
+      writeArray(group2,entry.first,0,&dummy);
+      for (size_t j = 0; j < prob->getNoFields(2); ++j)
+        writeArray(group2,prob->getField2Name(j),0,&dummy);
     }
-    Vector& psol = const_cast<Integrand*>(prob)->getSolution();
-    writeArray(group2, entry.first, loc > 0 ? psol.size() : 0, psol.ptr());
     H5Gclose(group2);
   }
 #endif
@@ -226,15 +258,15 @@ void HDF5Writer::writeSIM(int level, const DataEntry& entry)
 
 bool HDF5Writer::checkGroupExistence(int parent, const char* path)
 {
+  bool result = false;
 #ifdef HAS_HDF5
   // turn off errors to avoid cout spew
   herr_t (*old_func)(void*);
   void* old_client_data;
   H5Eget_auto(&old_func,&old_client_data);
   H5Eset_auto(NULL,NULL);
-  bool result =H5Gget_objinfo((hid_t)parent,path,0,NULL) == 0;
+  result = H5Gget_objinfo((hid_t)parent,path,0,NULL) == 0;
   H5Eset_auto(old_func,old_client_data);
-  return result;
 #endif
-  return false;
+  return result;
 }
diff --git a/src/Utility/HDF5Writer.h b/src/Utility/HDF5Writer.h
index bc39b7c6..03b2448a 100644
--- a/src/Utility/HDF5Writer.h
+++ b/src/Utility/HDF5Writer.h
@@ -3,35 +3,37 @@
 #pragma once
 
 #include "DataExporter.h"
+#include "MatVec.h"
+
+class SIMbase;
 
 
-class HDF5Writer : public DataWriter {
-  public:
+class HDF5Writer : public DataWriter
+{
+public:
+  HDF5Writer(const std::string& name, bool append = false);
+  virtual ~HDF5Writer() {}
 
-    HDF5Writer(const std::string& name);
-    virtual ~HDF5Writer();
+  virtual int getLastTimeLevel();
 
-    int getLastTimeLevel();
+  virtual void openFile(int level);
+  virtual void closeFile(int level);
 
-    void openFile(int level);
-    void closeFile(int level);
+  virtual void writeVector(int level, const DataEntry& entry);
+  virtual bool readVector(int level, const DataEntry& entry);
+  virtual void writeSIM(int level, const DataEntry& entry);
+  virtual bool readSIM(int level, const DataEntry& entry);
 
-    void writeVector(int level, const DataEntry& entry);
-    void readVector(int level, const DataEntry& entry);
-    void writeSIM(int level, const DataEntry& entry);
-    void readSIM(int level, const DataEntry& entry);
+  bool readField(int level, const std::string& name,
+                 Vector& vec, SIMbase* sim, int components);
 
-  protected:
-    void writeArray(int group, const std::string& name,
-                    int len, void* data);
-    void readArray(int group, const std::string& name,
-                   int& len, double*& data);
-    bool checkGroupExistence(int parent, const char* group);
+protected:
+  void writeArray(int group, const std::string& name,
+                  int len, void* data);
+  void readArray(int group, const std::string& name,
+                 int& len, double*& data);
+  bool checkGroupExistence(int parent, const char* group);
 
-    std::string m_hdf5;
-    int m_file;
-    unsigned int m_flag;
-
-    int m_rank; // MPI rank
-    int m_size; // number of MPI nodes
+  int          m_file;
+  unsigned int m_flag;
 };
diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h
index f12b3712..369bc0ee 100644
--- a/src/Utility/Utilities.h
+++ b/src/Utility/Utilities.h
@@ -23,10 +23,10 @@
 namespace utl
 {
   //! \brief Parses a character string into an integer or an integer range.
-  //! \param values The integer value(s) is(are) appended to this vector
-  //! \param argv Character string with integer data
+  //! \param values The integer value(s) is/are appended to this vector
+  //! \param[in] argv Character string with integer data
   //!
-  //! \details An integer range is recognised through the syntax <i>:<j>.
+  //! \details An integer range is recognised through the syntax \a i:j.
   void parseIntegers(std::vector<int>& values, const char* argv);
 
   //! \brief Reads one line, ignoring comment lines and leading blanks.
diff --git a/src/Utility/XMLWriter.C b/src/Utility/XMLWriter.C
index 5ad1f512..7a6e4857 100644
--- a/src/Utility/XMLWriter.C
+++ b/src/Utility/XMLWriter.C
@@ -3,31 +3,23 @@
 #include "XMLWriter.h"
 #include "SIMbase.h"
 #include "IntegrandBase.h"
-
-#ifdef PARALLEL_PETSC
-#include <mpi.h>
-#endif
+#include "StringUtils.h"
+#include "tinyxml.h"
+#include <fstream>
+#include <cstdio>
 
 
-XMLWriter::XMLWriter(const std::string& name) : m_xml(name+".xml")
+XMLWriter::XMLWriter(const std::string& name) : DataWriter(name+".xml")
 {
   m_doc = NULL;
-#ifdef PARALLEL_PETSC
-  MPI_Comm_rank(MPI_COMM_WORLD,&m_rank);
-  MPI_Comm_size(MPI_COMM_WORLD,&m_size);
-#else
-  m_rank = 0;
-#endif
+  m_node = NULL;
 }
 
-XMLWriter::~XMLWriter()
-{
-}
 
 int XMLWriter::getLastTimeLevel()
 {
   int result=-1;
-  TiXmlDocument doc(m_xml.c_str());
+  TiXmlDocument doc(m_name.c_str());
   doc.LoadFile();
   TiXmlHandle handle(&doc);
   TiXmlElement* levels = handle.FirstChild("info").
@@ -38,10 +30,12 @@ int XMLWriter::getLastTimeLevel()
   return result;
 }
 
+
 void XMLWriter::openFile(int level)
 {
   if (m_rank != 0)
     return;
+
   m_doc = new TiXmlDocument;
   TiXmlElement element("info");
   m_node = m_doc->InsertEndChild(element);
@@ -59,13 +53,37 @@ void XMLWriter::closeFile(int level)
   TiXmlText value(temp);
   pNewNode->InsertEndChild(value);
 
-  m_doc->SaveFile(m_xml);
+  m_doc->SaveFile(m_name);
   delete m_doc;
   m_doc = NULL;
 }
 
-void XMLWriter::readVector(int level, const DataEntry& entry)
+
+void XMLWriter::readInfo()
 {
+  TiXmlDocument doc(m_name);
+  doc.LoadFile();
+  TiXmlHandle handle(&doc);
+  TiXmlElement* elem = handle.FirstChild("info").
+                                         FirstChild("entry").ToElement();
+  while (elem) {
+    if (strcasecmp(elem->Attribute("type"),"field") == 0) {
+      Entry entry;
+      entry.name = elem->Attribute("name");
+      entry.description = elem->Attribute("description");
+      entry.patches = atoi(elem->Attribute("patches"));
+      entry.components = atoi(elem->Attribute("components"));
+      entry.patchfile = elem->Attribute("patchfile");
+      m_entry.push_back(entry);
+    }
+    elem = elem->NextSiblingElement("entry");
+  }
+}
+
+
+bool XMLWriter::readVector(int level, const DataEntry& entry)
+{
+  return true;
 }
 
 void XMLWriter::writeVector(int level, const DataEntry& entry)
@@ -81,8 +99,10 @@ void XMLWriter::writeVector(int level, const DataEntry& entry)
   m_node->InsertEndChild(element);
 }
 
-void XMLWriter::readSIM(int level, const DataEntry& entry)
+
+bool XMLWriter::readSIM(int level, const DataEntry& entry)
 {
+  return true;
 }
 
 void XMLWriter::writeSIM(int level, const DataEntry& entry)
@@ -90,25 +110,38 @@ void XMLWriter::writeSIM(int level, const DataEntry& entry)
   if (m_rank != 0)
     return;
 
+  // Assume that all fields use the same basis as the geometry.
+  // TODO: Not true for mixed methods, fix later...
+  std::string g2file(m_name);
+  std::ofstream os(replaceAll(g2file,".xml",".g2").c_str());
+  static_cast<SIMbase*>(entry.second.data)->dumpGeometry(os);
+
   int nPatch = static_cast<SIMbase*>(entry.second.data)->getNoPatches();
   const Integrand* prb = static_cast<SIMbase*>(entry.second.data)->getProblem();
 
   // primary solution
-  addField(entry.first,entry.second.description,"field",nPatch);
+  addField(entry.first,entry.second.description,"field",g2file,
+           prb->getNoFields(1),nPatch);
 
   // secondary solutions
   if (entry.second.size == -1)
     for (size_t j = 0; j < prb->getNoFields(2); j++)
-      addField(prb->getField2Name(j),"secondary","field",nPatch);
+      addField(prb->getField2Name(j),"secondary","field",g2file,1,nPatch);
 }
 
-void XMLWriter::addField(const std::string& name, const std::string& description,
-                         const std::string& type, int patches)
+
+void XMLWriter::addField (const std::string& name,
+                          const std::string& description,
+                          const std::string& type,
+                          const std::string& patchfile,
+                          int components, int patches)
 {
   TiXmlElement element("entry");
   element.SetAttribute("name",name.c_str());
   element.SetAttribute("description",description.c_str());
   element.SetAttribute("type",type.c_str());
   element.SetAttribute("patches",patches);
+  element.SetAttribute("components",components);
+  element.SetAttribute("patchfile",patchfile.c_str());
   m_node->InsertEndChild(element);
 }
diff --git a/src/Utility/XMLWriter.h b/src/Utility/XMLWriter.h
index 5a621002..2a92f1c8 100644
--- a/src/Utility/XMLWriter.h
+++ b/src/Utility/XMLWriter.h
@@ -3,32 +3,45 @@
 #pragma once
 
 #include "DataExporter.h"
-#include "tinyxml.h"
+
+class TiXmlDocument;
+class TiXmlNode;
 
 
-class XMLWriter : public DataWriter {
-  public:
-    XMLWriter(const std::string& name);
-    virtual ~XMLWriter();
+class XMLWriter : public DataWriter
+{
+public:
+  struct Entry {
+    std::string name;
+    std::string description;
+    std::string patchfile;
+    int patches;
+    int components;
+  };
 
-    int getLastTimeLevel();
+  XMLWriter(const std::string& name);
+  virtual ~XMLWriter() {}
 
-    void openFile(int level);
-    void closeFile(int level);
+  virtual int getLastTimeLevel();
 
-    void writeVector(int level, const DataEntry& entry);
-    void readVector(int level, const DataEntry& entry);
-    void writeSIM(int level, const DataEntry& entry);
-    void readSIM(int level, const DataEntry& entry);
+  void readInfo();
+  const std::vector<Entry>& getEntries() const { return m_entry; }
 
-  protected:
-    void addField(const std::string& name, const std::string& description,
-                  const std::string& type, int patches);
-    std::string m_xml;
+  virtual void openFile(int level);
+  virtual void closeFile(int level);
 
-    int m_rank; // MPI rank
-    int m_size; // number of MPI nodes
+  virtual void writeVector(int level, const DataEntry& entry);
+  virtual bool readVector(int level, const DataEntry& entry);
+  virtual void writeSIM(int level, const DataEntry& entry);
+  virtual bool readSIM(int level, const DataEntry& entry);
 
-    TiXmlDocument* m_doc;
-    TiXmlNode* m_node;
+protected:
+  void addField(const std::string& name, const std::string& description,
+                const std::string& type, const std::string& patchfile,
+                int components, int patches);
+
+  std::vector<Entry> m_entry;
+
+  TiXmlDocument* m_doc;
+  TiXmlNode* m_node;
 };