From 8442f85127ed9e9dbe3ada7e3aa94dc846d2af83 Mon Sep 17 00:00:00 2001 From: Achim Gsell Date: Thu, 14 May 2009 16:28:53 +0000 Subject: [PATCH] same basic functionality added --- src/h5_core/h5t_adjacencies.c | 461 +++++++++++++++++++++++++--------- 1 file changed, 345 insertions(+), 116 deletions(-) diff --git a/src/h5_core/h5t_adjacencies.c b/src/h5_core/h5t_adjacencies.c index 73ec377..efe21e7 100644 --- a/src/h5_core/h5t_adjacencies.c +++ b/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; }