src/h5_core/Makefile.am

@@ -29,20 +29,7 @@ HDFLIB = -L$(HDF5ROOT)/lib -lhdf5 -lz $(SZLIB) @LDFLAGS@
 # SZ LIBRARY
 SZLIB = @SZLIB@
 
-# Extra files that I wish to include in the dist tar ball.
-EXTRA_DIST = $(EXTRA_HEADERS)
-
-# Files that I don't want to include in the dist tar ball
-nodist_include_HEADERS =
-
-# What to build... Will be determined by configure script.
-lib_LIBRARIES = libH5_core.a
-
-# Listing of all possible targets that I may build.
-EXTRA_LIBRARIES = libH5_core.a
-
-# Header files that I wish to install in $(prefix)/include
-include_HEADERS = \
+EXTRA_HEADERS = \
 	h5_attribs.h \
 	h5_core.h \
 	h5_errorhandling.h \
@@ -58,6 +45,8 @@ include_HEADERS = \
 	h5_syscall_private.h \
 	h5_types.h \
 	h5b_errorhandling_private.h \
+	h5t_adjacencies.h \
+	h5t_adjacencies_private.h \
 	h5t_boundaries.h \
 	h5t_boundaries_private.h \
 	h5t_consts_private.h \
@@ -73,8 +62,20 @@ include_HEADERS = \
 	h5u_readwrite.h \
 	hdf5_misc.h
 
-# Listing of all possible headers that I may include
-EXTRA_HEADERS =
+# Extra files that I wish to include in the dist tar ball.
+EXTRA_DIST = $(EXTRA_HEADERS)
+
+# Files that I don't want to include in the dist tar ball
+nodist_include_HEADERS =
+
+# What to build... Will be determined by configure script.
+lib_LIBRARIES = libH5_core.a
+
+# Listing of all possible targets that I may build.
+EXTRA_LIBRARIES = libH5_core.a
+
+# Header files that I wish to install in $(prefix)/include
+include_HEADERS =
 
 # Listing of sources
 libH5_core_a_SOURCES = \
@@ -89,6 +90,7 @@ libH5_core_a_SOURCES = \
 	h5_qsort_r.c \
 	h5_readwrite.c \
 	h5_syscall.c \
+	h5t_adjacencies.c \
 	h5t_boundaries.c \
 	h5t_consts.c \
 	h5t_errorhandling.c \
@@ -100,6 +102,8 @@ libH5_core_a_SOURCES = \
 	h5u_readwrite.c \
 	hdf5_misc.c
 
+libH5_core_a_DEPENDENCIES = $(EXTRA_HEADERS)
+
 all: libH5_core.a
 
 libH5_core.a: $(libH5_core_a_OBJECTS)
@@ -109,7 +113,7 @@ libH5_core.a: $(libH5_core_a_OBJECTS)
 %.o : %.c
 	$(CC) $(CFLAGS) $(INC) -c $<
 
-$(libH5_core_a_OBJECTS): h5_core.h
+$(libH5_core_a_OBJECTS): $(libH5_core_a_DEPENDENCIES)
 
 #errorhandling.o: errorhandling.c
 

src/h5_core/h5_core.h

@@ -9,6 +9,7 @@
 #include "h5_maps.h"
 #include "h5_openclose.h"
 #include "h5_readwrite.h"
+#include "h5t_adjacencies.h"
 #include "h5t_boundaries.h"
 #include "h5t_inquiry.h"
 #include "h5t_map.h"

src/h5_core/h5_core_private.h

@@ -16,6 +16,7 @@
 
 #include "h5b_errorhandling_private.h"
 
+#include "h5t_adjacencies_private.h"
 #include "h5t_boundaries_private.h"
 #include "h5t_consts_private.h"
 #include "h5t_errorhandling_private.h"
@@ -47,17 +48,14 @@
    2E TT TT TT TT TT TT
 */
 #define H5_TET_MASK		( (h5_id_t) (ULLONG_MAX >> 8) )
-#define _h5t_build_triangle_id( idx, eid ) \
+#define _h5t_build_triangle_id( idx, eid )	 \
 	( (  1ull << (sizeof(eid)*8-4)) |	 \
-	  (idx << (sizeof(eid)*7)) |	 \
+	  (idx << (sizeof(eid)*7)) |		 \
 	  (eid & H5_TET_MASK))
 
-#define TRY(func)						\
+#define TRY( func )					\
 	if ( (int64_t)(ptrdiff_t)(func) < (int64_t)0 )	\
 		return H5_ERR;
