same basic functionality added

src/h5_core/h5t_adjacencies.c

@@ -17,6 +17,13 @@
 #include "h5_core/h5_core.h"
 #include "h5_core/h5_core_private.h"
 
+h5_err_t
+_compute_children_of_edge (
+	h5_file_t * const f,
+	h5_id_t local_kid,
+	h5_idlist_t *children
+	);
+
 static h5_err_t
 _calc_tets_of_vertices (
 	h5_file_t * const f
@@ -50,13 +57,13 @@ _cmp_te_node (
 	h5_te_node_key_t *key0 = (h5_te_node_key_t*)a;
 	h5_te_node_key_t *key1 = (h5_te_node_key_t*)b;
 
-	if ( key0->vid0 < key1->vid0 )
+	if ( key0->vids[0] < key1->vids[0] )
 		return -1;
-	if ( key0->vid0 > key1->vid0 )
+	if ( key0->vids[0] > key1->vids[0] )
 		return 1;
-	if ( key0->vid1 < key1->vid1 )
+	if ( key0->vids[1] < key1->vids[1] )
 		return -1;
-	if ( key0->vid1 > key1->vid1 )
+	if ( key0->vids[1] > key1->vids[1] )
 		return 1;
 	return 0;
 }
@@ -65,27 +72,30 @@ static h5_err_t
 _search_te (
 	h5_file_t * const f,
 	h5_te_node_t **node,
-	h5_id_t local_eid,
-	h5_id_t local_vid0,
-	h5_id_t local_vid1
+	h5_id_t local_kid
 	) {
 	h5t_fdata_t *t = f->t;
 	h5t_adjacencies_t *a = &t->adjacencies;
 	void *vnode;
 
+	h5_id_t local_eid = _h5t_get_elem_id ( local_kid );
+	h5_id_t face_id = _h5t_get_face_id ( local_kid );
+	h5_tet_data_t *tet_data = &t->elems_data.tets[local_eid];
+	h5_id_t *vids = (*node)->key.vids;
+	h5_id_t vid;
+	int map[6][2] = { { 0,1 }, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
+
+	vids[0] = tet_data->local_vids[map[face_id][0]];
+	vids[1] = tet_data->local_vids[map[face_id][1]];
+
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
+
 	if ( *node ==NULL ) {
 		TRY ( *node = _h5_alloc ( f, NULL, sizeof(**node) ) );
 		memset ( *node, 0, sizeof(**node) );
 	}
-
-	if ( local_vid0 < local_vid1 ) {
-		(*node)->key.vid0 = local_vid0; 
-		(*node)->key.vid1 = local_vid1;
-	} else {
-		(*node)->key.vid0 = local_vid1;
-		(*node)->key.vid1 = local_vid0;
-	}
-
 	TRY ( vnode = _h5_tsearch (
 		      f,
 		      *node,
@@ -106,22 +116,26 @@ static h5_err_t
 _find_te (
 	h5_file_t * const f,
 	h5_te_node_t **rnode,
-	h5_id_t local_vid0,
-	h5_id_t local_vid1
+	h5_id_t local_kid
 	) {
 	h5t_fdata_t *t = f->t;
 	h5t_adjacencies_t *a = &t->adjacencies;
 	h5_te_node_t node;
 	void *vnode = &node;
 
-	if ( local_vid0 < local_vid1 ) {
-		node.key.vid0 = local_vid0; 
-		node.key.vid1 = local_vid1;
-	} else {
-		node.key.vid0 = local_vid1;
-		node.key.vid1 = local_vid0;
-	}
+	h5_id_t local_eid = _h5t_get_elem_id ( local_kid );
+	h5_id_t face_id = _h5t_get_face_id ( local_kid );
+	h5_tet_data_t *tet_data = &t->elems_data.tets[local_eid];
+	h5_id_t *vids = node.key.vids;
+	h5_id_t vid;
+	int map[6][2] = { { 0,1 }, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
 
+	vids[0] = tet_data->local_vids[map[face_id][0]];
+	vids[1] = tet_data->local_vids[map[face_id][1]];
+
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
 	TRY ( vnode = _h5_tfind (
 		      f,
 		      vnode,
@@ -142,36 +156,27 @@ _calc_tets_of_edges (
 	h5_size_t num_tets = t->num_elems[t->num_levels-1];
 	h5_te_node_t *node = NULL;
 	h5_tet_data_t *tet = &t->elems_data.tets[0];
-	int i, j;
+	h5_id_t edge_id;
 
 	a->te_tree = NULL;
 
 	for ( local_tid = 0; local_tid < num_tets; local_tid++, tet++ ) {
-		/* for each edge of tet */
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[1] ) );
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[2] ) );
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[3] ) );
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[1],
-				   tet->local_vids[2] ) );
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[1],
-				   tet->local_vids[3] ) );
-		TRY ( _search_te ( f, &node, local_tid,
-				   tet->local_vids[2],
-				   tet->local_vids[3] ) );
+		for ( edge_id = 0; edge_id < 6; edge_id++ ) {
+			TRY (
+				_search_te (
+					f,
+					&node,
+					_h5t_build_edge_id (
+						edge_id,
+						local_tid ) )
+				);
+		}
 	}
 	if ( node && node->value.items == NULL ) {
 		_h5_free ( f, node );
 	}
 
