mapping of tets and triangles to global and local id added

.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

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

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

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

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()) ( \

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

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"

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

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;
 

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

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

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

src/h5/t_errorhandling.c

@@ -0,0 +1,41 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>	/* va_arg - System dependent ?! */
+#include <string.h>
+#include <errno.h>
+#include <fcntl.h>
+#include <hdf5.h>
+
+#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;
+}
+	

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

src/h5/t_map.c

@@ -0,0 +1,657 @@
+#include <hdf5.h>
+
+#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;
+}

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

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

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
 */

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] ) {

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

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