OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Working on writer

Browse Source
This commit is contained in:
Mark Berrill 2017-01-31 08:22:29 -05:00
parent 3434dcb43c
commit 099443a3ab
7 changed files with 568 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
IO/Mesh.cpp
View File

@ -29,6 +29,80 @@ Mesh::~Mesh( )
} }
/****************************************************
* MeshDataStruct *
****************************************************/
bool MeshDataStruct::check() const
{
enum VariableType { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
bool pass = mesh != nullptr;
for ( const auto& var : vars ) {
pass = pass && static_cast<int>(var->type)>=1 && static_cast<int>(var->type)<=3;
pass = pass && !var->data.empty();
}
if ( !pass )
return false;
const std::string& meshClass = mesh->className();
if ( meshClass == "PointList" ) {
const auto mesh2 = dynamic_cast<IO::PointList*>( mesh.get() );
if ( mesh2 == nullptr )
return false;
for ( const auto& var : vars ) {
if ( var->type == IO::NodeVariable ) {
pass = pass && var->data.size(0)==mesh2->points.size() && var->data.size(1)==var->dim;
} else if ( var->type == IO::EdgeVariable ) {
ERROR("Invalid type for PointList");
} else if ( var->type == IO::SurfaceVariable ) {
ERROR("Invalid type for PointList");
} else if ( var->type == IO::VolumeVariable ) {
ERROR("Invalid type for PointList");
} else {
ERROR("Invalid variable type");
}
}
} else if ( meshClass == "TriMesh" || meshClass == "TriList" ) {
const auto mesh2 = getTriMesh( mesh );
if ( mesh2 == nullptr )
return false;
for ( const auto& var : vars ) {
if ( var->type == IO::NodeVariable ) {
pass = pass && var->data.size(0)==mesh2->vertices->points.size() && var->data.size(1)==var->dim;
} else if ( var->type == IO::EdgeVariable ) {
ERROR("Not finished");
} else if ( var->type == IO::SurfaceVariable ) {
ERROR("Not finished");
} else if ( var->type == IO::VolumeVariable ) {
pass = pass && var->data.size(0)==mesh2->A.size() && var->data.size(1)==var->dim;
} else {
ERROR("Invalid variable type");
}
}
} else if ( meshClass == "DomainMesh" ) {
const auto mesh2 = dynamic_cast<IO::DomainMesh*>( mesh.get() );
if ( mesh2 == nullptr )
return false;
for ( const auto& var : vars ) {
if ( var->type == IO::NodeVariable ) {
pass = pass && (int) var->data.size(0)==(mesh2->nx+1) && (int) var->data.size(1)==(mesh2->ny+1)
&& (int) var->data.size(2)==(mesh2->nz+1) && var->data.size(3)==var->dim;
} else if ( var->type == IO::EdgeVariable ) {
ERROR("Not finished");
} else if ( var->type == IO::SurfaceVariable ) {
ERROR("Not finished");
} else if ( var->type == IO::VolumeVariable ) {
pass = pass && (int) var->data.size(0)==mesh2->nx && (int) var->data.size(1)==mesh2->ny
&& (int) var->data.size(2)==mesh2->nz && var->data.size(3)==var->dim;
} else {
ERROR("Invalid variable type");
}
}
} else {
ERROR("Unknown mesh class: "+mesh->className());
}
return pass;
}
/**************************************************** /****************************************************
* PointList * * PointList *
****************************************************/ ****************************************************/

12
IO/Mesh.h
View File

@ -64,6 +64,8 @@ public:
virtual std::pair<size_t,void*> pack( int level ) const; virtual std::pair<size_t,void*> pack( int level ) const;
//! Unpack the data //! Unpack the data
virtual void unpack( const std::pair<size_t,void*>& data ); virtual void unpack( const std::pair<size_t,void*>& data );
//! Access the points
const std::vector<Point>& getPoints() const { return points; }
public: public:
std::vector<Point> points; //!< List of points vertex std::vector<Point> points; //!< List of points vertex
}; };
@ -171,7 +173,13 @@ public:
std::string name; //!< Variable name std::string name; //!< Variable name
Array<double> data; //!< Variable data Array<double> data; //!< Variable data
//! Empty constructor //! Empty constructor
Variable(): type(NullVariable) {} Variable(): dim(0), type(NullVariable) {}
//! Constructor
Variable( int dim_, IO::VariableType type_, const std::string& name_ ):
dim(dim_), type(type_), name(name_) {}
//! Constructor
Variable( int dim_, IO::VariableType type_, const std::string& name_, const Array<double>& data_ ):
dim(dim_), type(type_), name(name_), data(data_) {}
//! Destructor //! Destructor
virtual ~Variable() {} virtual ~Variable() {}
protected: protected:
@ -189,6 +197,8 @@ struct MeshDataStruct {
std::string meshName; std::string meshName;
std::shared_ptr<Mesh> mesh; std::shared_ptr<Mesh> mesh;
std::vector<std::shared_ptr<Variable> > vars; std::vector<std::shared_ptr<Variable> > vars;
//! Check the data
bool check() const;
}; };

