upgrade MDAL to 0.3.0

Author: PeterPetrik

external/mdal/cmake/FindNetCDF.cmake

@@ -28,4 +28,3 @@ FIND_LIBRARY (NETCDF_LIBRARY
 INCLUDE (FindPackageHandleStandardArgs)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS (NetCDF
   DEFAULT_MSG NETCDF_LIBRARY NETCDF_INCLUDE_DIR)
-

external/mdal/cmake_templates/mdal_config.hpp.in

@@ -9,6 +9,8 @@
 
 #cmakedefine HAVE_NETCDF
 
+#cmakedefine HAVE_XML
+
 #endif // MDAL_CONFIG_HPP
 
 

external/mdal/frmts/mdal_2dm.cpp

@@ -12,6 +12,7 @@
 #include <vector>
 #include <map>
 #include <cassert>
+#include <limits>
 
 #include "mdal_2dm.hpp"
 #include "mdal.h"
@@ -63,7 +64,6 @@ size_t MDAL::Mesh2dm::vertexIndex( size_t vertexID ) const
   return vertexID;
 }
 
-
 MDAL::Driver2dm::Driver2dm():
   Driver( DRIVER_NAME,
           "2DM Mesh File",
@@ -135,6 +135,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 );
 
@@ -146,30 +149,40 @@ std::unique_ptr<MDAL::Mesh> MDAL::Driver2dm::load( const std::string &meshFile,
 
   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++;
@@ -236,6 +249,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() );
 }

external/mdal/frmts/mdal_2dm.hpp

@@ -42,6 +42,25 @@ namespace MDAL
       std::map<size_t, size_t> mVertexIDtoIndex;
   };
 
+  /**
+   * 2DM format specification used in TUFLOW 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)
+   */
   class Driver2dm: public Driver
   {
     public:

external/mdal/frmts/mdal_cf.cpp

@@ -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 (

external/mdal/frmts/mdal_cf.hpp

@@ -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;
 

external/mdal/frmts/mdal_flo2d.cpp

@@ -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

external/mdal/frmts/mdal_flo2d.hpp

@@ -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
       {

external/mdal/frmts/mdal_hec2d.cpp

@@ -18,51 +18,36 @@
 static HdfFile openHdfFile( const std::string &fileName )
 {
   HdfFile file( fileName );
-  if ( !file.isValid() )
-  {
-    throw MDAL_Status::Err_UnknownFormat;
-  }
+  if ( !file.isValid() ) throw MDAL_Status::Err_UnknownFormat;
   return file;
 }
 
 static HdfGroup openHdfGroup( const HdfFile &hdfFile, const std::string &name )
 {
   HdfGroup grp = hdfFile.group( name );
-  if ( !grp.isValid() )
-  {
-    throw MDAL_Status::Err_UnknownFormat;
-  }
+  if ( !grp.isValid() ) throw MDAL_Status::Err_UnknownFormat;
   return grp;
 }
 
 static HdfGroup openHdfGroup( const HdfGroup &hdfGroup, const std::string &name )
 {
   HdfGroup grp = hdfGroup.group( name );
-  if ( !grp.isValid() )
-  {
-    throw MDAL_Status::Err_UnknownFormat;
-  }
+  if ( !grp.isValid() ) throw MDAL_Status::Err_UnknownFormat;
   return grp;
 }
 
 static HdfDataset openHdfDataset( const HdfGroup &hdfGroup, const std::string &name )
 {
   HdfDataset dsFileType = hdfGroup.dataset( name );
-  if ( !dsFileType.isValid() )
-  {
-    throw MDAL_Status::Err_UnknownFormat;
-  }
+  if ( !dsFileType.isValid() ) throw MDAL_Status::Err_UnknownFormat;
   return dsFileType;
 }
 
 
 static std::string openHdfAttribute( const HdfFile &hdfFile, const std::string &name )
 {
   HdfAttribute attr = hdfFile.attribute( name );
-  if ( !attr.isValid() )
-  {
-    throw MDAL_Status::Err_UnknownFormat;
-  }
+  if ( !attr.isValid() ) throw MDAL_Status::Err_UnknownFormat;
   return attr.readString();
 }
 
@@ -360,10 +345,7 @@ std::vector<std::string> MDAL::DriverHec2D::read2DFlowAreasNamesOld( HdfGroup gG
 {
   HdfDataset dsNames = openHdfDataset( gGeom2DFlowAreas, "Names" );
   std::vector<std::string> names = dsNames.readArrayString();
-  if ( names.empty() )
-  {
-    throw MDAL_Status::Err_InvalidData;
-  }
+  if ( names.empty() ) throw MDAL_Status::Err_InvalidData;
   return names;
 }
 
@@ -432,10 +414,7 @@ std::vector<std::string> MDAL::DriverHec2D::read2DFlowAreasNames505( HdfGroup gG
   H5Tclose( attributeHID );
   H5Tclose( stringHID );
   std::vector<std::string> names;
-  if ( attributes.empty() )
-  {
-    throw MDAL_Status::Err_InvalidData;
-  }
+  if ( attributes.empty() ) throw MDAL_Status::Err_InvalidData;
 
   for ( const auto &attr : attributes )
   {