- files added

- implement traversal function with iterators
This commit is contained in:
2010-05-12 20:32:49 +00:00
parent 538b0989eb
commit c6cd77b523
11 changed files with 946 additions and 0 deletions
+10
View File
@@ -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
+38
View File
@@ -0,0 +1,38 @@
#include <stdlib.h>
#include <string.h>
#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
};
+5
View File
@@ -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
+246
View File
@@ -0,0 +1,246 @@
#include <stdlib.h>
#include <string.h>
#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
};
+5
View File
@@ -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
+383
View File
@@ -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 <time.h>
#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
};
@@ -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
+34
View File
@@ -0,0 +1,34 @@
#include <stdlib.h>
#include <string.h>
#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
};
+5
View File
@@ -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
+210
View File
@@ -0,0 +1,210 @@
#include <stdlib.h>
#include <string.h>
#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
};
+5
View File
@@ -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