-#define TRY2(func,exception)			\
-	if ( (int64_t)(ptrdiff_t)(func) < (int64_t)0 )	\
-		goto exception;
 
 /*!
   The functions declared here are not part of the API, but may be used

src/h5_core/h5_errorhandling.h

@@ -2,8 +2,10 @@
 #define __H5_ERRORHANDLING_H
 
 #define SET_FNAME( f, fname )	h5_set_funcname( f, fname );
-#define CHECK_FILEHANDLE( f ) \
-	TRY ( h5_check_filehandle ( f ) );
+#define CHECK_FILEHANDLE( f )					\
+	if ( h5_check_filehandle ( f ) != H5_SUCCESS )		\
+		return h5_get_errno( f );
+
 
 #define CHECK_WRITABLE_MODE( f )					\
 	if ( f->mode==H5_O_RDONLY )					\

src/h5_core/h5_maps.c

@@ -7,23 +7,43 @@
 #include "h5_core_private.h"
 
 h5_err_t
-_h5_alloc_smap (
+_h5_alloc_idlist (
 	h5_file_t * const f,
-	struct smap	*map,
+	h5_idlist_t	*list,
 	const h5_size_t	size
 	) {
-	int new = ( map->items == NULL );
-	size_t size_in_bytes = size * sizeof ( map->items[0] );
-	TRY ( map->items = _h5_alloc ( f, map->items, size_in_bytes ) );
-	map->size = size;
-	if ( new ) map->num_items = 0;
+	int new = ( list->items == NULL );
+	size_t size_in_bytes = size * sizeof ( list->items[0] );
+	TRY ( list->items = _h5_alloc ( f, list->items, size_in_bytes ) );
+	list->size = size;
+	if ( new ) list->num_items = 0;
 	return H5_SUCCESS;
 }
 
+h5_err_t
+_h5_append_to_idlist (
+	h5_file_t * const f,
+	h5_idlist_t *list,
+	h5_id_t id
+	) {
+	if ( list->num_items == list->size ) {
+		h5_size_t size = list->size;
+		if ( size == 0 ) {
+			size = 16 * sizeof(list->items[0]);
+		} else {
+			size *= 2;
+		}
+		TRY ( list->items = _h5_alloc ( f, list->items, size ) );
+	}
+	list->items[list->num_items++] = id;
+	return H5_SUCCESS;
+}
+
+
 h5_err_t
 _h5_alloc_idmap (
 	h5_file_t * const f,
-	struct idmap	*map,
+	h5_idmap_t	*map,
 	const h5_size_t	size
 	) {
 	int new = ( map->items == NULL );
@@ -38,7 +58,7 @@ _h5_alloc_idmap (
 h5_err_t
 _h5_insert_idmap (
 	h5_file_t * const f,
-	struct idmap *map,
+	h5_idmap_t *map,
 	h5_id_t global_id,
 	h5_id_t local_id
 	) {
@@ -75,7 +95,7 @@ _h5_insert_idmap (
  */
 h5_id_t
 _h5_search_idmap (
-	struct idmap *map,
+	h5_idmap_t *map,
 	h5_id_t value
 	) {
 
@@ -105,7 +125,7 @@ _cmp_idmap (
 
 h5_err_t
 _h5_sort_idmap (
-	struct idmap *map
+	h5_idmap_t *map
 	) {
 	qsort ( map->items, map->num_items, sizeof(map->items[0]), _cmp_idmap );
 	return H5_SUCCESS;

src/h5_core/h5_maps.h

@@ -2,36 +2,43 @@
 #define __H5_MAPS_H
 
 h5_err_t
-_h5_alloc_smap (
+_h5_alloc_idlist (
 	h5_file_t * const f,
-	struct smap	*map,
+	h5_idlist_t *list,
 	const h5_size_t	size
 	);
 
+h5_err_t
+_h5_append_to_idlist (
+	h5_file_t * const f,
+	h5_idlist_t *list,
+	h5_id_t id
+	);
+
 h5_err_t
 _h5_alloc_idmap (
 	h5_file_t * const f,
-	struct idmap	*map,
+	h5_idmap_t *map,
 	const h5_size_t	size
 	);
 
 h5_err_t
 _h5_insert_idmap (
 	h5_file_t * const f,
-	struct idmap *map,
+	h5_idmap_t *map,
 	h5_id_t global_id,
 	h5_id_t local_id
 	);
 
 h5_id_t
 _h5_search_idmap (
-	struct idmap *map,
+	h5_idmap_t *map,
 	h5_id_t value
 	);
 
 h5_err_t
 _h5_sort_idmap (
-	struct idmap *map
+	h5_idmap_t *map
 	);
 
 #endif

src/h5_core/h5_mpi.c

@@ -7,7 +7,7 @@
 h5_err_t
 _h5_mpi_allgather (
 	h5_file_t * const f,
-	const void * sendbuf,
+	void * sendbuf,
 	const int sendcount,
 	const MPI_Datatype sendtype,
 	void * recvbuf,
@@ -34,7 +34,7 @@ _h5_mpi_allgather (
 h5_err_t
 _h5_mpi_comm_size (
 	h5_file_t * const f,
-	MPI_comm comm,
+	MPI_Comm comm,
 	int *size
 	) {
 	int err = MPI_Comm_size ( comm, size );
@@ -51,7 +51,7 @@ _h5_mpi_comm_size (
 h5_err_t
 _h5_mpi_comm_rank (
 	h5_file_t * const f,
-	MPI_comm comm,
+	MPI_Comm comm,
 	int *rank
 	) {
 	int err = MPI_Comm_rank ( comm, rank );

src/h5_core/h5_mpi_private.h

@@ -5,7 +5,7 @@
 h5_err_t
 _h5_mpi_allgather (
 	h5_file_t * const f,
-	const void * sendbuf,
+	void * sendbuf,
 	const int sendcount,
 	const MPI_Datatype sendtype,
 	void * recvbuf,
@@ -17,14 +17,14 @@ _h5_mpi_allgather (
 h5_err_t
 _h5_mpi_comm_size (
 	h5_file_t * const f,
-	MPI_comm comm,
+	MPI_Comm comm,
 	int *size
 	);
 
 h5_err_t
 _h5_mpi_comm_rank (
 	h5_file_t * const f,
-	MPI_comm comm,
+	MPI_Comm comm,
 	int *rank
 	);
 #endif

src/h5_core/h5_openclose.c

@@ -63,7 +63,7 @@ static h5_int64_t
 _h5u_open_file (
 	h5_file_t *f			/*!< IN: file handle */
 	) {
-	TRY ( f->u = (struct h5u_fdata*) _h5_alloc( f, NULL, sizeof (*f->u) ) );
+	TRY ( f->u = (h5u_fdata_t*) _h5_alloc( f, NULL, sizeof (*f->u) ) );
 	h5u_fdata_t *u = f->u;
 
  	u->shape = 0;
@@ -146,11 +146,7 @@ _h5_open_file (
 #ifdef PARALLEL_IO
 	f->comm = comm;
 	TRY ( _h5_mpi_comm_size ( f, comm, &f->nprocs ) );
-
-	if (MPI_Comm_rank (comm, &f->myproc) != MPI_SUCCESS) {
-		return HANDLE_MPI_COMM_RANK_ERR;
-	}
-	
+	TRY ( _h5_mpi_comm_rank ( f, comm, &f->myproc ) );
 	
 	/* for the SP2... perhaps different for linux */
 	MPI_Info info = MPI_INFO_NULL;
@@ -158,7 +154,7 @@ _h5_open_file (
 	/* ks: IBM_large_block_io */
 	MPI_Info_create(&info);
 	MPI_Info_set(info, "IBM_largeblock_io", "true" );
-	TRY ( _h5_set_fapl_mpio_property ( f->access_prop, comm, info ) );
+	TRY ( _h5_set_fapl_mpio_property ( f, f->access_prop, comm, info ) );
 	MPI_Info_free(&info);
 	
 	TRY ( f->access_prop = _h5_create_property ( f, H5P_FILE_ACCESS ) );

src/h5_core/h5_syscall.c

@@ -1,5 +1,6 @@
 #include <stdlib.h>
 #include <hdf5.h>
+#include <search.h>
 #include "h5_core.h"
 #include "h5_core_private.h"
 
@@ -19,3 +20,43 @@ _h5_alloc (
 	}
 	return ptr;
 }
+
+void
+_h5_free (
+	h5_file_t * const f,
+	void *ptr
+	) {
+	if ( ptr ) free ( ptr );
+}
+	
+void *
+_h5_tsearch (
+	h5_file_t * const f,
+	const void *key,
+	void **rootp,
+	int (*compar) (const void *key1, const void *key2) 
+	) {
+	void *ptr = tsearch ( key, rootp, compar );
+	if ( ptr == NULL ) {
+		h5_error (
+			f,
+			H5_ERR_NOMEM,
+			"Out of memory." );
+		return (void*)(-1);
+	}
+	return ptr;
+}
+
+void *
+_h5_tfind (
+	h5_file_t * const f,
+	const void *key,
+	void *const *rootp,
+	int (*compar) (const void *key1, const void *key2) 
+	) {
+	void *ptr = tfind ( key, rootp, compar );
+	if ( ptr == NULL ) {
+		return (void*)(-1);
+	}
+	return ptr;
+}

src/h5_core/h5_syscall_private.h

@@ -5,6 +5,28 @@ void *
 _h5_alloc (
 	h5_file_t * const f,
 	void *ptr,
-	const size_t size );
+	const size_t size
+	);
 
+void
+_h5_free (
+	h5_file_t * const f,
+	void *ptr
+	);
+
+void *
+_h5_tsearch (
+	h5_file_t * const f,
+	const void *key,
+	void **rootp,
+	int (*compar) (const void *key1, const void *key2) 
+	);
+
+void *
+_h5_tfind (
+	h5_file_t * const f,
+	const void *key,
+	void *const *rootp,
+	int (*compar) (const void *key1, const void *key2) 
+	);
 #endif

src/h5_core/h5_types.h

@@ -54,9 +54,10 @@ typedef h5_err_t (*h5_errorhandler_t)(
 typedef unsigned long		MPI_Comm;
 #endif
 
-struct smap;
-typedef struct smap smap_t;
-struct idmap;
+struct h5_idlist;
+typedef struct h5_idlist h5_idlist_t;
+struct h5_idmap;
+typedef struct h5_idmap h5_idmap_t;
 
 /**
    \struct h5_file

src/h5_core/h5_types_private.h

@@ -1,7 +1,7 @@
 #ifndef __H5_TYPES_PRIVATE_H
 #define __H5_TYPES_PRIVATE_H
 
-struct smap {
+struct h5_idlist {
 	h5_size_t	size;		/* allocated space in number of items */
 	h5_size_t	num_items;	/* stored items	*/
 	h5_id_t		*items;
@@ -13,11 +13,9 @@ struct h5_idmap_el {
 };
 typedef struct h5_idmap_el h5_idmap_el_t;
 
-struct idmap {
+struct h5_idmap {
 	h5_size_t	size;		/* allocated space in number of items */
 	h5_size_t	num_items;	/* stored items	*/
 	h5_idmap_el_t	*items;
 };
-
-
 #endif

src/h5_core/h5t_errorhandling.c

@@ -8,22 +8,6 @@
 
 #include "h5_core.h"
 #include "h5_core_private.h"
-
-h5_err_t
-_h5t_handle_get_global_eid_err (
-	h5_file_t *f,
-	const h5_id_t * const global_vids
-	) {
-	struct h5t_fdata *t = f->t;
-	switch ( t->mesh_type ) {
-	case H5_OID_TETRAHEDRON:
-		return _h5t_error_global_tet_id_nexist ( f, global_vids );
-	case H5_OID_TRIANGLE:
-		return _h5t_error_global_tri_id_nexist ( f, global_vids );
-	default:
-		return h5_error_internal( f, __FILE__, __func__, __LINE__ );
-	}
-}
 	
 h5_err_t
 _h5t_error_illegal_object_type (

src/h5_core/h5t_errorhandling_private.h

@@ -1,13 +1,6 @@
 #ifndef __H5T_ERRORHANDLING_PRIVATE_H
 #define __H5T_ERRORHANDLING_PRIVATE_H
 
-
-h5_err_t
-_h5t_handle_get_global_eid_err (
-	h5_file_t *f,
-	const h5_id_t * const global_vids
-	);
-
 h5_err_t
 _h5t_error_illegal_object_type (
 	h5_file_t * const f,

src/h5_core/h5t_map.c

@@ -282,7 +282,8 @@ _qsort_cmp_elems1 (
 }
 
 /*!
-  Sort elems geometrically.
+  Sort elements. Store local id's in a sorted array so we can run a
+  binary search. 
 */
 h5_err_t
 _h5t_sort_elems (
@@ -296,7 +297,7 @@ _h5t_sort_elems (
 	int k;
 	h5_id_t i;
 	for ( k = 0; k < 2; k++ ) {
-		TRY( _h5_alloc_smap ( f, &t->sorted_elems[k], num_elems ) );
+		TRY( _h5_alloc_idlist ( f, &t->sorted_elems[k], num_elems ) );
 		for ( i = local_eid; i < num_elems; i++ ) {
 			t->sorted_elems[k].items[i] = i;
 		}
@@ -518,16 +519,14 @@ h5t_get_global_eid (
 	if ( t->vertices == NULL ) {
 		TRY ( _h5t_read_mesh ( f ) );
 	}
-	TRY2 ( h5t_map_global_vids2local (
+	TRY ( h5t_map_global_vids2local (
 		f,
 		global_vids,
 		t->mesh_type,
-		local_vids ), fail );
+		local_vids ) );
 	h5_id_t local_eid;
-	TRY2 ( local_eid = h5t_get_local_eid ( f, local_vids ), fail );
+	TRY ( local_eid = h5t_get_local_eid ( f, local_vids ) );
 	return h5t_map_local_eid2global ( f, local_eid );
-fail:
-	return _h5t_handle_get_global_eid_err ( f, global_vids );
 }
 
 static int
@@ -695,10 +694,11 @@ h5t_get_local_triangle_id  (
 	h5_id_t * const local_vids
 	) {
 	struct h5t_fdata *t = f->t;
-	
+	if ( t->vertices == NULL ) {
+		TRY ( _h5t_read_mesh ( f ) );
+	}
 	switch ( t->mesh_type ) {
 	case H5_OID_TETRAHEDRON: {
-		TRY ( _h5t_read_mesh ( f ) );
 		h5_id_t local_tid = _tetm_search_triangle ( f, local_vids );
 		if ( local_tid == -1 ) {
 			return _h5t_error_local_triangle_id_nexist( f, local_vids );
@@ -844,3 +844,79 @@ _h5t_rebuild_global_2_local_map_of_elems (
 
 	return H5_SUCCESS;
 }
+
+/*
+  setup structure "elems_data" with local ids for each element:
+  - translate the global vertex id's of each element to their
+    local id's
+  - translate the global parent id of each element to the
+    corresponding local id.
+*/
+h5_err_t
+_h5t_rebuild_elems_data (
+	h5_file_t * const f
+	) {
+	h5t_fdata_t *t = f->t;
+	h5_id_t local_eid = 0;
+	h5_id_t num_elems = t->num_elems[t->num_levels-1];
+	h5_id_t level_id = 0;
+	switch ( t->mesh_type ) {
+	case H5_OID_TETRAHEDRON: {
+		h5_tet_t *el;
+		h5_tet_data_t *el_data;
+
+		for ( local_eid=0, el=t->elems.tets, el_data=t->elems_data.tets;
+		      local_eid < num_elems;
+		      local_eid++, el++, el_data++ ) {
+			
+			TRY( h5t_map_global_vids2local (
+				     f,
+				     el->global_vids,
+				     t->mesh_type,
+				     el_data->local_vids
+				     ) );
+			if ( el->global_parent_eid >= 0 )
+				TRY ( el_data->local_parent_eid =
+				      h5t_map_global_eid2local (
+					      f,
+					      el->global_parent_eid ) );
+			if ( local_eid > t->num_elems_on_level[level_id] ) {
+				level_id++;
+			}
+			el_data->level_id = level_id;
+		}
+		break;
+	}
+	case H5_OID_TRIANGLE: {
+		h5_triangle_t *el;
+		h5_triangle_data_t *el_data;
+
+		for ( local_eid=0, el=t->elems.tris, el_data=t->elems_data.tris;
+		      local_eid < num_elems;
+		      local_eid++, el++, el_data++ ) {
+
+			TRY( h5t_map_global_vids2local (
+				     f,
+				     el->global_vids,
+				     t->mesh_type,
+				     el_data->local_vids
+				     ) );
+			if ( el->global_parent_eid >= 0 )
+				TRY ( el_data->local_parent_eid =
+				      h5t_map_global_eid2local (
+					      f,
+					      el->global_parent_eid ) );
+			if ( local_eid > t->num_elems_on_level[level_id] ) {
+				level_id++;
+			}
+			el_data->level_id = level_id;
+		}
+		break;
+	}
+	default:
+		return h5_error_internal (
+			f, __FILE__, __func__, __LINE__ );
+	}
+
+	return H5_SUCCESS;
+}

src/h5_core/h5t_map_private.h

@@ -34,4 +34,10 @@ h5_err_t
 _h5t_rebuild_global_2_local_map_of_elems (
 	h5_file_t * const f
 	);
+
+h5_err_t
+_h5t_rebuild_elems_data (
+	h5_file_t * const f
+	);
+
 #endif

src/h5_core/h5t_openclose.c

@@ -474,40 +474,23 @@ _release_memory (
 	) {
 	struct h5t_fdata *t = f->t;
 
-	if ( t->num_vertices ) {
-		free ( t->num_vertices );
-	}
-	t->num_vertices = NULL;
-
-	if ( t->num_elems ) {
-		free ( t->num_elems );
-	}
-	t->num_elems = NULL;
-
-	if ( t->num_elems_on_level ) {
-		free ( t->num_elems_on_level );
-	}
-	t->num_elems_on_level = NULL;
-
-	if ( t->map_vertex_g2l.items ) {
-		free ( t->map_vertex_g2l.items );
-	}
-	t->map_vertex_g2l.items = NULL;
-
-	if ( t->map_elem_g2l.items ) {
-		free ( t->map_elem_g2l.items );
-	}
-	t->map_elem_g2l.items = NULL;
-
-	if ( t->vertices ) {
-		free ( t->vertices );
-	}
+	_h5_free ( f, t->vertices );
 	t->vertices = NULL;
+	_h5_free ( f, t->vertices_data );
+	t->vertices_data = NULL;
+	_h5_free ( f, t->num_vertices );
+	t->num_vertices = NULL;
+	_h5_free ( f, t->map_vertex_g2l.items );
 
-	if ( t->elems.data ) {
-		free ( t->elems.data );
-	}
+	_h5_free ( f, t->elems.data );
 	t->elems.data = NULL;
+	_h5_free ( f, t->num_elems );
+	_h5_free ( f, t->elems_data.data );
+	t->elems_data.data = NULL;
+	_h5_free ( f, t->num_elems_on_level );
+	t->num_elems_on_level = NULL;
+	_h5_free ( f, t->map_elem_g2l.items );
+	t->map_elem_g2l.items = NULL;
 
 	return H5_SUCCESS;
 }

src/h5_core/h5t_readwrite.c

@@ -101,9 +101,6 @@ _h5t_write_mesh (
 	) {
 	h5t_fdata_t *t = f->t;
 	if ( ! t->mesh_changed ) return 0;
-	if ( t->level_changed ) {
-		TRY ( _h5t_close_level ( f ) );
-	}
 
 	TRY( _write_vertices( f ) );
 	TRY( _write_elems( f ) );
@@ -292,64 +289,11 @@ _h5t_read_elems (
 
 	TRY ( _h5t_sort_elems ( f ) );
 	TRY ( _h5t_rebuild_global_2_local_map_of_elems ( f ) );
-	/*
-	  setup structure with local vertex ids
-	 */
-	h5_id_t local_eid = 0;
-	h5_id_t num_elems = t->num_elems[t->num_levels-1];
-	switch ( t->mesh_type ) {
-	case H5_OID_TETRAHEDRON: {
-		h5_tet_t *el;
-		h5_tet_data_t *el_data;
-
-		for ( local_eid = 0, el = t->elems.tets, el_data = t->elems_data.tets;
-		      local_eid < num_elems;
-		      local_eid++, el++, el_data++ ) {
-
-			TRY( h5t_map_global_vids2local (
-				     f,
-				     el->global_vids,
-				     t->mesh_type,
-				     el_data->local_vids
-				     ) );
-			if ( el->global_parent_eid >= 0 )
-				TRY ( el_data->local_parent_eid =
-				      h5t_map_global_eid2local (
-					      f,
-					      el->global_parent_eid ) );
-		}
-		break;
-	}
-	case H5_OID_TRIANGLE: {
-		h5_triangle_t *el;
-		h5_triangle_data_t *el_data;
-
-		for ( local_eid = 0, el = t->elems.tris, el_data = t->elems_data.tris;
-		      local_eid < num_elems;
-		      local_eid++, el++, el_data++ ) {
-
-			TRY( h5t_map_global_vids2local (
-				     f,
-				     el->global_vids,
-				     t->mesh_type,
-				     el_data->local_vids
-				     ) );
-			if ( el->global_parent_eid >= 0 )
-				TRY ( el_data->local_parent_eid =
-				      h5t_map_global_eid2local (
-					      f,
-					      el->global_parent_eid ) );
-		}
-		break;
-	}
-	default:
-		return -1;
-	}
-	
+	TRY ( _h5t_rebuild_elems_data ( f ) );
+	TRY ( _h5t_rebuild_adj_data ( f ) );
 	return H5_SUCCESS;
 }
 
-
 h5_err_t
 h5t_start_traverse_elems (
 	h5_file_t * f
@@ -442,19 +386,10 @@ h5_err_t
 _h5t_read_mesh (
 	h5_file_t *f
 	) {
-	struct h5t_fdata *t = f->t;
 
-	if ( t->vertices == NULL ) {
-		TRY ( _h5t_read_vertices ( f ) );
-	}
+	TRY ( _h5t_read_vertices ( f ) );
+	TRY ( _h5t_read_elems ( f ) );
 
-	if ( t->elems.data == NULL ) {
-		TRY ( _h5t_read_elems ( f ) );
-	}
-
-	if ( t->sorted_elems[0].items == NULL ) {
-		_h5t_sort_elems ( f );
-	}
 	return H5_SUCCESS;
 }
 

src/h5_core/h5t_readwrite.h

@@ -46,7 +46,7 @@ h5t_add_num_triangles (
 	) ;
 
 h5_size_t
-h5t_add_num_elems (
+h5t_begin_store_elems (
 	h5_file_t * f,
 	const h5_size_t num
 	) ;

src/h5_core/h5t_storemesh.c

@@ -19,8 +19,8 @@ h5t_add_mesh (
 	) {
 	h5_id_t mesh_id = 0;
 
-	TRY( (mesh_id = h5t_open_mesh ( f, -1, mesh_type )) ); 
-	TRY( h5t_add_level ( f ) );
+	TRY ( (mesh_id = h5t_open_mesh ( f, -1, mesh_type )) ); 
+	TRY ( h5t_add_level ( f ) );
 	
 	f->t->mesh_changed = 1;
 
@@ -33,15 +33,9 @@ h5t_add_mesh (
    unique id's assigned by the mesher but this id's may not be
    consecutive numbered starting from 0.
  * Set the global vertex id's in element definitions.
- * Assign unique global id's to elements. This is a NOOP if we are
-   running serial.
-
-   Where are the global vertex id's stored?
-   - Vertex definitions
-
 */
-herr_t
-_h5t_assign_global_ids (
+h5_err_t
+_h5t_assign_global_vertex_ids (
 	h5_file_t * const f
 	) {
 	h5t_fdata_t *t = f->t;
@@ -49,41 +43,30 @@ _h5t_assign_global_ids (
 
 	if ( t->cur_level < 0 ) return H5_SUCCESS; /* no level defined */
 
-	/*
-	  assign global id to vertices
-	*/
 	/*
 	  simple in serial runs: global_id = local_id
 	*/
-	for ( local_id = 0; local_id < t->num_vertices[t->num_levels-1]; local_id++ ) {
+	for ( local_id = 0;
+	      local_id < t->num_vertices[t->num_levels-1];
+	      local_id++ ) {
 		t->vertices[local_id].global_vid = local_id;
 	}
 
 	return H5_SUCCESS;
 }
 
+/*!
+ Assign unique global id's to elements.
+
+ Nothing to do in serial case.
+*/
 h5_err_t
-_h5t_close_level (
+_h5t_assign_global_elem_ids (
 	h5_file_t * const f
 	) {
-	h5t_fdata_t *t = f->t;
-
-	if ( t->num_levels <= 0 ) return H5_SUCCESS;
-
-	t->num_vertices[t->cur_level] = t->last_stored_vid+1;
-	TRY ( _h5t_assign_global_ids ( f ) );
-
-	TRY ( _h5t_sort_vertices ( f ) );
-	TRY ( _h5t_rebuild_global_2_local_map_of_vertices ( f ) );
-	TRY ( _h5t_sort_elems ( f ) );
-	TRY ( _h5t_rebuild_global_2_local_map_of_elems ( f ) );
-
 	return H5_SUCCESS;
 }
 
-#if 0
-
-#endif
 h5_id_t
 h5t_add_level (
 	h5_file_t * const f
@@ -94,8 +77,6 @@ h5t_add_level (
 		return H5_ERR_INVAL;
 	}
 
-	TRY ( _h5t_close_level ( f ) );
-
 	/* t->num_levels will be set to zero on file creation(!) */
 	if ( t->num_levels == -1 ) {	/* unknown number of levels	*/
 		/* determine number of levels */
@@ -134,8 +115,10 @@ _h5t_alloc_num_vertices (
 
 	ssize_t size = num_vertices * sizeof ( t->vertices[0] );
 	TRY ( t->vertices = _h5_alloc (	f, t->vertices, size ) );
+	size = num_vertices * sizeof ( t->vertices_data[0] );
+	TRY ( t->vertices_data = _h5_alloc (	f, t->vertices_data, size ) );
 	TRY( _h5_alloc_idmap ( f, &t->map_vertex_g2l, num_vertices ) );
-	TRY( _h5_alloc_smap ( f, &t->sorted_lvertices, num_vertices ) );
+	TRY( _h5_alloc_idlist ( f, &t->sorted_lvertices, num_vertices ) );
 
 	return H5_SUCCESS;
 }
@@ -144,15 +127,16 @@ _h5t_alloc_num_vertices (
   Allocate memory for (more) vertices.
 */
 h5_err_t
-h5t_add_num_vertices (
+h5t_begin_store_vertices (
 	h5_file_t * const f,
 	const h5_size_t num
 	) {
-	struct h5t_fdata *t = f->t;
+	h5t_fdata_t *t = f->t;
 
 	if ( t->cur_level < 0 ) {
 		return _h5t_error_undef_level( f );
 	}
+	t->storing_data = 1;
 	h5_size_t cur_num_vertices = ( t->cur_level > 0 ?
 				       t->num_vertices[t->cur_level-1] : 0 );
 	t->num_vertices[t->cur_level] = cur_num_vertices+num;
@@ -190,13 +174,28 @@ h5t_store_vertex (
 	return local_id;
 }
 
+h5_err_t
+h5t_end_store_vertices (
+	h5_file_t * const f
+	) {
+	h5t_fdata_t *t = f->t;
+	t->storing_data = 0;
+
+	t->num_vertices[t->cur_level] = t->last_stored_vid+1;
+	TRY ( _h5t_assign_global_vertex_ids ( f ) );
+	TRY ( _h5t_sort_vertices ( f ) );
+	TRY ( _h5t_rebuild_global_2_local_map_of_vertices ( f ) );
+
+	return H5_SUCCESS;
+}
+
 h5_err_t
 _h5t_alloc_num_elems (
 	h5_file_t * const f,
 	const size_t cur_num_elems,
 	const size_t new_num_elems
 	) {
-	struct h5t_fdata *t = f->t;
+	h5t_fdata_t *t = f->t;
 	size_t sizeof_elem = 0;
 	size_t sizeof_lelem = 0;
 
@@ -228,12 +227,13 @@ _h5t_alloc_num_elems (
 }
 
 h5_err_t
-h5t_add_num_elems (
+h5t_begin_store_elems (
 	h5_file_t * const f,
 	const h5_size_t num
 	) {
 	struct h5t_fdata *t = f->t;
 
+	t->storing_data = 1;
 	size_t cur_num_elems = t->cur_level > 0 ?
 		t->num_elems[t->cur_level-1] : 0;
 	size_t new_num_elems = t->cur_level > 0 ?
@@ -257,7 +257,7 @@ h5t_store_elem    (
 	const h5_id_t local_parent_eid,
 	const h5_id_t local_vids[]
 	) {
-	struct h5t_fdata *t = f->t;
+	h5t_fdata_t *t = f->t;
 
 	/* level set? */
 	if ( t->cur_level < 0 )
@@ -317,7 +317,8 @@ _h5t_store_tet (
 
 	memcpy ( &tet->global_vids, local_vids, sizeof ( tet->global_vids ) );
 	_h5t_sort_local_vids ( f, tet->global_vids, 4 );
-	memcpy ( &tet_data->local_vids, &tet->global_vids, sizeof ( tet->global_vids ) );
+	memcpy ( &tet_data->local_vids, &tet->global_vids,
+		 sizeof ( tet->global_vids ) );
 
 	if ( local_parent_eid >= 0 ) {
 		if ( t->elems.tets[local_parent_eid].refined_on_level < 0 ) {
@@ -348,7 +349,8 @@ _h5t_store_triangle (
 
 	memcpy ( &tri->global_vids, vids, sizeof ( tri->global_vids ) );
 	_h5t_sort_local_vids ( f, tri->global_vids, 3 );
-	memcpy ( &tri_data->local_vids, &tri->global_vids, sizeof ( tri->global_vids ) );
+	memcpy ( &tri_data->local_vids, &tri->global_vids,
+		 sizeof ( tri->global_vids ) );
 
 	if ( local_parent_eid >= 0 ) {
 		if ( t->elems.tris[local_parent_eid].refined_on_level < 0 ) {
@@ -360,14 +362,33 @@ _h5t_store_triangle (
 	return local_eid;
 }
 
+h5_err_t
+h5t_end_store_elems (
+	h5_file_t * const f
+	) {
+	h5t_fdata_t *t = f->t;
+	t->storing_data = 0;
+
+	t->num_elems[t->cur_level] = t->last_stored_eid+1;
+	TRY ( _h5t_assign_global_elem_ids ( f ) );
+	TRY ( _h5t_sort_elems ( f ) );
+	TRY ( _h5t_rebuild_global_2_local_map_of_elems ( f ) );
+	TRY ( _h5t_rebuild_elems_data ( f ) );
+
+	return H5_SUCCESS;
+}
+
 
 h5_err_t
-h5t_refine_num_elems (
+h5t_begin_refine_elems (
 	h5_file_t * const f,
 	const h5_size_t num_elems_to_refine
 	) {
 	h5_size_t num_elems_to_add = 0;
-	switch ( f->t->mesh_type ) {
+	h5t_fdata_t *t = f->t;
+
+	t->storing_data = 1;
+	switch ( t->mesh_type ) {
 	case H5_OID_TETRAHEDRON:
 		num_elems_to_add = num_elems_to_refine*8;
 		break;
@@ -378,8 +399,21 @@ h5t_refine_num_elems (
 		return h5_error_internal ( f, __FILE__, __func__, __LINE__ );
 	}
 	h5_size_t num_vertices = (num_elems_to_add>>2)*3; /* this is an upper limit */
-	TRY ( h5t_add_num_vertices ( f, num_vertices ) );
-	TRY ( h5t_add_num_elems ( f, num_elems_to_add ) );
+	TRY ( h5t_begin_store_vertices ( f, num_vertices ) );
+	TRY ( h5t_begin_store_elems ( f, num_elems_to_add ) );
+
+	return H5_SUCCESS;
+}
+
+h5_err_t
+h5t_end_refine_elems (
+	h5_file_t * const f
+	) {
+	h5t_fdata_t *t = f->t;
+	t->storing_data = 0;
+
+	TRY ( h5t_end_store_vertices ( f ) );
+	TRY ( h5t_end_store_elems ( f ) );
 
 	return H5_SUCCESS;
 }

src/h5_core/h5t_storemesh.h

@@ -13,13 +13,7 @@ h5t_add_level (
 	);
 
 h5_size_t
-h5t_add_num_vertices (
-	h5_file_t * const f,
-	const h5_size_t num
-	);
-
-h5_err_t
-h5t_add_num_elems (
+h5t_begin_store_vertices (
 	h5_file_t * const f,
 	const h5_size_t num
 	);
@@ -31,6 +25,17 @@ h5t_store_vertex (
 	const h5_float64_t P[3]
 	);
 
+h5_err_t
+h5t_end_store_vertices (
+	h5_file_t * const f
+	);
+
+h5_err_t
+h5t_begin_store_elems (
+	h5_file_t * const f,
+	const h5_size_t num
+	);
+
 h5_id_t
 h5t_store_elem    (
 	h5_file_t * const f,
@@ -39,7 +44,12 @@ h5t_store_elem    (
 	);
 
 h5_err_t
-h5t_refine_num_elems (
+h5t_end_store_elems (
+	h5_file_t * const f
+	);
+
+h5_err_t
+h5t_begin_refine_elems (
 	h5_file_t * const f,
 	const h5_size_t num_elems_to_refine
 	);
@@ -50,4 +60,5 @@ h5t_refine_elem (
 	const h5_id_t local_eid
 	);
 
+
 #endif

src/h5_core/h5t_types_private.h

@@ -11,6 +11,11 @@ struct h5_vertex {
 };
 typedef struct h5_vertex	h5_vertex_t;
 
+struct h5_vertex_data {
+	h5_idlist_t	tv;
+};
+typedef struct h5_vertex_data h5_vertex_data_t;
+
 struct h5_triangle {
 	h5_id_t		global_eid;
 	h5_id_t		global_parent_eid;
@@ -40,12 +45,14 @@ typedef struct h5_element	h5_element_t;
 
 struct h5_triangle_data {
 	h5_id_t		local_parent_eid;
+	h5_id_t		level_id;
 	h5_3id_t	local_vids;
 };
 typedef struct h5_triangle_data h5_triangle_data_t;
 
 struct h5_tet_data {
 	h5_id_t		local_parent_eid;
+	h5_id_t		level_id;
 	h5_4id_t	local_vids;
 };
 typedef struct h5_tet_data	h5_tet_data_t;
@@ -118,6 +125,37 @@ struct h5_dtypes {
 };
 typedef struct h5_dtypes h5_dtypes_t;
 
+struct h5_te_node_key {
+	h5_id_t vid0;
+	h5_id_t vid1;
+};
+typedef struct h5_te_node_key h5_te_node_key_t;
+
+struct h5_te_node {
+	h5_te_node_key_t key;
+	h5_idlist_t value;
+};
+typedef struct h5_te_node h5_te_node_t;
+
+struct h5_td_node_key {
+	h5_id_t vid0;
+	h5_id_t vid1;
+	h5_id_t vid2;
+};
+typedef struct h5_td_node_key h5_td_node_key_t;
+
+struct h5_td_node {
+	h5_td_node_key_t key;
+	h5_idlist_t value;
+};
+typedef struct h5_td_node h5_td_node_t;
+
+struct h5t_adjacencies {
+	void * te_tree;
+	void * td_tree;
+};
+typedef struct h5t_adjacencies h5t_adjacencies_t;
+
 struct h5t_fdata {
 	/*** book-keeping ***/
 	char		mesh_name[16];
@@ -131,12 +169,14 @@ struct h5t_fdata {
 	h5_id_t		new_level;	/* idx of the first new level or -1 */
 	h5_size_t	num_levels;	/* number of levels */
 	h5_id_t		level_changed;
+	h5_id_t		storing_data;
 
 	/*** vertices ***/
 	h5_vertex_t	*vertices;
+	h5_vertex_data_t *vertices_data;
 	h5_size_t	*num_vertices;
-	struct idmap	map_vertex_g2l;	/* map global id to local id */
-	struct smap	sorted_lvertices;
+	h5_idmap_t	map_vertex_g2l;	/* map global id to local id */
+	h5_idlist_t	sorted_lvertices;
 	h5_id_t		last_retrieved_vid; 
 	h5_id_t		last_stored_vid;
 	h5_dataset_info_t	dsinfo_vertices;
@@ -148,7 +188,7 @@ struct h5t_fdata {
 	h5_elems_data_t	elems_data;	/* local, per element data */
 	h5_size_t	*num_elems;
 	h5_size_t	*num_elems_on_level;
-	struct idmap	map_elem_g2l;	/* map global id to local id */
+	h5_idmap_t	map_elem_g2l;	/* map global id to local id */
 
 	/*
 	  array with geometrically sorted local entitiy ids
@@ -156,7 +196,7 @@ struct h5t_fdata {
 	  [1]: 1,0,2,3 sorted
 	  ...
 	*/
-	struct smap	sorted_elems[H5_MAX_VERTICES_PER_ELEM];
+	h5_idlist_t	sorted_elems[H5_MAX_VERTICES_PER_ELEM];
 
 	h5_id_t		last_retrieved_eid; 
 	h5_id_t		last_stored_eid;
@@ -169,6 +209,7 @@ struct h5t_fdata {
 	h5_id_t		boundaries_gid;		/* hdf5 grp id container group */
 
 	boundary_t	boundary;
+	h5t_adjacencies_t adjacencies;
 
 	/*** HDF5 objects ***/
 	hid_t		topo_gid;		/* grp id of mesh in current

src/h5_core/h5u_readwrite.c

@@ -209,15 +209,17 @@ h5u_set_num_elements (
 	*/
 
 	TRY ( _h5_mpi_allgather (
-		&nparticles, 1, MPI_LONG_LONG,
-		f->pnparticles, 1, MPI_LONG_LONG,
-		f->comm ) );
+		      f,
+		      &nparticles, 1, MPI_LONG_LONG,
+		      f->u->pnparticles, 1, MPI_LONG_LONG,
+		      f->comm ) );
 
 	if ( f->myproc == 0 ) {
-		h5_debug ( "Particle offsets:" );
+		h5_debug ( f, "Particle offsets:" );
 		for(i=0;i<f->nprocs;i++) 
-			h5_debug ( "\tnp=%lld",
-					      (long long) f->pnparticles[i] );
+			h5_debug ( f, 
+				   "\tnp=%lld",
+				   (long long) f->u->pnparticles[i] );
 	}
 	/* should I create a selection here? */
 
@@ -225,29 +227,29 @@ h5u_set_num_elements (
 	stride[0] = 1;
 	start[0] = 0;
 	for (i=0; i<f->myproc; i++) {
-		start[0] += f->pnparticles[i];
+		start[0] += f->u->pnparticles[i];
 	}
 	
         /* compute total nparticles */
 	total = 0;
 	for (i=0; i < f->nprocs; i++) {
-		total += f->pnparticles[i];
+		total += f->u->pnparticles[i];
 	}
 
 	/* declare overall datasize */
-	TRY ( f->shape = _h5_create_space ( f, 1, &total, &total ) );
+	TRY ( f->u->shape = _h5_create_space ( f, 1, &total, &total ) );
 
 	/* declare overall data size  but then will select a subset */
-	TRY ( f->diskshape = _h5_create_space ( f, 1, &total, &total) );
+	TRY ( f->u->diskshape = _h5_create_space ( f, 1, &total, &total) );
 
 	/* declare local memory datasize */
-	TRY ( f->memshape = _h5_create_space (
-		      f, 1, &(f->nparticles), &dmax ) );
+	TRY ( f->u->memshape = _h5_create_space (
+		      f, 1, &(f->u->nparticles), &dmax ) );
 
 	count[0] = nparticles;
 	TRY ( _h5_select_hyperslab_of_space (
 		      f,
-		      f->diskshape,
+		      f->u->diskshape,
 		      H5S_SELECT_SET,
 		      start, stride, count,
 		      NULL ) );
@@ -330,8 +332,7 @@ h5u_set_view (
 	  For now, we interpret start=-1 to mean 0 and 
 	  end==-1 to mean end of file
 	*/
-	total = (hsize_t) h5u_get_num_elems ( f );
-	if ( total < 0 ) return HANDLE_H5_GET_NUM_PARTICLES_ERR ( f, total );
+	TRY ( total = (hsize_t) h5u_get_num_elems ( f ) );
 
 	if ( start == -1 ) start = 0;
 	if ( end == -1 )   end = total;
@@ -415,32 +416,31 @@ h5u_set_canonical_view (
 	h5_int64_t n = 0;
 	int i = 0;
 	
-	n = h5u_get_num_elems ( f );
-	if ( n < 0 ) return HANDLE_H5_GET_NUM_PARTICLES_ERR ( n );
+	TRY ( n = h5u_get_num_elems ( f ) );
 	/* 
 	   now lets query the attributes for this group to see if there
 	   is a 'pnparticles' group that contains the offsets for the
 	   processors.
 	*/
 	if ( h5_read_attrib (
+		     f,
 		     f->step_gid,
-		     "pnparticles", f->pnparticles ) < 0) {
+		     "pnparticles", f->u->pnparticles ) < 0) {
 		/*
 		  Attribute "pnparticles" is not available.  So
 		  subdivide the view into NP mostly equal pieces
 		*/
 		n /= f->nprocs;
 		for ( i=0; i<f->nprocs; i++ ) {
-			f->pnparticles[i] = n;
+			f->u->pnparticles[i] = n;
 		}
 	}
 
 	for ( i = 0; i < f->myproc; i++ ){
-		start += f->pnparticles[i];
+		start += f->u->pnparticles[i];
 	}
-	end = start + f->pnparticles[f->myproc] - 1;
-	herr = h5u_set_view ( f, start, end );
-	if ( herr < 0 ) return HANDLE_H5_SET_VIEW_ERR ( herr, start, end );
+	end = start + f->u->pnparticles[f->myproc] - 1;
+	TRY ( h5u_set_view ( f, start, end ) );
 
 #endif