diff --git a/.gitattributes b/.gitattributes index a4e5760..c086e80 100644 --- a/.gitattributes +++ b/.gitattributes @@ -431,6 +431,16 @@ src/h5core/h5t_tags.c -text src/h5core/h5t_tags_private.h -text src/h5core/h5t_tetm_adjacencies.c -text src/h5core/h5t_tetm_adjacencies_private.h -text +src/h5core/h5t_tetm_retrieve.c -text +src/h5core/h5t_tetm_retrieve_private.h -text +src/h5core/h5t_tetm_store.c -text +src/h5core/h5t_tetm_store_private.h -text +src/h5core/h5t_trim_adjacencies.c -text +src/h5core/h5t_trim_adjacencies_private.h -text +src/h5core/h5t_trim_retrieve.c -text +src/h5core/h5t_trim_retrieve_private.h -text +src/h5core/h5t_trim_store.c -text +src/h5core/h5t_trim_store_private.h -text src/h5core/h5t_types_private.h -text src/h5core/h5u_errorhandling_private.h -text src/h5core/h5u_readwrite.c -text diff --git a/src/h5core/h5t_tetm_retrieve.c b/src/h5core/h5t_tetm_retrieve.c new file mode 100644 index 0000000..f927487 --- /dev/null +++ b/src/h5core/h5t_tetm_retrieve.c @@ -0,0 +1,38 @@ +#include +#include + +#include "h5core/h5_core.h" +#include "h5_core_private.h" + +static h5_err_t +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_element = f->t->ref_element; + switch (iter->ref_element->dim - codim) { + case 0: // iterate vertices + iter->find = h5tpriv_find_tv2; + return h5tpriv_skip_to_next_elem_on_level (f, iter); + case 1: // iterate edges + iter->find = h5tpriv_find_te2; + return h5tpriv_skip_to_next_elem_on_level (f, iter); + case 2: // iterate faces + iter->find = h5tpriv_find_td2; + return h5tpriv_skip_to_next_elem_on_level (f, iter); + case 3: // iterate elems + iter->find = NULL; + return H5_SUCCESS; + default: + return h5_error_internal (f, __FILE__, __func__, __LINE__); + } +} + +struct h5t_retrieve_methods h5tpriv_tetm_retrieve_methods = { + begin_iterate_entities +}; + diff --git a/src/h5core/h5t_tetm_retrieve_private.h b/src/h5core/h5t_tetm_retrieve_private.h new file mode 100644 index 0000000..b1c3b6d --- /dev/null +++ b/src/h5core/h5t_tetm_retrieve_private.h @@ -0,0 +1,5 @@ +#ifndef __H5T_TETM_RETRIEVE_PRIVATE_H +#define __H5T_TETM_RETRIEVE_PRIVATE_H + +extern struct h5t_retrieve_methods h5tpriv_tetm_retrieve_methods; +#endif diff --git a/src/h5core/h5t_tetm_store.c b/src/h5core/h5t_tetm_store.c new file mode 100644 index 0000000..e8d99bc --- /dev/null +++ b/src/h5core/h5t_tetm_store.c @@ -0,0 +1,246 @@ +#include +#include + +#include "h5core/h5_core.h" +#include "h5_core_private.h" + +static h5_err_t +alloc_tets ( + h5_file_t* const f, + const size_t cur, + const size_t new + ) { + h5t_fdata_t *t = f->t; + const size_t nvertices = 4; + + /* alloc mem for elements */ + TRY ( t->elems.tets = h5priv_alloc ( + f, + t->elems.tets, + new * sizeof(t->elems.tets[0]) ) ); + memset ( + t->elems.tets + cur, + -1, + (new-cur) * sizeof(t->elems.tets[0]) ); + + /* alloc mem for local data of elements */ + TRY ( t->elems_ldta = h5priv_alloc ( + f, + t->elems_ldta, + new * sizeof (t->elems_ldta[0]) ) ); + memset ( + t->elems_ldta + cur, + -1, + (new-cur) * sizeof (t->elems_ldta[0]) ); + + /* alloc mem for local vertex IDs of elements */ + TRY ( t->elems_lvids = h5priv_alloc ( + f, + t->elems_lvids, + new * sizeof(t->elems_lvids[0]) * nvertices ) ); + memset ( + t->elems_lvids + cur * nvertices, + -1, + (new-cur) * sizeof(t->elems_lvids[0]) * nvertices ); + + /* re-init pointer to local vertex id in local data structure */ + h5_id_t *p = t->elems_lvids; + h5_id_t id; + for ( id = 0; id < new; id++, p+=nvertices ) { + t->elems_ldta[id].local_vids = p; + } + + /* alloc mem for global to local ID mapping */ + TRY ( h5priv_alloc_idmap ( f, &t->map_elem_g2l, new ) ); + + return H5_SUCCESS; +} + + +static h5_err_t +get_direct_children_of_edge ( + h5_file_t* const f, + h5_id_t face_idx, + h5_id_t elem_idx, + h5_id_t children[2] + ) { + /* + Please read the note about the offsets in the corresponding file + for triangle meshes. + */ + int offs[6][2] = { {0,1}, // edge 0 + {0,2}, // edge 1 + {1,2}, // edge 2 + {0,3}, // edge 3 + {1,3}, // edge 4 + {2,3} // edge 5 + }; + int num_faces = f->t->ref_element->num_faces[1]; + if ( ( face_idx < 0 ) || ( face_idx >= num_faces ) ) { + return h5_error_internal ( f, __FILE__, __func__, __LINE__ ); + } + children[0] = h5tpriv_build_edge_id ( face_idx, elem_idx+offs[face_idx][0] ); + children[1] = h5tpriv_build_edge_id ( face_idx, elem_idx+offs[face_idx][1] ); + return H5_SUCCESS; +} + +/*! + Refine edge. Store vertex, if new. + + Function can be used with tetrahedral and triangle meshes. + + \return local index of vertex +*/ +static h5_id_t +bisect_edge ( + h5_file_t* const f, + h5_id_t face_id, + h5_id_t elem_id + ) { + h5t_fdata_t *t = f->t; + h5_id_t vids[2]; + h5_idlist_t *retval; + /* + get all elements sharing the given edge + */ + TRY( h5tpriv_find_te2 (f, face_id, elem_id, &retval) ); + /* + check wether one of the found elements has been refined + */ + size_t i; + for ( i = 0; i < retval->num_items; i++ ) { + h5_id_t local_id = h5tpriv_get_elem_idx ( retval->items[i] ); + h5_elem_ldta_t *tet = &t->elems_ldta[local_id]; + if ( tet->local_child_eid >= 0 ) { + /* + this element has been refined! + return bisecting point + */ + h5_id_t face_id = h5tpriv_get_face_idx ( + retval->items[i] ); + h5_id_t kids[2], edge0[2], edge1[2]; + TRY ( get_direct_children_of_edge ( + f, + face_id, + tet->local_child_eid, + kids ) ); + TRY ( h5t_get_vertex_indices_of_edge ( f, kids[0], edge0 ) ); + TRY ( h5t_get_vertex_indices_of_edge ( f, kids[1], edge1 ) ); + if ( (edge0[0] == edge1[0]) || (edge0[0] == edge1[1]) ) + return edge0[0]; + else + return edge0[1]; + } + } + /* + None of the elements has been refined -> add new vertex. + */ + h5_float64_t *P0 = t->vertices[vids[0]].P; + h5_float64_t *P1 = t->vertices[vids[1]].P; + h5_float64_t P[3]; + + P[0] = ( P0[0] + P1[0] ) / 2.0; + P[1] = ( P0[1] + P1[1] ) / 2.0; + P[2] = ( P0[2] + P1[2] ) / 2.0; + + return h5t_store_vertex ( f, -1, P ); +} + +/*! + Refine tetrahedron \c elem_idx + + \return Local id of first new tetrahedron or \c -1 +*/ +static h5_id_t +refine_tet ( + h5_file_t* const f, + const h5_id_t elem_idx + ) { + h5t_fdata_t* t = f->t; + h5_id_t vertices[10]; + h5_id_t elem_idx_of_first_child; + h5_elem_ldta_t *tet = &t->elems_ldta[elem_idx]; + + if ( tet->local_child_eid >= 0 ) + return h5_error ( + f, + H5_ERR_INVAL, + "Tetrahedron %lld already refined.", + elem_idx ); + vertices[0] = tet->local_vids[0]; + vertices[1] = tet->local_vids[1]; + vertices[2] = tet->local_vids[2]; + vertices[3] = tet->local_vids[3]; + + vertices[4] = bisect_edge (f, 0, elem_idx); + vertices[5] = bisect_edge (f, 1, elem_idx); + vertices[6] = bisect_edge (f, 3, elem_idx); + vertices[7] = bisect_edge (f, 2, elem_idx); + vertices[8] = bisect_edge (f, 4, elem_idx); + vertices[9] = bisect_edge (f, 5, elem_idx); + + /* + add new tets + */ + h5_id_t new_elem[4]; + + new_elem[0] = vertices[0]; // child 0 + new_elem[1] = vertices[4]; + new_elem[2] = vertices[5]; + new_elem[3] = vertices[6]; + TRY( elem_idx_of_first_child = h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[4]; // child 1 + new_elem[1] = vertices[1]; + new_elem[2] = vertices[7]; + new_elem[3] = vertices[8]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[5]; // child 2 + new_elem[1] = vertices[7]; + new_elem[2] = vertices[2]; + new_elem[3] = vertices[9]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[6]; // child 3 + new_elem[1] = vertices[8]; + new_elem[2] = vertices[9]; + new_elem[3] = vertices[3]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[4]; // child 4 + new_elem[1] = vertices[5]; + new_elem[2] = vertices[6]; + new_elem[3] = vertices[8]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[4]; // child 5 + new_elem[1] = vertices[5]; + new_elem[2] = vertices[7]; + new_elem[3] = vertices[8]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[5]; // child 6 + new_elem[1] = vertices[7]; + new_elem[2] = vertices[8]; + new_elem[3] = vertices[9]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[5]; // child 7 + new_elem[1] = vertices[6]; + new_elem[2] = vertices[8]; + new_elem[3] = vertices[9]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + t->elems.tets[elem_idx].global_child_eid = elem_idx_of_first_child; + t->elems_ldta[elem_idx].local_child_eid = elem_idx_of_first_child; + t->num_elems_on_level[t->cur_level]--; + + return elem_idx_of_first_child; +} + +struct h5t_store_methods h5tpriv_tetm_store_methods = { + alloc_tets, + refine_tet, + get_direct_children_of_edge +}; diff --git a/src/h5core/h5t_tetm_store_private.h b/src/h5core/h5t_tetm_store_private.h new file mode 100644 index 0000000..8b4d785 --- /dev/null +++ b/src/h5core/h5t_tetm_store_private.h @@ -0,0 +1,5 @@ +#ifndef __H5T_TETM_STORE_PRIVATE_H +#define __H5T_TETM_STORE_PRIVATE_H + +extern struct h5t_store_methods h5tpriv_tetm_store_methods; +#endif diff --git a/src/h5core/h5t_trim_adjacencies.c b/src/h5core/h5t_trim_adjacencies.c new file mode 100644 index 0000000..b6c7dd5 --- /dev/null +++ b/src/h5core/h5t_trim_adjacencies.c @@ -0,0 +1,383 @@ +/* + Copyright 2006-2010 + Paul Scherrer Institut, Villigen, Switzerland; + Achim Gsell + All rights reserved. + + Authors + Achim Gsell + + Warning + This code is under development. +*/ + +#include +#include "h5core/h5_core.h" +#include "h5_core_private.h" + +/* + compute T(V) from current level up to highest levels. +*/ +static h5_err_t +compute_elems_of_vertices ( + h5_file_t * const f + ) { + h5t_fdata_t *t = f->t; + h5_id_t elem_idx = (t->cur_level <= 0) ? + 0 : t->num_elems[t->cur_level-1]; + h5_elem_ldta_t *el = &t->elems_ldta[elem_idx]; + h5_id_t num_elems = t->num_elems[t->num_levels-1]; + for (;elem_idx < num_elems; elem_idx++, el++) { + int face_idx; + int num_faces = t->ref_element->num_faces[0]; + for (face_idx = 0; face_idx < num_faces; face_idx++) { + h5_id_t vidx = el->local_vids[face_idx]; + TRY( h5priv_append_to_idlist ( + f, + &t->vertices_data[vidx].tv, + h5tpriv_build_vertex_id ( + face_idx, elem_idx)) ); + } + } + return H5_SUCCESS; +} + +/* + Compute T(E) from current level up to highest levels. + */ +static h5_err_t +compute_elems_of_edges ( + h5_file_t * const f + ) { + h5t_fdata_t *t = f->t; + h5_id_t elem_idx = (t->cur_level <= 0) ? + 0 : t->num_elems[t->cur_level-1]; + h5_elem_ldta_t *el = &t->elems_ldta[elem_idx]; + h5_id_t num_elems = t->num_elems[t->num_levels-1]; + h5_idlist_t *retval = NULL; + TRY( h5tpriv_resize_te_htab (f, 4*(num_elems-elem_idx)) ); + for (;elem_idx < num_elems; elem_idx++, el++) { + int face_idx; + int num_faces = t->ref_element->num_faces[1]; + for (face_idx = 0; face_idx < num_faces; face_idx++) { + TRY ( h5tpriv_search_te2 ( + f, face_idx, elem_idx, &retval ) ); + } + } + return H5_SUCCESS; +} + +static h5_err_t +rebuild_internal_structs ( + h5_file_t * const f + ) { + clock_t t1 = clock(); + TRY ( compute_elems_of_vertices ( f ) ); + clock_t t2 = clock(); + fprintf ( stderr, "compute_elems_of_vertices(): %f\n", + (float)(t2-t1)/CLOCKS_PER_SEC ); + t1 = clock(); + TRY ( compute_elems_of_edges ( f ) ); + t2 = clock(); + fprintf ( stderr, "compute_elems_of_edge(): %f\n", + (float)(t2-t1)/CLOCKS_PER_SEC ); + + return H5_SUCCESS; +} + +static h5_err_t +compute_children_of_edge ( + h5_file_t * const f, + h5_id_t kid, + h5_idlist_t *children + ) { + h5t_fdata_t *t = f->t; + h5_idlist_t *te; + + TRY ( h5tpriv_find_te2 ( + f, + h5tpriv_get_face_idx ( kid ), + h5tpriv_get_elem_idx ( kid ), + &te ) + ); + h5_id_t *edge = te->items; + h5_id_t *end = te->items+te->num_items; + for ( ; edge < end; edge++ ) { + h5_id_t elem_idx = h5tpriv_get_elem_idx ( *edge ); + h5_id_t face_idx = h5tpriv_get_face_idx ( *edge ); + h5_elem_ldta_t *el = &t->elems_ldta[elem_idx]; + if ( h5tpriv_elem_is_on_cur_level ( f, el ) == H5_OK ) { + TRY ( h5priv_append_to_idlist ( + f, children, *edge ) + ); + } else { + h5_id_t kids[2]; + TRY ( t->methods.store->get_direct_children_of_edge ( + f, + face_idx, + el->local_child_eid, + kids ) ); + TRY ( compute_children_of_edge ( + f, kids[0], children ) ); + TRY ( compute_children_of_edge ( + f, kids[1], children ) ); + } + } + return H5_SUCCESS; +} + +/* + Compute all sections of an edge. +*/ +static h5_err_t +compute_sections_of_edge ( + h5_file_t * const f, + h5_id_t kid, + h5_idlist_t *children + ) { + h5t_fdata_t *t = f->t; + h5_idlist_t *te; + + TRY ( h5tpriv_find_te2 ( + f, + h5tpriv_get_face_idx ( kid ), + h5tpriv_get_elem_idx ( kid ), + &te ) + ); + h5_id_t *edge = te->items; + h5_id_t *end = te->items+te->num_items; + int refined = 0; + for ( ; edge < end; edge++ ) { + h5_id_t eid = h5tpriv_get_elem_idx ( *edge ); + h5_id_t face_id = h5tpriv_get_face_idx ( *edge ); + h5_elem_ldta_t *el = &t->elems_ldta[eid]; + if ( ! h5tpriv_elem_is_on_cur_level ( f, el ) == H5_OK ) { + refined = 1; + h5_id_t kids[2]; + TRY ( t->methods.store->get_direct_children_of_edge ( + f, + face_id, + el->local_child_eid, + kids ) ); + TRY ( compute_sections_of_edge ( + f, kids[0], children ) ); + TRY ( compute_sections_of_edge ( + f, kids[1], children ) ); + } + } + if ( ! refined ) { + TRY ( h5priv_append_to_idlist ( f, children, te->items[0] ) ); + } + return H5_SUCCESS; +} + + +/* + map edge ID to unique ID + if unique ID not in list: add +*/ +static h5_err_t +add_edge ( + h5_file_t * const f, + h5_idlist_t *list, + h5_id_t face_idx, + h5_id_t elem_idx + ) { + h5_idlist_t *te; + TRY ( h5tpriv_find_te2 ( f, face_idx, elem_idx, &te ) ); + TRY ( h5priv_search_idlist ( f, list, te->items[0] ) ); + return H5_SUCCESS; +} + + +static h5_err_t +get_edges_upadjacent_to_vertex ( + h5_file_t * const f, + const h5_id_t entity_id, + h5_idlist_t **list + ) { + TRY ( h5priv_alloc_idlist ( f, list, 8 ) ); + + h5t_fdata_t *t = f->t; + h5_idlist_t *tv = &t->vertices_data[entity_id].tv; + h5_size_t i; + + h5_id_t *vertex_idp = tv->items; + for ( i = 0; i < tv->num_items; i++, vertex_idp++ ) { + h5_id_t elem_idx = h5tpriv_get_elem_idx ( *vertex_idp ); + h5_id_t face_idx = h5tpriv_get_face_idx ( *vertex_idp ); + h5_elem_ldta_t *el = &t->elems_ldta[elem_idx]; + + if ( h5tpriv_elem_is_on_cur_level ( f, el ) == H5_NOK ) { + continue; + } + /* + upward adjacend edges to vertices according reference element + */ + int map[3][2] = { {0,1}, // edge 0, 1 are adjacent to vertex 0 + {0,2}, // edge 0, 2 are adjacent to vertex 1 + {1,2} // edge 1, 2 are adjacent to vertex 2 + }; + TRY ( add_edge ( f, *list, map[face_idx][0], elem_idx ) ); + TRY ( add_edge ( f, *list, map[face_idx][1], elem_idx ) ); + } + return H5_SUCCESS; +} + +static h5_err_t +get_triangles_upadjacent_to_vertex ( + h5_file_t * const f, + const h5_id_t entity_id, + h5_idlist_t **list + ) { + TRY ( h5priv_alloc_idlist ( f, list, 8 ) ); + + h5t_fdata_t *t = f->t; + h5_idlist_t *tv = &t->vertices_data[entity_id].tv; + h5_size_t i; + h5_id_t *vertex_id = tv->items; + for ( i = 0; i < tv->num_items; i++, vertex_id++ ) { + h5_id_t elem_id = h5tpriv_get_elem_idx ( *vertex_id ); + h5_elem_ldta_t *el = &t->elems_ldta[elem_id]; + + if ( h5tpriv_elem_is_on_cur_level ( f, el ) == H5_NOK ) { + continue; + } + TRY ( h5priv_search_idlist ( f, *list, elem_id ) ); + } + return H5_SUCCESS; +} + +static h5_err_t +get_triangles_upadjacent_to_edge ( + h5_file_t * const f, + const h5_id_t entity_id, + h5_idlist_t **list + ) { + h5_idlist_t *children; + TRY ( h5priv_alloc_idlist ( f, &children, 8 ) ); + TRY ( compute_children_of_edge ( f, entity_id, children ) ); + TRY ( h5priv_alloc_idlist ( f, list, 8 ) ); + h5_id_t *edge_id = children->items; + h5_id_t *end = children->items+children->num_items; + for ( ; edge_id < end; edge_id++ ) { + h5_id_t elem_id = h5tpriv_get_elem_idx ( *edge_id ); + TRY ( h5priv_search_idlist ( f, *list, elem_id ) ); + } + TRY ( h5priv_free_idlist( f, &children ) ); + + return H5_SUCCESS; +} + +static h5_err_t +get_vertices_downadjacent_to_edge ( + h5_file_t * const f, + const h5_id_t entity_id, + h5_idlist_t **list + ) { + h5_idlist_t *children; + TRY ( h5priv_alloc_idlist ( f, &children, 8 ) ); + TRY( compute_sections_of_edge ( f, entity_id, children ) ); + TRY ( h5priv_alloc_idlist ( f, list, 8 ) ); + int i; + h5_id_t *edge_id = children->items; + for ( i = 0; i < children->num_items; i++, edge_id++ ) { + h5_id_t vids[2]; + TRY ( h5t_get_vertex_indices_of_edge ( f, *edge_id, vids ) ); + TRY ( h5priv_search_idlist ( f, *list, vids[0] ) ); + TRY ( h5priv_search_idlist ( f, *list, vids[1] ) ); + } + TRY ( h5priv_free_idlist( f, &children ) ); + return H5_SUCCESS; +} + +/* + Compute downward adjacent vertices of all edges of triangle. + */ +static h5_err_t +get_vertices_downadjacent_to_triangle ( + h5_file_t* const f, + const h5_id_t entity_id, + h5_idlist_t** list + ) { + h5_idlist_t* children; + TRY ( h5priv_alloc_idlist ( f, &children, 8 ) ); + + h5_id_t elem_idx = h5tpriv_get_elem_idx ( entity_id ); + // loop over all edges of triangle + h5_id_t face_idx; + h5_id_t num_faces = f->t->ref_element->num_faces[1]; + for (face_idx = 0; face_idx < num_faces; face_idx++) { + TRY( compute_sections_of_edge ( + f, + h5tpriv_build_edge_id (face_idx, elem_idx), + children) ); + } + TRY ( h5priv_alloc_idlist (f, list, 8) ); + h5_id_t *edge_idp = children->items; + int i; + for ( i = 0; i < children->num_items; i++, edge_idp++ ) { + h5_id_t vids[2]; + TRY ( h5t_get_vertex_indices_of_edge ( f, *edge_idp, vids ) ); + TRY ( h5priv_search_idlist ( f, *list, vids[0] ) ); + TRY ( h5priv_search_idlist ( f, *list, vids[1] ) ); + } + TRY ( h5priv_free_idlist(f, &children) ); + return H5_SUCCESS; +} + +static h5_err_t +get_edges_downadjacent_to_triangle ( + h5_file_t* const f, + const h5_id_t entity_id, + h5_idlist_t** list + ) { + h5_idlist_t* children; + TRY ( h5priv_alloc_idlist (f, &children, 8) ); + + h5_id_t elem_idx = h5tpriv_get_elem_idx (entity_id); + // loop over all edges of triangle + h5_id_t face_idx; + h5_id_t num_faces = f->t->ref_element->num_faces[1]; + for ( face_idx = 0; face_idx < num_faces; face_idx++ ) { + TRY( compute_sections_of_edge ( + f, + h5tpriv_build_edge_id (face_idx, elem_idx), + children) ); + } + TRY ( h5priv_alloc_idlist (f, list, 8) ); + h5_id_t *edge_idp = children->items; + int i; + for (i = 0; i < children->num_items; i++, edge_idp++) { + TRY ( h5priv_search_idlist (f, *list, *edge_idp) ); + } + TRY ( h5priv_free_idlist(f, &children) ); + return H5_SUCCESS; +} + + +static h5_err_t +release_internal_structs ( + h5_file_t * const f + ) { + /* TO BE WRITTEN @@@ */ + return H5_SUCCESS; +} + +struct h5t_adjacency_methods h5tpriv_trim_adjacency_methods = { + rebuild_internal_structs, + release_internal_structs, + get_edges_upadjacent_to_vertex, + get_triangles_upadjacent_to_vertex, + NULL, + get_triangles_upadjacent_to_edge, + NULL, + NULL, + get_vertices_downadjacent_to_edge, + get_vertices_downadjacent_to_triangle, + NULL, + get_edges_downadjacent_to_triangle, + NULL, + NULL +}; + diff --git a/src/h5core/h5t_trim_adjacencies_private.h b/src/h5core/h5t_trim_adjacencies_private.h new file mode 100644 index 0000000..bf06248 --- /dev/null +++ b/src/h5core/h5t_trim_adjacencies_private.h @@ -0,0 +1,5 @@ +#ifndef __H5T_TRIM_ADJACENCIES_PRIVATE_H +#define __H5T_TRIM_ADJACENCIES_PRIVATE_H + +extern struct h5t_adjacency_methods h5tpriv_trim_adjacency_methods; +#endif diff --git a/src/h5core/h5t_trim_retrieve.c b/src/h5core/h5t_trim_retrieve.c new file mode 100644 index 0000000..36b4bd1 --- /dev/null +++ b/src/h5core/h5t_trim_retrieve.c @@ -0,0 +1,34 @@ +#include +#include + +#include "h5core/h5_core.h" +#include "h5_core_private.h" + +static h5_err_t +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_element = f->t->ref_element; + switch (iter->ref_element->dim - codim) { + case 0: // iterate vertices + iter->find = h5tpriv_find_tv2; + return h5tpriv_skip_to_next_elem_on_level (f, iter); + case 1: // iterate edges + iter->find = h5tpriv_find_te2; + return h5tpriv_skip_to_next_elem_on_level (f, iter); + case 2: // iterate elems + iter->find = NULL; + return H5_SUCCESS; + default: + return h5_error_internal (f, __FILE__, __func__, __LINE__); + } +} + +struct h5t_retrieve_methods h5tpriv_trim_retrieve_methods = { + begin_iterate_entities +}; diff --git a/src/h5core/h5t_trim_retrieve_private.h b/src/h5core/h5t_trim_retrieve_private.h new file mode 100644 index 0000000..aeaa838 --- /dev/null +++ b/src/h5core/h5t_trim_retrieve_private.h @@ -0,0 +1,5 @@ +#ifndef __H5T_TRIM_RETRIEVE_PRIVATE_H +#define __H5T_TRIM_RETRIEVE_PRIVATE_H + +extern struct h5t_retrieve_methods h5tpriv_trim_retrieve_methods; +#endif diff --git a/src/h5core/h5t_trim_store.c b/src/h5core/h5t_trim_store.c new file mode 100644 index 0000000..c670bb3 --- /dev/null +++ b/src/h5core/h5t_trim_store.c @@ -0,0 +1,210 @@ +#include +#include + +#include "h5core/h5_core.h" +#include "h5_core_private.h" + +static h5_err_t +alloc_triangles ( + h5_file_t * const f, + const size_t cur, + const size_t new + ) { + h5t_fdata_t *t = f->t; + const size_t nvertices = 3; + + /* alloc mem for elements */ + TRY ( t->elems.tris = h5priv_alloc ( + f, + t->elems.tris, + new * sizeof(t->elems.tris[0]) ) ); + memset ( + t->elems.tris + cur, + -1, + (new-cur) * sizeof(t->elems.tris[0]) ); + + /* alloc mem for local data of elements */ + TRY ( t->elems_ldta = h5priv_alloc ( + f, + t->elems_ldta, + new * sizeof (t->elems_ldta[0]) ) ); + memset ( + t->elems_ldta + cur, + -1, + (new-cur) * sizeof (t->elems_ldta[0]) ); + + /* alloc mem for local vertex IDs of elements */ + TRY ( t->elems_lvids = h5priv_alloc ( + f, + t->elems_lvids, + new * sizeof(t->elems_lvids[0]) * nvertices ) ); + memset ( + t->elems_lvids + cur * sizeof(t->elems_lvids[0]) * nvertices, + -1, + (new-cur) * sizeof(t->elems_lvids[0]) * nvertices ); + + /* re-init pointer to local vertex id in local data structure */ + h5_id_t *p = t->elems_lvids; + h5_id_t id; + for ( id = 0; id < new; id++, p+=nvertices ) { + t->elems_ldta[id].local_vids = p; + } + + /* alloc mem for global to local ID mapping */ + TRY ( h5priv_alloc_idmap ( f, &t->map_elem_g2l, new ) ); + + return H5_SUCCESS; +} + +static h5_err_t +get_direct_children_of_edge ( + h5_file_t * const f, + const h5_id_t face_idx, + const h5_id_t elem_idx, // index of the first child! + h5_id_t children[2] + ) { + int num_faces = f->t->ref_element->num_faces[1]; + /* + Please note: The face index of the children and the father is + always the same. The only think we have to know, is the offset + to the element index of the first child. This is either 0, 1 or + 2. The third child is an "inner" child which doesn't superpose edges + of the parent. + + The direct children of edge 0 of an element are edge 0 of the + first child and edge 0 of the second child, giving the offset 0 + and 1 for this edge. + */ + int off[3][2] = { {0,1}, // edge 0 + {0,2}, // edge 1 + {1,2} // edge 2 + }; + + if ( ( face_idx < 0 ) || ( face_idx >= num_faces ) ) { + return h5_error_internal ( f, __FILE__, __func__, __LINE__ ); + } + children[0] = h5tpriv_build_edge_id ( face_idx, elem_idx+off[face_idx][0] ); + children[1] = h5tpriv_build_edge_id ( face_idx, elem_idx+off[face_idx][1] ); + return H5_SUCCESS; +} + +static h5_id_t +bisect_edge ( + h5_file_t * const f, + h5_id_t face_idx, + h5_id_t elem_idx + ) { + h5t_fdata_t *t = f->t; + h5_idlist_t *retval; + /* + get all elements sharing the given edge + */ + TRY( h5tpriv_find_te2 (f, face_idx, elem_idx, &retval) ); + /* + check wether one of the found elements has been refined + */ + size_t i; + for ( i = 0; i < retval->num_items; i++ ) { + h5_id_t local_id = h5tpriv_get_elem_idx ( retval->items[i] ); + h5_elem_ldta_t *tet = &t->elems_ldta[local_id]; + if ( tet->local_child_eid >= 0 ) { + /* + this element has been refined! + return bisecting point + */ + h5_id_t face_id = h5tpriv_get_face_idx ( + retval->items[i] ); + h5_id_t kids[2], edge0[2], edge1[2]; + TRY ( get_direct_children_of_edge ( + f, + face_id, + tet->local_child_eid, + kids ) ); + TRY ( h5t_get_vertex_indices_of_edge ( f, kids[0], edge0 ) ); + TRY ( h5t_get_vertex_indices_of_edge ( f, kids[1], edge1 ) ); + if ( (edge0[0] == edge1[0]) || (edge0[0] == edge1[1]) ) + return edge0[0]; + else + return edge0[1]; + } + } + /* + None of the elements has been refined -> add new vertex. + */ + h5_id_t vindices[2]; + TRY( h5t_get_vertex_indices_of_edge2 (f, face_idx, elem_idx, vindices) ); + h5_float64_t *P0 = t->vertices[vindices[0]].P; + h5_float64_t *P1 = t->vertices[vindices[1]].P; + h5_float64_t P[3]; + + P[0] = ( P0[0] + P1[0] ) / 2.0; + P[1] = ( P0[1] + P1[1] ) / 2.0; + P[2] = ( P0[2] + P1[2] ) / 2.0; + + return h5t_store_vertex ( f, -1, P ); +} + +/*! + Refine triangle \c local_eid + + \return Local id of first new triangle or \c -1 +*/ +static h5_id_t +refine_triangle ( + h5_file_t* const f, + const h5_id_t elem_idx + ) { + h5t_fdata_t *t = f->t; + h5_id_t vertices[6]; // local vertex indices + h5_id_t elem_idx_of_first_child; + h5_elem_ldta_t *el = &t->elems_ldta[elem_idx]; + + if ( el->local_child_eid >= 0 ) + return h5_error ( + f, + H5_ERR_INVAL, + "Element %lld already refined.", + elem_idx ); + + vertices[0] = el->local_vids[0]; + vertices[1] = el->local_vids[1]; + vertices[2] = el->local_vids[2]; + + vertices[3] = bisect_edge( f, 0, elem_idx ); + vertices[4] = bisect_edge( f, 1, elem_idx ); + vertices[5] = bisect_edge( f, 2, elem_idx ); + + h5_id_t new_elem[3]; + + new_elem[0] = vertices[0]; // V[0] < V[3] , V[4] + new_elem[1] = vertices[3]; + new_elem[2] = vertices[4]; + TRY( elem_idx_of_first_child = h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[3]; // V[3] < V[1] , V[5] + new_elem[1] = vertices[1]; + new_elem[2] = vertices[5]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[4]; // V[4] < V[5] , V[2] + new_elem[1] = vertices[5]; + new_elem[2] = vertices[2]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + new_elem[0] = vertices[3]; // V[3] < V[4] , V[5] + new_elem[1] = vertices[4]; + new_elem[2] = vertices[5]; + TRY( h5t_store_elem (f, elem_idx, new_elem) ); + + t->elems.tris[elem_idx].global_child_eid = elem_idx_of_first_child; + t->elems_ldta[elem_idx].local_child_eid = elem_idx_of_first_child; + t->num_elems_on_level[t->cur_level]--; + + return elem_idx_of_first_child; +} + +struct h5t_store_methods h5tpriv_trim_store_methods = { + alloc_triangles, + refine_triangle, + get_direct_children_of_edge +}; diff --git a/src/h5core/h5t_trim_store_private.h b/src/h5core/h5t_trim_store_private.h new file mode 100644 index 0000000..c51f1c8 --- /dev/null +++ b/src/h5core/h5t_trim_store_private.h @@ -0,0 +1,5 @@ +#ifndef __H5T_TRIM_STORE_PRIVATE_H +#define __H5T_TRIM_STORE_PRIVATE_H + +extern struct h5t_store_methods h5tpriv_trim_store_methods; +#endif