Skip to content
Permalink
Browse files

update MDAL 0.3.1

  • Loading branch information
PeterPetrik committed Apr 18, 2019
1 parent fd492ab commit f7420443925c38ca17c3883b8dd115f95beadd04
@@ -13,6 +13,7 @@
#include <map>
#include <cassert>
#include <limits>
#include <algorithm>

#include "mdal_2dm.hpp"
#include "mdal.h"
@@ -64,6 +65,19 @@ 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,
"2DM Mesh File",
@@ -146,6 +160,7 @@ 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 ) )
{
@@ -205,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];
@@ -29,21 +29,26 @@ 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 and BASEMENET solvers
* 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)
@@ -60,6 +65,13 @@ namespace MDAL
*
* 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
{
@@ -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,
@@ -317,7 +317,28 @@ void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
{
assert( raster_band );

double nodata = GDALGetRasterNoDataValue( raster_band, nullptr );
// nodata
int pbSuccess;
double nodata = GDALGetRasterNoDataValue( raster_band, &pbSuccess );
if ( pbSuccess == 0 ) nodata = std::numeric_limits<double>::quiet_NaN();
bool hasNoData = !std::isnan( nodata );

// offset and scale
double offset = 0.0;
double scale = GDALGetRasterScale( raster_band, &pbSuccess );
if ( ( pbSuccess == 0 ) || MDAL::equals( scale, 0.0 ) || std::isnan( scale ) )
{
scale = 1.0;
}
else
{
offset = GDALGetRasterOffset( raster_band, &pbSuccess );
if ( ( pbSuccess == 0 ) || std::isnan( offset ) )
{
offset = 0.0;
}
}

unsigned int mXSize = meshGDALDataset()->mXSize;
unsigned int mYSize = meshGDALDataset()->mYSize;

@@ -349,8 +370,11 @@ void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared
{
unsigned int idx = x + mXSize * y;
double val = mPafScanline[x];
if ( !MDAL::equals( val, nodata ) )
if ( !hasNoData || !MDAL::equals( val, nodata ) )
{
// Apply scale and offset
val = val * scale + offset;

// values is prepolulated with NODATA values, so store only legal values
if ( is_vector )
{
@@ -153,7 +153,7 @@ size_t MDAL::XdmfFunctionDataset::scalarData( size_t indexStart, size_t count, d
assert( mType != FunctionType::Join );

if ( mType == FunctionType::Subtract )
return substractFunction( indexStart, count, buffer );
return subtractFunction( indexStart, count, buffer );

if ( mType == FunctionType::Flow )
return flowFunction( indexStart, count, buffer );
@@ -176,7 +176,7 @@ size_t MDAL::XdmfFunctionDataset::activeData( size_t indexStart, size_t count, i
return count;
}

size_t MDAL::XdmfFunctionDataset::substractFunction( size_t indexStart, size_t count, double *buffer )
size_t MDAL::XdmfFunctionDataset::subtractFunction( size_t indexStart, size_t count, double *buffer )
{
std::vector<double> buf( 2 * count, std::numeric_limits<double>::quiet_NaN() );
size_t copyVals = extractRawData( indexStart, count, 2, buf );
@@ -86,7 +86,7 @@ namespace MDAL
* Currently we do not use any fancy bison/flex based
* expression parsing, just supporting few types of
* most common function types:
* - substraction (A-B)
* - subtraction (A-B)
* - join ( [A, B] vector)
* - magnitude
*
@@ -126,7 +126,7 @@ namespace MDAL
size_t activeData( size_t indexStart, size_t count, int *buffer ) override;

private:
size_t substractFunction( size_t indexStart, size_t count, double *buffer );
size_t subtractFunction( size_t indexStart, size_t count, double *buffer );
size_t flowFunction( size_t indexStart, size_t count, double *buffer );
size_t joinFunction( size_t indexStart, size_t count, double *buffer );
size_t extractRawData( size_t indexStart, size_t count, size_t nDatasets, std::vector<double> &buf );
@@ -22,7 +22,7 @@ static MDAL_Status sLastStatus;

const char *MDAL_Version()
{
return "0.3.0";
return "0.3.1";
}

MDAL_Status MDAL_LastStatus()

0 comments on commit f742044

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