Adding HDF5 support (credit Arne Morten) + some modifications related to VTF output

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@932 e10b68d5-8a6e-419e-a041-bce267b0401d

3rdparty/tinyxml/tinystr.C

@@ -0,0 +1,116 @@
+/*
+www.sourceforge.net/projects/tinyxml
+Original file by Yves Berquin.
+
+This software is provided 'as-is', without any express or implied
+warranty. In no event will the authors be held liable for any
+damages arising from the use of this software.
+
+Permission is granted to anyone to use this software for any
+purpose, including commercial applications, and to alter it and
+redistribute it freely, subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must
+not claim that you wrote the original software. If you use this
+software in a product, an acknowledgment in the product documentation
+would be appreciated but is not required.
+
+2. Altered source versions must be plainly marked as such, and
+must not be misrepresented as being the original software.
+
+3. This notice may not be removed or altered from any source
+distribution.
+*/
+
+/*
+ * THIS FILE WAS ALTERED BY Tyge Løvset, 7. April 2005.
+ */
+
+
+#ifndef TIXML_USE_STL
+
+#include "tinystr.h"
+
+// Error value for find primitive
+const TiXmlString::size_type TiXmlString::npos = static_cast< TiXmlString::size_type >(-1);
+
+
+// Null rep.
+TiXmlString::Rep TiXmlString::nullrep_ = { 0, 0, { '\0' } };
+
+
+void TiXmlString::reserve (size_type cap)
+{
+	if (cap > capacity())
+	{
+		TiXmlString tmp;
+		tmp.init(length(), cap);
+		memcpy(tmp.start(), data(), length());
+		swap(tmp);
+	}
+}
+
+
+TiXmlString& TiXmlString::assign(const char* str, size_type len)
+{
+	size_type cap = capacity();
+	if (len > cap || cap > 3*(len + 8))
+	{
+		TiXmlString tmp;
+		tmp.init(len);
+		memcpy(tmp.start(), str, len);
+		swap(tmp);
+	}
+	else
+	{
+		memmove(start(), str, len);
+		set_size(len);
+	}
+	return *this;
+}
+
+
+TiXmlString& TiXmlString::append(const char* str, size_type len)
+{
+	size_type newsize = length() + len;
+	if (newsize > capacity())
+	{
+		reserve (newsize + capacity());
+	}
+	memmove(finish(), str, len);
+	set_size(newsize);
+	return *this;
+}
+
+
+TiXmlString operator + (const TiXmlString & a, const TiXmlString & b)
+{
+	TiXmlString tmp;
+	tmp.reserve(a.length() + b.length());
+	tmp += a;
+	tmp += b;
+	return tmp;
+}
+
+TiXmlString operator + (const TiXmlString & a, const char* b)
+{
+	TiXmlString tmp;
+	TiXmlString::size_type b_len = static_cast<TiXmlString::size_type>( strlen(b) );
+	tmp.reserve(a.length() + b_len);
+	tmp += a;
+	tmp.append(b, b_len);
+	return tmp;
+}
+
+TiXmlString operator + (const char* a, const TiXmlString & b)
+{
+	TiXmlString tmp;
+	TiXmlString::size_type a_len = static_cast<TiXmlString::size_type>( strlen(a) );
+	tmp.reserve(a_len + b.length());
+	tmp.append(a, a_len);
+	tmp += b;
+	return tmp;
+}
+
+
+#endif	// TIXML_USE_STL

3rdparty/tinyxml/tinystr.h

@@ -0,0 +1,319 @@
+/*
+www.sourceforge.net/projects/tinyxml
+Original file by Yves Berquin.
+
+This software is provided 'as-is', without any express or implied
+warranty. In no event will the authors be held liable for any
+damages arising from the use of this software.
+
+Permission is granted to anyone to use this software for any
+purpose, including commercial applications, and to alter it and
+redistribute it freely, subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must
+not claim that you wrote the original software. If you use this
+software in a product, an acknowledgment in the product documentation
+would be appreciated but is not required.
+
+2. Altered source versions must be plainly marked as such, and
+must not be misrepresented as being the original software.
+
+3. This notice may not be removed or altered from any source
+distribution.
+*/
+
+/*
+ * THIS FILE WAS ALTERED BY Tyge Lovset, 7. April 2005.
+ *
+ * - completely rewritten. compact, clean, and fast implementation.
+ * - sizeof(TiXmlString) = pointer size (4 bytes on 32-bit systems)
+ * - fixed reserve() to work as per specification.
+ * - fixed buggy compares operator==(), operator<(), and operator>()
+ * - fixed operator+=() to take a const ref argument, following spec.
+ * - added "copy" constructor with length, and most compare operators.
+ * - added swap(), clear(), size(), capacity(), operator+().
+ */
+
+#ifndef TIXML_USE_STL
+
+#ifndef TIXML_STRING_INCLUDED
+#define TIXML_STRING_INCLUDED
+
+#include <assert.h>
+#include <string.h>
+
+/*	The support for explicit isn't that universal, and it isn't really
+	required - it is used to check that the TiXmlString class isn't incorrectly
+	used. Be nice to old compilers and macro it here:
+*/
+#if defined(_MSC_VER) && (_MSC_VER >= 1200 )
+	// Microsoft visual studio, version 6 and higher.
+	#define TIXML_EXPLICIT explicit
+#elif defined(__GNUC__) && (__GNUC__ >= 3 )
+	// GCC version 3 and higher.s
+	#define TIXML_EXPLICIT explicit
+#else
+	#define TIXML_EXPLICIT
+#endif
+
+
+/*
+   TiXmlString is an emulation of a subset of the std::string template.
+   Its purpose is to allow compiling TinyXML on compilers with no or poor STL support.
+   Only the member functions relevant to the TinyXML project have been implemented.
+   The buffer allocation is made by a simplistic power of 2 like mechanism : if we increase
+   a string and there's no more room, we allocate a buffer twice as big as we need.
+*/
+class TiXmlString
+{
+  public :
+	// The size type used
+  	typedef size_t size_type;
+
+	// Error value for find primitive
+	static const size_type npos; // = -1;
+
+
+	// TiXmlString empty constructor
+	TiXmlString () : rep_(&nullrep_)
+	{
+	}
+
+	// TiXmlString copy constructor
+	TiXmlString ( const TiXmlString & copy) : rep_(0)
+	{
+		init(copy.length());
+		memcpy(start(), copy.data(), length());
+	}
+
+	// TiXmlString constructor, based on a string
+	TIXML_EXPLICIT TiXmlString ( const char * copy) : rep_(0)
+	{
+		init( static_cast<size_type>( strlen(copy) ));
+		memcpy(start(), copy, length());
+	}
+
+	// TiXmlString constructor, based on a string
+	TIXML_EXPLICIT TiXmlString ( const char * str, size_type len) : rep_(0)
+	{
+		init(len);
+		memcpy(start(), str, len);
+	}
+
+	// TiXmlString destructor
+	~TiXmlString ()
+	{
+		quit();
+	}
+
+	// = operator
+	TiXmlString& operator = (const char * copy)
+	{
+		return assign( copy, (size_type)strlen(copy));
+	}
+
+	// = operator
+	TiXmlString& operator = (const TiXmlString & copy)
+	{
+		return assign(copy.start(), copy.length());
+	}
+
+
+	// += operator. Maps to append
+	TiXmlString& operator += (const char * suffix)
+	{
+		return append(suffix, static_cast<size_type>( strlen(suffix) ));
+	}
+
+	// += operator. Maps to append
+	TiXmlString& operator += (char single)
+	{
+		return append(&single, 1);
+	}
+
+	// += operator. Maps to append
+	TiXmlString& operator += (const TiXmlString & suffix)
+	{
+		return append(suffix.data(), suffix.length());
+	}
+
+
+	// Convert a TiXmlString into a null-terminated char *
+	const char * c_str () const { return rep_->str; }
+
+	// Convert a TiXmlString into a char * (need not be null terminated).
+	const char * data () const { return rep_->str; }
+
+	// Return the length of a TiXmlString
+	size_type length () const { return rep_->size; }
+
+	// Alias for length()
+	size_type size () const { return rep_->size; }
+
+	// Checks if a TiXmlString is empty
+	bool empty () const { return rep_->size == 0; }
+
+	// Return capacity of string
+	size_type capacity () const { return rep_->capacity; }
+
+
+	// single char extraction
+	const char& at (size_type index) const
+	{
+		assert( index < length() );
+		return rep_->str[ index ];
+	}
+
+	// [] operator
+	char& operator [] (size_type index) const
+	{
+		assert( index < length() );
+		return rep_->str[ index ];
+	}
+
+	// find a char in a string. Return TiXmlString::npos if not found
+	size_type find (char lookup) const
+	{
+		return find(lookup, 0);
+	}
+
+	// find a char in a string from an offset. Return TiXmlString::npos if not found
+	size_type find (char tofind, size_type offset) const
+	{
+		if (offset >= length()) return npos;
+
+		for (const char* p = c_str() + offset; *p != '\0'; ++p)
+		{
+		   if (*p == tofind) return static_cast< size_type >( p - c_str() );
+		}
+		return npos;
+	}
+
+	void clear ()
+	{
+		//Lee:
+		//The original was just too strange, though correct:
+		//	TiXmlString().swap(*this);
+		//Instead use the quit & re-init:
+		quit();
+		init(0,0);
+	}
+
+	/*	Function to reserve a big amount of data when we know we'll need it. Be aware that this
+		function DOES NOT clear the content of the TiXmlString if any exists.
+	*/
+	void reserve (size_type cap);
+
+	TiXmlString& assign (const char* str, size_type len);
+
+	TiXmlString& append (const char* str, size_type len);
+
+	void swap (TiXmlString& other)
+	{
+		Rep* r = rep_;
+		rep_ = other.rep_;
+		other.rep_ = r;
+	}
+
+  private:
+
+	void init(size_type sz) { init(sz, sz); }
+	void set_size(size_type sz) { rep_->str[ rep_->size = sz ] = '\0'; }
+	char* start() const { return rep_->str; }
+	char* finish() const { return rep_->str + rep_->size; }
+
+	struct Rep
+	{
+		size_type size, capacity;
+		char str[1];
+	};
+
+	void init(size_type sz, size_type cap)
+	{
+		if (cap)
+		{
+			// Lee: the original form:
+			//	rep_ = static_cast<Rep*>(operator new(sizeof(Rep) + cap));
+			// doesn't work in some cases of new being overloaded. Switching
+			// to the normal allocation, although use an 'int' for systems
+			// that are overly picky about structure alignment.
+			const size_type bytesNeeded = sizeof(Rep) + cap;
+			const size_type intsNeeded = ( bytesNeeded + sizeof(int) - 1 ) / sizeof( int ); 
+			rep_ = reinterpret_cast<Rep*>( new int[ intsNeeded ] );
+
+			rep_->str[ rep_->size = sz ] = '\0';
+			rep_->capacity = cap;
+		}
+		else
+		{
+			rep_ = &nullrep_;
+		}
+	}
+
+	void quit()
+	{
+		if (rep_ != &nullrep_)
+		{
+			// The rep_ is really an array of ints. (see the allocator, above).
+			// Cast it back before delete, so the compiler won't incorrectly call destructors.
+			delete [] ( reinterpret_cast<int*>( rep_ ) );
+		}
+	}
+
+	Rep * rep_;
+	static Rep nullrep_;
+
+} ;
+
+
+inline bool operator == (const TiXmlString & a, const TiXmlString & b)
+{
+	return    ( a.length() == b.length() )				// optimization on some platforms
+	       && ( strcmp(a.c_str(), b.c_str()) == 0 );	// actual compare
+}
+inline bool operator < (const TiXmlString & a, const TiXmlString & b)
+{
+	return strcmp(a.c_str(), b.c_str()) < 0;
+}
+
+inline bool operator != (const TiXmlString & a, const TiXmlString & b) { return !(a == b); }
+inline bool operator >  (const TiXmlString & a, const TiXmlString & b) { return b < a; }
+inline bool operator <= (const TiXmlString & a, const TiXmlString & b) { return !(b < a); }
+inline bool operator >= (const TiXmlString & a, const TiXmlString & b) { return !(a < b); }
+
+inline bool operator == (const TiXmlString & a, const char* b) { return strcmp(a.c_str(), b) == 0; }
+inline bool operator == (const char* a, const TiXmlString & b) { return b == a; }
+inline bool operator != (const TiXmlString & a, const char* b) { return !(a == b); }
+inline bool operator != (const char* a, const TiXmlString & b) { return !(b == a); }
+
+TiXmlString operator + (const TiXmlString & a, const TiXmlString & b);
+TiXmlString operator + (const TiXmlString & a, const char* b);
+TiXmlString operator + (const char* a, const TiXmlString & b);
+
+
+/*
+   TiXmlOutStream is an emulation of std::ostream. It is based on TiXmlString.
+   Only the operators that we need for TinyXML have been developped.
+*/
+class TiXmlOutStream : public TiXmlString
+{
+public :
+
+	// TiXmlOutStream << operator.
+	TiXmlOutStream & operator << (const TiXmlString & in)
+	{
+		*this += in;
+		return *this;
+	}
+
+	// TiXmlOutStream << operator.
+	TiXmlOutStream & operator << (const char * in)
+	{
+		*this += in;
+		return *this;
+	}
+
+} ;
+
+#endif	// TIXML_STRING_INCLUDED
+#endif	// TIXML_USE_STL

