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