added: VTU output

renamed HDF5toVTF to reflect that it can now output vtu files as well

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

Apps/HDF5toVTx/HDF5toVTx.C

@@ -7,7 +7,7 @@
 //!
 //! \author Arne Morten Kvarving / SINTEF
 //!
-//! \brief Convert a HDF5 results database to VTF for visualisation.
+//! \brief Convert a HDF5 results database to VTF/VTU for visualisation.
 //!
 //==============================================================================
 
@@ -22,6 +22,7 @@
 #include "ASMs3D.h"
 #include "ElementBlock.h"
 #include "VTF.h"
+#include "VTU.h"
 
 typedef std::map< std::string,std::vector<XMLWriter::Entry> > ProcessList;
 typedef std::map< std::string, std::vector<int> > VTFList;
@@ -111,6 +112,7 @@ void writePatchGeometry(ASMbase* patch, int id, VTF& myVtf, int* nViz)
   myVtf.writeGrid(lvb,str.str().c_str());
 }
 
+
 std::vector<RealArray*> generateFEModel(std::vector<ASMbase*> patches,
                                         int dims, int* n)
 {
@@ -133,6 +135,7 @@ std::vector<RealArray*> generateFEModel(std::vector<ASMbase*> patches,
   return result;
 }
 
+
 int main (int argc, char** argv)
 {
   int format = 0;
@@ -165,7 +168,7 @@ int main (int argc, char** argv)
 
   if (!infile) {
     std::cout <<"usage: "<< argv[0]
-              <<" <inputfile> [<vtffile>] [-nviz <nviz>] [-ndump <ndump>]"
+              <<" <inputfile> [<vtffile>] [<vtufile>] [-nviz <nviz>] [-ndump <ndump>]"
               <<" [-basis <basis>] [-1D|-2D]"<< std::endl;
     return 0;
   }
@@ -192,7 +195,11 @@ int main (int argc, char** argv)
   ProcessList processlist;
   std::map<std::string, std::vector<ASMbase*> > patches;
   std::vector<RealArray*> FEmodel;
-  VTF myVtf(vtffile,1);
+  VTF* myVtf;
+  if (strstr(vtffile,".vtf"))
+    myVtf = new VTF(vtffile,1);
+  else
+    myVtf = new VTU(vtffile,1);
 
   for (it = entry.begin(); it != entry.end(); ++it) {
     if (!it->basis.empty() && it->type != "restart") {
@@ -214,12 +221,13 @@ int main (int argc, char** argv)
   else
     gpatches = patches.begin()->second;
   for (int i=0;i<pit->second[0].patches;++i)
-    writePatchGeometry(gpatches[i],i+1,myVtf,n);
+    writePatchGeometry(gpatches[i],i+1,*myVtf,n);
   FEmodel = generateFEModel(gpatches,dims,n);
 
   bool ok = true;
   int block = 0;
   double time=0;
+  int k=1;
   for (int i = 0; i <= levels && ok; i += skip) {
     if (levels > 0) std::cout <<"\nTime level "<< i << " (t=" << time << ")" << std::endl;
     VTFList vlist, slist;
@@ -243,7 +251,7 @@ int main (int argc, char** argv)
               ok &= writeFieldPatch(tmp.getRow(r),1,*patches[pit->first][j],
                                     FEmodel[j],j+1,
                                     block,it->name.substr(pos,end),vlist,
-                                    slist,myVtf);
+                                    slist,*myVtf);
               pos = end+1;
             }
           }
@@ -251,21 +259,22 @@ int main (int argc, char** argv)
             ok &= writeFieldPatch(vec,it->components,
                                   *patches[pit->first][j],
                                   FEmodel[j],j+1,
-                                  block,it->name,vlist,slist,myVtf);
+                                  block,it->name,vlist,slist,*myVtf);
           }
         }
       }
     }
-    writeFieldBlocks(vlist,slist,myVtf,i+1);
+    writeFieldBlocks(vlist,slist,*myVtf,k);
 
     if (ok)
-      myVtf.writeState(i+1,"Time %g",time,0);
+      myVtf->writeState(k++,"Time %g",time,0);
     else
       return 3;
     pit = processlist.begin();
     time += pit->second.begin()->timestep*skip;
   }
   hdf.closeFile(levels,true);
+  delete myVtf;
 
   return 0;
 }

Apps/HDF5toVTx/VTU.C