3rdparty/tinyxml/tinyxmlerror.C

@@ -0,0 +1,53 @@
+/*
+www.sourceforge.net/projects/tinyxml
+Original code (2.0 and earlier )copyright (c) 2000-2006 Lee Thomason (www.grinninglizard.com)
+
+This software is provided 'as-is', without any express or implied 
+warranty. In no event will the authors be held liable for any 
+damages arising from the use of this software.
+
+Permission is granted to anyone to use this software for any 
+purpose, including commercial applications, and to alter it and 
+redistribute it freely, subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must
+not claim that you wrote the original software. If you use this
+software in a product, an acknowledgment in the product documentation
+would be appreciated but is not required.
+
+2. Altered source versions must be plainly marked as such, and
+must not be misrepresented as being the original software.
+
+3. This notice may not be removed or altered from any source
+distribution.
+*/
+
+#include "tinyxml.h"
+
+// The goal of the seperate error file is to make the first
+// step towards localization. tinyxml (currently) only supports
+// english error messages, but the could now be translated.
+//
+// It also cleans up the code a bit.
+//
+
+const char* TiXmlBase::errorString[ TIXML_ERROR_STRING_COUNT ] =
+{
+	"No error",
+	"Error",
+	"Failed to open file",
+	"Memory allocation failed.",
+	"Error parsing Element.",
+	"Failed to read Element name",
+	"Error reading Element value.",
+	"Error reading Attributes.",
+	"Error: empty tag.",
+	"Error reading end tag.",
+	"Error parsing Unknown.",
+	"Error parsing Comment.",
+	"Error parsing Declaration.",
+	"Error document empty.",
+	"Error null (0) or unexpected EOF found in input stream.",
+	"Error parsing CDATA.",
+	"Error when TiXmlDocument added to document, because TiXmlDocument can only be at the root.",
+};

Apps/FiniteDefElasticity/CMakeLists.txt

@@ -54,6 +54,7 @@ IF(${ENABLE_SAMG})
   FIND_PACKAGE(SAMG)
 ENDIF(${ENABLE_SAMG})
 FIND_PACKAGE(VTFWriter)
+FIND_PACKAGE(HDF5)
 
 # Required libraries
 SET(DEPLIBS ${IFEM_LIBRARIES}
@@ -87,6 +88,10 @@ IF(VTFWRITER_LIBRARIES)
   SET(DEPLIBS ${DEPLIBS} ${VTFWRITER_LIBRARIES})
 ENDIF(VTFWRITER_LIBRARIES)
 
+IF(HDF5_LIBRARIES)
+  SET(DEPLIBS ${DEPLIBS} ${HDF5_LIBRARIES})
+ENDIF(HDF5_LIBRARIES)
+
 INCLUDE_DIRECTORIES(${IFEM_INCLUDES} ${PROJECT_SOURCE_DIR}/../LinearElasticity)
 
 SET(EXECUTABLE_OUTPUT_PATH bin)

Apps/FiniteDefElasticity/Material/eigs3d.f

@@ -73,6 +73,7 @@ C
 C
 C         Entry section
 C
+      nrotd = 0
       ierr = 0
 C
       if (ipsw .gt. 0)                            then
@@ -94,13 +95,7 @@ C
          d(i) = v(i,i)
          b(i) = d(i)
          z(i) = zero
-C
-         do j = 1,3
-C
-          v(i,j) = zero
-C
-         end do ! j
-C
+         v(i,:) = zero
          v(i,i) = one
 C
       end do ! i
@@ -111,8 +106,6 @@ C
       g  = abs(d(1)) + abs(d(2)) + abs(d(3))
 C
       if (sm .lt. 1.d-13*g)                       go to 8000
-C
-      nrotd = 0
 C
       do its = 1, 50
 C

Apps/FiniteDefElasticity/Material/plas3d.f

@@ -344,7 +344,7 @@ C
      &                     f11,f22,f33,f12,f13,f23,Ypr,YY,gam,.true.)
 C
             if (ipsw .ge. 5 .and. it.lt.10)       then
-               write(iwr,"('   yield(',I1,') =')") it,yield
+               write(iwr,"('   yield(',I1,') =',1PE12.5)") it,yield
             end if
 C
             xx  = f1 - f3 * two3 * J2

Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C

@@ -470,6 +470,35 @@ bool NonlinearElasticityULMixed::evalBouMx (LocalIntegral*& elmInt,
 }
 
 