38
IO/Reader.cpp
View File

@ -4,6 +4,11 @@
#include "IO/IOHelpers.h" #include "IO/IOHelpers.h"
#include "common/Utilities.h" #include "common/Utilities.h"
#ifdef USE_SILO
#include "IO/silo.h"
#endif
#include <ProfilerApp.h> #include <ProfilerApp.h>
#include <iostream> #include <iostream>
#include <string.h> #include <string.h>
@ -44,6 +49,9 @@ std::vector<std::string> IO::readTimesteps( const std::string& filename )
while (fgets(buf,sizeof(buf),fid) != NULL) { while (fgets(buf,sizeof(buf),fid) != NULL) {
std::string line(buf); std::string line(buf);
line.resize(line.size()-1); line.resize(line.size()-1);
auto pos = line.find( "summary.silo" );
if ( pos != std::string::npos )
line.resize(pos);
if ( line.empty() ) if ( line.empty() )
continue; continue;
timesteps.push_back(line); timesteps.push_back(line);
@ -146,6 +154,36 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
} }
mesh->unpack( std::pair<size_t,void*>(bytes,data) ); mesh->unpack( std::pair<size_t,void*>(bytes,data) );
delete [] data; delete [] data;
} else if ( meshDatabase.format==4 ) {
// Reading a silo file
#ifdef USE_SILO
const DatabaseEntry& database = meshDatabase.domains[domain];
std::string filename = path + "/" + timestep + "/" + database.file;
int rank = std::stoi(database.file.substr(0,database.file.find(".silo")).c_str());
auto fid = silo::open( filename, silo::READ );
if ( meshDatabase.meshClass=="PointList" ) {
mesh.reset( new IO::PointList() );
ERROR("Not finished");
} else if ( meshDatabase.meshClass=="TriMesh" ) {
mesh.reset( new IO::TriMesh() );
ERROR("Not finished");
} else if ( meshDatabase.meshClass=="TriList" ) {
mesh.reset( new IO::TriList() );
ERROR("Not finished");
} else if ( meshDatabase.meshClass=="DomainMesh" ) {
std::vector<double> range;
std::vector<int> N;
silo::readUniformMesh( fid, database.name, range, N );
auto rankinfo = silo::read<int>( fid, database.name+"_rankinfo" );
RankInfoStruct rank_data( rankinfo[0], rankinfo[1], rankinfo[2], rankinfo[3] );
mesh.reset( new IO::DomainMesh( rank_data, N[0], N[1], N[2], range[1]-range[0], range[3]-range[2], range[5]-range[4] ) );
} else {
ERROR("Unknown mesh class");
}
silo::close( fid );
#else
ERROR("Build without silo support");
#endif
} else { } else {
ERROR("Unknown format"); ERROR("Unknown format");
} }

86
IO/Writer.cpp
View File

