From a390b30fcb126c78ab86104cb30d26b8ae672cfb Mon Sep 17 00:00:00 2001 From: Achim Gsell Date: Wed, 20 Aug 2008 14:29:26 +0000 Subject: [PATCH] mapping of tets and triangles to global and local id added --- .gitattributes | 5 + src/h5/Makefile.am | 15 +- src/h5/attribs.c | 6 +- src/h5/errorhandling.c | 2 +- src/h5/errorhandling_private.h | 10 +- src/h5/general.c | 24 +- src/h5/h5_core.h | 2 + src/h5/h5_private.h | 13 +- src/h5/h5_types.h | 7 + src/h5/maps.c | 42 +- src/h5/maps.h | 14 +- src/h5/readwrite.c | 4 +- src/h5/t_errorhandling.c | 41 ++ src/h5/t_errorhandling_private.h | 65 +++ src/h5/t_map.c | 657 +++++++++++++++++++++++++++++++ src/h5/t_map.h | 72 ++++ src/h5/t_map_private.h | 18 + src/h5/t_openclose.c | 4 - src/h5/t_readwrite.c | 352 +++++++---------- src/h5/t_readwrite.h | 17 +- src/h5/u_readwrite.c | 20 +- 21 files changed, 1108 insertions(+), 282 deletions(-) create mode 100644 src/h5/t_errorhandling.c create mode 100644 src/h5/t_errorhandling_private.h create mode 100644 src/h5/t_map.c create mode 100644 src/h5/t_map.h create mode 100644 src/h5/t_map_private.h diff --git a/.gitattributes b/.gitattributes index 5f6a7d8..077261f 100644 --- a/.gitattributes +++ b/.gitattributes @@ -387,6 +387,11 @@ src/h5/maps.c -text src/h5/maps.h -text src/h5/readwrite.c -text src/h5/readwrite.h -text +src/h5/t_errorhandling.c -text +src/h5/t_errorhandling_private.h -text +src/h5/t_map.c -text +src/h5/t_map.h -text +src/h5/t_map_private.h -text src/h5/t_openclose.c -text src/h5/t_openclose.h -text src/h5/t_readwrite.c -text diff --git a/src/h5/Makefile.am b/src/h5/Makefile.am index 8f5dfdd..15c5efb 100644 --- a/src/h5/Makefile.am +++ b/src/h5/Makefile.am @@ -47,8 +47,11 @@ include_HEADERS = \ h5_types.h \ attribs.h \ errorhandling.h \ + errorhandling_private.h \ maps.h \ readwrite.h \ + t_errorhandling_private.h \ + t_map.h \ t_openclose.h \ t_readwrite.h \ u_readwrite.h @@ -58,14 +61,16 @@ EXTRA_HEADERS = # Listing of sources libH5_a_SOURCES = \ - general.c \ attribs.c \ errorhandling.c \ + general.c \ maps.c \ - t_openclose.c \ readwrite.c \ - u_readwrite.c \ - t_readwrite.c + t_errorhandling.c \ + t_map.c \ + t_openclose.c \ + t_readwrite.c \ + u_readwrite.c all: libH5.a @@ -78,7 +83,7 @@ libH5.a: $(libH5_a_OBJECTS) $(libH5_a_OBJECTS): h5_core.h -errorhandling.o: errorhandling.c +#errorhandling.o: errorhandling.c clean: $(RM) -f *~ *.o *.a *.so diff --git a/src/h5/attribs.c b/src/h5/attribs.c index 2e7a3c3..d7e3f73 100644 --- a/src/h5/attribs.c +++ b/src/h5/attribs.c @@ -51,7 +51,7 @@ h5_read_attrib ( herr = H5Aclose ( attrib_id ); if ( herr < 0 ) return HANDLE_H5A_CLOSE_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -88,7 +88,7 @@ h5_write_attrib ( herr = H5Sclose ( space_id ); if ( herr < 0 ) return HANDLE_H5S_CLOSE_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -139,7 +139,7 @@ h5_get_attrib_info ( herr = H5Aclose ( attrib_id); if ( herr < 0 ) return HANDLE_H5A_CLOSE_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t diff --git a/src/h5/errorhandling.c b/src/h5/errorhandling.c index ed3230d..af50914 100644 --- a/src/h5/errorhandling.c +++ b/src/h5/errorhandling.c @@ -11,7 +11,7 @@ #include "H5Part.h" h5_error_handler _err_handler = h5_report_errorhandler; -h5_err_t _h5part_errno = H5PART_SUCCESS; +h5_err_t _h5part_errno = H5_SUCCESS; h5_id_t _debug = 0; static char *__funcname = "NONE"; diff --git a/src/h5/errorhandling_private.h b/src/h5/errorhandling_private.h index 51a5237..78293d8 100644 --- a/src/h5/errorhandling_private.h +++ b/src/h5/errorhandling_private.h @@ -1,6 +1,10 @@ #ifndef __ERRORHANDLING_PRIVATE_H #define __ERRORHANDLING_PRIVATE_H +extern h5part_error_handler _err_handler; +extern h5_err_t _h5part_errno; +extern h5_id_t _debug; + #define HANDLE_H5_ERR_NOT_IMPLEMENTED \ (*h5_get_errorhandler()) ( \ h5_get_funcname(), \ @@ -110,12 +114,6 @@ "No entry with index %lld and type %d in group %s!", \ (long long)idx, type, group_name ); -#define HANDLE_H5T_GID_NOT_EXIST_ERR( name, id ) \ - (*h5_get_errorhandler()) ( \ - h5_get_funcname(), \ - H5_ERR_NOENTRY, \ - "Global %s id %ld does not exist!", \ - name, (long)id ); #define HANDLE_H5_UNDEF_MESH_ERR \ (*h5_get_errorhandler()) ( \ diff --git a/src/h5/general.c b/src/h5/general.c index 35dc1dc..7d7dfd9 100644 --- a/src/h5/general.c +++ b/src/h5/general.c @@ -12,10 +12,6 @@ #include "H5Part.h" #include "H5Block.h" -extern h5part_error_handler _err_handler; -extern h5part_int64_t _h5part_errno; -extern unsigned _debug; - /*! \ingroup h5block_private @@ -23,7 +19,7 @@ extern unsigned _debug; Check whether \c f points to a valid file handle. - \return H5PART_SUCCESS or error code + \return H5_SUCCESS or error code */ h5part_int64_t h5_check_filehandle ( @@ -36,7 +32,7 @@ h5_check_filehandle ( return HANDLE_H5_BADFD_ERR; if ( f->block == NULL ) return HANDLE_H5_BADFD_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } @@ -62,7 +58,7 @@ _init ( void ) { if ( r5 < 0 ) return H5PART_ERR_INIT; } __init = 1; - return H5PART_SUCCESS; + return H5_SUCCESS; } /*! @@ -133,7 +129,7 @@ _h5b_open_file ( b->field_group_id = -1; b->have_layout = 0; - return H5PART_SUCCESS; + return H5_SUCCESS; } @@ -150,7 +146,7 @@ h5_open_file ( HANDLE_H5_INIT_ERR; return NULL; } - _h5part_errno = H5PART_SUCCESS; + _h5part_errno = H5_SUCCESS; h5_file *f = NULL; f = (h5_file*) malloc( sizeof (h5_file) ); @@ -301,7 +297,7 @@ h5_open_file ( De-initialize H5Block internal structure. Open HDF5 objects are closed and allocated memory freed. - \return H5PART_SUCCESS or error code + \return H5_SUCCESS or error code */ static h5part_int64_t _h5u_close_file ( @@ -338,7 +334,7 @@ _h5u_close_file ( De-initialize H5Block internal structure. Open HDF5 objects are closed and allocated memory freed. - \return H5PART_SUCCESS or error code + \return H5_SUCCESS or error code */ static h5part_int64_t _h5b_close_file ( @@ -371,7 +367,7 @@ _h5b_close_file ( free ( f->block ); f->block = NULL; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -379,7 +375,7 @@ h5_close_file ( h5_file *f ) { herr_t r = 0; - _h5part_errno = H5PART_SUCCESS; + _h5part_errno = H5_SUCCESS; CHECK_FILEHANDLE ( f ); @@ -439,7 +435,7 @@ h5_define_stepname_fmt ( } f->width_step_idx = (int)width; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5_err_t diff --git a/src/h5/h5_core.h b/src/h5/h5_core.h index 3ae880d..7892c17 100644 --- a/src/h5/h5_core.h +++ b/src/h5/h5_core.h @@ -1,3 +1,4 @@ + #ifndef __H5_CORE_H #define __H5_CORE_H @@ -55,6 +56,7 @@ h5_get_step ( #include "errorhandling.h" #include "maps.h" #include "readwrite.h" +#include "t_map.h" #include "t_openclose.h" #include "t_readwrite.h" #include "u_readwrite.h" diff --git a/src/h5/h5_private.h b/src/h5/h5_private.h index cfda1b1..20d9198 100644 --- a/src/h5/h5_private.h +++ b/src/h5/h5_private.h @@ -2,13 +2,21 @@ #define __H5_PRIVATE_H #include "errorhandling_private.h" +#include "t_map_private.h" +#include "t_errorhandling_private.h" -#define H5B_CONTAINER_GRPNAME "Block" +#define H5PART_GROUPNAME_STEP "Step" -#define H5T_CONTAINER_GRPNAME "Topo" +#define H5B_CONTAINER_GRPNAME "Block" + +#define H5T_CONTAINER_GRPNAME "Topo" #define H5BLOCK_GROUPNAME_BLOCK H5B_CONTAINER_GRPNAME +#define H5_TET_MASK ( (h5_id_t)(-1) >> 3 ) +#define _h5t_build_triangle_id( idx, entity_id ) \ + ( (idx << (sizeof(entity_id)*8 - 3)) | (entity_id & H5_TET_MASK)) + #define SET_FNAME( fname ) h5_set_funcname( fname ); #define CHECK_FILEHANDLE( f ) \ @@ -37,7 +45,6 @@ "Internal error: step_gid <= 0."); -#define H5PART_GROUPNAME_STEP "Step" /*! The functions declared here are not part of the API, but may be used diff --git a/src/h5/h5_types.h b/src/h5/h5_types.h index 4e7f93c..539145c 100644 --- a/src/h5/h5_types.h +++ b/src/h5/h5_types.h @@ -54,6 +54,8 @@ enum h5_mesh_types { /* enum with number of vertices(!) */ TETRAHEDRAL_MESH = 4 }; +#define H5_MAX_VERTICES_PER_ENTITY TETRAHEDRAL_MESH + struct h5_vertex { /* 32Byte */ h5_id_t id; h5_id_t unused; /* for right alignment */ @@ -209,6 +211,11 @@ struct h5t_fdata { h5_size_t *num_entities_on_level; struct idmap map_entity_g2l;/* map global id to local id */ + /* array with geometrically sorted local entitiy ids + [0]: 0,1,2,3 sorted + [1]: 1,0,2,3 sorted */ + struct smap sorted_lentities[H5_MAX_VERTICES_PER_ENTITY]; + h5_id_t last_retrieved_entity_id; h5_id_t last_stored_entity_id; diff --git a/src/h5/maps.c b/src/h5/maps.c index d0e9cc9..6e16d73 100644 --- a/src/h5/maps.c +++ b/src/h5/maps.c @@ -5,17 +5,33 @@ #include "h5_core.h" #include "h5_private.h" +h5_err_t +_h5_alloc_smap ( + struct smap *map, + h5_size_t size + ) { + int new = ( map == NULL ); + map->items = realloc ( map->items, size * sizeof ( map->items[0] ) ); + if ( map->items == NULL ) { + return HANDLE_H5_NOMEM_ERR; + } + map->size = size; + if ( new ) map->num_items = 0; + return H5_SUCCESS; +} h5_err_t _h5_alloc_idmap ( struct idmap *map, h5_size_t size ) { + int new = ( map == NULL ); map->items = realloc ( map->items, size * sizeof ( map->items[0] ) ); if ( map->items == NULL ) { return HANDLE_H5_NOMEM_ERR; } map->size = size; + if ( new ) map->num_items = 0; return H5_SUCCESS; } @@ -51,7 +67,7 @@ _h5_insert_idmap ( \ingroup h5_core - binary search in simple map + binary search in id map. \return index in array if found, othwise \c -(result+1) is the index where \c value must be inserted. @@ -78,15 +94,19 @@ _h5_search_idmap ( return -(low+1); // not found } -h5_id_t -h5t_map_vertex_id_global2local ( - h5_file *f, - h5_id_t global_id +int +_cmp_idmap ( + const void *id1, + const void *id2 ) { - struct h5t_fdata *t = &f->t; - - h5_id_t local_id = _h5_search_idmap ( &t->map_vertex_g2l, global_id ); - if ( local_id < 0 ) - return HANDLE_H5T_GID_NOT_EXIST_ERR ( "vertex", global_id ); - return local_id; + + return *(h5_id_t*)id1 - *(h5_id_t*)id2; +} + +h5_err_t +_h5_sort_idmap ( + struct idmap *map + ) { + qsort ( map->items, map->num_items, sizeof(map->items[0]), _cmp_idmap ); + return H5_SUCCESS; } diff --git a/src/h5/maps.h b/src/h5/maps.h index 344102f..0638fe9 100644 --- a/src/h5/maps.h +++ b/src/h5/maps.h @@ -1,6 +1,12 @@ #ifndef __MAPS_H #define __MAPS_H +h5_err_t +_h5_alloc_smap ( + struct smap *map, + h5_size_t size + ); + h5_err_t _h5_alloc_idmap ( struct idmap *map, @@ -20,9 +26,9 @@ _h5_search_idmap ( h5_id_t value ); -h5_id_t -h5t_map_vertex_id_global2local ( - h5_file *f, - h5_id_t global_id +h5_err_t +_h5_sort_idmap ( + struct idmap *map ); + #endif diff --git a/src/h5/readwrite.c b/src/h5/readwrite.c index b8a8f4a..60ac283 100644 --- a/src/h5/readwrite.c +++ b/src/h5/readwrite.c @@ -51,7 +51,7 @@ h5_write_data ( f->empty = 0; - return H5PART_SUCCESS; + return H5_SUCCESS; } @@ -178,7 +178,7 @@ h5_get_object_name ( if ( herr == 0 ) HANDLE_H5_NOENTRY_ERR( group_name, type, idx ); - return H5PART_SUCCESS; + return H5_SUCCESS; } h5_int64_t diff --git a/src/h5/t_errorhandling.c b/src/h5/t_errorhandling.c new file mode 100644 index 0000000..937dc21 --- /dev/null +++ b/src/h5/t_errorhandling.c @@ -0,0 +1,41 @@ +#include +#include +#include /* va_arg - System dependent ?! */ +#include +#include +#include +#include + +#include "h5_core.h" +#include "h5_private.h" + +h5_err_t +_h5t_handle_get_global_entity_id_err ( + h5_file *f, + const h5_id_t * const global_vids + ) { + struct h5t_fdata *t = &f->t; + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: + return _h5t_handle_get_global_tet_id_err ( global_vids ); + case TRIANGLE_MESH: + return _h5t_handle_get_global_tri_id_err ( global_vids ); + } + return -1; +} + +h5_err_t +_h5t_handle_get_local_entity_id_err ( + h5_file *f, + const h5_id_t * const local_vids + ) { + struct h5t_fdata *t = &f->t; + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: + return _h5t_handle_get_local_tet_id_err ( local_vids ); + case TRIANGLE_MESH: + return _h5t_handle_get_local_triangle_id_err ( local_vids ); + } + return -1; +} + diff --git a/src/h5/t_errorhandling_private.h b/src/h5/t_errorhandling_private.h new file mode 100644 index 0000000..e5f1107 --- /dev/null +++ b/src/h5/t_errorhandling_private.h @@ -0,0 +1,65 @@ +#ifndef __T_ERRORHANDLING_H +#define __T_ERRORHANDLING_H + +h5_err_t +_h5t_handle_get_global_entity_id_err ( + h5_file *f, + const h5_id_t * const global_vids + ); + +h5_err_t +_h5t_handle_get_local_entity_id_err ( + h5_file *f, + const h5_id_t * const local_vids + ); + +#define H5T_HANDLE_GID_NOT_EXIST_ERR( name, id ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Global %s id %ld does not exist!", \ + name, (long)id ); + +#define _h5t_handle_global_id_not_exist_err( name, id ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "%s with global id %ld does not exist!", \ + name, (long)id ); + +#define _h5t_handle_get_global_tet_id_err( vids ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Tetrahedron with global vertex ids (%d,%d,%d,%d) doesn't exist!", \ + vids[0], vids[1], vids[2], vids[3] ); + +#define _h5t_handle_get_global_tri_id_err( vids ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Triangle with global vertex ids (%d,%d,%d) doesn't exist!", \ + vids[0], vids[1], vids[2] ); + +#define _h5t_handle_get_local_tet_id_err( vids ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Tetrahedron with local vertex ids (%d,%d,%d,%d) doesn't exist!", \ + vids[0], vids[1], vids[2], vids[3] ); + +#define _h5t_handle_get_local_triangle_id_err( vids ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Triangle with local vertex ids (%d,%d,%d) doesn't exist!", \ + vids[0], vids[1], vids[2] ); + +#define _h5t_handle_get_global_triangle_id_err( vids ) \ + (*h5_get_errorhandler()) ( \ + h5_get_funcname(), \ + H5_ERR_NOENTRY, \ + "Triangle with global vertex ids (%d,%d,%d) doesn't exist!", \ + vids[0], vids[1], vids[2] ); + +#endif diff --git a/src/h5/t_map.c b/src/h5/t_map.c new file mode 100644 index 0000000..a36fc10 --- /dev/null +++ b/src/h5/t_map.c @@ -0,0 +1,657 @@ +#include + +#include "h5/h5_core.h" +#include "h5/h5_private.h" +#include "h5/t_map.h" + +/*! + Returns the local vertex id of the i-th vertex of an entity. For triangles + i is in [0,1,2], for tetraheda i is in [0,1,2,3]. +*/ + +static h5_id_t +_get_local_vertex_id_of_entity ( + h5_file * f, + int ith_vertex, + h5_id_t local_eid + ) { + struct h5t_fdata *t = &f->t; + + h5_id_t local_vid = -1; + h5_id_t global_vid = -1; + + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: { + local_vid = t->lentities.tets[local_eid].vertex_ids[ith_vertex]; + if ( local_vid == -1 ) { + global_vid = + t->entities.tets[local_eid].vertex_ids[ith_vertex]; + local_vid = _h5_search_idmap ( + &t->map_vertex_g2l, global_vid ); + t->lentities.tets[local_eid].vertex_ids[ith_vertex] = + local_vid; + } + break; + } + case TRIANGLE_MESH: { + local_vid = t->lentities.tris[local_eid].vertex_ids[ith_vertex]; + if ( local_vid == -1 ) { + global_vid = + t->entities.tris[local_eid].vertex_ids[ith_vertex]; + local_vid = _h5_search_idmap ( + &t->map_vertex_g2l, global_vid ); + t->lentities.tris[local_eid].vertex_ids[ith_vertex] = + local_vid; + } + break; + } + } + return local_vid; +} + +/*! + Return the coordinates of the i-th vertex of the entity given by its local + entity id. For triangles \c i is in \c [0,1,2], for tetraheda \c i is in + \c [0,1,2,3]. +*/ +static h5_float64_t* +_get_vertex_of_entity ( + h5_file * f, + int ith_vertex, + h5_id_t local_entity_id + ) { + struct h5t_fdata *t = &f->t; + + h5_id_t local_vid = _get_local_vertex_id_of_entity ( + f, ith_vertex, local_entity_id ); + if ( local_vid == -1 ) + return NULL; + return t->vertices[local_vid].P; +} + +/*! + Compare to vertices given by their coordinates +*/ + +static int +_cmp_vertices ( + h5_float64_t *P1, + h5_float64_t *P2 + ) { + int i; + for ( i = 0; i < 3; i++ ) { + if ( P1[i] < P2[i] ) return -1; + else if (P1[i] > P2[i] )return 1; + } + return 0; +} + +/*! + compare two entities given by their local vertex ids +*/ +static int +_vcmp_entities ( + h5_file * f, + h5_id_t *e1, + h5_id_t *e2 + ) { + struct h5t_fdata *t = &f->t; + int i; + + for ( i = 0; i < t->mesh_type; i++ ) { + int r = _cmp_vertices ( + t->vertices[e1[i]].P, + t->vertices[e2[i]].P ); + if ( r < 0 ) return -1; + else if ( r > 0 ) return 1; + } + return 0; +} + +/*! + compare two entities given by their local id +*/ +static int +_cmp_entities ( + h5_file * f, + h5_id_t local_eid1, + h5_id_t local_eid2 + ) { + struct h5t_fdata *t = &f->t; + int i; + for ( i = 0; i < t->mesh_type; i++ ) { + int r = _cmp_vertices ( + _get_vertex_of_entity ( f, i, local_eid1 ), + _get_vertex_of_entity ( f, i, local_eid2 ) ); + if ( r < 0 ) return -1; + else if ( r > 0 ) return 1; + } + return 0; +} + +static int +_cmp_entities1 ( + h5_file * f, + h5_id_t local_eid1, + h5_id_t local_eid2 + ) { + struct h5t_fdata *t = &f->t; + int imap[] = { 1, 0, 2, 3 }; + int i; + for ( i = 0; i < t->mesh_type; i++ ) { + int r = _cmp_vertices ( + _get_vertex_of_entity ( f, imap[i], local_eid1 ), + _get_vertex_of_entity ( f, imap[i], local_eid2 ) ); + if ( r < 0 ) return -1; + else if ( r > 0 ) return 1; + } + return 0; +} + +/* + The re-rentrant version of qsort(3) is not available on many systems, thus ... +*/ +static h5_file *_f; +static int +_qsort_cmp_entities0 ( const void* _local_eid1, const void* _local_eid2 ) { + h5_id_t local_eid1 = *(h5_id_t*)_local_eid1; + h5_id_t local_eid2 = *(h5_id_t*)_local_eid2; + return _cmp_entities ( _f, local_eid1, local_eid2 ); +} + +static int +_qsort_cmp_entities1 ( const void* _local_eid1, const void* _local_eid2 ) { + h5_id_t local_eid1 = *(h5_id_t*)_local_eid1; + h5_id_t local_eid2 = *(h5_id_t*)_local_eid2; + return _cmp_entities1 ( _f, local_eid1, local_eid2 ); +} + +/*! + Sort entities geometrically. +*/ +static h5_err_t +_sort_entities ( + h5_file *f + ) { + + struct h5t_fdata *t = &f->t; + h5_size_t num_entities = h5t_get_num_entities ( f ); + if ( num_entities < 0 ) return num_entities; + + int k; + h5_id_t i; + for ( k = 0; k < 2; k++ ) { + h5_err_t h5err = _h5_alloc_smap ( + &t->sorted_lentities[k], num_entities ); + if ( h5err < 0 ) return h5err; + for ( i = 0; i < num_entities; i++ ) { + t->sorted_lentities[k].items[i] = i; + } + t->sorted_lentities[k].num_items = num_entities; + } + + + _f = f; + qsort ( t->sorted_lentities[0].items, + num_entities, + sizeof(t->sorted_lentities[0].items[0]), + _qsort_cmp_entities0 ); + qsort ( t->sorted_lentities[1].items, + num_entities, + sizeof(t->sorted_lentities[1].items[0]), + _qsort_cmp_entities1 ); + + return H5_SUCCESS; +} + +/*! + Sort (small) array of local vertex ids geometrically. + */ +h5_err_t +_h5t_sort_local_vertex_ids ( + h5_file * const f, + h5_id_t * const local_vids, /* IN/OUT: local vertex ids */ + const h5_size_t size /* size of array */ + ) { + struct h5t_fdata *t = &f->t; + + h5_size_t i; + for ( i = 1; i < size; ++i ) { + h5_id_t local_vid = local_vids[i]; + + h5_id_t j = i; + while ( (j >= 1 ) && _cmp_vertices ( + t->vertices[local_vid].P, + t->vertices[local_vids[j-1]].P + ) < 0 ) { + local_vids[j] = local_vids[j-1]; + --j; + } + local_vids[j] = local_vid; + } + return H5_SUCCESS; +} + +h5_err_t +_h5t_sort_global_vertex_ids ( + h5_file * const f, + h5_id_t * const global_vids, /* IN/OUT: global vertex ids */ + const h5_size_t size /* size of array */ + ) { + + h5_id_t local_vids[H5_MAX_VERTICES_PER_ENTITY]; + const h5_id_t *global_vid = global_vids; + h5_id_t *local_vid = local_vids; + + h5_id_t i; + for ( i = 0; i < size; i++, local_vid++, global_vid++ ) { + *local_vid = h5t_map_global_vertex_id2local ( f, *global_vid ); + if ( *local_vid < 0 ) + return *local_vid; + } + h5_err_t h5err = _h5t_sort_local_vertex_ids ( f, local_vids, size ); + if ( h5err < 0 ) return h5err; + + for ( i = 0; i < size; i++ ) { + global_vids[i] = h5t_map_local_vertex_id2global ( f, local_vids[i] ); + } + return H5_SUCCESS; +} + + +/*! + Binary search an entity given by its local vertex ids. + + \result index in t->map_entity_s2l[0].items + */ +static h5_id_t +_search_entity ( + h5_file *f, + h5_id_t * const local_vids /* local vertex ids */ + ) { + struct h5t_fdata *t = &f->t; + + _h5t_sort_local_vertex_ids ( f, local_vids, t->mesh_type ); + + register h5_id_t low = 0; + register h5_id_t high = t->sorted_lentities[0].num_items - 1; + register h5_id_t *entity1 = local_vids; + while (low <= high) { + register int mid = (low + high) / 2; + + h5_id_t local_eid = t->sorted_lentities[0].items[mid]; + h5_id_t *entity2 = t->lentities.tets[local_eid].vertex_ids; + int diff = _vcmp_entities ( f, entity1, entity2 ); + if ( diff < 0 ) + high = mid - 1; + else if ( diff > 0 ) + low = mid + 1; + else + return mid; // found + } + return -(low+1); // not found +} + +/*! + Get local id of entity given by local vertex id's +*/ +h5_id_t +h5t_get_local_entity_id ( + h5_file *f, + h5_id_t * const local_vids /* IN/OUT: local vertex id's */ + ) { + struct h5t_fdata *t = &f->t; + + if ( t->vertices == NULL ) { + h5_err_t h5err = _h5t_read_vertices ( f ); + if ( h5err < 0 ) return h5err; + } + + if ( t->entities.data == NULL ) { + h5_err_t h5err = _h5t_read_entities ( f ); + if ( h5err < 0 ) return h5err; + } + + if ( t->sorted_lentities[0].items == NULL ) { + _sort_entities ( f ); + } + h5_id_t local_eid = _search_entity ( f, local_vids ); + if ( local_eid < 0 ) { + return _h5t_handle_get_local_entity_id_err ( f, local_vids ); + } + + return t->sorted_lentities[0].items[local_eid]; +} + +/*! + Map a global vertex id to corresponding local vertex id. +*/ +h5_id_t +h5t_map_global_vertex_id2local ( + h5_file *f, + const h5_id_t global_id + ) { + struct h5t_fdata *t = &f->t; + + h5_id_t local_id = _h5_search_idmap ( &t->map_vertex_g2l, global_id ); + if ( local_id < 0 ) + return _h5t_handle_global_id_not_exist_err ("vertex", global_id ); + return local_id; +} + +h5_err_t +h5t_map_global_vertex_ids2local ( + h5_file *f, + const h5_id_t * const global_vids, + const h5_id_t size, + h5_id_t * const local_vids + ) { + h5_id_t i; + + for ( i = 0; i < size; i++ ) { + local_vids[i] = h5t_map_global_vertex_id2local ( + f, global_vids[i] ); + if ( local_vids[i] < 0 ) + return _h5t_handle_global_id_not_exist_err ( + "vertex", global_vids[i] ); + } + return H5_SUCCESS; +} + +/*! + Map a local entity id to corresponding global id. +*/ +h5_id_t +h5t_map_local_entity_id2global ( + h5_file *f, + const h5_id_t local_eid + ) { + struct h5t_fdata *t = &f->t; + + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: + if ( local_eid < 0 || local_eid > t->num_entities[t->num_levels-1] ) + return HANDLE_H5_OUT_OF_RANGE_ERR ( "tet", local_eid ); + return t->entities.tets[local_eid].id; + case TRIANGLE_MESH: + if ( local_eid < 0 || local_eid > t->num_entities[t->num_levels-1] ) + return HANDLE_H5_OUT_OF_RANGE_ERR ( "triangle", local_eid ); + return t->entities.tris[local_eid].id; + } + return -1; +} + +/*! + Get global id of entity given by global vertex id's +*/ +h5_id_t +h5t_get_global_entity_id ( + h5_file *f, + const h5_id_t * const global_vids /* global vertex id's */ + ) { + struct h5t_fdata *t = &f->t; + h5_id_t local_vids[H5_MAX_VERTICES_PER_ENTITY]; + h5_err_t h5err = h5t_map_global_vertex_ids2local ( + f, + global_vids, + t->mesh_type, + local_vids ); + if ( h5err < 0 ) + return _h5t_handle_get_global_entity_id_err ( f, global_vids ); + + h5_id_t local_eid = h5t_get_local_entity_id ( f, local_vids ); + if ( local_eid < 0 ) + return _h5t_handle_get_global_entity_id_err ( f, global_vids ); + + return h5t_map_local_entity_id2global ( f, local_eid ); +} + +static int +_search_ith_vertex_in_entity ( + h5_file * const f, + const int i, + h5_id_t local_vid + ) { + struct h5t_fdata *t = &f->t; + + register h5_id_t low = 0; + register h5_id_t high = t->sorted_lentities[i].num_items - 1; + register h5_float64_t *vertex1 = t->vertices[local_vid].P; + + while (low <= high) { + register int mid = (low + high) / 2; + + h5_id_t local_eid = t->sorted_lentities[i].items[mid]; + h5_id_t local_vid2 = t->lentities.tets[local_eid].vertex_ids[i]; + h5_float64_t *vertex2 = t->vertices[local_vid2].P; + int diff = _cmp_vertices ( vertex1, vertex2 ); + + if ( diff < 0 ) + high = mid - 1; + else if ( diff > 0 ) + low = mid + 1; + else + return mid; // found + } + return -(low+1); // not found +} + +static h5_id_t +_tetm_contain_triangle ( + h5_file *f, + const h5_id_t * const local_vids, + int i, + h5_id_t local_eid + ) { + struct h5t_fdata *t = &f->t; + + h5_id_t *local_vids_of_entity = t->lentities.tets[local_eid].vertex_ids; + if ( i == 0 && + local_vids[1] == local_vids_of_entity[1] && + local_vids[2] == local_vids_of_entity[2] + ) return 0; + else if ( i == 0 && + local_vids[1] == local_vids_of_entity[1] && + local_vids[2] == local_vids_of_entity[3] + ) return 1; + else if ( i == 0 && + local_vids[1] == local_vids_of_entity[2] && + local_vids[2] == local_vids_of_entity[3] + ) return 2; + else if ( i == 1 && + local_vids[1] == local_vids_of_entity[2] && + local_vids[2] == local_vids_of_entity[3] + ) return 3; + return -1; +} + +/*! + \return local triangle id + */ +h5_id_t +_tetm_search_triangle ( + h5_file *f, + h5_id_t * const local_vids + ) { + struct h5t_fdata *t = &f->t; + + _h5t_sort_local_vertex_ids ( f, local_vids, t->mesh_type ); + + /* + search for Vertex_id(tri,0) in the tuple of all 0th vertices of all tets + if there is one: + take the smallest tet with Vertex_id(tri,0) == Vertex_id(tet,0) + loop over all tets with Vertex_id(tri,0) == Vertex_id(tet,0) + until we find a tet the triangle is belonging to. + else + search for Vertex_id(tri,0) in the tuple of all 1st vertices of + all tets + if there is one: + take the smallest tet with + Vertex_id(tri,0) == Vertex_id(tet,0) + loop over all tets with + Vertex_id(tri,0) == Vertex_id(tet,0) + until we find a tet the triangle is belonging to. + */ + + int i = 0; + h5_id_t idx; + for ( i = 0; i < 2; i++ ) { + idx = _search_ith_vertex_in_entity ( f, i, local_vids[0] ); + if ( idx >= 0 ) break; + } + if ( idx < 0 ) { + return -1; + } + + /* get leftmost */ + while ( local_vids[0] == _get_local_vertex_id_of_entity ( + f, i, t->sorted_lentities[i].items[idx] ) ) + idx--; + + h5_id_t local_eid; + do { + /* check whether triangle is in entity given by local id */ + local_eid = t->sorted_lentities[i].items[idx]; + h5_id_t tidx = _tetm_contain_triangle ( f, local_vids, i, local_eid ); + if ( tidx >= 0 ) { + return _h5t_build_triangle_id ( tidx, local_eid ); + } + idx++; + } while ( local_vids[0] == _get_local_vertex_id_of_entity ( + f, i, local_eid ) ); + + return -1; +} + +/*! + + */ +h5_id_t +h5t_get_global_triangle_id ( + h5_file * const f, + h5_id_t * const global_vids + ) { + struct h5t_fdata *t = &f->t; + + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: { + h5_id_t local_vids[4]; + h5_err_t h5err = h5t_map_global_vertex_ids2local ( + f, global_vids, 4, local_vids ); + if ( h5err < 0 ) return h5err; + h5_id_t local_tid = h5t_get_local_triangle_id ( f, local_vids ); + if ( local_tid < 0 ) + return _h5t_handle_get_global_triangle_id_err( global_vids ); + return h5t_map_local_triangle_id2global ( f, local_tid ); + } + case TRIANGLE_MESH: + return h5t_get_global_entity_id ( f, global_vids ); + } + return -1; +} + +/*! + + */ +h5_id_t +h5t_get_local_triangle_id ( + h5_file * const f, + h5_id_t * const local_vids + ) { + struct h5t_fdata *t = &f->t; + + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: { + struct h5t_fdata *t = &f->t; + + if ( t->vertices == NULL ) { + h5_err_t h5err = _h5t_read_vertices ( f ); + if ( h5err < 0 ) return h5err; + } + + if ( t->entities.data == NULL ) { + h5_err_t h5err = _h5t_read_entities ( f ); + if ( h5err < 0 ) return h5err; + } + + if ( t->sorted_lentities[0].items == NULL ) { + _sort_entities ( f ); + } + h5_id_t local_tid = _tetm_search_triangle ( f, local_vids ); + if ( local_tid == -1 ) { + return _h5t_handle_get_local_triangle_id_err ( + local_vids ); + } + return local_tid; + } + case TRIANGLE_MESH: + return h5t_get_local_entity_id ( f, local_vids ); + } + return -1; +} + +h5_id_t +h5t_map_global_entity_id2local ( + h5_file * const f, + const h5_id_t global_eid + ) { + struct h5t_fdata *t = &f->t; + h5_id_t local_eid = _h5_search_idmap ( &t->map_entity_g2l, global_eid ); + if ( local_eid < 0 ) + return _h5t_handle_global_id_not_exist_err ( "entity", global_eid ); + return local_eid; +} + +/*! + Map global triangle id to local id. + + \return local id of triangle + */ +h5_id_t +h5t_map_global_triangle_id2local ( + h5_file * const f, + const h5_id_t global_tri_id + ) { + struct h5t_fdata *t = &f->t; + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: { + h5_id_t global_tet_id = global_tri_id & H5_TET_MASK; + h5_id_t local_tet_id = h5t_map_global_entity_id2local ( + f, global_tet_id ); + if ( local_tet_id < 0 ) + return _h5t_handle_global_id_not_exist_err ( + "triangle", global_tri_id ); + return local_tet_id | (global_tri_id & ~H5_TET_MASK); + } + case TRIANGLE_MESH: + return h5t_map_global_entity_id2local ( f, global_tri_id ); + } + return -1; +} + +/*! + Map local triangle id to global id. + + \return global id of triangle + */ +h5_id_t +h5t_map_local_triangle_id2global ( + h5_file * const f, + const h5_id_t local_tri_id + ) { + struct h5t_fdata *t = &f->t; + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: { + h5_id_t local_tet_id = local_tri_id & H5_TET_MASK; + h5_id_t global_tet_id = h5t_map_local_entity_id2global ( + f, local_tet_id ); + if ( global_tet_id < 0 ) + return HANDLE_H5_OUT_OF_RANGE_ERR( "triangle", local_tri_id ); + return global_tet_id | (local_tri_id & ~H5_TET_MASK); + } + case TRIANGLE_MESH: + return h5t_map_local_entity_id2global ( f, local_tri_id ); + } + return -1; +} diff --git a/src/h5/t_map.h b/src/h5/t_map.h new file mode 100644 index 0000000..405d105 --- /dev/null +++ b/src/h5/t_map.h @@ -0,0 +1,72 @@ +#ifndef __T_MAP_H +#define __T_MAP_H + +h5_id_t +h5t_map_global_vertex_id2local ( + h5_file *f, + h5_id_t global_vid + ); + +h5_err_t +h5t_map_global_vertex_ids2local ( + h5_file *f, + const h5_id_t * const global_vids, + const h5_id_t size, + h5_id_t * const local_vids + ); + +h5_id_t +h5t_map_local_vertex_id2global ( + h5_file *f, + h5_id_t local_vid + ); + +h5_id_t +h5t_map_local_entity_id2global ( + h5_file *f, + h5_id_t local_eid + ); + +h5_id_t +h5t_map_global_entity_id2local ( + h5_file * const f, + const h5_id_t global_eid + ); + +h5_id_t +h5t_map_local_triangle_id2global ( + h5_file * const f, + const h5_id_t local_tid + ); + +h5_id_t +h5t_map_global_triangle_id2local ( + h5_file * const f, + const h5_id_t global_tid + ); + +h5_id_t +h5t_get_local_entity_id ( + h5_file *f, + h5_id_t * const local_vids + ); + +h5_id_t +h5t_get_global_entity_id ( + h5_file *f, + const h5_id_t * const global_vids + ); + +h5_id_t +h5t_get_local_triangle_id ( + h5_file * const f, + h5_id_t * const local_vids + ); + +h5_id_t +h5t_get_global_triangle_id ( + h5_file * const f, + h5_id_t * const global_vids + ); + +#endif diff --git a/src/h5/t_map_private.h b/src/h5/t_map_private.h new file mode 100644 index 0000000..6f4ea2f --- /dev/null +++ b/src/h5/t_map_private.h @@ -0,0 +1,18 @@ +#ifndef __T_MAP_PRIVATE_H +#define __T_MAP_PRIVATE_H + +h5_err_t +_h5t_sort_global_vertex_ids ( + h5_file * const f, + h5_id_t * const global_vids, + const h5_size_t size + ); + +h5_err_t +_h5t_sort_local_vertex_ids ( + h5_file * const f, + h5_id_t * const local_vids, + const h5_size_t size + ); + +#endif diff --git a/src/h5/t_openclose.c b/src/h5/t_openclose.c index 2ea66be..7797d8d 100644 --- a/src/h5/t_openclose.c +++ b/src/h5/t_openclose.c @@ -11,10 +11,6 @@ #include "h5/h5_private.h" #include "H5Part.h" -extern h5part_error_handler _err_handler; -extern h5part_int64_t _h5part_errno; -extern unsigned _debug; - /* create several HDF5 types */ diff --git a/src/h5/t_readwrite.c b/src/h5/t_readwrite.c index e116f0c..16d5896 100644 --- a/src/h5/t_readwrite.c +++ b/src/h5/t_readwrite.c @@ -8,6 +8,7 @@ #include "h5_core.h" #include "h5_private.h" + #include "H5Part.h" #include "H5Block.h" @@ -513,6 +514,24 @@ h5t_add_level ( return t->cur_level; } +static h5_err_t +_alloc_num_vertices ( + h5_file * f, + const h5_size_t num_vertices + ) { + struct h5t_fdata *t = &f->t; + + ssize_t num_bytes = num_vertices*sizeof ( t->vertices[0] ); + h5_debug ( "Allocating %ld bytes.", num_bytes ); + t->vertices = realloc ( t->vertices, num_bytes ); + if ( t->vertices == NULL ) { + return HANDLE_H5_NOMEM_ERR; + } + + return _h5_alloc_idmap (&t->map_vertex_g2l, num_vertices ); +} + + h5_err_t h5t_add_num_vertices ( h5_file * f, @@ -523,20 +542,11 @@ h5t_add_num_vertices ( if ( t->cur_level < 0 ) { return HANDLE_H5_UNDEF_LEVEL_ERR; } - ssize_t num_elems = (t->cur_level > 0 ? + ssize_t num_vertices = (t->cur_level > 0 ? t->num_vertices[t->cur_level-1] + num : num); - t->num_vertices[t->cur_level] = num_elems; - ssize_t num_bytes = num_elems*sizeof ( t->vertices[0] ); - h5_debug ( "Allocating %ld bytes.", num_bytes ); - t->vertices = realloc ( t->vertices, num_bytes ); - if ( t->vertices == NULL ) { - return HANDLE_H5_NOMEM_ERR; - } + t->num_vertices[t->cur_level] = num_vertices; - h5_err_t h5err = _h5_alloc_idmap (&t->map_vertex_g2l, num ); - if ( h5err < 0 ) return h5err; - - return H5_SUCCESS; + return _alloc_num_vertices ( f, num_vertices ); } /* @@ -616,9 +626,6 @@ _read_num_vertices ( h5_err_t h5err; struct h5t_fdata *t = &f->t; - if ( t->cur_level < 0 ) - return HANDLE_H5_UNDEF_LEVEL_ERR; - if ( t->mesh_gid < 0 ) { h5err = _open_mesh_group ( f ); if ( h5err < 0 ) return h5err; @@ -641,16 +648,13 @@ _read_num_vertices ( return H5_SUCCESS; } -static h5_err_t -_read_vertices ( +h5_err_t +_h5t_read_vertices ( h5_file * f ) { h5_err_t h5err; struct h5t_fdata *t = &f->t; - if ( t->cur_level < 0 ) - return HANDLE_H5_UNDEF_LEVEL_ERR; - if ( t->mesh_gid < 0 ) { h5err = _open_mesh_group ( f ); if ( h5err < 0 ) return h5err; @@ -661,12 +665,9 @@ _read_vertices ( if ( h5err < 0 ) return h5err; } - ssize_t num_elems = t->num_vertices[t->num_levels-1]; - ssize_t num_bytes = num_elems*sizeof ( t->vertices[0] ); - h5_debug ( "Allocating %ld bytes.", num_bytes ); - t->vertices = realloc ( t->vertices, num_bytes ); - if ( t->vertices == NULL ) - return HANDLE_H5_NOMEM_ERR; + h5err = _alloc_num_vertices ( f, t->num_vertices[t->num_levels-1] ); + if ( h5err < 0 ) return h5err; + h5err = _read_dataset ( f, t->mesh_gid, @@ -677,11 +678,19 @@ _read_vertices ( t->vertices ); if ( h5err < 0 ) return h5err; + h5_id_t local_vid = 0; + for ( ; local_vid < t->num_vertices[t->num_levels-1]; local_vid++ ) { + t->map_vertex_g2l.items[local_vid].global_id = + t->vertices[local_vid].id; + t->map_vertex_g2l.items[local_vid].local_id = local_vid; + t->map_vertex_g2l.num_items++; + } + _h5_sort_idmap ( &t->map_vertex_g2l ); return H5_SUCCESS; } h5_size_t -h5t_get_num_vertices ( +h5t_get_num_vertices_on_level ( h5_file * f ) { struct h5t_fdata *t = &f->t; @@ -762,7 +771,7 @@ h5t_traverse_vertices ( struct h5t_fdata *t = &f->t; if ( t->vertices == NULL ) { - h5_err_t h5err = _read_vertices ( f ); + h5_err_t h5err = _h5t_read_vertices ( f ); if ( h5err < 0 ) return h5err; } if ( t->last_retrieved_vertex_id+1 >= t->num_vertices[t->cur_level] ) { @@ -776,23 +785,15 @@ h5t_traverse_vertices ( return t->last_retrieved_vertex_id; } -h5_err_t -h5t_add_num_entities ( +static h5_err_t +_alloc_num_entities ( h5_file * f, - const h5_size_t num + size_t cur_num_entities, + size_t new_num_entities ) { struct h5t_fdata *t = &f->t; - ssize_t sizeof_entity = 0; - ssize_t sizeof_lentity = 0; - - ssize_t cur_num_entities = t->cur_level > 0 ? - t->num_entities[t->cur_level-1] : 0; - ssize_t num_entities = t->cur_level > 0 ? - num + t->num_entities[t->cur_level-1] : num; - t->num_entities[t->cur_level] = num_entities; - - t->num_entities_on_level[t->cur_level] = t->cur_level > 0 ? - num + t->num_entities_on_level[t->cur_level-1] : num; + size_t sizeof_entity = 0; + size_t sizeof_lentity = 0; switch ( t->mesh_type ) { case TETRAHEDRAL_MESH: @@ -808,160 +809,44 @@ h5t_add_num_entities ( } t->entities.data = realloc ( - t->entities.data, num_entities*sizeof_entity ); + t->entities.data, new_num_entities * sizeof_entity ); if ( t->entities.data == NULL ) { return H5_ERR_NOMEM; } + t->lentities.data = realloc ( - t->lentities.data, num_entities*sizeof_lentity ); + t->lentities.data, new_num_entities*sizeof_lentity ); if ( t->lentities.data == NULL ) { return H5_ERR_NOMEM; } memset ( t->lentities.data+cur_num_entities*sizeof_lentity, -1, - (num_entities-cur_num_entities) * sizeof_lentity ); + (new_num_entities-cur_num_entities) * sizeof_lentity ); - h5_err_t h5err = _h5_alloc_idmap (&t->map_entity_g2l, num_entities ); - if ( h5err < 0 ) return h5err; - - return H5_SUCCESS; + return _h5_alloc_idmap (&t->map_entity_g2l, new_num_entities ); } -static h5_id_t -_get_local_vertex_id_of_entity ( - h5_file * f, - int vertex, - h5_id_t local_eid - ) { - struct h5t_fdata *t = &f->t; - - h5_id_t local_vid = -1; - h5_id_t global_vid = -1; - - switch ( t->mesh_type ) { - case TETRAHEDRAL_MESH: { - local_vid = t->lentities.tets[local_eid].vertex_ids[vertex]; - if ( local_vid == -1 ) { - global_vid = - t->entities.tets[local_eid].vertex_ids[vertex]; - local_vid = - _h5_search_idmap ( &t->map_vertex_g2l, global_vid ); - t->lentities.tets[local_eid].vertex_ids[vertex] = local_vid; - } - break; - } - case TRIANGLE_MESH: { - local_vid = t->lentities.tris[local_eid].vertex_ids[vertex]; - if ( local_vid == -1 ) { - global_vid = - t->entities.tris[local_eid].vertex_ids[vertex]; - local_vid = - _h5_search_idmap ( &t->map_vertex_g2l, global_vid ); - t->lentities.tris[local_eid].vertex_ids[vertex] = local_vid; - } - break; - } - } - return local_vid; -} - -static h5_float64_t* -_get_vertex_of_entity ( - h5_file * f, - int vertex, - h5_id_t local_entity_id - ) { - struct h5t_fdata *t = &f->t; - - h5_id_t local_vid = _get_local_vertex_id_of_entity ( f, vertex, local_entity_id ); - if ( local_vid == -1 ) - return NULL; - return t->vertices[local_vid].P; -} - - -static int -_cmp_vertices ( - h5_float64_t *P1, - h5_float64_t *P2 - ) { - int i; - for ( i = 0; i < 3; i++ ) { - if ( P1[i] < P2[i] ) return -1; - else if (P1[i] > P2[i] )return 1; - } - return 0; -} - -static int -_cmp_entities ( - h5_file * f, - h5_id_t lid1, - h5_id_t lid2 - ) { - - int i; - for ( i = 0; i < 4; i++ ) { - int r = _cmp_vertices ( - _get_vertex_of_entity ( f, i, lid1 ), - _get_vertex_of_entity ( f, i, lid2 ) ); - if ( r < 0 ) return -1; - else if ( r > 0 ) return 1; - } - return 0; -} - - -/*! - Insertion sort of vertices - for small arrays only! -*/ h5_err_t -_h5t_isort_vertices ( - h5_file *f, - h5_id_t *global_ids, /* global ids */ - h5_size_t size /* size of array */ +h5t_add_num_entities ( + h5_file * f, + const h5_size_t num ) { struct h5t_fdata *t = &f->t; - h5_id_t *global_id = global_ids; + size_t cur_num_entities = t->cur_level > 0 ? + t->num_entities[t->cur_level-1] : 0; + size_t new_num_entities = t->cur_level > 0 ? + num + t->num_entities[t->cur_level-1] : num; + t->num_entities[t->cur_level] = new_num_entities; - h5_id_t *local_ids = malloc ( size * sizeof(global_ids[0]) ); - h5_id_t *local_id = local_ids; - h5_err_t h5err = H5_SUCCESS; - h5_size_t i; + t->num_entities_on_level[t->cur_level] = t->cur_level > 0 ? + num + t->num_entities_on_level[t->cur_level-1] : num; - if ( local_ids == NULL ) - return HANDLE_H5_NOMEM_ERR; - - for ( i = 0; i < size; i++, local_id++, global_id++ ) { - *local_id = _h5_search_idmap ( &t->map_vertex_g2l, *global_id ); - if ( *local_id < 0 ) { - h5err = *local_id; - goto cleanup; - } - } - - for ( i = 1; i < size; ++i ) { - h5_id_t gid = global_ids[i]; - h5_id_t lid = local_ids[i]; - h5_id_t j = i; - while ( _cmp_vertices ( - t->vertices[lid].P, - t->vertices[local_ids[j-1]].P - ) < 0 ) { - global_ids[j] = global_ids[j-1]; - local_ids[j] = local_ids[j-1]; - --j; - } - global_ids[j] = gid; - local_ids[j] = lid; - } -cleanup: - free ( local_ids ); - return h5err; + return _alloc_num_entities ( f, cur_num_entities, new_num_entities ); } + h5_id_t h5t_store_tet ( h5_file * f, @@ -977,7 +862,8 @@ h5t_store_tet ( more than allocated */ if ( t->last_stored_entity_id+1 >= t->num_entities[t->cur_level] ) - return HANDLE_H5_OVERFLOW_ERR( "tet", t->num_entities[t->cur_level] ); + return HANDLE_H5_OVERFLOW_ERR( + "tet", t->num_entities[t->cur_level] ); /* missing call to add the first level @@ -988,13 +874,16 @@ h5t_store_tet ( /* check parent id */ - if ( (t->cur_level == 0) && (parent_id != -1) ) { + if ( (t->cur_level == 0) && + (parent_id != -1) ) { return HANDLE_H5_PARENT_ID_ERR ( "tet", global_id, parent_id ); } - if ( (t->cur_level > 0) && (parent_id < 0) ) { + if ( (t->cur_level > 0) && + (parent_id < 0) ) { return HANDLE_H5_PARENT_ID_ERR ( "tet", global_id, parent_id ); } - if ( (t->cur_level > 0) && (parent_id >= t->num_entities[t->cur_level-1]) ) { + if ( (t->cur_level > 0) && + (parent_id >= t->num_entities[t->cur_level-1]) ) { return HANDLE_H5_PARENT_ID_ERR ( "tet", global_id, parent_id ); } /* @@ -1019,14 +908,15 @@ h5t_store_tet ( memcpy ( &tet->vertex_ids, vertex_ids, sizeof ( tet->vertex_ids ) ); - _h5t_isort_vertices ( f, tet->vertex_ids, 4 ); + _h5t_sort_global_vertex_ids ( f, tet->vertex_ids, 4 ); _h5_insert_idmap ( &t->map_entity_g2l, global_id, local_id ); if ( parent_id >= 0 ) { h5_id_t local_parent_id = _h5_search_idmap ( &t->map_entity_g2l, parent_id ); if ( t->entities.tets[local_parent_id].refined_on_level < 0 ) { - t->entities.tets[local_parent_id].refined_on_level = t->cur_level; + t->entities.tets[local_parent_id].refined_on_level = + t->cur_level; t->num_entities_on_level[t->cur_level]--; } } @@ -1049,7 +939,8 @@ h5t_store_triangle ( more than allocated */ if ( t->last_stored_entity_id+1 >= t->num_entities[t->cur_level] ) - return HANDLE_H5_OVERFLOW_ERR( "triangle", t->num_entities[t->cur_level] ); + return HANDLE_H5_OVERFLOW_ERR( + "triangle", t->num_entities[t->cur_level] ); /* missing call to add the first level @@ -1061,13 +952,16 @@ h5t_store_triangle ( check parent id */ if ( (t->cur_level == 0) && (parent_id != -1) ) { - return HANDLE_H5_PARENT_ID_ERR ( "triangle", global_id, parent_id ); + return HANDLE_H5_PARENT_ID_ERR ( + "triangle", global_id, parent_id ); } if ( (t->cur_level > 0) && (parent_id < 0) ) { - return HANDLE_H5_PARENT_ID_ERR ( "triangle", global_id, parent_id ); + return HANDLE_H5_PARENT_ID_ERR ( + "triangle", global_id, parent_id ); } - if ( (t->cur_level > 0) && (parent_id >= t->num_entities[t->cur_level-1]) ) { - return HANDLE_H5_PARENT_ID_ERR ( "triangle", global_id, parent_id ); + if ( (t->cur_level>0) && (parent_id >= t->num_entities[t->cur_level-1]) ) { + return HANDLE_H5_PARENT_ID_ERR ( + "triangle", global_id, parent_id ); } /* check id @@ -1095,7 +989,8 @@ h5t_store_triangle ( h5_id_t local_parent_id = _h5_search_idmap ( &t->map_entity_g2l, parent_id ); if ( t->entities.tris[local_parent_id].refined_on_level < 0 ) { - t->entities.tris[local_parent_id].refined_on_level = t->cur_level; + t->entities.tris[local_parent_id].refined_on_level = + t->cur_level; t->num_entities_on_level[t->cur_level]--; } } @@ -1111,9 +1006,6 @@ _read_num_entities ( h5_err_t h5err; struct h5t_fdata *t = &f->t; - if ( t->cur_level < 0 ) - return HANDLE_H5_UNDEF_LEVEL_ERR; - if ( t->mesh_gid < 0 ) { h5err = _open_mesh_group ( f ); if ( h5err < 0 ) return h5err; @@ -1185,17 +1077,13 @@ _open_file_space_entities ( return H5S_ALL; } -static h5_err_t -_read_entities ( +h5_err_t +_h5t_read_entities ( h5_file * f ) { h5_err_t h5err; - ssize_t num_bytes; struct h5t_fdata *t = &f->t; - if ( t->cur_level < 0 ) - return HANDLE_H5_UNDEF_LEVEL_ERR; - if ( t->mesh_gid < 0 ) { h5err = _open_mesh_group ( f ); if ( h5err < 0 ) return h5err; @@ -1206,24 +1094,8 @@ _read_entities ( if ( h5err < 0 ) return h5err; } - ssize_t num_elems = t->num_entities[t->num_levels-1]; + h5err = _alloc_num_entities ( f, 0, t->num_entities[t->num_levels-1] ); - switch ( t->mesh_type ) { - case TETRAHEDRAL_MESH: { - num_bytes = num_elems*sizeof ( t->entities.tets[0] ); - break; - } - case TRIANGLE_MESH: { - num_bytes = num_elems*sizeof ( t->entities.tris[0] ); - break; - } - default: - return -1; - } - h5_debug ( "Allocating %ld bytes.", num_bytes ); - t->entities.data = realloc ( t->entities.data, num_bytes ); - if ( t->entities.data == NULL ) - return HANDLE_H5_NOMEM_ERR; h5err = _read_dataset ( f, t->mesh_gid, @@ -1234,11 +1106,43 @@ _read_entities ( t->entities.data ); if ( h5err < 0 ) return h5err; + /* + setup structure with local vertex ids + */ + h5_id_t local_eid = 0; + h5_id_t num_entities = t->num_entities[t->num_levels-1]; + switch ( t->mesh_type ) { + case TETRAHEDRAL_MESH: + for ( local_eid = 0; local_eid < num_entities; local_eid++ ) { + h5err = h5t_map_global_vertex_ids2local ( + f, + t->entities.tets[local_eid].vertex_ids, + t->mesh_type, + t->lentities.tets[local_eid].vertex_ids + ); + if ( h5err < 0 ) return h5err; + } + break; + case TRIANGLE_MESH: + for ( local_eid = 0; local_eid < num_entities; local_eid++ ) { + h5err = h5t_map_global_vertex_ids2local ( + f, + t->entities.tris[local_eid].vertex_ids, + t->mesh_type, + t->lentities.tris[local_eid].vertex_ids + ); + if ( h5err < 0 ) return h5err; + } + break; + default: + return -1; + } + return H5_SUCCESS; } h5_size_t -h5t_get_num_entities ( +h5t_get_num_entities_on_level ( h5_file * f ) { struct h5t_fdata *t = &f->t; @@ -1256,6 +1160,22 @@ h5t_get_num_entities ( return t->num_entities_on_level[t->cur_level]; } +h5_size_t +h5t_get_num_entities ( + h5_file * f + ) { + struct h5t_fdata *t = &f->t; + + if ( t->cur_mesh < 0 ) { + return HANDLE_H5_UNDEF_MESH_ERR; + } + if ( t->num_entities == NULL ) { + h5_err_t h5err = _read_num_entities ( f ); + if ( h5err < 0 ) return h5err; + } + return t->num_entities[t->num_levels-1]; +} + h5_err_t h5t_start_traverse_tets ( h5_file * f @@ -1286,7 +1206,7 @@ _traverse_tets ( struct h5t_fdata *t = &f->t; if ( t->entities.data == NULL ) { - h5_err_t h5err = _read_entities ( f ); + h5_err_t h5err = _h5t_read_entities ( f ); if ( h5err < 0 ) return h5err; } if ( t->last_retrieved_entity_id+1 >= t->num_entities[t->cur_level] ) { @@ -1365,7 +1285,7 @@ _traverse_triangles ( struct h5t_fdata *t = &f->t; if ( t->entities.data == NULL ) { - h5_err_t h5err = _read_entities ( f ); + h5_err_t h5err = _h5t_read_entities ( f ); if ( h5err < 0 ) return h5err; } if ( t->last_retrieved_entity_id+1 >= t->num_entities[t->cur_level] ) { diff --git a/src/h5/t_readwrite.h b/src/h5/t_readwrite.h index 5199052..2eba30b 100644 --- a/src/h5/t_readwrite.h +++ b/src/h5/t_readwrite.h @@ -55,6 +55,11 @@ h5t_add_level ( h5_file * f ); +h5_err_t +_h5t_read_vertices ( + h5_file * f + ); + h5_size_t h5t_add_num_vertices ( h5_file * f, @@ -69,7 +74,7 @@ h5t_store_vertex ( ); h5_size_t -h5t_get_num_vertices ( +h5t_get_num_vertices_on_level ( h5_file * f ); @@ -96,6 +101,11 @@ h5t_get_num_entities ( h5_file * f ); +h5_size_t +h5t_get_num_entities_on_level ( + h5_file * f + ); + h5_id_t h5t_store_tet ( h5_file * f, @@ -112,6 +122,11 @@ h5t_store_triangle ( const h5_id_t vertex_ids[3] ); +h5_err_t +_h5t_read_entities ( + h5_file * f + ); + h5_err_t h5t_start_traverse_tets ( h5_file * f diff --git a/src/h5/u_readwrite.c b/src/h5/u_readwrite.c index 4e393e4..f8ca911 100644 --- a/src/h5/u_readwrite.c +++ b/src/h5/u_readwrite.c @@ -10,10 +10,6 @@ #include "h5_private.h" #include "H5Part.h" -extern h5part_error_handler _err_handler; -extern h5part_int64_t _h5part_errno; -extern unsigned _debug; - static hid_t _get_diskshape_for_reading ( h5_file *f, @@ -180,7 +176,7 @@ H5U_read_elems ( herr = H5Dclose ( dataset_id ); if ( herr < 0 ) return HANDLE_H5D_CLOSE_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -201,7 +197,7 @@ H5U_set_num_elements ( we don't know if things have changed globally */ if ( f->nparticles == nparticles ) { - return H5PART_SUCCESS; + return H5_SUCCESS; } #endif if ( f->diskshape != H5S_ALL ) { @@ -306,7 +302,7 @@ H5U_set_num_elements ( } #endif - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -357,7 +353,7 @@ H5U_reset_view ( if ( herr < 0 ) return HANDLE_H5S_CLOSE_ERR; f->memshape=H5S_ALL; } - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -378,7 +374,7 @@ H5U_set_view ( herr = H5U_reset_view ( f ); if ( herr < 0 ) return herr; - if ( start == -1 && end == -1 ) return H5PART_SUCCESS; + if ( start == -1 && end == -1 ) return H5_SUCCESS; /* View has been reset so H5PartGetNumParticles will tell @@ -434,7 +430,7 @@ H5U_set_view ( NULL ); if ( herr < 0 ) return HANDLE_H5S_SELECT_HYPERSLAB_ERR; - return H5PART_SUCCESS; + return H5_SUCCESS; } h5part_int64_t @@ -507,7 +503,7 @@ H5U_set_canonical_view ( #endif - return H5PART_SUCCESS; + return H5_SUCCESS; } @@ -544,5 +540,5 @@ H5U_get_dataset_info ( if ( *type < 0 ) return *type; } - return H5PART_SUCCESS; + return H5_SUCCESS; }