+size_t NonlinearElasticityULMixed::getNoFields (int fld) const
+{
+  if (fld < 2)
+    return nsd + 2;
+  else
+    return this->Elasticity::getNoFields(fld);
+}
+
+
+const char* NonlinearElasticityULMixed::getField1Name (size_t i,
+						       const char* prefix) const
+{
+  if (i < nsd)
+    return this->Elasticity::getField1Name(i,prefix);
+
+  static std::string name;
+  if (prefix) name = prefix + std::string(" ");
+
+  if (i == nsd)
+    name += "theta";
+  else if (i == (size_t)(nsd+1))
+    name += "p";
+  else
+    return 0;
+
+  return name.c_str();
+}
+
+
 NormBase* NonlinearElasticityULMixed::getNormIntegrand (AnaSol*) const
 {
   NonlinearElasticityULMixed* ulp;

Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h

@@ -90,6 +90,14 @@ public:
   //! returned pointer value.
   virtual NormBase* getNormIntegrand(AnaSol* = 0) const;
 
+  //! \brief Returns the number of primary/secondary solution field components.
+  //! \param[in] fld which field to consider (1=primary, 2=secondary)
+  virtual size_t getNoFields(int fld = 2) const;
+  //! \brief Returns the name of a primary solution field component.
+  //! \param[in] i Field component index
+  //! \param[in] prefix Name prefix for all components
+  virtual const char* getField1Name(size_t i, const char* prefix = 0) const;
+
 protected:
   mutable Tensor Fbar; //!< Mixed model deformation gradient
   mutable Matrix Dmat; //!< Projected mixed constitutive matrix

Apps/FiniteDefElasticity/main_NonLinEl.C

@@ -14,6 +14,8 @@
 #include "NonLinSIM.h"
 #include "SIMFiniteDefEl.h"
 #include "LinAlgInit.h"
+#include "HDF5Writer.h"
+#include "XMLWriter.h"
 #include "Profiler.h"
 #include "VTF.h"
 #include <fstream>
@@ -40,9 +42,10 @@
   \arg -nu \a nu : Number of visualization points per knot-span in u-direction
   \arg -nv \a nv : Number of visualization points per knot-span in v-direction
   \arg -nw \a nw : Number of visualization points per knot-span in w-direction
-  \arg -skip2nd : Skip VTF-output of secondary solution fields
+  \arg -hdf5 : Write primary and projected secondary solution to HDF5 file
+  \arg -skip2nd : Skip VTF/HDF5-output of secondary solution fields
   \arg -noEnergy : Skip integration of energy norms
-  \arg -saveInc \a dtSave : Time increment between each result save to VTF
+  \arg -saveInc \a dtSave : Time increment between each result save to VTF/HDF5
   \arg -dumpInc \a dtDump [raw] : Time increment between each solution dump
   \arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
   \arg -check : Data check only, read model and output to VTF (no solution)
@@ -69,6 +72,7 @@ int main (int argc, char** argv)
   int outPrec = 3;
   int n[3] = { 2, 2, 2 };
   std::vector<int> ignoredPatches;
+  bool doHDF5 = false;
   double dtSave = 0.0;
   double dtDump = 0.0;
   bool dumpWithID = true;
@@ -113,6 +117,8 @@ int main (int argc, char** argv)
       skip2nd = true;
     else if (!strcmp(argv[i],"-noEnergy"))
       energy = false;
+    else if (!strcmp(argv[i],"-hdf5"))
+      doHDF5 = true;
     else if (!strcmp(argv[i],"-saveInc") && i < argc-1)
       dtSave = atof(argv[++i]);
     else if (!strcmp(argv[i],"-outPrec") && i < argc-1)
@@ -185,10 +191,11 @@ int main (int argc, char** argv)
     std::cout <<"usage: "<< argv[0]
 	      <<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n      "
 	      <<" [-lag] [-spec] [-2D[pstrain]] [-UL] [-MX [<p>]|-mixed]"
-	      <<" [-nGauss <n>]\n       [-vtf <format>] [-nviz <nviz>]"
-	      <<" [-nu <nu>] [-nv <nv>] [-nw <nw>]\n      "
-	      <<" [-saveInc <dtSave>] [-skip2nd] [-dumpInc <dtDump> [raw]]\n   "
-	      <<"    [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
+	      <<" [-nGauss <n>]\n       [-vtf <format> [-nviz <nviz>]"
+	      <<" [-nu <nu>] [-nv <nv>] [-nw <nw>]]\n      "
+	      <<" [-saveInc <dtSave>] [-skip2nd] [-dumpInc <dtDump> [raw]]"
+	      <<" [-hdf5]\n      "
+	      <<" [-ignore <p1> <p2> ...] [-fixDup] [-checkRHS] [-check]\n";
     return 0;
   }
 
@@ -289,6 +296,20 @@ int main (int argc, char** argv)
   if (form >= 100)
     return 0; // model check
 
+  DataExporter* writer = 0;
+  if (doHDF5)
+  {
+    strtok(infile,".");
+    if (linalg.myPid == 0)
+      std::cout <<"\nWriting HDF5 file "<< infile <<".hdf5"<< std::endl;
+
+    writer = new DataExporter(true);
+    writer->registerField("u","solution",DataExporter::SIM, skip2nd ? 0 : -1);
+    writer->setFieldValue("u",model,(void*)&simulator.getSolution());
+    writer->registerWriter(new HDF5Writer(infile));
+    writer->registerWriter(new XMLWriter(infile));
+  }
+
   // Initialize the linear solver
   model->initSystem(solver,1,1);
   model->setAssociatedRHS(0,0);
@@ -321,10 +342,15 @@ int main (int argc, char** argv)
       if (!simulator.saveStep(++iStep,params.time.t,n,skip2nd))
 	return 6;
 
+      // Save solution variables to HDF5
+      if (writer)
+	writer->dumpTimeLevel();
+
       nextSave = params.time.t + dtSave;
     }
   }
 
+  if (writer) delete writer;
   if (oss) delete oss;
   return 0;
 }

Apps/LinearElasticity/main_LinEl3D.C

@@ -14,6 +14,8 @@
 #include "SIMLinEl3D.h"
 #include "SIMLinEl2D.h"
 #include "LinAlgInit.h"
+#include "HDF5Writer.h"
+#include "XMLWriter.h"
 #include "Profiler.h"
 #include <fstream>
 #include <stdlib.h>
@@ -39,6 +41,7 @@
   \arg -nu \a nu : Number of visualization points per knot-span in u-direction
   \arg -nv \a nv : Number of visualization points per knot-span in v-direction
   \arg -nw \a nw : Number of visualization points per knot-span in w-direction
+  \arg -hdf5 : Write primary and projected secondary solution to HDF5 file
   \arg -dumpASC : Dump model and solution to ASCII files for external processing
   \arg -ignore \a p1, \a p2, ... : Ignore these patches in the analysis
   \arg -eig \a iop : Eigenproblem solver to use (1...6)
@@ -74,6 +77,7 @@ int main (int argc, char** argv)
   bool checkRHS = false;
   bool vizRHS = false;
   bool fixDup = false;
+  bool dumpHDF5 = false;
   bool dumpASCII = false;
   bool twoD = false;
   char* infile = 0;
@@ -107,6 +111,8 @@ int main (int argc, char** argv)
       n[1] = atoi(argv[++i]);
     else if (!strcmp(argv[i],"-nw") && i < argc-1)
       n[2] = atoi(argv[++i]);
+    else if (!strcmp(argv[i],"-hdf5"))
+      dumpHDF5 = true;
     else if (!strcmp(argv[i],"-dumpASC"))
       dumpASCII = true;
     else if (!strcmp(argv[i],"-ignore"))
@@ -314,6 +320,19 @@ int main (int argc, char** argv)
 
   utl::profiler->start("Postprocessing");
 
