Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/libmesh_config.h"
22 :
23 : // only compile these functions if the user requests AMR support
24 : #ifdef LIBMESH_ENABLE_AMR
25 :
26 : #include "libmesh/elem.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/mesh_refinement.h"
29 : #include "libmesh/parallel.h"
30 : #include "libmesh/remote_elem.h"
31 :
32 : namespace libMesh
33 : {
34 :
35 :
36 :
37 : //-----------------------------------------------------------------
38 : // Mesh refinement methods
39 0 : bool MeshRefinement::limit_level_mismatch_at_node (const unsigned int max_mismatch)
40 : {
41 : // This function must be run on all processors at once
42 0 : parallel_object_only();
43 :
44 0 : bool flags_changed = false;
45 :
46 :
47 : // Vector holding the maximum element level that touches a node.
48 0 : std::vector<unsigned char> max_level_at_node (_mesh.n_nodes(), 0);
49 0 : std::vector<unsigned char> max_p_level_at_node (_mesh.n_nodes(), 0);
50 :
51 : // Loop over all the active elements & fill the vector
52 0 : for (auto & elem : _mesh.active_element_ptr_range())
53 : {
54 : const unsigned char elem_level =
55 0 : cast_int<unsigned char>(elem->level() +
56 0 : ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
57 : const unsigned char elem_p_level =
58 0 : cast_int<unsigned char>(elem->p_level() +
59 0 : ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
60 :
61 : // Set the max_level at each node
62 0 : for (const Node & node : elem->node_ref_range())
63 : {
64 0 : const dof_id_type node_number = node.id();
65 :
66 0 : libmesh_assert_less (node_number, max_level_at_node.size());
67 :
68 0 : max_level_at_node[node_number] =
69 0 : std::max (max_level_at_node[node_number], elem_level);
70 0 : max_p_level_at_node[node_number] =
71 0 : std::max (max_p_level_at_node[node_number], elem_p_level);
72 : }
73 0 : }
74 :
75 :
76 : // Now loop over the active elements and flag the elements
77 : // who violate the requested level mismatch. Alternatively, if
78 : // _enforce_mismatch_limit_prior_to_refinement is true, swap refinement flags
79 : // accordingly.
80 0 : for (auto & elem : _mesh.active_element_ptr_range())
81 : {
82 0 : const unsigned int elem_level = elem->level();
83 0 : const unsigned int elem_p_level = elem->p_level();
84 :
85 : // Skip the element if it is already fully flagged
86 : // unless we are enforcing mismatch prior to refinement and may need to
87 : // remove the refinement flag(s)
88 0 : if (elem->refinement_flag() == Elem::REFINE &&
89 0 : elem->p_refinement_flag() == Elem::REFINE
90 0 : && !_enforce_mismatch_limit_prior_to_refinement)
91 0 : continue;
92 :
93 : // Loop over the nodes, check for possible mismatch
94 0 : for (const Node & node : elem->node_ref_range())
95 : {
96 0 : const dof_id_type node_number = node.id();
97 :
98 : // Flag the element for refinement if it violates
99 : // the requested level mismatch
100 0 : if ((elem_level + max_mismatch) < max_level_at_node[node_number]
101 0 : && elem->refinement_flag() != Elem::REFINE)
102 : {
103 0 : elem->set_refinement_flag (Elem::REFINE);
104 0 : flags_changed = true;
105 : }
106 0 : if ((elem_p_level + max_mismatch) < max_p_level_at_node[node_number]
107 0 : && elem->p_refinement_flag() != Elem::REFINE)
108 : {
109 0 : elem->set_p_refinement_flag (Elem::REFINE);
110 0 : flags_changed = true;
111 : }
112 :
113 : // Possibly enforce limit mismatch prior to refinement
114 0 : flags_changed |= this->enforce_mismatch_limit_prior_to_refinement(elem, POINT, max_mismatch);
115 : }
116 0 : }
117 :
118 : // If flags changed on any processor then they changed globally
119 0 : this->comm().max(flags_changed);
120 :
121 0 : return flags_changed;
122 : }
123 :
124 :
125 :
126 0 : bool MeshRefinement::limit_level_mismatch_at_edge (const unsigned int max_mismatch)
127 : {
128 : // This function must be run on all processors at once
129 0 : parallel_object_only();
130 :
131 0 : bool flags_changed = false;
132 :
133 :
134 : // Maps holding the maximum element level that touches an edge
135 : std::map<std::pair<unsigned int, unsigned int>, unsigned char>
136 0 : max_level_at_edge;
137 : std::map<std::pair<unsigned int, unsigned int>, unsigned char>
138 0 : max_p_level_at_edge;
139 :
140 0 : std::unique_ptr<const Elem> edge, pedge;
141 : // Loop over all the active elements & fill the maps
142 0 : for (auto & elem : _mesh.active_element_ptr_range())
143 : {
144 : const unsigned char elem_level =
145 0 : cast_int<unsigned char>(elem->level() +
146 0 : ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
147 : const unsigned char elem_p_level =
148 0 : cast_int<unsigned char>(elem->p_level() +
149 0 : ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
150 :
151 : // Set the max_level at each edge
152 0 : for (auto n : elem->edge_index_range())
153 : {
154 0 : elem->build_edge_ptr(edge, n);
155 0 : dof_id_type childnode0 = edge->node_id(0);
156 0 : dof_id_type childnode1 = edge->node_id(1);
157 0 : if (childnode1 < childnode0)
158 0 : std::swap(childnode0, childnode1);
159 :
160 0 : for (const Elem * p = elem; p != nullptr; p = p->parent())
161 : {
162 0 : p->build_edge_ptr(pedge, n);
163 0 : dof_id_type node0 = pedge->node_id(0);
164 0 : dof_id_type node1 = pedge->node_id(1);
165 :
166 0 : if (node1 < node0)
167 0 : std::swap(node0, node1);
168 :
169 : // If elem does not share this edge with its ancestor
170 : // p, refinement levels of elements sharing p's edge
171 : // are not restricted by refinement levels of elem.
172 : // Furthermore, elem will not share this edge with any
173 : // of p's ancestors, so we can safely break out of the
174 : // for loop early.
175 0 : if (node0 != childnode0 && node1 != childnode1)
176 0 : break;
177 :
178 0 : childnode0 = node0;
179 0 : childnode1 = node1;
180 :
181 : std::pair<unsigned int, unsigned int> edge_key =
182 0 : std::make_pair(node0, node1);
183 :
184 : // Try emplacing (edge_key, elem_level) into max_level_at_edge map.
185 : // If emplacing fails, it means the edge_key already existed, so
186 : // we need to take the max of the previous and current values.
187 0 : if (auto [it, emplaced] = max_level_at_edge.emplace(edge_key, elem_level);
188 0 : !emplaced)
189 0 : it->second = std::max(it->second, elem_level);
190 :
191 : // Same thing for p_level
192 0 : if (auto [it, emplaced] = max_p_level_at_edge.emplace(edge_key, elem_p_level);
193 0 : !emplaced)
194 0 : it->second = std::max(it->second, elem_p_level);
195 : }
196 : }
197 0 : }
198 :
199 :
200 : // Now loop over the active elements and flag the elements
201 : // who violate the requested level mismatch
202 0 : for (auto & elem : _mesh.active_element_ptr_range())
203 : {
204 0 : const unsigned int elem_level = elem->level();
205 0 : const unsigned int elem_p_level = elem->p_level();
206 :
207 : // Skip the element if it is already fully flagged
208 0 : if (elem->refinement_flag() == Elem::REFINE &&
209 0 : elem->p_refinement_flag() == Elem::REFINE
210 0 : && !_enforce_mismatch_limit_prior_to_refinement)
211 0 : continue;
212 :
213 : // Loop over the nodes, check for possible mismatch
214 0 : for (auto n : elem->edge_index_range())
215 : {
216 0 : elem->build_edge_ptr(edge, n);
217 0 : dof_id_type node0 = edge->node_id(0);
218 0 : dof_id_type node1 = edge->node_id(1);
219 0 : if (node1 < node0)
220 0 : std::swap(node0, node1);
221 :
222 : std::pair<dof_id_type, dof_id_type> edge_key =
223 0 : std::make_pair(node0, node1);
224 :
225 : // Flag the element for refinement if it violates
226 : // the requested level mismatch
227 0 : if ((elem_level + max_mismatch) < max_level_at_edge[edge_key]
228 0 : && elem->refinement_flag() != Elem::REFINE)
229 : {
230 0 : elem->set_refinement_flag (Elem::REFINE);
231 0 : flags_changed = true;
232 : }
233 :
234 0 : if ((elem_p_level + max_mismatch) < max_p_level_at_edge[edge_key]
235 0 : && elem->p_refinement_flag() != Elem::REFINE)
236 : {
237 0 : elem->set_p_refinement_flag (Elem::REFINE);
238 0 : flags_changed = true;
239 : }
240 :
241 : // Possibly enforce limit mismatch prior to refinement
242 0 : flags_changed |= this->enforce_mismatch_limit_prior_to_refinement(elem, EDGE, max_mismatch);
243 : } // loop over edges
244 0 : } // loop over active elements
245 :
246 : // If flags changed on any processor then they changed globally
247 0 : this->comm().max(flags_changed);
248 :
249 0 : return flags_changed;
250 0 : }
251 :
252 :
253 :
254 21937 : bool MeshRefinement::limit_overrefined_boundary(const signed char max_mismatch)
255 : {
256 : // This function must be run on all processors at once
257 714 : parallel_object_only();
258 :
259 21937 : bool flags_changed = false;
260 :
261 : // Loop over all the active elements & look for mismatches to fix.
262 12949206 : for (auto & elem : _mesh.active_element_ptr_range())
263 : {
264 : // If we don't have an interior_parent then there's nothing to
265 : // be mismatched with.
266 13341038 : if ((elem->dim() >= LIBMESH_DIM) ||
267 6575367 : !elem->interior_parent())
268 6762691 : continue;
269 :
270 : const unsigned char elem_level =
271 5704 : cast_int<unsigned char>(elem->level() +
272 2980 : ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
273 : const unsigned char elem_p_level =
274 3492 : cast_int<unsigned char>(elem->p_level() +
275 2980 : ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
276 :
277 : // get all relevant interior elements
278 512 : std::set<Elem *> neighbor_set;
279 2980 : elem->find_interior_neighbors(neighbor_set);
280 :
281 7096 : for (auto & neighbor : neighbor_set)
282 4116 : if (max_mismatch >= 0)
283 : {
284 4372 : if ((elem_level > neighbor->level() + max_mismatch) &&
285 3428 : (neighbor->refinement_flag() != Elem::REFINE))
286 : {
287 108 : neighbor->set_refinement_flag(Elem::REFINE);
288 1320 : flags_changed = true;
289 : }
290 :
291 4404 : if ((elem_p_level > neighbor->p_level() + max_mismatch) &&
292 0 : (neighbor->p_refinement_flag() != Elem::REFINE))
293 : {
294 0 : neighbor->set_p_refinement_flag(Elem::REFINE);
295 0 : flags_changed = true;
296 : }
297 : }
298 20509 : }
299 :
300 : // If flags changed on any processor then they changed globally
301 21937 : this->comm().max(flags_changed);
302 :
303 21937 : return flags_changed;
304 : }
305 :
306 :
307 :
308 21369 : bool MeshRefinement::limit_underrefined_boundary(const signed char max_mismatch)
309 : {
310 : // This function must be run on all processors at once
311 698 : parallel_object_only();
312 :
313 21369 : bool flags_changed = false;
314 :
315 : // Loop over all the active elements & look for mismatches to fix.
316 12907546 : for (auto & elem : _mesh.active_element_ptr_range())
317 : {
318 : // If we don't have an interior_parent then there's nothing to
319 : // be mismatched with.
320 13301982 : if ((elem->dim() >= LIBMESH_DIM) ||
321 6558395 : !elem->interior_parent())
322 6741955 : continue;
323 :
324 : // get all relevant interior elements
325 264 : std::set<const Elem *> neighbor_set;
326 1632 : elem->find_interior_neighbors(neighbor_set);
327 :
328 3974 : for (const Elem * neighbor : neighbor_set)
329 : {
330 : const unsigned char neighbor_level =
331 4532 : cast_int<unsigned char>(neighbor->level() +
332 2342 : ((neighbor->refinement_flag() == Elem::REFINE) ? 1 : 0));
333 :
334 : const unsigned char neighbor_p_level =
335 2494 : cast_int<unsigned char>(neighbor->p_level() +
336 2342 : ((neighbor->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
337 :
338 2342 : if (max_mismatch >= 0)
339 : {
340 304 : if ((neighbor_level >
341 2478 : elem->level() + max_mismatch) &&
342 1998 : (elem->refinement_flag() != Elem::REFINE))
343 : {
344 8 : elem->set_refinement_flag(Elem::REFINE);
345 284 : flags_changed = true;
346 : }
347 :
348 304 : if ((neighbor_p_level >
349 2494 : elem->p_level() + max_mismatch) &&
350 0 : (elem->p_refinement_flag() != Elem::REFINE))
351 : {
352 0 : elem->set_p_refinement_flag(Elem::REFINE);
353 0 : flags_changed = true;
354 : }
355 : }
356 : } // loop over interior neighbors
357 19973 : }
358 :
359 : // If flags changed on any processor then they changed globally
360 21369 : this->comm().max(flags_changed);
361 :
362 21369 : return flags_changed;
363 : }
364 :
365 :
366 :
367 67418 : bool MeshRefinement::eliminate_unrefined_patches ()
368 : {
369 : // This function must be run on all processors at once
370 2248 : parallel_object_only();
371 :
372 67418 : bool flags_changed = false;
373 :
374 : // Quick return: if unrefined patches are allowed, then we are not
375 : // going to do anything here and can simply return that the flags
376 : // haven't changed.
377 67418 : if (_allow_unrefined_patches)
378 0 : return flags_changed;
379 :
380 : // Note: we *cannot* use a reference to the real pointer here, since
381 : // the pointer may be reseated below and we don't want to reseat
382 : // pointers held by the Mesh.
383 37335184 : for (Elem * elem : _mesh.active_element_ptr_range())
384 : {
385 : // First, see if there's any possibility we might have to flag
386 : // this element for h and p refinement - do we have any visible
387 : // neighbors? Next we'll check to see if any of those neighbors
388 : // are as coarse or coarser than us.
389 1005502 : bool h_flag_me = false,
390 1005502 : p_flag_me = false;
391 19673298 : for (auto neighbor : elem->neighbor_ptr_range())
392 : {
393 : // Quit if the element is not a local boundary
394 19670948 : if (neighbor != nullptr && neighbor != remote_elem)
395 : {
396 1005338 : h_flag_me = true;
397 1005338 : p_flag_me = true;
398 1005338 : break;
399 : }
400 : }
401 :
402 : // Skip the element if it is already fully flagged for refinement
403 19104049 : if (elem->p_refinement_flag() == Elem::REFINE)
404 1346 : p_flag_me = false;
405 20109555 : if (elem->refinement_flag() == Elem::REFINE)
406 : {
407 16722 : h_flag_me = false;
408 302209 : if (!p_flag_me)
409 2704 : continue;
410 : }
411 : // Test the parent if that is already flagged for coarsening
412 18801840 : else if (elem->refinement_flag() == Elem::COARSEN)
413 : {
414 46232 : libmesh_assert(elem->parent());
415 92464 : elem = elem->parent();
416 : // FIXME - this doesn't seem right - RHS
417 1186362 : if (elem->refinement_flag() != Elem::COARSEN_INACTIVE)
418 324 : continue;
419 46232 : p_flag_me = false;
420 : }
421 :
422 19100861 : const unsigned int my_level = elem->level();
423 1005354 : int my_p_adjustment = 0;
424 20106207 : if (elem->p_refinement_flag() == Elem::REFINE)
425 1260 : my_p_adjustment = 1;
426 19082723 : else if (elem->p_refinement_flag() == Elem::COARSEN)
427 : {
428 14 : libmesh_assert_greater (elem->p_level(), 0);
429 14 : my_p_adjustment = -1;
430 : }
431 : const unsigned int my_new_p_level =
432 19100861 : cast_int<unsigned int>(int(elem->p_level()) +
433 : my_p_adjustment);
434 :
435 : // Check all the element neighbors
436 19704194 : for (auto neighbor : elem->neighbor_ptr_range())
437 : {
438 : // Quit if the element is on a local boundary
439 19701205 : if (neighbor == nullptr || neighbor == remote_elem)
440 : {
441 29084 : h_flag_me = false;
442 29084 : p_flag_me = false;
443 29084 : break;
444 : }
445 : // if the neighbor will be equally or less refined than
446 : // we are, then we will not become an unrefined island.
447 : // So if we are still considering h refinement:
448 38083160 : if (h_flag_me &&
449 : // If our neighbor is already at a lower level,
450 : // it can't end up at a higher level even if it
451 : // is flagged for refinement once
452 18911519 : ((neighbor->level() < my_level) ||
453 : // If our neighbor is at the same level but isn't
454 : // flagged for refinement, it won't end up at a
455 : // higher level
456 16235779 : ((neighbor->active()) &&
457 2415291 : (neighbor->refinement_flag() != Elem::REFINE)) ||
458 : // If our neighbor is currently more refined but is
459 : // a parent flagged for coarsening, it will end up
460 : // at the same level.
461 71816 : (neighbor->refinement_flag() == Elem::COARSEN_INACTIVE)))
462 : {
463 : // We've proven we won't become an unrefined island,
464 : // so don't h refine to avoid that.
465 960964 : h_flag_me = false;
466 :
467 : // If we've also proven we don't need to p refine,
468 : // we don't need to check more neighbors
469 18311789 : if (!p_flag_me)
470 70040 : break;
471 : }
472 2641700 : if (p_flag_me)
473 : {
474 : // if active neighbors will have a p level
475 : // equal to or lower than ours, then we do not need to p
476 : // refine ourselves.
477 932208 : if (neighbor->active())
478 : {
479 907262 : int p_adjustment = 0;
480 17882759 : if (neighbor->p_refinement_flag() == Elem::REFINE)
481 50 : p_adjustment = 1;
482 16974437 : else if (neighbor->p_refinement_flag() == Elem::COARSEN)
483 : {
484 10 : libmesh_assert_greater (neighbor->p_level(), 0);
485 10 : p_adjustment = -1;
486 : }
487 16975493 : if (my_new_p_level >=
488 16975493 : cast_int<unsigned int>(int(neighbor->p_level())
489 : + p_adjustment))
490 : {
491 907020 : p_flag_me = false;
492 16971559 : if (!h_flag_me)
493 902748 : break;
494 : }
495 : }
496 : // If we have inactive neighbors, we need to
497 : // test all their active descendants which neighbor us
498 453710 : else if (neighbor->ancestor())
499 : {
500 453710 : if (neighbor->min_new_p_level_by_neighbor(elem,
501 24946 : my_new_p_level + 2) <= my_new_p_level)
502 : {
503 24930 : p_flag_me = false;
504 453004 : if (!h_flag_me)
505 3360 : break;
506 : }
507 : }
508 : }
509 : }
510 :
511 19100861 : if (h_flag_me)
512 : {
513 : // Parents that would create islands should no longer
514 : // coarsen
515 2783 : if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
516 : {
517 1420 : for (auto & child : elem->child_ref_range())
518 : {
519 32 : libmesh_assert_equal_to (child.refinement_flag(),
520 : Elem::COARSEN);
521 32 : child.set_refinement_flag(Elem::DO_NOTHING);
522 : }
523 8 : elem->set_refinement_flag(Elem::INACTIVE);
524 : }
525 : else
526 100 : elem->set_refinement_flag(Elem::REFINE);
527 2783 : flags_changed = true;
528 : }
529 19100861 : if (p_flag_me)
530 : {
531 34 : if (elem->p_refinement_flag() == Elem::COARSEN)
532 0 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
533 : else
534 6 : elem->set_p_refinement_flag(Elem::REFINE);
535 34 : flags_changed = true;
536 : }
537 62922 : }
538 :
539 : // If flags changed on any processor then they changed globally
540 67418 : this->comm().max(flags_changed);
541 :
542 67418 : return flags_changed;
543 : }
544 :
545 :
546 :
547 0 : bool MeshRefinement::enforce_mismatch_limit_prior_to_refinement(Elem * elem,
548 : NeighborType nt,
549 : unsigned max_mismatch)
550 : {
551 : // Eventual return value
552 0 : bool flags_changed = false;
553 :
554 : // If we are enforcing the limit prior to refinement then we
555 : // need to remove flags from any elements marked for refinement that
556 : // would cause a mismatch
557 0 : if (_enforce_mismatch_limit_prior_to_refinement
558 0 : && elem->refinement_flag() == Elem::REFINE)
559 : {
560 : // get all the relevant neighbors since we may have to refine
561 : // elements off edges or corners as well
562 0 : std::set<const Elem *> neighbor_set;
563 :
564 0 : if (nt == POINT)
565 0 : elem->find_point_neighbors(neighbor_set);
566 0 : else if (nt == EDGE)
567 0 : elem->find_edge_neighbors(neighbor_set);
568 : else
569 0 : libmesh_error_msg("Unrecognized NeighborType: " << nt);
570 :
571 : // Loop over the neighbors of element e
572 0 : for (const auto & neighbor : neighbor_set)
573 : {
574 0 : if ((elem->level() + 1 - max_mismatch) > neighbor->level())
575 : {
576 0 : elem->set_refinement_flag(Elem::DO_NOTHING);
577 0 : flags_changed = true;
578 : }
579 0 : if ((elem->p_level() + 1 - max_mismatch) > neighbor->p_level())
580 : {
581 0 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
582 0 : flags_changed = true;
583 : }
584 : } // loop over edge/point neighbors
585 : } // if _enforce_mismatch_limit_prior_to_refinement
586 :
587 0 : return flags_changed;
588 : }
589 :
590 : } // namespace libMesh
591 :
592 :
593 : #endif
|
