Skip to content

Commit

Permalink
fix for quad and triangle mesh. optimize storing of large datasets [u…
Browse files Browse the repository at this point in the history 
…grid]
  • Loading branch information
@PeterPetrik
PeterPetrik committed Oct 31, 2019
1 parent afc6901 commit cb97906
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 51 deletions.
54 changes: 44 additions & 10 deletions mdal/frmts/mdal_ugrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@

#include <netcdf.h>
#include <assert.h>
#include <algorithm>
#include <cmath>

MDAL::DriverUgrid::DriverUgrid()
: DriverCF(
Expand Down Expand Up @@ -470,7 +472,8 @@ void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
mNcFile.putAttrStr( mesh2dNodeZId, "standard_name", "altitude" );
mNcFile.putAttrStr( mesh2dNodeZId, "long_name", "z-coordinate of mesh nodes" );
mNcFile.putAttrStr( mesh2dNodeZId, "grid_mapping", "projected_coordinate_system" );
mNcFile.putAttrDouble( mesh2dNodeZId, "_FillValue", -999.0 );
double fillNodeZCoodVal = -999.0;
mNcFile.putAttrDouble( mesh2dNodeZId, "_FillValue", fillNodeZCoodVal );

// Faces 2D Variable
int mesh2FaceNodesId_dimIds [] { dimFaceCountId, dimMaxNodesPerFaceId };
Expand All @@ -480,7 +483,8 @@ void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
mNcFile.putAttrStr( mesh2FaceNodesId, "location", "face" );
mNcFile.putAttrStr( mesh2FaceNodesId, "long_name", "Mapping from every face to its corner nodes (counterclockwise)" );
mNcFile.putAttrInt( mesh2FaceNodesId, "start_index", 0 );
mNcFile.putAttrInt( mesh2FaceNodesId, "_FillValue", -999 );
int fillFace2DVertexValue = -999;
mNcFile.putAttrInt( mesh2FaceNodesId, "_FillValue", fillFace2DVertexValue );

// Projected Coordinate System
int pcsId = mNcFile.defineVar( "projected_coordinate_system", NC_INT, 0, nullptr );
Expand Down Expand Up @@ -518,22 +522,25 @@ void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
const size_t bufferSize = std::min( mesh->verticesCount(), maxBufferSize );
const size_t verticesCoordCount = bufferSize * 3;

std::unique_ptr<double[]> verticesCoordinates( new double[verticesCoordCount] );
std::vector<double> verticesCoordinates( verticesCoordCount );
std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();

size_t vertexIndex = 0;
size_t vertexFileIndex = 0;
while ( vertexIndex < mesh->verticesCount() )
{
size_t verticesRead = vertexIterator->next( bufferSize, verticesCoordinates.get() );
size_t verticesRead = vertexIterator->next( bufferSize, verticesCoordinates.data() );
if ( verticesRead == 0 )
break;

for ( size_t i = 0; i < verticesRead; i++ )
{
mNcFile.putDataDouble( mesh2dNodeXId, vertexFileIndex, verticesCoordinates[3 * i] );
mNcFile.putDataDouble( mesh2dNodeYId, vertexFileIndex, verticesCoordinates[3 * i + 1] );
mNcFile.putDataDouble( mesh2dNodeZId, vertexFileIndex, verticesCoordinates[3 * i + 2] );
if ( std::isnan( verticesCoordinates[3 * i + 2] ) )
mNcFile.putDataDouble( mesh2dNodeZId, vertexFileIndex, fillNodeZCoodVal );
else
mNcFile.putDataDouble( mesh2dNodeZId, vertexFileIndex, verticesCoordinates[3 * i + 2] );
vertexFileIndex++;
}
vertexIndex += verticesRead;
Expand All @@ -542,14 +549,41 @@ void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
// Write faces
std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
const size_t faceVerticesMax = mesh->faceVerticesMaximumCount();
const size_t facesCount = mesh->facesCount();
const size_t faceOffsetsBufferLen = std::min( facesCount, maxBufferSize );
const size_t vertexIndicesBufferLen = faceOffsetsBufferLen * faceVerticesMax;

int faceOffsets[1];
std::unique_ptr<int[]>vertexIndices( new int[faceVerticesMax] );
std::vector<int> faceOffsetsBuffer( faceOffsetsBufferLen );
std::vector<int> vertexIndicesBuffer( vertexIndicesBufferLen );

for ( size_t i = 0; i < mesh->facesCount(); ++i )
size_t faceIndex = 0;
while ( faceIndex < facesCount )
{
faceIterator->next( 1, faceOffsets, faceVerticesMax, vertexIndices.get() );
mNcFile.putDataArrayInt( mesh2FaceNodesId, i, faceVerticesMax, vertexIndices.get() );
size_t facesRead = faceIterator->next(
faceOffsetsBufferLen,
faceOffsetsBuffer.data(),
vertexIndicesBufferLen,
vertexIndicesBuffer.data() );
if ( facesRead == 0 )
break;

for ( size_t i = 0; i < facesRead; i++ )
{
std::vector<int> verticesFaceData( faceVerticesMax, fillFace2DVertexValue );
int startIndex = 0;
if ( i > 0 )
startIndex = faceOffsetsBuffer[ i - 1 ];
int endIndex = faceOffsetsBuffer[ i ];

size_t k = 0;
for ( int j = startIndex; j < endIndex; ++j )
{
int vertexIndex = vertexIndicesBuffer[ static_cast<size_t>( j ) ];
verticesFaceData[k++] = vertexIndex;
}
mNcFile.putDataArrayInt( mesh2FaceNodesId, faceIndex + i, faceVerticesMax, verticesFaceData.data() );
}
faceIndex += facesRead;
}

// Time values (not implemented)
Expand Down
38 changes: 38 additions & 0 deletions tests/mdal_testutils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,31 @@ bool compareVectors( const std::vector<double> &a, const std::vector<double> &b
return true;
}

bool compareMeshFrames( MeshH meshA, MeshH meshB )
{
// Vertices
int orignal_v_count = MDAL_M_vertexCount( meshA );
int saved_v_count = MDAL_M_vertexCount( meshB );
if ( orignal_v_count != saved_v_count ) return false;

std::vector<double> coordsA = getCoordinates( meshA, orignal_v_count );
std::vector<double> coordsB = getCoordinates( meshB, saved_v_count );
if ( !compareVectors( coordsA, coordsB ) )
return false;

// Faces
int orignal_f_count = MDAL_M_faceCount( meshA );
int saved_f_count = MDAL_M_faceCount( meshB );
if ( orignal_f_count != saved_f_count ) return false;

std::vector<int> verticesA = faceVertexIndices( meshA, orignal_f_count );
std::vector<int> verticesB = faceVertexIndices( meshB, saved_f_count );
if ( !compareVectors( verticesA, verticesB ) )
return false;

return true;
}

std::vector<double> getCoordinates( MeshH mesh, int verticesCount )
{
MeshVertexIteratorH iterator = MDAL_M_vertexIterator( mesh );
Expand Down Expand Up @@ -127,6 +152,19 @@ double getVertexZCoordinatesAt( MeshH mesh, int index )
return _getVertexCoordinatesAt( mesh, index, 2 );
}

std::vector<int> faceVertexIndices( MeshH mesh, int faceCount )
{
MeshFaceIteratorH iterator = MDAL_M_faceIterator( mesh );
int faceOffsetsBufferLen = faceCount + 1;
int vertexIndicesBufferLen = faceOffsetsBufferLen * MDAL_M_faceVerticesMaximumCount( mesh );
std::vector<int> faceOffsetsBuffer( static_cast<size_t>( faceOffsetsBufferLen ) );
std::vector<int> vertexIndicesBuffer( static_cast<size_t>( vertexIndicesBufferLen ) );
MDAL_FI_next( iterator, faceOffsetsBufferLen, faceOffsetsBuffer.data(),
vertexIndicesBufferLen, vertexIndicesBuffer.data() );
MDAL_FI_close( iterator );
return vertexIndicesBuffer;
}

int getFaceVerticesCountAt( MeshH mesh, int faceIndex )
{
MeshFaceIteratorH iterator = MDAL_M_faceIterator( mesh );
Expand Down
4 changes: 4 additions & 0 deletions tests/mdal_testutils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ double getVertexYCoordinatesAt( MeshH mesh, int index );
double getVertexZCoordinatesAt( MeshH mesh, int index );
int getFaceVerticesCountAt( MeshH mesh, int faceIndex );
int getFaceVerticesIndexAt( MeshH mesh, int faceIndex, int index );
std::vector<int> faceVertexIndices( MeshH mesh, int faceCount );

// Datasets
bool getActive( DatasetH dataset, int index );
Expand All @@ -37,4 +38,7 @@ double getValueY( DatasetH dataset, int index );
bool compareVectors( const std::vector<double> &a, const std::vector<double> &b );
bool compareVectors( const std::vector<int> &a, const std::vector<int> &b );

//! Same vertices (coords), faces and connectivity between them
bool compareMeshFrames( MeshH meshA, MeshH meshB );

#endif // MDAL_TESTUTILS_HPP
70 changes: 29 additions & 41 deletions tests/test_ugrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ TEST( MeshUgridTest, SaveDFlow11Manzese )

// Save the mesh as UGRID
std::string fileNameToSave = tmp_file( "/manzese_1d2d_small_map_saveTest.nc" );
std::string fileNameToSave2 = tmp_file( "/manzese_1d2d_small_map_saveTest2.nc" );
MDAL_SaveMesh( meshToSave, fileNameToSave.c_str(), "Ugrid" );
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );
Expand All @@ -35,57 +34,46 @@ TEST( MeshUgridTest, SaveDFlow11Manzese )
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

// Save again to make sure saved file can be loaded
MDAL_SaveMesh( savedMesh, fileNameToSave2.c_str(), "Ugrid" );
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

// Load second saved UGRID mesh
MeshH savedMesh2 = MDAL_LoadMesh( fileNameToSave2.c_str() );
EXPECT_NE( savedMesh2, nullptr );
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

// Compare second save with the original mesh

// Vertices
int orignal_v_count = MDAL_M_vertexCount( meshToSave );
int saved_v_count = MDAL_M_vertexCount( savedMesh2 );
EXPECT_EQ( orignal_v_count, saved_v_count );
// Compare saved with the original mesh
EXPECT_TRUE( compareMeshFrames( meshToSave, savedMesh ) );

// Faces
int orignal_f_count = MDAL_M_faceCount( meshToSave );
int saved_f_count = MDAL_M_faceCount( savedMesh2 );
EXPECT_EQ( orignal_f_count, saved_f_count );
// Close meshed and delete all the files
MDAL_CloseMesh( meshToSave );
MDAL_CloseMesh( savedMesh );
std::remove( fileNameToSave.c_str() );
}

// Test face 1
int orignal_f_v_count = getFaceVerticesCountAt( meshToSave, 1 );
int saved_f_v_count = getFaceVerticesCountAt( savedMesh2, 1 );
EXPECT_EQ( orignal_f_v_count, saved_f_v_count ); //quad
TEST( MeshUgridTest, SaveQuadAndTriangle )
{
//test driver capability
EXPECT_TRUE( MDAL_DR_SaveMeshCapability( MDAL_driverFromName( "Ugrid" ) ) );

// Test verticies
int orignal_f_v = getFaceVerticesIndexAt( meshToSave, 1, 0 );
int saved_f_v = getFaceVerticesIndexAt( savedMesh2, 1, 0 );
EXPECT_EQ( orignal_f_v, saved_f_v );
// Open mesh with both triangle and quad
std::string path = test_file( "/2dm/quad_and_triangle.2dm" );
MeshH meshToSave = MDAL_LoadMesh( path.c_str() );
EXPECT_NE( meshToSave, nullptr );
MDAL_Status s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

orignal_f_v = getFaceVerticesIndexAt( meshToSave, 1, 1 );
saved_f_v = getFaceVerticesIndexAt( savedMesh2, 1, 1 );
EXPECT_EQ( orignal_f_v, saved_f_v );
// Save the mesh as UGRID
std::string fileNameToSave = tmp_file( "/quad_and_triangle_saveTest.nc" );
MDAL_SaveMesh( meshToSave, fileNameToSave.c_str(), "Ugrid" );
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

orignal_f_v = getFaceVerticesIndexAt( meshToSave, 1, 2 );
saved_f_v = getFaceVerticesIndexAt( savedMesh2, 1, 2 );
EXPECT_EQ( orignal_f_v, saved_f_v );
// Load saved UGRID mesh
MeshH savedMesh = MDAL_LoadMesh( fileNameToSave.c_str() );
EXPECT_NE( savedMesh, nullptr );
s = MDAL_LastStatus();
ASSERT_EQ( MDAL_Status::None, s );

orignal_f_v = getFaceVerticesIndexAt( meshToSave, 1, 3 );
saved_f_v = getFaceVerticesIndexAt( savedMesh2, 1, 3 );
EXPECT_EQ( orignal_f_v, saved_f_v );
// Compare saved with the original mesh
EXPECT_TRUE( compareMeshFrames( meshToSave, savedMesh ) );

// Close meshed and delete all the files
MDAL_CloseMesh( meshToSave );
MDAL_CloseMesh( savedMesh );
MDAL_CloseMesh( savedMesh2 );
std::remove( fileNameToSave.c_str() );
std::remove( fileNameToSave2.c_str() );
}

TEST( MeshUgridTest, DFlow11Manzese )
Expand Down

0 comments on commit cb97906

Please sign in to comment.