- introducing types for local and globale IDs and indices

- different structurs for local elements and stored elements in file
- interator for border facets added

src/h5core/h5_errorhandling.c

@@ -37,8 +37,10 @@ h5_err_t
 h5_set_debuglevel (
 	const h5_id_t level	/*!< debug level */
 	) {
-	if (level < 0 || level > 5) return H5_ERR_INVAL;
-	h5priv_debug_level = level;
+	if (level < 0) 
+		h5priv_debug_level = 0;
+	else
+		h5priv_debug_level = level;
 	return H5_SUCCESS;
 }
 
@@ -168,15 +170,9 @@ h5priv_vprintf (
 	const char* fmt,
 	va_list ap
 	) {
-	size_t size = strlen (prefix);	// to avoid remark #981 with Intel CC
-	size += strlen (fmt);
-	size += strlen (__funcname) + 16;
-
-	char *fmt2 = (char*)malloc (size);
-	if (fmt2 == NULL) return;
-	sprintf (fmt2, "%s: %s: %s\n", prefix, __funcname, fmt); 
+	char fmt2[2048];
+	snprintf (fmt2, sizeof(fmt2), "%s: %s: %s\n", prefix, __funcname, fmt); 
 	vfprintf (f, fmt2, ap);
-	free (fmt2);
 }
 
 /*!
@@ -219,40 +215,6 @@ h5_verror (
 	h5priv_vprintf (stderr, "E", f->__funcname, fmt, ap);
 }
 
-
-/*!
-  \ingroup h5_core_errorhandling
-
-  Print a warning message to \c stderr.
-*/
-void
-h5_vwarn (
-	const h5_file_t* const f,
-	const char* fmt,
-	va_list ap
-	) {
-	if (h5priv_debug_level < 2) return;
-	h5priv_vprintf (stderr, "W", f->__funcname, fmt, ap);
-}
-
-/*!
-  \ingroup h5_core_errorhandling
-
-  Print a warning message to \c stderr.
-*/
-void
-h5_warn (
-	const h5_file_t* const f,
-	const char* fmt,
-	...
-	) {
-	va_list ap;
-	va_start (ap, fmt);
-	h5_vwarn (f, fmt, ap);
-	va_end (ap);
-}
-
-
 /*!
   \ingroup h5_core_errorhandling
 

src/h5core/h5_hdf5.c

@@ -97,13 +97,13 @@ h5priv_close_hdf5_group (
 	h5_file_t* const f,
 	const hid_t group_id
 	) {
+	if (group_id == 0 || group_id == -1) return H5_SUCCESS;
+
 	h5_debug (f, "%s (group_id=%lld, group_name=\"%s\")", 
 		  __func__,
 		  (long long)group_id,
 		  h5_get_objname (group_id));
 
-	h5_debug (f, "%s (group_id=%lld)", __func__, (long long)group_id);
-	if (group_id <= 0) return H5_SUCCESS; 
 	if (H5Gclose (group_id) < 0 ) {
 		return h5_error (
 			f,
@@ -1184,7 +1184,21 @@ h5priv_hdf5_link_exists (
 	const hid_t loc_id,
 	const char* name
 	) {
+	/* Save old error handler */
+	H5E_auto2_t old_func;
+	void *old_client_data;
+
+	H5Eget_auto2(H5E_DEFAULT, &old_func, &old_client_data);
+
+	/* Turn off error handling */
+	H5Eset_auto(H5E_DEFAULT, NULL, NULL);
+
+	/* Probe. Likely to fail, but that’s okay */
 	htri_t exists = H5Lexists ( loc_id, name, H5P_DEFAULT );
