Updating convertIO

IO/Writer.cpp

@@ -1,28 +1,69 @@
 #include "IO/Writer.h"
-#include "IO/MeshDatabase.h"
 #include "IO/IOHelpers.h"
+#include "IO/MeshDatabase.h"
 #include "IO/silo.h"
 #include "common/MPI.h"
 #include "common/Utilities.h"
 
-#include <sys/stat.h>
 #include <algorithm>
-#include <vector>
-#include <set>
 #include <memory>
+#include <set>
+#include <sys/stat.h>
+#include <vector>
 
 
-
-enum class Format { OLD, NEW, SILO, UNKNOWN }; 
-
+enum class Format { OLD, NEW, SILO, UNKNOWN };
 
 
 /****************************************************
-* Initialize the writer                             *
-****************************************************/
+ * Recursively create the subdirectory               *
+ ****************************************************/
+static void recursiveMkdir( const std::string &path, mode_t mode )
+{
+    // Iterate through the root directories until we create the desired path
+    for ( size_t pos = 0; pos < path.size(); ) {
+        // slide backwards in string until next slash found
+        pos++;
+        for ( ; pos < path.size(); pos++ ) {
+            if ( path[pos] == '/' || path[pos] == 92 )
+                break;
+        }
+        // Create the temporary path
+        auto path2 = path.substr( 0, pos );
+        // Check if the temporary path exists
+        struct stat status;
+        int result = stat( path2.data(), &status );
+        if ( result == 0 ) {
+            // if there is a part of the path that already exists make sure it is really a directory
+            if ( !S_ISDIR( status.st_mode ) ) {
+                ERROR(
+                    "Error in recursiveMkdir...\n"
+                    "    Cannot create directories in path = " +
+                    path +
+                    "\n    because some intermediate item in path exists and is NOT a directory" );
+            }
+            continue;
+        }
+        // Create the directory and test the result
+        result = mkdir( path2.data(), mode );
+        if ( result != 0 ) {
+            // Maybe another rank created the directory, check
+            int result = stat( path2.data(), &status );
+            if ( result != 0 && !S_ISDIR( status.st_mode ) )
+                ERROR( "Error in Utilities::recursiveMkdir...\n"
+                       "    Cannot create directory  = " +
+                       path2 );
+        }
+    }
+}
+
+
+/****************************************************
+ * Initialize the writer                             *
+ ****************************************************/
 static std::string global_IO_path;
 static Format global_IO_format = Format::UNKNOWN;
