From 9681dbf3b13384198b40c6dc072cb019b1a9186b Mon Sep 17 00:00:00 2001 From: Achim Gsell Date: Thu, 17 Jul 2008 14:01:20 +0000 Subject: [PATCH] mappings and bookkeeping functions we need for boundary triangles --- .gitattributes | 2 + src/h5/Makefile.am | 2 + src/h5/errorhandling_private.h | 7 + src/h5/h5_core.h | 1 + src/h5/h5_types.h | 61 ++++++-- src/h5/maps.c | 92 +++++++++++ src/h5/maps.h | 28 ++++ src/h5/t_readwrite.c | 273 ++++++++++++++++++++++++++------- 8 files changed, 394 insertions(+), 72 deletions(-) create mode 100644 src/h5/maps.c create mode 100644 src/h5/maps.h diff --git a/.gitattributes b/.gitattributes index 3eed191..dce6b3e 100644 --- a/.gitattributes +++ b/.gitattributes @@ -383,6 +383,8 @@ src/h5/general.c -text src/h5/h5_core.h -text src/h5/h5_private.h -text src/h5/h5_types.h -text +src/h5/maps.c -text +src/h5/maps.h -text src/h5/readwrite.c -text src/h5/readwrite.h -text src/h5/t_openclose.c -text diff --git a/src/h5/Makefile.am b/src/h5/Makefile.am index 14e33c3..8f5dfdd 100644 --- a/src/h5/Makefile.am +++ b/src/h5/Makefile.am @@ -47,6 +47,7 @@ include_HEADERS = \ h5_types.h \ attribs.h \ errorhandling.h \ + maps.h \ readwrite.h \ t_openclose.h \ t_readwrite.h \ @@ -60,6 +61,7 @@ libH5_a_SOURCES = \ general.c \ attribs.c \ errorhandling.c \ + maps.c \ t_openclose.c \ readwrite.c \ u_readwrite.c \ diff --git a/src/h5/errorhandling_private.h b/src/h5/errorhandling_private.h index 4e9c855..51a5237 100644 --- a/src/h5/errorhandling_private.h +++ b/src/h5/errorhandling_private.h @@ -110,6 +110,13 @@ "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()) ( \ h5_get_funcname(), \ diff --git a/src/h5/h5_core.h b/src/h5/h5_core.h index 2f07e2c..3ae880d 100644 --- a/src/h5/h5_core.h +++ b/src/h5/h5_core.h @@ -53,6 +53,7 @@ h5_get_step ( #include "attribs.h" #include "errorhandling.h" +#include "maps.h" #include "readwrite.h" #include "t_openclose.h" #include "t_readwrite.h" diff --git a/src/h5/h5_types.h b/src/h5/h5_types.h index 758233a..4e7f93c 100644 --- a/src/h5/h5_types.h +++ b/src/h5/h5_types.h @@ -1,4 +1,3 @@ - #ifndef __H5_TYPES_H #define __H5_TYPES_H @@ -50,9 +49,9 @@ struct h5_complex { }; typedef struct h5_complex h5_complex; -enum h5_mesh_types { - TETRAHEDRAL_MESH, - TRIANGLE_MESH +enum h5_mesh_types { /* enum with number of vertices(!) */ + TRIANGLE_MESH = 3, + TETRAHEDRAL_MESH = 4 }; struct h5_vertex { /* 32Byte */ @@ -76,6 +75,7 @@ struct h5_triangle { /*24Bytes*/ h5_id_t refined_on_level; }; + struct h5_tetrahedron { /* 24Bytes */ h5_id_t id; h5_id_t parent_id; @@ -84,6 +84,14 @@ struct h5_tetrahedron { /* 24Bytes */ h5_id_t vertex_ids[4]; }; +struct h5_ltriangle { + h5_id_t vertex_ids[3]; /* local(!) vertex ids */ +}; + +struct h5_ltetrahedron { + h5_id_t vertex_ids[4]; /* local(!) vertex ids */ +}; + typedef struct h5_vertex h5_vertex; typedef struct h5_edge h5_edge; typedef struct h5_triangle h5_triangle; @@ -151,9 +159,30 @@ struct h5t_fdata_level { }; union entities { - h5_tetrahedron *tets; - h5_triangle *tris; - void *data; + struct h5_tetrahedron *tets; + struct h5_triangle *tris; + void *data; +}; + +union lentities { + struct h5_ltetrahedron *tets; + struct h5_ltriangle *tris; + void *data; +}; + +struct smap { + h5_size_t size; /* allocated space in number of items */ + h5_size_t num_items; /* stored items */ + h5_id_t *items; +}; + +struct idmap { + h5_size_t size; /* allocated space in number of items */ + h5_size_t num_items; /* stored items */ + struct { + h5_id_t global_id; + h5_id_t local_id; + } *items; }; struct h5t_fdata { @@ -167,21 +196,22 @@ struct h5t_fdata { h5_id_t new_level; /* idx of the first new level or -1 */ h5_size_t num_levels; + h5_vertex *vertices; h5_size_t *num_vertices; - h5_size_t *num_entities; - h5_size_t *num_entities_on_level; - h5_size_t *map_entity_g2l;/* map global id to local id */ - + struct idmap map_vertex_g2l;/* map global id to local id */ h5_id_t last_retrieved_vertex_id; h5_id_t last_stored_vertex_id; - h5_vertex *vertices; + + union entities entities; + union lentities lentities; + hid_t entity_tid; /* type of mesh: tetrahedral, triangle ... */ + h5_size_t *num_entities; + h5_size_t *num_entities_on_level; + struct idmap map_entity_g2l;/* map global id to local id */ h5_id_t last_retrieved_entity_id; h5_id_t last_stored_entity_id; - hid_t entity_tid; - union entities entities; - /* HDF5 objects */ hid_t topo_gid; /* grp id of mesh in current level */ @@ -198,7 +228,6 @@ struct h5t_fdata { hid_t tet_tid; }; - /** \struct h5_file diff --git a/src/h5/maps.c b/src/h5/maps.c new file mode 100644 index 0000000..d0e9cc9 --- /dev/null +++ b/src/h5/maps.c @@ -0,0 +1,92 @@ +#include + +#include +#include "h5_types.h" +#include "h5_core.h" +#include "h5_private.h" + + +h5_err_t +_h5_alloc_idmap ( + struct idmap *map, + h5_size_t size + ) { + map->items = realloc ( map->items, size * sizeof ( map->items[0] ) ); + if ( map->items == NULL ) { + return HANDLE_H5_NOMEM_ERR; + } + map->size = size; + return H5_SUCCESS; +} + + +h5_err_t +_h5_insert_idmap ( + struct idmap *map, + h5_id_t global_id, + h5_id_t local_id + ) { + + if ( map->num_items == map->size ) + return HANDLE_H5_OVERFLOW_ERR( "g2lmap", map->size ); + + h5_id_t i = _h5_search_idmap ( map, global_id ); + if ( i >= 0 ) /* global id already in use ? */ + return -1; + + i = -(i+1); + + memmove ( + &map->items[i+1], + &map->items[i], + map->num_items - i ); + map->items[i].global_id = global_id; + map->items[i].local_id = local_id; + map->num_items++; + + return H5_SUCCESS; +} + +/*! + + \ingroup h5_core + + binary search in simple map + + \return index in array if found, othwise \c -(result+1) is the index + where \c value must be inserted. + + */ +h5_id_t +_h5_search_idmap ( + struct idmap *map, + h5_id_t value + ) { + + register int low = 0; + register int high = map->num_items - 1; + while (low <= high) { + register int mid = (low + high) / 2; + register h5_id_t diff = map->items[mid].global_id - value; + if ( diff > 0 ) + high = mid - 1; + else if ( diff < 0 ) + low = mid + 1; + else + return mid; // found + } + return -(low+1); // not found +} + +h5_id_t +h5t_map_vertex_id_global2local ( + h5_file *f, + 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 HANDLE_H5T_GID_NOT_EXIST_ERR ( "vertex", global_id ); + return local_id; +} diff --git a/src/h5/maps.h b/src/h5/maps.h new file mode 100644 index 0000000..344102f --- /dev/null +++ b/src/h5/maps.h @@ -0,0 +1,28 @@ +#ifndef __MAPS_H +#define __MAPS_H + +h5_err_t +_h5_alloc_idmap ( + struct idmap *map, + h5_size_t size + ); + +h5_err_t +_h5_insert_idmap ( + struct idmap *map, + h5_id_t global_id, + h5_id_t local_id + ); + +h5_id_t +_h5_search_idmap ( + struct idmap *map, + h5_id_t value + ); + +h5_id_t +h5t_map_vertex_id_global2local ( + h5_file *f, + h5_id_t global_id + ); +#endif diff --git a/src/h5/t_readwrite.c b/src/h5/t_readwrite.c index a8049ed..e116f0c 100644 --- a/src/h5/t_readwrite.c +++ b/src/h5/t_readwrite.c @@ -136,10 +136,6 @@ _write_obj ( return H5_SUCCESS; } -/* - store vertices -> Topo.VERTICES_COORD3D - levels -> TOPO.VERTICES_LEVELS3D - */ static h5_err_t _write_vertices ( h5_file * f @@ -280,10 +276,15 @@ _release_memory ( } t->num_entities_on_level = NULL; - if ( t->map_entity_g2l ) { - free ( t->map_entity_g2l ); + if ( t->map_vertex_g2l.items ) { + free ( t->map_vertex_g2l.items ); } - t->map_entity_g2l = NULL; + t->map_vertex_g2l.items = NULL; + + if ( t->map_entity_g2l.items ) { + free ( t->map_entity_g2l.items ); + } + t->map_entity_g2l.items = NULL; if ( t->vertices ) { free ( t->vertices ); @@ -532,6 +533,9 @@ h5t_add_num_vertices ( return HANDLE_H5_NOMEM_ERR; } + h5_err_t h5err = _h5_alloc_idmap (&t->map_vertex_g2l, num ); + if ( h5err < 0 ) return h5err; + return H5_SUCCESS; } @@ -698,7 +702,7 @@ h5t_get_num_vertices ( h5_id_t h5t_store_vertex ( h5_file * f, /*!< file handle */ - const h5_id_t id, /*!< global vertex id or -1 */ + const h5_id_t global_id, /*!< global vertex id or -1 */ const h5_float64_t P[3] /*!< coordinates */ ) { struct h5t_fdata *t = &f->t; @@ -720,20 +724,23 @@ h5t_store_vertex ( check id */ if ( (t->cur_level == 0) && ( - (id < 0) || (id >= t->num_vertices[0]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "vertex", id ); + (global_id < 0) || (global_id >= t->num_vertices[0]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "vertex", global_id ); } if ( (t->cur_level > 0) && ( - (id < t->num_vertices[t->cur_level-1]) || - (id >= t->num_vertices[t->cur_level]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "vertex", id ); + (global_id < t->num_vertices[t->cur_level-1]) || + (global_id >= t->num_vertices[t->cur_level]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "vertex", global_id ); } - h5_vertex *vertex = &t->vertices[++t->last_stored_vertex_id]; - vertex->id = id; + h5_id_t local_id = ++t->last_stored_vertex_id; + h5_vertex *vertex = &t->vertices[local_id]; + vertex->id = global_id; memcpy ( &vertex->P, P, sizeof ( vertex->P ) ); + + _h5_insert_idmap ( &t->map_vertex_g2l, global_id, local_id ); - return t->last_stored_vertex_id; + return local_id; } h5_err_t @@ -775,42 +782,190 @@ h5t_add_num_entities ( const h5_size_t num ) { struct h5t_fdata *t = &f->t; - ssize_t num_bytes = 0; + 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; + 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; + num + t->num_entities_on_level[t->cur_level-1] : num; switch ( t->mesh_type ) { case TETRAHEDRAL_MESH: - num_bytes = num_entities*sizeof ( t->entities.tets[0] ); + sizeof_entity = sizeof ( t->entities.tets[0] ); + sizeof_lentity = sizeof ( t->lentities.tets[0] ); break; case TRIANGLE_MESH: - num_bytes = num_entities*sizeof ( t->entities.tris[0] ); + sizeof_entity = sizeof ( t->entities.tris[0] ); + sizeof_lentity = sizeof ( t->lentities.tris[0] ); break; default: return -1; } - h5_debug ( "Allocating %ld bytes.", num_bytes ); - t->entities.data = realloc ( t->entities.data, num_bytes ); + + t->entities.data = realloc ( + t->entities.data, num_entities*sizeof_entity ); if ( t->entities.data == NULL ) { return H5_ERR_NOMEM; } + t->lentities.data = realloc ( + t->lentities.data, 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 ); - t->map_entity_g2l = realloc ( - t->map_entity_g2l, - num_entities*sizeof ( t->map_entity_g2l[0] ) ); + h5_err_t h5err = _h5_alloc_idmap (&t->map_entity_g2l, num_entities ); + if ( h5err < 0 ) return h5err; return H5_SUCCESS; } +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 */ + ) { + struct h5t_fdata *t = &f->t; + + h5_id_t *global_id = global_ids; + + 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; + + 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; +} h5_id_t h5t_store_tet ( h5_file * f, - const h5_id_t id, /*!< global tetrahedron id */ + const h5_id_t global_id, /*!< global tetrahedron id */ const h5_id_t parent_id, /*!< global parent id if level \c >0 else \c -1 */ const h5_id_t vertex_ids[4] /*!< tuple with vertex id's */ @@ -834,51 +989,55 @@ h5t_store_tet ( check parent id */ if ( (t->cur_level == 0) && (parent_id != -1) ) { - return HANDLE_H5_PARENT_ID_ERR ( "tet", id, parent_id ); + return HANDLE_H5_PARENT_ID_ERR ( "tet", global_id, parent_id ); } if ( (t->cur_level > 0) && (parent_id < 0) ) { - return HANDLE_H5_PARENT_ID_ERR ( "tet", id, parent_id ); + 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]) ) { - return HANDLE_H5_PARENT_ID_ERR ( "tet", id, parent_id ); + return HANDLE_H5_PARENT_ID_ERR ( "tet", global_id, parent_id ); } /* check id */ if ( (t->cur_level == 0) && ( - (id < 0) || (id >= t->num_entities[0]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "tet", id ); + (global_id < 0) || (global_id >= t->num_entities[0]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "tet", global_id ); } if ( (t->cur_level > 0) && ( - (id < t->num_entities[t->cur_level-1]) || - (id >= t->num_entities[t->cur_level]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "tet", id ); + (global_id < t->num_entities[t->cur_level-1]) || + (global_id >= t->num_entities[t->cur_level]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "tet", global_id ); } - - h5_tetrahedron *tet = &t->entities.tets[++t->last_stored_entity_id]; - tet->id = id; + h5_id_t local_id = ++t->last_stored_entity_id; + h5_tetrahedron *tet = &t->entities.tets[local_id]; + tet->id = global_id; tet->parent_id = parent_id; tet->refined_on_level = -1; tet->unused = 0; + memcpy ( &tet->vertex_ids, vertex_ids, sizeof ( tet->vertex_ids ) ); - t->map_entity_g2l[id] = t->last_stored_entity_id; + _h5t_isort_vertices ( 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 = t->map_entity_g2l[parent_id]; + 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->num_entities_on_level[t->cur_level]--; } } - return t->last_stored_entity_id; + return local_id; } h5_id_t h5t_store_triangle ( h5_file * f, - const h5_id_t id, /*!< global triangle id */ + const h5_id_t global_id, /*!< global triangle id */ const h5_id_t parent_id, /*!< global parent id if level \c >0 else \c -1 */ const h5_id_t vertex_ids[3] /*!< tuple with vertex id's */ @@ -902,43 +1061,45 @@ h5t_store_triangle ( check parent id */ if ( (t->cur_level == 0) && (parent_id != -1) ) { - return HANDLE_H5_PARENT_ID_ERR ( "triangle", 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", 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", id, parent_id ); + return HANDLE_H5_PARENT_ID_ERR ( "triangle", global_id, parent_id ); } /* check id */ if ( (t->cur_level == 0) && ( - (id < 0) || (id >= t->num_entities[0]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "triangle", id ); + (global_id < 0) || (global_id >= t->num_entities[0]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "triangle", global_id ); } if ( (t->cur_level > 0) && ( - (id < t->num_entities[t->cur_level-1]) || - (id >= t->num_entities[t->cur_level]) ) ) { - return HANDLE_H5_OUT_OF_RANGE_ERR( "triangle", id ); + (global_id < t->num_entities[t->cur_level-1]) || + (global_id >= t->num_entities[t->cur_level]) ) ) { + return HANDLE_H5_OUT_OF_RANGE_ERR( "triangle", global_id ); } - - h5_triangle *tri = &t->entities.tris[++t->last_stored_entity_id]; - tri->id = id; + h5_id_t local_id = ++t->last_stored_entity_id; + h5_triangle *tri = &t->entities.tris[local_id]; + tri->id = global_id; tri->parent_id = parent_id; tri->refined_on_level = -1; memcpy ( &tri->vertex_ids, vertex_ids, sizeof ( tri->vertex_ids ) ); - t->map_entity_g2l[id] = t->last_stored_entity_id; + _h5_insert_idmap ( &t->map_entity_g2l, global_id, local_id ); + if ( parent_id >= 0 ) { - h5_id_t local_parent_id = t->map_entity_g2l[parent_id]; + 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->num_entities_on_level[t->cur_level]--; } } - return t->last_stored_entity_id; + return local_id; } static h5_err_t