VTI (VTK XML Imagedata) Support

IO/Mesh.cpp

@@ -598,6 +598,8 @@ std::string getString( FileFormat type )
         return "silo";
     else if ( type == FileFormat::HDF5 )
         return "hdf5";
+    else if ( type == FileFormat::VTK )
+        return "vtk";
     else
         ERROR( "Invalid type" );
     return "";

IO/Mesh.h

@@ -24,7 +24,7 @@ enum class VariableType {
 };
 enum class DataType { Double, Float, Int, Null };
 enum class MeshType { PointMesh, SurfaceMesh, VolumeMesh, Unknown };
-enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 };
+enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 , VTK};
 
 
 //! Convert enums to/from strings (more future-proof than static_cast<int>)

IO/VtiWriter.cpp

@@ -0,0 +1,117 @@
+#include "IO/HDF5_IO.h"
+#include "IO/IOHelpers.h"
+#include "IO/MeshDatabase.h"
+#include "IO/Writer.h"
+#include "IO/silo.h"
+#include "IO/xmlvtk.h"
+#include "common/MPI.h"
+#include "common/Utilities.h"
+
+#include <algorithm>
+#include <memory>
+#include <set>
+#include <sys/stat.h>
+#include <vector>
+
+static void writeVti(
+    const std::string &fullpath, const IO::MeshDataStruct &meshData)
+{
+    const IO::DomainMesh &mesh = dynamic_cast<IO::DomainMesh &>( *meshData.mesh );
+    RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz );
+
+    VTIWriter vti = VTIWriter(std::string(fullpath));
+    vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, 
+                      ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz);    
+
+    vti.setSpacing( 1.0 , 1.0 , 1.0 );
+    vti.setOrigin(0,0,0);
+    vti.setCompress();                      
+
+    for ( size_t i = 0; i < meshData.vars.size(); i++ ) 
+    {
+        const auto &var = *meshData.vars[i];
+        if ( var.precision == IO::DataType::Double ) {
+            vti.addCellData( var.name , "Float64" , "binary" , var.dim , (unsigned char*) var.data.begin() );
+        } else if ( var.precision == IO::DataType::Float ) {
+            Array<float> data2( var.data.size() );
+            data2.copy( var.data );
+            vti.addCellData( var.name , "Float32" , "binary" , var.dim , (unsigned char*)  data2.begin());
+        } else if ( var.precision == IO::DataType::Int ) {
+            Array<int> data2( var.data.size() );
+            data2.copy( var.data );
+            vti.addCellData( var.name , "Int32" , "binary" , var.dim , (unsigned char*) var.data.begin() );
+        } else {
+            ERROR( "Unsupported format" );
+        }
+    }
+
+    vti.write();
+}
+
+void writeVtiSummary(
+    const std::vector<IO::MeshDatabase> &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename )
+{
+    const IO::DomainMesh &mesh = dynamic_cast<IO::DomainMesh &>( *meshData.mesh );
+    RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz );
+    PVTIWriter pvti = PVTIWriter( filename );
+    int rank = 0;
+    for ( const auto &data : meshes_written ) 
+    {        
+         for ( const auto &tmp : data.domains ) 
+            {  
+                RankInfoStruct info( rank, mesh.nprocx, mesh.nprocy, mesh.nprocz );
+                char filename[100];
+                sprintf( filename, "%05i.vti", rank );
+
+                VTIWriter vti = VTIWriter( filename ); 
+                vti.setWholeExtent( info.ix * mesh.nx, info.jy * mesh.ny , info.kz * mesh.nz, 
+                                  ( info.ix + 1 ) * mesh.nx, ( info.jy + 1 ) * mesh.ny, ( info.kz + 1 ) * mesh.nz);    
+                vti.setSpacing( 1.0 , 1.0 , 1.0 );
+                vti.setOrigin(0,0,0);
+                vti.setCompress();                             
+
+                for ( size_t i = 0; i < meshData.vars.size(); i++ ) 
+                {
+                    const auto &var = *meshData.vars[i];
+                    if ( var.precision == IO::DataType::Double ) {
+                        vti.addCellData( var.name , "Float64" , "binary" , var.dim , nullptr );
+                    } else if ( var.precision == IO::DataType::Float ) {
+                        vti.addCellData( var.name , "Float32" , "binary" , var.dim , nullptr );
+                    } else if ( var.precision == IO::DataType::Int ) {
+                        vti.addCellData( var.name , "Int32" , "binary"   , var.dim , nullptr );
+                    } else {
+                        ERROR( "Unsupported format" );
+                    }
+                }            
+
+                pvti.addVTIWriter(vti);           
+                rank++;
+
+            }
+    }    
+    pvti.write();
+}
+
+std::vector<IO::MeshDatabase> writeMeshesVti( const std::vector<IO::MeshDataStruct> &meshData,
+    const std::string &path, int rank )
+{
+    std::vector<IO::MeshDatabase> meshes_written;
+    char filename[100], fullpath[200];
+    sprintf( filename, "%05i.vti", rank );
+    sprintf( fullpath, "%s/%s", path.c_str(), filename );
+
+    for ( size_t i = 0; i < meshData.size(); i++ ) {
+//        auto mesh = meshData[i].mesh;
+    	auto database = getDatabase( fullpath , meshData[i], IO::FileFormat::VTK, rank );
+
+	if ( database.meshClass == "DomainMesh" ) {
+        	writeVti( fullpath, meshData[i] );
+	} else {
+	        ERROR( "Unknown mesh class or not implemented for vtk/vti output" );
+	}
+
+        meshes_written.push_back( database  );
+    }
+    return meshes_written;
+}
+