@ -199,17 +199,53 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
// Write a PointList mesh (and variables) to a file // Write a PointList mesh (and variables) to a file
static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database ) static void writeSiloPointList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{ {
ERROR("Not finished yet"); const IO::PointList& mesh = dynamic_cast<IO::PointList&>( *meshData.mesh );
const std::string meshname = database.domains[0].name;
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];
silo::writePointMeshVariable( fid, meshname, var.name, var.data );
}
} }
// Write a TriMesh mesh (and variables) to a file // Write a TriMesh mesh (and variables) to a file
static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct& meshData,
const IO::TriMesh& mesh, IO::MeshDatabase database )
{
const std::string meshname = database.domains[0].name;
const auto& points = mesh.vertices->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() };
const int *tri[] = { mesh.A.data(), mesh.B.data(), mesh.C.data() };
silo::writeTriMesh( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri );
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 );
silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, type );
printf("%s-%s: %i %i %i\n",meshname.c_str(),var.name.c_str(),points.size(),var.data.size(0),var.data.size(1));
}
}
static void writeSiloTriMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database ) static void writeSiloTriMesh( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{ {
ERROR("Not finished yet"); const IO::TriMesh& mesh = dynamic_cast<IO::TriMesh&>( *meshData.mesh );
writeSiloTriMesh2( fid, meshData, mesh, database );
} }
// Write a TriList mesh (and variables) to a file
static void writeSiloTriList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database ) static void writeSiloTriList( DBfile *fid, const IO::MeshDataStruct& meshData, IO::MeshDatabase database )
{ {
ERROR("Not finished yet"); auto mesh = getTriMesh( meshData.mesh );
writeSiloTriMesh2( fid, meshData, *mesh, database );
} }
// Write a DomainMesh mesh (and variables) to a file // 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 )
@ -223,6 +259,7 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData
std::array<int,3> baseindex = { info.ix, info.jy, info.kz }; std::array<int,3> baseindex = { info.ix, info.jy, info.kz };
const std::string meshname = database.domains[0].name; const std::string meshname = database.domains[0].name;
silo::writeUniformMesh<3>( fid, meshname, range, N ); 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++) { for (size_t i=0; i<meshData.vars.size(); i++) {
const IO::Variable& var = *meshData.vars[i]; const IO::Variable& var = *meshData.vars[i];
auto type = static_cast<silo::VariableType>( var.type ); auto type = static_cast<silo::VariableType>( var.type );
@ -251,21 +288,43 @@ static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string& filen
return database; return database;
} }
// Write the summary file for silo // Write the summary file for silo
std::pair<int,int> getSiloMeshType( const std::string& meshClass )
{
int meshType = 0;
int varType = 0;
if ( meshClass=="PointList" ) {
meshType = DB_POINTMESH;
varType = DB_POINTVAR;
} else if ( meshClass=="TriMesh" ) {
meshType = DB_UCDMESH;
varType = DB_UCDVAR;
} else if ( meshClass=="TriList" ) {
meshType = DB_UCDMESH;
varType = DB_UCDVAR;
} else if ( meshClass=="DomainMesh" ) {
meshType = DB_QUAD_RECT;
varType = DB_QUADVAR;
} else {
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 ); auto fid = silo::open( filename, silo::CREATE );
for ( const auto& data : meshes_written ) { for ( const auto& data : meshes_written ) {
std::vector<std::string> meshnames; 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 ); meshNames.push_back( tmp.file + ":" + tmp.name );
std::vector<int> meshTypes(meshnames.size(),DB_QUAD_RECT); // Not correct for all types silo::writeMultiMesh( fid, data.name, meshNames, meshTypes );
silo::writeMultiMesh( fid, data.name, meshnames, meshTypes );
for (const auto& variable : data.variables ) { for (const auto& variable : data.variables ) {
std::vector<int> varTypes(meshnames.size(),DB_QUADVAR); // Not correct for all types
std::vector<std::string> varnames; std::vector<std::string> varnames;
for ( const auto& tmp : data.domains ) for ( const auto& tmp : data.domains )
varnames.push_back( tmp.file + ":" + variable.name ); varnames.push_back( tmp.file + ":" + variable.name );
silo::writeMultiVar( fid, variable.name, varnames, varTypes, 3, variable.dim ); silo::writeMultiVar( fid, variable.name, varnames, varTypes );
} }
} }
silo::close( fid ); silo::close( fid );
@ -325,6 +384,11 @@ void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStr
IO::initialize( ); IO::initialize( );
PROFILE_START("writeData"); PROFILE_START("writeData");
int rank = comm_rank(comm); int rank = comm_rank(comm);
// Check the meshData before writing
for ( const auto& data : meshData ) {
if ( !data.check() )
ERROR("Error in meshData");
}
// Create the output directory // Create the output directory
std::string path = global_IO_path + "/" + subdir; std::string path = global_IO_path + "/" + subdir;
if ( rank == 0 ) { if ( rank == 0 ) {
@ -369,7 +433,7 @@ void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStr
} else if ( global_IO_format == Format::SILO ) { } else if ( global_IO_format == Format::SILO ) {
auto filename = global_IO_path+"/LBM.visit"; auto filename = global_IO_path+"/LBM.visit";
FILE *fid = fopen(filename.c_str(),"ab"); FILE *fid = fopen(filename.c_str(),"ab");
fprintf(fid,"%s/summary.silo\n",path.c_str()); fprintf(fid,"%s/summary.silo\n",subdir.c_str());
fclose(fid); fclose(fid);
} else { } else {
ERROR("Unknown format"); ERROR("Unknown format");

242
IO/silo.cpp
View File

@ -36,48 +36,34 @@ void close( DBfile* fid )
/**************************************************** /****************************************************
* Write a uniform mesh to silo * * Write/read an arbitrary vector *
****************************************************/ ****************************************************/
template<int NDIM> template<class TYPE> int getSiloType();
void writeUniformMesh( DBfile* fid, const std::string& meshname, template<> int getSiloType<int>() { return DB_INT; }
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N ) template<> int getSiloType<float>() { return DB_FLOAT; }
template<> int getSiloType<double>() { return DB_DOUBLE; }
template<class TYPE>
void write( DBfile* fid, const std::string& varname, const std::vector<TYPE>& data )
{ {
PROFILE_START("writeUniformMesh",2); int dims = data.size();
int dims[NDIM]; int err = DBWrite( fid, varname.c_str(), data.data(), &dims, 1, getSiloType<TYPE>() );
for (size_t d=0; d<N.size(); d++)
dims[d] = N[d]+1;
float *x = nullptr;
if ( NDIM >= 1 ) {
x = new float[dims[0]];
for (int i=0; i<N[0]; i++)
x[i] = range[0] + i*(range[1]-range[0])/N[0];
x[N[0]] = range[1];
}
float *y = nullptr;
if ( NDIM >= 2 ) {
y = new float[dims[1]];
for (int i=0; i<N[1]; i++)
y[i] = range[2] + i*(range[3]-range[2])/N[1];
y[N[1]] = range[3];
}
float *z = nullptr;
if ( NDIM >= 3 ) {
z = new float[dims[2]];
for (int i=0; i<N[2]; i++)
z[i] = range[4] + i*(range[5]-range[4])/N[2];
z[N[2]] = range[5];
}
float *coords[] = { x, y, z };
int err = DBPutQuadmesh( fid, meshname.c_str(), nullptr, coords, dims, NDIM, DB_FLOAT, DB_COLLINEAR, nullptr );
ASSERT( err == 0 ); ASSERT( err == 0 );
PROFILE_STOP("writeUniformMesh",2); }
template<class TYPE>
std::vector<TYPE> read( DBfile* fid, const std::string& varname )
{
int N = DBGetVarLength( fid, varname.c_str() );
std::vector<TYPE> data(N);
int err = DBReadVar( fid, varname.c_str(), data.data() );
ASSERT( err == 0 );
return data;
} }
/**************************************************** /****************************************************
* Write a vector/tensor quad variable * * Helper function to get variable suffixes *
****************************************************/ ****************************************************/
std::vector<std::string> getVarSuffix( int ndim, int nvars ) static std::vector<std::string> getVarSuffix( int ndim, int nvars )
{ {
std::vector<std::string> suffix(nvars); std::vector<std::string> suffix(nvars);
if ( nvars == 1 ) { if ( nvars == 1 ) {
@ -118,24 +104,93 @@ std::vector<std::string> getVarSuffix( int ndim, int nvars )
} }
return suffix; return suffix;
} }
/****************************************************
* Write/read a uniform mesh to silo *
****************************************************/
template<int NDIM>
void writeUniformMesh( DBfile* fid, const std::string& meshname,
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N )
{
PROFILE_START("writeUniformMesh",2);
int dims[NDIM];
for (size_t d=0; d<N.size(); d++)
dims[d] = N[d]+1;
float *x = nullptr;
if ( NDIM >= 1 ) {
x = new float[dims[0]];
for (int i=0; i<N[0]; i++)
x[i] = range[0] + i*(range[1]-range[0])/N[0];
x[N[0]] = range[1];
}
float *y = nullptr;
if ( NDIM >= 2 ) {
y = new float[dims[1]];
for (int i=0; i<N[1]; i++)
y[i] = range[2] + i*(range[3]-range[2])/N[1];
y[N[1]] = range[3];
}
float *z = nullptr;
if ( NDIM >= 3 ) {
z = new float[dims[2]];
for (int i=0; i<N[2]; i++)
z[i] = range[4] + i*(range[5]-range[4])/N[2];
z[N[2]] = range[5];
}
float *coords[] = { x, y, z };
int err = DBPutQuadmesh( fid, meshname.c_str(), nullptr, coords, dims, NDIM, DB_FLOAT, DB_COLLINEAR, nullptr );
ASSERT( err == 0 );
PROFILE_STOP("writeUniformMesh",2);
}
void readUniformMesh( DBfile* fid, const std::string& meshname,
std::vector<double>& range, std::vector<int>& N )
{
DBquadmesh* mesh = DBGetQuadmesh( fid, meshname.c_str() );
int ndim = mesh->ndims;
range.resize(2*ndim);
N.resize(ndim);
for (int d=0; d<ndim; d++) {
N[d] = mesh->dims[d];
range[2*d+0] = mesh->min_extents[d];
range[2*d+1] = mesh->max_extents[d];
}
delete mesh;
}
/****************************************************
* Write a vector/tensor quad variable *
****************************************************/
template<int NDIM> template<int NDIM>
void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N, void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
const std::string& varname, const Array<double>& data, VariableType type ) const std::string& varname, const Array<double>& data, VariableType type )
{ {
for (int d=0; d<NDIM; d++)
ASSERT(N[d]==(int)data.size(d));
PROFILE_START("writeUniformMeshVariable",2); PROFILE_START("writeUniformMeshVariable",2);
int nvars=1, dims[NDIM]={1}; int nvars=1, dims[NDIM]={1};
const double *vars[NDIM] = { nullptr }; const double *vars[NDIM] = { nullptr };
int vartype = 0;
if ( type == NodeVariable ) { if ( type == NodeVariable ) {
ERROR("Not finished"); ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
for (int d=0; d<NDIM; d++)
ASSERT(N[d]+1==(int)data.size(d));
vartype = DB_NODECENT;
nvars = data.size(NDIM);
size_t N = data.length()/nvars;
for (int d=0; d<NDIM; d++)
dims[d] = data.size(d);
for (int i=0; i<nvars; i++)
vars[i] = &data(i*N);
} else if ( type == EdgeVariable ) { } else if ( type == EdgeVariable ) {
ERROR("Not finished"); ERROR("Not finished");
} else if ( type == SurfaceVariable ) { } else if ( type == SurfaceVariable ) {
ERROR("Not finished"); ERROR("Not finished");
} else if ( type == VolumeVariable ) { } else if ( type == VolumeVariable ) {
ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 ); ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
nvars = data.size(NDIM+1); for (int d=0; d<NDIM; d++)
ASSERT(N[d]==(int)data.size(d));
vartype = DB_ZONECENT;
nvars = data.size(NDIM);
size_t N = data.length()/nvars; size_t N = data.length()/nvars;
for (int d=0; d<NDIM; d++) for (int d=0; d<NDIM; d++)
dims[d] = data.size(d); dims[d] = data.size(d);
@ -152,12 +207,97 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
for (int i=0; i<nvars; i++) for (int i=0; i<nvars; i++)
varnames[i] = const_cast<char*>(var_names[i].c_str()); varnames[i] = const_cast<char*>(var_names[i].c_str());
int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars, int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars,
varnames.data(), vars, dims, NDIM, nullptr, 0, DB_DOUBLE, DB_ZONECENT, nullptr ); varnames.data(), vars, dims, NDIM, nullptr, 0, DB_DOUBLE, vartype, nullptr );
ASSERT( err == 0 ); ASSERT( err == 0 );
PROFILE_STOP("writeUniformMeshVariable",2); PROFILE_STOP("writeUniformMeshVariable",2);
} }
/****************************************************
* Write a point mesh/variable to silo *
****************************************************/
void writePointMesh( DBfile* fid, const std::string& meshname,
int ndim, int N, const double *coords[] )
{
int err = DBPutPointmesh( fid, meshname.c_str(), ndim, coords, N, DB_DOUBLE, nullptr );
ASSERT( err == 0 );
}
void writePointMeshVariable( DBfile* fid, const std::string& meshname,
const std::string& varname, const Array<double>& data )
{
int N = data.size(0);
int nvars = data.size(1);
std::vector<const double*> vars(nvars);
for (int i=0; i<nvars; i++)
vars[i] = &data(0,i);
int err = DBPutPointvar( fid, varname.c_str(), meshname.c_str(), nvars, vars.data(), N, DB_DOUBLE, nullptr );
ASSERT( err == 0 );
}
/****************************************************
* Write a triangle mesh *
****************************************************/
void writeTriMesh( DBfile* fid, const std::string& meshName,
int ndim, int ndim_tri, int N, const double *coords[], int N_tri, const int *tri[] )
{
auto zoneName = meshName + "_zones";
std::vector<int> nodelist( (ndim_tri+1)*N_tri );
for (int i=0, j=0; i<N_tri; i++) {
for (int d=0; d<ndim_tri+1; d++, j++)
nodelist[j] = tri[d][i];
}
int shapetype = 0;
if ( ndim_tri==1 )
shapetype = DB_ZONETYPE_BEAM;
else if ( ndim_tri==2 )
shapetype = DB_ZONETYPE_TRIANGLE;
else if ( ndim_tri==3 )
shapetype = DB_ZONETYPE_PYRAMID;
else
ERROR("Unknown shapetype");
int shapesize = ndim_tri+1;
int shapecnt = N_tri;
DBPutZonelist2( fid, zoneName.c_str(), N_tri, ndim_tri, nodelist.data(),
nodelist.size(), 0, 0, 0, &shapetype, &shapesize, &shapecnt, 1, nullptr );
DBPutUcdmesh( fid, meshName.c_str(), ndim, nullptr, coords, N,
nodelist.size(), zoneName.c_str(), nullptr, DB_DOUBLE, nullptr );
}
void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
const std::string& varname, const Array<double>& data, VariableType type )
{
int nvars = 0;
int vartype = 0;
const double *vars[10] = { nullptr };
if ( type == NodeVariable ) {
vartype = DB_NODECENT;
nvars = data.size(1);
for (int i=0; i<nvars; i++)
vars[i] = &data(0,i);
} else if ( type == EdgeVariable ) {
ERROR("Not finished");
} else if ( type == SurfaceVariable ) {
ERROR("Not finished");
} else if ( type == VolumeVariable ) {
vartype = DB_ZONECENT;
nvars = data.size(1);
for (int i=0; i<nvars; i++)
vars[i] = &data(0,i);
} else {
ERROR("Invalid variable type");
}
auto suffix = getVarSuffix( ndim, nvars );
std::vector<std::string> var_names(nvars);
for (int i=0; i<nvars; i++)
var_names[i] = varname + suffix[i];
std::vector<char*> varnames(nvars,nullptr);
for (int i=0; i<nvars; i++)
varnames[i] = const_cast<char*>(var_names[i].c_str());
DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars,
varnames.data(), vars, data.size(0), nullptr, 0, DB_DOUBLE, vartype, nullptr );
}
/**************************************************** /****************************************************
* Write a multimesh * * Write a multimesh *
****************************************************/ ****************************************************/
@ -181,20 +321,12 @@ void writeMultiMesh( DBfile* fid, const std::string& meshname,
****************************************************/ ****************************************************/
void writeMultiVar( DBfile* fid, const std::string& varname, void writeMultiVar( DBfile* fid, const std::string& varname,
const std::vector<std::string>& varNames, const std::vector<std::string>& varNames,
const std::vector<int>& varTypes, int ndim, int nvar ) const std::vector<int>& varTypes )
{ {
ASSERT(varNames.size()==varTypes.size()); std::vector<char*> varnames(varNames.size(),nullptr);
auto suffix = getVarSuffix( ndim, nvar ); for (size_t j=0; j<varNames.size(); j++)
for (int i=0; i<nvar; i++) { varnames[j] = const_cast<char*>(varNames[j].c_str());
auto varname2 = varname + suffix[i]; DBPutMultivar( fid, varname.c_str(), varNames.size(), varnames.data(), (int*) varTypes.data(), nullptr );
auto varNames2 = varNames;
for ( auto& tmp : varNames2 )
tmp += suffix[i];
std::vector<char*> varnames(varNames.size(),nullptr);
for (size_t j=0; j<varNames.size(); j++)
varnames[j] = const_cast<char*>(varNames2[j].c_str());
DBPutMultivar( fid, varname2.c_str(), varNames.size(), varnames.data(), (int*) varTypes.data(), nullptr );
}
} }
@ -202,6 +334,12 @@ void writeMultiVar( DBfile* fid, const std::string& varname,
// Explicit instantiations // Explicit instantiations
template void silo::write<int>( DBfile*, const std::string&, const std::vector<int>& );
template void silo::write<float>( DBfile*, const std::string&, const std::vector<float>& );
template void silo::write<double>( DBfile*, const std::string&, const std::vector<double>& );
template std::vector<int> silo::read<int>( DBfile* fid, const std::string& varname );
template std::vector<float> silo::read<float>( DBfile* fid, const std::string& varname );
template std::vector<double> silo::read<double>( DBfile* fid, const std::string& varname );
template void silo::writeUniformMesh<1>( DBfile*, const std::string&, const std::array<double,2>&, const std::array<int,1>& ); template void silo::writeUniformMesh<1>( DBfile*, const std::string&, const std::array<double,2>&, const std::array<int,1>& );
template void silo::writeUniformMesh<2>( DBfile*, const std::string&, const std::array<double,4>&, const std::array<int,2>& ); template void silo::writeUniformMesh<2>( DBfile*, const std::string&, const std::array<double,4>&, const std::array<int,2>& );
template void silo::writeUniformMesh<3>( DBfile*, const std::string&, const std::array<double,6>&, const std::array<int,3>& ); template void silo::writeUniformMesh<3>( DBfile*, const std::string&, const std::array<double,6>&, const std::array<int,3>& );

110
IO/silo.h
View File

@ -44,10 +44,32 @@ DBfile* open( const std::string& filename, FileMode mode );
void close( DBfile* fid ); void close( DBfile* fid );
/*!
* @brief Write data to silo
* @detailed This function writes an arbitrary array to silo
* @param[in] fid Handle to the open file
* @param[in] varname Variable name
* @param[in] data Data to write
*/
template<class TYPE>
void write( DBfile* fid, const std::string& varname, const std::vector<TYPE>& data );
/*!
* @brief Write data to silo
* @detailed This function writes an arbitrary array to silo
* @param[in] fid Handle to the open file
* @param[in] varname Variable name
* @return Data read
*/
template<class TYPE>
std::vector<TYPE> read( DBfile* fid, const std::string& varname );
/*! /*!
* @brief Write a uniform grid * @brief Write a uniform grid
* @detailed This function writes a uniform grid to silo as a Quadmesh * @detailed This function writes a uniform grid to silo as a Quadmesh
* @return This function returns a handle to the file * @param[in] fid Handle to the open file
* @param[in] meshname Mesh name * @param[in] meshname Mesh name
* @param[in] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax } * @param[in] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax }
* @param[in] N Number of cells in each direction * @param[in] N Number of cells in each direction
@ -58,21 +80,20 @@ void writeUniformMesh( DBfile* fid, const std::string& meshname,
/*! /*!
* @brief Write a multimesh * @brief Read a uniform grid
* @detailed This function writes a multimesh to silo * @detailed This function reads a uniform grid from silo
* @return This function returns a handle to the file * @param[in] fid Handle to the open file
* @param[in] meshname Mesh name * @param[in] meshname Mesh name
* @param[in] subMeshNames Names of the sub meshes in the form "filename:meshname" * @param[out] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax }
* @param[in] subMeshTypes Type of each submesh * @param[out] N Number of cells in each direction
*/ */
void writeMultiMesh( DBfile* fid, const std::string& meshname, void readUniformMesh( DBfile* fid, const std::string& meshname,
const std::vector<std::string>& subMeshNames, std::vector<double>& range, std::vector<int>& N );
const std::vector<int>& subMeshTypes );
/*! /*!
* @brief Write a uniform grid * @brief Write a uniform grid variable
* @detailed This function writes a uniform grid to silo as a Quadmesh * @detailed This function writes a uniform grid variable to silo as a Quadmesh
* @param[in] fid Handle to the open file * @param[in] fid Handle to the open file
* @param[in] meshname Mesh name * @param[in] meshname Mesh name
* @param[in] N Number of cells in each direction * @param[in] N Number of cells in each direction
@ -85,6 +106,71 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
const std::string& varname, const Array<double>& data, VariableType type ); const std::string& varname, const Array<double>& data, VariableType type );
/*!
* @brief Write a pointmesh
* @detailed This function writes a pointmesh to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] ndim Number of dimensions
* @param[in] N Number of points
* @param[in] coords Coordinates of the points
*/
void writePointMesh( DBfile* fid, const std::string& meshname,
int ndim, int N, const double *coords[] );
/*!
* @brief Write a pointmesh grid variable
* @detailed This function writes a pointmesh variable to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] varname Variable name
* @param[in] data Variable data
*/
void writePointMeshVariable( DBfile* fid, const std::string& meshname,
const std::string& varname, const Array<double>& data );
/*!
* @brief Write a triangle mesh
* @detailed This function writes a triangle (or simplex) based mesh to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] ndim Number of dimensions for the coordinates
* @param[in] ndim_tri Number of dimensions for the triangles (2: surface, 3: volume)
* @param[in] N Number of points
* @param[in] coords Coordinates of the points
* @param[in] N_tri Number of triangles
* @param[in] tri Coordinates of the points
*/
void writeTriMesh( DBfile* fid, const std::string& meshname,
int ndim, int ndim_tri, int N, const double *coords[], int N_tri, const int *tri[] );
/*!
* @brief Write a pointmesh grid variable
* @detailed This function writes a pointmesh variable to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] varname Variable name
* @param[in] data Variable data
*/
void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
const std::string& varname, const Array<double>& data, VariableType type );
/*!
* @brief Write a multimesh
* @detailed This function writes a multimesh to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] subMeshNames Names of the sub meshes in the form "filename:meshname"
* @param[in] subMeshTypes Type of each submesh
*/
void writeMultiMesh( DBfile* fid, const std::string& meshname,
const std::vector<std::string>& subMeshNames,
const std::vector<int>& subMeshTypes );
/*! /*!
* @brief Write a multivariable * @brief Write a multivariable
* @detailed This function writes a multivariable to silo * @detailed This function writes a multivariable to silo
@ -98,7 +184,7 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
*/ */
void writeMultiVar( DBfile* fid, const std::string& varname, void writeMultiVar( DBfile* fid, const std::string& varname,
const std::vector<std::string>& subVarNames, const std::vector<std::string>& subVarNames,
const std::vector<int>& subVarTypes, int ndim, int nvar ); const std::vector<int>& subVarTypes );
}; // silo namespace }; // silo namespace
#endif #endif