+
+	/* Restore previous error handler */
+	H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);
+
 	if (exists < 0 )
 		return h5_error (f,
 			H5_ERR_HDF5,
@@ -1484,3 +1498,15 @@ h5_get_objname (
 
 	return objname;
 }
+
+#if 0
+herr_t
+h5priv_get_objinfo_by_name (
+	hid_t loc_id,
+	const char *object_name,
+	H5O_info_t *object_info
+	) {
+	herr_t herr = H5Oget_info_by_name (
+		loc_id, object_name, object_info, H5P_DEFAULT); 
+}
+#endif

src/h5core/h5_openclose.c

@@ -295,6 +295,7 @@ h5upriv_close_file (
 	h5_file_t* const f	/*!< file handle */
 	) {
 	struct h5u_fdata* u = f->u;
+	h5_debug (f, "%s ()", __func__);
 
 	f->__errno = H5_SUCCESS;
 	TRY( h5priv_close_hdf5_dataspace (f, u->shape) );
@@ -322,7 +323,7 @@ h5bpriv_close_file (
 	h5_file_t* const f	/*!< IN: file handle */
 	) {
 	struct h5b_fdata* b = f->b;
-
+	h5_debug (f, "%s ()", __func__);
 	TRY( h5priv_close_hdf5_group (f, b->block_gid) );
 	TRY( h5priv_close_hdf5_group (f, b->field_gid) );
 	TRY( h5priv_close_hdf5_dataspace (f, b->shape) );
@@ -362,7 +363,6 @@ h5_close_file (
 #ifndef PARALLEL_IO
 	TRY( h5tpriv_close_file (f) );
 #endif
-	TRY( h5priv_close_hdf5_group (f, f->step_gid) );
 	TRY( h5priv_close_hdf5_property (f, f->xfer_prop) );
 	TRY( h5priv_close_hdf5_property (f, f->access_prop) );
 	TRY( h5priv_close_hdf5_property (f, f->create_prop) );

src/h5core/h5_readwrite.c

@@ -24,14 +24,18 @@ h5priv_write_dataset_by_name (
 	h5_info (f, "Writing dataset %s/%s.",
 		 h5_get_objname (loc_id), dsinfo->name);
 
+#if 0
 	H5O_info_t obj_info;
 	herr_t herr = H5Oget_info_by_name(
 		loc_id,
 		dsinfo->name,
 		&obj_info,
 		H5P_DEFAULT);
-
-	if ((herr >= 0) && ((f->mode==H5_O_WRONLY) || (f->mode==H5_O_APPEND))) {
+#else
+	h5_err_t exists;
+	TRY( exists = h5priv_hdf5_link_exists (f, loc_id, dsinfo->name) );
+#endif
+	if ((exists > 0) && ((f->mode==H5_O_WRONLY) || (f->mode==H5_O_APPEND))) {
 		h5_warn (f,
 			 "Dataset %s/%s already exist.",
 			 h5_get_objname (loc_id), dsinfo->name);
@@ -46,7 +50,7 @@ h5priv_write_dataset_by_name (
 	hid_t diskspace_id;
 	hid_t memspace_id;
 
-	if (herr >= 0) {
+	if (exists) {
 		/* overwrite dataset */
 		TRY( dset_id = h5priv_open_hdf5_dataset (
 			     f, 
@@ -158,8 +162,8 @@ h5_err_t
 h5priv_close_step (
 	h5_file_t* const f
 	) {
-
-	if (f->step_gid < 0) return H5_SUCCESS;
+	h5_debug (f, "%s ()", __func__);
+	if (f->step_gid <= 0) return H5_SUCCESS;
 	TRY( h5tpriv_close_step (f) );
 	TRY( h5priv_close_hdf5_group (f, f->step_gid) );
 

src/h5core/h5_syscall.c

@@ -10,12 +10,13 @@ h5priv_alloc (
 	void* ptr,
 	const size_t size
 	) {
-	h5_debug (f, "Allocating %lu bytes.", size); 
+	h5_debug (f, "%s (ptr=%p, size=%lu)", __func__, ptr, size); 
 	ptr = realloc (ptr, size);
 	if (ptr == NULL) {
 		h5_error (f, H5_ERR_NOMEM, "Out of memory.");
 		return (void*)(H5_ERR);
 	}
+	h5_debug (f, "%s (): return address: 0x%p", __func__, ptr);
 	return ptr;
 }
 
@@ -25,12 +26,13 @@ h5priv_calloc (
 	const size_t count,
 	const size_t size
 	) {
-	h5_debug (f, "Allocating %lu * %lu bytes.", count, size); 
+	h5_debug (f, "%s (count=%lu , size=%lu)", __func__, count, size);
 	void* ptr = calloc (count, size);
 	if (ptr == NULL) {
 		h5_error (f, H5_ERR_NOMEM, "Out of memory.");
 		return (void*)(H5_ERR);
 	}
+	h5_debug (f, "%s (): return address: 0x%p", __func__, ptr);
 	return ptr;
 }
 
@@ -39,8 +41,10 @@ h5priv_free (
 	h5_file_t* const f,
 	void* ptr
 	) {
-#pragma unused f
-	if (ptr) free (ptr);
+	if (ptr) {
+		h5_debug (f, "%s (%p)", __func__,  ptr); 
+		free (ptr);
+	}
 	return H5_SUCCESS;
 }
 

src/h5core/h5t_access_private.h

@@ -48,7 +48,7 @@ h5tpriv_set_loc_elem_child_idx (
 		f, elem_idx, child_idx);
 }
 
-static inline h5_id_t
+static inline h5t_lvl_idx_t
 h5tpriv_get_loc_elem_level_idx (
 	h5_file_t* const f,
 	h5_loc_idx_t elem_idx
@@ -57,11 +57,11 @@ h5tpriv_get_loc_elem_level_idx (
 		f, elem_idx);
 }
 
-static inline h5_id_t
+static inline h5t_lvl_idx_t
 h5tpriv_set_loc_elem_level_idx (
 	h5_file_t* const f,
 	h5_loc_idx_t elem_idx,
-	h5_id_t lvl_idx
+	h5t_lvl_idx_t lvl_idx
 	) {
 	return (*f->t->methods.access->set_loc_elem_level_idx)(
 		f, elem_idx, lvl_idx);
@@ -253,4 +253,47 @@ h5tpriv_set_glb_elem_neighbor_idx (
 		f, elem_idx, face_idx, idx);
 }
 
+static inline h5_err_t
+h5tpriv_set_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	return (*f->t->methods.access->set_boundary_elem_flag)(f, elem_idx);
+}
+
+static inline h5_err_t
+h5tpriv_clear_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	return (*f->t->methods.access->clear_boundary_elem_flag)(f, elem_idx);
+}
+
+static inline int
+h5tpriv_is_boundary_elem (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	return (*f->t->methods.access->is_boundary_elem)(f, elem_idx);
+}
+
+static inline int
+h5tpriv_is_boundary_facet (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx,
+	h5_loc_idx_t facet_idx
+	) {
+	return (*f->t->methods.access->is_boundary_facet)(f, elem_idx, facet_idx);
+}
+
+static inline int
+h5tpriv_is_boundary_face (
+	h5_file_t* const f,
+	const int dim,
+	const h5_loc_idx_t elem_idx,
+	const h5_loc_idx_t facet_idx
+	) {
+	return (*f->t->methods.access->is_boundary_face)(f, dim, elem_idx, facet_idx);
+}
+
 #endif

src/h5core/h5t_access_tetm.c

@@ -46,21 +46,21 @@ set_loc_elem_child_idx (
 	return child_idx;
 }
 
-static h5_id_t
+static h5t_lvl_idx_t
 get_loc_elem_level_idx (
 	h5_file_t* const f,
 	const h5_loc_idx_t elem_idx
 	) {
-	return f->t->loc_elems.tets[elem_idx].idx;
+	return f->t->loc_elems.tets[elem_idx].level_idx;
 }
 
-static h5_id_t
+static h5t_lvl_idx_t
 set_loc_elem_level_idx (
 	h5_file_t* const f,
 	const h5_loc_idx_t elem_idx,
-	const h5_id_t level_idx
+	const h5t_lvl_idx_t level_idx
 	) {
-	f->t->loc_elems.tets[elem_idx].idx = level_idx;
+	f->t->loc_elems.tets[elem_idx].level_idx = level_idx;
 	return level_idx;
 }
 
@@ -240,6 +240,56 @@ set_glb_elem_neighbor_idx (
 	return neighbor_idx;
 }
 
+static h5_err_t
+set_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	f->t->loc_elems.tets[elem_idx].flags |= H5T_BOUNDARY_ELEM_FLAG;
+	return H5_SUCCESS;
+}
+
+static h5_err_t
+clear_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	f->t->loc_elems.tets[elem_idx].flags &= ~H5T_BOUNDARY_ELEM_FLAG;
+	return H5_SUCCESS;
+}
+
+
+static int
+is_boundary_elem  (
+	h5_file_t* const f,
+	const h5_loc_idx_t elem_idx
+	) {
+	return (f->t->loc_elems.tets[elem_idx].flags & H5T_BOUNDARY_ELEM_FLAG) ? 1 : 0;
+}
+
+static int
+is_boundary_facet (
+	h5_file_t* const f,
+	const h5_loc_idx_t elem_idx,
+	const h5_loc_idx_t facet_idx
+	) {
+	return (f->t->loc_elems.tets[elem_idx].neighbor_indices[facet_idx] == -1);
+}
+
+static int
+is_boundary_face (
+	h5_file_t* const f,
+	const int dim,
+	const h5_loc_idx_t elem_idx,
+	const h5_loc_idx_t facet_idx
+	) {
+#pragma unused f
+#pragma unused dim
+#pragma unused elem_idx
+#pragma unused facet_idx
+	return h5_error_internal (f, __FILE__, __func__, __LINE__);
+}
+
 struct h5t_access_methods h5tpriv_access_tetm_methods = {
 	get_loc_elem,
 	get_loc_elem_parent_idx,
@@ -267,5 +317,10 @@ struct h5t_access_methods h5tpriv_access_tetm_methods = {
 	get_glb_elem_neighbor_indices,
 	get_glb_elem_neighbor_idx,
 	set_glb_elem_neighbor_idx,
+	set_boundary_elem_flag,
+	clear_boundary_elem_flag,
+	is_boundary_elem,
+	is_boundary_facet,
+	is_boundary_face,
 };
 	

src/h5core/h5t_access_trim.c

@@ -46,21 +46,21 @@ set_loc_elem_child_idx (
 	return child_idx;
 }
 
-static h5_id_t
+static h5t_lvl_idx_t
 get_loc_elem_level_idx (
 	h5_file_t* const f,
 	const h5_loc_idx_t elem_idx
 	) {
-	return f->t->loc_elems.tris[elem_idx].idx;
+	return f->t->loc_elems.tris[elem_idx].level_idx;
 }
 
-static h5_id_t
+static h5t_lvl_idx_t
 set_loc_elem_level_idx (
 	h5_file_t* const f,
 	const h5_loc_idx_t elem_idx,
-	const h5_loc_idx_t level_idx
+	const h5t_lvl_idx_t level_idx
 	) {
-	f->t->loc_elems.tris[elem_idx].idx = level_idx;
+	f->t->loc_elems.tris[elem_idx].level_idx = level_idx;
 	return level_idx;
 }
 
@@ -240,6 +240,57 @@ set_glb_elem_neighbor_idx (
 	return neighbor_idx;
 }
 
+
+static h5_err_t
+set_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	f->t->loc_elems.tris[elem_idx].flags |= H5T_BOUNDARY_ELEM_FLAG;
+	return H5_SUCCESS;
+}
+
+static h5_err_t
+clear_boundary_elem_flag (
+	h5_file_t* const f,
+	h5_loc_idx_t elem_idx
+	) {
+	f->t->loc_elems.tris[elem_idx].flags &= ~H5T_BOUNDARY_ELEM_FLAG;
+	return H5_SUCCESS;
+}
+
+
+static int
+is_boundary_elem  (
+	h5_file_t* const f,
+	const h5_loc_idx_t elem_idx
+	) {
+	return (f->t->loc_elems.tris[elem_idx].flags & H5T_BOUNDARY_ELEM_FLAG) ? 1 : 0;
+}
+
+static int
+is_boundary_facet (
+	h5_file_t* const f,
+	const h5_loc_idx_t elem_idx,
+	const h5_loc_idx_t facet_idx
+	) {
+	return (f->t->loc_elems.tris[elem_idx].neighbor_indices[facet_idx] == -1);
+}
+
+static int
+is_boundary_face (
+	h5_file_t* const f,
+	const int dim,
+	const h5_loc_idx_t elem_idx,
+	const h5_loc_idx_t facet_idx
+	) {
+#pragma unused f
+#pragma unused dim
+#pragma unused elem_idx
+#pragma unused facet_idx
+	return h5_error_internal (f, __FILE__, __func__, __LINE__);
+}
+
 struct h5t_access_methods h5tpriv_access_trim_methods = {
 	get_loc_elem,
 	get_loc_elem_parent_idx,
@@ -267,5 +318,10 @@ struct h5t_access_methods h5tpriv_access_trim_methods = {
 	get_glb_elem_neighbor_indices,
 	get_glb_elem_neighbor_idx,
 	set_glb_elem_neighbor_idx,
+	set_boundary_elem_flag,
+	clear_boundary_elem_flag,
+	is_boundary_elem,
+	is_boundary_facet,
+	is_boundary_face,
 };
 	

src/h5core/h5t_adjacencies_private.h

@@ -31,9 +31,9 @@ h5tpriv_release_adjacency_structs (
 static inline h5_err_t
 h5tpriv_update_adjacency_structs (
 	h5_file_t* const f,
-	const h5_id_t level_id
+	const h5t_lvl_idx_t level_id
 	) {
-	h5_debug (f, "%s (%lld)", __func__, level_id);
+	h5_debug (f, "%s (%lld)", __func__, (long long)level_id);
 	return (*f->t->methods.adjacency->update_internal_structs)(f, level_id);
 }
 

src/h5core/h5t_adjacencies_tetm.c

@@ -56,7 +56,7 @@ release_tv (
 static inline h5_err_t
 compute_elems_of_vertices (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
 	/* expand structure */
 	TRY( alloc_tv (f) );
@@ -95,7 +95,7 @@ release_te (
 static inline h5_err_t
 compute_elems_of_edges (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
 	h5t_fdata_t *t = f->t;
 	h5_loc_idx_t elem_idx = (from_lvl <= 0) ? 0 : t->num_elems[from_lvl-1];
@@ -125,7 +125,7 @@ release_td (
 static inline h5_err_t
 compute_elems_of_triangles (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
 	h5t_fdata_t* t = f->t;
 	h5_loc_idx_t elem_idx = (from_lvl <= 0) ? 0 : t->num_elems[from_lvl-1];
@@ -266,7 +266,7 @@ compute_direct_children_of_triangle (
 		{{2,0},{2,2},{2,3},{1,7}},
 		{{3,1},{3,2},{3,3},{3,6}}
 	};
-	int num_faces = h5tpriv_ref_elem_get_num_triangles (f->t);
+	int num_faces = h5tpriv_ref_elem_get_num_facets (f->t);
 	if ((face_idx < 0) || (face_idx >= num_faces)) {
 		return h5_error_internal (f, __FILE__, __func__, __LINE__); 
 	}
@@ -456,11 +456,11 @@ get_triangles_uadj_to_vertex (
 			continue;
 		}
 		h5_loc_idx_t facet_idx;
-		facet_idx = h5tpriv_get_facet_connected_to_vertex (t->ref_elem, face_idx, 0);
+		facet_idx = h5tpriv_get_triangles_connected_to_vertex (t->ref_elem, face_idx, 0);
 		TRY( add_triangle (f, *list, facet_idx, elem_idx) );
-		facet_idx = h5tpriv_get_facet_connected_to_vertex (t->ref_elem, face_idx, 1);
+		facet_idx = h5tpriv_get_triangles_connected_to_vertex (t->ref_elem, face_idx, 1);
 		TRY( add_triangle (f, *list, facet_idx, elem_idx) );
-		facet_idx = h5tpriv_get_facet_connected_to_vertex (t->ref_elem, face_idx, 2);
+		facet_idx = h5tpriv_get_triangles_connected_to_vertex (t->ref_elem, face_idx, 2);
 		TRY( add_triangle (f, *list, facet_idx, elem_idx) );
 	}
 	return H5_SUCCESS;
@@ -842,7 +842,7 @@ get_adjacencies (
 static h5_err_t
 update_internal_structs (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
 	clock_t t1 = clock();
 	TRY( compute_elems_of_vertices (f, from_lvl) );

src/h5core/h5t_adjacencies_trim.c

@@ -55,9 +55,9 @@ release_tv (
 static inline h5_err_t
 compute_elems_of_vertices (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
-	h5_debug (f, "%s (%lld)", __func__, from_lvl);
+	h5_debug (f, "%s (%lld)", __func__, (long long)from_lvl);
 	/* expand structure */
 	TRY( alloc_tv (f) );
 
@@ -78,7 +78,7 @@ compute_elems_of_vertices (
 					     face_idx, idx)) );
 		}
 	}
-	h5_debug (f, "%s (%lld): done", __func__, from_lvl);
+	h5_debug (f, "%s (%lld): done", __func__, (long long)from_lvl);
 	return H5_SUCCESS;
 }
 
@@ -97,9 +97,9 @@ release_te (
 static inline h5_err_t
 compute_elems_of_edges (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
-	h5_debug (f, "%s (%lld)", __func__, from_lvl);
+	h5_debug (f, "%s (%lld)", __func__, (long long)from_lvl);
 	h5t_fdata_t *t = f->t;
 	h5_loc_idx_t elem_idx = (from_lvl <= 0) ? 0 : t->num_elems[from_lvl-1];
 	h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
@@ -113,7 +113,7 @@ compute_elems_of_edges (
 				      f, face_idx, elem_idx, &retval ) );
 		}
 	}
-	h5_debug (f, "%s (%lld): done", __func__, from_lvl);
+	h5_debug (f, "%s (): done", __func__);
 	return H5_SUCCESS;
 }
 
@@ -474,7 +474,7 @@ get_adjacencies (
 static inline h5_err_t
 update_internal_structs (
 	h5_file_t* const f,
-	const h5_id_t from_lvl
+	const h5t_lvl_idx_t from_lvl
 	) {
 	h5_debug (f, "%s (%lld)", __func__, (long long)from_lvl);
 	clock_t t1 = clock();

src/h5core/h5t_core_private.h

@@ -20,4 +20,76 @@
 
 #include "h5t_errorhandling_private.h"
 
+
+#if !(defined(ULLONG_MAX))
+#define ULLONG_MAX  (LLONG_MAX * 2ULL + 1ULL)
+#endif
+
+/*
+ ID's: 64bit
+
+ Vertices:
+   000100vv tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt
+   3V TT TT TT TT TT TT
+
+ Edges:
+   00100eee tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt
+   2E TT TT TT TT TT TT
+
+ Trinagles:
+   001100dd tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt
+   1D TT TT TT TT TT TT
+
+ Tets:
+   01000000 tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt tttttttt
+   00 TT TT TT TT TT TT
+
+*/
+#define H5T_ELEM_MASK		( (h5_id_t) (ULLONG_MAX >> 8) )
+#define H5T_ETYPE_MASK		( 7ull << (sizeof(h5_id_t)*8-4) )
+#define H5T_COMPONENT_MASK	(~H5T_ELEM_MASK)
+#define H5T_COMPONENT_ID_MASK	(15ull << (sizeof(h5_id_t)*7) )
+
+#define H5T_ETYPE_VERTEX	(1)
+#define H5T_ETYPE_EDGE		(2)
+#define H5T_ETYPE_TRIANGLE	(3)
+#define H5T_ETYPE_TET		(4)
+
+#define h5tpriv_set_entity_type( type, elem_idx )		\
+	(							\
+		((h5_id_t)(type) << (sizeof(h5_id_t)*8-4)) |	\
+		((h5_id_t)(elem_idx))				\
+		)
+#define h5tpriv_get_entity_type( entity_id )			\
+	((entity_id & H5T_ETYPE_MASK) >> (sizeof(h5_id_t)*8-4))
+
+
+#define h5tpriv_build_id( type, face_idx, elem_idx )		\
+	(							\
+		((h5_id_t)(type) << (sizeof(h5_id_t)*8-4)) |	\
+		((h5_id_t)(face_idx) << (sizeof(h5_id_t)*7)) |	\
+		((h5_id_t)(elem_idx) & H5T_ELEM_MASK)		\
+		)
+
+#define h5tpriv_build_vertex_id( face_idx, elem_idx )			\
+	( h5tpriv_build_id ( H5T_ETYPE_VERTEX, face_idx, elem_idx ) )
+
+#define h5tpriv_build_edge_id( face_idx, elem_idx )			\
+	( h5tpriv_build_id ( H5T_ETYPE_EDGE, face_idx, elem_idx ) )
+
+#define h5tpriv_build_triangle_id( face_idx, elem_idx )	 \
+	( h5tpriv_build_id ( H5T_ETYPE_TRIANGLE, face_idx, elem_idx ) )
+
+#define h5tpriv_build_elem_id( elem_idx )				\
+	( h5tpriv_build_id ( f->t->mesh_type, 0, elem_idx ) )
+
+#define h5tpriv_get_face_idx( entity_id )				\
+	( (entity_id & H5T_COMPONENT_ID_MASK) >> (sizeof(h5_id_t)*7) )
+
+#define h5tpriv_get_elem_idx( entity_id )		\
+	( entity_id & H5T_ELEM_MASK )
+
+
+#define H5T_BOUNDARY_ELEM_FLAG 1
+
 #endif

src/h5core/h5t_inquiry.c

@@ -60,7 +60,7 @@ h5t_get_num_levels (
 
   \return	Current level ID.
 */
-h5_id_t
+h5t_lvl_idx_t
 h5t_get_level (
 	h5_file_t* const f
 	) {

src/h5core/h5t_map.c

@@ -27,197 +27,6 @@ cmp_vertices (
 	return 0;
 }
 
-static int
-qsort_cmp_vertices (
-	void* _f,
-	const void* _vertex_idx1,
-	const void* _vertex_idx2
-	) {
-	h5_file_t* f = (h5_file_t*)_f;
-	h5_loc_idx_t vertex_idx1 = *(h5_loc_idx_t*)_vertex_idx1;
-	h5_loc_idx_t vertex_idx2 = *(h5_loc_idx_t*)_vertex_idx2;
-
-	return cmp_vertices (
-		f->t->vertices[vertex_idx1].P,
-		f->t->vertices[vertex_idx2].P );
-}
-
-/*!
-  Sort vertices. Store local id's in a sorted array for binary search. 
-*/
-h5_err_t
-h5tpriv_sort_vertices (
-	h5_file_t* const f
-	) {
-	h5t_fdata_t* t = f->t;
-	if (t->num_levels <= 0) return H5_SUCCESS;
-
-	h5_loc_idx_t vertex_idx = t->cur_level > 0 ?
-		t->num_vertices[t->cur_level-1] : 0;
-
-	h5_loc_idx_t num_vertices = t->num_vertices[t->num_levels-1];
-	for (; vertex_idx < num_vertices; vertex_idx++) {
-		t->sorted_lvertices.items[vertex_idx] = vertex_idx;
-	}
-	t->sorted_lvertices.num_items = num_vertices;
-
-	h5priv_qsort_r (
-		t->sorted_lvertices.items,
-		num_vertices,
-		sizeof (t->sorted_lvertices.items[0]),
-		f,
-		qsort_cmp_vertices);
-
-	return H5_SUCCESS;
-}
-
-
-/*!
-  Return local vertex index of a vertex given by its coordinates.
-
-  \return	local vertex idx if found
-  \return	else negativ value
- */
-h5_loc_idx_t
-h5tpriv_get_local_vid (
-	h5_file_t* const f,
-	h5_float64_t P[3]
-	) {
-	h5t_fdata_t*t = f->t;
-	register h5_loc_idx_t low = 0;
-	register h5_loc_idx_t high = t->sorted_lvertices.num_items - 1;
-	while (low <= high) {
-		register h5_loc_idx_t mid = (low + high) / 2;
-
-		h5_loc_idx_t vertex_idx = t->sorted_lvertices.items[mid];
-		h5_float64_t *P1 = t->vertices[vertex_idx].P;
-	        int diff = cmp_vertices ( P, P1 );
-           	if ( diff < 0 )
-               		high = mid - 1;
-           	else if ( diff > 0 )
-               		low = mid + 1;
-           	else
-               		return t->sorted_lvertices.items[mid]; // found
-       	}
-       	return -(low+1);  // not found
-}
-
-/*!
-  compare two elems given by their local id
-*/
-static inline int
-cmp_elems (
-	h5_file_t* const  f,
-	const h5_loc_idx_t elem_idx1,
-	const h5_loc_idx_t elem_idx2
-	) {
-	h5t_fdata_t* t = f->t;
-	int num_vertices = h5tpriv_ref_elem_get_num_vertices (t);
-	h5_loc_idx_t* indices1 = h5tpriv_get_loc_elem_vertex_indices (f, elem_idx1);
-	h5_loc_idx_t* indices2 = h5tpriv_get_loc_elem_vertex_indices (f, elem_idx2);
-
-	int i;
-	for (i = 0; i < num_vertices; i++) {
-		h5_float64_t* v1 = t->vertices[indices1[i]].P;
-		h5_float64_t* v2 = t->vertices[indices2[i]].P;
-		int r = cmp_vertices (v1, v2);
-		if (r < 0)
-			return -1;
-		else if (r > 0)
-			return 1;
-	}
-	return 0;
-}
-
-static inline int
-cmp_elems1 (
-	h5_file_t* f,
-	h5_loc_idx_t elem_idx1,
-	h5_loc_idx_t elem_idx2
-	) {
-	h5t_fdata_t *t = f->t;
-	int num_vertices = h5tpriv_ref_elem_get_num_vertices (t);
-	h5_loc_idx_t* indices1 = h5tpriv_get_loc_elem_vertex_indices (f, elem_idx1);
-	h5_loc_idx_t* indices2 = h5tpriv_get_loc_elem_vertex_indices (f, elem_idx2);
-
-	int imap[] = { 1, 0, 2, 3 };
-	int i;
-	for ( i = 0; i < num_vertices; i++ ) {
-		h5_float64_t* v1 = t->vertices[ indices1[imap[i]] ].P;
-		h5_float64_t* v2 = t->vertices[ indices2[imap[i]] ].P;
-		int r = cmp_vertices (v1, v2);
-		if (r < 0)
-			return -1;
-		else if (r > 0)
-			return 1;
-	}
-	return 0;
-}
-
-
-static int
-qsort_cmp_elems0 (
-	void* _f,
-	const void* _elem_idx1,
-	const void* _elem_idx2
-	) {
-	h5_file_t* f = (h5_file_t*)_f;
-	h5_loc_idx_t elem_idx1 = *(h5_loc_idx_t*)_elem_idx1;
-	h5_loc_idx_t elem_idx2 = *(h5_loc_idx_t*)_elem_idx2;
-	return cmp_elems (f, elem_idx1, elem_idx2);
-}
-
-static int
-qsort_cmp_elems1 (
-	void* _f,
-	const void* _elem_idx1,
-	const void* _elem_idx2
-	) {
-	h5_file_t* f = (h5_file_t*)_f;
-	h5_loc_idx_t elem_idx1 = *(h5_loc_idx_t*)_elem_idx1;
-	h5_loc_idx_t elem_idx2 = *(h5_loc_idx_t*)_elem_idx2;
-	return cmp_elems1 (f, elem_idx1, elem_idx2);
-}
-
-/*!
-  Sort elements. Store local id's in a sorted array so we can run a
-  binary search. 
-*/
-h5_err_t
-h5tpriv_sort_loc_elems (
-	h5_file_t* const f
-	) {
-	h5t_fdata_t* t = f->t;
-	if (t->num_levels <= 0) return H5_SUCCESS;
-	h5_loc_idx_t elem_idx = t->cur_level > 0 ? t->num_elems[t->cur_level-1] : 0;
-	h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
-
-	int k;
-	h5_loc_idx_t i;
-	for (k = 0; k < 2; k++) {
-		TRY( h5priv_alloc_idlist_items (f, &t->sorted_elems[k], num_elems) );
-		for (i = elem_idx; i < num_elems; i++) {
-			t->sorted_elems[k].items[i] = i;
-		}
-		t->sorted_elems[k].num_items = num_elems;
-	}
-
-	h5priv_qsort_r (
-		t->sorted_elems[0].items,
-		num_elems,
-		sizeof (t->sorted_elems[0].items[0]),
-		f,    
-		qsort_cmp_elems0);
-	h5priv_qsort_r (
-		t->sorted_elems[1].items,
-		num_elems,
-		sizeof (t->sorted_elems[1].items[0]),
-		f,
-		qsort_cmp_elems1);
-
-	return H5_SUCCESS;
-}
-
 /*!
   Sort (small) array of local vertex indices geometrically. 
  */
@@ -385,7 +194,19 @@ h5t_get_vertex_indices_of_entity (
 	h5_loc_id_t entity_type = h5tpriv_get_entity_type (entity_id);
 	h5_loc_idx_t face_idx = h5tpriv_get_face_idx (entity_id);
 	h5_loc_idx_t elem_idx = h5tpriv_get_elem_idx (entity_id);
-	int dim =  map_entity_type_to_dimension[entity_type];
+	int dim = map_entity_type_to_dimension[entity_type];
+
+	return h5t_get_vertex_indices_of_entity2 (f, dim, face_idx, elem_idx, vertex_indices);
+}
+
+h5_err_t
+h5t_get_vertex_indices_of_entity2 (
+	h5_file_t* const f,		// [in]
+	const int dim,			// [in] dimension
+	const h5_loc_idx_t face_idx,	// [in] vertex index according ref. element
+	const h5_loc_idx_t elem_idx,	// [in] local element index
+	h5_loc_idx_t* vertex_indices   	// [out]
+	) {
 	h5_loc_idx_t* indices = h5tpriv_get_loc_elem_vertex_indices (f, elem_idx);
 	const h5t_ref_elem_t* ref_elem = f->t->ref_elem;
 	int num_vertices = ref_elem->num_vertices_of_face[dim][face_idx];
@@ -436,8 +257,8 @@ h5t_get_vertex_indices_of_edge (
 }
 
 /*!
-  Get local vertex ID's of an edge. The edge is specified by the local 
-  element index and the sub-entity ID of the edge according the reference
+  Get local vertex indices of an edge. The edge is specified by the local 
+  element index and the face number of the edge according the reference
   element.
 
   This function can be used with tetrahedral and triangle meshes.

src/h5core/h5t_map_private.h

@@ -8,22 +8,12 @@ h5tpriv_sort_local_vertex_indices (
 	const h5_size_t size
 	);
 
-h5_err_t
-h5tpriv_sort_vertices (
-	h5_file_t * const f
-	);
-
 h5_loc_idx_t
 h5tpriv_get_local_vid (
 	h5_file_t * const f,
 	h5_float64_t P[3]
 	);
 
-h5_err_t
-h5tpriv_sort_loc_elems (
-	h5_file_t * const f
-	);
-
 h5_err_t
 h5tpriv_rebuild_vertex_indices_mapping (
 	h5_file_t * const f

src/h5core/h5t_openclose.c

@@ -265,7 +265,7 @@ init_fdata (
 	t->dsinfo_num_vertices.dims[0] = 0;
 	t->dsinfo_num_vertices.max_dims[0] = H5S_UNLIMITED;
 	t->dsinfo_num_vertices.chunk_dims[0] = 4096;
-	t->dsinfo_num_vertices.type_id = t->dtypes.h5_id_t;
+	t->dsinfo_num_vertices.type_id = t->dtypes.h5_glb_idx_t;
 	TRY( t->dsinfo_num_vertices.create_prop = h5priv_create_hdf5_property (
 		     f,
 		     H5P_DATASET_CREATE) );
@@ -298,7 +298,7 @@ init_fdata (
 	t->dsinfo_num_elems.dims[0] = 0;
 	t->dsinfo_num_elems.max_dims[0] = H5S_UNLIMITED;
 	t->dsinfo_num_elems.chunk_dims[0] = 4096;
-	t->dsinfo_num_elems.type_id = t->dtypes.h5_id_t;
+	t->dsinfo_num_elems.type_id = t->dtypes.h5_glb_idx_t;
 	TRY( t->dsinfo_num_elems.create_prop = h5priv_create_hdf5_property (
 		     f,
 		     H5P_DATASET_CREATE) );
@@ -315,7 +315,7 @@ init_fdata (
 	t->dsinfo_num_elems_on_level.dims[0] = 0;
 	t->dsinfo_num_elems_on_level.max_dims[0] = H5S_UNLIMITED;
 	t->dsinfo_num_elems_on_level.chunk_dims[0] = 4096;
-	t->dsinfo_num_elems_on_level.type_id = t->dtypes.h5_id_t;
+	t->dsinfo_num_elems_on_level.type_id = t->dtypes.h5_glb_idx_t;
 	TRY( t->dsinfo_num_elems_on_level.create_prop = h5priv_create_hdf5_property (
 		     f,
 		     H5P_DATASET_CREATE) );
@@ -347,7 +347,7 @@ h5tpriv_open_file (
 	TRY( (f->t = h5priv_calloc (f, 1, sizeof (*f->t))) );
 	h5t_fdata_t* t = f->t;
 
-	t->dtypes.h5_id_t = H5_INT64_T;
+	t->dtypes.h5_glb_idx_t = H5_INT64_T;
 	t->dtypes.h5_int64_t = H5_INT64_T;
 	t->dtypes.h5_float64_t = H5_FLOAT64_T;
 
@@ -549,14 +549,14 @@ h5t_close_mesh (
 h5_err_t
 h5t_set_level (
 	h5_file_t* const f,
-	const h5_id_t level_id
+	const h5t_lvl_idx_t level_id
 	) {
 	h5t_fdata_t* t = f->t;
 
 	if ((level_id < 0) || (level_id >= t->num_levels))
 		return HANDLE_H5_OUT_OF_RANGE_ERR (f, "Level", level_id);
 
-	h5_id_t prev_level = t->cur_level;
+	h5t_lvl_idx_t prev_level = t->cur_level;
 	t->cur_level = level_id;
 
 	if (level_id >= t->num_loaded_levels) {
@@ -578,7 +578,6 @@ h5tpriv_alloc_num_vertices (
 	ssize_t size = num * sizeof (t->vertices[0]);
 	TRY( t->vertices = h5priv_alloc (f, t->vertices, size) );
 	TRY( h5priv_alloc_idxmap (f, &t->map_vertex_g2l, num) );
-	TRY( h5priv_alloc_idlist_items (f, &t->sorted_lvertices, num) );
 
 	return H5_SUCCESS;
 }

src/h5core/h5t_readwrite.c

@@ -176,8 +176,6 @@ read_vertices (
 		     open_mem_space_vertices,
 		     open_file_space_vertices,
 		     t->vertices) );
-	TRY( h5tpriv_sort_vertices (f) );
-	TRY( h5tpriv_rebuild_vertex_indices_mapping (f) );
 
 	return H5_SUCCESS;
 }
@@ -235,8 +233,6 @@ open_file_space_elems (
 	return H5S_ALL;
 }
 
-
-
 static h5_err_t
 read_elems (
 	h5_file_t* const f
@@ -251,11 +247,6 @@ read_elems (
 		     open_mem_space_elems,
 		     open_file_space_elems,
 		     t->glb_elems.data) );
-
-	TRY( h5tpriv_sort_loc_elems (f) );
-	TRY( h5tpriv_rebuild_elem_indices_mapping (f) );
-	TRY( h5tpriv_init_loc_elems_struct (f) );
-
 	return H5_SUCCESS;
 }
 
@@ -284,9 +275,16 @@ h5tpriv_read_mesh (
 	TRY( read_num_levels (f) );
 	TRY( read_num_vertices (f) );
 	TRY( read_num_elems (f) );
+
 	TRY( read_vertices (f) );
+	TRY( h5tpriv_rebuild_vertex_indices_mapping (f) );
+
 	TRY( read_elems (f) );
+	TRY( h5tpriv_rebuild_elem_indices_mapping (f) );
+	TRY( h5tpriv_init_loc_elems_struct (f) );
 	TRY( h5tpriv_update_adjacency_structs (f, 0) );
+	TRY( h5tpriv_init_geom_boundary_info (f, 0) );
+
 	TRY( read_mtags (f) );
 	t->num_loaded_levels = t->num_levels;
 	return H5_SUCCESS;

src/h5core/h5t_readwrite_private.h

@@ -32,4 +32,12 @@ h5tpriv_init_loc_elems_struct (
 	return (*f->t->methods.read->init_loc_elems_struct) (f);
 }
 
+static inline h5_err_t
+h5tpriv_init_geom_boundary_info (
+	h5_file_t* const f,
+	const h5t_lvl_idx_t from_lvl
+	) {
+	return (*f->t->methods.read->init_geom_boundary_info) (f, from_lvl);
+}
+
 #endif

src/h5core/h5t_readwrite_tetm.c

@@ -11,13 +11,16 @@ init_loc_elems_struct (
 	h5t_fdata_t* const t = f->t;
 	h5_loc_idx_t elem_idx = 0;
 	const h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
-	h5_id_t level_idx = 0;
+	h5t_lvl_idx_t level_idx = 0;
 	int num_vertices = h5tpriv_ref_elem_get_num_vertices (t);
-	int num_edges = h5tpriv_ref_elem_get_num_edges (t);
+	int num_facets = h5tpriv_ref_elem_get_num_faces (t, 2);
 	h5_loc_tet_t* loc_elem = t->loc_elems.tets;
 	h5_glb_tet_t* glb_elem = t->glb_elems.tets;
 
-	for (elem_idx = 0; elem_idx < num_elems; elem_idx++, loc_elem++, glb_elem++) {
+	for (elem_idx = 0;
+	     elem_idx < num_elems; elem_idx++, loc_elem++, glb_elem++) {
+		// global element index
+		loc_elem->glb_idx = glb_elem->idx;
 		// local parent index
 		TRY( loc_elem->parent_idx =
 		     h5t_map_glb_elem_idx2loc (f, glb_elem->parent_idx) );
@@ -30,7 +33,7 @@ init_loc_elems_struct (
 		if (elem_idx >= t->num_elems[level_idx]) {
 			level_idx++;
 		}
-		loc_elem->idx = level_idx;
+		loc_elem->level_idx = level_idx;
 
 		// vertex indices
 		TRY( h5t_map_global_vertex_indices2local (
@@ -43,12 +46,93 @@ init_loc_elems_struct (
 		TRY( h5t_map_glb_elem_indices2loc (
 			     f,
 			     glb_elem->neighbor_indices,
-			     num_edges,
+			     num_facets,
 			     loc_elem->neighbor_indices) );
+		// on boundary?
+		int i;
+		for (i=0; i < num_facets; i++) {
+			if (loc_elem->neighbor_indices[i] == -1) {
+				loc_elem->flags |= H5T_BOUNDARY_ELEM_FLAG;
+				break;
+			}
+		}
+	}
+	return H5_SUCCESS;
+}
+
+static h5_err_t
+init_geom_boundary_info (
+	h5_file_t* const f,
+	const h5t_lvl_idx_t from_lvl
+	) {
+	h5t_fdata_t* const t = f->t;
+	h5_loc_idx_t elem_idx = 0;
+	const h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
+	int num_facets = h5tpriv_ref_elem_get_num_faces (t, 2);
+	h5_loc_tet_t* loc_elem = t->loc_elems.tets;
+	h5_glb_tet_t* glb_elem = t->glb_elems.tets;
+
+	for (elem_idx = (from_lvl <= 0) ? 0 : t->num_elems[from_lvl-1];
+	     elem_idx < num_elems; elem_idx++, loc_elem++, glb_elem++) {
+		// on boundary?
+		int i;
+		for (i=0; i < num_facets; i++) {
+			if (loc_elem->neighbor_indices[i] == -1) {
+				loc_elem->flags |= H5T_BOUNDARY_ELEM_FLAG;
+				break;
+			}
+		}
+		if (i == num_facets) {
+			continue; // no facet on boundary
+		}
+
+#if 0
+		// mark elements which are edge- or vertex- adjacent to the boundary
+		// are there more facets on the boundary?
+		int j;
+		for (j = i+1; j < num_facets; j++) {
+			if (loc_elem->neighbor_indices[j] == -1) {
+				break; // found another boundary facet
+			}
+		}
+		// get vertex ID's of vertices on boundary
+		h5_loc_idx_t vertex_indices[4];
+		if (j < num_facets) {
+			// all vertices on boundary
+			vertex_indices[0] = h5tpriv_build_vertex_id (0, elem_idx);
+			vertex_indices[1] = h5tpriv_build_vertex_id (1, elem_idx);
+			vertex_indices[2] = h5tpriv_build_vertex_id (2, elem_idx);
+			vertex_indices[3] = h5tpriv_build_vertex_id (3, elem_idx);
+			num_vertices = 4;
+		} else {
+			// three vertices on boundary
+			// get vertices of edge i
+			h5_loc_idx_t face_idx;
+			face_idx = t->ref_elem->map[2][i][0];
+			vertex_indices[0] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, face_idx);
+			face_idx = t->ref_elem->map[2][i][1];
+			vertex_indices[1] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, face_idx);
+			face_idx = t->ref_elem->map[2][i][2];
+			vertex_indices[2] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, face_idx);
+			num_vertices = 3;
+		}
+		// mark elements
+
+		for (i = 0; i < num_vertices; i++) {
+			const h5_loc_idx_t vertex_idx = vertex_indices[i];
+			h5_idlist_t* list = &t->adjacencies.tv.v[vertex_idx];
+			// set flag
+			for (j=0; j < list->num_items; j++) {
+				h5_loc_idx_t idx = h5tpriv_get_elem_idx (list->items[j]);
+				t->loc_elems.tris[idx].flags |= H5T_BOUNDARY_ELEM_FLAG;
+			}
+		}
+#endif
 	}
 	return H5_SUCCESS;
 }
 
 struct h5t_read_methods h5tpriv_read_tetm_methods = {
-	init_loc_elems_struct
+	init_loc_elems_struct,
+	init_geom_boundary_info,
 };

src/h5core/h5t_readwrite_trim.c

@@ -11,13 +11,15 @@ init_loc_elems_struct (
 	h5t_fdata_t* const t = f->t;
 	h5_loc_idx_t idx = 0;
 	const h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
-	h5_id_t level_idx = 0;
+	h5t_lvl_idx_t level_idx = 0;
 	int num_vertices = h5tpriv_ref_elem_get_num_vertices (t);
-	int num_edges = h5tpriv_ref_elem_get_num_edges (t);
+	int num_facets = h5tpriv_ref_elem_get_num_faces (t,1);
 	h5_loc_triangle_t* loc_elem = t->loc_elems.tris;
 	h5_glb_triangle_t* glb_elem = t->glb_elems.tris;
 
 	for (idx = 0; idx < num_elems; idx++, loc_elem++, glb_elem++) {
+		// global element index
+		loc_elem->glb_idx = glb_elem->idx;
 		// local parent index
 		TRY( loc_elem->parent_idx =
 		     h5t_map_glb_elem_idx2loc (f, glb_elem->parent_idx) );
@@ -30,7 +32,7 @@ init_loc_elems_struct (
 		if (idx >= t->num_elems[level_idx]) {
 			level_idx++;
 		}
-		loc_elem->idx = level_idx;
+		loc_elem->level_idx = level_idx;
 
 		// vertex indices
 		TRY( h5t_map_global_vertex_indices2local (
@@ -43,12 +45,90 @@ init_loc_elems_struct (
 		TRY( h5t_map_glb_elem_indices2loc (
 			     f,
 			     glb_elem->neighbor_indices,
-			     num_edges,
+			     num_facets,
 			     loc_elem->neighbor_indices) );
+		// on boundary?
+		int i;
+		for (i=0; i < num_facets; i++) {
+			if (loc_elem->neighbor_indices[i] == -1) {
+				loc_elem->flags |= H5T_BOUNDARY_ELEM_FLAG;
+				break;
+			}
+		}
+	}
+	return H5_SUCCESS;
+}
+
+static h5_err_t
+init_geom_boundary_info (
+	h5_file_t* const f,
+	const h5t_lvl_idx_t from_lvl
+	) {
+	h5t_fdata_t* const t = f->t;
+	h5_loc_idx_t elem_idx = 0;
+	const h5_loc_idx_t num_elems = t->num_elems[t->num_levels-1];
+	int num_facets = h5tpriv_ref_elem_get_num_faces (t, 1);
+	h5_loc_triangle_t* loc_elem = t->loc_elems.tris;
+	h5_glb_triangle_t* glb_elem = t->glb_elems.tris;
+
+	for (elem_idx = (from_lvl <= 0) ? 0 : t->num_elems[from_lvl-1];
+	     elem_idx < num_elems; elem_idx++, loc_elem++, glb_elem++) {
+		// on boundary?
+		int i;
+		for (i=0; i < num_facets; i++) {
+			if (loc_elem->neighbor_indices[i] == -1) {
+				loc_elem->flags |= H5T_BOUNDARY_ELEM_FLAG;
+				break;
+			}
+		}
+		if (i == num_facets) {
+			continue; // no facet on boundary
+		}
+
+#if 0
+		// mark elements which are edge- or vertex- adjacent to the boundary
+		// are there more facets on the boundary?
+		int j;
+		for (j = i+1; j < num_facets; j++) {
+			if (loc_elem->neighbor_indices[j] == -1) {
+				break; // found another boundary facet
+			}
+		}
+		// get vertex ID's of vertices on boundary
+		h5_loc_idx_t vertex_indices[4];
+		if (j < num_facets) {
+			// all vertices on boundary
+			vertex_indices[0] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, 0);
+			vertex_indices[1] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, 1);
+			vertex_indices[2] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, 2);
+			num_vertices = 3;
+		} else {
+			// two vertices on boundary
+			// get vertices of edge i
+			h5_loc_idx_t face_idx;
+			face_idx = t->ref_elem->map[1][i][0];
+			vertex_indices[0] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, face_idx);
+			face_idx = t->ref_elem->map[1][i][1];
+			vertex_indices[1] = h5tpriv_get_loc_elem_vertex_idx (f, elem_idx, face_idx);
+			num_vertices = 2;
+		}
+		// mark elements
+		for (i = 0; i < num_vertices; i++) {
+			// get upward adjacent elements
+			const h5_loc_idx_t vertex_idx = vertex_indices[i];
+			h5_idlist_t* list = &t->adjacencies.tv.v[vertex_idx];
+			// set flag
+			for (j=0; j < list->num_items; j++) {
+				h5_loc_idx_t idx = h5tpriv_get_elem_idx (list->items[j]);
+				t->loc_elems.tris[idx].flags |= H5T_BOUNDARY_ELEM_FLAG;
+			}
+		}
+#endif
 	}
 	return H5_SUCCESS;
 }
 
 struct h5t_read_methods h5tpriv_read_trim_methods = {
-	init_loc_elems_struct
+	init_loc_elems_struct,
+	init_geom_boundary_info,
 };

src/h5core/h5t_ref_elements.c

@@ -1,4 +1,5 @@
 #include "h5core/h5_core.h"
+#include "h5_core_private.h"
 
 /*
   
@@ -30,6 +31,11 @@
 
 const h5t_ref_elem_t h5t_tri_ref_elem = {
 	2,					// dimension
+	{
+		H5T_ETYPE_VERTEX,
+		H5T_ETYPE_EDGE,
+		H5T_ETYPE_TRIANGLE,
+	},
 	{
 		3,				// #vertices
 		3,				// #edges
@@ -87,6 +93,12 @@ const h5t_ref_elem_t h5t_tri_ref_elem = {
 
 const h5t_ref_elem_t h5t_tet_ref_elem = {
 	3,					// dimension
+	{
+		H5T_ETYPE_VERTEX,
+		H5T_ETYPE_EDGE,
+		H5T_ETYPE_TRIANGLE,
+		H5T_ETYPE_TET,
+	},
 	{
 		4,				// #vertices
 		6,				// #edges
@@ -99,7 +111,7 @@ const h5t_ref_elem_t h5t_tet_ref_elem = {
 		[2] = {3,3,3,3},	       	// #vertices of trinagles
 		[3] = {4}			// #vertices of tets
 	},
-	{	// map sub-elements to vertices
+	{	// map faces to vertices
 		[0] = {{0}, {1}, {2}, {3}},		 	  // 4 vertices
 		[1] = {{0,1}, {0,2}, {1,2}, {0,3}, {1,3}, {2,3}}, // 6 edges
 		[2] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,2,3}},	  // 4 triangles

src/h5core/h5t_ref_elements_private.h

@@ -3,11 +3,21 @@
 
 #define h5tpriv_ref_elem_get_num_vertices(this) (this->ref_elem->num_faces[0])
 #define h5tpriv_ref_elem_get_num_edges(this) (this->ref_elem->num_faces[1])
-#define h5tpriv_ref_elem_get_num_triangles(this) (this->ref_elem->num_faces[2])
+
+#define h5tpriv_ref_elem_get_num_facets(this) \
+	(this->ref_elem->num_faces[this->ref_elem->dim - 1])
+
 #define h5tpriv_ref_elem_get_num_faces(this, dim) (this->ref_elem->num_faces[dim])
 
 #define h5tpriv_ref_elem_get_dim(this) (this->ref_elem->dim)
 
-#define h5tpriv_get_edge_connected_to_vertex(this,face_idx, i) (this->edges_connected_to_vertex[face_idx][i])
-#define h5tpriv_get_facet_connected_to_vertex(this,face_idx, i) (this->edges_connected_to_vertex[face_idx][i])
+#define h5tpriv_get_edge_connected_to_vertex(this,face_idx, i) \
+	(this->edges_connected_to_vertex[face_idx][i])
+
+#define h5tpriv_get_triangles_connected_to_vertex(this,face_idx, i) \
+	(this->triangles_connected_to_vertex[face_idx][i])
+
+#define h5tpriv_ref_elem_get_entity_type(this,dim) \
+	(this->entity_types[dim])
+
 #endif

src/h5core/h5t_retrieve.c

@@ -36,13 +36,149 @@ h5tpriv_elem_is_on_cur_level (
 	h5_generic_loc_elem_t *el // ptr to local element
 	) {
 	h5t_fdata_t* t = f->t;
-	if ( (el->idx > t->cur_level) ||
+	if ( (el->level_idx > t->cur_level) ||
 	     (el->child_idx >= 0 && el->child_idx < num_elems_on_cur_level) ) {
 		return H5_NOK;
 	}
 	return H5_SUCCESS;
 }
 
+
+static h5_loc_id_t
+iterate_elems (
+	h5_file_t* const f,
+	h5t_entity_iterator_t* iter
+	) {
+	if ( h5tpriv_skip_to_next_elem_on_level (f, iter) == H5_NOK) {
+		h5_debug ( f, "Traversing done!" );
+		return H5_NOK;
+	}
+	return h5tpriv_build_elem_id ( iter->elem_idx );
+}
+
+static h5_loc_id_t
+iterate_boundary_elems (
+	h5_file_t* const f,
+	h5t_entity_iterator_t* iter
+	) {
+	do {
+		if ( h5tpriv_skip_to_next_elem_on_level (f, iter) == H5_NOK) {
+			h5_debug ( f, "Traversing done!" );
+			return H5_NOK;
+		}
+	} while (!h5tpriv_is_boundary_elem (f, iter->elem_idx));
+	return h5tpriv_build_elem_id ( iter->elem_idx );
+}
+
+
+
+/*
+  Iterate boundary facets (co-dim 1 entities).
+ */
+static h5_loc_id_t
+iterate_boundary_facets (
+	h5_file_t* const f,
+	h5t_entity_iterator_t* iter
+	) {
+	int num_facets = h5tpriv_ref_elem_get_num_facets (iter) - 1;
+	int dim = h5tpriv_ref_elem_get_dim (iter) - iter->codim;
+	do {
+		if (iter->face_idx >= num_facets) {
+			h5_loc_id_t elem_id;
+			TRY( elem_id = iterate_boundary_elems (f, iter) );
+			if (elem_id == H5_NOK) {
+				return H5_NOK;	// done!
+			}
+			iter->elem_idx = h5tpriv_get_elem_idx (elem_id);
+			iter->face_idx = 0;
+		} else {
+			iter->face_idx++;
+		}
+	} while (! h5tpriv_is_boundary_facet (f, iter->elem_idx, iter->face_idx));
+	int type = h5tpriv_ref_elem_get_entity_type (iter->ref_elem, dim);
+	return h5tpriv_build_id (type, iter->face_idx, iter->elem_idx );
+}
+
+/*!
+  Travere entities with co-dim > 0
+*/
+static h5_loc_id_t
+iterate_faces (
+	h5_file_t* const f,
+	h5t_entity_iterator_t* iter
+	) {
+	h5_idlist_t* entry;
+	int dim = h5tpriv_ref_elem_get_dim (iter) - iter->codim;
+	int num_faces = h5tpriv_ref_elem_get_num_faces (iter, dim) - 1;
+	int i = -1;
+	do {
+		if (iter->face_idx >= num_faces) {
+			if (h5tpriv_skip_to_next_elem_on_level (f, iter) == H5_NOK) {
+				h5_debug (f, "Traversing done!");
+				return H5_NOK;
+			}
+			iter->face_idx = 0;
+		} else {
+			iter->face_idx++;
+		}
+		// Skip already visited faces:
+
+		//  1. Get list of all elements with this face. Actually we
+		//     retrieve a list of entities sorted by the element index.
+		TRY( (iter->find)(f, iter->face_idx, iter->elem_idx, &entry) );
+
+		// 2. Go to first element in list which is on current level
+		i = -1;
+		h5_generic_loc_elem_t *el;
+		do {
+			i++;
+			h5_loc_idx_t idx = h5tpriv_get_elem_idx (entry->items[i]);
+			el = h5tpriv_get_loc_elem (f, idx);
+		} while (h5tpriv_elem_is_on_cur_level (f, el) == H5_NOK);
+
+		// 3. Face already visited if 
+	} while (iter->elem_idx > h5tpriv_get_elem_idx(entry->items[i]));
+	/*
+	  note: in above test iter->elem_idx is always greater or equal to the
+	  compared index. It cannot be smaller since iter->elem_idx is on the
+	  current level and the element index of entry->items[i] is the smallest
+	  element index with the given face on the current level.
+	*/
+	return entry->items[0];
+}
+
+/*
+  Iterate boundary faces
+ */
+static h5_loc_id_t
+iterate_boundary_faces (
+	h5_file_t* const f,
+	h5t_entity_iterator_t* iter
+	) {
+	// TODO!!!
+	int dim = h5tpriv_ref_elem_get_dim (iter) - iter->codim;
+	int num_faces = h5tpriv_ref_elem_get_num_faces (iter, dim) - 1;
+	do {
+		// iterate to next boundary face
+		do {
+			// first iterate over all faces of element, if done
+			// goto next element
+			if (iter->face_idx >= num_faces) {
+				h5_loc_id_t elem_id;
+				TRY( elem_id = iterate_boundary_elems (f, iter) );
+				if (elem_id == H5_NOK) {
+					return H5_NOK;	// done!
+				}
+				iter->face_idx = 0;
+			} else {
+				iter->face_idx++;
+			}
+		} while (! h5tpriv_is_boundary_face (f, dim, iter->elem_idx, iter->face_idx));
+		// Skip already visited faces
+	} while (0);
+	return h5_error_internal (f, __FILE__, __func__, __LINE__);
+}
+
 h5_err_t
 h5t_alloc_entity_iterator (
 	h5_file_t* f,
@@ -53,6 +189,16 @@ h5t_alloc_entity_iterator (
 	return h5t_begin_iterate_entities (f, *iter, codim);
 }
 
+h5_err_t
+h5t_alloc_boundary_face_iterator (
+	h5_file_t* f,
+	h5t_entity_iterator_t** iter,
+	int codim
+	) {
+	TRY( *iter = h5priv_alloc (f, NULL, sizeof (h5t_entity_iterator_t)) );
+	return h5t_begin_iterate_boundary_faces (f, *iter, codim);
+}
+
 h5_err_t
 h5t_release_entity_iterator (
 	h5_file_t* const f,
@@ -61,67 +207,47 @@ h5t_release_entity_iterator (
 	return h5priv_free (f, iter);
 }
 
-
 h5_err_t
 h5t_begin_iterate_entities (
 	h5_file_t* f,
 	h5t_entity_iterator_t* iter,
 	const int codim
 	) {
+	iter->face_idx = -1;
+	iter->elem_idx = -1;
+	iter->codim = codim;
+	iter->ref_elem = f->t->ref_elem;
+
+	if (iter->codim > 0) {
+		iter->iter = iterate_faces;
+	} else if (iter->codim == 0) {
+		iter->iter = iterate_elems;
+	}
 	return h5tpriv_init_iterator (f, iter, codim);
 }
 
-/*!
-  Travere entities with co-dim > 0
-*/
-static h5_loc_id_t
-iterate_faces (
-	h5_file_t* const f,
-	h5t_entity_iterator_t* iter
+h5_err_t
+h5t_begin_iterate_boundary_faces (
+	h5_file_t* f,
+	h5t_entity_iterator_t* iter,
+	const int codim
 	) {
-	h5_idlist_t* entry;
-	h5_size_t i;
-	int dim = h5tpriv_ref_elem_get_dim (iter) - iter->codim;
-	int num_faces = h5tpriv_ref_elem_get_num_faces (iter, dim) - 1;
-	do {
-		if (iter->face_idx >= num_faces) {
-			if (h5tpriv_skip_to_next_elem_on_level (f, iter) == H5_NOK) {
-				h5_debug (f, "Traversing done!");
-				return H5_NOK;
-			}
-			iter->face_idx = 0;
-		} else {
-			iter->face_idx++;
-		}
-		/*
-		  get list of all elements with this face
-		  and skip to first element in list which is on
-		  current level
-		*/
-		TRY( (iter->find)(f, iter->face_idx,
-				  iter->elem_idx, &entry) );
-		i = -1;
-		h5_generic_loc_elem_t *el;
-		do {
-			i++;
-			h5_loc_idx_t idx = h5tpriv_get_elem_idx (entry->items[i]);
-			el = h5tpriv_get_loc_elem (f, idx);
-		} while (h5tpriv_elem_is_on_cur_level (f, el) == H5_NOK);
-	} while (iter->elem_idx != h5tpriv_get_elem_idx(entry->items[i]));
-	
-	return entry->items[0];
-}
+	iter->face_idx = 999; // just a high enough number
+	iter->elem_idx = -1;
+	iter->codim = codim;
+	iter->ref_elem = f->t->ref_elem;
 
-static h5_loc_id_t
-iterate_elems (
-	h5_file_t* const f,
-	h5t_entity_iterator_t*iter
-	) {
-	if ( h5tpriv_skip_to_next_elem_on_level (f, iter) == H5_NOK) {
-		h5_debug ( f, "Traversing done!" );
-		return H5_NOK;
+	if (iter->codim <= 0 || iter->codim > iter->ref_elem->dim) {
+		return h5_error (f, H5_ERR_INVAL,
+				 "Co-dimension requested %d, but must be between %d and %d",
+				 codim, 1, iter->ref_elem->dim);
+	} else if (iter->codim == 1) {
+		iter->iter = iterate_boundary_facets;
 	}
-	return h5tpriv_build_elem_id ( iter->elem_idx );
+	else if (iter->codim > 1) {
+		iter->iter = iterate_boundary_faces;
+	}
+	return H5_SUCCESS;
 }
 
 h5_loc_id_t
@@ -129,13 +255,7 @@ h5t_iterate_entities (
 	h5_file_t* const f,
 	h5t_entity_iterator_t* iter
 	) {
-	if (iter->codim > 0) {
-		return iterate_faces (f, iter);
-	} else if (iter->codim == 0) {
-		return iterate_elems (f, iter);
-	} else {
-		return h5_error_internal (f, __FILE__, __func__, __LINE__);
-	}
+	return (iter->iter (f, iter));
 }
 		   
 h5_err_t
@@ -149,9 +269,10 @@ h5t_end_iterate_entities (
 	iter->codim = -1;
 	iter->ref_elem = NULL;
 	iter->find = NULL;
+	iter->iter = NULL;
 	return H5_SUCCESS;
 }
-		   
+
 h5_err_t
 h5t_get_vertex_coords_by_index (
 	h5_file_t* const f,

src/h5core/h5t_retrieve_tetm.c

@@ -9,10 +9,6 @@ begin_iterate_entities (
 	h5t_entity_iterator_t* const iter,
 	const int codim
 	) {
-	iter->face_idx = -1;
-	iter->elem_idx = -1;
-	iter->codim = codim;
-	iter->ref_elem = f->t->ref_elem;
 	switch (iter->ref_elem->dim - codim) {
 	case 0: // iterate vertices
 		iter->find = h5tpriv_find_tv2;

src/h5core/h5t_retrieve_trim.c

@@ -9,10 +9,6 @@ begin_iterate_entities (
 	h5t_entity_iterator_t* iter,
 	const int codim
 	) {
-	iter->face_idx = -1;
-	iter->elem_idx = -1;
-	iter->codim = codim;
-	iter->ref_elem = f->t->ref_elem;
 	switch (h5tpriv_ref_elem_get_dim (iter) - codim) {
 	case 0: // iterate vertices
 		iter->find = h5tpriv_find_tv2;

src/h5core/h5t_store.c

@@ -121,7 +121,7 @@ assign_glb_elem_data (
 }
 
 
-h5_id_t
+h5t_lvl_idx_t
 h5t_add_level (
 	h5_file_t* const f
 	) {
@@ -216,7 +216,6 @@ h5t_end_store_vertices (
 
 	t->num_vertices[t->cur_level] = t->last_stored_vid+1;
 	TRY( assign_global_vertex_indices (f) );
-	TRY( h5tpriv_sort_vertices (f) );
 	TRY( h5tpriv_rebuild_vertex_indices_mapping (f) );
 	return H5_SUCCESS;
 }
@@ -326,8 +325,6 @@ h5t_end_store_elems (
 
 	t->num_elems[t->cur_level] = t->last_stored_eid+1;
 
-	TRY( h5tpriv_sort_loc_elems (f) );
-
 	/* assign global indices to new indices */
 	TRY( assign_glb_elem_indices (f) );
 

src/h5core/h5t_store_tetm.c

@@ -244,7 +244,7 @@ compute_neighbor_of_face (
 				f, __FILE__, __func__, __LINE__);
 		}
 		if (td->num_items == 1) {
-			// neighbor is coarser or face is on the border
+			// neighbor is coarser or face is on the boundary
 			elem_idx = t->loc_elems.tets[elem_idx].parent_idx;
 			if (elem_idx == -1) {
 				// we are on the level of the macro grid
@@ -267,19 +267,21 @@ compute_neighbor_of_face (
 }
 
 /*
-  New level has been added, compute neighbores for new elements.
+  Compute neighbors for elements on given level.
  */
 static inline h5_err_t
-compute_neighbors_of_new_elems (
-	h5_file_t* const f
+compute_neighbors_of_elems (
+	h5_file_t* const f,
+	h5t_lvl_idx_t level
 	) {
 	h5t_fdata_t * const t = f->t;
-	if (t->cur_level < 0) {
-		// or should we consider this as an error?
-		return H5_SUCCESS;
+	if (level < 0 || level >= t->num_levels) {
+		return h5_error (f, H5_ERR_INVAL,
+				 "level idx %lld out of bound, must be in [%lld,%lld]",
+				 (long long)level, (long long)0, (long long)t->num_levels);
 	}
-	h5_loc_idx_t elem_idx = t->cur_level == 0 ? 0 : t->num_elems[t->cur_level-1];
-	const h5_loc_idx_t last_idx = t->num_elems[t->cur_level] - 1;
+	h5_loc_idx_t elem_idx = level == 0 ? 0 : t->num_elems[level-1];
+	const h5_loc_idx_t last_idx = t->num_elems[level] - 1;
 	h5_loc_tet_t *el = &t->loc_elems.tets[elem_idx];
 	while (elem_idx <= last_idx) {
 		h5_loc_idx_t face_idx = 0;
@@ -300,7 +302,8 @@ end_store_elems (
 	) {
 	h5t_fdata_t* const t = f->t;
 	TRY( h5tpriv_update_adjacency_structs (f, t->cur_level) );
-	TRY( compute_neighbors_of_new_elems (f) );
+	TRY( compute_neighbors_of_elems (f, t->cur_level) );
+	TRY( h5tpriv_init_geom_boundary_info (f, t->cur_level) );
 	return H5_SUCCESS;
 }
 

src/h5core/h5t_store_trim.c

@@ -27,7 +27,7 @@ alloc_triangles (
 		      t->loc_elems.tris,
 		      new * sizeof (t->loc_elems.tris[0]) ) );
 	memset (
-		t->glb_elems.tris + cur,
+		t->loc_elems.tris + cur,
 		-1,
 		(new-cur) * sizeof (t->loc_elems.tris[0]) );
 
@@ -208,7 +208,7 @@ compute_neighbor_of_face (
 				f, __FILE__, __func__, __LINE__);
 		}
 		if (te->num_items == 1) {
-			// neighbor is coarser or face is on the border
+			// neighbor is coarser or face is on the boundary
 			elem_idx = t->loc_elems.tris[elem_idx].parent_idx;
 			if (elem_idx == -1) {
 				// we are on the level of the macro grid
@@ -231,20 +231,22 @@ compute_neighbor_of_face (
 }
 
 /*
-  New level has been added, compute neighbores for new elements.
+  Compute neighbors for elements on given level.
  */
 static inline h5_err_t
-compute_neighbors_of_new_elems (
-	h5_file_t* const f
+compute_neighbors_of_elems (
+	h5_file_t* const f,
+	h5t_lvl_idx_t level
 	) {
 	h5_debug (f, "%s()", __func__);
 	h5t_fdata_t * const t = f->t;
-	if (t->cur_level < 0) {
-		// or should we consider this as an error?
-		return H5_SUCCESS;
+	if (level < 0 || level >= t->num_levels) {
+		return h5_error (f, H5_ERR_INVAL,
+				 "level idx %lld out of bound, must be in [%lld,%lld]",
+				 (long long)level, (long long)0, (long long)t->num_levels);
 	}
-	h5_loc_idx_t elem_idx = t->cur_level == 0 ? 0 : t->num_elems[t->cur_level-1];
-	const h5_loc_idx_t last_idx = t->num_elems[t->cur_level] - 1;
+	h5_loc_idx_t elem_idx = level == 0 ? 0 : t->num_elems[level-1];
+	const h5_loc_idx_t last_idx = t->num_elems[level] - 1;
 	h5_loc_triangle_t *el = &t->loc_elems.tris[elem_idx];
 	while (elem_idx <= last_idx) {
 		h5_loc_idx_t face_idx = 0;
@@ -267,7 +269,8 @@ end_store_elems (
 	h5t_fdata_t* t = f->t;
 
 	TRY( h5tpriv_update_adjacency_structs (f, t->cur_level) );
-	TRY( compute_neighbors_of_new_elems (f) );
+	TRY( compute_neighbors_of_elems (f, t->cur_level) );
+	TRY( h5tpriv_init_geom_boundary_info (f, t->cur_level) );
 	return H5_SUCCESS;
 }
 

src/h5core/h5t_types_private.h

@@ -20,9 +20,11 @@ typedef struct h5_glb_triangle {
 } h5_glb_triangle_t;
 
 typedef struct h5_loc_triangle {
-	h5_loc_idx_t	idx;
+	h5_glb_idx_t	glb_idx;
 	h5_loc_idx_t	parent_idx;
 	h5_loc_idx_t	child_idx;
+	h5t_lvl_idx_t	level_idx;
+	h5t_elem_flags_t	flags;
 	h5_loc_idx_t	vertex_indices[3];
 	h5_loc_idx_t	neighbor_indices[3];
 } h5_loc_triangle_t;
@@ -37,9 +39,11 @@ typedef struct h5_glb_tetrahedron {
 typedef h5_glb_tetrahedron_t	h5_glb_tet_t;
 
 typedef struct h5_loc_tetrahedron {
-	h5_loc_idx_t	idx;
+	h5_loc_idx_t	glb_idx;
 	h5_loc_idx_t	parent_idx;
 	h5_loc_idx_t	child_idx;
+	h5t_lvl_idx_t	level_idx;
+	h5t_elem_flags_t	flags;
 	h5_loc_idx_t	vertex_indices[4];
 	h5_loc_idx_t	neighbor_indices[4];
 } h5_loc_tetrahedron_t;
@@ -53,9 +57,11 @@ typedef struct h5_generic_glb_elem {
 } h5_generic_glb_elem_t;
 
 typedef struct h5_generic_loc_elem {
-	h5_loc_idx_t	idx;
+	h5_glb_idx_t	glb_idx;
 	h5_loc_idx_t	parent_idx;
 	h5_loc_idx_t	child_idx;
+	h5_int32_t	level_idx;
+	h5_int32_t	flags;
 	h5_loc_idx_t	indices[1];
 } h5_generic_loc_elem_t;
 
@@ -73,9 +79,6 @@ typedef union h5_loc_elems {
 
 /*** type ids' for compound types ***/
 typedef struct h5_dtypes {
-	hid_t		h5_id_t;		/* ID's */
-
-	hid_t		h5_glb_id_t;		/* ID's */
 	hid_t		h5_glb_idx_t;		/* ID's */
 
 	hid_t		h5_int64_t;		/* 64 bit signed integer */
@@ -123,10 +126,10 @@ struct h5t_access_methods {
 		h5_file_t* const, const h5_loc_idx_t);
 	h5_loc_id_t (*set_loc_elem_child_idx)(
 		h5_file_t* const, const h5_loc_idx_t, const h5_loc_idx_t);
-	h5_id_t (*get_loc_elem_level_idx)(
+	h5t_lvl_idx_t (*get_loc_elem_level_idx)(
 		h5_file_t* const, const h5_loc_idx_t);
-	h5_id_t (*set_loc_elem_level_idx)(
-		h5_file_t* const, const h5_loc_idx_t, const h5_id_t);
+	h5t_lvl_idx_t (*set_loc_elem_level_idx)(
+		h5_file_t* const, const h5_loc_idx_t, const h5t_lvl_idx_t);
 	h5_loc_idx_t* (*get_loc_elem_vertex_indices)(
 		h5_file_t* const, const h5_loc_idx_t);
 	h5_loc_idx_t (*get_loc_elem_vertex_idx)(
@@ -170,14 +173,20 @@ struct h5t_access_methods {
 	h5_glb_idx_t (*set_glb_elem_neighbor_idx)(
 		h5_file_t* const,
 		const h5_loc_idx_t, const h5_loc_idx_t, const h5_glb_idx_t);
+	h5_err_t (*set_boundary_elem_flag)(h5_file_t* const, const h5_loc_idx_t);
+	h5_err_t (*clear_boundary_elem_flag)(h5_file_t* const, const h5_loc_idx_t);
+	int (*is_boundary_elem)(h5_file_t* const, const h5_loc_idx_t);
+	int (*is_boundary_facet)(h5_file_t* const, const h5_loc_idx_t, const h5_loc_idx_t);
+	int (*is_boundary_face)(h5_file_t* const, const int, const h5_loc_idx_t, const h5_loc_idx_t);
 };
 
 struct h5t_read_methods {
 	h5_err_t (*init_loc_elems_struct)(h5_file_t* const);
+	h5_err_t (*init_geom_boundary_info)(h5_file_t* const, h5t_lvl_idx_t);
 };
 
 struct h5t_adjacency_methods {
-	h5_err_t (*update_internal_structs)(h5_file_t* const, h5_id_t);
+	h5_err_t (*update_internal_structs)(h5_file_t* const, h5t_lvl_idx_t);
 	h5_err_t (*release_internal_structs)(h5_file_t* const);
 	h5_err_t (*get_adjacencies)(
 		h5_file_t * const,
@@ -202,9 +211,9 @@ typedef struct h5t_fdata {
 	h5_id_t		cur_mesh;	/* id of current mesh */
 	h5_id_t		mesh_changed;	/* true if new or has been changed */
 	h5_id_t		num_meshes;	/* number of meshes */
-	h5_id_t		cur_level;	/* id of current level */
-	h5_size_t	num_levels;	/* number of levels */
-	h5_id_t		num_loaded_levels;
+	h5t_lvl_idx_t	cur_level;	/* idx of current level */
+	h5t_lvl_idx_t	num_levels;	/* number of levels */
+	h5t_lvl_idx_t	num_loaded_levels;
 
 	/*** HDF5 IDs ***/
 	hid_t		topo_gid;	/* grp id of mesh in current
@@ -222,8 +231,7 @@ typedef struct h5t_fdata {
 	h5_loc_vertex_t	*vertices;
 	h5_size_t	*num_vertices;
 	h5_idxmap_t	map_vertex_g2l;	/* map global to local idx */
-	h5_idlist_t	sorted_lvertices;
-	h5_loc_idx_t		last_stored_vid;
+	h5_loc_idx_t   	last_stored_vid;
 	h5_dsinfo_t	dsinfo_vertices;
 	h5_dsinfo_t	dsinfo_num_vertices;
 	
@@ -236,14 +244,6 @@ typedef struct h5t_fdata {
 	h5_size_t	*num_elems_on_level;
 	h5_idxmap_t	map_elem_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
-	  ...
-	*/
-	h5_idlist_t	sorted_elems[H5_MAX_VERTICES_PER_ELEM];
-
 	h5_loc_idx_t	last_stored_eid;
 	h5_dsinfo_t	dsinfo_elems;
 	h5_dsinfo_t	dsinfo_num_elems;