IO/Writer.cpp

@@ -15,7 +15,7 @@
 #include <vector>
 
 
-enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN };
+enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN, VTK };
 
 
 /****************************************************
@@ -26,7 +26,10 @@ std::vector<IO::MeshDatabase> writeMeshesSilo(
 void writeSiloSummary( const std::vector<IO::MeshDatabase> &, const std::string & );
 std::vector<IO::MeshDatabase> writeMeshesHDF5(
     const std::vector<IO::MeshDataStruct> &, const std::string &, IO::FileFormat, int, Xdmf & );
-
+std::vector<IO::MeshDatabase> writeMeshesVti( const std::vector<IO::MeshDataStruct> &meshData,
+    const std::string &path,  int rank );
+void writeVtiSummary(
+        const std::vector<IO::MeshDatabase> &meshes_written,const IO::MeshDataStruct &meshData, const std::string &filename );
 
 /****************************************************
  * Recursively create the subdirectory               *
@@ -90,6 +93,8 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap
         global_IO_format = Format::SILO;
     else if ( format == "hdf5" )
         global_IO_format = Format::HDF5;
+    else if ( format == "vtk" )
+        global_IO_format = Format::VTK;        
     else
         ERROR( "Unknown format" );
     int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank();
@@ -100,9 +105,13 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap
             filename = global_IO_path + "/summary.LBM";
         else if ( global_IO_format == Format::SILO || global_IO_format == Format::HDF5 )
             filename = global_IO_path + "/LBM.visit";
+        else if ( global_IO_format == Format::VTK)
+            filename = global_IO_path + "/LBM.pvd";
         else
             ERROR( "Unknown format" );
         auto fid = fopen( filename.c_str(), "wb" );
+        if ( global_IO_format == Format::VTK)
+            fprintf( fid, "<?xml version=\"1.0\"?>\n<VTKFile type=\"Collection\" version=\"0.1\">\n  <Collection>\n  </Collection>\n</VTKFile>\n" );
         fclose( fid );
     }
 }
@@ -278,13 +287,12 @@ static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
     return meshes_written;
 }
 
-
 /****************************************************
  * Write the mesh data                               *
  ****************************************************/
 void IO::writeData( const std::string &subdir, const std::vector<IO::MeshDataStruct> &meshData,
-    const Utilities::MPI &comm )
-{
+    const Utilities::MPI &comm , int timestep )    
+{      
     if ( global_IO_path.empty() )
         IO::initialize();
     PROFILE_START( "writeData" );
@@ -310,6 +318,9 @@ void IO::writeData( const std::string &subdir, const std::vector<IO::MeshDataStr
     } else if ( global_IO_format == Format::HDF5 ) {
         // Write hdf5
         meshes_written = writeMeshesHDF5( meshData, path, IO::FileFormat::HDF5, rank, xmf );
+    } else if ( global_IO_format == Format::VTK ) {
+        // Write vti (vtk/xml)
+         meshes_written = writeMeshesVti( meshData, path, rank );        
     } else {
         ERROR( "Unknown format" );
     }
@@ -328,6 +339,8 @@ void IO::writeData( const std::string &subdir, const std::vector<IO::MeshDataStr
             writeSiloSummary( meshes_written, path + "/summary.silo" );
         } else if ( global_IO_format == Format::HDF5 ) {
             xmf.write( path + "/summary.xmf" );
+        } else if  ( global_IO_format == Format::VTK ) {
+            writeVtiSummary( meshes_written, meshData[0], path + "/summary.pvti" );
         }
         // Add the timestep to the global summary file
         if ( global_IO_format == Format::OLD || global_IO_format == Format::NEW ) {
@@ -345,6 +358,32 @@ void IO::writeData( const std::string &subdir, const std::vector<IO::MeshDataStr
             FILE *fid     = fopen( filename.c_str(), "ab" );
             fprintf( fid, "%s/summary.xmf\n", subdir.c_str() );
             fclose( fid );
+        } else if ( global_IO_format == Format::VTK ) {
+            auto filename = global_IO_path + "/LBM.pvd";
+
+            FILE *fid     = fopen( filename.c_str(), "r" );
+            fseek(fid, 0, SEEK_END);
+
+            long filesize = ftell(fid);
+            char *buffer = (char *)malloc(filesize + 1);       
+            rewind(fid);
+            fread(buffer, 1, filesize, fid);
+
+            fclose(fid);
+            //std:cout << "Oi Estou aqui!!!! " << std::endl;
+            char newrow[1024];
+            snprintf(newrow, 1024, "\t<DataSet timestep=\"%d\" part=\"0\" file=\"%s/summary.pvti\"/>\n  </Collection>\n</VTKFile>\n", timestep, subdir.c_str());
+
+            char *pos = strstr(buffer, "  </Collection>\n</VTKFile>\n");
+            pos[0] = '\0';
+
+            fid     = fopen( filename.c_str(), "w" );
+            fputs(buffer,fid);
+            fputs(newrow,fid);
+            fclose(fid);
+
+            free(buffer);
+
         } else {
             ERROR( "Unknown format" );
         }

IO/Writer.h

@@ -38,7 +38,7 @@ void initialize(
  * @param[in] comm          The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof)
  */
 void writeData( const std::string &subdir, const std::vector<IO::MeshDataStruct> &meshData,
-    const Utilities::MPI &comm );
+    const Utilities::MPI &comm , int timestep = -1);
 
 
 /*!
@@ -48,12 +48,14 @@ void writeData( const std::string &subdir, const std::vector<IO::MeshDataStruct>
  * @param[in] meshData      The data to write
  * @param[in] comm          The comm to use for writing (usually MPI_COMM_WORLD or a dup thereof)
  */
+
+
 inline void writeData(
     int timestep, const std::vector<IO::MeshDataStruct> &meshData, const Utilities::MPI &comm )
 {
     char subdir[100];
     sprintf( subdir, "vis%03i", timestep );
-    writeData( subdir, meshData, comm );
+    writeData( subdir, meshData, comm , timestep);
 }