-void IO::initialize( const std::string& path, const std::string& format, bool append )
+void IO::initialize( const std::string &path, const std::string &format, bool append )
 {
     if ( path.empty() )
         global_IO_path = ".";
@@ -35,161 +76,170 @@ void IO::initialize( const std::string& path, const std::string& format, bool ap
     else if ( format == "silo" )
         global_IO_format = Format::SILO;
     else
-        ERROR("Unknown format");
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
-    if ( !append && rank==0 ) {
-        mkdir(path.c_str(),S_IRWXU|S_IRGRP);
+        ERROR( "Unknown format" );
+    int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank();
+    if ( !append && rank == 0 ) {
+        recursiveMkdir( path, S_IRWXU | S_IRGRP );
         std::string filename;
-        if ( global_IO_format==Format::OLD || global_IO_format==Format::NEW )
+        if ( global_IO_format == Format::OLD || global_IO_format == Format::NEW )
             filename = global_IO_path + "/summary.LBM";
-        else if ( global_IO_format==Format::SILO )
+        else if ( global_IO_format == Format::SILO )
             filename = global_IO_path + "/LBM.visit";
         else
-            ERROR("Unknown format");
-        auto fid = fopen(filename.c_str(),"wb");
-        fclose(fid);
+            ERROR( "Unknown format" );
+        auto fid = fopen( filename.c_str(), "wb" );
+        fclose( fid );
     }
 }
 
 
 // Write the mesh data in the original format
-static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO::MeshDataStruct>& meshData, const std::string& path )
+static std::vector<IO::MeshDatabase> writeMeshesOrigFormat(
+    const std::vector<IO::MeshDataStruct> &meshData, const std::string &path, int rank )
 {
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
     std::vector<IO::MeshDatabase> meshes_written;
-    for (size_t i=0; i<meshData.size(); i++) {
+    for ( size_t i = 0; i < meshData.size(); i++ ) {
         char domainname[100], filename[100], fullpath[200];
-        sprintf(domainname,"%05i",rank);
-        sprintf(filename,"%s.%05i",meshData[i].meshName.c_str(),rank);
-	    sprintf(fullpath,"%s/%s",path.c_str(),filename);
-        FILE *fid = fopen(fullpath,"wb");
-        INSIST(fid!=NULL,std::string("Error opening file: ")+fullpath);
+        sprintf( domainname, "%05i", rank );
+        sprintf( filename, "%s.%05i", meshData[i].meshName.c_str(), rank );
+        sprintf( fullpath, "%s/%s", path.c_str(), filename );
+        FILE *fid = fopen( fullpath, "wb" );
+        INSIST( fid != NULL, std::string( "Error opening file: " ) + fullpath );
         std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
         IO::MeshDatabase mesh_entry;
-        mesh_entry.name = meshData[i].meshName;
-        mesh_entry.type = meshType(*mesh);
+        mesh_entry.name      = meshData[i].meshName;
+        mesh_entry.type      = meshType( *mesh );
         mesh_entry.meshClass = meshData[i].mesh->className();
-        mesh_entry.format = 1;
+        mesh_entry.format    = 1;
         IO::DatabaseEntry domain;
-        domain.name = domainname;
-        domain.file = filename;
+        domain.name   = domainname;
+        domain.file   = filename;
         domain.offset = 0;
-        mesh_entry.domains.push_back(domain);
+        mesh_entry.domains.push_back( domain );
         if ( !meshData[i].vars.empty() ) {
-            printf("Warning: variables are not supported with this format (original)\n");
-            //for (size_t j=0; j<meshData[i].vars.size(); j++)
+            printf( "Warning: variables are not supported with this format (original)\n" );
+            // for (size_t j=0; j<meshData[i].vars.size(); j++)
             //    mesh_entry.variables.push_back( meshData[i].vars[j]->name );
         }
         const std::string meshClass = mesh->className();
-        if ( meshClass=="PointList" ) {
+        if ( meshClass == "PointList" ) {
             // List of points
-            std::shared_ptr<IO::PointList> pointlist = std::dynamic_pointer_cast<IO::PointList>(mesh);
-            const std::vector<Point>& P = pointlist->points;
-            for (size_t i=0; i<P.size(); i++) {
+            std::shared_ptr<IO::PointList> pointlist =
+                std::dynamic_pointer_cast<IO::PointList>( mesh );
+            const std::vector<Point> &P = pointlist->points;
+            for ( size_t i = 0; i < P.size(); i++ ) {
                 double x[3];
-                x[0] = P[i].x;  x[1] = P[i].y;  x[2] = P[i].z;
-                fwrite(x,sizeof(double),3,fid);
+                x[0] = P[i].x;
+                x[1] = P[i].y;
+                x[2] = P[i].z;
+                fwrite( x, sizeof( double ), 3, fid );
             }
-        } else if ( meshClass=="TriList" || meshClass=="TriMesh" ) {
+        } else if ( meshClass == "TriList" || meshClass == "TriMesh" ) {
             // Triangle mesh
-            std::shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
-            const std::vector<Point>& A = trilist->A;
-            const std::vector<Point>& B = trilist->B;
-            const std::vector<Point>& C = trilist->C;
-            for (size_t i=0; i<A.size(); i++) {
+            std::shared_ptr<IO::TriList> trilist = IO::getTriList( mesh );
+            const std::vector<Point> &A          = trilist->A;
+            const std::vector<Point> &B          = trilist->B;
+            const std::vector<Point> &C          = trilist->C;
+            for ( size_t i = 0; i < A.size(); i++ ) {
                 double tri[9];
-                tri[0] = A[i].x;  tri[1] = A[i].y;  tri[2] = A[i].z;
-                tri[3] = B[i].x;  tri[4] = B[i].y;  tri[5] = B[i].z;
-                tri[6] = C[i].x;  tri[7] = C[i].y;  tri[8] = C[i].z;
-                fwrite(tri,sizeof(double),9,fid);
+                tri[0] = A[i].x;
+                tri[1] = A[i].y;
+                tri[2] = A[i].z;
+                tri[3] = B[i].x;
+                tri[4] = B[i].y;
+                tri[5] = B[i].z;
+                tri[6] = C[i].x;
+                tri[7] = C[i].y;
+                tri[8] = C[i].z;
+                fwrite( tri, sizeof( double ), 9, fid );
             }
-        } else if ( meshClass=="DomainMesh" ) {
+        } else if ( meshClass == "DomainMesh" ) {
             // This format was never supported with the old format
         } else {
-            ERROR("Unknown mesh");
+            ERROR( "Unknown mesh" );
         }
-        fclose(fid);
+        fclose( fid );
         std::sort( mesh_entry.variables.begin(), mesh_entry.variables.end() );
-        mesh_entry.variables.erase( std::unique( mesh_entry.variables.begin(), mesh_entry.variables.end() ), mesh_entry.variables.end() );
-        meshes_written.push_back(mesh_entry);
+        mesh_entry.variables.erase(
+            std::unique( mesh_entry.variables.begin(), mesh_entry.variables.end() ),
+            mesh_entry.variables.end() );
+        meshes_written.push_back( mesh_entry );
     }
     return meshes_written;
 }
 
 
 // Create the database entry for the mesh data
-static IO::MeshDatabase getDatabase( const std::string& filename, const IO::MeshDataStruct& mesh, int format )
+static IO::MeshDatabase getDatabase(
+    const std::string &filename, const IO::MeshDataStruct &mesh, int format, int rank )
 {
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
     char domainname[100];
-    sprintf(domainname,"%s_%05i",mesh.meshName.c_str(),rank);
+    sprintf( domainname, "%s_%05i", mesh.meshName.c_str(), rank );
     // Create the MeshDatabase
     IO::MeshDatabase database;
-    database.name = mesh.meshName;
-    database.type = meshType(*(mesh.mesh));
+    database.name      = mesh.meshName;
+    database.type      = meshType( *( mesh.mesh ) );
     database.meshClass = mesh.mesh->className();
-    database.format = format;
+    database.format    = format;
     // Write the mesh
     IO::DatabaseEntry domain;
-    domain.name = domainname;
-    domain.file = filename;
+    domain.name   = domainname;
+    domain.file   = filename;
     domain.offset = -1;
-    database.domains.push_back(domain);
+    database.domains.push_back( domain );
     // Write the variables
-    for (size_t i=0; i<mesh.vars.size(); i++) {
+    for ( size_t i = 0; i < mesh.vars.size(); i++ ) {
         // Add basic variable info
         IO::VariableDatabase info;
         info.name = mesh.vars[i]->name;
         info.type = mesh.vars[i]->type;
-        info.dim = mesh.vars[i]->dim;
-        database.variables.push_back(info);
+        info.dim  = mesh.vars[i]->dim;
+        database.variables.push_back( info );
         // Add domain variable info
         IO::DatabaseEntry variable;
-        variable.name = mesh.vars[i]->name;
-        variable.file = filename;
+        variable.name   = mesh.vars[i]->name;
+        variable.file   = filename;
         variable.offset = -1;
-        std::pair<std::string,std::string> key(domain.name,mesh.vars[i]->name);
-        database.variable_data.insert( 
-            std::pair<std::pair<std::string,std::string>,IO::DatabaseEntry>(key,variable) );
+        std::pair<std::string, std::string> key( domain.name, mesh.vars[i]->name );
+        database.variable_data.insert(
+            std::pair<std::pair<std::string, std::string>, IO::DatabaseEntry>( key, variable ) );
     }
     return database;
 }
 
 
 // Write a mesh (and variables) to a file
-static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
-    const IO::MeshDataStruct& mesh, int format )
+static IO::MeshDatabase write_domain(
+    FILE *fid, const std::string &filename, const IO::MeshDataStruct &mesh, int format, int rank )
 {
     const int level = 0;
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
     // Create the MeshDatabase
-    IO::MeshDatabase database = getDatabase( filename, mesh, format );
+    IO::MeshDatabase database = getDatabase( filename, mesh, format, rank );
     // Write the mesh
-    IO::DatabaseEntry& domain = database.domains[0];
-    domain.offset = ftell(fid);
-    std::pair<size_t,void*> data = mesh.mesh->pack(level);
-    fprintf(fid,"Mesh: %s-%05i: %lu\n",mesh.meshName.c_str(),rank,data.first);
-    fwrite(data.second,1,data.first,fid);
-    fprintf(fid,"\n");
-    delete [] (char*) data.second;
+    IO::DatabaseEntry &domain      = database.domains[0];
+    domain.offset                  = ftell( fid );
+    std::pair<size_t, void *> data = mesh.mesh->pack( level );
+    fprintf( fid, "Mesh: %s-%05i: %lu\n", mesh.meshName.c_str(), rank, data.first );
+    fwrite( data.second, 1, data.first, fid );
+    fprintf( fid, "\n" );
+    delete[]( char * ) data.second;
     // Write the variables
-    for (size_t i=0; i<mesh.vars.size(); i++) {
-        std::pair<std::string,std::string> key(domain.name,mesh.vars[i]->name);
-        IO::DatabaseEntry& variable = database.variable_data[key];
-        variable.offset = ftell(fid);
-        int dim = mesh.vars[i]->dim;
-        int type = static_cast<int>(mesh.vars[i]->type);
-        size_t N = mesh.vars[i]->data.length();
-        if ( type == static_cast<int>(IO::VariableType::NullVariable) ) {
-            ERROR("Variable type not set");
+    for ( size_t i = 0; i < mesh.vars.size(); i++ ) {
+        std::pair<std::string, std::string> key( domain.name, mesh.vars[i]->name );
+        IO::DatabaseEntry &variable = database.variable_data[key];
+        variable.offset             = ftell( fid );
+        int dim                     = mesh.vars[i]->dim;
+        int type                    = static_cast<int>( mesh.vars[i]->type );
+        size_t N                    = mesh.vars[i]->data.length();
+        if ( type == static_cast<int>( IO::VariableType::NullVariable ) ) {
+            ERROR( "Variable type not set" );
         }
-        size_t N_mesh = mesh.mesh->numberPointsVar(mesh.vars[i]->type);
-        ASSERT(N==dim*N_mesh);
-        fprintf(fid,"Var: %s-%05i-%s: %i, %i, %lu, %lu, double\n",
-            database.name.c_str(), rank, variable.name.c_str(),
-            dim, type, N_mesh, N*sizeof(double) );
-        fwrite(mesh.vars[i]->data.data(),sizeof(double),N,fid);
-        fprintf(fid,"\n");
+        size_t N_mesh = mesh.mesh->numberPointsVar( mesh.vars[i]->type );
+        ASSERT( N == dim * N_mesh );
+        fprintf( fid, "Var: %s-%05i-%s: %i, %i, %lu, %lu, double\n", database.name.c_str(), rank,
+            variable.name.c_str(), dim, type, N_mesh, N * sizeof( double ) );
+        fwrite( mesh.vars[i]->data.data(), sizeof( double ), N, fid );
+        fprintf( fid, "\n" );
     }
     return database;
 }
@@ -198,11 +248,12 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
 #ifdef USE_SILO
 // Write a PointList mesh (and variables) to a file
 template<class TYPE>
-static void writeSiloPointMesh( DBfile *fid, const IO::PointList& mesh, const std::string& meshname )
+static void writeSiloPointMesh(
+    DBfile *fid, const IO::PointList &mesh, const std::string &meshname )
 {
-    const auto& points = mesh.getPoints();
-    std::vector<TYPE> x(points.size()), y(points.size()), z(points.size());
-    for (size_t i=0; i<x.size(); i++) {
+    const auto &points = mesh.getPoints();
+    std::vector<TYPE> x( points.size() ), y( points.size() ), z( points.size() );
+    for ( size_t i = 0; i < x.size(); i++ ) {
         x[i] = points[i].x;
         y[i] = points[i].y;
         z[i] = points[i].z;
@@ -210,28 +261,29 @@ static void writeSiloPointMesh( DBfile *fid, const IO::PointList& mesh, const st
     const TYPE *coords[] = { x.data(), y.data(), z.data() };
     silo::writePointMesh<TYPE>( fid, meshname, 3, points.size(), coords );
 }
-static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
+static void writeSiloPointList(
+    DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database )
 {
-    const IO::PointList& mesh = dynamic_cast<IO::PointList&>( *meshData.mesh );
+    const IO::PointList &mesh  = dynamic_cast<IO::PointList &>( *meshData.mesh );
     const std::string meshname = database.domains[0].name;
     if ( meshData.precision == IO::DataType::Double ) {
         writeSiloPointMesh<double>( fid, mesh, meshname );
     } else if ( meshData.precision == IO::DataType::Float ) {
         writeSiloPointMesh<float>( fid, mesh, meshname );
     } else {
-        ERROR("Unsupported format");
+        ERROR( "Unsupported format" );
     }
-    const auto& points = mesh.getPoints();
-    std::vector<double> x(points.size()), y(points.size()), z(points.size());
-    for (size_t i=0; i<x.size(); i++) {
+    const auto &points = mesh.getPoints();
+    std::vector<double> x( points.size() ), y( points.size() ), z( points.size() );
+    for ( size_t i = 0; i < x.size(); i++ ) {
         x[i] = points[i].x;
         y[i] = points[i].y;
         z[i] = points[i].z;
     }
     const double *coords[] = { x.data(), y.data(), z.data() };
     silo::writePointMesh( fid, meshname, 3, points.size(), coords );
-    for (size_t i=0; i<meshData.vars.size(); i++) {
-        const IO::Variable& var = *meshData.vars[i];
+    for ( size_t i = 0; i < meshData.vars.size(); i++ ) {
+        const IO::Variable &var = *meshData.vars[i];
         if ( var.precision == IO::DataType::Double ) {
             silo::writePointMeshVariable( fid, meshname, var.name, var.data );
         } else if ( var.precision == IO::DataType::Float ) {
@@ -243,27 +295,27 @@ static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData,
             data2.copy( var.data );
             silo::writePointMeshVariable( fid, meshname, var.name, data2 );
         } else {
-            ERROR("Unsupported format");
+            ERROR( "Unsupported format" );
         }
     }
 }
 // Write a TriMesh mesh (and variables) to a file
 template<class TYPE>
-static void writeSiloTriMesh( DBfile *fid, const IO::TriMesh& mesh, const std::string& meshname )
+static void writeSiloTriMesh( DBfile *fid, const IO::TriMesh &mesh, const std::string &meshname )
 {
-    const auto& points = mesh.vertices->getPoints();
-    std::vector<TYPE> x(points.size()), y(points.size()), z(points.size());
-    for (size_t i=0; i<x.size(); i++) {
+    const auto &points = mesh.vertices->getPoints();
+    std::vector<TYPE> x( points.size() ), y( points.size() ), z( points.size() );
+    for ( size_t i = 0; i < x.size(); i++ ) {
         x[i] = points[i].x;
         y[i] = points[i].y;
         z[i] = points[i].z;
     }
     const TYPE *coords[] = { x.data(), y.data(), z.data() };
-    const int *tri[] = { mesh.A.data(), mesh.B.data(), mesh.C.data() };
+    const int *tri[]     = { mesh.A.data(), mesh.B.data(), mesh.C.data() };
     silo::writeTriMesh<TYPE>( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri );
 }
-static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
-    const IO::TriMesh& mesh, IO::MeshDatabase database )
+static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct &meshData,
+    const IO::TriMesh &mesh, IO::MeshDatabase database )
 {
     const std::string meshname = database.domains[0].name;
     if ( meshData.precision == IO::DataType::Double ) {
@@ -271,11 +323,11 @@ static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
     } else if ( meshData.precision == IO::DataType::Float ) {
         writeSiloTriMesh<float>( fid, mesh, meshname );
     } else {
-        ERROR("Unsupported format");
+        ERROR( "Unsupported format" );
     }
-    for (size_t i=0; i<meshData.vars.size(); i++) {
-        const IO::Variable& var = *meshData.vars[i];
-        auto type = static_cast<silo::VariableType>( var.type );
+    for ( size_t i = 0; i < meshData.vars.size(); i++ ) {
+        const IO::Variable &var = *meshData.vars[i];
+        auto type               = static_cast<silo::VariableType>( var.type );
         if ( var.precision == IO::DataType::Double ) {
             silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, type );
         } else if ( var.precision == IO::DataType::Float ) {
@@ -287,35 +339,40 @@ static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
             data2.copy( var.data );
             silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, type );
         } else {
-            ERROR("Unsupported format");
+            ERROR( "Unsupported format" );
         }
     }
 }
-static void writeSiloTriMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
+static void writeSiloTriMesh(
+    DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database )
 {
-    const IO::TriMesh& mesh = dynamic_cast<IO::TriMesh&>( *meshData.mesh );
+    const IO::TriMesh &mesh = dynamic_cast<IO::TriMesh &>( *meshData.mesh );
     writeSiloTriMesh2( fid, meshData, mesh, database );
 }
-static void writeSiloTriList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
+static void writeSiloTriList(
+    DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database )
 {
     auto mesh = getTriMesh( meshData.mesh );
     writeSiloTriMesh2( fid, meshData, *mesh, database );
 }
 // Write a DomainMesh mesh (and variables) to a file
-static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
+static void writeSiloDomainMesh(
+    DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database )
 {
-    const IO::DomainMesh& mesh = dynamic_cast<IO::DomainMesh&>( *meshData.mesh );
+    const IO::DomainMesh &mesh = dynamic_cast<IO::DomainMesh &>( *meshData.mesh );
     RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz );
-    std::array<double,6> range = { info.ix*mesh.Lx/info.nx, (info.ix+1)*mesh.Lx/info.nx, 
-                                   info.jy*mesh.Ly/info.ny, (info.jy+1)*mesh.Ly/info.ny,
-                                   info.kz*mesh.Lz/info.nz, (info.kz+1)*mesh.Lz/info.nz };
-    std::array<int,3> N = { mesh.nx, mesh.ny, mesh.nz };
-    auto meshname = database.domains[0].name;
+    std::array<double, 6> range = { info.ix * mesh.Lx / info.nx,
+        ( info.ix + 1 ) * mesh.Lx / info.nx, info.jy * mesh.Ly / info.ny,
+        ( info.jy + 1 ) * mesh.Ly / info.ny, info.kz * mesh.Lz / info.nz,
+        ( info.kz + 1 ) * mesh.Lz / info.nz };
+    std::array<int, 3> N        = { mesh.nx, mesh.ny, mesh.nz };
+    auto meshname               = database.domains[0].name;
     silo::writeUniformMesh<3>( fid, meshname, range, N );
-    silo::write<int>( fid, meshname+"_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } );
-    for (size_t i=0; i<meshData.vars.size(); i++) {
-        const auto& var = *meshData.vars[i];
-        auto type = static_cast<silo::VariableType>( var.type );
+    silo::write<int>(
+        fid, meshname + "_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } );
+    for ( size_t i = 0; i < meshData.vars.size(); i++ ) {
+        const auto &var = *meshData.vars[i];
+        auto type       = static_cast<silo::VariableType>( var.type );
         if ( var.precision == IO::DataType::Double ) {
             silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, type );
         } else if ( var.precision == IO::DataType::Float ) {
@@ -327,65 +384,66 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData
             data2.copy( var.data );
             silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, type );
         } else {
-            ERROR("Unsupported format");
+            ERROR( "Unsupported format" );
         }
     }
 }
 // Write a mesh (and variables) to a file
-static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string& filename,
-    const IO::MeshDataStruct& mesh, int format )
+static IO::MeshDatabase write_domain_silo(
+    DBfile *fid, const std::string &filename, const IO::MeshDataStruct &mesh, int format, int rank )
 {
     // Create the MeshDatabase
-    auto database = getDatabase( filename, mesh, format );
-    if ( database.meshClass=="PointList" ) {
+    auto database = getDatabase( filename, mesh, format, rank );
+    if ( database.meshClass == "PointList" ) {
         writeSiloPointList( fid, mesh, database );
-    } else if ( database.meshClass=="TriMesh" ) {
+    } else if ( database.meshClass == "TriMesh" ) {
         writeSiloTriMesh( fid, mesh, database );
-    } else if ( database.meshClass=="TriList" ) {
+    } else if ( database.meshClass == "TriList" ) {
         writeSiloTriList( fid, mesh, database );
-    } else if ( database.meshClass=="DomainMesh" ) {
+    } else if ( database.meshClass == "DomainMesh" ) {
         writeSiloDomainMesh( fid, mesh, database );
     } else {
-        ERROR("Unknown mesh class");
+        ERROR( "Unknown mesh class" );
     }
     return database;
 }
 // Write the summary file for silo
-std::pair<int,int> getSiloMeshType( const std::string& meshClass )
+std::pair<int, int> getSiloMeshType( const std::string &meshClass )
 {
     int meshType = 0;
-    int varType = 0;
-    if ( meshClass=="PointList" ) {
+    int varType  = 0;
+    if ( meshClass == "PointList" ) {
         meshType = DB_POINTMESH;
         varType  = DB_POINTVAR;
-    } else if ( meshClass=="TriMesh" ) {
+    } else if ( meshClass == "TriMesh" ) {
         meshType = DB_UCDMESH;
         varType  = DB_UCDVAR;
-    } else if ( meshClass=="TriList" ) {
+    } else if ( meshClass == "TriList" ) {
         meshType = DB_UCDMESH;
         varType  = DB_UCDVAR;
-    } else if ( meshClass=="DomainMesh" ) {
+    } else if ( meshClass == "DomainMesh" ) {
         meshType = DB_QUAD_RECT;
         varType  = DB_QUADVAR;
     } else {
-        ERROR("Unknown mesh class");
+        ERROR( "Unknown mesh class" );
     }
     return std::make_pair( meshType, varType );
 }
-void writeSiloSummary( const std::vector<IO::MeshDatabase>& meshes_written, const std::string& filename )
+void writeSiloSummary(
+    const std::vector<IO::MeshDatabase> &meshes_written, const std::string &filename )
 {
     auto fid = silo::open( filename, silo::CREATE );
-    for ( const auto& data : meshes_written ) {
+    for ( const auto &data : meshes_written ) {
         auto type = getSiloMeshType( data.meshClass );
         std::vector<int> meshTypes( data.domains.size(), type.first );
         std::vector<int> varTypes( data.domains.size(), type.second );
         std::vector<std::string> meshNames;
-        for ( const auto& tmp : data.domains )
+        for ( const auto &tmp : data.domains )
             meshNames.push_back( tmp.file + ":" + tmp.name );
         silo::writeMultiMesh( fid, data.name, meshNames, meshTypes );
-        for (const auto& variable : data.variables ) {
+        for ( const auto &variable : data.variables ) {
             std::vector<std::string> varnames;
-            for ( const auto& tmp : data.domains )
+            for ( const auto &tmp : data.domains )
                 varnames.push_back( tmp.file + ":" + variable.name );
             silo::writeMultiVar( fid, variable.name, varnames, varTypes );
         }
@@ -396,113 +454,111 @@ void writeSiloSummary( const std::vector<IO::MeshDatabase>& meshes_written, cons
 
 
 // Write the mesh data in the new format
-static std::vector<IO::MeshDatabase> writeMeshesNewFormat( 
-    const std::vector<IO::MeshDataStruct>& meshData, const std::string& path, int format )
+static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
+    const std::vector<IO::MeshDataStruct> &meshData, const std::string &path, int format, int rank )
 {
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
     std::vector<IO::MeshDatabase> meshes_written;
     char filename[100], fullpath[200];
-    sprintf(filename,"%05i",rank);
-    sprintf(fullpath,"%s/%s",path.c_str(),filename);
-    FILE *fid = fopen(fullpath,"wb");
-    for (size_t i=0; i<meshData.size(); i++) {
+    sprintf( filename, "%05i", rank );
+    sprintf( fullpath, "%s/%s", path.c_str(), filename );
+    FILE *fid = fopen( fullpath, "wb" );
+    for ( size_t i = 0; i < meshData.size(); i++ ) {
         std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
-        meshes_written.push_back( write_domain(fid,filename,meshData[i],format) );
+        meshes_written.push_back( write_domain( fid, filename, meshData[i], format, rank ) );
     }
-    fclose(fid);
+    fclose( fid );
     return meshes_written;
 }
 
 
 // Write the mesh data to silo
-static std::vector<IO::MeshDatabase> writeMeshesSilo( 
-    const std::vector<IO::MeshDataStruct>& meshData, const std::string& path, int format )
+static std::vector<IO::MeshDatabase> writeMeshesSilo(
+    const std::vector<IO::MeshDataStruct> &meshData, const std::string &path, int format, int rank )
 {
 #ifdef USE_SILO
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
     std::vector<IO::MeshDatabase> meshes_written;
     char filename[100], fullpath[200];
-    sprintf(filename,"%05i.silo",rank);
-    sprintf(fullpath,"%s/%s",path.c_str(),filename);
+    sprintf( filename, "%05i.silo", rank );
+    sprintf( fullpath, "%s/%s", path.c_str(), filename );
     auto fid = silo::open( fullpath, silo::CREATE );
-    for (size_t i=0; i<meshData.size(); i++) {
+    for ( size_t i = 0; i < meshData.size(); i++ ) {
         auto mesh = meshData[i].mesh;
-        meshes_written.push_back( write_domain_silo(fid,filename,meshData[i],format) );
+        meshes_written.push_back( write_domain_silo( fid, filename, meshData[i], format, rank ) );
     }
     silo::close( fid );
     return meshes_written;
 #else
-    ERROR("Application built without silo support");
+    NULL_USE( meshData );
+    NULL_USE( path );
+    NULL_USE( format );
+    NULL_USE( rank );
+    ERROR( "Application built without silo support" );
     return std::vector<IO::MeshDatabase>();
 #endif
-}        
+}
 
 
 /****************************************************
-* Write the mesh data                               *
-****************************************************/
-void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStruct>& meshData, const Utilities::MPI& comm )
+ * Write the mesh data                               *
+ ****************************************************/
+void IO::writeData( const std::string &subdir, const std::vector<IO::MeshDataStruct> &meshData,
+    const Utilities::MPI &comm )
 {
     if ( global_IO_path.empty() )
-        IO::initialize( );
-    PROFILE_START("writeData");
-    int rank = Utilities::MPI(MPI_COMM_WORLD).getRank();
+        IO::initialize();
+    PROFILE_START( "writeData" );
+    int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank();
     // Check the meshData before writing
-    for ( const auto& data : meshData ) {
+    for ( const auto &data : meshData ) {
         if ( !data.check() )
-            ERROR("Error in meshData");
+            ERROR( "Error in meshData" );
     }
     // Create the output directory
     std::string path = global_IO_path + "/" + subdir;
-    if ( rank == 0 ) {
-        mkdir(path.c_str(),S_IRWXU|S_IRGRP);
-    }
-    comm.barrier();
+    recursiveMkdir( path, S_IRWXU | S_IRGRP );
     // Write the mesh files
     std::vector<IO::MeshDatabase> meshes_written;
     if ( global_IO_format == Format::OLD ) {
         // Write the original triangle format
-        meshes_written = writeMeshesOrigFormat( meshData, path );
+        meshes_written = writeMeshesOrigFormat( meshData, path, rank );
     } else if ( global_IO_format == Format::NEW ) {
         // Write the new format (double precision)
-        meshes_written = writeMeshesNewFormat( meshData, path, 2 );
+        meshes_written = writeMeshesNewFormat( meshData, path, 2, rank );
     } else if ( global_IO_format == Format::SILO ) {
         // Write silo
-        meshes_written = writeMeshesSilo( meshData, path, 4 );
+        meshes_written = writeMeshesSilo( meshData, path, 4, rank );
     } else {
-        ERROR("Unknown format");
+        ERROR( "Unknown format" );
     }
     // Gather a complete list of files on rank 0
-    meshes_written = gatherAll(meshes_written,comm);
+    meshes_written = gatherAll( meshes_written, comm );
     // Write the summary files
     if ( rank == 0 ) {
         // Write the summary file for the current timestep
         char filename[200];
-        sprintf(filename,"%s/LBM.summary",path.c_str());
-        write(meshes_written,filename);
-        // Write summary silo file if needed
-        #ifdef USE_SILO
+        sprintf( filename, "%s/LBM.summary", path.c_str() );
+        write( meshes_written, filename );
+// Write summary silo file if needed
+#ifdef USE_SILO
         if ( global_IO_format == Format::SILO ) {
-            sprintf(filename,"%s/summary.silo",path.c_str());
-            writeSiloSummary(meshes_written,filename);
+            sprintf( filename, "%s/summary.silo", path.c_str() );
+            writeSiloSummary( meshes_written, filename );
         }
-        #endif
+#endif
         // Add the timestep to the global summary file
         if ( global_IO_format == Format::OLD || global_IO_format == Format::NEW ) {
-            auto filename = global_IO_path+"/summary.LBM";
-            FILE *fid = fopen(filename.c_str(),"ab");
-            fprintf(fid,"%s/\n",subdir.c_str());
-            fclose(fid);
+            auto filename = global_IO_path + "/summary.LBM";
+            FILE *fid     = fopen( filename.c_str(), "ab" );
+            fprintf( fid, "%s/\n", subdir.c_str() );
+            fclose( fid );
         } else if ( global_IO_format == Format::SILO ) {
-            auto filename = global_IO_path+"/LBM.visit";
-            FILE *fid = fopen(filename.c_str(),"ab");
-            fprintf(fid,"%s/summary.silo\n",subdir.c_str());
-            fclose(fid);
+            auto filename = global_IO_path + "/LBM.visit";
+            FILE *fid     = fopen( filename.c_str(), "ab" );
+            fprintf( fid, "%s/summary.silo\n", subdir.c_str() );
+            fclose( fid );
         } else {
-            ERROR("Unknown format");
+            ERROR( "Unknown format" );
         }
     }
-    PROFILE_STOP("writeData");
+    PROFILE_STOP( "writeData" );
 }
-
-