h5t_store_trim.c 8.52 KB
Newer Older
1
/*
gsell's avatar
gsell committed
2
  Copyright (c) 2006-2016, The Regents of the University of California,
3 4 5 6 7 8 9
  through Lawrence Berkeley National Laboratory (subject to receipt of any
  required approvals from the U.S. Dept. of Energy) and the Paul Scherrer
  Institut (Switzerland).  All rights reserved.

  License: see file COPYING in top level of source distribution.
*/

gsell's avatar
gsell committed
10
#include "private/h5t_types.h"
11
#include "private/h5t_err.h"
gsell's avatar
gsell committed
12 13 14 15 16 17 18
#include "private/h5t_access.h"
#include "private/h5t_adjacencies.h"
#include "private/h5t_core.h"
#include "private/h5t_map.h"
#include "private/h5t_model.h"
#include "private/h5t_io.h"
#include "private/h5t_store.h"
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

#include "h5core/h5t_map.h"

static h5_err_t
alloc_loc_elems (
        h5t_mesh_t* const m,
        const size_t cur,
        const size_t new
        ) {
	H5_PRIV_FUNC_ENTER (h5_err_t, "m=%p, cur=%zu, new=%zu", m, cur, new);

	/* alloc mem for local data of elements */
	TRY (m->loc_elems = h5_alloc (
	             m->loc_elems,
	             new * sizeof (h5_loc_tri_t)));
	memset (
	        (h5_loc_tri_t*)m->loc_elems + cur,
	        -1,
	        (new-cur) * sizeof (h5_loc_tri_t));

gsell's avatar
gsell committed
39
	H5_RETURN (H5_SUCCESS);
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
}

/*
   Bisect edge and return local vertex index of the bisecting point.
 */