@@ -0,0 +1,159 @@
+// $Id$
+//==============================================================================
+//!
+//! \file VTU.C
+//!
+//! \date Jun 14 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Basic VTU file writer class
+//!
+//==============================================================================
+
+
+#include "ElementBlock.h"
+#include "VTU.h"
+
+#include <fstream>
+#include <iomanip>
+#include <sstream>
+
+
+VTU::VTU(const char* base, int format)
+  : VTF(NULL,format), m_base(base)
+{
+  m_base = m_base.substr(0,m_base.rfind('.'));
+}
+
+
+VTU::~VTU()
+{
+  for (size_t i=0;i<m_geom.size();++i)
+    delete m_geom[i];
+  m_geom.clear();
+}
+
+
+bool VTU::writeGrid(const ElementBlock* block, const char* name)
+{
+  m_geom.push_back(block);
+  return true;
+}
+
+
+bool VTU::writeVres(const std::vector<real>& field, int blockID, int geomID,
+                    int components)
+{
+  m_field[blockID].data = new Vector(field.data(),field.size());
+  m_field[blockID].components = components;
+  m_field[blockID].patch = geomID;
+  return true;
+}
+
+
+bool VTU::writeNres(const std::vector<real>& field, int blockID, int geomID)
+{
+  m_field[blockID].data = new Vector(field.data(),field.size());
+  m_field[blockID].components = 1;
+  m_field[blockID].patch = geomID;
+  return true;
+}
+
+
+bool VTU::writeVblk(const std::vector<int>& vBlockIDs, 
+                    const char* resultName, int idBlock, int iStep)
+{
+  for (size_t i=0;i<vBlockIDs.size();++i)
+    m_field[vBlockIDs[i]].name = resultName;
+  return true;
+}
+
+
+bool VTU::writeSblk(const std::vector<int>& sBlockIDs, 
+                    const char* resultName, int idBlock, int iStep)
+{
+  for (size_t i=0;i<sBlockIDs.size();++i)
+    m_field[sBlockIDs[i]].name = resultName;
+  return true;
+}
+
+
+bool VTU::writeState(int iStep, const char* fmt, real refValue, int refType)
+{
+  std::ofstream file;
+  std::stringstream str;
+  str << m_base << "-" << std::setfill('0') << std::setw(5) << iStep-1 << ".vtu";
+  file.open(str.str().c_str());
+
+  file << "<?xml version=\"1.0\"?>" << std::endl;
+  file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
+  file << "\t<UnstructuredGrid>" << std::endl;
+  for (size_t i=0;i<m_geom.size();++i) {
+    file << "\t\t<Piece NumberOfCells=\"" << m_geom[i]->getNoElms()
+         << "\" NumberOfPoints=\"" << m_geom[i]->getNoNodes()
+         << "\">" << std::endl;
+
+    // dump geometry
+    file << "\t\t\t<Points>" << std::endl;
+    file << "\t\t\t\t<DataArray type=\"Float32\" Name=\"Coordinates\""
+         << " NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
+    file << "\t\t\t\t\t";
+    for (std::vector<Vec3>::const_iterator it  = m_geom[i]->begin_XYZ();
+                                           it != m_geom[i]->end_XYZ();++it) {
+      it->print(file);
+      file << " ";
+    }
+    file << std::endl << "\t\t\t\t</DataArray>" << std::endl;
+    file << "\t\t\t</Points>" << std::endl;
+
+    file << "\t\t\t<Cells>" << std::endl;
+    file << "\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\""
+         << " NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
+    file << "\t\t\t\t\t";
+    for (size_t j=0;j<m_geom[i]->getNoElms()*m_geom[i]->getNoElmNodes();++j)
+      file << m_geom[i]->getElements()[j] << " ";
+    file << std::endl << "\t\t\t\t</DataArray>" << std::endl;
+
+    file << "\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\""
+         << " NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
+    file << "\t\t\t\t\t";
+    std::string type = "9";
+    if (m_geom[i]->getNoElmNodes() == 8)
+      type="12";
+    for (size_t k=0;k<m_geom[i]->getNoElms();++k)
+      file << type << " ";
+    file << std::endl << "\t\t\t\t</DataArray>" << std::endl;
+
+    file << "\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\""
+         << " NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
+    file << "\t\t\t\t\t";
+    for (size_t k=0;k<m_geom[i]->getNoElms();++k)
+      file << (k+1)*m_geom[i]->getNoElmNodes() << " ";
+    file << std::endl << "\t\t\t\t</DataArray>" << std::endl;
+    file << "\t\t\t</Cells>" << std::endl;
+
+    // now add point datas
+    file << "\t\t\t<PointData Scalars=\"scalars\">" << std::endl;
+    for (std::map<int,FieldInfo>::iterator it2  = m_field.begin();
+                                           it2 != m_field.end();++it2) {
+      if (it2->second.patch == (int)i+1) {
+        file << "\t\t\t\t<DataArray type=\"Float32\" Name=\"" << it2->second.name << "\""
+             << " NumberOfComponents=\"" << it2->second.components
+             << "\" format=\"ascii\">" << std::endl;
+        file << "\t\t\t\t\t";
+        for (size_t k=0;k<it2->second.data->size();++k)
+          file << (*it2->second.data)[k] << " ";
+        delete it2->second.data;
+        file << std::endl << "\t\t\t\t</DataArray>" << std::endl;
+      }
+    }
+    file << "\t\t\t</PointData>" << std::endl;
+    file << "\t\t</Piece>" << std::endl;
+  }
+  file << "\t</UnstructuredGrid>" << std::endl;
+  file << "</VTKFile>" << std::endl;
+  file.close();
+  m_field.clear();
+  return true;
+}