-#define DEBUG
+#undef DEBUG
 #ifdef DEBUG
 	node = NULL;
 	tet = &t->elems_data.tets[0];
@@ -207,17 +212,17 @@ _cmp_td_node (
 	h5_td_node_key_t *key0 = (h5_td_node_key_t*)a;
 	h5_td_node_key_t *key1 = (h5_td_node_key_t*)b;
 
-	if ( key0->vid0 < key1->vid0 )
+	if ( key0->vids[0] < key1->vids[0] )
 		return -1;
-	if ( key0->vid0 > key1->vid0 )
+	if ( key0->vids[0] > key1->vids[0] )
 		return 1;
-	if ( key0->vid1 < key1->vid1 )
+	if ( key0->vids[1] < key1->vids[1] )
 		return -1;
-	if ( key0->vid1 > key1->vid1 )
+	if ( key0->vids[1] > key1->vids[1] )
 		return 1;
-	if ( key0->vid2 < key1->vid2 )
+	if ( key0->vids[2] < key1->vids[2] )
 		return -1;
-	if ( key0->vid2 > key1->vid2 )
+	if ( key0->vids[2] > key1->vids[2] )
 		return 1;
 	return 0;
 }
@@ -226,67 +231,43 @@ static h5_err_t
 _search_td (
 	h5_file_t * const f,
 	h5_td_node_t **node,
-	h5_id_t local_eid,
-	h5_id_t local_vid0,
-	h5_id_t local_vid1,
-	h5_id_t local_vid2
+	h5_id_t local_did
 	) {
 	h5t_fdata_t *t = f->t;
 	h5t_adjacencies_t *a = &t->adjacencies;
 	void *vnode;
 
+	h5_id_t local_eid = _h5t_get_elem_id ( local_did );
+	h5_id_t face_id = _h5t_get_face_id ( local_did );
+	h5_tet_data_t *tet_data = &t->elems_data.tets[local_eid];
+	h5_id_t *vids = (*node)->key.vids;
+	h5_id_t vid;
+	int map[4][3] = { { 1,2,3 }, {0,2,3}, {0,1,3}, {0,1,2} };
+
+	vids[0] = tet_data->local_vids[map[face_id][0]];
+	vids[1] = tet_data->local_vids[map[face_id][1]];
+	vids[2] = tet_data->local_vids[map[face_id][2]];
+
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
+	if ( vids[1] > vids[2] ) {
+		vid = vids[1]; vids[1] = vids[2]; vids[2] = vid;
+	}
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
 	if ( *node ==NULL ) {
 		TRY ( *node = _h5_alloc ( f, NULL, sizeof(**node) ) );
 		memset ( *node, 0, sizeof(**node) );
 	}
 
-	if ( local_vid0 <= local_vid1 ) {
-		if ( local_vid0 <= local_vid2 ) {
-			(*node)->key.vid0 = local_vid0;
-			if ( local_vid1 <= local_vid2 ) {
-				(*node)->key.vid1 = local_vid1;
-				(*node)->key.vid2 = local_vid2;
-			} else {
-				(*node)->key.vid1 = local_vid2;
-				(*node)->key.vid2 = local_vid1;
-			}
-		} else { // v2 < v0 && v0 <= v1
-			(*node)->key.vid0 = local_vid2;
-			if ( local_vid0 <= local_vid1 ) {
-				(*node)->key.vid1 = local_vid0;
-				(*node)->key.vid2 = local_vid1;
-			} else {
-				(*node)->key.vid1 = local_vid1;
-				(*node)->key.vid2 = local_vid0;
-			}
-		}
-	} else { // v1 < v0
-		if ( local_vid1 < local_vid2 ) {
-			(*node)->key.vid0 = local_vid1;
-			if ( local_vid0 <= local_vid2 ) {
-				(*node)->key.vid1 = local_vid0;
-				(*node)->key.vid2 = local_vid2;
-			} else {
-				(*node)->key.vid1 = local_vid2;
-				(*node)->key.vid2 = local_vid0;
-			}
-		} else {
-			(*node)->key.vid0 = local_vid2;
-			if ( local_vid0 <= local_vid1 ) {
-				(*node)->key.vid1 = local_vid0;
-				(*node)->key.vid2 = local_vid1;
-			} else {
-				(*node)->key.vid1 = local_vid1;
-				(*node)->key.vid2 = local_vid0;
-			}
-		}
-	}
 	TRY ( vnode = _h5_tsearch (
 		      f,
 		      *node,
 		      (void**)&a->td_tree,
 		      _cmp_td_node ) );
-	h5_te_node_t *rnode = *(h5_te_node_t **)vnode;
+	h5_td_node_t *rnode = *(h5_td_node_t **)vnode;
 	TRY ( _h5_append_to_idlist (
 		      f,
 		      &rnode->value,
@@ -297,6 +278,47 @@ _search_td (
 	return H5_SUCCESS;
 }
 
+static h5_err_t
+_find_td (
+	h5_file_t * const f,
+	h5_td_node_t **rnode,
+	h5_id_t local_did
+	) {
+	h5t_fdata_t *t = f->t;
+	h5t_adjacencies_t *a = &t->adjacencies;
+	h5_td_node_t node;
+	void *vnode = &node;
+
+	h5_id_t local_eid = _h5t_get_elem_id ( local_did );
+	h5_id_t face_id = _h5t_get_face_id ( local_did );
+	h5_tet_data_t *tet_data = &t->elems_data.tets[local_eid];
+	h5_id_t *vids = node.key.vids;
+	h5_id_t vid;
+	int map[4][3] = { { 1,2,3 }, {0,2,3}, {0,1,3}, {0,1,2} };
+
+	vids[0] = tet_data->local_vids[map[face_id][0]];
+	vids[1] = tet_data->local_vids[map[face_id][1]];
+	vids[2] = tet_data->local_vids[map[face_id][2]];
+
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
+	if ( vids[1] > vids[2] ) {
+		vid = vids[1]; vids[1] = vids[2]; vids[2] = vid;
+	}
+	if ( vids[0] > vids[1] ) {
+		vid = vids[0]; vids[0] = vids[1]; vids[1] = vid;
+	}
+	TRY ( vnode = _h5_tfind (
+		      f,
+		      vnode,
+		      (void**)&a->td_tree,
+		      _cmp_td_node ) );
+	*rnode = *(h5_td_node_t **)vnode;
+
+	return H5_SUCCESS;
+}
+
 static h5_err_t
 _calc_tets_of_triangles (
 	h5_file_t * const f
@@ -307,28 +329,21 @@ _calc_tets_of_triangles (
 	h5_size_t num_tets = t->num_elems[t->num_levels-1];
 	h5_td_node_t *node = NULL;
 	h5_tet_data_t *tet = &t->elems_data.tets[0];
+	h5_id_t face_id;
 
 	a->td_tree = NULL;
 
 	for ( local_tid = 0; local_tid < num_tets; local_tid++, tet++ ) {
-		/* for each triangle of tet */
-		TRY ( _search_td ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[1],
-				   tet->local_vids[2] ) );
-		TRY ( _search_td ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[1],
-				   tet->local_vids[3] ) );
-		TRY ( _search_td ( f, &node, local_tid,
-				   tet->local_vids[0],
-				   tet->local_vids[2],
-				   tet->local_vids[3] ) );
-		TRY ( _search_td ( f, &node,
-				   local_tid,
-				   tet->local_vids[1], 
-				   tet->local_vids[2],
-				   tet->local_vids[3] ) );
+		for ( face_id = 0; face_id < 4; face_id++ ) {
+			TRY (
+				_search_td (
+					f,
+					&node,
+					_h5t_build_triangle_id (
+						face_id,
+						local_tid) )
+				);
+		}
 	}
 	if ( node && node->value.items == NULL ) {
 		_h5_free ( f, node );
@@ -345,5 +360,219 @@ _h5t_rebuild_adj_data (
 	TRY ( _calc_tets_of_edges ( f ) );
 	TRY ( _calc_tets_of_triangles ( f ) );
 
+	h5_id_t local_kid = _h5t_build_edge_id ( 0, 0 );
+	h5_idlist_t children;
+
+	f->t->cur_level = 2;
+	memset ( &children, 0, sizeof(children) );
+	_compute_children_of_edge (
+		f,
+		local_kid,
+		&children
+	);
+	f->t->cur_level = 0;
+	return H5_SUCCESS;
+}
+
+h5_err_t
+_compute_direct_children_of_edge (
+	h5_file_t * const f,
+	h5_id_t local_kid,
+	h5_id_t	kids[4]
+	) {
+	h5_id_t local_eid = _h5t_get_elem_id ( local_kid );
+	h5_id_t face_id = _h5t_get_face_id ( local_kid );
+
+	if ( face_id == 0 ) {
+		kids[0] = _h5t_build_edge_id ( 3, local_eid+0 );
+		kids[1] = _h5t_build_edge_id ( 3, local_eid+1 );
+	} else if ( face_id == 1 ) {
+		kids[0] = _h5t_build_edge_id ( 5, local_eid+1 );
+		kids[1] = _h5t_build_edge_id ( 5, local_eid+2 );
+	} else if ( face_id == 2 ) {
+		kids[0] = _h5t_build_edge_id ( 2, local_eid+0 );
+		kids[1] = _h5t_build_edge_id ( 2, local_eid+2 );
+	} else if ( face_id == 3 ) {
+		kids[0] = _h5t_build_edge_id ( 0, local_eid+0 );
+		kids[1] = _h5t_build_edge_id ( 0, local_eid+3 );
+	} else if ( face_id == 4 ) {
+		kids[0] = _h5t_build_edge_id ( 4, local_eid+1 );
+		kids[1] = _h5t_build_edge_id ( 4, local_eid+3 );
+	} else if ( face_id == 5 ) {
+		kids[0] = _h5t_build_edge_id ( 1, local_eid+2 );
+		kids[1] = _h5t_build_edge_id ( 1, local_eid+3 );
+	} else {
+		return  h5_error_internal ( f, __FILE__, __func__, __LINE__ ); 
+	}
+	return H5_SUCCESS;
+}
+
+/*
+def _compute_children_of_edge ( E, children ):
+	for E' = in =T(E):
+		T' = tetrahedron bound to E'
+		if T' has no children on level ≤ L: next
+		(E'0, E'1) = _compute_direct_children_of edge ( E' )
+		for i in [0,1]:
+			T'i = tetrahedron bound to D'i
+			if T'i has children on level ≤ L:
+				_compute_children_of_edge ( E'i, children )
+			else:
+				children.append ( E'i ) 
+*/
+h5_err_t
+_compute_children_of_edge (
+	h5_file_t * const f,
+	h5_id_t local_kid,
+	h5_idlist_t *children
+	) {
+	h5t_fdata_t *t = f->t;
+	h5_te_node_t *te;
+	int i,k;
+
+	TRY ( _find_te ( f, &te, local_kid ) );
+	for ( k = 0; k < te->value.num_items; k++ ) {
+		h5_id_t local_eid = _h5t_get_elem_id ( te->value.items[k] );
+		h5_tet_data_t *tet = &t->elems_data.tets[local_eid];
+
+		/* no children? */
+		if ( tet->local_child_eid < 0 )
+			continue;
+		/* no children on current level? */
+		h5_tet_data_t *child_tet = &t->elems_data.tets[tet->local_child_eid];
+		if ( (tet->local_child_eid >= 0) &&
+		     (child_tet->level_id > t->cur_level) )
+			continue;
+
+		h5_id_t local_kids[2];
+		TRY ( _compute_direct_children_of_edge (
+			f,
+			te->value.items[k],	/* parent triangle */
+			local_kids ) );		/* result */
+		for ( i = 0; i < 4; i++ ) {
+			local_eid = _h5t_get_elem_id ( local_kids[i] );
+			tet = &t->elems_data.tets[local_eid];
+
+			/* no children? */
+			if ( tet->local_child_eid < 0 ) {
+				TRY ( _h5_append_to_idlist (
+					      f, children, local_kids[i] )
+					);
+			}
+			/* no children on current level? */
+			child_tet = &t->elems_data.tets[tet->local_child_eid];
+			if ( child_tet->level_id > t->cur_level) {
+				TRY ( _h5_append_to_idlist (
+					      f, children, local_kids[i] ) );
+			} else {
+				TRY ( _compute_children_of_edge (
+					      f, local_kids[i], children ) );
+			}
+		}
+	}
+	return H5_SUCCESS;
+}
+
+
+h5_err_t
+_compute_direct_children_of_triangle (
+	h5_file_t * const f,
+	h5_id_t local_did,
+	h5_id_t	children_dids[4]
+	) {
+
+	h5_id_t local_eid = _h5t_get_elem_id ( local_did );
+	h5_id_t face_id = _h5t_get_face_id ( local_did );
+	if ( face_id == 0 ) {
+		children_dids[0] = _h5t_build_triangle_id ( 0, local_eid+1 );
+		children_dids[1] = _h5t_build_triangle_id ( 0, local_eid+2 );
+		children_dids[2] = _h5t_build_triangle_id ( 0, local_eid+3 );
+		children_dids[3] = _h5t_build_triangle_id ( 0, local_eid+7 );
+	} else if ( face_id == 1 ) {
+		children_dids[0] = _h5t_build_triangle_id ( 3, local_eid+0 );
+		children_dids[1] = _h5t_build_triangle_id ( 3, local_eid+2 );
+		children_dids[2] = _h5t_build_triangle_id ( 3, local_eid+3 );
+		children_dids[3] = _h5t_build_triangle_id ( 3, local_eid+6 );
+	} else if ( face_id == 2 ) {
+		children_dids[0] = _h5t_build_triangle_id ( 2, local_eid+0 );
+		children_dids[1] = _h5t_build_triangle_id ( 2, local_eid+1 );
+		children_dids[2] = _h5t_build_triangle_id ( 2, local_eid+3 );
+		children_dids[3] = _h5t_build_triangle_id ( 1, local_eid+4 );
+	} else if ( face_id == 3 ) {
+		children_dids[0] = _h5t_build_triangle_id ( 1, local_eid+0 );
+		children_dids[1] = _h5t_build_triangle_id ( 1, local_eid+1 );
+		children_dids[2] = _h5t_build_triangle_id ( 1, local_eid+2 );
+		children_dids[3] = _h5t_build_triangle_id ( 2, local_eid+5 );
+	} else {
+		return  h5_error_internal ( f, __FILE__, __func__, __LINE__ ); 
+	}
+	return H5_SUCCESS;
+}
+
+/*
+def _compute_children_of_triangle ( D, children ):
+	for D' in T(D):
+		T' = tetrahedron bound to D'
+		if T' has no children on level ≤ L: next
+		(D'0,D'1,D'2,D'3) = _compute_direct_children_of_triangle ( D' )
+		for i in [0,1,2,3]:
+			T'i = tetrahedron bound to D'i
+			if T'i has children on level ≤ L:
+				_compute_children_of_triangle ( D'i, children )
+			else:
+				children.append ( D'i) 
+*/
+h5_err_t
+_compute_children_of_triangle (
+	h5_file_t * const f,
+	h5_id_t local_did,
+	h5_idlist_t *children 
+	) {
+
+	h5t_fdata_t *t = f->t;
+	h5_td_node_t *td;
+	int i, k;
+
+	TRY ( _find_td ( f, &td, local_did ) );
+	for ( k = 0; k < td->value.num_items; k++ ) {
+		h5_id_t local_eid = _h5t_get_elem_id ( td->value.items[k] );
+		h5_tet_data_t *tet = &t->elems_data.tets[local_eid];
+
+		/* no children? */
+		if ( tet->local_child_eid < 0 )
+			continue;
+		/* no children on current level? */
+		h5_tet_data_t *child_tet = &t->elems_data.tets[tet->local_child_eid];
+		if ( (tet->local_child_eid >= 0) &&
+		     (child_tet->level_id > t->cur_level) )
+			continue;
+
+		h5_id_t local_dids[4];
+		TRY ( _compute_direct_children_of_triangle (
+			f,
+			td->value.items[k],	/* parent triangle */
+			local_dids ) );		/* result */
+		for ( i = 0; i < 4; i++ ) {
+			local_eid = _h5t_get_elem_id ( local_dids[i] );
+			tet = &t->elems_data.tets[local_eid];
+
+			/* no children? */
+			if ( tet->local_child_eid < 0 ) {
+				TRY ( _h5_append_to_idlist (
+					      f, children, local_dids[i] )
+					);
+			}
+			/* no children on current level? */
+			child_tet = &t->elems_data.tets[tet->local_child_eid];
+			if ( child_tet->level_id > t->cur_level) {
+				TRY ( _h5_append_to_idlist (
+					      f, children, local_dids[i] ) );
+			} else {
+				TRY ( _compute_children_of_triangle (
+					      f, local_dids[i], children ) );
+			}
+		}
+	}
+
 	return H5_SUCCESS;
 }