Skip to content
Permalink
Browse files

update MDAL 0.3.1

  • Loading branch information
PeterPetrik committed Apr 18, 2019
1 parent 1fac5a6 commit da7f655393922aaef7b444f951b64401e2ef4a7e
@@ -28,4 +28,3 @@ FIND_LIBRARY (NETCDF_LIBRARY
INCLUDE (FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS (NetCDF
DEFAULT_MSG NETCDF_LIBRARY NETCDF_INCLUDE_DIR)

@@ -9,6 +9,8 @@

#cmakedefine HAVE_NETCDF

#cmakedefine HAVE_XML

#endif // MDAL_CONFIG_HPP


@@ -12,6 +12,8 @@
#include <vector>
#include <map>
#include <cassert>
#include <limits>
#include <algorithm>

#include "mdal_2dm.hpp"
#include "mdal.h"
@@ -63,6 +65,18 @@ size_t MDAL::Mesh2dm::vertexIndex( size_t vertexID ) const
return vertexID;
}

size_t MDAL::Mesh2dm::maximumVertexId() const
{
size_t maxIndex = verticesCount() - 1;
if ( mVertexIDtoIndex.empty() )
return maxIndex;
else
{
// std::map is sorted!
size_t maxID = mVertexIDtoIndex.rbegin()->first;
return std::max( maxIndex, maxID );
}
}

MDAL::Driver2dm::Driver2dm():
Driver( DRIVER_NAME,
@@ -135,6 +149,9 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
std::vector<Vertex> vertices( vertexCount );
std::vector<Face> faces( faceCount );

// Basement 3.x supports definition of elevation for cell centers
std::vector<double> elementCenteredElevation;

in.clear();
in.seekg( 0, std::ios::beg );

@@ -143,33 +160,44 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
size_t faceIndex = 0;
size_t vertexIndex = 0;
std::map<size_t, size_t> vertexIDtoIndex;
size_t lastVertexID = 0;

while ( std::getline( in, line ) )
{
if ( startsWith( line, "E4Q" ) )
if ( startsWith( line, "E4Q" ) ||
startsWith( line, "E3T" )
)
{
chunks = split( line, ' ' );
assert( faceIndex < faceCount );

const size_t faceVertexCount = MDAL::toSizeT( line[1] );
assert( ( faceVertexCount == 3 ) || ( faceVertexCount == 4 ) );

Face &face = faces[faceIndex];
face.resize( 4 );
face.resize( faceVertexCount );

// chunks format here
// E** id vertex_id1, vertex_id2, ... material_id (elevation - optional)
// vertex ids are numbered from 1
// Right now we just store node IDs here - we will convert them to node indices afterwards
for ( size_t i = 0; i < 4; ++i )
face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1
assert( chunks.size() > faceVertexCount + 1 );

faceIndex++;
}
else if ( startsWith( line, "E3T" ) )
{
chunks = split( line, ' ' );
assert( faceIndex < faceCount );
for ( size_t i = 0; i < faceVertexCount; ++i )
face[i] = MDAL::toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1

Face &face = faces[faceIndex];
face.resize( 3 );
// Right now we just store node IDs here - we will convert them to node indices afterwards
for ( size_t i = 0; i < 3; ++i )
// OK, now find out if there is optional cell elevation (BASEMENT 3.x)
if ( chunks.size() == faceVertexCount + 4 )
{
face[i] = toSizeT( chunks[i + 2] ) - 1; // 2dm is numbered from 1

// initialize dataset if it is still empty
if ( elementCenteredElevation.empty() )
{
elementCenteredElevation = std::vector<double>( faceCount, std::numeric_limits<double>::quiet_NaN() );
}

// add Bed Elevation (Face) value
elementCenteredElevation[faceIndex] = MDAL::toDouble( chunks[ faceVertexCount + 3 ] );
}

faceIndex++;
@@ -192,7 +220,22 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
else if ( startsWith( line, "ND" ) )
{
chunks = split( line, ' ' );
size_t nodeID = toSizeT( chunks[1] ) - 1; // 2dm is numbered from 1
size_t nodeID = toSizeT( chunks[1] );

if ( nodeID != 0 )
{
// specification of 2DM states that ID should be positive integer numbered from 1
// but it seems some formats do not respect that
if ( ( lastVertexID != 0 ) && ( nodeID <= lastVertexID ) )
{
// the algorithm requires that the file has NDs orderer by index
if ( status ) *status = MDAL_Status::Err_InvalidData;
return nullptr;
}
lastVertexID = nodeID;
}
nodeID -= 1; // 2dm is numbered from 1

_parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID, status );
assert( vertexIndex < vertexCount );
Vertex &vertex = vertices[vertexIndex];
@@ -236,6 +279,10 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
);
mesh->faces = faces;
mesh->vertices = vertices;

// Add Bed Elevations
MDAL::addFaceScalarDatasetGroup( mesh.get(), elementCenteredElevation, "Bed Elevation (Face)" );
MDAL::addBedElevationDatasetGroup( mesh.get(), vertices );

return std::unique_ptr<Mesh>( mesh.release() );
}
@@ -29,19 +29,50 @@ namespace MDAL
~Mesh2dm() override;


//! Some formats supports gaps in the vertex indexing, but we return continuos array from MDAL
//! for most of the formats this returns
//! HYDRO_AS-2D supports gaps in the vertex indexing,
//! but we use continuos array of vertices in MDAL
//! \param vertexID internal index/ID of the vertex that native format uses
//! \returns index of the vertex in the continuous array of vertices we returned by readVertices()
//! \returns index of the vertex in the continuous array of vertices we returned by readVertices().
//! For invalid vertexID it is returned index that is out of vertices array bounds.
virtual size_t vertexIndex( size_t vertexID ) const;

//! Returns maximum vertex ID.
//! For meshes without gaps in vertex indexing, it is vertex count - 1
virtual size_t maximumVertexId() const;

private:
// 2dm supports "gaps" in the mesh indexing
// Store only the indices that have different index and ID
// https://github.com/lutraconsulting/MDAL/issues/51
//! 2dm supports "gaps" in the mesh indexing
//! Store only the indices that have different index and ID
//! https://github.com/lutraconsulting/MDAL/issues/51
std::map<size_t, size_t> mVertexIDtoIndex;
};

/**
* 2DM format specification used in TUFLOW, HYDRO_AS-2D and BASEMENET solvers
* Text file format representing mesh vertices (ND) and faces (E**)
* ND id x y z
* Supports lines, triangles and polygons up to 9 vertices (implemented triangles and quads)
* E3T id 1 2 3 mat_id -> face type, id, vertex indices ..., material index
*
* full specification here: https://www.xmswiki.com/wiki/SMS:2D_Mesh_Files_*.2dm
*
* Exception for the official specification is for recognition of cell-centered
* elevation values supported by BASEMENT 3.x releases
* If face definition has extra column, it is parsed and recognized as
* elevation, e.g. format for triangle
* E3T id 1 2 3 mat_id elevation
* and added automatically as "Bed Elevation (Face)"
*
* Note that some 2dm formats do have some extra columns after mat_id column with
* data with unknown origin/name (e.g. tests/data/2dm/regular_grid.2dm)
*
* HYDRO_AS-2D also allows gaps in vertex indexing. In this case we support only files
* where the vertices are sorted by ID in the source file (limitation of the implementation)
*
* Vertex/Face IDs should be indexed from 1. We support indexing from 0 for datasets in xmdf format,
* but not for ascii dat format (since the loop is from 0 to maximumVertexId() which is in this case
* numberical_limits<size_t>::max(); (limitation of the implementation)
*/
class Driver2dm: public Driver
{
public:
@@ -13,6 +13,7 @@
#include <map>
#include <cassert>
#include <memory>
#include <limits>

#include "mdal_ascii_dat.hpp"
#include "mdal.h"
@@ -103,7 +104,8 @@ void MDAL::DriverAsciiDat::loadOldFormat( std::ifstream &in,
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->verticesCount() != fileNodeCount )
size_t meshIdCount = maximumId( mesh ) + 1;
if ( meshIdCount != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( cardType == "SCALAR" || cardType == "VECTOR" )
@@ -167,7 +169,8 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
if ( cardType == "ND" && items.size() >= 2 )
{
size_t fileNodeCount = toSizeT( items[1] );
if ( mesh->verticesCount() != fileNodeCount )
size_t meshIdCount = maximumId( mesh ) + 1;
if ( meshIdCount != fileNodeCount )
EXIT_WITH_ERROR( MDAL_Status::Err_IncompatibleMesh );
}
else if ( cardType == "NC" && items.size() >= 2 )
@@ -247,6 +250,15 @@ void MDAL::DriverAsciiDat::loadNewFormat( std::ifstream &in,
}
}

size_t MDAL::DriverAsciiDat::maximumId( const MDAL::Mesh *mesh ) const
{
const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
if ( m2dm )
return m2dm->maximumVertexId();
else
return mesh->verticesCount() - 1;
}

/**
* The DAT format contains "datasets" and each dataset has N-outputs. One output
* represents data for all vertices/faces for one timestep
@@ -265,6 +277,14 @@ void MDAL::DriverAsciiDat::load( const std::string &datFile, MDAL::Mesh *mesh, M
return;
}

size_t mID = maximumId( mesh );
if ( mID == std::numeric_limits<size_t>::max() )
{
// This happens when mesh is 2DM and vertices are numbered from 0
if ( status ) *status = MDAL_Status::Err_IncompatibleMesh;
return;
}

std::ifstream in( mDatFile, std::ifstream::in );
std::string line;
if ( !std::getline( in, line ) )
@@ -316,17 +336,22 @@ void MDAL::DriverAsciiDat::readVertexTimestep(

const Mesh2dm *m2dm = dynamic_cast<const Mesh2dm *>( mesh );
double *values = dataset->values();
for ( size_t i = 0; i < mesh->verticesCount(); ++i )
size_t meshIdCount = maximumId( mesh ) + 1; // these are native format indexes (IDs). For formats without gaps it equals vertex array index
size_t vertexCount = mesh->verticesCount();

for ( size_t id = 0; id < meshIdCount; ++id )
{
std::string line;
std::getline( stream, line );
std::vector<std::string> tsItems = split( line, ' ' );

size_t index;
if ( m2dm )
index = m2dm->vertexIndex( i );
index = m2dm->vertexIndex( id ); //this index may be out of values array
else
index = i;
index = id;

if ( index >= vertexCount ) continue;

if ( isVector )
{
@@ -35,6 +35,9 @@ namespace MDAL
* such dataset, the dataset name contains "_els_" substring
* (e.g. depth_els_1.dat)
*
* HYDRO_AS-2D solver can have mesh that has numbering gaps, but
* speficies values for even missing indexes in dataset file
*
* In one file, there is always one dataset group stored.
*
* Sometime the "older" datasets may have some part of the
@@ -60,6 +63,12 @@ namespace MDAL
void loadOldFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;
void loadNewFormat( std::ifstream &in, Mesh *mesh, MDAL_Status *status ) const;

//! Gets maximum (native) index.
//! For meshes without indexing gap it is vertexCount - 1
//! For some HYDRO_AS-2D meshes with indexing gaps, it returns
//! maximum native index of the vertex in defined in the mesh
size_t maximumId( const Mesh *mesh ) const;

void readVertexTimestep(
const Mesh *mesh,
std::shared_ptr<DatasetGroup> group,
@@ -125,10 +125,7 @@ MDAL::cfdataset_info_map MDAL::DriverCF::parseDatasetGroupInfo()
}
while ( true );

if ( dsinfo_map.size() == 0 )
{
throw MDAL_Status::Err_InvalidData;
}
if ( dsinfo_map.size() == 0 ) throw MDAL_Status::Err_InvalidData;

return dsinfo_map;
}
@@ -388,16 +385,6 @@ void MDAL::CFDimensions::setDimension( MDAL::CFDimensions::Type type,
mCount[type] = count;
}

size_t MDAL::CFDimensions::faceCount() const
{
return size( Face2D ) + size( Line1D );
}

size_t MDAL::CFDimensions::vertexCount() const
{
return size( Vertex1D ) + size( Vertex2D );
}

bool MDAL::CFDimensions::isDatasetType( MDAL::CFDimensions::Type type ) const
{
return (
@@ -41,10 +41,6 @@ namespace MDAL
size_t size( Type type ) const;
//! Sets a dimension
void setDimension( Type type, size_t count, int ncid = -1 );
//! Returns number of faces (1D + 2D)
size_t faceCount() const;
//! Returns number of vertices (1D + 2D)
size_t vertexCount() const;
//! Returns whether the type is one that case be used for datasets definition
bool isDatasetType( Type type ) const;

@@ -502,23 +502,6 @@ void MDAL::DriverFlo2D::createMesh( const std::vector<CellCenter> &cells, double
mMesh->vertices = vertices;
}

bool MDAL::DriverFlo2D::isFlo2DFile( const std::string &fileName )
{
std::vector<std::string> required_files =
{
"CADPTS.DAT",
"FPLAIN.DAT"
};

for ( const std::string &str : required_files )
{
std::string fn( fileNameFromDir( fileName, str ) );
if ( !fileExists( fn ) )
return false;
}
return true;
}

bool MDAL::DriverFlo2D::parseHDF5Datasets( const std::string &datFileName )
{
//return true on error
@@ -25,8 +25,6 @@ namespace MDAL
bool canRead( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &resultsFile, MDAL_Status *status ) override;

static bool isFlo2DFile( const std::string &fileName );

private:
struct CellCenter
{

0 comments on commit da7f655

Please sign in to comment.
You can’t perform that action at this time.