static h5_loc_idx_t
bisect_edge (
        h5t_mesh_t* const m,
        const h5_loc_idx_t face_idx,
        const h5_loc_idx_t elem_idx
        ) {
	H5_PRIV_FUNC_ENTER (h5_loc_idx_t,
	                    "m=%p, face_idx=%lld, elem_idx=%lld",
	                    m, (long long)face_idx, (long long)elem_idx);
	h5_loc_idlist_t* retval;
	// get all elements sharing the given edge
	TRY (h5tpriv_find_te2 (m, face_idx, elem_idx, &retval));
	// check weather one of the found elements has been refined
	size_t i;
	for (i = 0; i < retval->num_items; i++) {
//		// check if it is shared boundary edge that has been refined already
//		if (num_b_edges > 0) {
//			h5_glb_idx_t my_glb_idx = m->loc_elems[elem_idx].glb_idx;
//			for (int j = 0; j < num_b_edges; j++) {
//				if (b_edges[j].idx == my_glb_idx &&
//						b_edges[j].face_idx == face_idx) {
gsell's avatar
gsell committed
66
//					H5_LEAVE (b_edges[j].vtx);
67 68 69 70 71 72 73 74 75 76 77
//				}
//			}
//		}
		h5_loc_id_t kids[2] = {-1,-1};
		TRY (h5tpriv_get_loc_entity_children (m, retval->items[i], kids));
		if (kids[0] >= 0) {
			// element has been refined, return bisecting point
			h5_loc_idx_t edge0[2], edge1[2];
			TRY (h5t_get_loc_vertex_indices_of_edge (m, kids[0], edge0));
			TRY (h5t_get_loc_vertex_indices_of_edge (m, kids[1], edge1));
			if ((edge0[0] == edge1[0]) || (edge0[0] == edge1[1])) {
gsell's avatar
gsell committed
78
				H5_LEAVE (edge0[0]);
79
			} else {
gsell's avatar
gsell committed
80
				H5_LEAVE (edge0[1]);
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
			}
		}
	}
	/*
	   None of the elements has been refined -> add new vertex.
	 */
	h5_loc_idx_t indices[2];
	TRY( h5t_get_loc_vertex_indices_of_edge2 (m, face_idx, elem_idx, indices) );
	h5_float64_t* P0 = m->vertices[indices[0]].P;
	h5_float64_t* P1 = m->vertices[indices[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;

gsell's avatar
gsell committed
97
	H5_RETURN (h5t_store_vertex (m, -1, P));
98 99 100 101 102 103 104 105 106 107 108 109 110
}

/*
   Please read note about number of new vertices in tetrahedral
   mesh implementation.
 */
static h5_err_t
pre_refine_triangle (
        h5t_mesh_t* const m
        ) {
	H5_PRIV_FUNC_ENTER (h5_err_t, "m=%p", m);
	unsigned int num_interior_elems_to_refine = m->marked_entities->num_items;
	TRY (h5t_begin_store_vertices (m, num_interior_elems_to_refine*3 + 64));
111
	TRY (h5t_begin_store_elems (m, num_interior_elems_to_refine*4));
gsell's avatar
gsell committed
112
	H5_RETURN (H5_SUCCESS);
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
}

/*!
   Refine triangle \c local_eid

   \return Local index of first new triangle or \c H5_ERR
 */
static h5_loc_idx_t
refine_triangle (
        h5t_mesh_t* const m,
        const h5_loc_idx_t elem_idx
        ) {
	H5_PRIV_FUNC_ENTER (h5_loc_idx_t,
	                    "m=%p, elem_idx=%lld",
	                    m, (long long)elem_idx);
	h5_loc_idx_t vertices[6];       // local vertex indices
	h5_loc_idx_t elem_idx_of_first_child;
	h5_loc_tri_t* el = (h5_loc_tri_t*)m->loc_elems + elem_idx;

	if (el->child_idx >= 0)
133 134 135 136
		H5_RETURN_ERROR (
			H5_ERR_INVAL,
			"Element %lld already refined.",
			(long long)elem_idx);
137 138 139 140 141 142 143 144 145 146 147 148 149 150

	vertices[0] = el->vertex_indices[0];
	vertices[1] = el->vertex_indices[1];
	vertices[2] = el->vertex_indices[2];

	vertices[3] = bisect_edge (m, 0, elem_idx);
	vertices[4] = bisect_edge (m, 1, elem_idx);
	vertices[5] = bisect_edge (m, 2, elem_idx);

	h5_loc_idx_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];
151
	TRY (elem_idx_of_first_child = h5tpriv_add_cell (m, elem_idx, new_elem, NULL));
152 153 154 155

	new_elem[0] = vertices[3];  // V[3] < V[1] , V[5]
	new_elem[1] = vertices[1];
	new_elem[2] = vertices[5];
156
	TRY (h5tpriv_add_cell (m, elem_idx, new_elem, NULL));
157 158 159 160

	new_elem[0] = vertices[4];  // V[4] < V[5] , V[2]
	new_elem[1] = vertices[5];
	new_elem[2] = vertices[2];
161
	TRY (h5tpriv_add_cell (m, elem_idx, new_elem, NULL));
162 163 164 165

	new_elem[0] = vertices[3];  // V[3] < V[4] , V[5]
	new_elem[1] = vertices[5];
	new_elem[2] = vertices[4]; // TODO check if that ordering is correct!
166
	TRY (h5tpriv_add_cell (m, elem_idx, new_elem, NULL));
167 168 169 170

	((h5_loc_tri_t*)m->loc_elems)[elem_idx].child_idx = elem_idx_of_first_child;
	m->num_interior_leaf_elems[m->leaf_level]--;

gsell's avatar
gsell committed
171
	H5_RETURN (elem_idx_of_first_child);
172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
}

static inline h5_loc_idx_t
compute_neighbor_of_face (
        h5t_mesh_t* const m,
        h5_loc_idx_t elem_idx,
        const h5_loc_idx_t face_idx
        ) {
	H5_PRIV_FUNC_ENTER (h5_loc_idx_t
	                    , "m=%p, elem_idx=%lld, face_idx=%lld",
	                    m, (long long)elem_idx, (long long)face_idx);
	h5_loc_idlist_t* te;
	h5_loc_idx_t neighbor_idx = -2;

	do {
		TRY( h5tpriv_find_te2 (
		             m,
		             face_idx,
		             elem_idx,
		             &te) );
		if (te == NULL) {
gsell's avatar
gsell committed
193
			H5_LEAVE (h5_error_internal ());
194 195 196 197 198 199 200 201 202 203 204
		}
		if (te->num_items == 1) {
			h5_loc_idx_t old_elem_idx = elem_idx;
			// neighbor is coarser or face is on the boundary
			elem_idx = ((h5_loc_tri_t*)m->loc_elems)[elem_idx].parent_idx;
			if (elem_idx == -1) {
				// we are on the level of the macro grid
				neighbor_idx = -1;
			}
			if (elem_idx < -1) { // this should only happen if we are on the boarder
				// of a loaded chunk and the parent is on a different chunk
205
				if (h5_debug_mask >= 6) {
206 207 208
					h5_debug ("Elem %d is on different proc than its parent %d \n"
						"therefore neighborhood idx is not correct resolved", old_elem_idx, elem_idx);
				}
gsell's avatar
gsell committed
209 210
				assert (m->f->myproc != h5priv_find_proc_to_write (m, old_elem_idx));
				H5_LEAVE (~0); //TODO what is a resonable output here?
211 212 213 214 215 216 217 218 219 220 221
			}
		} else if (te->num_items == 2) {
			// neighbor has same level of coarsness
			if (h5tpriv_get_elem_idx(te->items[0]) == elem_idx) {
				neighbor_idx = h5tpriv_get_elem_idx (te->items[1]);
			} else {
				neighbor_idx = h5tpriv_get_elem_idx (te->items[0]);
			}

		} else {
			printf ("elem %d face %d num_items %d", elem_idx, face_idx, te->num_items);
gsell's avatar
gsell committed
222
			H5_LEAVE (h5_error_internal ());
223 224
		}
	} while (neighbor_idx < -1);
gsell's avatar
gsell committed
225
	H5_RETURN (neighbor_idx);
226 227 228 229 230 231 232 233 234 235 236 237
}

/*
   Compute neighbors for elements on given level.
 */
static inline h5_err_t
compute_neighbors_of_elems (
        h5t_mesh_t* const m,
        h5_lvl_idx_t level
        ) {
	H5_PRIV_FUNC_ENTER (h5_err_t, "m=%p, level=%d", m, level);
	if (level < 0 || level >= m->num_leaf_levels) {
238 239 240 241 242 243
		H5_RETURN_ERROR (
			H5_ERR_INVAL,
			"level idx %lld out of bound, must be in [%lld,%lld]",
			(long long)level,
			(long long)0,
			(long long)m->num_leaf_levels);
244 245 246 247 248 249 250 251 252 253 254 255 256 257
	}
	h5_loc_idx_t elem_idx = level == 0 ? 0 : m->num_interior_elems[level-1];
	const h5_loc_idx_t last_idx = m->num_interior_elems[level] - 1;
	h5_loc_tri_t *el = (h5_loc_tri_t*)m->loc_elems + elem_idx;
	while (elem_idx <= last_idx) {
		h5_loc_idx_t face_idx = 0;
		for (; face_idx < 3; face_idx++) {
			el->neighbor_indices[face_idx] =
			        compute_neighbor_of_face (m, elem_idx, face_idx);
		}
		elem_idx++;
		el++;
	}

gsell's avatar
gsell committed
258
	H5_RETURN (H5_SUCCESS);
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
}
/*
 * returns number of newly created triangles when refining a triangle
 */
static int
get_num_new_triangles (
		void
		) {
	return 4;
}
static h5_err_t
end_store_elems (
        h5t_mesh_t* const m
        ) {
	H5_PRIV_FUNC_ENTER (h5_err_t, "m=%p", m);

	h5_loc_idx_t start_idx = (m->leaf_level > 0) ? m->num_interior_elems[m->leaf_level-1] : 0;
	h5_loc_idx_t count = m->num_interior_elems[m->leaf_level] - start_idx;

	TRY( h5tpriv_update_internal_structs (m, m->leaf_level) );
	TRY( compute_neighbors_of_elems (m, m->leaf_level) );
	TRY( h5tpriv_init_elem_flags (m, start_idx, count) );
gsell's avatar
gsell committed
281
	H5_RETURN (H5_SUCCESS);
282 283 284 285 286 287 288 289 290
}

struct h5t_store_methods h5tpriv_trim_store_methods = {
	alloc_loc_elems,
	pre_refine_triangle,
	refine_triangle,
	get_num_new_triangles,
	end_store_elems,
};