Apps/HDF5toVTx/VTU.h

@@ -0,0 +1,53 @@
+#pragma once
+
+// $Id$
+//==============================================================================
+//!
+//! \file VTU.h
+//!
+//! \date Jun 14 2011
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Basic VTU file writer class
+//!
+//==============================================================================
+
+
+#include "MatVec.h"
+#include "VTF.h"
+#include <string>
+
+
+class ElementBlock;
+
+
+class VTU : public VTF {
+  public:
+    VTU(const char* base, int format);
+    virtual ~VTU();
+
+    bool writeGrid(const ElementBlock* lvb, const char* name);
+
+    bool writeVres(const std::vector<double>& field, int blockID,
+                   int geomID, int components);
+    bool writeNres(const std::vector<double>& vec, int blockID, int geomID);
+
+    bool writeVblk(const std::vector<int>& vBlockIDs, 
+                   const char* resultName, int idBlock, int iStep=1);
+
+    bool writeSblk(const std::vector<int>& sBlockIDs,
+                   const char* resultName, int idBlock, int iStep=1);
+
+    bool writeState(int iStep, const char* fmt, real refValue, int refType=0);
+  protected:
+    std::string m_base;
+    std::vector<const ElementBlock*> m_geom;
+    struct FieldInfo {
+      Vector* data;
+      int components;
+      int patch;
+      std::string name;
+    };
+    std::map<int,FieldInfo> m_field;
+};

CMakeLists.txt

@@ -145,9 +145,9 @@ ADD_EXECUTABLE(LinEl ${LinEl_SRCS})
 TARGET_LINK_LIBRARIES(LinEl IFEM ${DEPLIBS})
 
 IF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