+  if (dumpHDF5 && !displ.front().empty())
+  {
+    strtok(infile,".");
+    if (linalg.myPid == 0)
+      std::cout <<"\nWriting HDF5 file "<< infile <<".hdf5"<< std::endl;
+    DataExporter exporter(true);
+    exporter.registerField("u","displacement",DataExporter::SIM,-1);
+    exporter.setFieldValue("u",model,&displ.front());
+    exporter.registerWriter(new HDF5Writer(infile));
+    exporter.registerWriter(new XMLWriter(infile));
+    exporter.dumpTimeLevel();
+  }
+
   if (format >= 0)
   {
     // Write VTF-file with model geometry

CMakeLists.txt

@@ -43,6 +43,7 @@ IF(${ENABLE_SAMG})
   FIND_PACKAGE(SAMG)
 ENDIF(${ENABLE_SAMG})
 FIND_PACKAGE(VTFWriter)
+FIND_PACKAGE(HDF5)
 
 # Required libraries
 SET(DEPLIBS ${GoTrivariate_LIBRARIES} ${GoTools_LIBRARIES}
@@ -61,6 +62,7 @@ SET(INCLUDES
   ${PROJECT_SOURCE_DIR}/src/LinAlg
   ${PROJECT_SOURCE_DIR}/src/SIM
   ${PROJECT_SOURCE_DIR}/src/Utility
+  ${PROJECT_SOURCE_DIR}/3rdparty/tinyxml
 )
 
 IF(SuperLU_MT_LIBRARIES AND SuperLU_MT_INCLUDES AND "${ENABLE_SUPERLU_MT}")
@@ -104,6 +106,12 @@ IF(VTFWRITER_LIBRARIES AND VTFWRITER_INCLUDES)
   SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAS_VTFAPI=1")
 ENDIF(VTFWRITER_LIBRARIES AND VTFWRITER_INCLUDES)
 
+IF(HDF5_LIBRARIES AND HDF5_INCLUDE_DIR)
+  SET(DEPLIBS ${DEPLIBS} ${HDF5_LIBRARIES})
+  SET(INCLUDES ${INCLUDES} ${HDF5_INCLUDE_DIR})
+  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAS_HDF5=1")
+ENDIF(HDF5_LIBRARIES AND HDF5_INCLUDE_DIR)
+
 INCLUDE_DIRECTORIES(${INCLUDES})
 
 SET(EXECUTABLE_OUTPUT_PATH bin)
@@ -118,7 +126,8 @@ ENDIF(NOT WIN32)
 
 # Make the IFEM library
 
-FILE(GLOB_RECURSE IFEM_SRCS ${PROJECT_SOURCE_DIR}/src/*.[Cf])
+FILE(GLOB_RECURSE IFEM_SRCS ${PROJECT_SOURCE_DIR}/src/*.[Cf]
+                            ${PROJECT_SOURCE_DIR}/3rdparty/*.C)
 ADD_LIBRARY(IFEM ${IFEM_SRCS})
 
 # Make the Apps

cmake/Modules/FindIFEM.cmake

@@ -16,8 +16,9 @@ IF(IFEM_INCLUDES)
   FIND_PATH(SIM_INCL SIMbase.h ${PROJECT_SOURCE_DIR}/../../src/SIM)
   FIND_PATH(Alg_INCL MatVec.h  ${PROJECT_SOURCE_DIR}/../../src/LinAlg)
   FIND_PATH(Utl_INCL Vec3.h    ${PROJECT_SOURCE_DIR}/../../src/Utility)
+  FIND_PATH(Utl_3rd  tinyxml.h ${PROJECT_SOURCE_DIR}/../../3rdparty/tinyxml)
   SET(IFEM_INCLUDES ${IFEM_INCLUDES}
-    ${Int_INCL} ${SIM_INCL} ${Alg_INCL} ${Utl_INCL})
+    ${Int_INCL} ${SIM_INCL} ${Alg_INCL} ${Utl_INCL} ${Utl_3rd})
 
   FIND_LIBRARY(IFEM_LIBRARIES
     NAMES IFEM

src/ASM/ASMbase.C

@@ -479,17 +479,26 @@ void ASMbase::extractNodeVec (const Vector& globRes, Vector& nodeVec,
 }
 
 
-void ASMbase::injectNodeVec (const Vector& nodeVec, Vector& globRes,
+bool ASMbase::injectNodeVec (const Vector& nodeVec, Vector& globRes,
 			     unsigned char nndof) const
 {
   if (nndof == 0) nndof = nf;
 
+  if (nodeVec.size() < MLGN.size()*nndof)
+  {
+    std::cerr <<" *** ASMbase::injectNodeVec:: Invalid patch vector, size = "
+	      << nodeVec.size() <<" < "<< MLGN.size()*nndof << std::endl;
+    return false;
+  }
+
   for (size_t i = 0; i < MLGN.size(); i++)
   {
     int n = MLGN[i]-1;
     for (unsigned char j = 0; j < nndof; j++)
       globRes[nndof*n+j] = nodeVec[nndof*i+j];
   }
+
+  return true;
 }
 
 

src/ASM/ASMbase.h

@@ -320,7 +320,7 @@ public:
   //! \param[in] nodeVec Nodal result vector for this patch
   //! \param globVec Global solution vector in DOF-order
   //! \param[in] nndof Number of DOFs per node (the default is \a nf)
-  virtual void injectNodeVec(const Vector& nodeVec, Vector& globVec,
+  virtual bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
 			     unsigned char nndof = 0) const;
 
 protected:

src/ASM/IntegrandBase.h

@@ -333,15 +333,16 @@ public:
   //! \brief Returns whether there are any traction/flux values to write to VTF.
   virtual bool hasTractionValues() const { return false; }
 
-  //! \brief Returns the number of (secondary) solution field components.
-  virtual size_t getNoFields() const { return 0; }
-
-  //! \brief Returns the name of a (secondary) solution field component.
-  virtual const char* getFieldLabel(size_t, const char* = 0) const { return 0; }
-
   //! \brief Returns a pointer to an Integrand for solution norm evaluation.
   virtual NormBase* getNormIntegrand(AnaSol* = 0) const { return 0; }
 
+  //! \brief Returns the number of primary/secondary solution field components.
+  virtual size_t getNoFields(int = 2) const { return 0; }
+  //! \brief Returns the name of a primary solution field component.
+  virtual const char* getField1Name(size_t, const char* = 0) const { return 0; }
+  //! \brief Returns the name of a secondary solution field component.
+  virtual const char* getField2Name(size_t, const char* = 0) const { return 0; }
+
   //! \brief Returns the number of solution vectors.
   size_t getNoSolutions() const { return primsol.size(); }
 

src/ASM/SAMpatchPara.C

@@ -139,7 +139,6 @@ bool SAMpatchPara::initForAssembly (SystemVector& sysRHS,
 }
 
 
-
 bool SAMpatchPara::assembleSystem (SystemVector& sysRHS,
 				   const Matrix& eK, int iel,
 				   Vector* reactionForces) const
@@ -189,11 +188,6 @@ bool SAMpatchPara::getElmEqns (IntVec& meen, int iel, int nedof) const
   int nenod = mpmnpc[iel] - ip;
   if (nedof < 1) nedof = nenod*ndof/nnod;
 
-#ifdef USE_F77SAM
-  int neldof, neslv, neprd;
-  meen.resize(nedof,0);
-  elmeq_(madof,mmnpc+ip-1,mpmceq,meqn,nenod,&meen.front(),neldof,neslv,neprd);
-#else
   meen.clear();
   meen.reserve(nedof);
   for (int i = 0; i < nenod; i++, ip++)
@@ -202,12 +196,10 @@ bool SAMpatchPara::getElmEqns (IntVec& meen, int iel, int nedof) const
     for (int j = madof[node-1]; j < madof[node]; j++)
       meen.push_back(j);
   }
-  int neldof = meen.size();
-#endif
-  if (neldof == nedof) return true;
+  if (meen.size() == nedof) return true;
 
   std::cerr <<"SAMpatchPara::getElmEqns: Invalid element matrix dimension "
-	    << nedof <<" (should have been "<< neldof <<")"<< std::endl;
+	    << nedof <<" (should have been "<< meen.size() <<")"<< std::endl;
   return false;
 }
 

src/Integrands/AxSymmElasticity.C

@@ -254,15 +254,29 @@ bool AxSymmElasticity::evalSol (Vector& s, const Vector& N,
 }
 
 
-const char* AxSymmElasticity::getFieldLabel (size_t i, const char* prefix) const
+const char* AxSymmElasticity::getField1Name (size_t i, const char* prefix) const
 {
-  static const char* s[5] = { "s_r","s_z","s_t","s_zr", "von Mises stress" };
+  if (i > 2)
+    return 0;
+  else if (!prefix)
+    return i == 0 ? "u_r" : "u_z";
 
-  static std::string label;
-  if (prefix)
-    label = std::string(prefix) + " " + std::string(s[i]);
-  else
-    label = s[i];
+  static std::string name;
+  name = prefix + std::string(i == 0 ? " u_r" : " u_z");
 
-  return label.c_str();
+  return name.c_str();
+}
+
+
+const char* AxSymmElasticity::getField2Name (size_t i, const char* prefix) const
+{
+  if (i > 4) return 0;
+
+  static const char* s[5] = { "s_rr","s_zz","s_tt","s_zr", "von Mises stress" };
+  if (!prefix) return s[i];
+
+  static std::string name;
+  name = prefix + std::string(" ") + s[i];
+
+  return name.c_str();
 }

src/Integrands/AxSymmElasticity.h

@@ -67,13 +67,17 @@ public:
   virtual bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
 		       const Vec3& X) const;
 
-  //! \brief Returns the number of secondary solution field components.
-  virtual size_t getNoFields() const { return 5; }
-
+  //! \brief Returns the number of primary/secondary solution field components.
+  //! \param[in] fld which field set to consider (1=primary, 2=secondary)
+  virtual size_t getNoFields(int fld = 2) const { return fld > 1 ? 5 : 2; }
   //! \brief Returns the name of a secondary solution field component.
   //! \param[in] i Field component index
   //! \param[in] prefix Name prefix for all components
-  virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
+  virtual const char* getField1Name(size_t i, const char* prefix = 0) const;
+  //! \brief Returns the name of a secondary solution field component.
+  //! \param[in] i Field component index
+  //! \param[in] prefix Name prefix for all components
+  virtual const char* getField2Name(size_t i, const char* prefix = 0) const;
 
   //! \brief Returns \e true if this is an axial-symmetric problem.
   virtual bool isAxiSymmetric() const { return true; }

src/Integrands/Elasticity.C

@@ -512,35 +512,53 @@ bool Elasticity::evalSol (Vector& s, const STensorFunc& asol,
 }
 
 
-size_t Elasticity::getNoFields () const
+size_t Elasticity::getNoFields (int fld) const
 {
-  if (nsd == 2 && material->isPlaneStrain())
+  if (fld < 2)
+    return nsd;
+  else if (nsd == 2 && material->isPlaneStrain())
     return 5;
   else
     return nsd*(nsd+1)/2 + 1;
 }
 
 
-const char* Elasticity::getFieldLabel (size_t i, const char* prefix) const
+const char* Elasticity::getField1Name (size_t i, const char* prefix) const
 {
-  static const char* s[6] = { "s_xx","s_yy","s_zz","s_xy","s_yz","s_xz" };
+  if (i >= nsd) return 0;
 
-  static std::string label;
+  static const char* s[3] = { "u_x", "u_y", "u_z" };
+  if (!prefix) return s[i];
+
+  static std::string name;
+  name = prefix + std::string(" ") + s[i];
+
+  return name.c_str();
+}
+
+
+const char* Elasticity::getField2Name (size_t i, const char* prefix) const
+{
+  if (i >= this->getNoFields(2)) return 0;
+
+  static const char* s[6] = { "s_xx", "s_yy", "s_zz", "s_xy", "s_yz", "s_xz" };
+
+  static std::string name;
   if (prefix)
-    label = std::string(prefix) + " ";
+    name = std::string(prefix) + " ";
   else
-    label.clear();
+    name.clear();
 
   if (nsd == 1)
-    label += "Axial stress";
+    name += "Axial stress";
   else if (i+1 >= this->getNoFields())
-    label += "von Mises stress";
+    name += "von Mises stress";
   else if (i == 2 && this->getNoFields() == 4)
-    label += s[3];
+    name += s[3];
   else
-    label += s[i];
+    name += s[i];
 
-  return label.c_str();
+  return name.c_str();
 }
 
 
@@ -572,7 +590,7 @@ bool ElasticityNorm::initElementBou (const std::vector<int>& MNPC)
 }
 
 
-size_t ElasticityNorm::getNoFields () const
+size_t ElasticityNorm::getNoFields (int) const
 {
   return anasol ? 4 : 2;
 }

src/Integrands/Elasticity.h

@@ -126,13 +126,17 @@ public:
   //! \param[in] asol Pointer to analytical solution fields (optional)
   virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
 
-  //! \brief Returns the number of secondary solution field components.
-  virtual size_t getNoFields() const;
-
+  //! \brief Returns the number of primary/secondary solution field components.
+  //! \param[in] fld which field set to consider (1=primary, 2=secondary)
+  virtual size_t getNoFields(int fld = 2) const;
+  //! \brief Returns the name of a primary solution field component.
+  //! \param[in] i Field component index
+  //! \param[in] prefix Name prefix for all components
+  virtual const char* getField1Name(size_t i, const char* prefix = 0) const;
   //! \brief Returns the name of a secondary solution field component.
   //! \param[in] i Field component index
   //! \param[in] prefix Name prefix for all components
-  virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
+  virtual const char* getField2Name(size_t i, const char* prefix = 0) const;
 
   //! \brief Returns \e true if this is an axial-symmetric problem.
   virtual bool isAxiSymmetric() const { return false; }
@@ -261,7 +265,7 @@ public:
 		       const Vec3& X, const Vec3& normal) const;
 
   //! \brief Returns the number of norm quantities.
-  virtual size_t getNoFields() const;
+  virtual size_t getNoFields(int = 0) const;
 
 protected:
   //! \brief Returns the element norm object to use in the integration.

src/Integrands/Poisson.C

@@ -245,19 +245,28 @@ bool Poisson::evalSol (Vector& s, const VecFunc& asol, const Vec3& X) const
 }
 
 
-const char* Poisson::getFieldLabel (size_t i, const char* prefix) const
+const char* Poisson::getField1Name (size_t, const char* prefix) const
 {
-  static const char* s3[3] = { "q_x","q_y","q_z" };
+  if (!prefix) return "u";
 
-  static std::string label;
-  if (prefix)
-    label = std::string(prefix) + " ";
-  else
-    label.clear();
+  static std::string name;
+  name = prefix + std::string(" u");
 
-  label += s3[i];
+  return name.c_str();
+}
 
-  return label.c_str();
+
+const char* Poisson::getField2Name (size_t i, const char* prefix) const
+{
+  if (i >= nsd) return 0;
+
+  static const char* s[3] = { "q_x","q_y","q_z" };
+  if (!prefix) return s[i];
+
+  static std::string name;
+  name = prefix + std::string(" ") + s[i];
+
+  return name.c_str();
 }
 
 

src/Integrands/Poisson.h

@@ -112,13 +112,16 @@ public:
   //! \param[in] asol Pointer to analytical solution (optional)
   virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
 
-  //! \brief Returns the number of secondary solution field components.
-  virtual size_t getNoFields() const { return nsd; }
-
+  //! \brief Returns the number of primary/secondary solution field components.
+  //! \param[in] fld which field set to consider (1=primary, 2=secondary)
+  virtual size_t getNoFields(int fld = 2) const { return fld > 1 ? nsd : 1; }
+  //! \brief Returns the name of the primary solution field.
+  //! \param[in] prefix Name prefix
+  virtual const char* getField1Name(size_t, const char* prefix = 0) const;
   //! \brief Returns the name of a secondary solution field component.
   //! \param[in] i Field component index
   //! \param[in] prefix Name prefix for all components
-  virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
+  virtual const char* getField2Name(size_t i, const char* prefix = 0) const;
 
   //! \brief Sets up the constitutive matrix at current point.
   //! \param[out] C \f$ nsd\times nsd\f$-matrix
@@ -178,7 +181,7 @@ public:
 		       const Vec3& X) const;
 
   //! \brief Returns the number of norm quantities.
-  virtual size_t getNoFields() const { return anasol ? 3 : 1; }
+  virtual size_t getNoFields(int = 0) const { return anasol ? 3 : 1; }
 
 private:
   Poisson& problem; //!< The problem-specific data

src/SIM/NonLinSIM.C

@@ -394,7 +394,8 @@ bool NonLinSIM::saveModel (char* fileName, int format, int* nViz)
 }
 
 
-bool NonLinSIM::saveStep (int iStep, double time, int* nViz, bool psolOnly)
+bool NonLinSIM::saveStep (int iStep, double time, int* nViz,
+			  bool psolOnly, const char* vecName)
 {
   PROFILE1("NonLinSIM::saveStep");
 
@@ -416,7 +417,8 @@ bool NonLinSIM::saveStep (int iStep, double time, int* nViz, bool psolOnly)
 
   // Write solution fields
   if (!solution.empty())
-    if (!model->writeGlvS(solution.front(),nViz,iStep,nBlock,time,psolOnly))
+    if (!model->writeGlvS(solution.front(),nViz,iStep,nBlock,time,
+			  psolOnly,vecName))
       return false;
 
   // Write element norms (only when no additional visualization points are used)

src/SIM/NonLinSIM.h

@@ -99,7 +99,9 @@ public:
   //! \param[in] time Current time/load parameter
   //! \param[in] nViz Number of visualization points over each knot-span
   //! \param[in] psolOnly If \e true, skip secondary solution field output
-  bool saveStep(int iStep, double time, int* nViz, bool psolOnly = false);
+  //! \param[in] vecName Optional name of primary solution vector field
+  bool saveStep(int iStep, double time, int* nViz, bool psolOnly = false,
+		const char* vecName = 0);
 
   //! \brief Dumps the primary solution for inspection.
   //! \param[in] iStep Load/time step identifier

src/SIM/SIMbase.C

@@ -657,6 +657,10 @@ bool SIMbase::solveSystem (Vector& solution, int printSol,
     }
     std::cout << std::endl;
   }
+#if SP_DEBUG > 2
+  else
+    std::cout <<"\nSolution vector:"<< *b;
+#endif
 
   return true;
 }
@@ -1031,7 +1035,8 @@ 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)
+			 int iStep, int& nBlock, double time,
+			 bool psolOnly, const char* vecName)
 {
   if (psol.empty())
     return true;
@@ -1044,10 +1049,12 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
   Matrix field;
   size_t i, j;
   int geomID = 0;
+  const size_t pMAX = 6;
   const size_t sMAX = 21;
-  std::vector<int> vID, dID[3], sID[sMAX];
-  bool haveAsol = false;
+  std::vector<int> vID[2], dID[pMAX], sID[sMAX];
   bool scalarEq = myModel.front()->getNoFields() == 1;
+  size_t nVcomp = scalarEq ? 1 : this->getNoSpaceDim();
+  bool haveAsol = false;
   if (mySol)
     if (scalarEq)
       haveAsol = mySol->hasScalarSol() > 1;
@@ -1067,26 +1074,16 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
     if (!myModel[i]->evalSolution(field,myProblem->getSolution(),nViz))
       return false;
 
-    if (scalarEq)
-    {
-      if (!myVtf->writeNres(field.getRow(1),++nBlock,++geomID))
-	return false;
-      else
-	dID[0].push_back(nBlock);
-    }
+    if (!myVtf->writeVres(field,++nBlock,++geomID,nVcomp))
+      return false;
     else
-    {
-      if (!myVtf->writeVres(field,++nBlock,++geomID,this->getNoSpaceDim()))
+      vID[0].push_back(nBlock);
+
+    for (j = 0; j < field.rows() && j < pMAX; j++)
+      if (!myVtf->writeNres(field.getRow(1+j),++nBlock,geomID))
 	return false;
       else
-	vID.push_back(nBlock);
-
-      for (j = 0; j < field.rows() && j < 3; j++)
-	if (!myVtf->writeNres(field.getRow(1+j),++nBlock,geomID))
-	  return false;
-	else
-	  dID[j].push_back(nBlock);
-    }
+	dID[j].push_back(nBlock);
 
     if (psolOnly) continue; // skip secondary solution
 
@@ -1100,7 +1097,7 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
       if (!myVtf->writeVres(field,++nBlock,geomID))
 	return false;
       else
-	vID.push_back(nBlock);
+	vID[1].push_back(nBlock);
 
     size_t k = 0;
     for (j = 1; j <= field.rows() && k < sMAX; j++)
@@ -1156,44 +1153,37 @@ bool SIMbase::writeGlvS (const Vector& psol, const int* nViz,
 
   // Write result block identifications
 
+  bool ok = true;
   int idBlock = 10;
-  if (scalarEq)
-  {
-    if (!psolOnly)
-      if (!myVtf->writeVblk(vID,"Flux",idBlock,iStep))
-	return false;
-
-    if (!myVtf->writeSblk(dID[0],"u",++idBlock,iStep))
-      return false;
-  }
+  if (vecName)
+    ok = myVtf->writeVblk(vID[0],vecName,idBlock,iStep);
   else
-  {
-    if (!myVtf->writeDblk(vID,"Solution",idBlock,iStep))
+    ok = myVtf->writeDblk(vID[0],"Solution",idBlock,iStep);
+  if (!ok) return false;
+
+  if (scalarEq && !psolOnly)
+    if (!myVtf->writeVblk(vID[1],"Flux",idBlock+1,iStep))
       return false;
 
-    std::string label("u_x");
-    for (j = 0; j < 3 && !dID[j].empty(); j++)
-      if (!myVtf->writeSblk(dID[j],label.c_str(),++idBlock,iStep))
-	return false;
-      else
-	label[2]++;
-  }
+  for (j = 0; j < pMAX && !dID[j].empty(); j++)
+    if (!myVtf->writeSblk(dID[j],myProblem->getField1Name(j),++idBlock,iStep))
+      return false;
 
   if (psolOnly) return true;
 
-  size_t nf = myProblem->getNoFields();
+  size_t nf = myProblem->getNoFields(2);
   for (i = j = 0; i < nf && !sID[j].empty(); i++, j++)
-    if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,haveAsol?"FE":0),
+    if (!myVtf->writeSblk(sID[j],myProblem->getField2Name(i,haveAsol?"FE":0),
 			  ++idBlock,iStep)) return false;
 
   if (discretization == Spline)
     for (i = 0; i < nf && !sID[j].empty(); i++, j++)
-      if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,"Projected"),
+      if (!myVtf->writeSblk(sID[j],myProblem->getField2Name(i,"Projected"),
 			    ++idBlock,iStep)) return false;
 
   if (haveAsol)
     for (i = 0; i < nf && !sID[j].empty(); i++, j++)
-      if (!myVtf->writeSblk(sID[j],myProblem->getFieldLabel(i,"Exact"),
+      if (!myVtf->writeSblk(sID[j],myProblem->getField2Name(i,"Exact"),
 			    ++idBlock,iStep)) return false;
 
   return true;
@@ -1374,10 +1364,7 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
     myModel[i]->extractNodeVec(psol,patchSol);
     for (k = 1; k <= nf; k++)
     {
-      os <<"# FE u";
-      if (nf > 1)
-	os <<'_'<< char('w'+k);
-
+      os << myProblem->getField1Name(k,"# FE");
       for (j = 1; j <= myModel[i]->getNoNodes(); j++)
       {
 	std::pair<int,int> dofs = mySam->getNodeDOFs(j);
@@ -1396,7 +1383,7 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
     // Write the secondary solution
     for (j = 1; j <= field.rows(); j++)
     {
-      os <<"# "<< myProblem->getFieldLabel(j-1,"FE");
+      os << myProblem->getField2Name(j-1,"# FE");
       for (k = 1; k <= field.cols(); k++)
 	os <<"\n"<< field(j,k);
       os << std::endl;
@@ -1511,9 +1498,40 @@ bool SIMbase::project (Vector& psol)
     myModel[i]->extractNodeVec(psol,myProblem->getSolution());
     if (!myModel[i]->evalSolution(values,*myProblem))
       return false;
-    myModel[i]->injectNodeVec(values,ssol,values.rows());
+    if (!myModel[i]->injectNodeVec(values,ssol,values.rows()))
+      return false;
   }
   psol = ssol;
 
   return true;
 }
+
+
+bool SIMbase::extractPatchSolution (const Vector& sol,
+				    int pindx, unsigned char nndof)
+{
+  if (pindx < 0 || (size_t)pindx >= myModel.size() || sol.empty())
+    return false;
+
+  myModel[pindx]->extractNodeVec(sol,myProblem->getSolution(),nndof);
+  return true;
+}
+
+
+bool SIMbase::injectPatchSolution (Vector& sol, const Vector& locSol,
+				   int pindx, unsigned char nndof)
+{
+  if (pindx < 0 || (size_t)pindx >= myModel.size())
+    return false;
+
+  return myModel[pindx]->injectNodeVec(locSol,sol,nndof);
+}
+
+
+bool SIMbase::evalSecondarySolution (Matrix& field, int pindx) const
+{
+  if (pindx < 0 || (size_t)pindx >= myModel.size())
+    return false;
+
+  return myModel[pindx]->evalSolution(field,*myProblem);
+}

src/SIM/SIMbase.h

@@ -118,6 +118,9 @@ public:
   //! \brief Prints out problem-specific data to the given stream.
   void printProblem(std::ostream& os) const;
 
+  //! \brief Returns a pointer to the problem-specific data object.
+  const Integrand* getProblem() const { return myProblem; }
+
   //! \brief Returns the number of spatial dimensions in the model.
   size_t getNoSpaceDim() const;
   //! \brief Returns the model size in terms of number of DOFs.
@@ -126,6 +129,8 @@ public:
   size_t getNoNodes(bool unique = false) const;
   //! \brief Returns the number of solution vectors.
   size_t getNoSolutions() const;
+  //! \brief Returns the total number of patches in the model.
+  int getNoPatches() const { return nGlPatches; }
 
   //! \brief Initializes time-dependent in-homogeneous Dirichlet coefficients.
   //! \param[in] time Current time
@@ -231,14 +236,24 @@ public:
 		   int nev, int ncv, int iop, double shift,
 		   size_t iA = 0, size_t iB = 1);
 
-  //! \brief Project the secondary solution associated with \a psol.
+  //! \brief Projects the secondary solution associated with \a psol.
   //! \param psol Control point values of primary(in)/secondary(out) solution
   //!
-  //! \details The secondary solution, defined through the Integrand class,
+  //! \details The secondary solution, defined through the Integrand object,
   //! corresponding to the primary solution \a psol is projected onto the
   //! spline basis to obtain the control point values of the secondary solution.
   bool project(Vector& psol);
 
+  //! \brief Evaluates the secondary solution field for specified patch.
+  //! \param[out] field Control point values of the secondary solution field
+  //! \param[in] pindx Local patch index to evaluate solution field for
+  //!
+  //! \details The secondary solution is derived from the primary solution,
+  //! which is assumed to be stored within the Integrand object.
+  //! The solution is evaluated at the Greville points and then projected onto
+  //! the spline basis to obtain the control point values.
+  bool evalSecondarySolution(Matrix& field, int pindx) const;
+
 
   // Post-processing methods
   // =======================
@@ -271,20 +286,25 @@ public:
   //! \param[in] nViz Number of visualization points over each knot-span
   //! \param[in] iStep Load/time step identifier
   //! \param nBlock Running result block counter
-  virtual bool writeGlvV(const Vector& vec, const char* fieldName,
-			 const int* nViz, int iStep, int& nBlock) const;
+  bool writeGlvV(const Vector& vec, const char* fieldName,
+		 const int* nViz, int iStep, int& nBlock) const;
 
   //! \brief Writes solution fields for a given load/time step to the VTF-file.
-  //! \details If an analytical solution is provided, the exact stress fields
-  //! are written to the VTF-file as well.
   //! \param[in] psol Primary solution vector
   //! \param[in] nViz Number of visualization points over each knot-span
   //! \param[in] iStep Load/time step identifier
   //! \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
+  //!
+  //! \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
+  //! 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);
+			 int& nBlock, double time = 0.0, bool psolOnly = false,
+			 const char* vecName = 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
@@ -352,6 +372,7 @@ protected:
   //! \brief Creates the computational FEM model from the spline patches.
   bool createFEMmodel();
 
+public:
   //! \brief Returns the local patch index for the given global patch number.
   //! \details For serial applications this is an identity mapping only, whereas
   //! for parallel applications the local (1-based) patch index on the current
@@ -359,6 +380,31 @@ protected:
   //! If \a patchNo is not on current processor, 0 is returned.
   int getLocalPatchIndex(int patchNo) const;
 
+  //! \brief Extracts the local solution vector for a specified patch.
+  //! \param[in] sol Global primary solution vectors in DOF-order
+  //! \param[in] pindx Local patch index to extract solution vectors for
+  //! \param[in] nndof Number of DOFs per node (optional)
+  //!
+  //! \details The extracted patch-level solution vector is stored within the
+  //! Integrand \a *myProblem such that \a evalSolution can be invoked to get
+  //! the secondary solution field values within the same patch afterwards,
+  bool extractPatchSolution(const Vector& sol,
+			    int pindx, unsigned char nndof = 0);
+
+  //! \brief Injects a patch-wise solution vector into the global vector.
+  //! \param sol Global primary solution vector in DOF-order
+  //! \param[in] locSol Local solution vector associated with specified patch
+  //! \param[in] pindx Local patch index to inject solution vector for
+  //! \param[in] nndof Number of DOFs per node (optional)
+  bool injectPatchSolution(Vector& sol, const Vector& locSol,
+			   int pindx, unsigned char nndof = 0);
+
+protected:
+  //! \brief Extracts local solution vector(s) for the given patch.
+  //! \param[in] patch Pointer to the patch to extract solution vectors for
+  //! \param[in] sol Global primary solution vectors in DOF-order
+  void extractPatchSolution(const ASMbase* patch, const Vectors& sol);
+
   //! \brief Initializes material properties for integration of interior terms.
   virtual bool initMaterial(size_t) { return true; }
 
@@ -372,11 +418,6 @@ protected:
   //! \brief Finalizes the system assembly.
   bool finalizeAssembly(bool newLHSmatrix);
 
-  //! \brief Extracts local solution vector(s) for the given patch.
-  //! \param[in] patch Pointer to the patch to extract solution vectors for
-  //! \param[in] sol Global primary solution vectors in DOF-order
-  void extractPatchSolution(const ASMbase* patch, const Vectors& sol);
-
 public:
   //! \brief Enum defining the available discretization methods.
   enum Discretization { Spline, Lagrange, Spectral };

src/Utility/DataExporter.C

@@ -0,0 +1,139 @@
+// $Id$
+
+#include "DataExporter.h"
+#include <iostream>
+#include <algorithm>
+
+
+DataExporter::DataExporter(bool dynWrts) : deleteWriters(dynWrts), m_level(-1)
+{
+}
+
+DataExporter::~DataExporter()
+{
+  if (deleteWriters)
+    for (size_t i = 0; i < m_writers.size(); i++)
+      delete m_writers[i];
+}
+
+
+bool DataExporter::registerField(const std::string& name,
+                                 const std::string& description,
+                                 FieldType field, size_t size)
+{
+  if (m_entry.find(name) != m_entry.end())
+    return false;
+
+  FileEntry entry;
+  entry.description = description;
+  entry.field = field;
+  entry.size = size;
+  entry.data = NULL;
+  m_entry.insert(make_pair(name,entry));
+
+  return true;
+}
+
+bool DataExporter::registerWriter(DataWriter* writer)
+{
+  m_writers.push_back(writer);
+  return true;
+}
+
+bool DataExporter::setFieldValue(const std::string& name,
+				 void* data, void* data2)
+{
+  std::map<std::string,FileEntry>::iterator it = m_entry.find(name);
+  if (it == m_entry.end())
+    return false;
+
+  it->second.data = data;
+  it->second.data2 = data2;
+  return true;
+}
+
+bool DataExporter::dumpTimeLevel()
+{
+  if (m_level == -1)
+    m_level = getWritersTimeLevel()+1;
+
+  std::map<std::string,FileEntry>::iterator it;
+  std::vector<DataWriter*>::iterator it2;
+  for (it2 = m_writers.begin(); it2 != m_writers.end(); ++it2) {
+    (*it2)->openFile(m_level);
+    for (it = m_entry.begin(); it != m_entry.end(); ++it) {
+      if (!it->second.data)
+        return false;
+      switch (it->second.field) {
+        case VECTOR:
+          (*it2)->writeVector(m_level,*it);
+          break;
+        case SIM:
+          (*it2)->writeSIM(m_level,*it);
+          break;
+        default:
+	  std::cout <<"DataExporter: Invalid field type registered, skipping"
+		    << std::endl;
+          break;
+      }
+    }
+    (*it2)->closeFile(m_level);
+  }
+  m_level++;
+
+  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;
+
+  input->openFile(level);
+  std::map<std::string,FileEntry>::iterator it;
+  for (it = m_entry.begin(); it != m_entry.end(); ++it)
+  {
+    if (!it->second.data)
+      return false;
+    switch (it->second.field)
+    {
+      case VECTOR:
+        input->readVector(level,*it);
+        break;
+      case SIM:
+        input->readSIM(level,*it);
+        break;
+      default:
+	std::cout << "DataExporter::loadTimeLevel: Invalid field type "
+		  << "registered, skipping" << std::endl;
+        break;
+    }
+  }
+  input->closeFile(level);
+  m_level++;
+
+  return true;
+}
+
+int DataExporter::getTimeLevel()
+{
+  if (m_level == -1)
+    m_level = getWritersTimeLevel();
+
+  return m_level;
+}
+
+int DataExporter::getWritersTimeLevel()
+{
+  std::vector<int> levels;
+  std::vector<DataWriter*>::const_iterator it2;
+  for (it2 = m_writers.begin(); it2 != m_writers.end(); ++it2)
+    levels.push_back((*it2)->getLastTimeLevel());
+  return *min_element(levels.begin(),levels.end());
+}

src/Utility/DataExporter.h

@@ -0,0 +1,79 @@
+// $Id$
+
+#pragma once
+
+#include <map>
+#include <string>
+#include <vector>
+
+
+class DataWriter;
+
+class DataExporter {
+  public:
+    enum FieldType {
+      VECTOR,
+      SIM
+    };
+
+    struct FileEntry {
+      std::string description;
+      FieldType field;
+      int size;
+      void* data;
+      void* data2;
+    };
+
+    DataExporter(bool dynamicWriters = false);
+    virtual ~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);
+
+    bool registerWriter(DataWriter* writer);
+
+    bool setFieldValue(const std::string& name, void* data, void* data2=NULL);
+
+    bool dumpTimeLevel();
+
+    //! \brief Loads last time level with first registered writer by default.
+    bool loadTimeLevel(int level=-1, DataWriter* input=NULL);
+    int getTimeLevel();
+
+  protected:
+    int getWritersTimeLevel();
+
+    std::map<std::string,FileEntry> m_entry;
+    std::vector<DataWriter*> m_writers;
+    bool deleteWriters;
+    int m_level;
+};
+
+
+typedef std::pair<std::string,DataExporter::FileEntry> DataEntry;
+
+class DataWriter
+{
+protected:
+  DataWriter() {}
+
+public:
+  virtual ~DataWriter() {}
+
+  virtual int getLastTimeLevel() = 0;
+
+  virtual void openFile(int level) = 0;
+  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 void writeSIM(int level, const DataEntry& entry) = 0;
+  virtual void readSIM(int level, const DataEntry& entry) = 0;
+};

src/Utility/HDF5Writer.C

@@ -0,0 +1,240 @@
+// $Id$
+
+#include "HDF5Writer.h"
+#include "SIMbase.h"
+#include "IntegrandBase.h"
+#include <sstream>
+#include <numeric>
+#include <sys/stat.h>
+
+#ifdef HAS_HDF5
+#include <hdf5.h>
+#endif
+
+#ifdef PARALLEL_PETSC
+#include <mpi.h>
+#endif
+
+
+HDF5Writer::HDF5Writer (const std::string& name) : m_hdf5(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)
+    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()
+{
+#ifdef HAS_HDF5
+  if (m_flag == H5F_ACC_TRUNC)
+    return -1;
+
+  hid_t acc_tpl = H5P_DEFAULT;
+#ifdef PARALLEL_PETSC
+  MPI_Info info = MPI_INFO_NULL;
+  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;
+  while (1) {
+    std::stringstream str;
+    str << '/' << result;
+    if (!checkGroupExistence(m_file,str.str().c_str()))
+      break;
+    result++;
+  }
+  H5Fclose(m_file);
+
+  return result-1;
+#else
+  return -1;
+#endif
+}
+
+void HDF5Writer::openFile(int level)
+{
+#ifdef HAS_HDF5
+  hid_t acc_tpl = H5P_DEFAULT;
+#ifdef PARALLEL_PETSC
+  MPI_Info info = MPI_INFO_NULL;
+  acc_tpl = H5Pcreate (H5P_FILE_ACCESS);
+  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);
+  else
+    m_file = H5Fopen(m_hdf5.c_str(),m_flag,acc_tpl);
+
+  std::stringstream str;
+  str << '/' << level;
+  if (!checkGroupExistence(m_file,str.str().c_str())) {
+    hid_t group = H5Gcreate2(m_file,str.str().c_str(),0,H5P_DEFAULT,H5P_DEFAULT);
+    H5Gclose(group);
+  }
+#ifdef PARALLEL_PETSC
+  H5Pclose(acc_tpl);
+#endif
+#endif
+}
+
+void HDF5Writer::closeFile(int level)
+{
+#ifdef HAS_HDF5
+  H5Fflush(m_file,H5F_SCOPE_GLOBAL);
+  H5Fclose(m_file);
+  m_flag = H5F_ACC_RDWR;
+#endif
+}
+
+void HDF5Writer::readArray(int group, const std::string& name,
+                           int& len, double*& data)
+{
+#ifdef HAS_HDF5
+  hid_t set = H5Dopen2(group,name.c_str(),H5P_DEFAULT);
+  hsize_t siz = H5Dget_storage_size(set);
+  siz /= 8;
+  len = siz;
+  data = new double[siz];
+  H5Dread(set,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT,data);
+  H5Dclose(set);
+#else
+  std::cout << "HDF5Writer: compiled without HDF5 support, no data read" << std::endl;
+#endif
+}
+
+void HDF5Writer::writeArray(int group, const std::string& name,
+                            int len, void* data)
+{
+#ifdef HAS_HDF5
+#ifdef PARALLEL_PETSC
+  int lens[m_size];
+  std::fill(lens,lens+m_size,len);
+  MPI_Alltoall(lens,1,MPI_INT,lens,1,MPI_INT,MPI_COMM_WORLD);
+  int total_len=std::accumulate(lens,lens+m_size,0);
+  hsize_t start = (hsize_t)std::accumulate(lens,lens+m_rank,0);
+#else
+  int total_len = len;
+  hsize_t start = 0;
+#endif
+  hsize_t stride = 1;
+  hsize_t siz = total_len;
+  hid_t space = H5Screate_simple(1,&siz,NULL);
+  hid_t set = H5Dcreate2(group,name.c_str(),
+                         H5T_NATIVE_DOUBLE,space,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
+  hid_t file_space = H5Dget_space(set);
+  if (len ) {
+    siz = len;
+    H5Sselect_hyperslab(file_space,H5S_SELECT_SET,&start,&stride,&siz,NULL);
+    hid_t mem_space = H5Screate_simple(1,&siz,NULL);
+    H5Dwrite(set,H5T_NATIVE_DOUBLE,mem_space,file_space,H5P_DEFAULT,data);
+    H5Sclose(mem_space);
+  }
+  H5Sclose(file_space);
+  H5Dclose(set);
+#else
+  std::cout << "HDF5Writer: compiled without HDF5 support, no data written" << std::endl;
+#endif
+}
+
+void HDF5Writer::readVector(int level, const DataEntry& entry)
+{
+//  readArray(level,entry.first,entry.second.size,entry.second.data);
+}
+
+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)
+{
+#ifdef HAS_HDF5
+  SIMbase* sim = static_cast<SIMbase*>(entry.second.data);
+  Vector* sol = static_cast<Vector*>(entry.second.data2);
+  if (!sol) return;
+
+  for (int i = 0; i < sim->getNoPatches(); ++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,entry.first,siz,tmp);
+      sim->injectPatchSolution(*sol,Vector(tmp,siz),loc-1,entry.second.size);
+      delete[] tmp;
+    }
+    H5Gclose(group2);
+  }
+#endif
+}
+
+void HDF5Writer::writeSIM(int level, const DataEntry& entry)
+{
+#ifdef HAS_HDF5
+  SIMbase* sim = static_cast<SIMbase*>(entry.second.data);
+  Vector* sol = static_cast<Vector*>(entry.second.data2);
+  if (!sol) return;
+
+  const Integrand* prob = sim->getProblem();
+  for (int i = 0; i < sim->getNoPatches(); ++i) {
+    std::stringstream str;
+    str << level;
+    str << '/';
+    str << i+1;
+    hid_t group2;
+    if (checkGroupExistence(m_file,str.str().c_str()))
+      group2 = H5Gopen2(m_file,str.str().c_str(),H5P_DEFAULT);
+    else
+      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, loc > 0 ? psol.size() : 0, psol.ptr());
+    H5Gclose(group2);
+  }
+#endif
+}
+
+bool HDF5Writer::checkGroupExistence(int parent, const char* path)
+{
+#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;
+  H5Eset_auto(old_func,old_client_data);
+  return result;
+#endif
+  return false;
+}

src/Utility/HDF5Writer.h

@@ -0,0 +1,37 @@
+// $Id$
+
+#pragma once
+
+#include "DataExporter.h"
+
+
+class HDF5Writer : public DataWriter {
+  public:
+
+    HDF5Writer(const std::string& name);
+    virtual ~HDF5Writer();
+
+    int getLastTimeLevel();
+
+    void openFile(int level);
+    void closeFile(int level);
+
+    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);
+
+  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
+};

src/Utility/VTF.C

@@ -127,7 +127,7 @@ bool VTF::writeVres (const std::vector<real>& nodeResult,
   const size_t ncmp = nres/(nnod > 0 ? nnod : 1);
   if (nres != ncmp*nnod)
     return showError("Invalid size of result array",nres);
-  else if (nvc < 1)
+  else if (nvc < 1 || nvc > ncmp)
     nvc = ncmp;
 
 #ifdef HAS_VTFAPI
@@ -136,6 +136,13 @@ bool VTF::writeVres (const std::vector<real>& nodeResult,
   if (nres == 3*nnod && nvc == 3)
     for (size_t i = 0; i < nres; i++)
       resVec[i] = nodeResult[i];
+  else if (nres == nnod && nvc == 1)
+    for (size_t i = 0; i < nnod; i++)
+    {
+      // Writing a scalar solution as Z-deflection
+      resVec[3*i] = resVec[3*i+1] = 0.0f;
+      resVec[3*i+2] = nodeResult[i];
+    }
   else
     for (size_t i = 0; i < nnod; i++)
       for (size_t j = 0; j < 3; j++)

src/Utility/XMLWriter.C

@@ -0,0 +1,114 @@
+// $Id$
+
+#include "XMLWriter.h"
+#include "SIMbase.h"
+#include "IntegrandBase.h"
+
+#ifdef PARALLEL_PETSC
+#include <mpi.h>
+#endif
+
+
+XMLWriter::XMLWriter(const std::string& name) : m_xml(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
+}
+
+XMLWriter::~XMLWriter()
+{
+}
+
+int XMLWriter::getLastTimeLevel()
+{
+  int result=-1;
+  TiXmlDocument doc(m_xml.c_str());
+  doc.LoadFile();
+  TiXmlHandle handle(&doc);
+  TiXmlElement* levels = handle.FirstChild("info").
+                                            FirstChild("levels").ToElement();
+  if (levels && levels->FirstChild())
+    result = atoi(levels->FirstChild()->Value());
+
+  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);
+}
+
+void XMLWriter::closeFile(int level)
+{
+  if (!m_doc || m_rank != 0)
+    return;
+
+  TiXmlElement element2("levels");
+  TiXmlNode *pNewNode = m_node->InsertEndChild(element2);
+  char temp[8];
+  sprintf(temp,"%i",level);
+  TiXmlText value(temp);
+  pNewNode->InsertEndChild(value);
+
+  m_doc->SaveFile(m_xml);
+  delete m_doc;
+  m_doc = NULL;
+}
+
+void XMLWriter::readVector(int level, const DataEntry& entry)
+{
+}
+
+void XMLWriter::writeVector(int level, const DataEntry& entry)
+{
+  if (m_rank != 0)
+    return;
+
+  TiXmlElement element("entry");
+  element.SetAttribute("name",entry.first.c_str());
+  element.SetAttribute("description",entry.second.description.c_str());
+  element.SetAttribute("type","vector");
+  element.SetAttribute("size",entry.second.size);
+  m_node->InsertEndChild(element);
+}
+
+void XMLWriter::readSIM(int level, const DataEntry& entry)
+{
+}
+
+void XMLWriter::writeSIM(int level, const DataEntry& entry)
+{
+  if (m_rank != 0)
+    return;
+
+  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);
+
+  // secondary solutions
+  if (entry.second.size == -1)
+    for (size_t j = 0; j < prb->getNoFields(2); j++)
+      addField(prb->getField2Name(j),"secondary","field",nPatch);
+}
+
+void XMLWriter::addField(const std::string& name, const std::string& description,
+                         const std::string& type, 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);
+  m_node->InsertEndChild(element);
+}

src/Utility/XMLWriter.h

@@ -0,0 +1,34 @@
+// $Id$
+
+#pragma once
+
+#include "DataExporter.h"
+#include "tinyxml.h"
+
+
+class XMLWriter : public DataWriter {
+  public:
+    XMLWriter(const std::string& name);
+    virtual ~XMLWriter();
+
+    int getLastTimeLevel();
+
+    void openFile(int level);
+    void closeFile(int level);
+
+    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);
+
+  protected:
+    void addField(const std::string& name, const std::string& description,
+                  const std::string& type, int patches);
+    std::string m_xml;
+
+    int m_rank; // MPI rank
+    int m_size; // number of MPI nodes
+
+    TiXmlDocument* m_doc;
+    TiXmlNode* m_node;
+};