Skip to content

Commit

Permalink
[FEATURE] Mesh 1D Renderer (#34848)
Browse files Browse the repository at this point in the history 
* MDAL 0.5.90 : support for custom Logger and 1D meshes
* [FEATURE] [MESH] Support rendering of 1D meshes, see qgis/QGIS-Enhancement-Proposals#164

1D mesh consist of edges (edge is straight line segment with 2 vertices) and the data that is defined on either
vertices or edges. Such data can be loaded by MDAL and rendered as mesh layer in QGIS.
  • Loading branch information
@PeterPetrik
PeterPetrik committed Mar 9, 2020
1 parent a4b9774 commit 0675c0b
Show file tree
Hide file tree
Showing 113 changed files with 3,558 additions and 1,032 deletions.
262 changes: 202 additions & 60 deletions external/mdal/api/mdal.h
View file

Large diffs are not rendered by default.

85 changes: 52 additions & 33 deletions external/mdal/frmts/mdal_2dm.cpp
View file
Expand Up @@ -18,17 +18,20 @@
#include "mdal_2dm.hpp"
#include "mdal.h"
#include "mdal_utils.hpp"
#include "mdal_logger.hpp"

#define DRIVER_NAME "2DM"

MDAL::Mesh2dm::Mesh2dm( size_t verticesCount,
size_t edgesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
MDAL::BBox extent,
const std::string &uri,
const std::map<size_t, size_t> vertexIDtoIndex )
: MemoryMesh( DRIVER_NAME,
verticesCount,
edgesCount,
facesCount,
faceVerticesMaximumCount,
extent,
Expand All @@ -39,15 +42,15 @@ MDAL::Mesh2dm::Mesh2dm( size_t verticesCount,

MDAL::Mesh2dm::~Mesh2dm() = default;

bool _parse_vertex_id_gaps( std::map<size_t, size_t> &vertexIDtoIndex, size_t vertexIndex, size_t vertexID, MDAL_Status *status )
bool _parse_vertex_id_gaps( std::map<size_t, size_t> &vertexIDtoIndex, size_t vertexIndex, size_t vertexID )
{
if ( vertexIndex == vertexID )
return false;

std::map<size_t, size_t>::iterator search = vertexIDtoIndex.find( vertexID );
if ( search != vertexIDtoIndex.end() )
{
if ( status ) *status = MDAL_Status::Warn_ElementNotUnique;
MDAL::Log::warning( Warn_ElementNotUnique, DRIVER_NAME, "could not find vertex" );
return true;
}

Expand Down Expand Up @@ -105,22 +108,23 @@ bool MDAL::Driver2dm::canReadMesh( const std::string &uri )
return true;
}

std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile, MDAL_Status *status )
std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile )
{
mMeshFile = meshFile;

if ( status ) *status = MDAL_Status::None;
MDAL::Log::resetLastStatus();

std::ifstream in( mMeshFile, std::ifstream::in );
std::string line;
if ( !std::getline( in, line ) || !startsWith( line, "MESH2D" ) )
{
if ( status ) *status = MDAL_Status::Err_UnknownFormat;
MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), meshFile + " could not be opened" );
return nullptr;
}

size_t faceCount = 0;
size_t vertexCount = 0;
size_t edgesCount = 0;

// Find out how many nodes and elements are contained in the .2dm mesh file
while ( std::getline( in, line ) )
Expand All @@ -134,20 +138,24 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
{
vertexCount++;
}
else if ( startsWith( line, "E2L" ) ||
startsWith( line, "E3L" ) ||
else if ( startsWith( line, "E2L" ) )
{
edgesCount++;
}
else if ( startsWith( line, "E3L" ) ||
startsWith( line, "E6T" ) ||
startsWith( line, "E8Q" ) ||
startsWith( line, "E9Q" ) )
{
if ( status ) *status = MDAL_Status::Warn_UnsupportedElement;
faceCount += 1; // We still count them as elements
MDAL::Log::warning( MDAL_Status::Err_UnsupportedElement, name(), "found unsupported element" );
return nullptr;
}
}

