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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure
2017-02-15 16:38:51 -05:00
parent 4acd3efafb 3c590b56b4
commit 1227aaf482
10 changed files with 860 additions and 166 deletions
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 *
****************************************************/

12
IO/Mesh.h
View File

@@ -64,6 +64,8 @@ public:
virtual std::pair<size_t,void*> pack( int level ) const;
//! Unpack the data
virtual void unpack( const std::pair<size_t,void*>& data );
//! Access the points
const std::vector<Point>& getPoints() const { return points; }
public:
std::vector<Point> points; //!< List of points vertex
};
@@ -171,7 +173,13 @@ public:
std::string name; //!< Variable name
Array<double> data; //!< Variable data
//! 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
virtual ~Variable() {}
protected:
@@ -189,6 +197,8 @@ struct MeshDataStruct {
std::string meshName;
std::shared_ptr<Mesh> mesh;
std::vector<std::shared_ptr<Variable> > vars;
//! Check the data
bool check() const;
};

8
IO/MeshDatabase.cpp
View File

@@ -203,6 +203,14 @@ MeshDatabase& MeshDatabase::operator=(const MeshDatabase& rhs)
this->variable_data = rhs.variable_data;
return *this;
}
VariableDatabase MeshDatabase::getVariableDatabase( const std::string& varname ) const
{
for (size_t i=0; i<variables.size(); i++) {
if ( variables[i].name == varname )
return variables[i];
}
return VariableDatabase();
}
/****************************************************

1
IO/MeshDatabase.h
View File

@@ -60,6 +60,7 @@ struct MeshDatabase {
std::vector<DatabaseEntry> domains; //!< List of the domains
std::vector<VariableDatabase> variables; //!< List of the variables
std::map<variable_id,DatabaseEntry> variable_data; //!< Data for the variables
VariableDatabase getVariableDatabase( const std::string& varname ) const;
public:
MeshDatabase();
~MeshDatabase();

179
IO/Reader.cpp
View File

@@ -4,6 +4,11 @@
#include "IO/IOHelpers.h"
#include "common/Utilities.h"
#ifdef USE_SILO
#include "IO/silo.h"
#endif
#include <ProfilerApp.h>
#include <iostream>
#include <string.h>
@@ -39,14 +44,21 @@ std::vector<std::string> IO::readTimesteps( const std::string& filename )
FILE *fid= fopen(filename.c_str(),"rb");
if ( fid==NULL )
ERROR("Error opening file");
auto pos = std::min(filename.find_last_of(47),filename.find_last_of(90));
std::string path = getPath( filename );
if ( !path.empty() )
path += "/";
std::vector<std::string> timesteps;
char buf[1000];
while (fgets(buf,sizeof(buf),fid) != NULL) {
std::string line(buf);
line.resize(line.size()-1);
auto pos = line.find( "summary.silo" );
if ( pos != std::string::npos )
line.resize(pos);
if ( line.empty() )
continue;
timesteps.push_back(line);
timesteps.push_back(path+line);
}
fclose(fid);
PROFILE_STOP("readTimesteps");
@@ -146,6 +158,61 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
}
mesh->unpack( std::pair<size_t,void*>(bytes,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" ) {
Array<double> coords = silo::readPointMesh( fid, database.name );
ASSERT(coords.size(1)==3);
std::shared_ptr<IO::PointList> mesh2( new IO::PointList( coords.size(0) ) );
for (size_t i=0; i<coords.size(1); i++) {
mesh2->points[i].x = coords(i,0);
mesh2->points[i].y = coords(i,1);
mesh2->points[i].z = coords(i,2);
}
mesh = mesh2;
} else if ( meshDatabase.meshClass=="TriMesh" || meshDatabase.meshClass=="TriList" ) {
Array<double> coords;
Array<int> tri;
silo::readTriMesh( fid, database.name, coords, tri );
ASSERT( tri.size(1)==3 && coords.size(1)==3 );
int N_tri = tri.size(0);
int N_point = coords.size(0);
std::shared_ptr<IO::TriMesh> mesh2( new IO::TriMesh( N_tri, N_point ) );
for (int i=0; i<N_point; i++) {
mesh2->vertices->points[i].x = coords(i,0);
mesh2->vertices->points[i].y = coords(i,1);
mesh2->vertices->points[i].z = coords(i,2);
}
for (int i=0; i<N_tri; i++) {
mesh2->A[i] = tri(i,0);
mesh2->B[i] = tri(i,1);
mesh2->C[i] = tri(i,2);
}
if ( meshDatabase.meshClass=="TriMesh" ) {
mesh = mesh2;
} else if ( meshDatabase.meshClass=="TriList" ) {
auto trilist = IO::getTriList( std::dynamic_pointer_cast<IO::Mesh>( mesh2 ) );
mesh = trilist;
}
} 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 {
ERROR("Unknown format");
}
@@ -163,33 +230,61 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
it = meshDatabase.variable_data.find(key);
if ( it==meshDatabase.variable_data.end() )
return std::shared_ptr<IO::Variable>();
const DatabaseEntry& database = it->second;
std::string filename = path + "/" + timestep + "/" + database.file;
FILE *fid = fopen(filename.c_str(),"rb");
fseek(fid,database.offset,SEEK_SET);
char line[1000];
fgetl(line,1000,fid);
size_t i1 = find(line,':');
size_t i2 = find(&line[i1+1],':')+i1+1;
std::vector<std::string> values = splitList(&line[i2+1],',');
ASSERT(values.size()==5);
int dim = atoi(values[0].c_str());
int type = atoi(values[1].c_str());
size_t N = atol(values[2].c_str());
size_t bytes = atol(values[3].c_str());
std::string precision = values[4];
std::shared_ptr<IO::Variable> var( new IO::Variable() );
var->dim = dim;
var->type = static_cast<IO::VariableType>(type);
var->name = variable;
var->data.resize(N*dim);
if ( precision=="double" ) {
size_t count = fread(var->data.data(),sizeof(double),N*dim,fid);
ASSERT(count*sizeof(double)==bytes);
std::shared_ptr<IO::Variable> var;
if ( meshDatabase.format == 2 ) {
const DatabaseEntry& database = it->second;
std::string filename = path + "/" + timestep + "/" + database.file;
FILE *fid = fopen(filename.c_str(),"rb");
fseek(fid,database.offset,SEEK_SET);
char line[1000];
fgetl(line,1000,fid);
size_t i1 = find(line,':');
size_t i2 = find(&line[i1+1],':')+i1+1;
std::vector<std::string> values = splitList(&line[i2+1],',');
ASSERT(values.size()==5);
int dim = atoi(values[0].c_str());
int type = atoi(values[1].c_str());
size_t N = atol(values[2].c_str());
size_t bytes = atol(values[3].c_str());
std::string precision = values[4];
var = std::shared_ptr<IO::Variable>( new IO::Variable() );
var->dim = dim;
var->type = static_cast<IO::VariableType>(type);
var->name = variable;
var->data.resize(N*dim);
if ( precision=="double" ) {
size_t count = fread(var->data.data(),sizeof(double),N*dim,fid);
ASSERT(count*sizeof(double)==bytes);
} else {
ERROR("Format not implimented");
}
fclose(fid);
} else if ( meshDatabase.format == 4 ) {
// Reading a silo file
#ifdef USE_SILO
const auto& database = meshDatabase.domains[domain];
auto variableDatabase = meshDatabase.getVariableDatabase( variable );
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 );
var.reset( new Variable( variableDatabase.dim, variableDatabase.type, variable ) );
if ( meshDatabase.meshClass=="PointList" ) {
var->data = silo::readPointMeshVariable( fid, variable );
} else if ( meshDatabase.meshClass=="TriMesh" || meshDatabase.meshClass=="TriList" ) {
var->data = silo::readTriMeshVariable( fid, variable );
} else if ( meshDatabase.meshClass=="DomainMesh" ) {
var->data = silo::readUniformMeshVariable( fid, variable );
} else {
ERROR("Unknown mesh class");
}
silo::close( fid );
#else
ERROR("Build without silo support");
#endif
} else {
ERROR("Format not implimented");
ERROR("Unknown format");
}
fclose(fid);
return var;
}
@@ -202,7 +297,9 @@ void IO::reformatVariable( const IO::Mesh& mesh, IO::Variable& var )
if ( mesh.className() == "DomainMesh" ) {
const IO::DomainMesh& mesh2 = dynamic_cast<const IO::DomainMesh&>( mesh );
if ( var.type == NodeVariable ) {
ERROR("Not finished");
size_t N2 = var.data.length() / ((mesh2.nx+1)*(mesh2.ny+1)*(mesh2.nz+1));
ASSERT( (mesh2.nx+1)*(mesh2.ny+1)*(mesh2.nz+1)*N2 == var.data.length() );
var.data.reshape( { (size_t) mesh2.nx+1, (size_t) mesh2.ny+1, (size_t) mesh2.nz+1, N2 } );
} else if ( var.type == EdgeVariable ) {
ERROR("Not finished");
} else if ( var.type == SurfaceVariable ) {
@@ -215,9 +312,31 @@ void IO::reformatVariable( const IO::Mesh& mesh, IO::Variable& var )
ERROR("Invalid variable type");
}
} else if ( mesh.className() == "PointList" ) {
ERROR("Not finished");
} else if ( mesh.className() == "TriMesh" ) {
ERROR("Not finished");
const IO::PointList& mesh2 = dynamic_cast<const IO::PointList&>( mesh );
size_t N = mesh2.points.size();
size_t N_var = var.data.length()/N;
ASSERT( N*N_var == var.data.length() );
var.data.reshape( { N, N_var } );
} else if ( mesh.className()=="TriMesh" || mesh.className() == "TriList" ) {
std::shared_ptr<Mesh> mesh_ptr( const_cast<Mesh*>(&mesh), []( void* ) {} );
std::shared_ptr<TriMesh> mesh2 = getTriMesh( mesh_ptr );
if ( var.type == NodeVariable ) {
size_t N = mesh2->vertices->points.size();
size_t N_var = var.data.length()/N;
ASSERT( N*N_var == var.data.length() );
var.data.reshape( { N, N_var } );
} else if ( var.type == EdgeVariable ) {
ERROR("Not finished");
} else if ( var.type == SurfaceVariable ) {
ERROR("Not finished");
} else if ( var.type == VolumeVariable ) {
size_t N = mesh2->A.size();
size_t N_var = var.data.length()/N;
ASSERT( N*N_var == var.data.length() );
var.data.reshape( { N, N_var } );
} else {
ERROR("Invalid variable type");
}
} else {
ERROR("Unknown mesh type");
}

87
IO/Writer.cpp
View File

@@ -199,17 +199,52 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
// Write a PointList mesh (and variables) to a file
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
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 );
}
}
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 )
{
ERROR("Not finished yet");
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 )
@@ -223,6 +258,7 @@ static void writeSiloDomainMesh( DBfile *fid, const IO::MeshDataStruct& meshData
std::array<int,3> baseindex = { info.ix, info.jy, info.kz };
const std::string 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 IO::Variable& var = *meshData.vars[i];
auto type = static_cast<silo::VariableType>( var.type );
@@ -251,21 +287,43 @@ static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string& filen
return database;
}
// 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 )
{
auto fid = silo::open( filename, silo::CREATE );
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 )
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 );
meshNames.push_back( tmp.file + ":" + tmp.name );
silo::writeMultiMesh( fid, data.name, meshNames, meshTypes );
for (const auto& variable : data.variables ) {
std::vector<int> varTypes(meshnames.size(),DB_QUADVAR); // Not correct for all types
std::vector<std::string> varnames;
for ( const auto& tmp : data.domains )
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 );
@@ -325,6 +383,11 @@ void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStr
IO::initialize( );
PROFILE_START("writeData");
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
std::string path = global_IO_path + "/" + subdir;
if ( rank == 0 ) {
@@ -364,12 +427,12 @@ void IO::writeData( const std::string& subdir, const std::vector<IO::MeshDataStr
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",path.c_str());
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",path.c_str());
fprintf(fid,"%s/summary.silo\n",subdir.c_str());
fclose(fid);
} else {
ERROR("Unknown format");

328
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>
void writeUniformMesh( DBfile* fid, const std::string& meshname,
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N )
template<class TYPE> int getSiloType();
template<> int getSiloType<int>() { return DB_INT; }
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[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 );
int dims = data.size();
int err = DBWrite( fid, varname.c_str(), (void*) data.data(), &dims, 1, getSiloType<TYPE>() );
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);
if ( nvars == 1 ) {
@@ -118,24 +104,93 @@ std::vector<std::string> getVarSuffix( int ndim, int nvars )
}
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]-1;
range[2*d+0] = mesh->min_extents[d];
range[2*d+1] = mesh->max_extents[d];
}
DBFreeQuadmesh( mesh );
}
/****************************************************
* Write a vector/tensor quad variable *
****************************************************/
template<int NDIM>
void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
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);
int nvars=1, dims[NDIM]={1};
const double *vars[NDIM] = { nullptr };
int vartype = 0;
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 ) {
ERROR("Not finished");
} else if ( type == SurfaceVariable ) {
ERROR("Not finished");
} else if ( type == VolumeVariable ) {
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;
for (int d=0; d<NDIM; d++)
dims[d] = data.size(d);
@@ -152,10 +207,181 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
for (int i=0; i<nvars; i++)
varnames[i] = const_cast<char*>(var_names[i].c_str());
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 );
PROFILE_STOP("writeUniformMeshVariable",2);
}
Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname )
{
auto var = DBGetQuadvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<double> data( var->nels, var->nvals );
if ( var->datatype == DB_DOUBLE ) {
for (int i=0; i<var->nvals; i++)
memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
} else {
ERROR("Unsupported format");
}
DBFreeQuadvar( var );
std::vector<size_t> dims( var->ndims+1, var->nvals );
for (int d=0; d<var->ndims; d++)
dims[d] = var->dims[d];
data.reshape( dims );
return data;
}
/****************************************************
* Read/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 );
}
Array<double> readPointMesh( DBfile* fid, const std::string& meshname )
{
auto mesh = DBGetPointmesh( fid, meshname.c_str() );
int N = mesh->nels;
int ndim = mesh->ndims;
Array<double> coords(N,ndim);
if ( mesh->datatype == DB_DOUBLE ) {
for (int d=0; d<ndim; d++)
memcpy(&coords(0,d),mesh->coords[d],N*sizeof(double));
} else {
ERROR("Unsupported format");
}
DBFreePointmesh( mesh );
return coords;
}
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 );
}
Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname )
{
auto var = DBGetPointvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<double> data( var->nels, var->nvals );
if ( var->datatype == DB_DOUBLE ) {
for (int i=0; i<var->nvals; i++)
memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
} else {
ERROR("Unsupported format");
}
DBFreeMeshvar( var );
return data;
}
/****************************************************
* Read/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 readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coords, Array<int>& tri )
{
auto mesh = DBGetUcdmesh( fid, meshname.c_str() );
int ndim = mesh->ndims;
int N_nodes = mesh->nnodes;
coords.resize(N_nodes,ndim);
if ( mesh->datatype == DB_DOUBLE ) {
for (int d=0; d<ndim; d++)
memcpy(&coords(0,d),mesh->coords[d],N_nodes*sizeof(double));
} else {
ERROR("Unsupported format");
}
auto zones = mesh->zones;
int N_zones = zones->nzones;
int ndim_zones = zones->ndims;
ASSERT( zones->nshapes==1 );
int type = zones->shapetype[0];
int shapesize = zones->shapesize[0];
tri.resize(N_zones,shapesize);
for (int i=0; i<N_zones; i++) {
for (int j=0; j<shapesize; j++)
tri(i,j) = zones->nodelist[i*shapesize+j];
}
DBFreeUcdmesh( mesh );
}
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 );
}
Array<double> readTriMeshVariable( DBfile* fid, const std::string& varname )
{
auto var = DBGetUcdvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<double> data( var->nels, var->nvals );
if ( var->datatype == DB_DOUBLE ) {
for (int i=0; i<var->nvals; i++)
memcpy( &data(0,i), var->vals[i], var->nels*sizeof(double) );
} else {
ERROR("Unsupported format");
}
DBFreeUcdvar( var );
return data;
}
/****************************************************
@@ -181,20 +407,12 @@ void writeMultiMesh( DBfile* fid, const std::string& meshname,
****************************************************/
void writeMultiVar( DBfile* fid, const std::string& varname,
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());
auto suffix = getVarSuffix( ndim, nvar );
for (int i=0; i<nvar; i++) {
auto varname2 = varname + suffix[i];
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 );
}
std::vector<char*> varnames(varNames.size(),nullptr);
for (size_t j=0; j<varNames.size(); j++)
varnames[j] = const_cast<char*>(varNames[j].c_str());
DBPutMultivar( fid, varname.c_str(), varNames.size(), varnames.data(), (int*) varTypes.data(), nullptr );
}
@@ -202,6 +420,12 @@ void writeMultiVar( DBfile* fid, const std::string& varname,
// 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<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>& );

162
IO/silo.h
View File

@@ -44,10 +44,32 @@ DBfile* open( const std::string& filename, FileMode mode );
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
* @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] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax }
* @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
* @detailed This function writes a multimesh to silo
* @return This function returns a handle to the file
* @brief Read a uniform grid
* @detailed This function reads a uniform grid from 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
* @param[out] range Range of mesh { xmin, xmax, ymin, ymax, zmin, zmax }
* @param[out] N Number of cells in each direction
*/
void writeMultiMesh( DBfile* fid, const std::string& meshname,
const std::vector<std::string>& subMeshNames,
const std::vector<int>& subMeshTypes );
void readUniformMesh( DBfile* fid, const std::string& meshname,
std::vector<double>& range, std::vector<int>& N );
/*!
* @brief Write a uniform grid
* @detailed This function writes a uniform grid to silo as a Quadmesh
* @brief Write a uniform grid variable
* @detailed This function writes a uniform grid variable to silo as a Quadmesh
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] N Number of cells in each direction
@@ -85,6 +106,123 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
const std::string& varname, const Array<double>& data, VariableType type );
/*!
* @brief Read a uniform mesh grid variable
* @detailed This function read a uniform mesh variable to silo
* @param[in] fid Handle to the open file
* @param[in] varname Variable name
* @return Variable data
*/
Array<double> readUniformMeshVariable( DBfile* fid, const std::string& varname );
/*!
* @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 Read a pointmesh
* @detailed This function reads a pointmesh from silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @return Returns the coordinates as a N x ndim array
*/
Array<double> readPointMesh( DBfile* fid, const std::string& meshname );
/*!
* @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 Read a pointmesh grid variable
* @detailed This function reads a pointmesh variable from silo
* @param[in] fid Handle to the open file
* @param[in] varname Variable name
* @return Variable data
*/
Array<double> readPointMeshVariable( DBfile* fid, const std::string& varname );
/*!
* @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 Read a triangle mesh
* @detailed This function reads a triangle (or simplex) based mesh to silo
* @param[in] fid Handle to the open file
* @param[in] meshname Mesh name
* @param[in] coords Coordinates of the points
* @param[in] tri Coordinates of the points
*/
void readTriMesh( DBfile* fid, const std::string& meshname, Array<double>& coords, Array<int>& tri );
/*!
* @brief Write a triangle mesh grid variable
* @detailed This function writes a triangle mesh 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 Read a triangle mesh grid variable
* @detailed This function read a triangle mesh variable to silo
* @param[in] fid Handle to the open file
* @param[in] varname Variable name
* @return Variable data
*/
Array<double> readTriMeshVariable( DBfile* fid, const std::string& varname );
/*!
* @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
* @detailed This function writes a multivariable to silo
@@ -98,7 +236,7 @@ void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const s
*/
void writeMultiVar( DBfile* fid, const std::string& varname,
const std::vector<std::string>& subVarNames,
const std::vector<int>& subVarTypes, int ndim, int nvar );
const std::vector<int>& subVarTypes );
}; // silo namespace
#endif

