Branch data Line data Source code
1 : : #ifndef AMR_refinement_h
2 : : #define AMR_refinement_h
3 : :
4 : : #include <algorithm>
5 : :
6 : : #include "Macro.hpp"
7 : : #include "tet_store.hpp"
8 : : #include "node_connectivity.hpp"
9 : :
10 : : // TODO: make this have a base class to support multiple generator schemes
11 : : // using the policy design pattern
12 : :
13 : : #if defined(STRICT_GNUC)
14 : : #pragma GCC diagnostic push
15 : : #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
16 : : #endif
17 : :
18 : : namespace AMR {
19 : :
20 [ + - ]: 2128 : class refinement_t {
21 : : private:
22 : :
23 : : size_t DEFAULT_REFINEMENT_LEVEL = 0; //TODO: Is this in the right place?
24 : : size_t MIN_REFINEMENT_LEVEL = DEFAULT_REFINEMENT_LEVEL;
25 : : // list of "intermediate" edges to be deleted
26 : : std::vector< edge_t > delete_list;
27 : :
28 : : public:
29 : :
30 : : size_t MAX_REFINEMENT_LEVEL = 3;
31 : :
32 : : // TODO: Document this
33 : : child_id_list_t generate_child_ids( tet_store_t& tet_store, size_t parent_id, size_t count = MAX_CHILDREN)
34 : : {
35 : : //return morton_id_generator_t::get_children_ids(parent_id);
36 : : return tet_store.generate_child_ids(parent_id, count);
37 : : }
38 : :
39 : : /**
40 : : * @brief function to detect when an invalid refinement is
41 : : * invoked
42 : : *
43 : : * @param tet_store Tet store to use
44 : : * @param tet_id Id the of the tet which will be refined
45 : : *
46 : : * @return A bool stating if the tet can be validly refined
47 : : */
48 : : bool check_allowed_refinement( tet_store_t& tet_store, size_t tet_id)
49 : : {
50 : : Refinement_State& master_element = tet_store.data(tet_id);
51 : :
52 : : // These asserts mean we never actually try refine a 1:2 or 1:4
53 : : assert( master_element.refinement_case !=
54 : : Refinement_Case::one_to_two);
55 : : assert( master_element.refinement_case !=
56 : : Refinement_Case::one_to_four);
57 : :
58 : : // cppcheck-suppress assertWithSideEffect
59 : : assert( tet_store.is_active(tet_id) );
60 : :
61 : : // Check this won't take us past the max refinement level
62 [ + - ][ + - ]: 99847 : if (master_element.refinement_level >= MAX_REFINEMENT_LEVEL)
[ + - ]
63 : : {
64 : : return false;
65 : : }
66 : :
67 : : // If we got here, we didn't detect anything which tells us not
68 : : // to refine
69 : : return true;
70 : : }
71 : :
72 : : /**
73 : : * @brief Method which takes a tet id, and deduces the other
74 : : * parameters needed to perform a 1:2
75 : : *
76 : : * @param tet_store Tet store to use
77 : : * @param node_connectivity Mesh node connectivity (graph)
78 : : * @param tet_id The id to refine 1:2
79 : : */
80 : 511 : void refine_one_to_two( tet_store_t& tet_store, node_connectivity_t& node_connectivity, size_t tet_id)
81 : : {
82 : 511 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
83 : 511 : node_pair_t nodes = find_single_refinement_nodes(tet_store,edge_list);
84 : 511 : refine_one_to_two( tet_store, node_connectivity, tet_id, nodes[0], nodes[1]);
85 : 511 : }
86 : :
87 : : /*
88 : : //! @brief Method which takes a tet id, and transforms arguments
89 : : //! into the form needed for the main 1:2 refinement method
90 : : //! @param tet_id The id to refine 1:2
91 : : void refine_one_to_two(
92 : : size_t tet_id,
93 : : std::string edge_key
94 : : )
95 : : {
96 : : std::vector<std::string> nodes = util::split(edge_key,KEY_DELIM);
97 : : size_t edge_node_A_id = std::stoul (nodes[0],nullptr,0);
98 : : size_t edge_node_B_id = std::stoul (nodes[1],nullptr,0);
99 : : refine_one_to_two( tet_id, edge_node_A_id, edge_node_B_id);
100 : : }
101 : : */
102 : : /**
103 : : * @brief Refine a given tet id into 2 children.
104 : : * NOTE: Does not do any validity checking (currently?)
105 : : *
106 : : * @param tet_store Tet store to use
107 : : * @param node_connectivity Mesh node connectivity (graph)
108 : : * @param tet_id Id of tet to refine
109 : : * @param edge_node_A_id The first node of id of the edge which
110 : : * will be split
111 : : * @param edge_node_B_id The second node of id of the
112 : : * edge which will be split
113 : : */
114 : 511 : void refine_one_to_two(
115 : : tet_store_t& tet_store,
116 : : node_connectivity_t& node_connectivity,
117 : : size_t tet_id,
118 : : size_t edge_node_A_id,
119 : : size_t edge_node_B_id
120 : : )
121 : : {
122 : :
123 : : trace_out << "refine_one_to_two" << std::endl;
124 : 0 : if (!check_allowed_refinement(tet_store,tet_id)) return;
125 : :
126 : 511 : tet_t original_tet = tet_store.get(tet_id);
127 : :
128 : : //coordinate_t original_tet_c = node_connectivity->id_to_coordinate(id);
129 : :
130 : 511 : size_t new_node_id = node_connectivity.add( edge_node_A_id, edge_node_B_id );
131 : :
132 : : /// Split existing tet into two new tets
133 : :
134 : : // The two new tets will be the same, but for each an edge will
135 : : // be cut, losing an edge replaced by E
136 : :
137 : : tet_t new_tet1;
138 : : tet_t new_tet2;
139 : :
140 : : // Create a new tet that is based on the original
141 : : copy_tet(&new_tet1, &original_tet);
142 : :
143 : : // Replace all node ids in tet that were pointing to A with new_node_id
144 : : replace_node(&new_tet1, edge_node_A_id, new_node_id);
145 : :
146 : : // Create a new tet that is based on the original
147 : : copy_tet(&new_tet2, &original_tet);
148 : :
149 : : // Replace all node ids in tet that were pointing to B with new_node_id
150 : : replace_node(&new_tet2, edge_node_B_id, new_node_id);
151 : :
152 : : // Now, update the edge list
153 : :
154 : : // Generate edges for split
155 : 511 : tet_store.edge_store.split(edge_node_A_id, edge_node_B_id, new_node_id,
156 : : Edge_Lock_Case::intermediate);
157 : :
158 : : child_id_list_t child_list = generate_child_ids(tet_store,tet_id, 2);
159 : :
160 : 511 : size_t first_child_id = child_list[0];
161 : 511 : size_t second_child_id = child_list[1];
162 : :
163 : : // Add the two new tets to the system
164 : : size_t new_tet_id = first_child_id;
165 [ + - ]: 511 : tet_store.add(
166 : : first_child_id,
167 : : new_tet1,
168 : : Refinement_Case::one_to_two,
169 : : tet_id
170 : : );
171 : :
172 : : //size_t new_tet_id2 = second_child_id;
173 [ + - ]: 511 : tet_store.add(
174 : : second_child_id,
175 : : new_tet2,
176 : : Refinement_Case::one_to_two,
177 : : tet_id
178 : : );
179 : :
180 : : //trace_out << "1:2 DOING REFINE OF " << tet_id << ". Adding " << child_list[0] << " and " << child_list[1] << std::endl;
181 : :
182 : : // This call is only needed to add a single edge, from the new
183 : : // node to the node on the normal to that face, but avoids
184 : : // directly calculating which nodes that is
185 [ + - ]: 511 : tet_store.generate_edges(new_tet_id);
186 : :
187 : : // Currently we lock one per tet, around the split node. We
188 : : // also need to lock the two "arms" which come out from it
189 : : //lock_edges_from_node(new_tet_id, new_node_id, Edge_Lock_Case::intermediate);
190 : : //lock_edges_from_node(new_tet_id2, new_node_id, Edge_Lock_Case::intermediate);
191 : :
192 : : // Deactivate parent tet?
193 : : tet_store.deactivate(tet_id);
194 : : //lock_edges_from_node(new_node_id, Edge_Lock_Case::intermediate);
195 : : trace_out << "Adding " << new_node_id << " to intermediate list " << std::endl;
196 : : tet_store.intermediate_list.insert(new_node_id);
197 : : }
198 : :
199 : : /**
200 : : * @brief Method which takes a tet id, and deduces the other
201 : : * parameters needed to perform a 1:4
202 : : *
203 : : * @param tet_store Tet store to use
204 : : * @param node_connectivity Mesh node connectivity (graph)
205 : : * @param tet_id The id to refine 1:4
206 : : */
207 : 706 : void refine_one_to_four( tet_store_t& tet_store,
208 : : node_connectivity_t& node_connectivity, size_t tet_id)
209 : : {
210 : : trace_out << "do refine 1:4 " << std::endl;
211 : : //bool face_refine = false;
212 : : size_t face_refine_id = 0; // FIXME: Does this need a better default
213 : 706 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
214 : :
215 : : // Iterate over each face
216 [ + - ]: 1811 : for (size_t face = 0; face < NUM_TET_FACES; face++)
217 : : {
218 : : int num_face_refine_edges = 0;
219 : :
220 : 1811 : face_ids_t face_ids = face_list[face];
221 : : trace_out << "face ids " <<
222 : : face_ids[0] << ", " <<
223 : : face_ids[1] << ", " <<
224 : : face_ids[2] << ", " <<
225 : : std::endl;
226 : :
227 : 1811 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
228 : : // For this face list, see which ones need refining
229 : : trace_out << "Looping to " << NUM_FACE_NODES << std::endl;
230 [ + + ]: 5754 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
231 : : {
232 : : trace_out << "nodes " << k << std::endl;
233 : :
234 : 4531 : edge_t edge = face_edge_list[k];
235 [ + + ]: 4531 : if (tet_store.edge_store.get(edge).needs_refining == 1)
236 : : {
237 : 2909 : num_face_refine_edges++;
238 : : trace_out << "Ref " << edge << " Num face => " << num_face_refine_edges << std::endl;
239 : : }
240 : :
241 : : // Check for locked edges
242 : : // This case only cares about faces with no locks
243 [ + + ]: 4531 : if (tet_store.edge_store.lock_case(edge) != Edge_Lock_Case::unlocked)
244 : : {
245 : : // Abort this face
246 : : trace_out << "Face has lock it's not this one " << face << std::endl;
247 : : num_face_refine_edges = 0;
248 : 588 : break;
249 : : }
250 : : trace_out << "Num face => " << num_face_refine_edges << std::endl;
251 : : }
252 [ + + ]: 1223 : if (num_face_refine_edges >= 2)
253 : : {
254 : : assert(num_face_refine_edges < 4);
255 : : //face_refine = true;
256 : : trace_out << "Accepting face " << face << std::endl;
257 : : face_refine_id = face;
258 : 706 : break;
259 : : }
260 : : }
261 : :
262 : 706 : tet_t tet = tet_store.get(tet_id);
263 : : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list, face_refine_id);
264 : 706 : size_t opposite_id = tet[opposite_offset];
265 : :
266 : : trace_out << "1:4 tet mark id " << tet_id << std::endl;
267 : : trace_out << "opposite offset " << opposite_offset << std::endl;
268 : : trace_out << "opposite id " << opposite_id << std::endl;
269 : : trace_out << "face refine id " << face_refine_id << std::endl;
270 : : trace_out << "face list 0 " << face_list[face_refine_id][0] << std::endl;
271 : : trace_out << "face list 1 " << face_list[face_refine_id][1] << std::endl;
272 : : trace_out << "face list 2 " << face_list[face_refine_id][2] << std::endl;
273 : :
274 : 706 : refine_one_to_four(tet_store, node_connectivity, tet_id, face_list[face_refine_id], opposite_id);
275 : 706 : }
276 : :
277 : : /**
278 : : * @brief Method which takes a tet id, and deduces the other
279 : : * parameters needed to perform a 1:4, as a part of an 8:4 deref
280 : : *
281 : : * @param tet_store Tet store to use
282 : : * @param node_connectivity Mesh node connectivity (graph)
283 : : * @param tet_id The id to refine 1:4
284 : : */
285 : 0 : void deref_refine_one_to_four( tet_store_t& tet_store,
286 : : node_connectivity_t& node_connectivity, size_t tet_id)
287 : : {
288 : : trace_out << "do refine 1:4 " << std::endl;
289 : : //bool face_refine = false;
290 : : size_t face_refine_id = 0; // FIXME: Does this need a better default
291 : 0 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
292 : :
293 : : // Iterate over each face
294 [ - - ]: 0 : for (size_t face = 0; face < NUM_TET_FACES; face++)
295 : : {
296 : : int num_face_refine_edges = 0;
297 : :
298 : 0 : face_ids_t face_ids = face_list[face];
299 : : trace_out << "face ids " <<
300 : : face_ids[0] << ", " <<
301 : : face_ids[1] << ", " <<
302 : : face_ids[2] << ", " <<
303 : : std::endl;
304 : :
305 : 0 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
306 : : // For this face list, see which ones need refining
307 : : trace_out << "Looping to " << NUM_FACE_NODES << std::endl;
308 [ - - ]: 0 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
309 : : {
310 : : trace_out << "nodes " << k << std::endl;
311 : :
312 : 0 : edge_t edge = face_edge_list[k];
313 [ - - ]: 0 : if (tet_store.edge_store.get(edge).needs_refining == 2)
314 : : {
315 : 0 : num_face_refine_edges++;
316 : : trace_out << "Ref " << edge << " Num face => " << num_face_refine_edges << std::endl;
317 : : }
318 : :
319 : : // Check for locked edges
320 : : // This case only cares about faces with no locks
321 [ - - ]: 0 : if (tet_store.edge_store.lock_case(edge) != Edge_Lock_Case::unlocked)
322 : : {
323 : : // Abort this face
324 : : trace_out << "Face has lock it's not this one " << face << std::endl;
325 : : num_face_refine_edges = 0;
326 : 0 : break;
327 : : }
328 : : trace_out << "Num face => " << num_face_refine_edges << std::endl;
329 : : }
330 [ - - ]: 0 : if (num_face_refine_edges >= 2)
331 : : {
332 : : assert(num_face_refine_edges < 4);
333 : : //face_refine = true;
334 : : trace_out << "Accepting face " << face << std::endl;
335 : : face_refine_id = face;
336 : 0 : break;
337 : : }
338 : : }
339 : :
340 : 0 : tet_t tet = tet_store.get(tet_id);
341 : : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list, face_refine_id);
342 : 0 : size_t opposite_id = tet[opposite_offset];
343 : :
344 : : trace_out << "1:4 tet mark id " << tet_id << std::endl;
345 : : trace_out << "opposite offset " << opposite_offset << std::endl;
346 : : trace_out << "opposite id " << opposite_id << std::endl;
347 : : trace_out << "face refine id " << face_refine_id << std::endl;
348 : : trace_out << "face list 0 " << face_list[face_refine_id][0] << std::endl;
349 : : trace_out << "face list 1 " << face_list[face_refine_id][1] << std::endl;
350 : : trace_out << "face list 2 " << face_list[face_refine_id][2] << std::endl;
351 : :
352 : 0 : refine_one_to_four(tet_store, node_connectivity, tet_id, face_list[face_refine_id], opposite_id);
353 : 0 : }
354 : :
355 : : /**
356 : : * @brief Refine a given tet id into 4 children.
357 : : * NOTE: Does not do any validity checking (currently?)
358 : : *
359 : : * @param tet_store Tet store to use
360 : : * @param node_connectivity Mesh node connectivity (graph)
361 : : * @param tet_id The id of the tet to refine
362 : : * @param face_ids The ids which make the face to be split
363 : : * @param opposite_id The remaining id which is "opposite" the
364 : : * split face
365 : : */
366 : 706 : void refine_one_to_four(
367 : : tet_store_t& tet_store,
368 : : node_connectivity_t& node_connectivity,
369 : : size_t tet_id,
370 : : std::array<size_t, NUM_FACE_NODES> face_ids,
371 : : size_t opposite_id
372 : : )
373 : : {
374 : :
375 : : trace_out << "refine_one_to_four" << std::endl;
376 : 0 : if (!check_allowed_refinement(tet_store,tet_id)) return;
377 : :
378 : : trace_out << "Refining tet_id " << tet_id <<
379 : : " 1:4 opposite edge " << opposite_id << std::endl;
380 : :
381 : 706 : tet_t t = tet_store.get(tet_id);
382 : : trace_out << "Tet has nodes " <<
383 : : t[0] << ", " <<
384 : : t[1] << ", " <<
385 : : t[2] << ", " <<
386 : : t[3] << ", " <<
387 : : std::endl;
388 : :
389 : : trace_out << "face_ids " <<
390 : : face_ids[0] << ", " <<
391 : : face_ids[1] << ", " <<
392 : : face_ids[2] << ", " <<
393 : : std::endl;
394 : :
395 : 706 : size_t A = face_ids[0];
396 : 706 : size_t B = face_ids[1];
397 : 706 : size_t C = face_ids[2];
398 : : size_t D = opposite_id;
399 : :
400 : : trace_out <<
401 : : " A " << A <<
402 : : " B " << B <<
403 : : " C " << C <<
404 : : " D " << D <<
405 : : std::endl;
406 : :
407 : : // Make new nodes
408 : : //coordinate_t AB_mid = node_connectivity->find_mid_point(A, B);
409 : 706 : size_t AB = node_connectivity.add(A,B);
410 : :
411 : : //coordinate_t AC_mid = node_connectivity->find_mid_point(A, C);
412 : 706 : size_t AC = node_connectivity.add(A,C);
413 : :
414 : : //coordinate_t BC_mid = node_connectivity->find_mid_point(B, C);
415 : 706 : size_t BC = node_connectivity.add(B,C);
416 : :
417 : : // Use nodes to update edges
418 : : // All added edges will be locked due to containing intermediate points
419 : : // Split Outer face edges
420 : 706 : tet_store.edge_store.split(A, C, AC, Edge_Lock_Case::intermediate);
421 : 706 : tet_store.edge_store.split(A, B, AB, Edge_Lock_Case::intermediate);
422 : 706 : tet_store.edge_store.split(B, C, BC, Edge_Lock_Case::intermediate);
423 : :
424 : : // Connect D to intermediate points
425 : 706 : tet_store.edge_store.generate(D, AC, Edge_Lock_Case::intermediate);
426 : 706 : tet_store.edge_store.generate(D, BC, Edge_Lock_Case::intermediate);
427 : 706 : tet_store.edge_store.generate(D, AB, Edge_Lock_Case::intermediate);
428 : : // Connect inner edges
429 : 706 : tet_store.edge_store.generate(AC, BC, Edge_Lock_Case::intermediate);
430 : 706 : tet_store.edge_store.generate(AC, AB, Edge_Lock_Case::intermediate);
431 : 706 : tet_store.edge_store.generate(AB, BC, Edge_Lock_Case::intermediate);
432 : :
433 : : // Make new Tets
434 : : // This is just the node opposite the face plus each pair
435 : : // of the news nodes, and the old corner
436 : : // FIXME: How to find that near corner programatically?
437 : :
438 : : // Hard coded solution
439 : : // A AC AB D
440 : : // AC AB BC D
441 : : // AC BC C D
442 : : // AB B BC D
443 : :
444 : : size_t num_children = 4;
445 : : child_id_list_t child = generate_child_ids(tet_store,tet_id, num_children);
446 : :
447 : : // Outsides
448 [ + - ]: 706 : tet_store.add(child[0], A, AB, AC, D, Refinement_Case::one_to_four, tet_id);
449 [ + - ]: 706 : tet_store.add(child[2], AC, BC, C, D, Refinement_Case::one_to_four, tet_id);
450 [ + - ]: 706 : tet_store.add(child[3], AB, B, BC, D, Refinement_Case::one_to_four, tet_id);
451 : :
452 : : // Center
453 : 706 : size_t center_id = child[1]; // 1 to preserve Jacobian order
454 [ + - ]: 706 : tet_store.add(center_id, AC, AB, BC, D, Refinement_Case::one_to_four, tet_id);
455 : :
456 : :
457 : : // TODO: replace this with a more concise way to lock the correct edges
458 : :
459 : 706 : tet_store.add_center(center_id);
460 : : /*
461 : : lock_edges_from_node(child[0], AB, Edge_Lock_Case::intermediate);
462 : : lock_edges_from_node(child[0], AC, Edge_Lock_Case::intermediate);
463 : : lock_edges_from_node(child[2], AC, Edge_Lock_Case::intermediate);
464 : : lock_edges_from_node(child[2], BC, Edge_Lock_Case::intermediate);
465 : : lock_edges_from_node(child[3], AB, Edge_Lock_Case::intermediate);
466 : : lock_edges_from_node(child[3], BC, Edge_Lock_Case::intermediate);
467 : : lock_edges_from_node(center_id, AC, Edge_Lock_Case::intermediate);
468 : : lock_edges_from_node(center_id, AB, Edge_Lock_Case::intermediate);
469 : : lock_edges_from_node(center_id, BC, Edge_Lock_Case::intermediate);
470 : : */
471 : :
472 : :
473 : : tet_store.deactivate(tet_id);
474 : :
475 : : //trace_out << "1:4 DOING REFINE OF " << tet_id << ". Adding "
476 : : // << child[0] << ", "
477 : : // << child[1] << ", "
478 : : // << child[2] << ", "
479 : : // << child[3]
480 : : // << std::endl;
481 : :
482 : : /*
483 : : lock_edges_from_node(AB, Edge_Lock_Case::intermediate);
484 : : lock_edges_from_node(AC, Edge_Lock_Case::intermediate);
485 : : lock_edges_from_node(BC, Edge_Lock_Case::intermediate);
486 : : */
487 : :
488 : : trace_out << "Adding " << AB << " to intermediate list " << std::endl;
489 : : tet_store.intermediate_list.insert(AB);
490 : : trace_out << "Adding " << AC << " to intermediate list " << std::endl;
491 : : tet_store.intermediate_list.insert(AC);
492 : : trace_out << "Adding " << BC << " to intermediate list " << std::endl;
493 : : tet_store.intermediate_list.insert(BC);
494 : :
495 : : }
496 : :
497 : : /**
498 : : * @brief Refine a given tet id into 8 children.
499 : : * NOTE: Does not do any validity checking (currently?)
500 : : *
501 : : * @param tet_store Tet store to use
502 : : * @param node_connectivity Mesh node connectivity (graph)
503 : : * @param tet_id Id of tet to refine
504 : : */
505 : 98630 : void refine_one_to_eight( tet_store_t& tet_store,
506 : : node_connectivity_t& node_connectivity, size_t tet_id)
507 : : {
508 : :
509 : : trace_out << "refine_one_to_eight" << std::endl;
510 : 0 : if (!check_allowed_refinement(tet_store,tet_id)) return;
511 : :
512 : : // Split every edge into two
513 : : // Makes 4 tets out of the old corners and 3 near mid-points
514 : : // Make 4 out of the midpoints
515 : :
516 : : // For Tet {ABCD} need to know all (non-repeating) node pairs
517 : : // {AB} {AC} {AD} {BC} {BD} {CD}
518 : : // This can either be hard coded, or generated with a 2d loop
519 : : // The loop would just be i=0..4, j=i..4
520 : : //
521 : :
522 : :
523 : 98630 : tet_t tet = tet_store.get(tet_id);
524 : :
525 : : size_t A = tet[0];
526 : : size_t B = tet[1];
527 : : size_t C = tet[2];
528 : : size_t D = tet[3];
529 : :
530 : : trace_out << "A " << A << " B " << B << " C " << C << " D " << D
531 : : << std::endl;
532 : :
533 : : // Generate pairs of nodes (i.e edges)
534 : : // Hard coding for now, can swap out for loop
535 : : //coordinate_t AB_mid = node_connectivity->find_mid_point(A,B);
536 : 98630 : size_t AB = node_connectivity.add(A,B);
537 : :
538 : : //coordinate_t AC_mid = node_connectivity->find_mid_point(A,C);
539 : 98630 : size_t AC = node_connectivity.add(A,C);
540 : :
541 : : //coordinate_t AD_mid = node_connectivity->find_mid_point(A,D);
542 : 98630 : size_t AD = node_connectivity.add(A,D);
543 : :
544 : : //coordinate_t BC_mid = node_connectivity->find_mid_point(B,C);
545 : 98630 : size_t BC = node_connectivity.add(B,C);
546 : :
547 : : //coordinate_t BD_mid = node_connectivity->find_mid_point(B,D);
548 : 98630 : size_t BD = node_connectivity.add(B,D);
549 : :
550 : : //coordinate_t CD_mid = node_connectivity->find_mid_point(C,D);
551 : 98630 : size_t CD = node_connectivity.add(C,D);
552 : :
553 : : // Update edges
554 : :
555 : 98630 : tet_store.edge_store.split(A, C, AC, Edge_Lock_Case::unlocked);
556 : : tet_store.edge_store.split(A, B, AB, Edge_Lock_Case::unlocked);
557 : : tet_store.edge_store.split(A, D, AD, Edge_Lock_Case::unlocked);
558 : : tet_store.edge_store.split(B, C, BC, Edge_Lock_Case::unlocked);
559 : : tet_store.edge_store.split(B, D, BD, Edge_Lock_Case::unlocked);
560 : : tet_store.edge_store.split(C, D, CD, Edge_Lock_Case::unlocked);
561 : :
562 : :
563 : : // Outside edges for face ABC
564 : 98630 : tet_store.edge_store.generate(AC, BC, Edge_Lock_Case::unlocked);
565 : 98630 : tet_store.edge_store.generate(AC, AB, Edge_Lock_Case::unlocked);
566 : 98630 : tet_store.edge_store.generate(AB, BC, Edge_Lock_Case::unlocked);
567 : :
568 : : // Outside edges for face ACD
569 : 98630 : tet_store.edge_store.generate(AC, AD, Edge_Lock_Case::unlocked);
570 : 98630 : tet_store.edge_store.generate(AD, CD, Edge_Lock_Case::unlocked);
571 : 98630 : tet_store.edge_store.generate(AC, CD, Edge_Lock_Case::unlocked);
572 : :
573 : : // Outside edges for face BCD
574 : 98630 : tet_store.edge_store.generate(BD, CD, Edge_Lock_Case::unlocked);
575 : 98630 : tet_store.edge_store.generate(BD, BC, Edge_Lock_Case::unlocked);
576 : 98630 : tet_store.edge_store.generate(CD, BC, Edge_Lock_Case::unlocked);
577 : :
578 : : // Outside edges for face ABD
579 : 98630 : tet_store.edge_store.generate(AD, BD, Edge_Lock_Case::unlocked);
580 : 98630 : tet_store.edge_store.generate(AB, AD, Edge_Lock_Case::unlocked);
581 : 98630 : tet_store.edge_store.generate(AB, BD, Edge_Lock_Case::unlocked);
582 : :
583 : : // Interior Edges
584 : 98630 : tet_store.edge_store.generate(AC, BD, Edge_Lock_Case::unlocked);
585 : 98630 : tet_store.edge_store.generate(CD, AD, Edge_Lock_Case::unlocked);
586 : :
587 : : // Add the new tets
588 : : //
589 : : // External
590 : : // A AB AC AD - A
591 : : // B BA BC BD - B
592 : : // C CA CB CD - C
593 : : // D DA DB DC - D
594 : : // -
595 : : // Internal (for a face BDC, it's the intermediate and mid opposite)
596 : : // BC CD DB AC - BDC
597 : : // AB BD AD AC - ABD
598 : : // AB AC BC BD - ABC
599 : : // AC AD CD BD - ACD
600 : : //
601 : :
602 : : // TODO: This is actually generating IDs not trying to get them
603 : : child_id_list_t child = generate_child_ids(tet_store,tet_id);
604 : :
605 : : // This order should give a positive Jacobian
606 [ + - ]: 98630 : tet_store.add(child[0], A, AB, AC, AD, Refinement_Case::one_to_eight, tet_id);
607 [ + - ]: 98630 : tet_store.add(child[1], B, BC, AB, BD, Refinement_Case::one_to_eight, tet_id);
608 [ + - ]: 98630 : tet_store.add(child[2], C, AC, BC, CD, Refinement_Case::one_to_eight, tet_id);
609 [ + - ]: 98630 : tet_store.add(child[3], D, AD, CD, BD, Refinement_Case::one_to_eight, tet_id);
610 : :
611 [ + - ]: 98630 : tet_store.add(child[4], BC, CD, AC, BD, Refinement_Case::one_to_eight, tet_id);
612 [ + - ]: 98630 : tet_store.add(child[5], AB, BD, AC, AD, Refinement_Case::one_to_eight, tet_id);
613 [ + - ]: 98630 : tet_store.add(child[6], AB, BC, AC, BD, Refinement_Case::one_to_eight, tet_id);
614 [ + - ]: 98630 : tet_store.add(child[7], AC, BD, CD, AD, Refinement_Case::one_to_eight, tet_id);
615 : :
616 : : tet_store.deactivate(tet_id);
617 : :
618 : : //trace_out << "1:8 DOING REFINE OF " << tet_id << ". "
619 : : // << child[0] << ", "
620 : : // << child[1] << ", "
621 : : // << child[2] << ", "
622 : : // << child[3] << ", "
623 : : // << child[4] << ", "
624 : : // << child[5] << ", "
625 : : // << child[6] << ", "
626 : : // << child[7]
627 : : // << std::endl;
628 : :
629 : : }
630 : :
631 : : // This is just a simple assignment, but I wanted to abstract it
632 : : // for if we change the underlying type to something which a simple
633 : : // assignment is no longer safe
634 : : /**
635 : : * @brief Function to duplicate (deep copy) a tet. Useful for when
636 : : * you want to make a tet that's very similar to an existing one
637 : : *
638 : : * @param out The tet to store the copy
639 : : * @param original The tet to copy the data from
640 : : */
641 : : void copy_tet(tet_t* out, tet_t* original)
642 : : {
643 : : // NOTE: This will do a deep copy, so is safer than it may look
644 : 1022 : *out = *original;
645 : : }
646 : :
647 : : // This is just a std::replace, but may need to be more complicated
648 : : // in the future?
649 : : /**
650 : : * @brief function to take an existing list of tet ids and
651 : : * replace one. This can be useful for when you want to build very
652 : : * similar tets which share nodes
653 : : *
654 : : * @param tet Tet to perform operation on
655 : : * @param remove Element to be replaced
656 : : * @param add Element to replace with
657 : : */
658 : : void replace_node(tet_t* tet, size_t remove, size_t add)
659 : : {
660 : : std::replace(tet->begin(), tet->end(), remove, add);
661 : : }
662 : :
663 : : /**
664 : : * @brief Function to find out slot in the x,y,z data arrays a tet lives
665 : : *
666 : : * // NOTE: this is _currently_ trivial, but will be nice if we for
667 : : * example swap data stores to a map
668 : : *
669 : : * @param tet_store Tet store to use
670 : : * @param tet tet of the tet to look for
671 : : * @param element offset into that tet to look at
672 : : *
673 : : * @return tet into data arrays the tet lives
674 : : */
675 : : // TODO: Move this (or rename?)
676 : : size_t tet_id_to_node_id( tet_store_t& tet_store, size_t tet, size_t element) {
677 : : return tet_store.get(tet)[element];
678 : : }
679 : :
680 : : /**
681 : : * @brief Function to find the nodes which make up the
682 : : * single (or first?º edge which needs to be refined in an given
683 : : * edge_list
684 : : *
685 : : * @param tet_store Tet store to use
686 : : * @param edge_list The edge list to search for a refinement edge
687 : : *
688 : : * @return The node pair which represent the edge which needs
689 : : * refining
690 : : */
691 : : node_pair_t find_single_refinement_nodes( tet_store_t& tet_store, edge_list_t edge_list)
692 : : {
693 : : node_pair_t returned_nodes;
694 : : bool found_break = false;
695 : : for (size_t k = 0; k < NUM_TET_EDGES; k++)
696 : : {
697 : : edge_t edge = edge_list[k];
698 : :
699 : : if (tet_store.edge_store.get(edge).needs_refining == 1)
700 : : {
701 : : returned_nodes[0] = edge.first();
702 : : returned_nodes[1] = edge.second();
703 : :
704 : : trace_out << "1:2 needs to be split on " <<
705 : : returned_nodes[0] << " and " <<
706 : : returned_nodes[1] << std::endl;
707 : :
708 : : found_break = true;
709 : : break;
710 : : }
711 : : }
712 : :
713 : : assert(found_break);
714 : :
715 : : return returned_nodes;
716 : : }
717 : :
718 : 684 : void lock_intermediates(
719 : : tet_store_t& tet_store,
720 : : std::unordered_set<size_t> intermediate_list,
721 : : Edge_Lock_Case lock_case
722 : : )
723 : : {
724 : : // Loop over all edges
725 : : // If the edge is in the intermediate_list, deal with it
726 [ + + ]: 2214247 : for (const auto& p : tet_store.edge_store.edges)
727 : : {
728 : 2213563 : auto e = p.first;
729 : 2213563 : size_t k1 = e.first();
730 : 2213563 : size_t k2 = e.second();
731 : : // Can we make this double search cheaper?
732 : : if (
733 [ + + ][ + + ]: 4420383 : (intermediate_list.count(k1)) ||
734 : : (intermediate_list.count(k2))
735 : : )
736 : : {
737 : : trace_out << "Locking intermediate " << e << " from " << k1 << " and " << k2 << std::endl;
738 [ + - ]: 15541 : tet_store.edge_store.get(e).lock_case = lock_case;
739 [ + - ]: 15541 : tet_store.edge_store.get(e).needs_refining = 0;
740 : : }
741 : : }
742 : :
743 : 684 : }
744 : :
745 : : // TODO: remove this, it's horrible and not efficient.
746 : : // WARNING: THIS GOES OVER ALL TETS!!!!
747 : : void lock_edges_from_node(
748 : : tet_store_t& tet_store,
749 : : size_t node_id,
750 : : Edge_Lock_Case lock_case
751 : : )
752 : : {
753 : : // Iterate over edges of ALL tet
754 : : for (const auto& kv : tet_store.tets)
755 : : {
756 : : size_t tet_id = kv.first;
757 : : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
758 : : for (size_t k = 0; k < NUM_TET_EDGES; k++)
759 : : {
760 : : // If it contains that node id, mark it using lock_case
761 : : edge_t edge = edge_list[k];
762 : :
763 : : size_t edge_node_A_id = edge.first();
764 : : size_t edge_node_B_id = edge.second();
765 : :
766 : : if ((edge_node_A_id == node_id) || (edge_node_B_id == node_id)) {
767 : : trace_out << " found node in " << edge_node_A_id << " - " << edge_node_B_id << " set to " << lock_case << std::endl;
768 : : tet_store.edge_store.get(edge).lock_case = lock_case;
769 : : tet_store.edge_store.get(edge).needs_refining = 0;
770 : : }
771 : : }
772 : : }
773 : : }
774 : : // void lock_edges_from_node(
775 : : // tet_store_t& tet_store,
776 : : // size_t tet_id,
777 : : // size_t node_id,
778 : : // Edge_Lock_Case lock_case
779 : : // )
780 : : // {
781 : : // // Iterate over edges of of tet
782 : : // edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
783 : : // for (size_t k = 0; k < NUM_TET_EDGES; k++)
784 : : // {
785 : : // // If it contains that node id, mark it using lock_case
786 : : // edge_t edge = edge_list[k];
787 : : //
788 : : // size_t edge_node_A_id = edge.first();
789 : : // size_t edge_node_B_id = edge.second();
790 : : //
791 : : // if ((edge_node_A_id == node_id) || (edge_node_B_id == node_id)) {
792 : : // tet_store.edge_store.get(edge).lock_case = lock_case;
793 : : // }
794 : : // }
795 : : // }
796 : :
797 : :
798 : : ///// DEREFINEMENT STARTS HERE /////
799 : : /**
800 : : * @brief Function to iterate over children and remove them
801 : : *
802 : : * @param tet_store Tet store to use
803 : : * @param parent_id Id of the parent for whom you will delete the
804 : : * children
805 : : */
806 : : void derefine_children(tet_store_t& tet_store, size_t parent_id)
807 : : {
808 : : // For a given tet_id, find and delete its children
809 : : Refinement_State& parent = tet_store.data(parent_id);
810 : : for (auto c : parent.children)
811 : : {
812 : : tet_store.erase(c);
813 : : //tet_store.deactivate(c);
814 : :
815 : : /*
816 : : auto children = tet_store.data(c).children;
817 : : // Debug printing
818 : : std::cout << "tet " << c << "has ";
819 : : for (auto child : children)
820 : : {
821 : : std::cout << " _child " << child;
822 : : }
823 : : std::cout << std::endl;
824 : : */
825 : : }
826 : : parent.children.clear();
827 : : }
828 : :
829 : : /**
830 : : * @brief Common code for derefinement. Deactives the children and
831 : : * actives the parent
832 : : *
833 : : * @param tet_store Tet store to use
834 : : * @param parent_id The id of the parent
835 : : */
836 : : void generic_derefine(tet_store_t& tet_store, size_t parent_id)
837 : : {
838 : 40417 : derefine_children(tet_store,parent_id);
839 : 40417 : tet_store.activate(parent_id);
840 : : }
841 : :
842 : : /**
843 : : * @brief Perform 2->1 derefinement on tet
844 : : *
845 : : * @param tet_store Tet store to use
846 : : * @param parent_id The id of the parent
847 : : */
848 : : void derefine_two_to_one(tet_store_t& tet_store, node_connectivity_t&, size_t parent_id)
849 : : {
850 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
851 : : // build a delete-list of edges/intermediates first, mesh_adapter
852 : : // deletes edges from this list later
853 : : determine_deletelist_of_intermediates(tet_store, parent_id);
854 : : generic_derefine(tet_store,parent_id);
855 : : }
856 : :
857 : : /**
858 : : * @brief Perform 4->1 derefinement on tet
859 : : *
860 : : * @param tet_store Tet store to use
861 : : * @param parent_id The id of the parent
862 : : */
863 : : void derefine_four_to_one(tet_store_t& tet_store, node_connectivity_t&, size_t parent_id)
864 : : {
865 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
866 : : // build a delete-list of edges/intermediates first, mesh_adapter
867 : : // deletes edges from this list later
868 : : determine_deletelist_of_intermediates(tet_store, parent_id);
869 : : generic_derefine(tet_store,parent_id);
870 : : }
871 : :
872 : : /**
873 : : * @brief Perform 8->1 derefinement on tet
874 : : *
875 : : * @param tet_store Tet store to use
876 : : * @param parent_id The id of the parent
877 : : */
878 : 40417 : void derefine_eight_to_one(tet_store_t& tet_store, node_connectivity_t&, size_t parent_id)
879 : : {
880 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
881 : :
882 : : generic_derefine(tet_store,parent_id);
883 : :
884 : : // TODO: Do we delete the nodes? Do we even have nodes?
885 : :
886 : : // Delete the center edges
887 : : // If edge isn't in the parent, delete it? Is there a better way?
888 : 40417 : edge_list_t parent_edges = tet_store.generate_edge_keys(parent_id);
889 : :
890 : : Refinement_State& parent = tet_store.data(parent_id);
891 [ - + ]: 40417 : for (auto c : parent.children)
892 : : {
893 [ - - ]: 0 : edge_list_t child_edges = tet_store.generate_edge_keys(c);
894 : : // build a delete-list of non-matching edges first, then
895 : : // mesh_adapter deletes edges from this list later
896 [ - - ]: 0 : determine_deletelist_of_non_matching_edges(child_edges, parent_edges);
897 : : }
898 : 40417 : }
899 : :
900 : : // TODO: Document This.
901 : 0 : void derefine_four_to_two(tet_store_t& tet_store, node_connectivity_t& node_connectivity, size_t parent_id)
902 : : {
903 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
904 : : auto edge = find_edge_not_derefined(tet_store,
905 : 0 : node_connectivity, parent_id);
906 : 0 : derefine_four_to_one(tet_store, node_connectivity, parent_id);
907 : 0 : refine_one_to_two( tet_store, node_connectivity, parent_id,
908 : : edge.first(), edge.second() );
909 : 0 : }
910 : :
911 : : // TODO: Document This.
912 : 0 : void derefine_eight_to_two(tet_store_t& tet_store, node_connectivity_t& node_connectivity, size_t parent_id)
913 : : {
914 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
915 : :
916 : : auto edge = find_edge_not_derefined(tet_store,
917 : 0 : node_connectivity, parent_id);
918 : 0 : derefine_eight_to_one(tet_store, node_connectivity, parent_id);
919 : 0 : refine_one_to_two( tet_store, node_connectivity, parent_id,
920 : : edge.first(), edge.second() );
921 : 0 : }
922 : :
923 : : // TODO: Document This.
924 : 0 : void derefine_eight_to_four(tet_store_t& tet_store, node_connectivity_t& node_connectivity, size_t parent_id)
925 : : {
926 : : //if (!check_allowed_derefinement(tet_store,parent_id)) return;
927 : : // TODO: think about if the logic for these derefs are right
928 : 0 : derefine_eight_to_one(tet_store, node_connectivity, parent_id);
929 : 0 : deref_refine_one_to_four( tet_store, node_connectivity, parent_id);
930 : 0 : }
931 : :
932 : : /**
933 : : * @brief Loop over children and determine delete-list of all intermediate edges
934 : : *
935 : : * @param tet_store Tet store to use
936 : : * @param parent_id Id of parent
937 : : */
938 : 343 : void determine_deletelist_of_intermediates(tet_store_t& tet_store, size_t parent_id)
939 : : {
940 : : Refinement_State& parent = tet_store.data(parent_id);
941 : 343 : auto parent_edges = tet_store.generate_edge_keys(parent_id);
942 : : std::set< edge_t > parent_edge_set;
943 [ + + ][ + - ]: 2401 : for (auto pe:parent_edges) parent_edge_set.insert(pe);
944 [ + + ]: 1425 : for (auto c : parent.children)
945 : : {
946 [ + - ]: 1082 : edge_list_t edge_list = tet_store.generate_edge_keys(c);
947 [ + + ]: 7574 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
948 : : {
949 : 6492 : edge_t edge = edge_list[k];
950 : : // accept this code may try delete an edge which has already gone
951 : 6492 : if (tet_store.edge_store.exists(edge)) {
952 : : if (parent_edge_set.count(edge) == 0)
953 : : {
954 : : trace_out << "child " << c << " adding to delete list: "
955 : : << edge.first() << " - " << edge.second() << std::endl;
956 [ + - ]: 4936 : delete_list.push_back(edge);
957 : : }
958 : : }
959 : : }
960 : : }
961 : 343 : }
962 : :
963 : : /**
964 : : * @brief Deletes the intermediate edge in the delete list for derefinement
965 : : *
966 : : * @param tet_store Tet store to use
967 : : */
968 : 326 : void delete_intermediates_of_children(tet_store_t& tet_store)
969 : : {
970 [ + + ]: 5262 : for (const auto& edge : delete_list) {
971 : 4936 : tet_store.edge_store.erase(edge);
972 : : }
973 : :
974 : : delete_list.clear();
975 : 326 : }
976 : :
977 : : /**
978 : : * @brief If edge in candidate is not present in basis, add edge
979 : : * (candidate) to delete list
980 : : *
981 : : * @param candidate The edge list which is to be searched and deleted
982 : : * @param basis The edge list to check against
983 : : */
984 : 0 : void determine_deletelist_of_non_matching_edges(edge_list_t candidate,
985 : : edge_list_t basis)
986 : : {
987 : : trace_out << "Looking for edges to delete" << std::endl;
988 : :
989 : : // TODO: Sanity check this now we changed to edge_t
990 : :
991 : : // Loop over the edges in each child. Look over the basis and
992 : : // if we can't find it, delete it
993 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
994 : : {
995 : 0 : edge_t search_key = candidate[k];
996 : :
997 : : // Search the basis for it
998 : : bool found_it = false;
999 : :
1000 [ - - ]: 0 : for (size_t l = 0; l < NUM_TET_EDGES; l++)
1001 : : {
1002 [ - - ]: 0 : edge_t key = basis[l];
1003 : : if (search_key == key)
1004 : : {
1005 : : found_it = true;
1006 : : }
1007 : : }
1008 : :
1009 : : // If we didn't find it, delete it
1010 [ - - ]: 0 : if (!found_it)
1011 : : {
1012 : : // Delete it
1013 : : //tet_store.edge_store.erase(search_key);
1014 : 0 : delete_list.push_back(search_key);
1015 : : }
1016 : : }
1017 : 0 : }
1018 : :
1019 : :
1020 : : /**
1021 : : * @brief function to detect when an invalid derefinement is
1022 : : * invoked
1023 : : *
1024 : : * @param tet_store Tet store to use
1025 : : * @param tet_id Id the of the tet which will be de-refined
1026 : : *
1027 : : * @return A bool stating if the tet can be validly de-refined
1028 : : */
1029 : : bool check_allowed_derefinement( tet_store_t& tet_store, size_t tet_id)
1030 : : {
1031 : : Refinement_State& master_element = tet_store.data(tet_id);
1032 : :
1033 : : // Check this won't take us past the max refinement level
1034 : : if (master_element.refinement_level <= MIN_REFINEMENT_LEVEL)
1035 : : {
1036 : : return false;
1037 : : }
1038 : :
1039 : : // If we got here, we didn't detect anything which tells us not
1040 : : // to
1041 : : return true;
1042 : : }
1043 : :
1044 : :
1045 : : // HERE BE DRAGONS! THIS IS DANGEROUS IF YOU USE IT WRONG
1046 : : // For every child of parent_id, set his children to our won
1047 : : // TODO: set a flag for the curious user to know we trashed the children
1048 : : void overwrite_children(
1049 : : tet_store_t& tet_store,
1050 : : const child_id_list_t& to_be_replaced,
1051 : : const child_id_list_t& replace_with
1052 : : )
1053 : : {
1054 : : for (auto c : to_be_replaced)
1055 : : {
1056 : : tet_store.data(c).children = replace_with;
1057 : : }
1058 : : }
1059 : :
1060 : : /**
1061 : : * @brief function to detect which edge should not get derefined
1062 : : *
1063 : : * @param tet_store Tet store to use
1064 : : * @param node_connectivity Node connectivity to use
1065 : : * @param tet_id Id the of the tet which will be de-refined
1066 : : *
1067 : : * @return Array of size two containing nodes of required edge
1068 : : */
1069 : 0 : edge_t find_edge_not_derefined(
1070 : : tet_store_t& tet_store,
1071 : : node_connectivity_t& node_connectivity,
1072 : : size_t tet_id)
1073 : : {
1074 : : // 2 nonparent nodes set to derefine for a 4:2
1075 : : // 5 nonparent nodes set to derefine for an 8:2
1076 : : // will have 2 or 5 nonparent nodes set to deref; Figure out which
1077 : : // edge is the one that is not set to deref
1078 : :
1079 : 0 : auto derefine_node_set = find_derefine_node_set(tet_store, tet_id);
1080 : :
1081 : : //// Do number of points
1082 : : //std::unordered_set<size_t> derefine_node_set;
1083 : :
1084 : : // Find the set of nodes which are not in the parent
1085 : : std::unordered_set<size_t> non_parent_nodes =
1086 [ - - ]: 0 : child_exclusive_nodes(tet_store, tet_id);
1087 : :
1088 : : // from the above non_parent_nodes set and derefine_node_set,
1089 : : // figureout which node should be removed
1090 : : std::size_t ed_A(0), ed_B(0);
1091 [ - - ]: 0 : for (auto npn:non_parent_nodes) {
1092 [ - - ]: 0 : if (derefine_node_set.count(npn) == 0) {
1093 : : // we've found the node that should not be removed, now we
1094 : : // need to find the edge it belongs to
1095 : 0 : auto nonderef_edge = node_connectivity.get(npn);
1096 : 0 : ed_A = nonderef_edge[0];
1097 : 0 : ed_B = nonderef_edge[1];
1098 : : //std::cout << "do-not-deref-APAN " << "A " << nd_edge[0]
1099 : : // << " B " << nd_edge[1] << std::endl;
1100 : : }
1101 : : }
1102 : :
1103 : : assert(ed_A!=ed_B);
1104 : 0 : edge_t nd_edge(ed_A, ed_B);
1105 : 0 : return nd_edge;
1106 : : }
1107 : :
1108 : : /**
1109 : : * @brief function to detect what intermediate/non-parent nodes are
1110 : : * marked for derefinement
1111 : : *
1112 : : * @param tet_store Tet store to use
1113 : : * @param tet_id Id of the tet which will be de-refined
1114 : : *
1115 : : * @return Set of nodes of marked for derefinement
1116 : : */
1117 : 83268 : std::unordered_set< size_t > find_derefine_node_set(
1118 : : tet_store_t& tet_store,
1119 : : size_t tet_id)
1120 : : {
1121 : : // Set of nodes which are not in the parent
1122 : : std::unordered_set<size_t> non_parent_nodes =
1123 : 83268 : child_exclusive_nodes(tet_store, tet_id);
1124 : : std::unordered_set<size_t> derefine_node_set, unmarked_deref_node_set,
1125 : : final_deref_node_set;
1126 : :
1127 [ + - ]: 83268 : child_id_list_t children = tet_store.data(tet_id).children;
1128 : :
1129 : : // Look at children
1130 : : trace_out << tet_id << " Looping over " << children.size() << "children" << std::endl;
1131 [ + + ]: 744300 : for (size_t i = 0; i < children.size(); i++)
1132 : : {
1133 : : trace_out << "child: " << children[i] << std::endl;
1134 : : // TODO: Is this in element or tet ids?
1135 [ + - ]: 661032 : edge_list_t edge_list = tet_store.generate_edge_keys(children[i]);
1136 [ + + ]: 4627224 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1137 : : {
1138 : 3966192 : edge_t edge = edge_list[k];
1139 : : // TODO: where do we makr the edges that need to be derefed? parent of child?
1140 : : // Check each node, see if its an intermediate
1141 : 3966192 : size_t A = edge.first();
1142 : 3966192 : size_t B = edge.second();
1143 : : trace_out << "checking edge for deref " << A << " - " << B << std::endl;
1144 : :
1145 : : //if (tet_store.is_intermediate(A))
1146 [ + + ]: 3966192 : if (non_parent_nodes.count(A) )
1147 : : {
1148 [ + - ][ + + ]: 2963412 : if (tet_store.edge_store.get(edge).needs_derefining) {
1149 : : trace_out << "Adding " << A << std::endl;
1150 : : derefine_node_set.insert(A);
1151 : : }
1152 : : else {
1153 : : unmarked_deref_node_set.insert(A);
1154 : : //trace_out << "NOT added " << A << std::endl;
1155 : : }
1156 : : }
1157 : :
1158 : : //if (tet_store.is_intermediate(B))
1159 [ + + ]: 3966192 : if (non_parent_nodes.count(B))
1160 : : {
1161 [ + - ][ + + ]: 3961674 : if (tet_store.edge_store.get(edge).needs_derefining) {
1162 : : trace_out << "Adding " << B << std::endl;
1163 : : derefine_node_set.insert(B);
1164 : : }
1165 : : else {
1166 : : unmarked_deref_node_set.insert(B);
1167 : : //trace_out << "NOT added " << B << std::endl;
1168 : : }
1169 : : }
1170 : : }
1171 : : }
1172 : :
1173 : : //trace_out << "marked for deref: " << derefine_node_set.size() << std::endl;
1174 : : //trace_out << "NOT marked for deref: " << unmarked_deref_node_set.size() << std::endl;
1175 : :
1176 : : // remove nodes that are unmarked for derefinement
1177 [ + + ]: 577852 : for (auto drnode : derefine_node_set) {
1178 [ + + ]: 494584 : if (unmarked_deref_node_set.count(drnode) == 0)
1179 : : final_deref_node_set.insert(drnode);
1180 : : }
1181 : : derefine_node_set = final_deref_node_set;
1182 : 83268 : return derefine_node_set;
1183 : : }
1184 : :
1185 : :
1186 [ + - ]: 166536 : std::unordered_set<size_t> child_exclusive_nodes(tet_store_t& tet_store,
1187 : : size_t tet_id)
1188 : : {
1189 : : std::unordered_set<size_t> non_parent_nodes;
1190 : :
1191 : : // array
1192 [ + - ]: 166536 : auto parent_tet = tet_store.get(tet_id);
1193 : :
1194 : : // convert to set
1195 [ + - ]: 166536 : std::unordered_set<size_t> parent_set(begin(parent_tet), end(parent_tet));
1196 : :
1197 [ + - ]: 166536 : child_id_list_t children = tet_store.data(tet_id).children;
1198 [ + + ]: 1488600 : for (size_t i = 0; i < children.size(); i++)
1199 : : {
1200 [ + - ]: 1322064 : auto child_tet = tet_store.get( children[i] );
1201 : :
1202 : : // Look at nodes, if not present add to set
1203 [ + + ]: 6610320 : for (std::size_t j = 0; j < NUM_TET_NODES; j++)
1204 : : {
1205 : 5288256 : auto node = child_tet[j];
1206 [ + + ]: 5288256 : if (parent_set.count(node) == 0)
1207 : : {
1208 : : non_parent_nodes.insert(node);
1209 : : }
1210 : : }
1211 : : }
1212 : :
1213 : : trace_out <<" Found " << non_parent_nodes.size() << " non parent nodes " << std::endl;
1214 : 166536 : return non_parent_nodes;
1215 : :
1216 : : }
1217 : : };
1218 : : }
1219 : :
1220 : : #if defined(STRICT_GNUC)
1221 : : #pragma GCC diagnostic pop
1222 : : #endif
1223 : :
1224 : : #endif // guard
|
