Skip to content

Commit

Permalink
update MDAL 0.7.3
Browse files Browse the repository at this point in the history
  • Loading branch information
@PeterPetrik @nyalldawson
PeterPetrik authored and nyalldawson committed Nov 30, 2020
1 parent 50fd353 commit a34ff39
Show file tree
Hide file tree
Showing 12 changed files with 154 additions and 86 deletions.
4 changes: 2 additions & 2 deletions external/mdal/api/mdal.h
View file
Expand Up @@ -234,10 +234,10 @@ MDAL_EXPORT void MDAL_CloseMesh( MDAL_MeshH mesh );
/**
* Creates a empty mesh in memory
*
* \since MDAL 0.7
*
* \note the mesh is editable (vertices and faces can be added, see MDAL_M_addVertices() and MDAL_M_addFaces()),
* and can be saved with MDAL_SaveMesh()
*
* \since MDAL 0.7
*/
MDAL_EXPORT MDAL_MeshH MDAL_CreateMesh( MDAL_DriverH driver );

Expand Down
2 changes: 1 addition & 1 deletion external/mdal/frmts/mdal_ascii_dat.cpp
View file
Expand Up @@ -485,7 +485,7 @@ bool MDAL::DriverAsciiDat::persist( MDAL::DatasetGroup *group )
if ( !MDAL::contains( uri, "_els" ) && group->dataLocation() != MDAL_DataLocation::DataOnVertices )
{
// Should contain _els in name for edges/faces dataset but it does not
int pos = uri.size() - 4;
int pos = MDAL::toInt( uri.size() ) - 4;
uri.insert( std::max( 0, pos ), "_els" );
group->replaceUri( uri );
}
Expand Down
18 changes: 9 additions & 9 deletions external/mdal/frmts/mdal_flo2d.cpp
View file
Expand Up @@ -144,8 +144,8 @@ void MDAL::DriverFlo2D::parseCHANBANKFile( const std::string &datFileName,
{
throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while loading CHANBANK file, wrong lineparts count (2)" );
}
int leftBank = MDAL::toSizeT( lineParts[0] ) - 1; //numbered from 1
int rightBank = MDAL::toSizeT( lineParts[1] ) - 1;
int leftBank = MDAL::toInt( MDAL::toSizeT( lineParts[0] ) ) - 1; //numbered from 1
int rightBank = MDAL::toInt( MDAL::toSizeT( lineParts[1] ) ) - 1;

std::map<size_t, size_t>::const_iterator it = cellIdToVertices.find( rightBank );
if ( it != cellIdToVertices.end() )
Expand Down Expand Up @@ -197,7 +197,7 @@ void MDAL::DriverFlo2D::parseCHANFile( const std::string &datFileName, const std
{
throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while loading CHAN file, wrong chanel element line" );
}
int currentCellId = MDAL::toSizeT( lineParts[1] ) - 1;
int currentCellId = MDAL::toInt( MDAL::toSizeT( lineParts[1] ) ) - 1;
if ( previousCellId >= 0 )
{
std::map<size_t, size_t>::const_iterator it1 = cellIdToVertices.find( previousCellId );
Expand Down Expand Up @@ -818,8 +818,8 @@ void MDAL::DriverFlo2D::createMesh2d( const std::vector<CellCenter> &cells, cons
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;
size_t width = MDAL::toSizeT( ( vertexExtent.maxX - vertexExtent.minX ) / cell_size + 1 );
size_t heigh = MDAL::toSizeT( ( vertexExtent.maxY - vertexExtent.minY ) / cell_size + 1 );
std::vector<std::vector<size_t>> vertexGrid( width, std::vector<size_t>( heigh, INVALID_INDEX ) );

Vertices vertices;
Expand All @@ -828,13 +828,13 @@ void MDAL::DriverFlo2D::createMesh2d( const std::vector<CellCenter> &cells, cons
{
Face &e = faces[i];

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

for ( size_t position = 0; position < 4; ++position )
{
size_t xPos;
size_t yPos;
size_t xPos = 0;
size_t yPos = 0;

switch ( position )
{
Expand Down
122 changes: 85 additions & 37 deletions external/mdal/frmts/mdal_hec2d.cpp
View file
Expand Up @@ -202,21 +202,6 @@ static std::vector<MDAL::RelativeTimestamp> readTimes( const HdfFile &hdfFile )
return convertTimeData( times, dataTimeUnits );
}

static std::vector<int> readFace2Cells( const HdfFile &hdfFile, const std::string &flowAreaName, size_t *nFaces )
{
// First read face to node mapping
HdfGroup gGeom = openHdfGroup( hdfFile, "Geometry" );
HdfGroup gGeom2DFlowAreas = openHdfGroup( gGeom, "2D Flow Areas" );
HdfGroup gArea = openHdfGroup( gGeom2DFlowAreas, flowAreaName );
HdfDataset dsFace2Cells = openHdfDataset( gArea, "Faces Cell Indexes" );

std::vector<hsize_t> fdims = dsFace2Cells.dims();
std::vector<int> face2Cells = dsFace2Cells.readArrayInt(); //2x nFaces

*nFaces = fdims[0];
return face2Cells;
}

void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
const HdfGroup &rootGroup,
const std::vector<size_t> &areaElemStartIndex,
Expand All @@ -226,16 +211,14 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
const std::vector<RelativeTimestamp> &times,
const DateTime &referenceTime )
{
double eps = std::numeric_limits<double>::min();

std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
name(),
mMesh.get(),
mFileName,
datasetName
);
group->setDataLocation( MDAL_DataLocation::DataOnFaces );
group->setIsScalar( true );
group->setIsScalar( false );
group->setReferenceTime( referenceTime );

std::vector<std::shared_ptr<MDAL::MemoryDataset2D>> datasets;
Expand All @@ -253,34 +236,99 @@ void MDAL::DriverHec2D::readFaceOutput( const HdfFile &hdfFile,
{
std::string flowAreaName = flowAreaNames[nArea];

size_t nFaces;
std::vector<int> face2Cells = readFace2Cells( hdfFile, flowAreaName, &nFaces );

HdfGroup gFlowAreaRes = openHdfGroup( rootGroup, flowAreaName );
HdfDataset dsVals = openHdfDataset( gFlowAreaRes, rawDatasetName );
std::vector<float> vals = dsVals.readArray();

HdfGroup gGeom = openHdfGroup( hdfFile, "Geometry" );
HdfGroup gGeom2DFlowAreas = openHdfGroup( gGeom, "2D Flow Areas" );
HdfGroup gArea = openHdfGroup( gGeom2DFlowAreas, flowAreaName );
HdfDataset dsCellFaceInfo = openHdfDataset( gArea, "Cells Face and Orientation Info" );
std::vector<hsize_t> celldims = dsCellFaceInfo.dims();
std::vector<int> cellFaceInfo = dsCellFaceInfo.readArrayInt();

HdfDataset dsCellFaceOrValues = openHdfDataset( gArea, "Cells Face and Orientation Values" );
std::vector<int> cellFaceOrValues = dsCellFaceOrValues.readArrayInt();

HdfDataset dsFacePointIndex = openHdfDataset( gArea, "Faces FacePoint Indexes" );
std::vector<hsize_t> fdims = dsFacePointIndex.dims();
size_t nFaces = static_cast<size_t>( fdims[0] );
std::vector<int> facePointIndex = dsFacePointIndex.readArrayInt();

HdfDataset dsCoords = openHdfDataset( gArea, "FacePoints Coordinate" );
std::vector<hsize_t> cdims = dsCoords.dims();
std::vector<double> coords = dsCoords.readArrayDouble(); //2xnNodes matrix in array

for ( size_t tidx = 0; tidx < times.size(); ++tidx )
{
std::shared_ptr<MDAL::MemoryDataset2D> dataset = datasets[tidx];
double *values = dataset->values();

for ( size_t i = 0; i < nFaces; ++i )
for ( size_t cell_idx = 0; cell_idx < celldims.at( 0 ); ++cell_idx )
{
size_t idx = tidx * nFaces + i;
double val = static_cast<double>( vals[idx] ); // This is value on face!

if ( !std::isnan( val ) && fabs( val ) > eps ) //not nan and not 0
size_t cell_mdal_idx = cell_idx + areaElemStartIndex[nArea];
double valx = 0;
double valy = 0;
size_t consideredValueCount = 0;
size_t firstPosition = static_cast<size_t>( cellFaceInfo[cell_idx * 2] );
size_t faceCount = static_cast<size_t>( cellFaceInfo[cell_idx * 2 + 1] );
for ( size_t f = 0; f < faceCount; ++f )
{
for ( size_t c = 0; c < 2; ++c )
//get face indexes
size_t faceIndex1 = static_cast<size_t>( cellFaceOrValues[( firstPosition + f ) * 2] );
size_t faceIndex2 = static_cast<size_t>( cellFaceOrValues[( firstPosition + ( f + 1 ) % faceCount ) * 2] );
double val1 = static_cast<double>( vals[faceIndex1 + tidx * nFaces] );
double val2 = static_cast<double>( vals[faceIndex2 + tidx * nFaces] );
if ( std::isnan( val1 ) || std::isnan( val2 ) )
continue;

size_t indexPoint11 = facePointIndex[faceIndex1 * 2];
size_t indexPoint12 = facePointIndex[faceIndex1 * 2 + 1];
size_t indexPoint21 = facePointIndex[faceIndex2 * 2];
size_t indexPoint22 = facePointIndex[faceIndex2 * 2 + 1];
bool commonIndex = ( indexPoint11 == indexPoint21 ||
indexPoint11 == indexPoint22 ||
indexPoint12 == indexPoint21 ||
indexPoint12 == indexPoint22 );
if ( !commonIndex )
{
size_t cell_idx = static_cast<size_t>( face2Cells[2 * i + c] ) + areaElemStartIndex[nArea];
// Take just maximum
if ( std::isnan( values[cell_idx] ) || values[cell_idx] < val )
{
values[cell_idx] = val;
}
// should not happen, but better to prevent
continue;
}

double dx1 = coords[indexPoint11 * 2] - coords[indexPoint12 * 2];
double dy1 = coords[indexPoint11 * 2 + 1] - coords[indexPoint12 * 2 + 1];
double dx2 = coords[indexPoint21 * 2] - coords[indexPoint22 * 2];
double dy2 = coords[indexPoint21 * 2 + 1] - coords[indexPoint22 * 2 + 1];
double l1 = sqrt( dx1 * dx1 + dy1 * dy1 );
double l2 = sqrt( dx2 * dx2 + dy2 * dy2 );
if ( l1 == 0 || l2 == 0 )
{
continue;
}
double nx1 = -dy1 / l1;
double ny1 = dx1 / l1;
double nx2 = -dy2 / l2;
double ny2 = dx2 / l2;

double deter = nx1 * ny2 - nx2 * ny1;
if ( deter == 0 ) //colinear face, forbidden by hecras, but better to prevent
continue;
valx += ( ny2 * val1 - ny1 * val2 ) / deter;
valy += ( nx1 * val2 - nx2 * val1 ) / deter;
consideredValueCount++;
}
if ( consideredValueCount != 0 )
{
valx /= consideredValueCount;
valy /= consideredValueCount;
values[cell_mdal_idx * 2] = valx;
values[cell_mdal_idx * 2 + 1] = valy;
}
else
{
values[cell_mdal_idx * 2] = std::numeric_limits<double>::quiet_NaN();
values[cell_mdal_idx * 2 + 1] = std::numeric_limits<double>::quiet_NaN();
}
}
}
Expand All @@ -302,15 +350,15 @@ void MDAL::DriverHec2D::readFaceResults( const HdfFile &hdfFile,
// UNSTEADY
HdfGroup flowGroup = get2DFlowAreasGroup( hdfFile, "Unsteady Time Series" );
MDAL::DateTime referenceDateTime = readReferenceDateTime( hdfFile );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Shear Stress", "Face Shear Stress", mTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Velocity", "Face Velocity", mTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Shear Stress", "Shear Stress", mTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Face Velocity", "Velocity", mTimes, referenceDateTime );

// SUMMARY
flowGroup = get2DFlowAreasGroup( hdfFile, "Summary Output" );
std::vector<MDAL::RelativeTimestamp> dummyTimes( 1, MDAL::RelativeTimestamp() );

readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Shear Stress", "Face Shear Stress/Maximums", dummyTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Velocity", "Face Velocity/Maximums", dummyTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Shear Stress", "Shear Stress/Maximums", dummyTimes, referenceDateTime );
readFaceOutput( hdfFile, flowGroup, areaElemStartIndex, flowAreaNames, "Maximum Face Velocity", "Velocity/Maximums", dummyTimes, referenceDateTime );
}


Expand Down

0 comments on commit a34ff39

Please sign in to comment.