// Allocate memory
std::vector<Vertex> vertices( vertexCount );
std::vector<Face> faces( faceCount );
Vertices vertices( vertexCount );
Edges edges( edgesCount );
Faces faces( faceCount );

// Basement 3.x supports definition of elevation for cell centers
std::vector<double> elementCenteredElevation;
Expand All @@ -159,6 +167,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,

size_t faceIndex = 0;
size_t vertexIndex = 0;
size_t edgeIndex = 0;
std::map<size_t, size_t> vertexIDtoIndex;
size_t lastVertexID = 0;

Expand Down Expand Up @@ -202,20 +211,18 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,

faceIndex++;
}
else if ( startsWith( line, "E2L" ) ||
startsWith( line, "E3L" ) ||
startsWith( line, "E6T" ) ||
startsWith( line, "E8Q" ) ||
startsWith( line, "E9Q" ) )
else if ( startsWith( line, "E2L" ) )
{
// We do not yet support these elements
// format: E2L id n1 n2 matid
chunks = split( line, ' ' );
assert( faceIndex < faceCount );

//size_t elemID = toSizeT( chunks[1] );
assert( false ); //TODO mark element as unusable

faceIndex++;
assert( edgeIndex < edgesCount );
assert( chunks.size() > 4 );
size_t startVertexIndex = MDAL::toSizeT( chunks[2] ) - 1; // 2dm is numbered from 1
size_t endVertexIndex = MDAL::toSizeT( chunks[3] ) - 1; // 2dm is numbered from 1
Edge &edge = edges[edgeIndex];
edge.startVertex = startVertexIndex;
edge.endVertex = endVertexIndex;
edgeIndex++;
}
else if ( startsWith( line, "ND" ) )
{
Expand All @@ -229,14 +236,14 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
if ( ( lastVertexID != 0 ) && ( nodeID <= lastVertexID ) )
{
// the algorithm requires that the file has NDs orderer by index
if ( status ) *status = MDAL_Status::Err_InvalidData;
MDAL::Log::error( MDAL_Status::Err_InvalidData, name(), "nodes are not ordered by index" );
return nullptr;
}
lastVertexID = nodeID;
}
nodeID -= 1; // 2dm is numbered from 1

_parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID, status );
_parse_vertex_id_gaps( vertexIDtoIndex, vertexIndex, nodeID );
assert( vertexIndex < vertexCount );
Vertex &vertex = vertices[vertexIndex];
vertex.x = toDouble( chunks[2] );
Expand All @@ -260,7 +267,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
}
else if ( vertices.size() < nodeID )
{
if ( status ) *status = MDAL_Status::Warn_ElementWithInvalidNode;
MDAL::Log::warning( MDAL_Status::Warn_ElementWithInvalidNode, name(), "found invalid node" );
}
}
//TODO check validity of the face
Expand All @@ -270,6 +277,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
std::unique_ptr< Mesh2dm > mesh(
new Mesh2dm(
vertices.size(),
edges.size(),
faces.size(),
MAX_VERTICES_PER_FACE_2DM,
computeExtent( vertices ),
Expand All @@ -279,6 +287,7 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
);
mesh->faces = faces;
mesh->vertices = vertices;
mesh->edges = edges;

// Add Bed Elevations
MDAL::addFaceScalarDatasetGroup( mesh.get(), elementCenteredElevation, "Bed Elevation (Face)" );
Expand All @@ -287,21 +296,21 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
return std::unique_ptr<Mesh>( mesh.release() );
}

void MDAL::Driver2dm::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Status *status )
void MDAL::Driver2dm::save( const std::string &uri, MDAL::Mesh *mesh )
{
if ( status ) *status = MDAL_Status::None;
MDAL::Log::resetLastStatus();

std::ofstream file( uri, std::ofstream::out );

if ( !file.is_open() )
{
if ( status ) *status = MDAL_Status::Err_FailToWriteToDisk;
MDAL::Log::error( MDAL_Status::Err_FailToWriteToDisk, name(), "Could not open file " + uri );
}

std::string line = "MESH2D";
file << line << std::endl;

//write vertices
// write vertices
std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();
double vertex[3];
for ( size_t i = 0; i < mesh->verticesCount(); ++i )
Expand All @@ -320,12 +329,12 @@ void MDAL::Driver2dm::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Statu
file << line << std::endl;
}