109
tests/TestWriter.cpp
View File

@ -257,32 +257,77 @@ int main(int argc, char **argv)
std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) ); std::shared_ptr<IO::DomainMesh> domain( new IO::DomainMesh(rank_data,6,7,8,1.0,1.0,1.0) );
// Create the variables // Create the variables
std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() ); std::shared_ptr<IO::Variable> set_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_set_mag") );
std::shared_ptr<IO::Variable> dist_list( new IO::Variable() ); std::shared_ptr<IO::Variable> set_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_set_vec") );
std::shared_ptr<IO::Variable> domain_var( new IO::Variable() ); std::shared_ptr<IO::Variable> list_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_list_mag") );
dist_set1->dim = 1; std::shared_ptr<IO::Variable> list_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_list_vec") );
dist_list->dim = 1; std::shared_ptr<IO::Variable> point_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_point_mag") );
domain_var->dim = 1; std::shared_ptr<IO::Variable> point_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_point_vec") );
dist_set1->name = "Distance"; std::shared_ptr<IO::Variable> domain_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_domain_mag") );
dist_list->name = "Distance"; std::shared_ptr<IO::Variable> domain_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_domain_vec") );
domain_var->name = "Distance"; std::shared_ptr<IO::Variable> set_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_set_mag") );
dist_set1->type = IO::NodeVariable; std::shared_ptr<IO::Variable> set_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_set_vec") );
dist_list->type = IO::NodeVariable; std::shared_ptr<IO::Variable> list_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_list_mag") );
domain_var->type = IO::VolumeVariable; std::shared_ptr<IO::Variable> list_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_list_vec") );
dist_set1->data.resize( N_points ); std::shared_ptr<IO::Variable> domain_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_domain_mag") );
for (int i=0; i<N_points; i++) std::shared_ptr<IO::Variable> domain_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_domain_vec") );
dist_set1->data(i) = distance(set1->points[i]); point_node_mag->data.resize( N_points );
dist_list->data.resize( 3, N_tri ); point_node_vec->data.resize( N_points, 3 );
for (int i=0; i<N_tri; i++) { for (int i=0; i<N_points; i++) {
dist_list->data(0,i) = distance(trilist->A[i]); point_node_mag->data(i) = distance(set1->points[i]);
dist_list->data(1,i) = distance(trilist->B[i]); point_node_vec->data(i,0) = set1->points[i].x;
dist_list->data(2,i) = distance(trilist->C[i]); point_node_vec->data(i,1) = set1->points[i].y;
point_node_vec->data(i,2) = set1->points[i].z;
} }
domain_var->data.resize(domain->nx,domain->ny,domain->nz); set_node_mag->data = point_node_mag->data;
set_node_vec->data = point_node_vec->data;
list_node_mag->data.resize( 3*N_tri );
list_node_vec->data.resize( 3*N_tri, 3 );
for (int i=0; i<N_points; i++) {
list_node_mag->data(3*i+0) = distance(trilist->A[i]);
list_node_mag->data(3*i+1) = distance(trilist->B[i]);
list_node_mag->data(3*i+2) = distance(trilist->C[i]);
list_node_vec->data(3*i+0,0) = trilist->A[i].x;
list_node_vec->data(3*i+0,1) = trilist->A[i].y;
list_node_vec->data(3*i+0,2) = trilist->A[i].z;
list_node_vec->data(3*i+1,0) = trilist->B[i].x;
list_node_vec->data(3*i+1,1) = trilist->B[i].y;
list_node_vec->data(3*i+1,2) = trilist->B[i].z;
list_node_vec->data(3*i+2,0) = trilist->C[i].x;
list_node_vec->data(3*i+2,1) = trilist->C[i].y;
list_node_vec->data(3*i+2,2) = trilist->C[i].z;
}
domain_node_mag->data.resize(domain->nx+1,domain->ny+1,domain->nz+1);
domain_node_vec->data.resize({(size_t)domain->nx+1,(size_t)domain->ny+1,(size_t)domain->nz+1,3});
for (int i=0; i<domain->nx+1; i++) {
for (int j=0; j<domain->ny+1; j++) {
for (int k=0; k<domain->nz+1; k++) {
domain_node_mag->data(i,j,k) = distance(Point(i,j,k));
domain_node_vec->data(i,j,k,0) = Point(i,j,k).x;
domain_node_vec->data(i,j,k,1) = Point(i,j,k).y;
domain_node_vec->data(i,j,k,2) = Point(i,j,k).z;
}
}
}
set_cell_mag->data.resize( N_tri );
set_cell_vec->data.resize( N_tri, 3 );
for (int i=0; i<N_tri; i++) {
set_cell_mag->data(i) = i;
set_cell_vec->data(i,0) = 3*i+0;
set_cell_vec->data(i,1) = 3*i+1;
set_cell_vec->data(i,2) = 3*i+2;
}
list_cell_mag->data = set_cell_mag->data;
list_cell_vec->data = set_cell_vec->data;
domain_cell_mag->data.resize(domain->nx,domain->ny,domain->nz);
domain_cell_vec->data.resize({(size_t)domain->nx,(size_t)domain->ny,(size_t)domain->nz,3});
for (int i=0; i<domain->nx; i++) { for (int i=0; i<domain->nx; i++) {
for (int j=0; j<domain->ny; j++) { for (int j=0; j<domain->ny; j++) {
for (int k=0; k<domain->nz; k++) { for (int k=0; k<domain->nz; k++) {
domain_var->data(i,j,k) = distance(Point(i,j,k)); domain_cell_mag->data(i,j,k) = distance(Point(i,j,k));
domain_cell_vec->data(i,j,k,0) = Point(i,j,k).x;
domain_cell_vec->data(i,j,k,1) = Point(i,j,k).y;
domain_cell_vec->data(i,j,k,2) = Point(i,j,k).z;
} }
} }
} }
@ -291,21 +336,31 @@ int main(int argc, char **argv)
std::vector<IO::MeshDataStruct> meshData(4); std::vector<IO::MeshDataStruct> meshData(4);
meshData[0].meshName = "pointmesh"; meshData[0].meshName = "pointmesh";
meshData[0].mesh = set1; meshData[0].mesh = set1;
meshData[0].vars.push_back(dist_set1); meshData[0].vars.push_back(point_node_mag);
meshData[0].vars.push_back(point_node_vec);
meshData[1].meshName = "trimesh"; meshData[1].meshName = "trimesh";
meshData[1].mesh = trimesh; meshData[1].mesh = trimesh;
meshData[1].vars.push_back(dist_set1); meshData[1].vars.push_back(set_node_mag);
meshData[1].vars.push_back(set_node_vec);
meshData[1].vars.push_back(set_cell_mag);
meshData[1].vars.push_back(set_cell_vec);
meshData[2].meshName = "trilist"; meshData[2].meshName = "trilist";
meshData[2].mesh = trilist; meshData[2].mesh = trilist;
meshData[2].vars.push_back(dist_list); meshData[2].vars.push_back(list_node_mag);
meshData[2].vars.push_back(list_node_vec);
meshData[2].vars.push_back(list_cell_mag);
meshData[2].vars.push_back(list_cell_vec);
meshData[3].meshName = "domain"; meshData[3].meshName = "domain";
meshData[3].mesh = domain; meshData[3].mesh = domain;
meshData[3].vars.push_back(domain_var); meshData[3].vars.push_back(domain_node_mag);
meshData[3].vars.push_back(domain_node_vec);
meshData[3].vars.push_back(domain_cell_mag);
meshData[3].vars.push_back(domain_cell_vec);
// Run the tests // Run the tests
testWriter( "old", meshData, ut ); testWriter( "old", meshData, ut );
testWriter( "new", meshData, ut ); testWriter( "new", meshData, ut );
//testWriter( "silo", meshData, ut ); testWriter( "silo", meshData, ut );
// Finished // Finished
ut.report(); ut.report();