-  FILE(GLOB_RECURSE HDF2VTF_SRCS ${PROJECT_SOURCE_DIR}/Apps/HDF5toVTF/*.C)
+  FILE(GLOB_RECURSE HDF2VTF_SRCS ${PROJECT_SOURCE_DIR}/Apps/HDF5toVTx/*.C)
-  ADD_EXECUTABLE(HDF5toVTF ${HDF2VTF_SRCS})
+  ADD_EXECUTABLE(HDF5toVTx ${HDF2VTF_SRCS})
-  TARGET_LINK_LIBRARIES(HDF5toVTF IFEM ${DEPLIBS})
+  TARGET_LINK_LIBRARIES(HDF5toVTx IFEM ${DEPLIBS})
 ENDIF(HDF5_LIBRARIES AND VTFWRITER_LIBRARIES)
 
 # Regression tests

src/Utility/VTF.C

@@ -33,6 +33,11 @@ real VTF::vecOffset[3] = { 0.0, 0.0, 0.0 };
 
 VTF::VTF (const char* filename, int type)
 {
+  if (!filename) {
+    myFile = 0;
+    return;
+  }
+
   myState = 0;
 #if HAS_VTFAPI == 1
   // Create the VTF file object
@@ -71,6 +76,8 @@ VTF::VTF (const char* filename, int type)
 
 VTF::~VTF ()
 {
+  if (!myFile) return;
+
   std::vector<int> geomID(myBlocks.size());
   for (size_t i = 0; i < myBlocks.size(); i++)
   {
@@ -82,8 +89,6 @@ VTF::~VTF ()
     showError("Error writing geometry");
 
 #if HAS_VTFAPI == 1
-  if (!myFile) return;
-
   size_t i;
   for (i = 0; i < myDBlock.size(); i++)
     if (myDBlock[i])
@@ -120,8 +125,6 @@ VTF::~VTF ()
 
   delete myFile;
 #elif HAS_VTFAPI == 2
-  if (!myFile) return;
-
   size_t i;
   for (i = 0; i < myDBlock.size(); i++)
     if (myDBlock[i])

src/Utility/VTF.h

@@ -50,15 +50,15 @@ public:
   //! \brief Writes the FE geometry to the VTF-file.
   //! \param[in] g The FE grid that all results written are referred to
   //! \param[in] partname Name of the geometry being written
-  bool writeGrid(const ElementBlock* g, const char* partname);
+  virtual bool writeGrid(const ElementBlock* g, const char* partname);
 
   //! \brief Writes a block of scalar nodal results to the VTF-file.
   //! \param[in] nodeResult Vector of nodal results,
   //!            length must equal the number of nodes in the geometry block
   //! \param[in] idBlock Result block identifier
   //! \param[in] gID Geometry block identifier
-  bool writeNres(const std::vector<real>& nodeResult,
+  virtual bool writeNres(const std::vector<real>& nodeResult,
-		 int idBlock = 1, int gID = 1);
+		         int idBlock = 1, int gID = 1);
   //! \brief Writes a block of scalar element results to the VTF-file.
   //! \param[in] elmResult Vector of element results (one per element)
   //!            length must equal the number of elements in the geometry block
@@ -72,8 +72,8 @@ public:
   //! \param[in] idBlock Result block identifier
   //! \param[in] gID Geometry block identifier
   //! \param[in] nvc Number of components per node in \a nodeResult
-  bool writeVres(const std::vector<real>& nodeResult,
+  virtual bool writeVres(const std::vector<real>& nodeResult,
-		 int idBlock = 1, int gID = 1, size_t nvc = 0);
+		         int idBlock = 1, int gID = 1, size_t nvc = 0);
   //! \brief Writes a block of point vector results to the VTF-file.
   //! \details This method creates a separate geometry block consisting of the
   //! attack points of the result vectors, since they are independent of the
@@ -94,8 +94,8 @@ public:
   //! \param[in] resultName Name of the result quantity
   //! \param[in] idBlock Scalar block identifier
   //! \param[in] iStep Load/Time step identifier
-  bool writeSblk(const std::vector<int>& sBlockIDs,
+  virtual bool writeSblk(const std::vector<int>& sBlockIDs,
-                 const char* resultName = 0, int idBlock = 1, int iStep = 1);
+                         const char* resultName = 0, int idBlock = 1, int iStep = 1);
   //! \brief Writes a vector block definition to the VTF-file.
   //! \param[in] vBlockID The result block that makes up this vector block
   //! \param[in] resultName Name of the result quantity
@@ -108,8 +108,8 @@ public:
   //! \param[in] resultName Name of the result quantity
   //! \param[in] idBlock Vector block identifier
   //! \param[in] iStep Load/Time step identifier
-  bool writeVblk(const std::vector<int>& vBlockIDs,
+  virtual bool writeVblk(const std::vector<int>& vBlockIDs,
-                 const char* resultName = 0, int idBlock = 1, int iStep = 1);
+                         const char* resultName = 0, int idBlock = 1, int iStep = 1);
   //! \brief Writes a displacement block definition to the VTF-file.
   //! \param[in] dBlockIDs All result blocks that make up the displacement block
   //! \param[in] resultName Name of the result quantity
@@ -123,7 +123,7 @@ public:
   //! \param[in] fmt Format string for step name
   //! \param[in] refValue Reference value for the step (time, frequency, etc.)
   //! \param[in] refType Reference value type (0=Time, 1=Frequency, 2=Load case)
-  bool writeState(int iStep, const char* fmt, real refValue, int refType = 0);
+  virtual bool writeState(int iStep, const char* fmt, real refValue, int refType = 0);
 
   //! \brief Returns the pointer to a geometry block.
   const ElementBlock* getBlock(int geomID) const { return myBlocks[geomID-1]; }