//write faces
// write faces
std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
for ( size_t i = 0; i < mesh->facesCount(); ++i )
{
int faceOffsets[1];
int vertexIndices[4]; //max 4 vertices for a face
int vertexIndices[MAX_VERTICES_PER_FACE_2DM];
faceIterator->next( 1, faceOffsets, 4, vertexIndices );

if ( faceOffsets[0] > 2 && faceOffsets[0] < 5 )
Expand All @@ -343,9 +352,19 @@ void MDAL::Driver2dm::save( const std::string &uri, MDAL::Mesh *mesh, MDAL_Statu
line.append( std::to_string( vertexIndices[j] + 1 ) );
}
}

file << line << std::endl;
}

// write edges
std::unique_ptr<MDAL::MeshEdgeIterator> edgeIterator = mesh->readEdges();
for ( size_t i = 0; i < mesh->edgesCount(); ++i )
{
int startIndex;
int endIndex;
edgeIterator->next( 1, &startIndex, &endIndex );
line = "E2L " + std::to_string( mesh->facesCount() + i + 1 ) + " " + std::to_string( startIndex + 1 ) + " " + std::to_string( endIndex + 1 ) + " 1";
file << line << std::endl;
}

file.close();
}
10 changes: 7 additions & 3 deletions external/mdal/frmts/mdal_2dm.hpp
View file
Expand Up @@ -23,6 +23,7 @@ namespace MDAL
{
public:
Mesh2dm( size_t verticesCount,
size_t edgesCount,
size_t facesCount,
size_t faceVerticesMaximumCount,
BBox extent,
Expand Down Expand Up @@ -54,7 +55,10 @@ namespace MDAL
* 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)
* The format supports lines, triangles and quads, where the elememts could be
* stored with some intermediate points (e.g. triangle defined by 6 vertices, vertices
* and the middle of the edges) We do support only simple definition, so E6T, E8Q and E9Q
* are not supported.
* 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
Expand Down Expand Up @@ -87,8 +91,8 @@ namespace MDAL
{return MAX_VERTICES_PER_FACE_2DM;}

bool canReadMesh( const std::string &uri ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile, MDAL_Status *status ) override;
void save( const std::string &uri, Mesh *mesh, MDAL_Status *status ) override;
std::unique_ptr< Mesh > load( const std::string &meshFile ) override;
void save( const std::string &uri, Mesh *mesh ) override;

private:
std::string mMeshFile;
Expand Down
6 changes: 3 additions & 3 deletions external/mdal/frmts/mdal_3di.cpp
View file
Expand Up @@ -57,13 +57,13 @@ void MDAL::Driver3Di::populateFacesAndVertices( Vertices &vertices, Faces &faces
int ncidX = mNcFile->getVarId( "Mesh2DContour_x" );
double fillX = mNcFile->getFillValue( ncidX );
std::vector<double> faceVerticesX( arrsize );
if ( nc_get_var_double( mNcFile->handle(), ncidX, faceVerticesX.data() ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_get_var_double( mNcFile->handle(), ncidX, faceVerticesX.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );

// Y coordinate
int ncidY = mNcFile->getVarId( "Mesh2DContour_y" );
double fillY = mNcFile->getFillValue( ncidY );
std::vector<double> faceVerticesY( arrsize );
if ( nc_get_var_double( mNcFile->handle(), ncidY, faceVerticesY.data() ) ) throw MDAL_Status::Err_UnknownFormat;
if ( nc_get_var_double( mNcFile->handle(), ncidY, faceVerticesY.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );

// now populate create faces and backtrack which vertices
// are used in multiple faces
Expand Down Expand Up @@ -135,7 +135,7 @@ void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
"Bed Elevation"
);

group->setDataLocation( MDAL_DataLocation::DataOnFaces2D );
group->setDataLocation( MDAL_DataLocation::DataOnFaces );
group->setIsScalar( true );

std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
Expand Down

0 comments on commit 0675c0b

Please sign in to comment.