171
tests/TestWriter.cpp
View File

@@ -74,19 +74,19 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
ut.failure("Incorrent number of timesteps");
// Check the mesh lists
for (size_t i=0; i<timesteps.size(); i++) {
for ( const auto& timestep : timesteps ) {
// Load the list of meshes and check its size
PROFILE_START(format+"-read-getMeshList");
std::vector<IO::MeshDatabase> list = IO::getMeshList(".",timesteps[i]);
auto databaseList = IO::getMeshList(".",timestep);
PROFILE_STOP(format+"-read-getMeshList");
if ( list.size()==meshData.size() )
if ( databaseList.size()==meshData.size() )
ut.passes("Corrent number of meshes found");
else
ut.failure("Incorrent number of meshes found");
// Check the number of domains for each mesh
bool pass = true;
for (size_t j=0; j<list.size(); j++)
pass = pass && (int)list[j].domains.size()==nprocs;
for ( const auto& database : databaseList )
pass = pass && (int)database.domains.size()==nprocs;
if ( pass ) {
ut.passes("Corrent number of domains for mesh");
} else {
@@ -94,17 +94,18 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
continue;
}
// For each domain, load the mesh and check its data
for (size_t j=0; j<list.size(); j++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
for ( const auto& database : databaseList ) {
pass = true;
for (size_t k=0; k<database.domains.size(); k++) {
PROFILE_START(format+"-read-getMesh");
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timestep,database,k);
PROFILE_STOP(format+"-read-getMesh");
if ( mesh.get()==NULL ) {
printf("Failed to load %s\n",list[j].name.c_str());
printf("Failed to load %s\n",database.name.c_str());
pass = false;
break;
}
if ( list[j].name=="pointmesh" ) {
if ( database.name=="pointmesh" ) {
// Check the pointmesh
std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
if ( pmesh.get()==NULL ) {
@@ -116,7 +117,7 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
break;
}
}
if ( list[j].name=="trimesh" || list[j].name=="trilist" ) {
if ( database.name=="trimesh" || database.name=="trilist" ) {
// Check the trimesh/trilist
std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
@@ -147,7 +148,7 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
pass = false;
}
}
if ( list[j].name=="domain" && format!="old" ) {
if ( database.name=="domain" && format!="old" ) {
// Check the domain mesh
const IO::DomainMesh& mesh1 = *std::dynamic_pointer_cast<IO::DomainMesh>(mesh);
if ( mesh1.nprocx!=domain->nprocx || mesh1.nprocy!=domain->nprocy || mesh1.nprocz!=domain->nprocz )
@@ -158,29 +159,30 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
pass = false;
}
}
}
if ( pass ) {
ut.passes("Meshes loaded correctly");
} else {
ut.failure("Meshes did not load correctly");
continue;
}
// For each domain, load the variables and check their data
if ( format=="old" )
continue; // Old format does not support variables
for (size_t j=0; j<list.size(); j++) {
if ( pass ) {
ut.passes("Mesh \"" + database.name + "\" loaded correctly");
} else {
ut.failure("Mesh \"" + database.name + "\" did not load correctly");
continue;
}
// Load the variables and check their data
if ( format=="old" )
continue; // Old format does not support variables
const IO::MeshDataStruct* mesh0 = NULL;
for (size_t k=0; k<meshData.size(); k++) {
if ( meshData[k].meshName == list[j].name ) {
if ( meshData[k].meshName == database.name ) {
mesh0 = &meshData[k];
break;
}
}
for (size_t v=0; v<list[j].variables.size(); v++) {
for (size_t k=0; k<list[j].domains.size(); k++) {
for (size_t k=0; k<database.domains.size(); k++) {
auto mesh = IO::getMesh(".",timestep,database,k);
for (size_t v=0; v<mesh0->vars.size(); v++) {
PROFILE_START(format+"-read-getVariable");
std::shared_ptr<const IO::Variable> variable =
IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
std::shared_ptr<IO::Variable> variable =
IO::getVariable(".",timestep,database,k,mesh0->vars[v]->name);
if ( format=="new" )
IO::reformatVariable( *mesh, *variable );
PROFILE_STOP(format+"-read-getVariable");
const IO::Variable& var1 = *mesh0->vars[v];
const IO::Variable& var2 = *variable;
@@ -193,9 +195,9 @@ void testWriter( const std::string& format, const std::vector<IO::MeshDataStruct
pass = pass && approx_equal(var1.data(m),var2.data(m));
}
if ( pass ) {
ut.passes("Variables match");
ut.passes("Variable \"" + variable->name + "\" matched");
} else {
ut.failure("Variables did not match");
ut.failure("Variable \"" + variable->name + "\" did not match");
break;
}
}
@@ -257,32 +259,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) );
// Create the variables
std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
std::shared_ptr<IO::Variable> dist_list( new IO::Variable() );
std::shared_ptr<IO::Variable> domain_var( new IO::Variable() );
dist_set1->dim = 1;
dist_list->dim = 1;
domain_var->dim = 1;
dist_set1->name = "Distance";
dist_list->name = "Distance";
domain_var->name = "Distance";
dist_set1->type = IO::NodeVariable;
dist_list->type = IO::NodeVariable;
domain_var->type = IO::VolumeVariable;
dist_set1->data.resize( N_points );
for (int i=0; i<N_points; i++)
dist_set1->data(i) = distance(set1->points[i]);
dist_list->data.resize( 3, N_tri );
for (int i=0; i<N_tri; i++) {
dist_list->data(0,i) = distance(trilist->A[i]);
dist_list->data(1,i) = distance(trilist->B[i]);
dist_list->data(2,i) = distance(trilist->C[i]);
std::shared_ptr<IO::Variable> set_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_set_mag") );
std::shared_ptr<IO::Variable> set_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_set_vec") );
std::shared_ptr<IO::Variable> list_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_list_mag") );
std::shared_ptr<IO::Variable> list_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_list_vec") );
std::shared_ptr<IO::Variable> point_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_point_mag") );
std::shared_ptr<IO::Variable> point_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_point_vec") );
std::shared_ptr<IO::Variable> domain_node_mag( new IO::Variable(1,IO::NodeVariable,"Node_domain_mag") );
std::shared_ptr<IO::Variable> domain_node_vec( new IO::Variable(3,IO::NodeVariable,"Node_domain_vec") );
std::shared_ptr<IO::Variable> set_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_set_mag") );
std::shared_ptr<IO::Variable> set_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_set_vec") );
std::shared_ptr<IO::Variable> list_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_list_mag") );
std::shared_ptr<IO::Variable> list_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_list_vec") );
std::shared_ptr<IO::Variable> domain_cell_mag( new IO::Variable(1,IO::VolumeVariable,"Cell_domain_mag") );
std::shared_ptr<IO::Variable> domain_cell_vec( new IO::Variable(3,IO::VolumeVariable,"Cell_domain_vec") );
point_node_mag->data.resize( N_points );
point_node_vec->data.resize( N_points, 3 );
for (int i=0; i<N_points; i++) {
point_node_mag->data(i) = distance(set1->points[i]);
point_node_vec->data(i,0) = set1->points[i].x;
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 j=0; j<domain->ny; j++) {
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 +338,31 @@ int main(int argc, char **argv)
std::vector<IO::MeshDataStruct> meshData(4);
meshData[0].meshName = "pointmesh";
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].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].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].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
testWriter( "old", meshData, ut );
testWriter( "new", meshData, ut );
//testWriter( "silo", meshData, ut );
testWriter( "silo", meshData, ut );
// Finished
ut.report();

4
tests/lbpm_morphdrain_pp.cpp
View File

@@ -285,7 +285,7 @@ int main(int argc, char **argv)
double sw_diff_old = 1.0;
double sw_diff_new = 1.0;
while (sw_new > SW){
while (sw_new > SW && Rcrit_new > 0.99){
Rcrit_old = Rcrit_new;
Rcrit_new -= deltaR;// decrease critical radius
@@ -298,7 +298,7 @@ int main(int argc, char **argv)
while (GlobalNumber != 0){
if (rank==0) printf("GlobalNumber=%i \n",GlobalNumber);
int LocalNumber=0;
int LocalNumber=GlobalNumber=0;
for(k=0; k<Nz; k++){
for(j=0; j<Ny; j++){
for(i=0; i<Nx; i++){