Skip to content
Permalink
Browse files

update to MDAL 0.7.1, fixes bug in FLO2D reader

  • Loading branch information
PeterPetrik authored and nyalldawson committed Oct 1, 2020
1 parent b23288b commit 9c9b070d4753aa6e1b598a8a75149c8c15e360b9
@@ -22,21 +22,7 @@

#define FLO2D_NAN 0.0

struct VertexCompare
{
bool operator()( const MDAL::Vertex &lhs, const MDAL::Vertex &rhs ) const
{
double resX = 0;
resX += lhs.x * 1000000;
resX += lhs.y * 1000;

double resY = 0;
resY += rhs.x * 1000000;
resY += rhs.y * 1000;

return resX < resY;
}
};
#define INVALID_INDEX std::numeric_limits<size_t>::max()

static std::string fileNameFromDir( const std::string &mainFileName, const std::string &name )
{
@@ -99,7 +85,8 @@ void MDAL::DriverFlo2D::addStaticDataset(
mMesh->datasetGroups.push_back( group );
}

void MDAL::DriverFlo2D::parseCADPTSFile( const std::string &datFileName, std::vector<CellCenter> &cells )

void MDAL::DriverFlo2D::parseCADPTSFile( const std::string &datFileName, std::vector<CellCenter> &cells, MDAL::BBox &cellCenterExtent )
{
std::string cadptsFile( fileNameFromDir( datFileName, "CADPTS.DAT" ) );
if ( !MDAL::fileExists( cadptsFile ) )
@@ -122,8 +109,16 @@ void MDAL::DriverFlo2D::parseCADPTSFile( const std::string &datFileName, std::ve
cc.id = MDAL::toSizeT( lineParts[0] ) - 1; //numbered from 1
cc.x = MDAL::toDouble( lineParts[1] );
cc.y = MDAL::toDouble( lineParts[2] );
cc.conn.resize( 4 );
cells.push_back( cc );

if ( cc.x > cellCenterExtent.maxX )
cellCenterExtent.maxX = cc.x;
if ( cc.x < cellCenterExtent.minX )
cellCenterExtent.minX = cc.x;
if ( cc.y > cellCenterExtent.maxY )
cellCenterExtent.maxY = cc.y;
if ( cc.y < cellCenterExtent.minY )
cellCenterExtent.minY = cc.y;
}
}

@@ -481,7 +476,8 @@ void MDAL::DriverFlo2D::createMesh1d( const std::string &datFileName, const std:

void MDAL::DriverFlo2D::parseFPLAINFile( std::vector<double> &elevations,
const std::string &datFileName,
std::vector<CellCenter> &cells )
std::vector<CellCenter> &cells,
double &cellSize )
{
elevations.clear();
// FPLAIN.DAT - CONNECTIVITY (ELEM NUM, ELEM N, ELEM E, ELEM S, ELEM W, MANNING-N, BED ELEVATION)
@@ -494,21 +490,41 @@ void MDAL::DriverFlo2D::parseFPLAINFile( std::vector<double> &elevations,
std::ifstream fplainStream( fplainFile, std::ifstream::in );
std::string line;

bool cellSizeCalculated = false;

while ( std::getline( fplainStream, line ) )
{
line = MDAL::rtrim( line );
std::vector<std::string> lineParts = MDAL::split( line, ' ' );
if ( lineParts.size() != 7 )
{
throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while loading Fplain file, wrong lineparts count (7)" );
throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while loading FPLAIN.DAT file, wrong lineparts count (7)" );
}
size_t cc_i = MDAL::toSizeT( lineParts[0] ) - 1; //numbered from 1
for ( size_t j = 0; j < 4; ++j )

if ( !cellSizeCalculated )
{
cells[cc_i].conn[j] = MDAL::toInt( lineParts[j + 1] ) - 1; //numbered from 1, 0 boundary Vertex
size_t cc_i = MDAL::toSizeT( lineParts[0] ) - 1; //numbered from 1
for ( int i = 1; i < 5; ++i ) //search the first cell that have a neighbor to calculate cell size
{
int neighborCell = MDAL::toInt( lineParts[i] );
if ( neighborCell != 0 )
{
if ( i % 2 == 1 ) //North or South neighbor
cellSize = fabs( cells[neighborCell - 1].y - cells[cc_i].y );
else // East or West
cellSize = fabs( cells[neighborCell - 1].x - cells[cc_i].x );

cellSizeCalculated = true;
break;
}
}
}

elevations.push_back( MDAL::toDouble( lineParts[6] ) );
}

if ( !cellSizeCalculated )
throw MDAL::Error( MDAL_Status::Err_IncompatibleMesh, "Only isolated cell(s), not possible to calculate cell size" );
}

static void addDatasetToGroup( std::shared_ptr<MDAL::DatasetGroup> group, std::shared_ptr<MDAL::MemoryDataset2D> dataset )
@@ -756,31 +772,6 @@ void MDAL::DriverFlo2D::parseVELFPVELOCFile( const std::string &datFileName )
addStaticDataset( maxVel, "Velocity/Maximums", datFileName );
}

double MDAL::DriverFlo2D::calcCellSize( const std::vector<CellCenter> &cells )
{
// find first cell that is not izolated from the others
// and return its distance to the neighbor's cell center
for ( size_t i = 0; i < cells.size(); ++i )
{
for ( size_t j = 0; j < 4; ++j )
{
int idx = cells[i].conn[0];
if ( idx > -1 )
{
if ( ( j == 0 ) || ( j == 2 ) )
{
return fabs( cells[static_cast<size_t>( idx )].y - cells[i].y );
}
else
{
return fabs( cells[static_cast<size_t>( idx )].x - cells[i].x );
}
}
}
}
throw MDAL::Error( MDAL_Status::Err_IncompatibleMesh, "Did not find izolated cell" );
}

MDAL::Vertex MDAL::DriverFlo2D::createVertex( size_t position, double half_cell_size, const CellCenter &cell )
{
MDAL::Vertex n;
@@ -813,38 +804,69 @@ MDAL::Vertex MDAL::DriverFlo2D::createVertex( size_t position, double half_cell_
return n;
}

void MDAL::DriverFlo2D::createMesh2d( const std::vector<CellCenter> &cells, double half_cell_size )
void MDAL::DriverFlo2D::createMesh2d( const std::vector<CellCenter> &cells, const BBox &cellCenterExtent, double cell_size )
{
// Create all Faces from cell centers.
// Vertexs must be also created, they are not stored in FLO-2D files
// try to reuse Vertexs already created for other Faces by usage of unique_Vertexs set.
Faces faces;

double half_cell_size = cell_size / 2;
Faces faces( cells.size(), Face( 4 ) );

BBox vertexExtent( cellCenterExtent.minX - half_cell_size,
cellCenterExtent.maxX + half_cell_size,
cellCenterExtent.minY - half_cell_size,
cellCenterExtent.maxY + half_cell_size );

size_t width = ( vertexExtent.maxX - vertexExtent.minX ) / cell_size + 1;
size_t heigh = ( vertexExtent.maxY - vertexExtent.minY ) / cell_size + 1;
std::vector<std::vector<size_t>> vertexGrid( width, std::vector<size_t>( heigh, INVALID_INDEX ) );

Vertices vertices;
std::map<Vertex, size_t, VertexCompare> unique_vertices; //vertex -> id
size_t vertex_idx = 0;

for ( size_t i = 0; i < cells.size(); ++i )
{
Face e( 4 );
Face &e = faces[i];

size_t xVertexIdx = ( cells[i].x - vertexExtent.minX ) / cell_size;
size_t yVertexIdx = ( cells[i].y - vertexExtent.minY ) / cell_size;

for ( size_t position = 0; position < 4; ++position )
{
Vertex n = createVertex( position, half_cell_size, cells[i] );
const auto iter = unique_vertices.find( n );
if ( iter == unique_vertices.end() )
size_t xPos;
size_t yPos;

switch ( position )
{
unique_vertices[n] = vertex_idx;
vertices.push_back( n );
e[position] = vertex_idx;
++vertex_idx;
case 0:
xPos = 1;
yPos = 0;
break;

case 1:
xPos = 1;
yPos = 1;
break;

case 2:
xPos = 0;
yPos = 1;
break;

case 3:
xPos = 0;
yPos = 0;
break;
}
else

if ( vertexGrid[xVertexIdx + xPos][yVertexIdx + yPos] == INVALID_INDEX )
{
e[position] = iter->second;
vertices.push_back( createVertex( position, half_cell_size, cells.at( i ) ) );
vertexGrid[xVertexIdx + xPos][yVertexIdx + yPos] = vertices.size() - 1;
}
}

faces.push_back( e );
e[position] = vertexGrid[xVertexIdx + xPos][yVertexIdx + yPos];
}
}

mMesh.reset(
@@ -1074,11 +1096,11 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverFlo2D::loadMesh1d()
{
std::vector<CellCenter> cells;
std::map<size_t, size_t> cellsIdToVertex;

MDAL::BBox cellCenterExtent;
try
{
// Parse cells position
parseCADPTSFile( mDatFileName, cells );
parseCADPTSFile( mDatFileName, cells, cellCenterExtent );
createMesh1d( mDatFileName, cells, cellsIdToVertex );
parseHYCHANFile( mDatFileName, cellsIdToVertex );
}
@@ -1100,13 +1122,14 @@ std::unique_ptr<MDAL::Mesh> MDAL::DriverFlo2D::loadMesh2d()
try
{
// Parse mMesh info
parseCADPTSFile( mDatFileName, cells );
MDAL::BBox cellCenterExtent;
parseCADPTSFile( mDatFileName, cells, cellCenterExtent );
std::vector<double> elevations;
parseFPLAINFile( elevations, mDatFileName, cells );
double cell_size = calcCellSize( cells );
double cell_size;
parseFPLAINFile( elevations, mDatFileName, cells, cell_size );

// Create mMesh
createMesh2d( cells, cell_size / 2.0 );
createMesh2d( cells, cellCenterExtent, cell_size );

// create output for bed elevation
addStaticDataset( elevations, "Bed Elevation", mDatFileName );
@@ -66,28 +66,26 @@ namespace MDAL
size_t id;
double x;
double y;
std::vector<int> conn; // north, east, south, west cell center index, -1 boundary Vertex
};

std::unique_ptr< MDAL::MemoryMesh > mMesh;
std::string mDatFileName;

void createMesh2d( const std::vector<CellCenter> &cells, double half_cell_size );
void createMesh2d( const std::vector<CellCenter> &cells, const MDAL::BBox &cellCenterExtent, double cell_size );
void createMesh1d( const std::string &datFileName, const std::vector<CellCenter> &cells, std::map<size_t, size_t> &cellsIdToVertex );

void parseOUTDatasets( const std::string &datFileName, const std::vector<double> &elevations );
bool parseHDF5Datasets( MDAL::MemoryMesh *mesh, const std::string &timedepFileName );
void parseVELFPVELOCFile( const std::string &datFileName );
void parseDEPTHFile( const std::string &datFileName, const std::vector<double> &elevations );
void parseTIMDEPFile( const std::string &datFileName, const std::vector<double> &elevations );
void parseFPLAINFile( std::vector<double> &elevations, const std::string &datFileName, std::vector<CellCenter> &cells );
void parseCADPTSFile( const std::string &datFileName, std::vector<CellCenter> &cells );
void parseFPLAINFile( std::vector<double> &elevations, const std::string &datFileName, std::vector<CellCenter> &cells, double &cellZize );
void parseCADPTSFile( const std::string &datFileName, std::vector<CellCenter> &cells, BBox &cellCenterExtent );
void parseCHANBANKFile( const std::string &datFileName, std::map<size_t, size_t> &cellIdToVertices, std::map<size_t, std::vector<size_t> > &duplicatedRightBankToVertex, size_t &verticesCount );
void parseCHANFile( const std::string &datFileName, const std::map<size_t, size_t> &cellIdToVertices, std::vector<Edge> &edges );
void parseHYCHANFile( const std::string &datFileName, const std::map<size_t, size_t> &cellIdToVertices );
void addStaticDataset( std::vector<double> &vals, const std::string &groupName, const std::string &datFileName );
static MDAL::Vertex createVertex( size_t position, double half_cell_size, const CellCenter &cell );
static double calcCellSize( const std::vector<CellCenter> &cells );

std::unique_ptr< Mesh > loadMesh2d();
std::unique_ptr<Mesh> loadMesh1d();
@@ -21,7 +21,7 @@ static const char *EMPTY_STR = "";

const char *MDAL_Version()
{
return "0.7.0";
return "0.7.1";
}

MDAL_Status MDAL_LastStatus()
@@ -25,10 +25,10 @@ namespace MDAL
BBox() {}
BBox( double lx, double ux, double ly, double uy ): minX( lx ), maxX( ux ), minY( ly ), maxY( uy ) {}

double minX;
double maxX;
double minY;
double maxY;
double minX = std::numeric_limits<double>::max();
double maxX = -std::numeric_limits<double>::max();
double minY = std::numeric_limits<double>::max();
double maxY = -std::numeric_limits<double>::max();
};

typedef struct

0 comments on commit 9c9b070

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