LBPM/IO/Xdmf.cpp
2021-06-02 15:36:44 -04:00

621 lines
23 KiB
C++

#include "IO/Xdmf.h"
#include "common/Array.h"
#include "common/UtilityMacros.h"
ArraySize squeeze( const ArraySize &x )
{
int Nd = 0;
size_t N[5] = { 1 };
for ( size_t i = 0; i < x.maxDim(); i++ ) {
if ( x[i] != 1 )
N[Nd++] = x[i];
}
return ArraySize( std::max( Nd, 1 ), N );
}
// Helper functions
static void addDataItem(
FILE *xmf, const std::string &indent, ArraySize size, const std::string &location )
{
size = squeeze( size );
if ( size.ndim() == 1 ) {
fprintf( xmf, "%s<DataItem Dimensions=\"%lu\"", indent.data(), size[0] );
} else if ( size.ndim() == 2 ) {
fprintf( xmf, "%s<DataItem Dimensions=\"%lu %lu\"", indent.data(), size[1], size[0] );
} else if ( size.ndim() == 3 ) {
fprintf( xmf, "%s<DataItem Dimensions=\"%lu %lu %lu\"", indent.data(), size[2], size[1],
size[0] );
} else if ( size.ndim() == 4 ) {
fprintf( xmf, "%s<DataItem Dimensions=\"%lu %lu %lu %lu\"", indent.data(), size[3], size[2],
size[1], size[0] );
} else {
ERROR( "Invalid number of dimensions" );
}
fprintf( xmf, " Format=\"HDF\">\n" );
fprintf( xmf, "%s %s\n", indent.data(), location.data() );
fprintf( xmf, "%s</DataItem>\n", indent.data() );
}
template<class TYPE>
static void addVariable( FILE *xmf, const std::string &indent, const std::string &name,
const std::string &type, const std::string &center, ArraySize size,
const std::string &location )
{
fprintf( xmf, "%s<Attribute Name=\"%s\" AttributeType=\"%s\" Center=\"%s\">\n", indent.data(),
name.data(), type.data(), center.data() );
addDataItem( xmf, indent + " ", size, location );
fprintf( xmf, "%s</Attribute>\n", indent.data() );
}
/****************************************************************
* Enum functions *
****************************************************************/
/*template<class TYPE>
static Xdmf::DataType getDataType()
{
if ( std::is_same<TYPE, char>::value )
return Xdmf::DataType::Char;
else if ( std::is_same<TYPE, int32_t>::value )
return Xdmf::DataType::Int32;
else if ( std::is_same<TYPE, int64_t>::value )
return Xdmf::DataType::Int64;
else if ( std::is_same<TYPE, uint32_t>::value )
return Xdmf::DataType::Uint32;
else if ( std::is_same<TYPE, uint64_t>::value )
return Xdmf::DataType::Uint64;
else if ( std::is_same<TYPE, float>::value )
return Xdmf::DataType::Float;
else if ( std::is_same<TYPE, double>::value )
return Xdmf::DataType::Double;
else
ERROR( "Invalid type" );
}*/
static const char *TopologyTypeNames[] = { "", "Polyvertex", "Polyline", "Polygon", "Triangle",
"Quadrilateral", "Tetrahedron", "Pyramid", "Wedge", "Hexahedron", "Edge_3", "Triangle_6",
"Quadrilateral_8", "Tetrahedron_10", "Pyramid_13", "Wedge_15", "Hexahedron_20", "Mixed",
"CurvilinearMesh2D", "CurvilinearMesh3D", "RectangularMesh2D", "RectangularMesh3D",
"UniformMesh2D", "UniformMesh3D" };
static const uint8_t TopologyTypeDOFs[] = { 0, 1, 2, 0, 3, 4, 4, 5, 6, 8, 3, 6, 8, 10, 13, 15, 20,
0, 0, 0, 0, 0, 0, 0 };
/****************************************************************
* Create a mesh *
****************************************************************/
Xdmf::MeshData Xdmf::createPointMesh( const std::string &name, uint8_t NDIM, size_t N,
const std::string &x, const std::string &y, const std::string &z )
{
return createUnstructuredMesh( name, NDIM, TopologyType::Polyvertex, N, "", N, x, y, z );
}
Xdmf::MeshData Xdmf::createUniformMesh(
const std::string &name, const std::vector<double> &range, ArraySize size )
{
ASSERT( range.size() == 2 * size.ndim() );
MeshData data;
data.name = name;
data.size = size;
if ( size.ndim() == 2 )
data.type = TopologyType::UniformMesh2D;
else if ( size.ndim() == 3 )
data.type = TopologyType::UniformMesh3D;
else
ERROR( "# of dimensions != 2 or 3" );
for ( int i = 0; i < 2 * size.ndim(); i++ )
data.range[i] = range[i];
return data;
}
Xdmf::MeshData Xdmf::createCurvilinearMesh( const std::string &name, ArraySize size,
const std::string &x, const std::string &y, const std::string &z )
{
MeshData data;
data.name = name;
if ( size.ndim() == 2 )
data.type = TopologyType::CurvilinearMesh2D;
else if ( size.ndim() == 3 )
data.type = TopologyType::CurvilinearMesh3D;
else
ERROR( "Invalid size for Curvilinear mesh" );
data.size = size;
data.x = x;
data.y = y;
data.z = z;
return data;
}
Xdmf::MeshData Xdmf::createUnstructuredMesh( const std::string &name, uint8_t NDIM,
TopologyType type, size_t NumElements, const std::string &dofMap, size_t NumNodes,
const std::string &x, const std::string &y, const std::string &z )
{
ASSERT( type != TopologyType::Null );
MeshData data;
data.name = name;
data.type = type;
data.size = { NDIM, NumElements, NumNodes };
data.dofMap = dofMap;
data.x = x;
data.y = y;
data.z = z;
return data;
}
/****************************************************************
* Add a variable *
****************************************************************/
void Xdmf::MeshData::addVariable( const std::string &meshName, const std::string &varName,
ArraySize varSize, RankType rank, Center center, const std::string &varData )
{
VarData var;
var.name = varName;
var.size = varSize;
var.data = varData;
var.rankType = rank;
var.center = center;
vars.push_back( std::move( var ) );
}
/****************************************************************
* Add a mesh domain *
****************************************************************/
void Xdmf::addMesh( const std::string &meshName, const MeshData &domain )
{
auto &domains = d_meshData[meshName];
for ( const auto &domain2 : domains )
ASSERT( domain2.name != domain.name );
domains.push_back( domain );
}
/****************************************************************
* Write a variable *
****************************************************************/
static void writeVariable( FILE *fid, const Xdmf::VarData &var, const std::string &indent )
{
// Write the variable name
fprintf( fid, "%s<Attribute Name=\"%s\"", indent.data(), var.name.data() );
// Write the variable type
if ( var.rankType == Xdmf::RankType::Scalar ) {
fprintf( fid, " AttributeType=\"Scalar\"" );
} else if ( var.rankType == Xdmf::RankType::Vector ) {
fprintf( fid, " AttributeType=\"Vector\"" );
} else if ( var.rankType == Xdmf::RankType::Tensor ) {
fprintf( fid, " AttributeType=\"Tensor\"" );
} else if ( var.rankType == Xdmf::RankType::Tensor6 ) {
fprintf( fid, " AttributeType=\"Tensor6\"" );
} else if ( var.rankType == Xdmf::RankType::Matrix ) {
fprintf( fid, " AttributeType=\"Matrix\"" );
} else if ( var.rankType == Xdmf::RankType::GlobalID ) {
fprintf( fid, " AttributeType=\"GlobalID\"" );
} else {
ERROR( "Unknown center type" );
}
// Write the variable centering
if ( var.center == Xdmf::Center::Node ) {
fprintf( fid, " Center=\"Node\">\n" );
} else if ( var.center == Xdmf::Center::Cell ) {
fprintf( fid, " Center=\"Cell\">\n" );
} else if ( var.center == Xdmf::Center::Grid ) {
fprintf( fid, " Center=\"Grid\">\n" );
} else if ( var.center == Xdmf::Center::Face ) {
fprintf( fid, " Center=\"Face\">\n" );
} else if ( var.center == Xdmf::Center::Edge ) {
fprintf( fid, " Center=\"Edge\">\n" );
} else if ( var.center == Xdmf::Center::Other ) {
fprintf( fid, " Center=\"Other\">\n" );
} else {
ERROR( "Unknown center type" );
}
// Write the data item
addDataItem( fid, indent + " ", var.size, var.data );
// Finished
fprintf( fid, "%s</Attribute>\n", indent.data() );
}
/****************************************************************
* Write the mesh grid *
****************************************************************/
static void writeMeshGrid( FILE *fid, const Xdmf::MeshData &mesh, const std::string &indent )
{
const char *s = indent.data();
double x0[3] = { mesh.range[0], mesh.range[2], mesh.range[4] };
double dx[3] = { ( mesh.range[1] - mesh.range[0] ) / mesh.size[0],
( mesh.range[3] - mesh.range[2] ) / mesh.size[1],
( mesh.range[5] - mesh.range[4] ) / mesh.size[2] };
switch ( mesh.type ) {
case Xdmf::TopologyType::UniformMesh2D:
// Write a uniform 2d mesh
fprintf( fid, "%s<Grid Name=\"%s\" GridType=\"Uniform\">\n", s, mesh.name.data() );
fprintf( fid,
"%s <Topology TopologyType=\"2DCoRectMesh\" NumberOfElements=\"%lu %lu\"/>\n", s,
mesh.size[1] + 1, mesh.size[0] + 1 );
fprintf( fid, "%s <Geometry GeometryType=\"ORIGIN_DXDY\">\n", s );
fprintf(
fid, "%s <DataItem Format=\"XML\" NumberType=\"float\" Dimensions=\"2\">\n", s );
fprintf( fid, "%s %0.12e %0.12e\n", s, x0[0], x0[1] );
fprintf( fid, "%s </DataItem>\n", s );
fprintf(
fid, "%s <DataItem Format=\"XML\" NumberType=\"float\" Dimensions=\"2\">\n", s );
fprintf( fid, "%s %0.12e %0.12e\n", s, dx[0], dx[1] );
fprintf( fid, "%s </DataItem>\n", s );
fprintf( fid, "%s </Geometry>\n", s );
break;
case Xdmf::TopologyType::UniformMesh3D:
// Write a uniform 3d mesh
fprintf( fid, "%s<Grid Name=\"%s\" GridType=\"Uniform\">\n", s, mesh.name.data() );
fprintf( fid,
"%s <Topology TopologyType=\"3DCoRectMesh\" NumberOfElements=\"%lu %lu\"/>\n", s,
mesh.size[1] + 1, mesh.size[0] + 1 );
fprintf( fid, "%s <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n", s );
fprintf(
fid, "%s <DataItem Format=\"XML\" NumberType=\"float\" Dimensions=\"3\">\n", s );
fprintf( fid, "%s %0.12e %0.12e %0.12e\n", s, x0[0], x0[1], x0[2] );
fprintf( fid, "%s </DataItem>\n", s );
fprintf(
fid, "%s <DataItem Format=\"XML\" NumberType=\"float\" Dimensions=\"3\">\n", s );
fprintf( fid, "%s %0.12e %0.12e %0.12e\n", s, dx[0], dx[1], dx[2] );
fprintf( fid, "%s </DataItem>\n", s );
fprintf( fid, "%s </Geometry>\n", s );
break;
case Xdmf::TopologyType::CurvilinearMesh2D:
// Write a 2D curvillinear mesh
fprintf( fid, "%s<Grid Name=\"%s\" GridType=\"Uniform\">\n", s, mesh.name.data() );
fprintf( fid, "%s <Topology TopologyType=\"2DSMesh\" NumberOfElements=\"%lu %lu\"/>\n", s,
mesh.size[1] + 1, mesh.size[0] + 1 );
fprintf( fid, "%s <Geometry GeometryType=\"X_Y\">\n", s );
addDataItem( fid, indent + " ", mesh.size + 1, mesh.x );
addDataItem( fid, indent + " ", mesh.size + 1, mesh.y );
fprintf( fid, "%s </Geometry>\n", s );
break;
case Xdmf::TopologyType::CurvilinearMesh3D:
// Write a 3D curvillinear mesh
fprintf( fid, "%s<Grid Name=\"%s\" GridType=\"Uniform\">\n", s, mesh.name.data() );
fprintf( fid, "%s <Topology TopologyType=\"3DSMesh\" NumberOfElements=\"%lu %lu %lu\"/>\n",
s, mesh.size[2] + 1, mesh.size[1] + 1, mesh.size[0] + 1 );
fprintf( fid, "%s <Geometry GeometryType=\"X_Y_Z\">\n", s );
addDataItem( fid, indent + " ", mesh.size + 1, mesh.x );
addDataItem( fid, indent + " ", mesh.size + 1, mesh.y );
addDataItem( fid, indent + " ", mesh.size + 1, mesh.z );
fprintf( fid, "%s </Geometry>\n", s );
break;
case Xdmf::TopologyType::Polyvertex:
case Xdmf::TopologyType::Polyline:
case Xdmf::TopologyType::Polygon:
case Xdmf::TopologyType::Triangle:
case Xdmf::TopologyType::Quadrilateral:
case Xdmf::TopologyType::Tetrahedron:
case Xdmf::TopologyType::Pyramid:
case Xdmf::TopologyType::Wedge:
case Xdmf::TopologyType::Hexahedron:
case Xdmf::TopologyType::Edge_3:
case Xdmf::TopologyType::Triangle_6:
case Xdmf::TopologyType::Quadrilateral_8:
case Xdmf::TopologyType::Tetrahedron_10:
case Xdmf::TopologyType::Pyramid_13:
case Xdmf::TopologyType::Wedge_15:
case Xdmf::TopologyType::Hexahedron_20:
// Write an unstructured mesh
{
int NDIM = mesh.size[0];
size_t Nelem = mesh.size[1];
size_t Nnode = mesh.size[2];
uint8_t Ndofs = TopologyTypeDOFs[static_cast<int>( mesh.type )];
auto type = TopologyTypeNames[static_cast<int>( mesh.type )];
fprintf( fid, "%s<Grid Name=\"%s\">\n", s, mesh.name.data() );
fprintf( fid, "%s <Topology TopologyType=\"%s\"", s, type );
fprintf( fid, " NumberOfElements=\"%lu\">\n", Nelem );
if ( !mesh.dofMap.empty() )
addDataItem( fid, indent + " ", { Ndofs, Nelem }, mesh.dofMap );
fprintf( fid, "%s </Topology>\n", s );
if ( NDIM == 2 ) {
if ( mesh.y.empty() ) {
fprintf( fid, "%s <Geometry GeometryType=\"XY\">\n", s );
addDataItem( fid, indent + " ", { 2, Nnode }, mesh.x );
} else {
fprintf( fid, "%s <Geometry GeometryType=\"X_Y\">\n", s );
addDataItem( fid, indent + " ", Nnode, mesh.x );
addDataItem( fid, indent + " ", Nnode, mesh.y );
}
} else if ( NDIM == 3 ) {
if ( mesh.y.empty() ) {
fprintf( fid, "%s <Geometry GeometryType=\"XYZ\">\n", s );
addDataItem( fid, indent + " ", { 2, Nnode }, mesh.x );
} else {
fprintf( fid, "%s <Geometry GeometryType=\"X_Y_Z\">\n", s );
addDataItem( fid, indent + " ", Nnode, mesh.x );
addDataItem( fid, indent + " ", Nnode, mesh.y );
addDataItem( fid, indent + " ", Nnode, mesh.z );
}
} else {
ERROR( "Dimensions other than 2 or 3 are not supported" );
}
fprintf( fid, "%s </Geometry>\n", s );
}
break;
default: {
auto msg = "Invalid mesh type: " + std::to_string( static_cast<int>( mesh.type ) ) + " - " +
mesh.name;
ERROR( msg );
}
}
// Write the variables
for ( const auto &var : mesh.vars )
writeVariable( fid, var, indent + " " );
fprintf( fid, "%s </Grid>\n", s );
}
/****************************************************************
* Write the XDMF xml file *
****************************************************************/
void Xdmf::write( const std::string &filename ) const
{
// Create XDMF file
auto fid = fopen( filename.data(), "w" );
fprintf( fid, "<?xml version=\"1.0\" ?>\n" );
fprintf( fid, "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n" );
fprintf( fid, "<Xdmf Version=\"2.0\">\n" );
fprintf( fid, "<Domain>\n" );
// Write an empty mesh to enable collections to work properly
fprintf( fid, " <Grid Name=\"\" GridType=\"Uniform\"></Grid>\n\n" );
// Write each mesh
for ( const auto &data : d_meshData ) {
auto name = data.first;
auto domains = data.second;
if ( domains.empty() )
continue;
if ( domains.size() == 1u && name == domains[0].name ) {
writeMeshGrid( fid, domains[0], " " );
} else {
fprintf( fid, " <Grid Name=\"%s\" GridType=\"Collection\">\n", name.data() );
for ( const auto &domain : domains )
writeMeshGrid( fid, domain, " " );
fprintf( fid, " </Grid>\n\n" );
}
}
fprintf( fid, "</Domain>\n" );
fprintf( fid, "</Xdmf>\n" );
fclose( fid );
}
/****************************************************************
* Pack/Unpack data *
****************************************************************/
template<class T>
typename std::enable_if<std::is_trivially_copyable<T>::value, size_t>::type size( const T & )
{
return sizeof( T );
}
template<class T>
typename std::enable_if<std::is_trivially_copyable<T>::value, char *>::type pack(
char *ptr, const T &x )
{
memcpy( ptr, &x, sizeof( T ) );
return ptr + sizeof( T );
}
template<class T>
typename std::enable_if<std::is_trivially_copyable<T>::value, char *>::type unpack(
char *ptr, T &x )
{
memcpy( &x, ptr, sizeof( T ) );
return ptr + sizeof( T );
}
static size_t size( const std::string &str ) { return sizeof( int ) + str.size(); }
static char *pack( char *ptr, const std::string &str )
{
int N = str.size();
memcpy( ptr, &N, sizeof( int ) );
ptr += sizeof( int );
memcpy( ptr, str.data(), str.size() );
ptr += str.size();
return ptr;
}
static char *unpack( char *ptr, std::string &str )
{
int N = 0;
memcpy( &N, ptr, sizeof( int ) );
ASSERT( N >= 0 && N < 1000 );
ptr += sizeof( int );
str = std::string( ptr, N );
ptr += N;
return ptr;
}
static size_t size( const Xdmf::VarData &data )
{
size_t bytes = 0;
bytes += size( data.name );
bytes += size( data.size );
bytes += size( data.rankType );
bytes += size( data.center );
bytes += size( data.data );
return bytes;
}
static char *pack( char *ptr, const Xdmf::VarData &data )
{
ptr = pack( ptr, data.name );
ptr = pack( ptr, data.size );
ptr = pack( ptr, data.rankType );
ptr = pack( ptr, data.center );
ptr = pack( ptr, data.data );
return ptr;
}
static char *unpack( char *ptr, Xdmf::VarData &data )
{
int rankType = 0, center = 0;
ptr = unpack( ptr, data.name );
ptr = unpack( ptr, data.size );
ptr = unpack( ptr, rankType );
ptr = unpack( ptr, center );
ptr = unpack( ptr, data.data );
data.rankType = static_cast<Xdmf::RankType>( rankType );
data.center = static_cast<Xdmf::Center>( center );
return ptr;
}
static size_t size( const Xdmf::MeshData &data )
{
int N_vars = data.vars.size();
size_t bytes = 0;
bytes += size( data.name );
bytes += size( data.type );
bytes += size( data.size );
bytes += size( data.range );
bytes += size( data.x );
bytes += size( data.y );
bytes += size( data.z );
bytes += size( N_vars );
for ( int i = 0; i < N_vars; i++ )
bytes += size( data.vars[i] );
return bytes;
}
static char *pack( char *ptr, const Xdmf::MeshData &data )
{
int N_vars = data.vars.size();
ptr = pack( ptr, data.name );
ptr = pack( ptr, data.type );
ptr = pack( ptr, data.size );
ptr = pack( ptr, data.range );
ptr = pack( ptr, data.x );
ptr = pack( ptr, data.y );
ptr = pack( ptr, data.z );
ptr = pack( ptr, N_vars );
for ( int i = 0; i < N_vars; i++ )
ptr = pack( ptr, data.vars[i] );
return ptr;
}
static char *unpack( char *ptr, Xdmf::MeshData &data )
{
int N_vars = 0;
ptr = unpack( ptr, data.name );
ptr = unpack( ptr, data.type );
ptr = unpack( ptr, data.size );
ptr = unpack( ptr, data.range );
ptr = unpack( ptr, data.x );
ptr = unpack( ptr, data.y );
ptr = unpack( ptr, data.z );
ptr = unpack( ptr, N_vars );
data.vars.resize( N_vars );
for ( int i = 0; i < N_vars; i++ )
ptr = unpack( ptr, data.vars[i] );
return ptr;
}
static size_t size( const std::vector<Xdmf::MeshData> &data )
{
size_t bytes = 0;
int N = data.size();
bytes += size( N );
for ( int i = 0; i < N; i++ )
bytes += size( data[i] );
return bytes;
}
static char *pack( char *ptr, const std::vector<Xdmf::MeshData> &data )
{
int N = data.size();
ptr = pack( ptr, N );
for ( int i = 0; i < N; i++ )
ptr = pack( ptr, data[i] );
return ptr;
}
static char *unpack( char *ptr, std::vector<Xdmf::MeshData> &data )
{
data.clear();
int N = data.size();
ptr = unpack( ptr, N );
data.resize( N );
for ( int i = 0; i < N; i++ )
ptr = unpack( ptr, data[i] );
return ptr;
}
static size_t size( const std::map<std::string, std::vector<Xdmf::MeshData>> &data )
{
size_t bytes = 0;
int N_map = data.size();
bytes += size( N_map );
for ( const auto &tmp : data ) {
bytes += size( tmp.first );
bytes += size( tmp.second );
}
return bytes;
}
static char *pack( char *ptr, const std::map<std::string, std::vector<Xdmf::MeshData>> &data )
{
int N_map = data.size();
ptr = pack( ptr, N_map );
for ( const auto &tmp : data ) {
ptr = pack( ptr, tmp.first );
ptr = pack( ptr, tmp.second );
}
return ptr;
}
static char *unpack( char *ptr, std::map<std::string, std::vector<Xdmf::MeshData>> &data )
{
data.clear();
int N_map = data.size();
ptr = unpack( ptr, N_map );
for ( int i = 0; i < N_map; i++ ) {
std::string name;
std::vector<Xdmf::MeshData> data2;
ptr = unpack( ptr, name );
ptr = unpack( ptr, data2 );
data[name] = std::move( data2 );
}
return ptr;
}
/****************************************************************
* Gather all data to rank 0 *
****************************************************************/
void Xdmf::gather( const Utilities::MPI &comm )
{
if ( comm.getRank() == 0 ) {
for ( int i = 1; i < comm.getSize(); i++ ) {
// Recieve the data
size_t N_meshes = 0, N_bytes = 0;
comm.recv( &N_meshes, 1, i, 717 );
comm.recv( &N_bytes, 1, i, 718 );
auto buf = new char[N_bytes];
comm.recv( buf, N_bytes, i, 719 );
// Unpack the data
std::map<std::string, std::vector<MeshData>> data;
unpack( buf, data );
delete[] buf;
// Add the meshes
for ( auto tmp : data ) {
const auto &name = tmp.first;
const auto &domains = tmp.second;
if ( domains.size() == 1u && domains[0].name == name ) {
// We are dealing with a single mesh
ASSERT( d_meshData.find( name ) == d_meshData.end() );
d_meshData.insert( tmp );
} else {
// Add the domains
auto &meshes = d_meshData[name];
for ( auto domain : domains ) {
for ( const auto &tmp : meshes )
ASSERT( tmp.name != domain.name );
meshes.push_back( domain );
}
}
}
}
} else {
// Send the number of meshes
size_t N_meshes = d_meshData.size();
comm.send( &N_meshes, 1, 0, 717 );
// Pack the send data
size_t N_bytes = size( d_meshData );
comm.send( &N_bytes, 1, 0, 718 );
auto buf = new char[N_bytes];
pack( buf, d_meshData );
// Send the data to rank 0
comm.send( buf, N_bytes, 0, 719 );
delete[] buf;
// Clear the internal data
d_meshData.clear();
}
}