Branch data Line data Source code
1 : : #include "mesh_adapter.hpp"
2 : :
3 : : #include <assert.h> // for assert
4 : : #include <cstddef> // for size_t
5 : : #include <iostream> // for operator<<, endl, basic_os...
6 : : #include <set> // for set
7 : : #include <utility> // for pair
8 : : #include "AMR/AMR_types.hpp" // for Edge_Refinement, edge_list_t
9 : : #include "AMR/Loggers.hpp" // for trace_out
10 : : #include "AMR/Refinement_State.hpp" // for Refinement_Case, Refinemen...
11 : : #include "AMR/edge.hpp" // for operator<<, edge_t
12 : : #include "AMR/edge_store.hpp" // for edge_store_t
13 : : #include "AMR/marked_refinements_store.hpp" // for marked_refinements_store_t
14 : : #include "AMR/node_connectivity.hpp" // for node_connectivity_t
15 : : #include "AMR/refinement.hpp" // for refinement_t
16 : : #include "AMR/tet_store.hpp" // for tet_store_t
17 : :
18 : : #if defined(__clang__)
19 : : #pragma clang diagnostic push
20 : : #pragma clang diagnostic ignored "-Wunreachable-code"
21 : : #pragma clang diagnostic ignored "-Wdocumentation"
22 : : #endif
23 : :
24 : : namespace AMR {
25 : :
26 : :
27 : : #ifdef ENABLE_NODE_STORE
28 : : /**
29 : : * @brief This accepts external coord arrays and allows the node_store to
30 : : * track the new node positions as they are added
31 : : *
32 : : * @param m_x X coodinates
33 : : * @param m_y Y coodinates
34 : : * @param m_z Z coodinates
35 : : * @param graph_size Total number of nodes
36 : : */
37 : : // TODO: remove graph size and use m.size()
38 : : // TODO: remove these pointers
39 : : //void mesh_adapter_t::init_node_store(coord_type* m_x, coord_type* m_y, coord_type* m_z)
40 : : //{
41 : : // assert( m_x->size() == m_y->size() );
42 : : // assert( m_x->size() == m_z->size() );
43 : :
44 : : // node_store.set_x(*m_x);
45 : : // node_store.set_y(*m_y);
46 : : // node_store.set_z(*m_z);
47 : : //}
48 : : #endif
49 : :
50 : 22 : std::pair< bool, std::size_t > mesh_adapter_t::check_same_face(
51 : : std::size_t tet_id,
52 : : const std::unordered_set<std::size_t>& inactive_nodes)
53 : : {
54 : 22 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
55 : :
56 : : Assert(inactive_nodes.size()==3 || inactive_nodes.size()==2,
57 : : "Incorrectly sized inactive nodes set");
58 : :
59 : : // for a tet ABCD, the keys (edges) are ordered
60 : : // A-B, A-C, A-D, B-C, B-D, C-D
61 : : // 0-1, 0-2, 0-3, 1-2, 1-3, 2-3
62 : :
63 : : std::array< std::array< std::size_t, 3 >, 4 >
64 : : edges_on_face;
65 : :
66 : : // A-B-C
67 : 22 : edges_on_face[0][0] =
68 : 22 : tk::cref_find(node_connectivity.data(),edge_list[0].get_data());
69 : 22 : edges_on_face[0][1] =
70 : 22 : tk::cref_find(node_connectivity.data(),edge_list[1].get_data());
71 : 22 : edges_on_face[0][2] =
72 : 22 : tk::cref_find(node_connectivity.data(),edge_list[3].get_data());
73 : :
74 : : // A-B-D
75 : 22 : edges_on_face[1][0] =
76 : 22 : tk::cref_find(node_connectivity.data(),edge_list[0].get_data());
77 : 22 : edges_on_face[1][1] =
78 : 22 : tk::cref_find(node_connectivity.data(),edge_list[2].get_data());
79 : 22 : edges_on_face[1][2] =
80 : 22 : tk::cref_find(node_connectivity.data(),edge_list[4].get_data());
81 : :
82 : : // B-C-D
83 : 22 : edges_on_face[2][0] =
84 : 22 : tk::cref_find(node_connectivity.data(),edge_list[3].get_data());
85 : 22 : edges_on_face[2][1] =
86 : 22 : tk::cref_find(node_connectivity.data(),edge_list[4].get_data());
87 : 22 : edges_on_face[2][2] =
88 : 22 : tk::cref_find(node_connectivity.data(),edge_list[5].get_data());
89 : :
90 : : // A-C-D
91 : 22 : edges_on_face[3][0] =
92 : 22 : tk::cref_find(node_connectivity.data(),edge_list[1].get_data());
93 : 22 : edges_on_face[3][1] =
94 : 22 : tk::cref_find(node_connectivity.data(),edge_list[2].get_data());
95 : 22 : edges_on_face[3][2] =
96 : 22 : tk::cref_find(node_connectivity.data(),edge_list[5].get_data());
97 : :
98 : : //Iterate over edges to determine if inactive_nodes are all part of a face
99 : : bool same_face(false);
100 : : [[maybe_unused]] bool tnode_set(false);
101 : : std::size_t third_node = 0;
102 [ + + ]: 110 : for(const auto& face : edges_on_face)
103 : : {
104 : : std::size_t icount = 0;
105 [ + + ]: 352 : for (const auto& np_node : face) {
106 : 128 : if (inactive_nodes.count(np_node)) ++icount;
107 : : }
108 [ + + ]: 88 : if (inactive_nodes.size() == icount) {
109 : : same_face = true;
110 : : // if the two inactive_nodes being checked are on the same parent
111 : : // face, determine the third node on that face
112 [ + + ]: 16 : if (inactive_nodes.size() == 2) {
113 [ - + ][ + - ]: 2 : for (auto fn:face) {
114 : : if (inactive_nodes.count(fn) == 0) {
115 : : third_node = fn;
116 : : tnode_set = true;
117 : 2 : break;
118 : : }
119 : : }
120 : : }
121 : : }
122 : : }
123 : :
124 : : if (same_face && inactive_nodes.size() == 2)
125 : : Assert(tnode_set, "Third node on face not set in derefine");
126 : :
127 : 22 : return {same_face, third_node};
128 : : }
129 : :
130 : : /** @brief Consume an existing mesh, and turn it into the AMRs
131 : : * representations of tets and nodes
132 : : *
133 : : * @param tetinpoel Vector of nodes grouped together in blocks of 4 to
134 : : * represent tets
135 : : */
136 : 2484 : void mesh_adapter_t::consume_tets(const std::vector< std::size_t >& tetinpoel )
137 : : {
138 [ + + ]: 1055790 : for (size_t i = 0; i < tetinpoel.size(); i+=4)
139 : : {
140 : : tet_t t = {
141 : : {
142 : : tetinpoel[i],
143 : 1053306 : tetinpoel[i+1],
144 : 1053306 : tetinpoel[i+2],
145 : 1053306 : tetinpoel[i+3]
146 : : }
147 : 1053306 : };
148 : :
149 : : trace_out << "Consume tet " << i << std::endl;
150 : 1053306 : tet_store.add(t, AMR::Refinement_Case::initial_grid);
151 : : }
152 : 2484 : }
153 : :
154 : : /**
155 : : * @brief Place holder function to evaluate error estimate at
156 : : * nodes, and therefor mark things as needing to be refined
157 : : */
158 : : //void mesh_adapter_t::evaluate_error_estimate() {
159 : : // for (auto& kv : tet_store.edge_store.edges)
160 : : // {
161 : : // // Mark them as needing refinement
162 : : // if (kv.second.refinement_criteria > refinement_cut_off)
163 : : // {
164 : : // kv.second.needs_refining = 1;
165 : : // }
166 : : // else
167 : : // {
168 : : // // TODO: Check this won't be overwriting valuable
169 : : // // information from last iteration
170 : : // kv.second.needs_refining = 0;
171 : : // }
172 : : // }
173 : : //}
174 : :
175 : : /**
176 : : * @brief Helper function to apply uniform refinement to all tets
177 : : */
178 : 31 : void mesh_adapter_t::mark_uniform_refinement()
179 : : {
180 [ + + ]: 33650 : for (auto& kv : tet_store.edge_store.edges) {
181 : : auto& local = kv.second;
182 [ + - ]: 33619 : if (local.lock_case == Edge_Lock_Case::unlocked)
183 : 33619 : local.needs_refining = 1;
184 : : }
185 : 31 : mark_refinement();
186 : 31 : }
187 : :
188 : : /**
189 : : * @brief Helper function to apply uniform derefinement to all tets
190 : : */
191 : 24 : void mesh_adapter_t::mark_uniform_derefinement()
192 : : {
193 : 24 : const auto& inp = tet_store.get_active_inpoel();
194 : 24 : auto& edge_store = tet_store.edge_store;
195 [ + + ]: 49409 : for (std::size_t t=0; t<inp.size()/4; ++t) {
196 : : const auto edges =
197 : 49385 : edge_store.generate_keys(
198 : 49385 : {inp[t*4+0], inp[t*4+1], inp[t*4+2], inp[t*4+3]});
199 [ + + ]: 345695 : for (const auto& tetedge : edges) {
200 : : auto e = edge_store.edges.find(tetedge);
201 [ + - ]: 296310 : if (e != end(edge_store.edges)) {
202 : : auto& local = e->second;
203 : 296310 : local.needs_derefining = 1;
204 : : }
205 : : }
206 : : }
207 : 24 : mark_derefinement();
208 : 24 : }
209 : :
210 : : /**
211 : : * @brief For a given set of edges, set their refinement criteria for
212 : : * refinement
213 : : *
214 : : * @param remote Vector of edges and edge tags
215 : : */
216 : 16 : void mesh_adapter_t::mark_error_refinement(
217 : : const std::vector< std::pair< edge_t, edge_tag > >& remote )
218 : : {
219 [ + + ]: 5360 : for (const auto& r : remote) {
220 : 5344 : auto& local = tet_store.edge_store.get( r.first );
221 [ + + ]: 5344 : if (r.second == edge_tag::REFINE) {
222 [ - + ]: 3612 : if (local.lock_case > Edge_Lock_Case::unlocked) {
223 : 0 : local.needs_refining = 0;
224 : : } else {
225 : 3612 : local.needs_refining = 1;
226 : : // an edge deemed to be 'needing refinement', cannot be derefined
227 : 3612 : local.needs_derefining = 0;
228 : : }
229 [ + - ]: 1732 : } else if (r.second == edge_tag::DEREFINE) {
230 : 1732 : local.needs_derefining = 1;
231 : : }
232 : : }
233 : :
234 : 16 : mark_refinement();
235 : 16 : mark_derefinement();
236 : 16 : }
237 : :
238 : 0 : void mesh_adapter_t::mark_error_refinement_corr( const EdgeData& edges )
239 : : {
240 [ - - ]: 0 : for (const auto& r : edges)
241 : : {
242 : 0 : auto& edgeref = tet_store.edge_store.get( edge_t(r.first) );
243 : 0 : edgeref.needs_refining = std::get<0>(r.second);
244 : 0 : edgeref.needs_derefining = std::get<1>(r.second);
245 : : Assert(edgeref.lock_case == Edge_Lock_Case::unlocked ?
246 : : edgeref.lock_case <= std::get<2>(r.second) : true,
247 : : "Edge " + std::to_string(r.first[0]) +
248 : : "-" + std::to_string(r.first[1]) +
249 : : " : current edge-lock " + std::to_string(edgeref.lock_case) +
250 : : " stricter than received edge-lock " +
251 : : std::to_string(std::get<2>(r.second)));
252 : 0 : edgeref.lock_case = std::get<2>(r.second);
253 : : }
254 : 0 : mark_refinement();
255 : 0 : mark_derefinement();
256 : 0 : }
257 : :
258 : : /**
259 : : * @brief Function to detect the compatibility class (1,
260 : : * 2, or 3) based on the number of locked edges and the existence
261 : : * of intermediate edges
262 : : *
263 : : * @param num_locked_edges The number of locked edges
264 : : * @param num_intermediate_edges The number of intermediate edges
265 : : * @param refinement_case The refinement case of the tet
266 : : * @param normal TODO: Document this!
267 : : *
268 : : * @return The compatibili4y class of the current scenario
269 : : */
270 : 50634 : int mesh_adapter_t::detect_compatibility(
271 : : int num_locked_edges,
272 : : int num_intermediate_edges,
273 : : AMR::Refinement_Case refinement_case,
274 : : int normal
275 : : )
276 : : {
277 : : int compatibility = 0;
278 : 50634 : num_locked_edges += num_intermediate_edges;
279 : :
280 : : /*
281 : : // Split this into three categories
282 : : // 1. Normal elements without locked edges. => 1
283 : : //if (normal) {
284 : :
285 : : // 3. Intermediate elements with at least one edge marked for refinement => 3
286 : : if (num_intermediate_edges > 0)
287 : : {
288 : : compatibility = 3;
289 : : }
290 : : else if (num_locked_edges == 0) {
291 : : compatibility = 1;
292 : : }
293 : : // 2. Normal elements with locked edges. => 2
294 : : else {
295 : : compatibility = 2;
296 : : }
297 : : //}
298 : : */
299 : :
300 : : //else {
301 : : //if (num_intermediate_edges > 0) { compatibility = 3; }
302 : : //}
303 : :
304 : :
305 : :
306 : : // Only 1:2 and 1:4 are intermediates and eligible for class3 // NOT TRUE!
307 : : /*
308 : : if (num_intermediate_edges > 0)
309 : : {
310 : : if (!normal) {
311 : : trace_out << " not normal 3 " << std::endl;
312 : : compatibility = 3;
313 : : }
314 : : else { // Attempt to allow for "normal" 1:4 and 1:8
315 : : compatibility = 2;
316 : : trace_out << " normal 3 " << std::endl;
317 : : }
318 : :
319 : : }
320 : : else {
321 : : if (num_locked_edges == 0) {
322 : : trace_out << " no lock 1 " << std::endl;
323 : : compatibility = 1;
324 : : }
325 : : else {
326 : : trace_out << " lock 2 " << std::endl;
327 : : compatibility = 2;
328 : : }
329 : : }
330 : : */
331 : :
332 : :
333 : : // Old implementation
334 : : // Only 1:2 and 1:4 are intermediates and eligible for class3 // NOT TRUE!
335 : 50634 : if (
336 [ - + ]: 50634 : (refinement_case == AMR::Refinement_Case::one_to_two) or
337 : : (refinement_case == AMR::Refinement_Case::one_to_four)
338 : : )
339 : : {
340 [ - - ]: 0 : if (!normal) {
341 : : trace_out << " not normal 3 " << std::endl;
342 : : compatibility = 3;
343 : : }
344 : : else { // Attempt to allow for "normal" 1:4 and 1:8
345 : : compatibility = 2;
346 : : trace_out << " normal 3 " << std::endl;
347 : : }
348 : :
349 : : }
350 : : else {
351 [ - + ]: 50634 : if (num_locked_edges == 0) {
352 : : trace_out << " no lock 1 " << std::endl;
353 : : compatibility = 1;
354 : : }
355 : : else {
356 : : trace_out << " lock 2 " << std::endl;
357 : : compatibility = 2;
358 : : }
359 : : }
360 : :
361 : : assert(compatibility > 0);
362 : : assert(compatibility < 4);
363 : 50634 : return compatibility;
364 : : }
365 : :
366 : : /**
367 : : * @brief Function which implements the main refinement algorithm from
368 : : * the paper Iterating over the cells, deciding which refinement and
369 : : * compatibility types are appropriate etc
370 : : */
371 : 47 : void mesh_adapter_t::mark_refinement() {
372 : :
373 : : #ifndef AMR_MAX_ROUNDS
374 : : // Paper says the average actual num rounds will be 5-15
375 : : #define AMR_MAX_ROUNDS 50
376 : : #endif
377 : : const size_t max_num_rounds = AMR_MAX_ROUNDS;
378 : :
379 : : // Mark refinements
380 : : size_t iter;
381 : : //Iterate until convergence
382 [ + - ]: 107 : for (iter = 0; iter < max_num_rounds; iter++)
383 : : {
384 : :
385 : 107 : tet_store.marked_refinements.get_state_changed() = false;
386 : :
387 : : // Loop over Tets.
388 [ + + ]: 56949 : for (const auto& kv : tet_store.tets)
389 : : {
390 : 56842 : size_t tet_id = kv.first;
391 : :
392 : : trace_out << "Process tet " << tet_id << std::endl;
393 : :
394 : : // Only apply checks to tets on the active list
395 : : if (tet_store.is_active(tet_id)) {
396 : : int num_locked_edges = 0;
397 : : int num_intermediate_edges = 0;
398 : :
399 : : // Loop over nodes and count the number which need refining
400 : : int num_to_refine = 0;
401 : :
402 : : // This is useful for later inspection
403 [ + - ]: 51555 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
404 : :
405 : : //Iterate over edges
406 [ + + ]: 360885 : for(auto & key : edge_list)
407 : : {
408 : :
409 : : trace_out << "Edge " << key << std::endl;
410 : :
411 : : //Count locked edges and edges in need of
412 : : // refinement
413 : : // Count Locked Edges
414 [ + - ][ - + ]: 309330 : if(tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::locked)
415 : : {
416 : : trace_out << "Found locked edge " << key << std::endl;
417 : : trace_out << "Locked :" << tet_store.edge_store.get(key).lock_case << std::endl;
418 : 0 : num_locked_edges++;
419 : : }
420 [ + - ]: 309330 : else if(tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::intermediate
421 [ + - ][ + - ]: 309330 : || tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::temporary)
[ + - ]
422 : : {
423 : : trace_out << "Found intermediate edge " << key << std::endl;
424 : 0 : num_intermediate_edges++;
425 : : }
426 : : else
427 : : {
428 : : // Count edges which need refining
429 : : // We check in here as we won't refine a
430 : : // locked edge and will thus ignore it
431 [ + - ][ + + ]: 309330 : if (tet_store.edge_store.get(key).needs_refining == 1)
432 : : {
433 : 298129 : num_to_refine++;
434 : : trace_out << "key needs ref " << key << std::endl;
435 : : }
436 : : }
437 : : }
438 : :
439 : : // TODO: Should this be a reference?
440 : : AMR::Refinement_Case refinement_case = tet_store.get_refinement_case(tet_id);
441 : 51555 : int normal = tet_store.is_normal(tet_id);
442 : :
443 : : trace_out << "Checking " << tet_id <<
444 : : " ref case " << refinement_case <<
445 : : " num ref " << num_to_refine <<
446 : : " normal " << normal <<
447 : : std::endl;
448 : :
449 : :
450 : :
451 : : //If we have any tets to refine
452 [ + + ]: 51555 : if (num_to_refine > 0)
453 : : {
454 : : //Determine compatibility case
455 : :
456 [ + - ]: 50634 : int compatibility = detect_compatibility(num_locked_edges,
457 : : num_intermediate_edges, refinement_case, normal);
458 : :
459 : : trace_out << "Compat " << compatibility << std::endl;
460 : :
461 : : // Now check num_to_refine against situations
462 [ + - ]: 50634 : if (compatibility == 1)
463 : : {
464 [ + - ]: 50634 : refinement_class_one(num_to_refine, tet_id);
465 : : }
466 [ - - ]: 0 : else if (compatibility == 2)
467 : : {
468 [ - - ]: 0 : refinement_class_two(edge_list, tet_id);
469 : : }
470 [ - - ]: 0 : else if (compatibility == 3)
471 : : {
472 [ - - ]: 0 : refinement_class_three(tet_id);
473 : : }
474 : :
475 : : /*
476 : : // Write temp mesh out
477 : : std::string temp_file = "temp." +
478 : : std::to_string(iter) + "." +
479 : : std::to_string(tet_id) + ".exo";
480 : :
481 : : std::cout << "Writing " << temp_file << std::endl;
482 : : Adaptive_UnsMesh outmesh(
483 : : get_active_inpoel(), x(), y(), z()
484 : : );
485 : : tk::ExodusIIMeshWriter( temp_file, tk::ExoWriter::CREATE ).
486 : : writeMesh(outmesh);
487 : : */
488 : :
489 : : } // if num_to_refine
490 : : else {
491 : : // If we got here, we don't want to refine this guy
492 [ + - ]: 921 : tet_store.marked_refinements.add(tet_id, AMR::Refinement_Case::none);
493 : : }
494 : : } // if active
495 : : else {
496 : : trace_out << "Inactive" << std::endl;
497 : : }
498 : : } // For
499 : :
500 : : // If nothing changed during that round, break
501 [ + + ]: 107 : if (!tet_store.marked_refinements.get_state_changed())
502 : : {
503 : : trace_out << "Terminating loop at iter " << iter << std::endl;
504 : : break;
505 : : }
506 : : trace_out << "End iter " << iter << std::endl;
507 : : }
508 : : trace_out << "Loop took " << iter << " rounds." << std::endl;
509 : :
510 : : //std::cout << "Print Tets" << std::endl;
511 : : //print_tets();
512 : 47 : }
513 : :
514 : : /**
515 : : * @brief Helper function to print tet information when needed
516 : : */
517 : 0 : void mesh_adapter_t::print_tets() {
518 : : tet_store.print_tets();
519 : 0 : }
520 : :
521 : : /**
522 : : * @brief Function to call refinement after each tet has had it's
523 : : * refinement case marked and calculated
524 : : */
525 : 47 : void mesh_adapter_t::perform_refinement()
526 : : {
527 : : // Track tets which need to be deleted this iteration
528 : : std::set<size_t> round_two;
529 : :
530 : : trace_out << "Perform ref" << std::endl;
531 : :
532 : : // Do refinements
533 [ + + ]: 137891 : for (const auto& kv : tet_store.tets)
534 : : {
535 [ + + ]: 137844 : size_t tet_id = kv.first;
536 : :
537 : : trace_out << "Do refine of " << tet_id << std::endl;
538 : : if (tet_store.has_refinement_decision(tet_id))
539 : : {
540 [ + + ][ + - ]: 23419 : switch(tet_store.marked_refinements.get(tet_id))
[ - + ]
541 : : {
542 : 130 : case AMR::Refinement_Case::one_to_two:
543 [ + - ]: 130 : refiner.refine_one_to_two(tet_store,node_connectivity,tet_id);
544 : : break;
545 : 176 : case AMR::Refinement_Case::one_to_four:
546 [ + - ]: 176 : refiner.refine_one_to_four(tet_store,node_connectivity,tet_id);
547 : : break;
548 : 22907 : case AMR::Refinement_Case::one_to_eight:
549 [ + - ]: 22907 : refiner.refine_one_to_eight(tet_store,node_connectivity,tet_id);
550 : : break;
551 [ - - ]: 0 : case AMR::Refinement_Case::two_to_eight:
552 [ - - ]: 0 : round_two.insert( tet_store.get_parent_id(tet_id) );
553 : : //std::cout << "2->8\n";
554 : 0 : break;
555 [ - - ]: 0 : case AMR::Refinement_Case::four_to_eight:
556 [ - - ]: 0 : round_two.insert( tet_store.get_parent_id(tet_id));
557 : : //std::cout << "4->8\n";
558 : 0 : break;
559 : : case AMR::Refinement_Case::initial_grid:
560 : : // Do nothing
561 : : case AMR::Refinement_Case::none:
562 : : // Do nothing
563 : : break;
564 : : // No need for default as enum is explicitly covered
565 : : }
566 : : // Mark tet as not needing refinement
567 : 23419 : tet_store.marked_refinements.erase(tet_id);
568 : : }
569 : : }
570 : :
571 : : trace_out << "round_two size " << round_two.size() << std::endl;
572 [ - - ][ - + ]: 47 : for (const auto i : round_two)
573 : : {
574 : : trace_out << "round two i " << i << std::endl;
575 : :
576 : : // Cache children as we're about to change this data
577 [ - - ]: 0 : auto former_children = tet_store.data(i).children;
578 : :
579 : : AMR::Refinement_State& element = tet_store.data(i);
580 : :
581 [ - - ]: 0 : if (element.children.size() == 2)
582 : : {
583 : : trace_out << "perform 2:8" << std::endl;
584 [ - - ]: 0 : refiner.derefine_two_to_one(tet_store,node_connectivity,i);
585 : : }
586 [ - - ]: 0 : else if (element.children.size() == 4)
587 : : {
588 : : trace_out << "perform 4:8" << std::endl;
589 [ - - ]: 0 : refiner.derefine_four_to_one(tet_store,node_connectivity,i);
590 : : }
591 : : else {
592 : : std::cout << "num children " << element.children.size() << std::endl;
593 : : assert(0);
594 : : }
595 : : // remove tets and edges marked for deletion above
596 : 0 : refiner.delete_intermediates_of_children(tet_store);
597 : 0 : tet_store.process_delete_list();
598 : :
599 [ - - ]: 0 : refiner.refine_one_to_eight(tet_store,node_connectivity,i);
600 : :
601 : : // Grab children after it has been updated
602 [ - - ]: 0 : auto current_children = tet_store.data(i).children;
603 : :
604 : : // I want to set the children stored in *my* own children, to be
605 : : // the value of my new children....
606 : : //refiner.overwrite_children(tet_store, former_children, current_children);
607 : :
608 [ - - ]: 0 : tet_store.unset_marked_children(i); // FIXME: This will not work well in parallel
609 [ - - ]: 0 : element.refinement_case = AMR::Refinement_Case::one_to_eight;
610 : : }
611 : :
612 : : // Clean up dead edges
613 : : // clean_up_dead_edges(); // Nothing get's marked as "dead" atm?
614 : :
615 : : //std::cout << "Total Edges : " << tet_store.edge_store.size() << std::endl;
616 : : //std::cout << "Total Tets : " << tet_store.size() << std::endl;
617 : : //std::cout << "Total Nodes : " << m_x.size() << std::endl;
618 : :
619 : : trace_out << "Done ref" << std::endl;
620 : : node_connectivity.print();
621 : : node_connectivity.print();
622 [ + - ]: 47 : tet_store.print_node_types();
623 : : tet_store.print_tets();
624 : : //node_connectivity.print();
625 : :
626 : : //reset_intermediate_edges();
627 [ + - ]: 47 : remove_edge_locks(1);
628 [ + - ]: 47 : remove_normals();
629 : :
630 [ + - ]: 47 : lock_intermediates();
631 : :
632 [ + + ]: 194921 : for (auto& kv : tet_store.edge_store.edges) {
633 : : auto& local = kv.second;
634 [ + + ]: 194874 : if (local.needs_refining == 1) local.needs_refining = 0;
635 : : }
636 : 47 : }
637 : :
638 : 71 : void mesh_adapter_t::lock_intermediates()
639 : : {
640 : : /*
641 : : for (auto k : tet_store.intermediate_list)
642 : : {
643 : : refiner.lock_edges_from_node(tet_store,k, Edge_Lock_Case::intermediate);
644 : : }
645 : : */
646 : : // TODO: Passing tet_store twice probably isn't the best
647 : 71 : refiner.lock_intermediates(tet_store, tet_store.intermediate_list, Edge_Lock_Case::intermediate);
648 : 71 : }
649 : :
650 : : /**
651 : : * @brief A method implementing "Algorithm 1" from the paper
652 : : *
653 : : * @param num_to_refine Number of edges to refine
654 : : * @param tet_id The id of the given tet
655 : : */
656 : 50634 : void mesh_adapter_t::refinement_class_one(int num_to_refine, size_t tet_id)
657 : : {
658 : : trace_out << "Refinement Class One" << std::endl;
659 : :
660 : : // "If nrefine = 1
661 : : // Accept as a 1:2 refinement"
662 [ + + ]: 50634 : if (num_to_refine == 1)
663 : : {
664 : : tet_store.mark_one_to_two(tet_id);
665 : : }
666 : :
667 : : // "Else if nrefine = 2 OR nrefine = 3"
668 [ + + ]: 50105 : else if (num_to_refine > 1 && num_to_refine < 4)
669 : : {
670 : :
671 : : // We need to detect if the edges which need to refine are
672 : : // on the same face
673 : : // and if so which face so we know how to 1:4
674 : :
675 : 686 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
676 : : bool edges_on_same_face = false;
677 : : size_t face_refine_id = 0;
678 : :
679 : : // Iterate over each face
680 [ + + ]: 2099 : for (size_t face = 0; face < NUM_TET_FACES; face++)
681 : : {
682 : : int num_face_refine_edges = 0;
683 : 2018 : face_ids_t face_ids = face_list[face];
684 : :
685 : : trace_out << "Face is " <<
686 : : face_ids[0] << ", " <<
687 : : face_ids[1] << ", " <<
688 : : face_ids[2] << ", " <<
689 : : std::endl;
690 : :
691 : 2018 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
692 : : // For this face list, see which ones need refining
693 [ + + ]: 8072 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
694 : : {
695 : 6054 : edge_t key = face_edge_list[k];
696 [ + + ]: 6054 : if (tet_store.edge_store.get(key).needs_refining == 1)
697 : : {
698 : 3346 : num_face_refine_edges++;
699 : : }
700 : : }
701 [ + + ]: 2018 : if (num_face_refine_edges == num_to_refine)
702 : : {
703 : : edges_on_same_face = true;
704 : : face_refine_id = face;
705 : : trace_out << "Breaking with face value " << face << std::endl;
706 : 605 : break;
707 : : }
708 : : }
709 : :
710 : : // "If active edges are on the same face
711 : : // Activate any inactive edges of the face
712 : : // Accept as a 1:4 // refinement"
713 : : if (edges_on_same_face)
714 : : {
715 : : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list,
716 : : face_refine_id);
717 : :
718 : 605 : tet_t tet = tet_store.get(tet_id);
719 : : size_t opposite_id = tet[opposite_offset];
720 : :
721 : : trace_out << "face_refine_id " << face_refine_id << std::endl;
722 : : trace_out << "opposite_offset " << opposite_offset << std::endl;
723 : : trace_out << "opposite_id " << opposite_id << std::endl;
724 : :
725 : : // Activate edges on this face
726 : 605 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_list[face_refine_id]);
727 : :
728 [ + + ]: 2420 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
729 : : {
730 : 1815 : edge_t key = face_edge_list[k];
731 : 1815 : tet_store.edge_store.mark_for_refinement(key);
732 : : }
733 : :
734 : : //refiner.refine_one_to_four(tet_id, face_list[face_refine_id],
735 : : //opposite_id);
736 : : tet_store.mark_one_to_four(tet_id);
737 : : }
738 : : // "Else if active edges are not on the same face
739 : : // Activate all edges
740 : : // Accept as a 1:8 refinement"
741 : : else {
742 : : //refiner.refine_one_to_eight(tet_id);
743 : 81 : tet_store.mark_edges_for_refinement(tet_id);
744 : : tet_store.mark_one_to_eight(tet_id);
745 : : }
746 : :
747 : : }
748 : :
749 : : // "Else if nrefine > 3
750 : : // Activate any inactive edges
751 : : // Accept as a 1:8 refinement"
752 [ + - ]: 49419 : else if (num_to_refine > 3)
753 : : {
754 : : //refiner.refine_one_to_eight(tet_id);
755 : 49419 : tet_store.mark_edges_for_refinement(tet_id);
756 : : tet_store.mark_one_to_eight(tet_id);
757 : : }
758 : 50634 : }
759 : :
760 : : // TODO: Document this
761 : 0 : void mesh_adapter_t::lock_tet_edges(size_t tet_id) {
762 : 0 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
763 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
764 : : {
765 : 0 : edge_t key = edge_list[k];
766 [ - - ]: 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
767 : : {
768 : : trace_out << "LOCKING! " << key << std::endl;
769 : 0 : tet_store.edge_store.get(key).lock_case = AMR::Edge_Lock_Case::locked;
770 : : }
771 : : }
772 : 0 : }
773 : :
774 : : // TODO: Document this
775 : : // TODO: This has too similar a name to deactivate_tet
776 : 0 : void mesh_adapter_t::deactivate_tet_edges(size_t tet_id) {
777 : 0 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
778 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
779 : : {
780 : 0 : edge_t key = edge_list[k];
781 : :
782 : 0 : tet_store.edge_store.unmark_for_refinement(key);
783 : : trace_out << "Deactivating " << key << std::endl;
784 : 0 : tet_store.edge_store.get(key).needs_derefining = false;
785 : : }
786 : 0 : }
787 : :
788 : : /**
789 : : * @brief Unmarks edges of given tet for derefinement only
790 : : *
791 : : * @param tet_id The id of the given tet
792 : : */
793 : 25898 : void mesh_adapter_t::deactivate_deref_tet_edges(size_t tet_id) {
794 : 25898 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
795 [ + + ]: 181286 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
796 : : {
797 : 155388 : edge_t key = edge_list[k];
798 : :
799 : : trace_out << "Deactivating " << key << std::endl;
800 : 155388 : tet_store.edge_store.get(key).needs_derefining = false;
801 : : }
802 : 25898 : }
803 : :
804 : : /**
805 : : * @brief An implementation of "Algorithm 2" from the paper
806 : : *
807 : : * @param edge_list The list of edges for the given tet
808 : : * @param tet_id The id of the given tet
809 : : */
810 : 0 : void mesh_adapter_t::refinement_class_two(edge_list_t edge_list, size_t tet_id)
811 : : {
812 : : trace_out << "Refinement Class Two" << std::endl;
813 : :
814 : :
815 : : // "Deactivate all locked edges"
816 : :
817 : : // count number of active edges
818 : : int num_active_edges = 0;
819 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
820 : : {
821 : 0 : edge_t key = edge_list[k];
822 [ - - ]: 0 : if (tet_store.edge_store.get(key).lock_case != AMR::Edge_Lock_Case::unlocked)
823 : : {
824 : : tet_store.edge_store.unmark_for_refinement(key);
825 : : }
826 : : // "Count number of active edges"
827 [ - - ]: 0 : if (tet_store.edge_store.get(key).needs_refining == 1) {
828 : 0 : num_active_edges++;
829 : : }
830 : : }
831 : :
832 : : // Find out of two active edges live on the same face
833 : : bool face_refine = false;
834 : : size_t face_refine_id = 0; // FIXME: Does this need a better default
835 : 0 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
836 : :
837 : : // Iterate over each face
838 [ - - ]: 0 : for (size_t face = 0; face < NUM_TET_FACES; face++)
839 : : {
840 : : trace_out << "face " << face << std::endl;
841 : : int num_face_refine_edges = 0;
842 : : int num_face_locked_edges = 0;
843 : :
844 : 0 : face_ids_t face_ids = face_list[face];
845 : 0 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
846 : : // For this face list, see which ones need refining
847 [ - - ]: 0 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
848 : : {
849 : 0 : edge_t key = face_edge_list[k];
850 : : trace_out << "Checking " << key << std::endl;
851 [ - - ][ - - ]: 0 : if (tet_store.edge_store.get(key).needs_refining == 1)
852 : : {
853 : 0 : num_face_refine_edges++;
854 : : trace_out << "ref! " << key << std::endl;
855 : : }
856 : :
857 : : // Check for locked edges
858 : : // This case only cares about faces with no locks
859 [ - - ][ - - ]: 0 : if (tet_store.edge_store.get(key).lock_case != AMR::Edge_Lock_Case::unlocked)
860 : : {
861 : 0 : num_face_locked_edges++;
862 : : trace_out << "locked! " << key << std::endl;
863 : : }
864 : : }
865 : :
866 : :
867 : : // Decide if we want to process this face
868 [ - - ]: 0 : if (num_face_refine_edges >= 2 && num_face_locked_edges == 0)
869 : : {
870 : : // We can refine this face
871 : : face_refine = true;
872 : : face_refine_id = face;
873 : 0 : break;
874 : : }
875 : : }
876 : :
877 : : // "If nrefine = 1
878 : : // Accept as 1:2 refinement"
879 : : // TODO: can we hoist this higher
880 [ - - ]: 0 : if (num_active_edges == 1)
881 : : {
882 : : //node_pair_t nodes = find_single_refinement_nodes(edge_list);
883 : : //refine_one_to_two( tet_id, nodes[0], nodes[1]);
884 : : tet_store.mark_one_to_two(tet_id);
885 : : }
886 : : // "Else if any face has nrefine >= 2 AND no locked edges
887 : : // Active any inactive edges of the face
888 : : // Accept as a 1:4 refinement"
889 [ - - ]: 0 : else if (face_refine)
890 : : {
891 : : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list, face_refine_id);
892 : :
893 : 0 : tet_t tet = tet_store.get(tet_id);
894 : : size_t opposite_id = tet[opposite_offset];
895 : :
896 : : trace_out << "Tet ID " << tet_id << std::endl;
897 : : trace_out << "Opposite offset " << opposite_offset << std::endl;
898 : : trace_out << "Opposite id " << opposite_id << std::endl;
899 : : trace_out << "Face refine id " << face_refine_id << std::endl;
900 : :
901 : : edge_list_t face_edge_list =
902 : 0 : AMR::edge_store_t::generate_keys_from_face_ids(face_list[face_refine_id]);
903 : :
904 [ - - ]: 0 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
905 : : {
906 : 0 : edge_t key = face_edge_list[k];
907 : 0 : tet_store.edge_store.mark_for_refinement(key);
908 : : }
909 : :
910 : : //refiner.refine_one_to_four(tet_id, face_list[face_refine_id], opposite_id);
911 : : tet_store.mark_one_to_four(tet_id);
912 : : }
913 : :
914 : : // "Else
915 : : // Deactivate all edges
916 : : // Mark all edges as locked"
917 : : else {
918 : : trace_out << "Class 2 causes some locking.." << std::endl;
919 : 0 : deactivate_tet_edges(tet_id);
920 : 0 : lock_tet_edges(tet_id);
921 : : }
922 : :
923 : 0 : }
924 : :
925 : : /**
926 : : * @brief Based on a tet_id, decide if it's current state of locked
927 : : * and marked edges maps to a valid refinement case. The logic for
928 : : * this was derived from talking to JW and reading Chicoma.
929 : : *
930 : : * It basically just checks if something a 1:2 and has 3
931 : : * intermediates and 3 makred edges, or is a 1:4 and has 5/6
932 : : * intermediates
933 : : *
934 : : * @param child_id the id of the tet to check
935 : : *
936 : : * @return A bool saying if the tet is in a valid state to be refined
937 : : */
938 : 0 : bool mesh_adapter_t::check_valid_refinement_case(size_t child_id) {
939 : :
940 : : trace_out << "check valid ref " << child_id << std::endl;
941 : 0 : edge_list_t edge_list = tet_store.generate_edge_keys(child_id);
942 : :
943 : : size_t num_to_refine = 0;
944 : : size_t num_intermediate = 0;
945 : : size_t unlocked = 0;
946 : : size_t locked = 0;
947 : :
948 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
949 : : {
950 : 0 : edge_t key = edge_list[k];
951 : : trace_out << "Key " << key << std::endl;
952 : :
953 : : // Count intermediate edges
954 : 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::intermediate
955 [ - - ][ - - ]: 0 : || tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::temporary)
956 : : {
957 : : trace_out << "found intermediate" << std::endl;
958 : 0 : num_intermediate++;
959 : : }
960 : :
961 : : // Count number of marked for refinement
962 [ - - ]: 0 : if (tet_store.edge_store.get(key).needs_refining == 1)
963 : : {
964 : : trace_out << "found refine" << std::endl;
965 : 0 : num_to_refine++;
966 : : }
967 : :
968 : 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
969 : : {
970 : : trace_out << "found unlocked" << std::endl;
971 : : unlocked++;
972 : : }
973 : :
974 : 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::locked)
975 : : {
976 : : trace_out << "found locked" << std::endl;
977 : : locked++;
978 : : }
979 : :
980 : : }
981 : :
982 : : AMR::Refinement_State& element = tet_store.data(child_id);
983 : :
984 : : trace_out <<
985 : : "Intermediates " << num_intermediate <<
986 : : " num to refine " << num_to_refine <<
987 : : " unlocked " << unlocked <<
988 : : " locked " << locked <<
989 : : " Case " << element.refinement_case <<
990 : : std::endl;
991 : :
992 : : // check if element is 1:2
993 [ - - ]: 0 : if (element.refinement_case == AMR::Refinement_Case::one_to_two)
994 : : {
995 : : // If so check it has 3 intermediates and 3 which need refining
996 [ - - ]: 0 : if (num_intermediate != 3 || num_to_refine != 3) {
997 : 0 : return false;
998 : : }
999 : : else {
1000 : : trace_out << "True " <<
1001 : : "Intermediates " << num_intermediate <<
1002 : : " num to refine " << num_to_refine <<
1003 : : " Case " << element.refinement_case <<
1004 : : " 2:1 " << AMR::Refinement_Case::one_to_two <<
1005 : : std::endl;
1006 : : }
1007 : : }
1008 : :
1009 : : // check if element is 1:4
1010 [ - - ]: 0 : else if (element.refinement_case == AMR::Refinement_Case::one_to_four)
1011 : : {
1012 : : // TODO: Check if it's a center tet for a 1:4
1013 : : // FIXME: Is this even needed? How else would you get these
1014 : : // combinations? Can't we just combine these two checks?
1015 : :
1016 : : bool is_center_tet = tet_store.is_center(child_id);
1017 : :
1018 : : if (is_center_tet)
1019 : : {
1020 [ - - ]: 0 : if (num_to_refine != 0 || num_intermediate != 6)
1021 : : {
1022 : : trace_out << "Fail compat 1:4 center" << std::endl;
1023 : 0 : return false;
1024 : : }
1025 : : }
1026 : : else { // Is one of the outsides (not center)
1027 [ - - ]: 0 : if (num_to_refine != 1 || num_intermediate != 5)
1028 : : {
1029 : : trace_out << "Fail compat 1:4 non center" << std::endl;
1030 : 0 : return false;
1031 : : }
1032 : : }
1033 : : }
1034 : :
1035 : : // If it makes it here, it's compatible
1036 : : return true;
1037 : :
1038 : : }
1039 : :
1040 : : /**
1041 : : * @brief Place holder method for the implementation of "Algorithm
1042 : : * 3" from the paper
1043 : : */
1044 : : // TODO: Does this parse a childs siblings multiple times?
1045 : 0 : void mesh_adapter_t::refinement_class_three(size_t tet_id) {
1046 : :
1047 : : trace_out << "Refinement Class Three" << std::endl;
1048 : :
1049 : : // "Identify parent element iparent"
1050 : : // TODO: WE should either always use the id to fetch, or always do the data lookup
1051 : : //size_t parent_id = master_elements.get_parent(tet_id);
1052 : : size_t parent_id = tet_store.get_parent_id(tet_id);
1053 : :
1054 : : trace_out << "Parent id = " << parent_id << std::endl;
1055 : :
1056 : : // NOTE: This implies comms when we use these ids?
1057 : 0 : child_id_list_t children = tet_store.data(parent_id).children;
1058 : :
1059 : : // "Do for each child element ielement
1060 : : // Activate all non-locked edges
1061 : : // Deactivate all locked edges"
1062 [ - - ]: 0 : for (size_t i = 0; i < children.size(); i++)
1063 : : {
1064 : : // TODO: Is this in element or tet ids?
1065 : : trace_out << "Checking child " << children[i] << std::endl;
1066 [ - - ]: 0 : edge_list_t edge_list = tet_store.generate_edge_keys(children[i]);
1067 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1068 : : {
1069 : 0 : edge_t key = edge_list[k];
1070 : : trace_out << "Compat 3 " << key << std::endl;
1071 [ - - ][ - - ]: 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
1072 : : {
1073 : : trace_out << "Compat 3 marking edge " << key << std::endl;
1074 : : tet_store.edge_store.mark_for_refinement(key);
1075 : : }
1076 : : else {
1077 : : tet_store.edge_store.unmark_for_refinement(key);
1078 : : }
1079 : : }
1080 : : }
1081 : :
1082 : : // "Set compatible = TRUE
1083 : : bool compatible = true;
1084 : :
1085 : : // Do for each child element ielement
1086 : : // If ielement is not a valid refinement case
1087 : : // compatible = FALSE"
1088 [ - - ]: 0 : for (size_t i = 0; i < children.size(); i++)
1089 : : {
1090 : 0 : size_t child = children[i];
1091 [ - - ][ - - ]: 0 : if ( !check_valid_refinement_case(child) )
1092 : : {
1093 : : trace_out << "Compat 3 Marking compatible false because of invalid refinement case" << std::endl;
1094 : :
1095 : : compatible = false;
1096 : : }
1097 : : else {
1098 : : trace_out << "Is compatible" << std::endl;
1099 : : }
1100 : : }
1101 : :
1102 : : // "If compatible = FALSE
1103 : : // Do for each child element ielement
1104 : : // Deactive all edges of ielement
1105 : : // Mark all edges of ielement as locked
1106 : : // Mark ielement as normal"
1107 [ - - ]: 0 : if (compatible == false)
1108 : : {
1109 [ - - ]: 0 : for (size_t i = 0; i < children.size(); i++)
1110 : : {
1111 : 0 : size_t child = children[i];
1112 [ - - ]: 0 : deactivate_tet_edges(child);
1113 [ - - ]: 0 : lock_tet_edges(child);
1114 : : trace_out << "Compat 3 locking edges of " << child << std::endl;
1115 : : // Here we interpret normal to mean "don't treat it like it has intermediates"
1116 : : tet_store.mark_normal(child);
1117 : : trace_out << "Compat 3 " << child << std::endl;
1118 : : }
1119 : : }
1120 : : else {
1121 : : trace_out << "TIME TO 2:8 " << tet_id << std::endl;
1122 : : // Accept as 2:8 or 4:8
1123 : : AMR::Refinement_State& element = tet_store.data(tet_id);
1124 [ - - ]: 0 : if (element.refinement_case == AMR::Refinement_Case::one_to_two)
1125 : : {
1126 : : tet_store.mark_two_to_eight(tet_id);
1127 : : }
1128 [ - - ]: 0 : else if (element.refinement_case == AMR::Refinement_Case::one_to_four)
1129 : : {
1130 : : tet_store.mark_four_to_eight(tet_id);
1131 : : }
1132 : : else {
1133 : : trace_out << " I don't know what to do with this..it looks like you're trying to 2/4:8 an 8... " << std::endl;
1134 : : }
1135 : :
1136 : : }
1137 : 0 : }
1138 : :
1139 : 47 : void mesh_adapter_t::remove_normals()
1140 : : {
1141 [ + + ]: 137891 : for (const auto& kv : tet_store.tets)
1142 : : {
1143 : 137844 : size_t tet_id = kv.first;
1144 : : tet_store.set_normal(tet_id, 0);
1145 : : }
1146 : 47 : }
1147 : :
1148 : 118 : void mesh_adapter_t::remove_edge_locks(int intermediate)
1149 : : {
1150 [ + + ]: 283420 : for (const auto& kv : tet_store.tets)
1151 : : {
1152 : 283302 : size_t tet_id = kv.first;
1153 : :
1154 : : trace_out << "Process tet removelock " << tet_id << std::endl;
1155 : :
1156 : : // Only apply checks to tets on the active list
1157 : : if (tet_store.is_active(tet_id)) {
1158 : : // change it from intermediate to locked
1159 : 249425 : update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::locked, AMR::Edge_Lock_Case::unlocked);
1160 [ + - ]: 249425 : if (intermediate) {
1161 : 249425 : update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::intermediate, AMR::Edge_Lock_Case::unlocked);
1162 : : }
1163 : : }
1164 : : }
1165 : :
1166 : 118 : }
1167 : : //void mesh_adapter_t::reset_intermediate_edges()
1168 : : //{
1169 : : // for (const auto& kv : tet_store.tets)
1170 : : // {
1171 : : // size_t tet_id = kv.first;
1172 : :
1173 : : // trace_out << "Process tet reset " << tet_id << std::endl;
1174 : :
1175 : : // // Only apply checks to tets on the active list
1176 : : // if (tet_store.is_active(tet_id)) {
1177 : : // // change it from intermediate to locked
1178 : : // update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::intermediate, AMR::Edge_Lock_Case::locked);
1179 : : // }
1180 : : // }
1181 : : //}
1182 : :
1183 : 627017 : void mesh_adapter_t::update_tet_edges_lock_type(size_t tet_id, AMR::Edge_Lock_Case check, AMR::Edge_Lock_Case new_case) {
1184 : 627017 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
1185 [ + + ]: 4389119 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1186 : : {
1187 : 3762102 : edge_t key = edge_list[k];
1188 [ + + ]: 3762102 : if (tet_store.edge_store.get(key).lock_case == check)
1189 : : {
1190 : 3796 : tet_store.edge_store.get(key).lock_case = new_case;
1191 : : }
1192 : : }
1193 : 627017 : }
1194 : :
1195 : : /**
1196 : : * @brief This unlocks edges that were previously locked with a `temporary'
1197 : : * lock, indicating a parallel compatibility induced locking
1198 : : */
1199 : 71 : void mesh_adapter_t::remove_edge_temp_locks()
1200 : : {
1201 [ + + ]: 145529 : for (const auto& kv : tet_store.tets)
1202 : : {
1203 : 145458 : size_t tet_id = kv.first;
1204 : :
1205 : : trace_out << "Process tet remove temp lock " << tet_id << std::endl;
1206 : :
1207 : : // Only apply checks to tets on the active list
1208 : : if (tet_store.is_active(tet_id)) {
1209 : : // change it from temporary to unlocked
1210 : 128167 : update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::temporary,
1211 : : AMR::Edge_Lock_Case::unlocked);
1212 : : }
1213 : : }
1214 : 71 : }
1215 : :
1216 : 40 : void mesh_adapter_t::mark_derefinement()
1217 : : {
1218 : : const size_t max_num_rounds = AMR_MAX_ROUNDS;
1219 : :
1220 : : // Mark refinements
1221 : : size_t iter;
1222 : : //Iterate until convergence
1223 [ + - ]: 65 : for (iter = 0; iter < max_num_rounds; iter++)
1224 : : {
1225 : 65 : tet_store.marked_derefinements.get_state_changed() = false;
1226 : :
1227 : : // set of elements which have been considered for derefinement
1228 : : std::unordered_set< size_t > done_deref_marking;
1229 : :
1230 : : // Loop over tets
1231 [ + + ]: 132964 : for (const auto& kv : tet_store.tets)
1232 : : {
1233 : : // this loop only runs for active tets
1234 : 132899 : if (!tet_store.is_active(kv.first)) {
1235 [ + - ]: 16442 : deactivate_deref_tet_edges(kv.first);
1236 : 118687 : continue;
1237 : : }
1238 : : size_t activetet_id = kv.first;
1239 : :
1240 : : // check if activetet_id has a parent (assign to tet_id)
1241 : : // if it does not, activetet_id is not a derefinement candidate
1242 : : size_t tet_id;
1243 : : const auto& activetet_data = tet_store.data(activetet_id);
1244 [ + + ]: 116457 : if (!activetet_data.has_parent) {
1245 [ + - ]: 1551 : deactivate_deref_tet_edges(activetet_id);
1246 : 1551 : continue;
1247 : : }
1248 : : else {
1249 [ + + ]: 114906 : tet_id = activetet_data.parent_id;
1250 : : }
1251 : :
1252 : : // if already considered for deref, do not reconsider
1253 : 99949 : if (done_deref_marking.count(tet_id) > 0) {
1254 : 99949 : continue;
1255 : : }
1256 : : done_deref_marking.insert(tet_id);
1257 : :
1258 [ + - ]: 14957 : child_id_list_t children = tet_store.data(tet_id).children;
1259 : :
1260 : : // check if any child of tet_id (i.e. any active tet) is marked
1261 : : // for refinement
1262 : : bool is_child_ref(false);
1263 [ + + ]: 130337 : for (size_t i=0; i<children.size(); i++) {
1264 [ + - ]: 115380 : edge_list_t chedge_list = tet_store.generate_edge_keys(children[i]);
1265 : : // Check each edge, see if it is marked for refinement
1266 [ + + ]: 807660 : for (size_t k=0; k<NUM_TET_EDGES; k++) {
1267 : 692280 : edge_t edge = chedge_list[k];
1268 : :
1269 [ + - ][ + + ]: 692280 : if (tet_store.edge_store.get(edge).needs_refining == 1) {
1270 : : is_child_ref = true;
1271 : 30120 : continue;
1272 : : }
1273 : : }
1274 : : }
1275 : : // deactivate from deref if marked for ref
1276 [ + + ]: 14957 : if (is_child_ref) {
1277 : : trace_out << tet_id << " Looping cancelled since child marked for refinement." << std::endl;
1278 [ + + ]: 6705 : for (auto child_id : children) {
1279 [ + - ]: 5960 : deactivate_deref_tet_edges(child_id);
1280 : : }
1281 [ + - ]: 745 : deactivate_deref_tet_edges(tet_id);
1282 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1283 : : continue;
1284 [ + - ]: 745 : }
1285 : :
1286 : : // check if tet_id has been marked for deref-ref
1287 [ + - ]: 14212 : edge_list_t pedge_list = tet_store.generate_edge_keys(tet_id);
1288 : : // Check each edge, see if it is marked for refinement
1289 [ + + ]: 99484 : for (size_t k=0; k<NUM_TET_EDGES; k++) {
1290 : 85272 : edge_t edge = pedge_list[k];
1291 : :
1292 : : // deactivate child-edges from deref if marked '2'
1293 [ + - ][ + + ]: 85272 : if (tet_store.edge_store.get(edge).needs_refining == 2) {
1294 [ + - ]: 236 : auto edge_nodes = edge.get_data();
1295 : 236 : auto ch_node = node_connectivity.data().at(edge_nodes);
1296 : : std::array< edge_t, 2> ch_edge;
1297 [ - + ]: 236 : ch_edge[0] = {edge_nodes[0], ch_node};
1298 [ - + ]: 236 : ch_edge[1] = {edge_nodes[1], ch_node};
1299 [ + - ]: 236 : tet_store.edge_store.get(ch_edge[0]).needs_derefining = 0;
1300 [ + - ]: 236 : tet_store.edge_store.get(ch_edge[1]).needs_derefining = 0;
1301 : : }
1302 : : }
1303 : :
1304 : : // This is useful for later inspection
1305 : : //edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
1306 : : std::size_t num_to_derefine = 0; // Nodes
1307 : :
1308 : : AMR::Refinement_Case refinement_case = tet_store.get_refinement_case(tet_id);
1309 : :
1310 [ + - ]: 14212 : auto derefine_node_set = refiner.find_derefine_node_set(tet_store, tet_id);
1311 : : // Find the set of nodes which are not in the parent
1312 : : std::unordered_set<size_t> non_parent_nodes =
1313 [ + - ]: 14212 : refiner.child_exclusive_nodes(tet_store, tet_id);
1314 : :
1315 : : //for (auto drnode: derefine_node_set)
1316 : : // trace_out << "derefine node: " << drnode << std::endl;
1317 : :
1318 : : num_to_derefine = derefine_node_set.size();
1319 : :
1320 : : if (num_to_derefine > 0) {
1321 : : trace_out << "num_to_derefine " << num_to_derefine << std::endl;
1322 : : trace_out << "ref_case " << refinement_case << std::endl;
1323 : : trace_out << "num children " << children.size() << std::endl;
1324 : : }
1325 : :
1326 : : //num_to_derefine = convert_derefine_edges_to_points(tet_store, tet_id, num_edges_to_derefine, refinement_case);
1327 : :
1328 : : // "If nderefine = 1
1329 [ + + ]: 14212 : if (num_to_derefine == 1)
1330 : : {
1331 : : // If icase = 1:2
1332 : :
1333 : : //if (refinement_case == AMR::Refinement_Case::one_to_two)
1334 [ + + ]: 390 : if (children.size() == 2)
1335 : : {
1336 : : // Accept as 2:1 derefine"
1337 : : trace_out << "Accept as 2:1" << std::endl;
1338 : : //refiner.derefine_two_to_one(tet_store, node_connectivity, tet_id);
1339 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::two_to_one);
1340 : : }
1341 : : // "Else
1342 : : else {
1343 : : // Deactivate all points"
1344 [ + + ]: 108 : for (auto child_id : children) {
1345 [ + - ]: 96 : deactivate_deref_tet_edges(child_id);
1346 : : }
1347 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1348 : : trace_out << "giving up on deref decision. deactivate near 2:1 ntd = 1" << std::endl;
1349 : : }
1350 : : }
1351 : :
1352 : :
1353 : : // "If nderefine = 2
1354 : : else if (num_to_derefine == 2)
1355 : : {
1356 : : // If icase = 1:4
1357 : : //if (refinement_case == AMR::Refinement_Case::one_to_four)
1358 [ - + ]: 14 : if (children.size() == 4)
1359 : : {
1360 : : // Accept as 4:2 derefine"
1361 : : trace_out << "Accept as 4:2" << std::endl;
1362 : : //refiner.derefine_four_to_two(tet_store, node_connectivity, tet_id);
1363 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::four_to_two);
1364 : : }
1365 : : // "Else
1366 : : else {
1367 : : // Deactivate all points"
1368 [ + + ]: 126 : for (auto child_id : children) {
1369 [ + - ]: 112 : deactivate_deref_tet_edges(child_id);
1370 : : }
1371 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1372 : : trace_out << "giving up on deref decision. deactivate near 4:2 ntd = 2" << std::endl;
1373 : : }
1374 : : }
1375 : :
1376 : : // "If nderefine = 3
1377 : : else if (num_to_derefine == 3)
1378 : : {
1379 : : // If icase = 1:4
1380 : : //if (refinement_case == AMR::Refinement_Case::one_to_four)
1381 [ + + ]: 522 : if (children.size() == 4)
1382 : : {
1383 : : // Accept as 4:1 derefine"
1384 : : trace_out << "Accept as 4:1" << std::endl;
1385 : : //refiner.derefine_four_to_one(tet_store, node_connectivity, tet_id);
1386 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::four_to_one);
1387 : : }
1388 : : // "Else if icase = 1:8
1389 : : //else if (refinement_case == AMR::Refinement_Case::one_to_eight)
1390 [ + - ]: 20 : else if (children.size() == 8)
1391 : : {
1392 : : // we have a list of (non-parent) nodes that is marked
1393 : : // for derefinement. First, determine the nodes that are
1394 : : // unmarked for derefinement (or inactive_nodes). Then,
1395 : : // determine if these are on a single face.
1396 : : std::unordered_set<size_t> inactive_node_set;
1397 [ + + ][ + + ]: 140 : for (auto npn : non_parent_nodes) {
1398 : : if (derefine_node_set.count(npn) == 0)
1399 : : inactive_node_set.insert(npn);
1400 : : }
1401 : : Assert(inactive_node_set.size() == 3, "Incorrectly "
1402 : : "sized inactive-node set");
1403 [ + - ]: 20 : auto same_face = check_same_face(tet_id, inactive_node_set);
1404 : : // If inactive points lie on same face
1405 [ + + ]: 20 : if (same_face.first == true)
1406 : : {
1407 : : // Accept as 8:4 derefinement
1408 : : trace_out << "Accept as 8:4" << std::endl;
1409 : :
1410 : : // create a vector of node-array-pairs to mark edges
1411 : : // for refinement 1:4
1412 : : std::vector< std::array< std::size_t, 2 > > ref_edges;
1413 : : trace_out << "inactive nodes on same face: ";
1414 [ + - ][ + + ]: 56 : for (auto n:inactive_node_set) {
1415 : : trace_out << n << ", ";
1416 [ + - ]: 42 : ref_edges.push_back(node_connectivity.get(n));
1417 : : }
1418 : : trace_out << std::endl;
1419 : :
1420 [ + - ][ + - ]: 28 : tet_store.edge_store.mark_edges_for_deref_ref(ref_edges);
[ + - ][ - - ]
1421 : : //refiner.derefine_eight_to_four(tet_store, node_connectivity, tet_id);
1422 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_four);
1423 : : }
1424 : : // "Else
1425 : : else {
1426 : : // Deactivate all points"
1427 [ + + ]: 54 : for (auto child_id : children) {
1428 [ + - ]: 48 : deactivate_deref_tet_edges(child_id);
1429 : : }
1430 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1431 : : trace_out << "giving up on deref decision. deactivate near 8:4 ntd = 3" << std::endl;
1432 : : }
1433 : :
1434 : : }
1435 : : }
1436 : :
1437 : : // "If nderefine = 4
1438 : : else if (num_to_derefine == 4)
1439 : : //else if (children.size() == 4)
1440 : : {
1441 : : // we have a list of (non-parent) nodes that is marked
1442 : : // for derefinement. First, determine the nodes that are
1443 : : // unmarked for derefinement (or inactive_nodes). Then,
1444 : : // determine if these are on a single face.
1445 : : std::unordered_set<size_t> inactive_node_set;
1446 [ + + ][ + + ]: 14 : for (auto npn : non_parent_nodes) {
1447 : : if (derefine_node_set.count(npn) == 0)
1448 : : inactive_node_set.insert(npn);
1449 : : }
1450 : : Assert(inactive_node_set.size() == 2, "Incorrectly "
1451 : : "sized inactive-node set");
1452 : : // Check if the inactive point belong to the same parent
1453 : : // face and deactivate the third point on that face
1454 [ + - ]: 2 : auto same_face = check_same_face(tet_id, inactive_node_set);
1455 : :
1456 [ + - ]: 2 : if (same_face.first == true)
1457 : : {
1458 : : // deactivate the edges associated with same_face.second
1459 [ + + ]: 18 : for (size_t i = 0; i < children.size(); i++)
1460 : : {
1461 [ + - ]: 16 : edge_list_t edge_list = tet_store.generate_edge_keys(children[i]);
1462 [ + + ]: 112 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1463 : : {
1464 : 96 : edge_t edge = edge_list[k];
1465 : : size_t A = edge.first();
1466 : : size_t B = edge.second();
1467 [ + + ][ + + ]: 96 : if (A == same_face.second || B == same_face.second)
1468 [ + - ]: 36 : tet_store.edge_store.get(edge).needs_derefining = false;
1469 : : }
1470 : : }
1471 : :
1472 : : // create a vector of node-array-pairs to mark edges
1473 : : // for refinement 1:4
1474 : : inactive_node_set.insert(same_face.second);
1475 : : std::vector< std::array< std::size_t, 2 > > ref_edges;
1476 [ + - ][ + + ]: 8 : for (auto n:inactive_node_set) {
1477 [ + - ]: 6 : ref_edges.push_back(node_connectivity.get(n));
1478 : : }
1479 : :
1480 [ + - ][ + - ]: 4 : tet_store.edge_store.mark_edges_for_deref_ref(ref_edges);
[ + - ][ - - ]
1481 : :
1482 : : // Accept as 8:4 derefinement
1483 : : trace_out << "Accept as 8:4" << std::endl;
1484 : : //refiner.derefine_eight_to_four(tet_store, node_connectivity, tet_id);
1485 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_four);
1486 : : }
1487 : : // "Else
1488 : : else {
1489 : : // Deactivate all points"
1490 [ - - ]: 0 : for (auto child_id : children) {
1491 [ - - ]: 0 : deactivate_deref_tet_edges(child_id);
1492 : : }
1493 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1494 : : trace_out << "giving up on deref decision. deactivate near 8:4 ntd = 4" << std::endl;
1495 : : }
1496 : : }
1497 : :
1498 : : // "If nderefine = 5
1499 : : else if (num_to_derefine == 5)
1500 : : {
1501 : : // Accept as 8:2 derefine"
1502 : : trace_out << "Accept as 8:2 " << std::endl;
1503 : : //refiner.derefine_eight_to_two(tet_store, node_connectivity, tet_id);
1504 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_two);
1505 : : }
1506 : :
1507 : : // "If nderefine = 6
1508 : : else if (num_to_derefine == 6)
1509 : : {
1510 : : // Accept as 8:1 derefine"
1511 : : trace_out << "Accept as 8:1" << std::endl;
1512 : : //refiner.derefine_eight_to_one(tet_store, node_connectivity, tet_id);
1513 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_one);
1514 : : }
1515 : :
1516 : : // "If nderefine = 0
1517 : : else {
1518 : : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1519 : : // Deactivate all points"
1520 [ + + ]: 1062 : for (auto child_id : children) {
1521 [ + - ]: 944 : deactivate_deref_tet_edges(child_id);
1522 : : }
1523 : : trace_out << "giving up with no deref decision because nderefine = 0" << std::endl;
1524 : : }
1525 : : }
1526 : :
1527 : : // If nothing changed during that round, break
1528 [ + + ]: 65 : if (!tet_store.marked_derefinements.get_state_changed())
1529 : : {
1530 : : trace_out << "Terminating loop at iter " << iter << std::endl;
1531 : : break;
1532 : : }
1533 : : trace_out << "End iter " << iter << std::endl;
1534 : : // clear out set of elements considered during this iteration
1535 : : done_deref_marking.clear();
1536 : : }
1537 : : trace_out << "Deref Loop took " << iter << " rounds." << std::endl;
1538 : 40 : }
1539 : :
1540 : : // TODO: document
1541 : 24 : void mesh_adapter_t::perform_derefinement()
1542 : : {
1543 : : trace_out << "Perform deref" << std::endl;
1544 : :
1545 : : // Do derefinements
1546 [ + + ]: 7638 : for (const auto& kv : tet_store.tets)
1547 : : {
1548 [ + + ]: 7614 : size_t tet_id = kv.first;
1549 : : //size_t parent_id = 0;
1550 : :
1551 : : // TODO: Do I really want to loop all tets?
1552 : :
1553 : : // TODO: is this doing a double lookup?
1554 : : if (tet_store.has_derefinement_decision(tet_id))
1555 : : {
1556 : : trace_out << "Do derefine of " << tet_id << std::endl;
1557 : : //size_t parent_id = tet_store.get_parent_id(tet_id);
1558 : : //trace_out << "Parent = " << parent_id << std::endl;
1559 [ + + ][ - + ]: 6562 : switch(tet_store.marked_derefinements.get(tet_id))
[ - - ][ + ]
1560 : : {
1561 : 130 : case AMR::Derefinement_Case::two_to_one:
1562 : 130 : refiner.derefine_two_to_one(tet_store,node_connectivity,tet_id);
1563 : : trace_out << "Completed derefine 2:1 of " << tet_id << std::endl;
1564 : : break;
1565 : 176 : case AMR::Derefinement_Case::four_to_one:
1566 : 176 : refiner.derefine_four_to_one(tet_store,node_connectivity,tet_id);
1567 : : trace_out << "Completed derefine 4:1 of " << tet_id << std::endl;
1568 : : break;
1569 : 0 : case AMR::Derefinement_Case::four_to_two:
1570 : 0 : refiner.derefine_four_to_two(tet_store,node_connectivity,tet_id);
1571 : : trace_out << "Completed derefine 4:2 of " << tet_id << std::endl;
1572 : : break;
1573 : 5974 : case AMR::Derefinement_Case::eight_to_one:
1574 : 5974 : refiner.derefine_eight_to_one(tet_store,node_connectivity,tet_id);
1575 : : trace_out << "Completed derefine 8:1 of " << tet_id << std::endl;
1576 : : break;
1577 : 0 : case AMR::Derefinement_Case::eight_to_two:
1578 : 0 : refiner.derefine_eight_to_two(tet_store,node_connectivity,tet_id);
1579 : : trace_out << "Completed derefine 8:2 of " << tet_id << std::endl;
1580 : : break;
1581 : 0 : case AMR::Derefinement_Case::eight_to_four:
1582 : 0 : refiner.derefine_eight_to_four(tet_store,node_connectivity,tet_id);
1583 : : trace_out << "Completed derefine 8:4 of " << tet_id << std::endl;
1584 : : break;
1585 : : case AMR::Derefinement_Case::skip:
1586 : : // What do we do with skip?
1587 : : break;
1588 : : }
1589 : : // Mark tet as not needing derefinement
1590 : 6562 : tet_store.marked_derefinements.erase(tet_id);
1591 : : }
1592 : : }
1593 : :
1594 : : node_connectivity.print();
1595 : 24 : refiner.delete_intermediates_of_children(tet_store);
1596 : 24 : tet_store.process_delete_list();
1597 : 24 : tet_store.print_node_types();
1598 : :
1599 : 24 : lock_intermediates();
1600 : :
1601 [ + + ]: 130086 : for (auto& kv : tet_store.edge_store.edges) {
1602 : : auto& local = kv.second;
1603 : 130062 : local.needs_derefining = 0;
1604 [ + + ]: 130062 : if (local.needs_refining == 2) local.needs_refining = 0;
1605 : : }
1606 : 24 : }
1607 : :
1608 : : }
1609 : :
1610 : : #if defined(__clang__)
1611 : : #pragma clang diagnostic pop
1612 : : #endif
|
