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