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 : 17 : 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 [ + - ]: 17 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
55 : :
56 [ + + ][ - + ]: 17 : 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 : 17 : edges_on_face[0][0] =
68 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[0].get_data());
69 : 17 : edges_on_face[0][1] =
70 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[1].get_data());
71 : 17 : edges_on_face[0][2] =
72 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[3].get_data());
73 : :
74 : : // A-B-D
75 : 17 : edges_on_face[1][0] =
76 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[0].get_data());
77 : 17 : edges_on_face[1][1] =
78 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[2].get_data());
79 : 17 : edges_on_face[1][2] =
80 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[4].get_data());
81 : :
82 : : // B-C-D
83 : 17 : edges_on_face[2][0] =
84 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[3].get_data());
85 : 17 : edges_on_face[2][1] =
86 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[4].get_data());
87 : 17 : edges_on_face[2][2] =
88 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[5].get_data());
89 : :
90 : : // A-C-D
91 : 17 : edges_on_face[3][0] =
92 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[1].get_data());
93 : 17 : edges_on_face[3][1] =
94 [ + - ]: 17 : tk::cref_find(node_connectivity.data(),edge_list[2].get_data());
95 : 17 : edges_on_face[3][2] =
96 [ + - ]: 17 : 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 : 17 : bool same_face(false);
100 : 17 : [[maybe_unused]] bool tnode_set(false);
101 : 17 : std::size_t third_node = 0;
102 [ + + ]: 85 : for(const auto& face : edges_on_face)
103 : : {
104 : 68 : std::size_t icount = 0;
105 [ + + ]: 272 : for (const auto& np_node : face) {
106 [ + - ][ + + ]: 204 : if (inactive_nodes.count(np_node)) ++icount;
107 : : }
108 [ + + ]: 68 : if (inactive_nodes.size() == icount) {
109 : 12 : 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 [ + + ]: 12 : if (inactive_nodes.size() == 2) {
113 [ + - ]: 6 : for (auto fn:face) {
114 [ + - ][ + + ]: 6 : if (inactive_nodes.count(fn) == 0) {
115 : 4 : third_node = fn;
116 : 4 : tnode_set = true;
117 : 4 : break;
118 : : }
119 : : }
120 : : }
121 : : }
122 : : }
123 : :
124 [ + + ][ + + ]: 17 : if (same_face && inactive_nodes.size() == 2)
[ + + ]
125 [ - + ][ - - ]: 4 : Assert(tnode_set, "Third node on face not set in derefine");
[ - - ][ - - ]
126 : :
127 : 17 : 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 : 2128 : void mesh_adapter_t::consume_tets(const std::vector< std::size_t >& tetinpoel )
137 : : {
138 [ + + ]: 2341405 : for (size_t i = 0; i < tetinpoel.size(); i+=4)
139 : : {
140 : : tet_t t = {
141 : : {
142 : 2339277 : tetinpoel[i],
143 : 2339277 : tetinpoel[i+1],
144 : 2339277 : tetinpoel[i+2],
145 : 2339277 : tetinpoel[i+3]
146 : : }
147 : 2339277 : };
148 : :
149 : : trace_out << "Consume tet " << i << std::endl;
150 [ + - ]: 2339277 : tet_store.add(t, AMR::Refinement_Case::initial_grid);
151 : : }
152 : 2128 : }
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 : 287 : void mesh_adapter_t::mark_uniform_refinement()
179 : : {
180 [ + + ]: 399906 : for (auto& kv : tet_store.edge_store.edges) {
181 : 399619 : auto& local = kv.second;
182 [ + + ]: 399619 : if (local.lock_case == Edge_Lock_Case::unlocked)
183 : 399208 : local.needs_refining = 1;
184 : : }
185 : 287 : mark_refinement();
186 : 287 : }
187 : :
188 : : /**
189 : : * @brief Helper function to apply uniform derefinement to all tets
190 : : */
191 : 72 : void mesh_adapter_t::mark_uniform_derefinement()
192 : : {
193 : 72 : const auto& inp = tet_store.get_active_inpoel();
194 : 72 : auto& edge_store = tet_store.edge_store;
195 [ + + ]: 324860 : for (std::size_t t=0; t<inp.size()/4; ++t) {
196 : : const auto edges =
197 : : edge_store.generate_keys(
198 [ + - ]: 324788 : {inp[t*4+0], inp[t*4+1], inp[t*4+2], inp[t*4+3]});
199 [ + + ]: 2273516 : for (const auto& tetedge : edges) {
200 [ + - ]: 1948728 : auto e = edge_store.edges.find(tetedge);
201 [ + - ]: 1948728 : if (e != end(edge_store.edges)) {
202 : 1948728 : auto& local = e->second;
203 : : //if (local.lock_case == Edge_Lock_Case::unlocked) {
204 : 1948728 : local.needs_derefining = 1;
205 : : // trace_out << "edge marked for deref: " << local.A << " - "
206 : : // << local.B << std::endl;
207 : : //}
208 : : }
209 : : }
210 : : }
211 : 72 : mark_derefinement();
212 : 72 : }
213 : :
214 : : /**
215 : : * @brief For a given set of edges, set their refinement criteria for
216 : : * refinement
217 : : *
218 : : * @param remote Vector of edges and edge tags
219 : : */
220 : 24 : void mesh_adapter_t::mark_error_refinement(
221 : : const std::vector< std::pair< edge_t, edge_tag > >& remote )
222 : : {
223 [ + + ]: 11856 : for (const auto& r : remote) {
224 [ + - ]: 11832 : auto& local = tet_store.edge_store.get( r.first );
225 [ + + ]: 11832 : if (r.second == edge_tag::REFINE) {
226 [ + + ]: 8844 : if (local.lock_case > Edge_Lock_Case::unlocked) {
227 : 254 : local.needs_refining = 0;
228 : : } else {
229 : 8590 : local.needs_refining = 1;
230 : : // an edge deemed to be 'needing refinement', cannot be derefined
231 : 8590 : local.needs_derefining = 0;
232 : : }
233 [ + - ]: 2988 : } else if (r.second == edge_tag::DEREFINE) {
234 [ + + ]: 2988 : if (local.lock_case > Edge_Lock_Case::unlocked) {
235 : 354 : local.needs_derefining = 0;
236 : : } else {
237 : 2634 : local.needs_derefining = 1;
238 : : }
239 : : }
240 : : }
241 : :
242 : 24 : mark_refinement();
243 : 24 : mark_derefinement();
244 : 24 : }
245 : :
246 : 3 : void mesh_adapter_t::mark_error_refinement_corr( const EdgeData& edges )
247 : : {
248 [ + + ]: 29 : for (const auto& r : edges)
249 : : {
250 [ + - ]: 26 : auto& edgeref = tet_store.edge_store.get( edge_t(r.first) );
251 : 26 : edgeref.needs_refining = r.second.first;
252 [ - + ]: 26 : assert(edgeref.lock_case <= r.second.second);
253 : 26 : edgeref.lock_case = r.second.second;
254 : : }
255 : 3 : mark_refinement();
256 : 3 : }
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 : 276621 : 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 : 276621 : int compatibility = 0;
278 : 276621 : 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 [ + + ]: 276621 : if (
336 [ + + ]: 276441 : (refinement_case == AMR::Refinement_Case::one_to_two) or
337 : : (refinement_case == AMR::Refinement_Case::one_to_four)
338 : : )
339 : : {
340 [ + - ]: 576 : if (!normal) {
341 : : trace_out << " not normal 3 " << std::endl;
342 : 576 : compatibility = 3;
343 : : }
344 : : else { // Attempt to allow for "normal" 1:4 and 1:8
345 : 0 : compatibility = 2;
346 : : trace_out << " normal 3 " << std::endl;
347 : : }
348 : :
349 : : }
350 : : else {
351 [ + + ]: 276045 : if (num_locked_edges == 0) {
352 : : trace_out << " no lock 1 " << std::endl;
353 : 271709 : compatibility = 1;
354 : : }
355 : : else {
356 : : trace_out << " lock 2 " << std::endl;
357 : 4336 : compatibility = 2;
358 : : }
359 : : }
360 : :
361 [ - + ]: 276621 : assert(compatibility > 0);
362 [ - + ]: 276621 : assert(compatibility < 4);
363 : 276621 : 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 : 654 : 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 : 654 : const size_t max_num_rounds = AMR_MAX_ROUNDS;
378 : :
379 : : // Mark refinements
380 : : size_t iter;
381 : : //Iterate until convergence
382 [ + - ]: 1046 : for (iter = 0; iter < max_num_rounds; iter++)
383 : : {
384 : :
385 : 1046 : tet_store.marked_refinements.get_state_changed() = false;
386 : :
387 : : // Loop over Tets.
388 [ + + ]: 668677 : for (const auto& kv : tet_store.tets)
389 : : {
390 : 667631 : 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 [ + - ][ + + ]: 667631 : if (tet_store.is_active(tet_id)) {
396 : 602964 : int num_locked_edges = 0;
397 : 602964 : int num_intermediate_edges = 0;
398 : :
399 : : // Loop over nodes and count the number which need refining
400 : 602964 : int num_to_refine = 0;
401 : :
402 : : // This is useful for later inspection
403 [ + - ]: 602964 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
404 : :
405 : : //Iterate over edges
406 [ + + ]: 4220748 : 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 [ + - ][ - + ]: 3617784 : 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 [ + - ][ + + ]: 3617784 : else if(tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::intermediate)
421 : : {
422 : : trace_out << "Found intermediate edge " << key << std::endl;
423 : 45592 : num_intermediate_edges++;
424 : : }
425 : : else
426 : : {
427 : : // Count edges which need refining
428 : : // We check in here as we won't refine a
429 : : // locked edge and will thus ignore it
430 [ + - ][ + + ]: 3572192 : if (tet_store.edge_store.get(key).needs_refining == 1)
431 : : {
432 : 1629446 : num_to_refine++;
433 : : trace_out << "key needs ref " << key << std::endl;
434 : : }
435 : : }
436 : : }
437 : :
438 : : // TODO: Should this be a reference?
439 [ + - ]: 602964 : AMR::Refinement_Case refinement_case = tet_store.get_refinement_case(tet_id);
440 [ + - ]: 602964 : int normal = tet_store.is_normal(tet_id);
441 : :
442 : : trace_out << "Checking " << tet_id <<
443 : : " ref case " << refinement_case <<
444 : : " num ref " << num_to_refine <<
445 : : " normal " << normal <<
446 : : std::endl;
447 : :
448 : :
449 : :
450 : : //If we have any tets to refine
451 [ + + ]: 602964 : if (num_to_refine > 0)
452 : : {
453 : : //Determine compatibility case
454 : :
455 [ + - ]: 276621 : int compatibility = detect_compatibility(num_locked_edges,
456 : : num_intermediate_edges, refinement_case, normal);
457 : :
458 : : trace_out << "Compat " << compatibility << std::endl;
459 : :
460 : : // Now check num_to_refine against situations
461 [ + + ]: 276621 : if (compatibility == 1)
462 : : {
463 [ + - ]: 271709 : refinement_class_one(num_to_refine, tet_id);
464 : : }
465 [ + + ]: 4912 : else if (compatibility == 2)
466 : : {
467 [ + - ]: 4336 : refinement_class_two(edge_list, tet_id);
468 : : }
469 [ + - ]: 576 : else if (compatibility == 3)
470 : : {
471 [ + - ]: 576 : refinement_class_three(tet_id);
472 : : }
473 : :
474 : : /*
475 : : // Write temp mesh out
476 : : std::string temp_file = "temp." +
477 : : std::to_string(iter) + "." +
478 : : std::to_string(tet_id) + ".exo";
479 : :
480 : : std::cout << "Writing " << temp_file << std::endl;
481 : : Adaptive_UnsMesh outmesh(
482 : : get_active_inpoel(), x(), y(), z()
483 : : );
484 : : tk::ExodusIIMeshWriter( temp_file, tk::ExoWriter::CREATE ).
485 : : writeMesh(outmesh);
486 : : */
487 : :
488 : : } // if num_to_refine
489 : : else {
490 : : // If we got here, we don't want to refine this guy
491 [ + - ]: 326343 : tet_store.marked_refinements.add(tet_id, AMR::Refinement_Case::none);
492 : : }
493 : : } // if active
494 : : else {
495 : : trace_out << "Inactive" << std::endl;
496 : : }
497 : : } // For
498 : :
499 : : // If nothing changed during that round, break
500 [ + + ]: 1046 : if (!tet_store.marked_refinements.get_state_changed())
501 : : {
502 : : trace_out << "Terminating loop at iter " << iter << std::endl;
503 : 654 : break;
504 : : }
505 : : trace_out << "End iter " << iter << std::endl;
506 : : }
507 : : trace_out << "Loop took " << iter << " rounds." << std::endl;
508 : :
509 : : //std::cout << "Print Tets" << std::endl;
510 : : //print_tets();
511 : 654 : }
512 : :
513 : : /**
514 : : * @brief Helper function to print tet information when needed
515 : : */
516 : 0 : void mesh_adapter_t::print_tets() {
517 : 0 : tet_store.print_tets();
518 : 0 : }
519 : :
520 : : /**
521 : : * @brief Function to call refinement after each tet has had it's
522 : : * refinement case marked and calculated
523 : : */
524 : 344 : void mesh_adapter_t::perform_refinement()
525 : : {
526 : : // Track tets which need to be deleted this iteration
527 : 688 : std::set<size_t> round_two;
528 : :
529 : : trace_out << "Perform ref" << std::endl;
530 : :
531 : : // Do refinements
532 [ + + ]: 1077732 : for (const auto& kv : tet_store.tets)
533 : : {
534 : 1077388 : size_t tet_id = kv.first;
535 : :
536 : : trace_out << "Do refine of " << tet_id << std::endl;
537 [ + - ][ + + ]: 1077388 : if (tet_store.has_refinement_decision(tet_id))
538 : : {
539 [ + - ][ + + ]: 135942 : switch(tet_store.marked_refinements.get(tet_id))
[ + + ][ + + ]
[ - ]
540 : : {
541 : 511 : case AMR::Refinement_Case::one_to_two:
542 [ + - ]: 511 : refiner.refine_one_to_two(tet_store,node_connectivity,tet_id);
543 : 511 : break;
544 : 706 : case AMR::Refinement_Case::one_to_four:
545 [ + - ]: 706 : refiner.refine_one_to_four(tet_store,node_connectivity,tet_id);
546 : 706 : break;
547 : 98593 : case AMR::Refinement_Case::one_to_eight:
548 [ + - ]: 98593 : refiner.refine_one_to_eight(tet_store,node_connectivity,tet_id);
549 : 98593 : break;
550 : 30 : case AMR::Refinement_Case::two_to_eight:
551 [ + - ][ + - ]: 30 : round_two.insert( tet_store.get_parent_id(tet_id) );
552 : : //std::cout << "2->8\n";
553 : 30 : break;
554 : 66 : case AMR::Refinement_Case::four_to_eight:
555 [ + - ][ + - ]: 66 : round_two.insert( tet_store.get_parent_id(tet_id));
556 : : //std::cout << "4->8\n";
557 : 66 : break;
558 : 36036 : case AMR::Refinement_Case::initial_grid:
559 : : // Do nothing
560 : : case AMR::Refinement_Case::none:
561 : : // Do nothing
562 : 36036 : break;
563 : : // No need for default as enum is explicitly covered
564 : : }
565 : : // Mark tet as not needing refinement
566 [ + - ]: 135942 : tet_store.marked_refinements.erase(tet_id);
567 : : }
568 : : }
569 : :
570 : : trace_out << "round_two size " << round_two.size() << std::endl;
571 [ + + ]: 381 : for (const auto i : round_two)
572 : : {
573 : : trace_out << "round two i " << i << std::endl;
574 : :
575 : : // Cache children as we're about to change this data
576 [ + - ][ + - ]: 74 : auto former_children = tet_store.data(i).children;
577 : :
578 [ + - ]: 37 : AMR::Refinement_State& element = tet_store.data(i);
579 : :
580 [ + + ]: 37 : if (element.children.size() == 2)
581 : : {
582 : : trace_out << "perform 2:8" << std::endl;
583 [ + - ]: 15 : refiner.derefine_two_to_one(tet_store,node_connectivity,i);
584 : : }
585 [ + - ]: 22 : else if (element.children.size() == 4)
586 : : {
587 : : trace_out << "perform 4:8" << std::endl;
588 [ + - ]: 22 : refiner.derefine_four_to_one(tet_store,node_connectivity,i);
589 : : }
590 : : else {
591 [ - - ][ - - ]: 0 : std::cout << "num children " << element.children.size() << std::endl;
[ - - ]
592 : 0 : assert(0);
593 : : }
594 : : // remove tets and edges marked for deletion above
595 [ + - ]: 37 : refiner.delete_intermediates_of_children(tet_store);
596 [ + - ]: 37 : tet_store.process_delete_list();
597 : :
598 [ + - ]: 37 : refiner.refine_one_to_eight(tet_store,node_connectivity,i);
599 : :
600 : : // Grab children after it has been updated
601 [ + - ][ + - ]: 37 : auto current_children = tet_store.data(i).children;
602 : :
603 : : // I want to set the children stored in *my* own children, to be
604 : : // the value of my new children....
605 : : //refiner.overwrite_children(tet_store, former_children, current_children);
606 : :
607 [ + - ]: 37 : tet_store.unset_marked_children(i); // FIXME: This will not work well in parallel
608 : 37 : element.refinement_case = AMR::Refinement_Case::one_to_eight;
609 : : }
610 : :
611 : : // Clean up dead edges
612 : : // clean_up_dead_edges(); // Nothing get's marked as "dead" atm?
613 : :
614 : : //std::cout << "Total Edges : " << tet_store.edge_store.size() << std::endl;
615 : : //std::cout << "Total Tets : " << tet_store.size() << std::endl;
616 : : //std::cout << "Total Nodes : " << m_x.size() << std::endl;
617 : :
618 : : trace_out << "Done ref" << std::endl;
619 : 344 : node_connectivity.print();
620 : 344 : node_connectivity.print();
621 [ + - ]: 344 : tet_store.print_node_types();
622 [ + - ]: 344 : tet_store.print_tets();
623 : : //node_connectivity.print();
624 : :
625 : : //reset_intermediate_edges();
626 [ + - ]: 344 : remove_edge_locks(1);
627 [ + - ]: 344 : remove_normals();
628 : :
629 [ + - ]: 344 : lock_intermediates();
630 : :
631 [ + + ]: 1546468 : for (auto& kv : tet_store.edge_store.edges) {
632 : 1546124 : auto& local = kv.second;
633 : 1546124 : local.needs_refining = 0;
634 : : }
635 : 344 : }
636 : :
637 : 684 : void mesh_adapter_t::lock_intermediates()
638 : : {
639 : : /*
640 : : for (auto k : tet_store.intermediate_list)
641 : : {
642 : : refiner.lock_edges_from_node(tet_store,k, Edge_Lock_Case::intermediate);
643 : : }
644 : : */
645 : : // TODO: Passing tet_store twice probably isn't the best
646 [ + - ]: 684 : refiner.lock_intermediates(tet_store, tet_store.intermediate_list, Edge_Lock_Case::intermediate);
647 : 684 : }
648 : :
649 : : /**
650 : : * @brief A method implementing "Algorithm 1" from the paper
651 : : *
652 : : * @param num_to_refine Number of edges to refine
653 : : * @param tet_id The id of the given tet
654 : : */
655 : 271709 : void mesh_adapter_t::refinement_class_one(int num_to_refine, size_t tet_id)
656 : : {
657 : : trace_out << "Refinement Class One" << std::endl;
658 : :
659 : : // "If nrefine = 1
660 : : // Accept as a 1:2 refinement"
661 [ + + ]: 271709 : if (num_to_refine == 1)
662 : : {
663 : 1048 : tet_store.mark_one_to_two(tet_id);
664 : : }
665 : :
666 : : // "Else if nrefine = 2 OR nrefine = 3"
667 [ + - ][ + + ]: 270661 : else if (num_to_refine > 1 && num_to_refine < 4)
668 : : {
669 : :
670 : : // We need to detect if the edges which need to refine are
671 : : // on the same face
672 : : // and if so which face so we know how to 1:4
673 : :
674 [ + - ]: 1410 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
675 : 1410 : bool edges_on_same_face = false;
676 : 1410 : size_t face_refine_id = 0;
677 : :
678 : : // Iterate over each face
679 [ + + ]: 4110 : for (size_t face = 0; face < NUM_TET_FACES; face++)
680 : : {
681 : 3999 : int num_face_refine_edges = 0;
682 : 3999 : face_ids_t face_ids = face_list[face];
683 : :
684 : : trace_out << "Face is " <<
685 : : face_ids[0] << ", " <<
686 : : face_ids[1] << ", " <<
687 : : face_ids[2] << ", " <<
688 : : std::endl;
689 : :
690 [ + - ]: 3999 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
691 : : // For this face list, see which ones need refining
692 [ + + ]: 15996 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
693 : : {
694 : 11997 : edge_t key = face_edge_list[k];
695 [ + - ][ + + ]: 11997 : if (tet_store.edge_store.get(key).needs_refining == 1)
696 : : {
697 : 6723 : num_face_refine_edges++;
698 : : }
699 : : }
700 [ + + ]: 3999 : if (num_face_refine_edges == num_to_refine)
701 : : {
702 : 1299 : edges_on_same_face = true;
703 : 1299 : face_refine_id = face;
704 : : trace_out << "Breaking with face value " << face << std::endl;
705 : 1299 : break;
706 : : }
707 : : }
708 : :
709 : : // "If active edges are on the same face
710 : : // Activate any inactive edges of the face
711 : : // Accept as a 1:4 // refinement"
712 [ + + ]: 1410 : if (edges_on_same_face)
713 : : {
714 : 1299 : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list,
715 : : face_refine_id);
716 : :
717 [ + - ]: 1299 : tet_t tet = tet_store.get(tet_id);
718 : 1299 : size_t opposite_id = tet[opposite_offset];
719 : :
720 : : trace_out << "face_refine_id " << face_refine_id << std::endl;
721 : : trace_out << "opposite_offset " << opposite_offset << std::endl;
722 : : trace_out << "opposite_id " << opposite_id << std::endl;
723 : :
724 : : // Activate edges on this face
725 [ + - ]: 1299 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_list[face_refine_id]);
726 : :
727 [ + + ]: 5196 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
728 : : {
729 : 3897 : edge_t key = face_edge_list[k];
730 [ + - ]: 3897 : tet_store.edge_store.mark_for_refinement(key);
731 : : }
732 : :
733 : : //refiner.refine_one_to_four(tet_id, face_list[face_refine_id],
734 : : //opposite_id);
735 [ + - ]: 1299 : tet_store.mark_one_to_four(tet_id);
736 : : }
737 : : // "Else if active edges are not on the same face
738 : : // Activate all edges
739 : : // Accept as a 1:8 refinement"
740 : : else {
741 : : //refiner.refine_one_to_eight(tet_id);
742 [ + - ]: 111 : tet_store.mark_edges_for_refinement(tet_id);
743 [ + - ]: 111 : tet_store.mark_one_to_eight(tet_id);
744 : 1410 : }
745 : :
746 : : }
747 : :
748 : : // "Else if nrefine > 3
749 : : // Activate any inactive edges
750 : : // Accept as a 1:8 refinement"
751 [ + - ]: 269251 : else if (num_to_refine > 3)
752 : : {
753 : : //refiner.refine_one_to_eight(tet_id);
754 : 269251 : tet_store.mark_edges_for_refinement(tet_id);
755 : 269251 : tet_store.mark_one_to_eight(tet_id);
756 : : }
757 : 271709 : }
758 : :
759 : : // TODO: Document this
760 : 0 : void mesh_adapter_t::lock_tet_edges(size_t tet_id) {
761 [ - - ]: 0 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
762 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
763 : : {
764 : 0 : edge_t key = edge_list[k];
765 [ - - ][ - - ]: 0 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
766 : : {
767 : : trace_out << "LOCKING! " << key << std::endl;
768 [ - - ]: 0 : tet_store.edge_store.get(key).lock_case = AMR::Edge_Lock_Case::locked;
769 : : }
770 : : }
771 : 0 : }
772 : :
773 : : // TODO: Document this
774 : : // TODO: This has too similar a name to deactivate_tet
775 : 0 : void mesh_adapter_t::deactivate_tet_edges(size_t tet_id) {
776 [ - - ]: 0 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
777 [ - - ]: 0 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
778 : : {
779 : 0 : edge_t key = edge_list[k];
780 : :
781 [ - - ]: 0 : tet_store.edge_store.unmark_for_refinement(key);
782 : : trace_out << "Deactivating " << key << std::endl;
783 [ - - ]: 0 : tet_store.edge_store.get(key).needs_derefining = false;
784 : : }
785 : 0 : }
786 : :
787 : : /**
788 : : * @brief Unmarks edges of given tet for derefinement only
789 : : *
790 : : * @param tet_id The id of the given tet
791 : : */
792 : 103675 : void mesh_adapter_t::deactivate_deref_tet_edges(size_t tet_id) {
793 [ + - ]: 103675 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
794 [ + + ]: 725725 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
795 : : {
796 : 622050 : edge_t key = edge_list[k];
797 : :
798 : : trace_out << "Deactivating " << key << std::endl;
799 [ + - ]: 622050 : tet_store.edge_store.get(key).needs_derefining = false;
800 : : }
801 : 103675 : }
802 : :
803 : : /**
804 : : * @brief An implementation of "Algorithm 2" from the paper
805 : : *
806 : : * @param edge_list The list of edges for the given tet
807 : : * @param tet_id The id of the given tet
808 : : */
809 : 4336 : void mesh_adapter_t::refinement_class_two(edge_list_t edge_list, size_t tet_id)
810 : : {
811 : : trace_out << "Refinement Class Two" << std::endl;
812 : :
813 : :
814 : : // "Deactivate all locked edges"
815 : :
816 : : // count number of active edges
817 : 4336 : int num_active_edges = 0;
818 [ + + ]: 30352 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
819 : : {
820 : 26016 : edge_t key = edge_list[k];
821 [ + - ][ + + ]: 26016 : if (tet_store.edge_store.get(key).lock_case != AMR::Edge_Lock_Case::unlocked)
822 : : {
823 [ + - ]: 16056 : tet_store.edge_store.unmark_for_refinement(key);
824 : : }
825 : : // "Count number of active edges"
826 [ + - ][ + + ]: 26016 : if (tet_store.edge_store.get(key).needs_refining == 1) {
827 : 9101 : num_active_edges++;
828 : : }
829 : : }
830 : :
831 : : // Find out of two active edges live on the same face
832 : 4336 : bool face_refine = false;
833 : 4336 : size_t face_refine_id = 0; // FIXME: Does this need a better default
834 [ + - ]: 4336 : face_list_t face_list = tet_store.generate_face_lists(tet_id);
835 : :
836 : : // Iterate over each face
837 [ + + ]: 15235 : for (size_t face = 0; face < NUM_TET_FACES; face++)
838 : : {
839 : : trace_out << "face " << face << std::endl;
840 : 13298 : int num_face_refine_edges = 0;
841 : 13298 : int num_face_locked_edges = 0;
842 : :
843 : 13298 : face_ids_t face_ids = face_list[face];
844 [ + - ]: 13298 : edge_list_t face_edge_list = AMR::edge_store_t::generate_keys_from_face_ids(face_ids);
845 : : // For this face list, see which ones need refining
846 [ + + ]: 53192 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
847 : : {
848 : 39894 : edge_t key = face_edge_list[k];
849 : : trace_out << "Checking " << key << std::endl;
850 [ + - ][ + + ]: 39894 : if (tet_store.edge_store.get(key).needs_refining == 1)
851 : : {
852 : 14167 : num_face_refine_edges++;
853 : : trace_out << "ref! " << key << std::endl;
854 : : }
855 : :
856 : : // Check for locked edges
857 : : // This case only cares about faces with no locks
858 [ + - ][ + + ]: 39894 : if (tet_store.edge_store.get(key).lock_case != AMR::Edge_Lock_Case::unlocked)
859 : : {
860 : 24020 : num_face_locked_edges++;
861 : : trace_out << "locked! " << key << std::endl;
862 : : }
863 : : }
864 : :
865 : :
866 : : // Decide if we want to process this face
867 [ + + ][ + - ]: 13298 : if (num_face_refine_edges >= 2 && num_face_locked_edges == 0)
868 : : {
869 : : // We can refine this face
870 : 2399 : face_refine = true;
871 : 2399 : face_refine_id = face;
872 : 2399 : break;
873 : : }
874 : : }
875 : :
876 : : // "If nrefine = 1
877 : : // Accept as 1:2 refinement"
878 : : // TODO: can we hoist this higher
879 [ + + ]: 4336 : if (num_active_edges == 1)
880 : : {
881 : : //node_pair_t nodes = find_single_refinement_nodes(edge_list);
882 : : //refine_one_to_two( tet_id, nodes[0], nodes[1]);
883 [ + - ]: 1937 : tet_store.mark_one_to_two(tet_id);
884 : : }
885 : : // "Else if any face has nrefine >= 2 AND no locked edges
886 : : // Active any inactive edges of the face
887 : : // Accept as a 1:4 refinement"
888 [ + - ]: 2399 : else if (face_refine)
889 : : {
890 : 2399 : size_t opposite_offset = AMR::node_connectivity_t::face_list_opposite(face_list, face_refine_id);
891 : :
892 [ + - ]: 2399 : tet_t tet = tet_store.get(tet_id);
893 : 2399 : size_t opposite_id = tet[opposite_offset];
894 : :
895 : : trace_out << "Tet ID " << tet_id << std::endl;
896 : : trace_out << "Opposite offset " << opposite_offset << std::endl;
897 : : trace_out << "Opposite id " << opposite_id << std::endl;
898 : : trace_out << "Face refine id " << face_refine_id << std::endl;
899 : :
900 : : edge_list_t face_edge_list =
901 [ + - ]: 2399 : AMR::edge_store_t::generate_keys_from_face_ids(face_list[face_refine_id]);
902 : :
903 [ + + ]: 9596 : for (size_t k = 0; k < NUM_FACE_NODES; k++)
904 : : {
905 : 7197 : edge_t key = face_edge_list[k];
906 [ + - ]: 7197 : tet_store.edge_store.mark_for_refinement(key);
907 : : }
908 : :
909 : : //refiner.refine_one_to_four(tet_id, face_list[face_refine_id], opposite_id);
910 [ + - ]: 2399 : tet_store.mark_one_to_four(tet_id);
911 : : }
912 : :
913 : : // "Else
914 : : // Deactivate all edges
915 : : // Mark all edges as locked"
916 : : else {
917 : : trace_out << "Class 2 causes some locking.." << std::endl;
918 [ - - ]: 0 : deactivate_tet_edges(tet_id);
919 [ - - ]: 0 : lock_tet_edges(tet_id);
920 : : }
921 : :
922 : 4336 : }
923 : :
924 : : /**
925 : : * @brief Based on a tet_id, decide if it's current state of locked
926 : : * and marked edges maps to a valid refinement case. The logic for
927 : : * this was derived from talking to JW and reading Chicoma.
928 : : *
929 : : * It basically just checks if something a 1:2 and has 3
930 : : * intermediates and 3 makred edges, or is a 1:4 and has 5/6
931 : : * intermediates
932 : : *
933 : : * @param child_id the id of the tet to check
934 : : *
935 : : * @return A bool saying if the tet is in a valid state to be refined
936 : : */
937 : 1944 : bool mesh_adapter_t::check_valid_refinement_case(size_t child_id) {
938 : :
939 : : trace_out << "check valid ref " << child_id << std::endl;
940 [ + - ]: 1944 : edge_list_t edge_list = tet_store.generate_edge_keys(child_id);
941 : :
942 : 1944 : size_t num_to_refine = 0;
943 : 1944 : size_t num_intermediate = 0;
944 : 1944 : size_t unlocked = 0;
945 : 1944 : size_t locked = 0;
946 : :
947 [ + + ]: 13608 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
948 : : {
949 : 11664 : edge_t key = edge_list[k];
950 : : trace_out << "Key " << key << std::endl;
951 : :
952 : : // Count intermediate edges
953 [ + - ][ + + ]: 11664 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::intermediate)
954 : : {
955 : : trace_out << "found intermediate" << std::endl;
956 : 9396 : num_intermediate++;
957 : : }
958 : :
959 : : // Count number of marked for refinement
960 [ + - ][ + + ]: 11664 : if (tet_store.edge_store.get(key).needs_refining == 1)
961 : : {
962 : : trace_out << "found refine" << std::endl;
963 : 2268 : num_to_refine++;
964 : : }
965 : :
966 [ + - ][ + + ]: 11664 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
967 : : {
968 : : trace_out << "found unlocked" << std::endl;
969 : 2268 : unlocked++;
970 : : }
971 : :
972 [ + - ][ - + ]: 11664 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::locked)
973 : : {
974 : : trace_out << "found locked" << std::endl;
975 : 0 : locked++;
976 : : }
977 : :
978 : : }
979 : :
980 [ + - ]: 1944 : AMR::Refinement_State& element = tet_store.data(child_id);
981 : :
982 : : trace_out <<
983 : : "Intermediates " << num_intermediate <<
984 : : " num to refine " << num_to_refine <<
985 : : " unlocked " << unlocked <<
986 : : " locked " << locked <<
987 : : " Case " << element.refinement_case <<
988 : : std::endl;
989 : :
990 : : // check if element is 1:2
991 [ + + ]: 1944 : if (element.refinement_case == AMR::Refinement_Case::one_to_two)
992 : : {
993 : : // If so check it has 3 intermediates and 3 which need refining
994 [ + - ][ - + ]: 360 : if (num_intermediate != 3 || num_to_refine != 3) {
995 : 0 : return false;
996 : : }
997 : : else {
998 : : trace_out << "True " <<
999 : : "Intermediates " << num_intermediate <<
1000 : : " num to refine " << num_to_refine <<
1001 : : " Case " << element.refinement_case <<
1002 : : " 2:1 " << AMR::Refinement_Case::one_to_two <<
1003 : : std::endl;
1004 : : }
1005 : : }
1006 : :
1007 : : // check if element is 1:4
1008 [ + - ]: 1584 : else if (element.refinement_case == AMR::Refinement_Case::one_to_four)
1009 : : {
1010 : : // TODO: Check if it's a center tet for a 1:4
1011 : : // FIXME: Is this even needed? How else would you get these
1012 : : // combinations? Can't we just combine these two checks?
1013 : :
1014 [ + - ]: 1584 : bool is_center_tet = tet_store.is_center(child_id);
1015 : :
1016 [ + + ]: 1584 : if (is_center_tet)
1017 : : {
1018 [ + - ][ - + ]: 396 : if (num_to_refine != 0 || num_intermediate != 6)
1019 : : {
1020 : : trace_out << "Fail compat 1:4 center" << std::endl;
1021 : 0 : return false;
1022 : : }
1023 : : }
1024 : : else { // Is one of the outsides (not center)
1025 [ + - ][ - + ]: 1188 : if (num_to_refine != 1 || num_intermediate != 5)
1026 : : {
1027 : : trace_out << "Fail compat 1:4 non center" << std::endl;
1028 : 0 : return false;
1029 : : }
1030 : : }
1031 : : }
1032 : :
1033 : : // If it makes it here, it's compatible
1034 : 1944 : return true;
1035 : :
1036 : : }
1037 : :
1038 : : /**
1039 : : * @brief Place holder method for the implementation of "Algorithm
1040 : : * 3" from the paper
1041 : : */
1042 : : // TODO: Does this parse a childs siblings multiple times?
1043 : 576 : void mesh_adapter_t::refinement_class_three(size_t tet_id) {
1044 : :
1045 : : trace_out << "Refinement Class Three" << std::endl;
1046 : :
1047 : : // "Identify parent element iparent"
1048 : : // TODO: WE should either always use the id to fetch, or always do the data lookup
1049 : : //size_t parent_id = master_elements.get_parent(tet_id);
1050 [ + - ]: 576 : size_t parent_id = tet_store.get_parent_id(tet_id);
1051 : :
1052 : : trace_out << "Parent id = " << parent_id << std::endl;
1053 : :
1054 : : // NOTE: This implies comms when we use these ids?
1055 [ + - ][ + - ]: 1152 : child_id_list_t children = tet_store.data(parent_id).children;
1056 : :
1057 : : // "Do for each child element ielement
1058 : : // Activate all non-locked edges
1059 : : // Deactivate all locked edges"
1060 [ + + ]: 2520 : for (size_t i = 0; i < children.size(); i++)
1061 : : {
1062 : : // TODO: Is this in element or tet ids?
1063 : : trace_out << "Checking child " << children[i] << std::endl;
1064 [ + - ]: 1944 : edge_list_t edge_list = tet_store.generate_edge_keys(children[i]);
1065 [ + + ]: 13608 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1066 : : {
1067 : 11664 : edge_t key = edge_list[k];
1068 : : trace_out << "Compat 3 " << key << std::endl;
1069 [ + - ][ + + ]: 11664 : if (tet_store.edge_store.get(key).lock_case == AMR::Edge_Lock_Case::unlocked)
1070 : : {
1071 : : trace_out << "Compat 3 marking edge " << key << std::endl;
1072 [ + - ]: 2268 : tet_store.edge_store.mark_for_refinement(key);
1073 : : }
1074 : : else {
1075 [ + - ]: 9396 : tet_store.edge_store.unmark_for_refinement(key);
1076 : : }
1077 : : }
1078 : : }
1079 : :
1080 : : // "Set compatible = TRUE
1081 : 576 : bool compatible = true;
1082 : :
1083 : : // Do for each child element ielement
1084 : : // If ielement is not a valid refinement case
1085 : : // compatible = FALSE"
1086 [ + + ]: 2520 : for (size_t i = 0; i < children.size(); i++)
1087 : : {
1088 : 1944 : size_t child = children[i];
1089 [ + - ][ - + ]: 1944 : if ( !check_valid_refinement_case(child) )
1090 : : {
1091 : : trace_out << "Compat 3 Marking compatible false because of invalid refinement case" << std::endl;
1092 : :
1093 : 0 : compatible = false;
1094 : : }
1095 : : else {
1096 : : trace_out << "Is compatible" << std::endl;
1097 : : }
1098 : : }
1099 : :
1100 : : // "If compatible = FALSE
1101 : : // Do for each child element ielement
1102 : : // Deactive all edges of ielement
1103 : : // Mark all edges of ielement as locked
1104 : : // Mark ielement as normal"
1105 [ - + ]: 576 : if (compatible == false)
1106 : : {
1107 [ - - ]: 0 : for (size_t i = 0; i < children.size(); i++)
1108 : : {
1109 : 0 : size_t child = children[i];
1110 [ - - ]: 0 : deactivate_tet_edges(child);
1111 [ - - ]: 0 : lock_tet_edges(child);
1112 : : trace_out << "Compat 3 locking edges of " << child << std::endl;
1113 : : // Here we interpret normal to mean "don't treat it like it has intermediates"
1114 [ - - ]: 0 : tet_store.mark_normal(child);
1115 : : trace_out << "Compat 3 " << child << std::endl;
1116 : : }
1117 : : }
1118 : : else {
1119 : : trace_out << "TIME TO 2:8 " << tet_id << std::endl;
1120 : : // Accept as 2:8 or 4:8
1121 [ + - ]: 576 : AMR::Refinement_State& element = tet_store.data(tet_id);
1122 [ + + ]: 576 : if (element.refinement_case == AMR::Refinement_Case::one_to_two)
1123 : : {
1124 [ + - ]: 180 : tet_store.mark_two_to_eight(tet_id);
1125 : : }
1126 [ + - ]: 396 : else if (element.refinement_case == AMR::Refinement_Case::one_to_four)
1127 : : {
1128 [ + - ]: 396 : tet_store.mark_four_to_eight(tet_id);
1129 : : }
1130 : : else {
1131 : : 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;
1132 : : }
1133 : :
1134 : : }
1135 : 576 : }
1136 : :
1137 : 344 : void mesh_adapter_t::remove_normals()
1138 : : {
1139 [ + + ]: 1077910 : for (const auto& kv : tet_store.tets)
1140 : : {
1141 : 1077566 : size_t tet_id = kv.first;
1142 [ + - ]: 1077566 : tet_store.set_normal(tet_id, 0);
1143 : : }
1144 : 344 : }
1145 : :
1146 : 344 : void mesh_adapter_t::remove_edge_locks(int intermediate)
1147 : : {
1148 [ + + ]: 1077910 : for (const auto& kv : tet_store.tets)
1149 : : {
1150 : 1077566 : size_t tet_id = kv.first;
1151 : :
1152 : : trace_out << "Process tet removelock " << tet_id << std::endl;
1153 : :
1154 : : // Only apply checks to tets on the active list
1155 [ + - ][ + + ]: 1077566 : if (tet_store.is_active(tet_id)) {
1156 : : // change it from intermediate to locked
1157 [ + - ]: 951372 : update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::locked, AMR::Edge_Lock_Case::unlocked);
1158 [ + - ]: 951372 : if (intermediate) {
1159 [ + - ]: 951372 : update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::intermediate, AMR::Edge_Lock_Case::unlocked);
1160 : : }
1161 : : }
1162 : : }
1163 : :
1164 : 344 : }
1165 : : //void mesh_adapter_t::reset_intermediate_edges()
1166 : : //{
1167 : : // for (const auto& kv : tet_store.tets)
1168 : : // {
1169 : : // size_t tet_id = kv.first;
1170 : :
1171 : : // trace_out << "Process tet reset " << tet_id << std::endl;
1172 : :
1173 : : // // Only apply checks to tets on the active list
1174 : : // if (tet_store.is_active(tet_id)) {
1175 : : // // change it from intermediate to locked
1176 : : // update_tet_edges_lock_type(tet_id, AMR::Edge_Lock_Case::intermediate, AMR::Edge_Lock_Case::locked);
1177 : : // }
1178 : : // }
1179 : : //}
1180 : :
1181 : 1902744 : void mesh_adapter_t::update_tet_edges_lock_type(size_t tet_id, AMR::Edge_Lock_Case check, AMR::Edge_Lock_Case new_case) {
1182 [ + - ]: 1902744 : edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
1183 [ + + ]: 13319208 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1184 : : {
1185 : 11416464 : edge_t key = edge_list[k];
1186 [ + - ][ + + ]: 11416464 : if (tet_store.edge_store.get(key).lock_case == check)
1187 : : {
1188 [ + - ]: 4274 : tet_store.edge_store.get(key).lock_case = new_case;
1189 : : }
1190 : : }
1191 : 1902744 : }
1192 : :
1193 : 96 : void mesh_adapter_t::mark_derefinement()
1194 : : {
1195 : 96 : const size_t max_num_rounds = AMR_MAX_ROUNDS;
1196 : :
1197 : : // Mark refinements
1198 : : size_t iter;
1199 : : //Iterate until convergence
1200 [ + - ]: 173 : for (iter = 0; iter < max_num_rounds; iter++)
1201 : : {
1202 : 173 : tet_store.marked_derefinements.get_state_changed() = false;
1203 : :
1204 : : // set of elements which have been considered for derefinement
1205 : 173 : std::unordered_set< size_t > done_deref_marking;
1206 : :
1207 : : // Loop over tets
1208 [ + + ]: 761464 : for (const auto& kv : tet_store.tets)
1209 : : {
1210 : : // this loop only runs for active tets
1211 [ + - ][ + + ]: 761291 : if (!tet_store.is_active(kv.first)) {
1212 [ + - ]: 88525 : deactivate_deref_tet_edges(kv.first);
1213 : 678023 : continue;
1214 : : }
1215 : 672766 : size_t activetet_id = kv.first;
1216 : :
1217 : : // check if activetet_id has a parent (assign to tet_id)
1218 : : // if it does not, activetet_id is not a derefinement candidate
1219 : : size_t tet_id;
1220 [ + - ]: 672766 : const auto& activetet_data = tet_store.data(activetet_id);
1221 [ + + ]: 672766 : if (!activetet_data.has_parent) {
1222 [ + - ]: 2336 : deactivate_deref_tet_edges(activetet_id);
1223 : 2336 : continue;
1224 : : }
1225 : : else {
1226 : 670430 : tet_id = activetet_data.parent_id;
1227 : : }
1228 : :
1229 : : // if already considered for deref, do not reconsider
1230 [ + - ][ + + ]: 670430 : if (done_deref_marking.count(tet_id) > 0) {
1231 : 585928 : continue;
1232 : : }
1233 [ + - ]: 84502 : done_deref_marking.insert(tet_id);
1234 : :
1235 [ + - ][ + - ]: 84502 : child_id_list_t children = tet_store.data(tet_id).children;
1236 : :
1237 : : // check if any child of tet_id (i.e. any active tet) is marked
1238 : : // for refinement
1239 : 84502 : bool is_child_ref(false);
1240 [ + + ]: 755406 : for (size_t i=0; i<children.size(); i++) {
1241 [ + - ]: 670904 : edge_list_t chedge_list = tet_store.generate_edge_keys(children[i]);
1242 : : // Check each edge, see if it is marked for refinement
1243 [ + + ]: 4696328 : for (size_t k=0; k<NUM_TET_EDGES; k++) {
1244 : 4025424 : edge_t edge = chedge_list[k];
1245 : :
1246 [ + - ][ + + ]: 4025424 : if (tet_store.edge_store.get(edge).needs_refining == 1) {
1247 : 47834 : is_child_ref = true;
1248 : 47834 : continue;
1249 : : }
1250 : : }
1251 : : }
1252 : : // deactivate from deref if marked for ref
1253 [ + + ]: 84502 : if (is_child_ref) {
1254 [ + + ]: 11106 : for (auto child_id : children) {
1255 [ + - ]: 9872 : deactivate_deref_tet_edges(child_id);
1256 : : }
1257 [ + - ]: 1234 : deactivate_deref_tet_edges(tet_id);
1258 : 1234 : continue;
1259 : : }
1260 : :
1261 : : // This is useful for later inspection
1262 : : //edge_list_t edge_list = tet_store.generate_edge_keys(tet_id);
1263 : 83268 : std::size_t num_to_derefine = 0; // Nodes
1264 : :
1265 [ + - ]: 83268 : AMR::Refinement_Case refinement_case = tet_store.get_refinement_case(tet_id);
1266 : :
1267 [ + - ]: 166536 : auto derefine_node_set = refiner.find_derefine_node_set(tet_store, tet_id);
1268 : : // Find the set of nodes which are not in the parent
1269 : : std::unordered_set<size_t> non_parent_nodes =
1270 [ + - ]: 166536 : refiner.child_exclusive_nodes(tet_store, tet_id);
1271 : :
1272 : : //for (auto drnode: derefine_node_set)
1273 : : // trace_out << "derefine node: " << drnode << std::endl;
1274 : :
1275 : 83268 : num_to_derefine = derefine_node_set.size();
1276 : :
1277 : : if (num_to_derefine > 0) {
1278 : : trace_out << "num_to_derefine " << num_to_derefine << std::endl;
1279 : : trace_out << "ref_case " << refinement_case << std::endl;
1280 : : trace_out << "num children " << children.size() << std::endl;
1281 : : }
1282 : :
1283 : : //num_to_derefine = convert_derefine_edges_to_points(tet_store, tet_id, num_edges_to_derefine, refinement_case);
1284 : :
1285 : : // "If nderefine = 1
1286 [ + + ]: 83268 : if (num_to_derefine == 1)
1287 : : {
1288 : : // If icase = 1:2
1289 : :
1290 : : //if (refinement_case == AMR::Refinement_Case::one_to_two)
1291 [ + + ]: 393 : if (children.size() == 2)
1292 : : {
1293 : : // Accept as 2:1 derefine"
1294 : : trace_out << "Accept as 2:1" << std::endl;
1295 : : //refiner.derefine_two_to_one(tet_store, node_connectivity, tet_id);
1296 [ + - ]: 378 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::two_to_one);
1297 : : }
1298 : : // "Else
1299 : : else {
1300 : : // Deactivate all points"
1301 [ + + ]: 135 : for (auto child_id : children) {
1302 [ + - ]: 120 : deactivate_deref_tet_edges(child_id);
1303 : : }
1304 [ + - ]: 15 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1305 : : trace_out << "giving up on deref decision. deactivate near 2:1 ntd = 1" << std::endl;
1306 : : }
1307 : : }
1308 : :
1309 : :
1310 : : // "If nderefine = 2
1311 [ + + ]: 82875 : else if (num_to_derefine == 2)
1312 : : {
1313 : : // If icase = 1:4
1314 : : //if (refinement_case == AMR::Refinement_Case::one_to_four)
1315 [ - + ]: 14 : if (children.size() == 4)
1316 : : {
1317 : : // Accept as 4:2 derefine"
1318 : : trace_out << "Accept as 4:2" << std::endl;
1319 : : //refiner.derefine_four_to_two(tet_store, node_connectivity, tet_id);
1320 [ - - ]: 0 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::four_to_two);
1321 : : }
1322 : : // "Else
1323 : : else {
1324 : : // Deactivate all points"
1325 [ + + ]: 126 : for (auto child_id : children) {
1326 [ + - ]: 112 : deactivate_deref_tet_edges(child_id);
1327 : : }
1328 [ + - ]: 14 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1329 : : trace_out << "giving up on deref decision. deactivate near 4:2 ntd = 2" << std::endl;
1330 : : }
1331 : : }
1332 : :
1333 : : // "If nderefine = 3
1334 [ + + ]: 82861 : else if (num_to_derefine == 3)
1335 : : {
1336 : : // If icase = 1:4
1337 : : //if (refinement_case == AMR::Refinement_Case::one_to_four)
1338 [ + + ]: 515 : if (children.size() == 4)
1339 : : {
1340 : : // Accept as 4:1 derefine"
1341 : : trace_out << "Accept as 4:1" << std::endl;
1342 : : //refiner.derefine_four_to_one(tet_store, node_connectivity, tet_id);
1343 [ + - ]: 502 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::four_to_one);
1344 : : }
1345 : : // "Else if icase = 1:8
1346 : : //else if (refinement_case == AMR::Refinement_Case::one_to_eight)
1347 [ + - ]: 13 : else if (children.size() == 8)
1348 : : {
1349 : : // we have a list of (non-parent) nodes that is marked
1350 : : // for derefinement. First, determine the nodes that are
1351 : : // unmarked for derefinement (or inactive_nodes). Then,
1352 : : // determine if these are on a single face.
1353 : 26 : std::unordered_set<size_t> inactive_node_set;
1354 [ + + ]: 91 : for (auto npn : non_parent_nodes) {
1355 [ + - ][ + + ]: 78 : if (derefine_node_set.count(npn) == 0)
1356 [ + - ]: 39 : inactive_node_set.insert(npn);
1357 : : }
1358 [ - + ][ - - ]: 13 : Assert(inactive_node_set.size() == 3, "Incorrectly "
[ - - ][ - - ]
1359 : : "sized inactive-node set");
1360 [ + - ]: 13 : auto same_face = check_same_face(tet_id, inactive_node_set);
1361 : : // If inactive points lie on same face
1362 [ + + ]: 13 : if (same_face.first == true)
1363 : : {
1364 : : // Accept as 8:4 derefinement
1365 : : trace_out << "Accept as 8:4" << std::endl;
1366 : :
1367 : : // create a vector of node-array-pairs to mark edges
1368 : : // for refinement 1:4
1369 : 16 : std::vector< std::array< std::size_t, 2 > > ref_edges;
1370 [ + + ]: 32 : for (auto n:inactive_node_set) {
1371 [ + - ][ + - ]: 24 : ref_edges.push_back(node_connectivity.get(n));
1372 : : }
1373 : :
1374 [ + - ][ + - ]: 8 : tet_store.edge_store.mark_edges_for_deref_ref(ref_edges);
1375 : : //refiner.derefine_eight_to_four(tet_store, node_connectivity, tet_id);
1376 [ + - ]: 8 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_four);
1377 : : }
1378 : : // "Else
1379 : : else {
1380 : : // Deactivate all points"
1381 [ + + ]: 45 : for (auto child_id : children) {
1382 [ + - ]: 40 : deactivate_deref_tet_edges(child_id);
1383 : : }
1384 [ + - ]: 5 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1385 : : trace_out << "giving up on deref decision. deactivate near 8:4 ntd = 3" << std::endl;
1386 : : }
1387 : :
1388 : : }
1389 : : }
1390 : :
1391 : : // "If nderefine = 4
1392 [ + + ]: 82346 : else if (num_to_derefine == 4)
1393 : : //else if (children.size() == 4)
1394 : : {
1395 : : // we have a list of (non-parent) nodes that is marked
1396 : : // for derefinement. First, determine the nodes that are
1397 : : // unmarked for derefinement (or inactive_nodes). Then,
1398 : : // determine if these are on a single face.
1399 : 8 : std::unordered_set<size_t> inactive_node_set;
1400 [ + + ]: 28 : for (auto npn : non_parent_nodes) {
1401 [ + - ][ + + ]: 24 : if (derefine_node_set.count(npn) == 0)
1402 [ + - ]: 8 : inactive_node_set.insert(npn);
1403 : : }
1404 [ - + ][ - - ]: 4 : Assert(inactive_node_set.size() == 2, "Incorrectly "
[ - - ][ - - ]
1405 : : "sized inactive-node set");
1406 : : // Check if the inactive point belong to the same parent
1407 : : // face and deactivate the third point on that face
1408 [ + - ]: 4 : auto same_face = check_same_face(tet_id, inactive_node_set);
1409 : :
1410 [ + - ]: 4 : if (same_face.first == true)
1411 : : {
1412 : : // deactivate the edges associated with same_face.second
1413 [ + + ]: 36 : for (size_t i = 0; i < children.size(); i++)
1414 : : {
1415 [ + - ]: 32 : edge_list_t edge_list = tet_store.generate_edge_keys(children[i]);
1416 [ + + ]: 224 : for (size_t k = 0; k < NUM_TET_EDGES; k++)
1417 : : {
1418 : 192 : edge_t edge = edge_list[k];
1419 : 192 : size_t A = edge.first();
1420 : 192 : size_t B = edge.second();
1421 [ + + ][ + + ]: 192 : if (A == same_face.second || B == same_face.second)
1422 [ + - ]: 54 : tet_store.edge_store.get(edge).needs_derefining = false;
1423 : : }
1424 : : }
1425 : :
1426 : : // create a vector of node-array-pairs to mark edges
1427 : : // for refinement 1:4
1428 [ + - ]: 4 : inactive_node_set.insert(same_face.second);
1429 : 8 : std::vector< std::array< std::size_t, 2 > > ref_edges;
1430 [ + + ]: 16 : for (auto n:inactive_node_set) {
1431 [ + - ][ + - ]: 12 : ref_edges.push_back(node_connectivity.get(n));
1432 : : }
1433 : :
1434 [ + - ][ + - ]: 4 : tet_store.edge_store.mark_edges_for_deref_ref(ref_edges);
1435 : :
1436 : : // Accept as 8:4 derefinement
1437 : : trace_out << "Accept as 8:4" << std::endl;
1438 : : //refiner.derefine_eight_to_four(tet_store, node_connectivity, tet_id);
1439 [ + - ]: 4 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_four);
1440 : : }
1441 : : // "Else
1442 : : else {
1443 : : // Deactivate all points"
1444 [ - - ]: 0 : for (auto child_id : children) {
1445 [ - - ]: 0 : deactivate_deref_tet_edges(child_id);
1446 : : }
1447 [ - - ]: 0 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1448 : : trace_out << "giving up on deref decision. deactivate near 8:4 ntd = 4" << std::endl;
1449 : : }
1450 : : }
1451 : :
1452 : : // "If nderefine = 5
1453 [ - + ]: 82342 : else if (num_to_derefine == 5)
1454 : : {
1455 : : // Accept as 8:2 derefine"
1456 : : trace_out << "Accept as 8:2 " << std::endl;
1457 : : //refiner.derefine_eight_to_two(tet_store, node_connectivity, tet_id);
1458 [ - - ]: 0 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_two);
1459 : : }
1460 : :
1461 : : // "If nderefine = 6
1462 [ + + ]: 82342 : else if (num_to_derefine == 6)
1463 : : {
1464 : : // Accept as 8:1 derefine"
1465 : : trace_out << "Accept as 8:1" << std::endl;
1466 : : //refiner.derefine_eight_to_one(tet_store, node_connectivity, tet_id);
1467 [ + - ]: 82058 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::eight_to_one);
1468 : : }
1469 : :
1470 : : // "If nderefine = 0
1471 : : else {
1472 [ + - ]: 284 : tet_store.mark_derefinement_decision(tet_id, AMR::Derefinement_Case::skip);
1473 : : // Deactivate all points"
1474 [ + + ]: 1720 : for (auto child_id : children) {
1475 [ + - ]: 1436 : deactivate_deref_tet_edges(child_id);
1476 : : }
1477 : : trace_out << "giving up with no deref decision because nderefine = 0" << std::endl;
1478 : : }
1479 : : }
1480 : :
1481 : : // If nothing changed during that round, break
1482 [ + + ]: 173 : if (!tet_store.marked_derefinements.get_state_changed())
1483 : : {
1484 : : trace_out << "Terminating loop at iter " << iter << std::endl;
1485 : 96 : break;
1486 : : }
1487 : : trace_out << "End iter " << iter << std::endl;
1488 : : // clear out set of elements considered during this iteration
1489 : 77 : done_deref_marking.clear();
1490 : : }
1491 : : trace_out << "Deref Loop took " << iter << " rounds." << std::endl;
1492 : 96 : }
1493 : :
1494 : : // TODO: document
1495 : 289 : void mesh_adapter_t::perform_derefinement()
1496 : : {
1497 : : trace_out << "Perform deref" << std::endl;
1498 : :
1499 : : // Do derefinements
1500 [ + + ]: 576562 : for (const auto& kv : tet_store.tets)
1501 : : {
1502 : 576273 : size_t tet_id = kv.first;
1503 : : //size_t parent_id = 0;
1504 : :
1505 : : // TODO: Do I really want to loop all tets?
1506 : :
1507 : : // TODO: is this doing a double lookup?
1508 [ + - ][ + + ]: 576273 : if (tet_store.has_derefinement_decision(tet_id))
1509 : : {
1510 : : trace_out << "Do derefine of " << tet_id << std::endl;
1511 : : //size_t parent_id = tet_store.get_parent_id(tet_id);
1512 : : //trace_out << "Parent = " << parent_id << std::endl;
1513 [ + - ][ + + ]: 40769 : switch(tet_store.marked_derefinements.get(tet_id))
[ - + ][ - - ]
[ + - ]
1514 : : {
1515 : 130 : case AMR::Derefinement_Case::two_to_one:
1516 [ + - ]: 130 : refiner.derefine_two_to_one(tet_store,node_connectivity,tet_id);
1517 : : trace_out << "Completed derefine 2:1 of " << tet_id << std::endl;
1518 : 130 : break;
1519 : 176 : case AMR::Derefinement_Case::four_to_one:
1520 [ + - ]: 176 : refiner.derefine_four_to_one(tet_store,node_connectivity,tet_id);
1521 : : trace_out << "Completed derefine 4:1 of " << tet_id << std::endl;
1522 : 176 : break;
1523 : 0 : case AMR::Derefinement_Case::four_to_two:
1524 [ - - ]: 0 : refiner.derefine_four_to_two(tet_store,node_connectivity,tet_id);
1525 : : trace_out << "Completed derefine 4:2 of " << tet_id << std::endl;
1526 : 0 : break;
1527 : 40417 : case AMR::Derefinement_Case::eight_to_one:
1528 [ + - ]: 40417 : refiner.derefine_eight_to_one(tet_store,node_connectivity,tet_id);
1529 : : trace_out << "Completed derefine 8:1 of " << tet_id << std::endl;
1530 : 40417 : break;
1531 : 0 : case AMR::Derefinement_Case::eight_to_two:
1532 [ - - ]: 0 : refiner.derefine_eight_to_two(tet_store,node_connectivity,tet_id);
1533 : : trace_out << "Completed derefine 8:2 of " << tet_id << std::endl;
1534 : 0 : break;
1535 : 0 : case AMR::Derefinement_Case::eight_to_four:
1536 [ - - ]: 0 : refiner.derefine_eight_to_four(tet_store,node_connectivity,tet_id);
1537 : : trace_out << "Completed derefine 8:4 of " << tet_id << std::endl;
1538 : 0 : break;
1539 : 46 : case AMR::Derefinement_Case::skip:
1540 : : // What do we do with skip?
1541 : 46 : break;
1542 : : }
1543 : : // Mark tet as not needing derefinement
1544 [ + - ]: 40769 : tet_store.marked_derefinements.erase(tet_id);
1545 : : }
1546 : : }
1547 : :
1548 : 289 : node_connectivity.print();
1549 : 289 : refiner.delete_intermediates_of_children(tet_store);
1550 : 289 : tet_store.process_delete_list();
1551 : 289 : tet_store.print_node_types();
1552 : :
1553 [ + + ]: 1507277 : for (auto& kv : tet_store.edge_store.edges) {
1554 : 1506988 : auto& local = kv.second;
1555 : 1506988 : local.needs_derefining = 0;
1556 : : }
1557 : 289 : }
1558 : :
1559 : : }
1560 : :
1561 : : #if defined(__clang__)
1562 : : #pragma clang diagnostic pop
1563 : : #endif
|
