Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
20 : #include "libmesh/libmesh_config.h"
21 :
22 : // only compile these functions if the user requests AMR support
23 : #ifdef LIBMESH_ENABLE_AMR
24 :
25 : #include "libmesh/boundary_info.h"
26 : #include "libmesh/elem_range.h"
27 : #include "libmesh/error_vector.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/mesh_communication.h"
31 : #include "libmesh/mesh_refinement.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/parallel_ghost_sync.h"
34 : #include "libmesh/partitioner.h"
35 : #include "libmesh/remote_elem.h"
36 : #include "libmesh/sync_refinement_flags.h"
37 : #include "libmesh/int_range.h"
38 :
39 : #ifdef DEBUG
40 : // Some extra validation for DistributedMesh
41 : #include "libmesh/mesh_tools.h"
42 : #endif // DEBUG
43 :
44 : #ifdef LIBMESH_ENABLE_PERIODIC
45 : #include "libmesh/periodic_boundaries.h"
46 : #endif
47 :
48 : // C++ includes
49 : #include <atomic>
50 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
51 : #include <cmath> // for isnan(), when it's defined
52 : #include <limits>
53 :
54 :
55 : // Anonymous namespace for helper functions
56 : // namespace {
57 : //
58 : // using namespace libMesh;
59 : //
60 : // struct SyncCoarsenInactive
61 : // {
62 : // bool operator() (const Elem * elem) const {
63 : // // If we're not an ancestor, there's no chance our coarsening
64 : // // settings need to be changed.
65 : // if (!elem->ancestor())
66 : // return false;
67 : //
68 : // // If we don't have any remote children, we already know enough to
69 : // // determine the correct refinement flag ourselves.
70 : // //
71 : // // If we know we have a child that isn't being coarsened, that
72 : // // also forces a specific flag.
73 : // //
74 : // // Either way there's nothing we need to communicate.
75 : // bool found_remote_child = false;
76 : // for (auto & child : elem->child_ref_range())
77 : // {
78 : // if (child.refinement_flag() != Elem::COARSEN)
79 : // return false;
80 : // if (&child == remote_elem)
81 : // found_remote_child = true;
82 : // }
83 : // return found_remote_child;
84 : // }
85 : // };
86 : // }
87 :
88 :
89 :
90 : namespace libMesh
91 : {
92 :
93 : //-----------------------------------------------------------------
94 : // Mesh refinement methods
95 425047 : MeshRefinement::MeshRefinement (MeshBase & m) :
96 : ParallelObject(m),
97 400571 : _mesh(m),
98 400571 : _use_member_parameters(false),
99 400571 : _coarsen_by_parents(false),
100 400571 : _refine_fraction(0.3),
101 400571 : _coarsen_fraction(0.0),
102 400571 : _max_h_level(libMesh::invalid_uint),
103 400571 : _coarsen_threshold(10),
104 400571 : _nelem_target(0),
105 400571 : _absolute_global_tolerance(0.0),
106 400571 : _face_level_mismatch_limit(1),
107 400571 : _edge_level_mismatch_limit(0),
108 400571 : _node_level_mismatch_limit(0),
109 400571 : _overrefined_boundary_limit(0),
110 400571 : _underrefined_boundary_limit(0),
111 400571 : _allow_unrefined_patches(false),
112 400571 : _enforce_mismatch_limit_prior_to_refinement(false)
113 : #ifdef LIBMESH_ENABLE_PERIODIC
114 437285 : , _periodic_boundaries(nullptr)
115 : #endif
116 : {
117 425047 : }
118 :
119 :
120 :
121 : #ifdef LIBMESH_ENABLE_PERIODIC
122 10 : void MeshRefinement::set_periodic_boundaries_ptr(PeriodicBoundaries * pb_ptr)
123 : {
124 10 : _periodic_boundaries = pb_ptr;
125 10 : }
126 : #endif
127 :
128 :
129 :
130 402038 : MeshRefinement::~MeshRefinement () = default;
131 :
132 :
133 :
134 63051 : void MeshRefinement::clear ()
135 : {
136 2100 : _new_nodes_map.clear();
137 63051 : }
138 :
139 :
140 :
141 47308501 : Node * MeshRefinement::add_node(Elem & parent,
142 : unsigned int child,
143 : unsigned int node,
144 : processor_id_type proc_id)
145 : {
146 8013324 : LOG_ONE_ENTITY_SCOPE("add_node()", "MeshRefinement");
147 :
148 47308501 : unsigned int parent_n = parent.as_parent_node(child, node);
149 :
150 47308501 : if (parent_n != libMesh::invalid_uint)
151 14854118 : return parent.node_ptr(parent_n);
152 :
153 : const std::vector<std::pair<dof_id_type, dof_id_type>>
154 35125491 : bracketing_nodes = parent.bracketing_nodes(child, node);
155 :
156 : // If we're not a parent node, we *must* be bracketed by at least
157 : // one pair of parent nodes
158 1902432 : libmesh_assert(bracketing_nodes.size());
159 :
160 : // Return the node if it already exists.
161 : //
162 : // We'll leave the processor_id untouched in this case - if we're
163 : // repartitioning later or if this is a new unpartitioned node,
164 : // we'll update it then, and if not then we don't want to update it.
165 33223059 : if (const auto new_node_id = _new_nodes_map.find(bracketing_nodes);
166 : new_node_id != DofObject::invalid_id)
167 23729721 : return _mesh.node_ptr(new_node_id);
168 :
169 : // Otherwise we need to add a new node.
170 : //
171 : // Figure out where to add the point:
172 :
173 525099 : Point p; // defaults to 0,0,0
174 :
175 116385164 : for (auto n : parent.node_index_range())
176 : {
177 : // The value from the embedding matrix
178 106891826 : const Real em_val = parent.embedding_matrix(child,node,n);
179 :
180 106891826 : if (em_val != 0.)
181 : {
182 53398110 : p.add_scaled (parent.point(n), em_val);
183 :
184 : // If we'd already found the node we shouldn't be here
185 2738954 : libmesh_assert_not_equal_to (em_val, 1);
186 : }
187 : }
188 :
189 : // Although we're leaving new nodes unpartitioned at first, with a
190 : // DistributedMesh we would need a default id based on the numbering
191 : // scheme for the requested processor_id.
192 9493338 : Node * new_node = _mesh.add_point (p, DofObject::invalid_id, proc_id);
193 :
194 525099 : libmesh_assert(new_node);
195 :
196 : // But then we'll make sure this node is marked as unpartitioned.
197 9493338 : new_node->processor_id() = DofObject::invalid_processor_id;
198 :
199 : // Add the node to the map.
200 9493338 : _new_nodes_map.add_node(*new_node, bracketing_nodes);
201 :
202 : // Return the address of the new node
203 525099 : return new_node;
204 : }
205 :
206 :
207 :
208 0 : Elem * MeshRefinement::add_elem (Elem * elem)
209 : {
210 0 : libmesh_assert(elem);
211 0 : return _mesh.add_elem (elem);
212 : }
213 :
214 6516497 : Elem * MeshRefinement::add_elem (std::unique_ptr<Elem> elem)
215 : {
216 389790 : libmesh_assert(elem);
217 7296077 : return _mesh.add_elem(std::move(elem));
218 : }
219 :
220 :
221 1590 : void MeshRefinement::create_parent_error_vector(const ErrorVector & error_per_cell,
222 : ErrorVector & error_per_parent,
223 : Real & parent_error_min,
224 : Real & parent_error_max)
225 : {
226 : // This function must be run on all processors at once
227 82 : parallel_object_only();
228 :
229 : // Make sure the input error vector is valid
230 : #ifdef DEBUG
231 11590 : for (const auto & val : error_per_cell)
232 : {
233 11508 : libmesh_assert_greater_equal (val, 0);
234 : // isnan() isn't standard C++ yet
235 : #ifdef isnan
236 : libmesh_assert(!isnan(val));
237 : #endif
238 : }
239 :
240 : // Use a reference to std::vector to avoid confusing
241 : // this->comm().verify
242 82 : std::vector<ErrorVectorReal> & epc = error_per_parent;
243 82 : libmesh_assert(this->comm().verify(epc));
244 : #endif // #ifdef DEBUG
245 :
246 : // error values on uncoarsenable elements will be left at -1
247 1590 : error_per_parent.clear();
248 1672 : error_per_parent.resize(error_per_cell.size(), 0.0);
249 :
250 : {
251 : // Find which elements are uncoarsenable
252 78930 : for (auto & elem : _mesh.active_local_element_ptr_range())
253 : {
254 42409 : Elem * parent = elem->parent();
255 :
256 : // Active elements are uncoarsenable
257 46902 : error_per_parent[elem->id()] = -1.0;
258 :
259 : // Grandparents and up are uncoarsenable
260 121584 : while (parent)
261 : {
262 22384 : parent = parent->parent();
263 107426 : if (parent)
264 : {
265 14446 : const dof_id_type parentid = parent->id();
266 7223 : libmesh_assert_less (parentid, error_per_parent.size());
267 78460 : error_per_parent[parentid] = -1.0;
268 : }
269 : }
270 1426 : }
271 :
272 : // Sync between processors.
273 : // Use a reference to std::vector to avoid confusing
274 : // this->comm().min
275 82 : std::vector<ErrorVectorReal> & epp = error_per_parent;
276 1590 : this->comm().min(epp);
277 : }
278 :
279 : // The parent's error is defined as the square root of the
280 : // sum of the children's errors squared, so errors that are
281 : // Hilbert norms remain Hilbert norms.
282 : //
283 : // Because the children may be on different processors, we
284 : // calculate local contributions to the parents' errors squared
285 : // first, then sum across processors and take the square roots
286 : // second.
287 : Threads::parallel_for_each
288 1590 : (_mesh.active_local_element_stored_range(),
289 63579 : [&error_per_cell, &error_per_parent](const Elem * elem)
290 : {
291 8986 : const Elem * parent = elem->parent();
292 :
293 : // Calculate each contribution to parent cells
294 42409 : if (parent)
295 : {
296 7938 : const dof_id_type parentid = parent->id();
297 3969 : libmesh_assert_less (parentid, error_per_parent.size());
298 :
299 : // If the parent has grandchildren we won't be able to
300 : // coarsen it, so forget it. Otherwise, add this child's
301 : // contribution to the sum of the squared child errors
302 40158 : if (error_per_parent[parentid] != -1.0)
303 38608 : error_per_parent[parentid] += (error_per_cell[elem->id()] *
304 3496 : error_per_cell[elem->id()]);
305 : }
306 42409 : });
307 :
308 : // Sum the vector across all processors
309 1590 : this->comm().sum(static_cast<std::vector<ErrorVectorReal> &>(error_per_parent));
310 :
311 : // Calculate the min and max as we loop
312 1590 : parent_error_min = std::numeric_limits<double>::max();
313 1590 : parent_error_max = 0.;
314 :
315 302728 : for (auto i : index_range(error_per_parent))
316 : {
317 : // If this element isn't a coarsenable parent with error, we
318 : // have nothing to do. Just flag it as -1 and move on
319 : // Note that this->comm().sum might have left uncoarsenable
320 : // elements with error_per_parent=-n_proc, so reset it to
321 : // error_per_parent=-1
322 312646 : if (error_per_parent[i] < 0.)
323 : {
324 257098 : error_per_parent[i] = -1.;
325 257098 : continue;
326 : }
327 :
328 : // The error estimator might have already given us an
329 : // estimate on the coarsenable parent elements; if so then
330 : // we want to retain that estimate
331 45804 : if (error_per_cell[i])
332 : {
333 0 : error_per_parent[i] = error_per_cell[i];
334 0 : continue;
335 : }
336 : // if not, then e_parent = sqrt(sum(e_child^2))
337 : else
338 44040 : error_per_parent[i] = std::sqrt(error_per_parent[i]);
339 :
340 44040 : parent_error_min = std::min (parent_error_min,
341 1764 : error_per_parent[i]);
342 50759 : parent_error_max = std::max (parent_error_max,
343 1764 : error_per_parent[i]);
344 : }
345 1590 : }
346 :
347 :
348 :
349 63051 : void MeshRefinement::update_nodes_map ()
350 : {
351 63051 : this->_new_nodes_map.init(_mesh);
352 63051 : }
353 :
354 :
355 :
356 1400 : bool MeshRefinement::test_level_one (bool libmesh_dbg_var(libmesh_assert_pass)) const
357 : {
358 : // This function must be run on all processors at once
359 1400 : parallel_object_only();
360 :
361 : // We may need a PointLocator for topological_neighbor() tests
362 : // later, which we need to make sure gets constructed on all
363 : // processors at once.
364 1400 : std::unique_ptr<PointLocatorBase> point_locator;
365 :
366 : #ifdef LIBMESH_ENABLE_PERIODIC
367 : bool has_periodic_boundaries =
368 1400 : _periodic_boundaries && !_periodic_boundaries->empty();
369 1400 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
370 :
371 1400 : if (has_periodic_boundaries)
372 200 : point_locator = _mesh.sub_point_locator();
373 : #endif
374 :
375 1400 : bool failure = false;
376 :
377 : #ifndef NDEBUG
378 1400 : Elem * failed_elem = nullptr;
379 1400 : Elem * failed_neighbor = nullptr;
380 : #endif // !NDEBUG
381 :
382 281099 : for (auto & elem : _mesh.active_local_element_ptr_range())
383 1355733 : for (auto n : elem->side_index_range())
384 : {
385 : Elem * neighbor =
386 1076034 : topological_neighbor(elem, point_locator.get(), n);
387 :
388 2046699 : if (!neighbor || !neighbor->active() ||
389 970665 : neighbor == remote_elem)
390 105369 : continue;
391 :
392 970665 : if ((neighbor->level() + 1 < elem->level()) ||
393 1941330 : (neighbor->p_level() + 1 < elem->p_level()) ||
394 970665 : (neighbor->p_level() > elem->p_level() + 1))
395 : {
396 0 : failure = true;
397 : #ifndef NDEBUG
398 0 : failed_elem = elem;
399 0 : failed_neighbor = neighbor;
400 : #endif // !NDEBUG
401 0 : break;
402 : }
403 0 : }
404 :
405 : // If any processor failed, we failed globally
406 1400 : this->comm().max(failure);
407 :
408 1400 : if (failure)
409 : {
410 : // We didn't pass the level one test, so libmesh_assert that
411 : // we're allowed not to
412 : #ifndef NDEBUG
413 0 : if (libmesh_assert_pass)
414 : {
415 0 : libMesh::out << "MeshRefinement Level one failure, element: "
416 0 : << *failed_elem
417 0 : << std::endl;
418 0 : libMesh::out << "MeshRefinement Level one failure, neighbor: "
419 0 : << *failed_neighbor
420 0 : << std::endl;
421 : }
422 : #endif // !NDEBUG
423 0 : libmesh_assert(!libmesh_assert_pass);
424 0 : return false;
425 : }
426 1400 : return true;
427 0 : }
428 :
429 :
430 :
431 654 : bool MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass)) const
432 : {
433 : // This function must be run on all processors at once
434 654 : parallel_object_only();
435 :
436 654 : bool found_flag = false;
437 :
438 : #ifndef NDEBUG
439 654 : Elem * failed_elem = nullptr;
440 : #endif
441 :
442 : // Search for local flags
443 133347 : for (auto & elem : _mesh.active_local_element_ptr_range())
444 132693 : if (elem->refinement_flag() == Elem::REFINE ||
445 132693 : elem->refinement_flag() == Elem::COARSEN ||
446 398079 : elem->p_refinement_flag() == Elem::REFINE ||
447 132693 : elem->p_refinement_flag() == Elem::COARSEN)
448 : {
449 0 : found_flag = true;
450 : #ifndef NDEBUG
451 0 : failed_elem = elem;
452 : #endif
453 0 : break;
454 0 : }
455 :
456 : // If we found a flag on any processor, it counts
457 654 : this->comm().max(found_flag);
458 :
459 654 : if (found_flag)
460 : {
461 : #ifndef NDEBUG
462 0 : if (libmesh_assert_pass)
463 : {
464 : libMesh::out <<
465 0 : "MeshRefinement test_unflagged failure, element: " <<
466 0 : *failed_elem << std::endl;
467 : }
468 : #endif
469 : // We didn't pass the "elements are unflagged" test,
470 : // so libmesh_assert that we're allowed not to
471 0 : libmesh_assert(!libmesh_assert_pass);
472 0 : return false;
473 : }
474 654 : return true;
475 : }
476 :
477 :
478 :
479 19406 : bool MeshRefinement::refine_and_coarsen_elements ()
480 : {
481 : // This function must be run on all processors at once
482 682 : parallel_object_only();
483 :
484 : // We can't yet turn a non-level-one mesh into a level-one mesh
485 682 : if (_face_level_mismatch_limit)
486 682 : libmesh_assert(test_level_one(true));
487 :
488 : // Possibly clean up the refinement flags from
489 : // a previous step. While we're at it, see if this method should be
490 : // a no-op.
491 19406 : std::atomic<bool> elements_flagged{false};
492 :
493 : Threads::parallel_for_each
494 19406 : (_mesh.element_stored_range(),
495 14342997 : [&elements_flagged](Elem * elem)
496 : {
497 : // This might be left over from the last step
498 344078 : const Elem::RefinementState flag = elem->refinement_flag();
499 :
500 : // Set refinement flag to INACTIVE if the
501 : // element isn't active
502 344078 : if (!elem->active())
503 : {
504 84278 : elem->set_refinement_flag(Elem::INACTIVE);
505 84278 : elem->set_p_refinement_flag(Elem::INACTIVE);
506 : }
507 6148503 : else if (flag == Elem::JUST_REFINED)
508 0 : elem->set_refinement_flag(Elem::DO_NOTHING);
509 6148503 : else if (!elements_flagged)
510 : {
511 47464 : if (flag == Elem::REFINE || flag == Elem::COARSEN)
512 606 : elements_flagged = true;
513 : else
514 : {
515 : const Elem::RefinementState pflag =
516 2574 : elem->p_refinement_flag();
517 32689 : if (pflag == Elem::REFINE || pflag == Elem::COARSEN)
518 48 : elements_flagged = true;
519 : }
520 : }
521 8179086 : });
522 :
523 : // Did *any* processor find elements flagged for AMR/C?
524 19406 : bool any_elements_flagged = elements_flagged;
525 19406 : _mesh.comm().max(any_elements_flagged);
526 :
527 : // If we have nothing to do, let's not bother verifying that nothing
528 : // is compatible with nothing.
529 19406 : if (!any_elements_flagged)
530 28 : return false;
531 :
532 : // Parallel consistency has to come first, or coarsening
533 : // along processor boundaries might occasionally be falsely
534 : // prevented
535 : #ifdef DEBUG
536 654 : bool flags_were_consistent = this->make_flags_parallel_consistent();
537 :
538 654 : libmesh_assert (flags_were_consistent);
539 : #endif
540 :
541 : // Smooth refinement and coarsening flags
542 18412 : _smooth_flags(true, true);
543 :
544 : // First coarsen the flagged elements.
545 : const bool coarsening_changed_mesh =
546 18412 : this->_coarsen_elements ();
547 :
548 : // First coarsen the flagged elements.
549 : // FIXME: test_level_one now tests consistency across periodic
550 : // boundaries, which requires a point_locator, which just got
551 : // invalidated by _coarsen_elements() and hasn't yet been cleared by
552 : // prepare_for_use().
553 :
554 : // libmesh_assert(this->make_coarsening_compatible());
555 : // libmesh_assert(this->make_refinement_compatible());
556 :
557 : // FIXME: This won't pass unless we add a redundant find_neighbors()
558 : // call or replace find_neighbors() with on-the-fly neighbor updating
559 : // libmesh_assert(!this->eliminate_unrefined_patches());
560 :
561 : // We can't contract the mesh ourselves anymore - a System might
562 : // need to restrict old coefficient vectors first
563 : // _mesh.contract();
564 :
565 : // First coarsen the flagged elements.
566 : // Now refine the flagged elements. This will
567 : // take up some space, maybe more than what was freed.
568 : const bool refining_changed_mesh =
569 18412 : this->_refine_elements();
570 :
571 : // First coarsen the flagged elements.
572 : // Finally, the new mesh needs to be prepared for use
573 18412 : if (coarsening_changed_mesh || refining_changed_mesh)
574 : {
575 : #ifdef DEBUG
576 646 : _mesh.libmesh_assert_valid_parallel_ids();
577 : #endif
578 :
579 18160 : _mesh.prepare_for_use ();
580 :
581 646 : if (_face_level_mismatch_limit)
582 646 : libmesh_assert(test_level_one(true));
583 646 : libmesh_assert(test_unflagged(true));
584 646 : libmesh_assert(this->make_coarsening_compatible());
585 646 : libmesh_assert(this->make_refinement_compatible());
586 : // FIXME: This won't pass unless we add a redundant find_neighbors()
587 : // call or replace find_neighbors() with on-the-fly neighbor updating
588 : // libmesh_assert(!this->eliminate_unrefined_patches());
589 :
590 18160 : return true;
591 : }
592 : else
593 : {
594 8 : if (_face_level_mismatch_limit)
595 8 : libmesh_assert(test_level_one(true));
596 8 : libmesh_assert(test_unflagged(true));
597 8 : libmesh_assert(this->make_coarsening_compatible());
598 8 : libmesh_assert(this->make_refinement_compatible());
599 : }
600 :
601 : // Otherwise there was no change in the mesh,
602 : // let the user know. Also, there is no need
603 : // to prepare the mesh for use since it did not change.
604 8 : return false;
605 :
606 : }
607 :
608 :
609 :
610 :
611 :
612 :
613 :
614 25759 : bool MeshRefinement::coarsen_elements ()
615 : {
616 : // This function must be run on all processors at once
617 856 : parallel_object_only();
618 :
619 : // We can't yet turn a non-level-one mesh into a level-one mesh
620 856 : if (_face_level_mismatch_limit)
621 4 : libmesh_assert(test_level_one(true));
622 :
623 : // Possibly clean up the refinement flags from
624 : // a previous step
625 : Threads::parallel_for_each
626 25759 : (_mesh.element_stored_range(),
627 541574 : [](Elem * elem)
628 : {
629 : // Set refinement flag to INACTIVE if the
630 : // element isn't active
631 1083148 : if (!elem->active())
632 : {
633 184474 : elem->set_refinement_flag(Elem::INACTIVE);
634 184474 : elem->set_p_refinement_flag(Elem::INACTIVE);
635 : }
636 :
637 : // This might be left over from the last step
638 10268273 : if (elem->refinement_flag() == Elem::JUST_REFINED)
639 135776 : elem->set_refinement_flag(Elem::DO_NOTHING);
640 541574 : });
641 :
642 : // Parallel consistency has to come first, or coarsening
643 : // along processor boundaries might occasionally be falsely
644 : // prevented
645 25759 : bool flags_were_consistent = this->make_flags_parallel_consistent();
646 :
647 : // In theory, we should be able to remove the above call, which can
648 : // be expensive and should be unnecessary. In practice, doing
649 : // consistent flagging in parallel is hard, it's impossible to
650 : // verify at the library level if it's being done by user code, and
651 : // we don't want to abort large parallel runs in opt mode... but we
652 : // do want to warn that they should be fixed.
653 856 : libmesh_assert(flags_were_consistent);
654 :
655 856 : if (!flags_were_consistent)
656 : libmesh_warning("Warning: Refinement flags were not consistent between processors! "
657 : "Correcting and continuing.\n");
658 :
659 : // Smooth coarsening flags
660 25759 : _smooth_flags(false, true);
661 :
662 : // Coarsen the flagged elements.
663 : const bool mesh_changed =
664 25759 : this->_coarsen_elements ();
665 :
666 856 : if (_face_level_mismatch_limit)
667 4 : libmesh_assert(test_level_one(true));
668 856 : libmesh_assert(this->make_coarsening_compatible());
669 : // FIXME: This won't pass unless we add a redundant find_neighbors()
670 : // call or replace find_neighbors() with on-the-fly neighbor updating
671 : // libmesh_assert(!this->eliminate_unrefined_patches());
672 :
673 : // We can't contract the mesh ourselves anymore - a System might
674 : // need to restrict old coefficient vectors first
675 : // _mesh.contract();
676 :
677 : // Finally, the new mesh may need to be prepared for use
678 25759 : if (mesh_changed)
679 213 : _mesh.prepare_for_use ();
680 :
681 25759 : return mesh_changed;
682 : }
683 :
684 :
685 :
686 :
687 :
688 :
689 :
690 26745 : bool MeshRefinement::refine_elements ()
691 : {
692 : // This function must be run on all processors at once
693 882 : parallel_object_only();
694 :
695 : // Prevent refinement when the mesh has disjoint neighbor boundary pairs.
696 : //
697 : // Disjoint neighbor boundary interfaces violate the topological continuity
698 : // assumed by the refinement algorithm. Refinement on such meshes can
699 : // produce invalid neighbor relationships (for example reverse indices
700 : // equal to invalid_uint), which may trigger assertions like
701 : // `rev < neigh->n_neighbors()` or corrupt neighbor data.
702 : //
703 : // For safety, refinement is disabled while disjoint neighbor boundary pairs are
704 : // present.
705 26745 : if (_mesh.get_disjoint_neighbor_boundary_pairs())
706 : {
707 211 : libmesh_not_implemented_msg(
708 : "Mesh refinement is not yet implemented for meshes with disjoint neighbor boundary interfaces.\n");
709 : }
710 :
711 880 : if (_face_level_mismatch_limit)
712 28 : libmesh_assert(test_level_one(true));
713 :
714 : // Possibly clean up the refinement flags from
715 : // a previous step
716 : Threads::parallel_for_each
717 26674 : (_mesh.element_stored_range(),
718 480538 : [](Elem * elem)
719 : {
720 : // Set refinement flag to INACTIVE if the
721 : // element isn't active
722 961076 : if (!elem->active())
723 : {
724 112958 : elem->set_refinement_flag(Elem::INACTIVE);
725 112958 : elem->set_p_refinement_flag(Elem::INACTIVE);
726 : }
727 :
728 : // This might be left over from the last step
729 8877291 : if (elem->refinement_flag() == Elem::JUST_REFINED)
730 6 : elem->set_refinement_flag(Elem::DO_NOTHING);
731 480538 : });
732 :
733 : // Parallel consistency has to come first, or coarsening
734 : // along processor boundaries might occasionally be falsely
735 : // prevented
736 26674 : bool flags_were_consistent = this->make_flags_parallel_consistent();
737 :
738 : // In theory, we should be able to remove the above call, which can
739 : // be expensive and should be unnecessary. In practice, doing
740 : // consistent flagging in parallel is hard, it's impossible to
741 : // verify at the library level if it's being done by user code, and
742 : // we don't want to abort large parallel runs in opt mode... but we
743 : // do want to warn that they should be fixed.
744 880 : libmesh_assert(flags_were_consistent);
745 :
746 880 : if (!flags_were_consistent)
747 : libmesh_warning("Warning: Refinement flags were not consistent between processors! "
748 : "Correcting and continuing.\n");
749 :
750 : // Smooth refinement flags
751 26674 : _smooth_flags(true, false);
752 :
753 : // Now refine the flagged elements. This will
754 : // take up some space, maybe more than what was freed.
755 : const bool mesh_changed =
756 26674 : this->_refine_elements();
757 :
758 880 : if (_face_level_mismatch_limit)
759 28 : libmesh_assert(test_level_one(true));
760 880 : libmesh_assert(this->make_refinement_compatible());
761 :
762 : // FIXME: This won't pass unless we add a redundant find_neighbors()
763 : // call or replace find_neighbors() with on-the-fly neighbor updating
764 : // libmesh_assert(!this->eliminate_unrefined_patches());
765 :
766 : // Finally, the new mesh needs to be prepared for use
767 26674 : if (mesh_changed)
768 2264 : _mesh.prepare_for_use ();
769 :
770 26674 : return mesh_changed;
771 : }
772 :
773 :
774 :
775 97112 : bool MeshRefinement::make_flags_parallel_consistent()
776 : {
777 : // This function must be run on all processors at once
778 2432 : parallel_object_only();
779 :
780 4210 : LOG_SCOPE ("make_flags_parallel_consistent()", "MeshRefinement");
781 :
782 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
783 97112 : &Elem::set_refinement_flag);
784 : Parallel::sync_dofobject_data_by_id
785 191792 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
786 :
787 : SyncRefinementFlags psync(_mesh, &Elem::p_refinement_flag,
788 97112 : &Elem::set_p_refinement_flag);
789 : Parallel::sync_dofobject_data_by_id
790 191792 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), psync);
791 :
792 : // If we weren't consistent in both h and p on every processor then
793 : // we weren't globally consistent
794 192207 : bool parallel_consistent = hsync.parallel_consistent &&
795 96862 : psync.parallel_consistent;
796 97112 : this->comm().min(parallel_consistent);
797 :
798 99544 : return parallel_consistent;
799 : }
800 :
801 :
802 :
803 53943 : bool MeshRefinement::make_coarsening_compatible()
804 : {
805 : // This function must be run on all processors at once
806 3214 : parallel_object_only();
807 :
808 : // We may need a PointLocator for topological_neighbor() tests
809 : // later, which we need to make sure gets constructed on all
810 : // processors at once.
811 52239 : std::unique_ptr<PointLocatorBase> point_locator;
812 :
813 : #ifdef LIBMESH_ENABLE_PERIODIC
814 : bool has_periodic_boundaries =
815 53943 : _periodic_boundaries && !_periodic_boundaries->empty();
816 3214 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
817 :
818 3214 : if (has_periodic_boundaries)
819 532 : point_locator = _mesh.sub_point_locator();
820 : #endif
821 :
822 8132 : LOG_SCOPE ("make_coarsening_compatible()", "MeshRefinement");
823 :
824 : // Unless we encounter a specific situation level-one
825 : // will be satisfied after executing this loop just once
826 3214 : bool level_one_satisfied = true;
827 :
828 :
829 : // Unless we encounter a specific situation we will be compatible
830 : // with any selected refinement flags
831 53943 : bool compatible_with_refinement = true;
832 :
833 :
834 : // find the maximum h and p levels in the mesh
835 53943 : unsigned int max_level = 0;
836 53943 : unsigned int max_p_level = 0;
837 :
838 : {
839 : // First we look at all the active level-0 elements. Since it doesn't make
840 : // sense to coarsen them we must un-set their coarsen flags if
841 : // they are set.
842 30848992 : for (auto & elem : _mesh.active_element_ptr_range())
843 : {
844 16416349 : max_level = std::max(max_level, elem->level());
845 16416349 : max_p_level =
846 16416349 : std::max(max_p_level,
847 17149313 : static_cast<unsigned int>(elem->p_level()));
848 :
849 16429901 : if ((elem->level() == 0) &&
850 123471 : (elem->refinement_flag() == Elem::COARSEN))
851 1040 : elem->set_refinement_flag(Elem::DO_NOTHING);
852 :
853 17766601 : if ((elem->p_level() == 0) &&
854 1350252 : (elem->p_refinement_flag() == Elem::COARSEN))
855 10 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
856 49025 : }
857 : }
858 :
859 : // Even if there are no refined elements on this processor then
860 : // there may still be work for us to do on e.g. ancestor elements.
861 : // At the very least we need to be in the loop if a distributed mesh
862 : // needs to synchronize data.
863 : #if 0
864 : if (max_level == 0 && max_p_level == 0)
865 : {
866 : // But we still have to check with other processors
867 : this->comm().min(compatible_with_refinement);
868 :
869 : return compatible_with_refinement;
870 : }
871 : #endif
872 :
873 : // Loop over all the active elements. If an element is marked
874 : // for coarsening we better check its neighbors. If ANY of these neighbors
875 : // are marked for refinement AND are at the same level then there is a
876 : // conflict. By convention refinement wins, so we un-mark the element for
877 : // coarsening. Level-one would be violated in this case so we need to re-run
878 : // the loop.
879 53943 : if (_face_level_mismatch_limit)
880 : {
881 :
882 27403 : repeat:
883 1882 : level_one_satisfied = true;
884 :
885 376 : do
886 : {
887 2258 : level_one_satisfied = true;
888 :
889 39533586 : for (auto & elem : _mesh.active_element_ptr_range())
890 : {
891 1095900 : bool my_flag_changed = false;
892 :
893 21515557 : if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
894 : // the coarsen flag is set
895 : {
896 4937841 : const unsigned int my_level = elem->level();
897 :
898 21103732 : for (auto n : elem->side_index_range())
899 : {
900 : const Elem * neighbor =
901 16204505 : topological_neighbor(elem, point_locator.get(), n);
902 :
903 16204505 : if (neighbor != nullptr && // I have a
904 15878468 : neighbor != remote_elem) // neighbor here
905 : {
906 1265988 : if (neighbor->active()) // and it is active
907 : {
908 16204938 : if ((neighbor->level() == my_level) &&
909 543220 : (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
910 : // and wants to be refined
911 : {
912 24168 : elem->set_refinement_flag(Elem::DO_NOTHING);
913 1094 : my_flag_changed = true;
914 24168 : break;
915 : }
916 : }
917 : else // I have a neighbor and it is not active. That means it has children.
918 : { // While it _may_ be possible to coarsen us if all the children of
919 : // that element want to be coarsened, it is impossible to know at this
920 : // stage. Forget about it for the moment... This can be handled in
921 : // two steps.
922 201334 : elem->set_refinement_flag(Elem::DO_NOTHING);
923 8904 : my_flag_changed = true;
924 201334 : break;
925 : }
926 : }
927 : }
928 : }
929 21515557 : if (elem->p_refinement_flag() == Elem::COARSEN) // If
930 : // the element is active and the order reduction flag is set
931 : {
932 164 : const unsigned int my_p_level = elem->p_level();
933 :
934 7129 : for (auto n : elem->side_index_range())
935 : {
936 : const Elem * neighbor =
937 5896 : topological_neighbor(elem, point_locator.get(), n);
938 :
939 5896 : if (neighbor != nullptr && // I have a
940 4637 : neighbor != remote_elem) // neighbor here
941 : {
942 388 : if (neighbor->active()) // and it is active
943 : {
944 348 : if ((neighbor->p_level() > my_p_level &&
945 0 : neighbor->p_refinement_flag() != Elem::COARSEN)
946 3772 : || (neighbor->p_level() == my_p_level &&
947 116 : neighbor->p_refinement_flag() == Elem::REFINE))
948 : {
949 217 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
950 14 : my_flag_changed = true;
951 14 : break;
952 : }
953 : }
954 : else // I have a neighbor and it is not active.
955 : { // We need to find which of its children
956 : // have me as a neighbor, and maintain
957 : // level one p compatibility with them.
958 : // Because we currently have level one h
959 : // compatibility, we don't need to check
960 : // grandchildren
961 :
962 20 : libmesh_assert(neighbor->has_children());
963 1616 : for (auto & subneighbor : neighbor->child_ref_range())
964 1513 : if (&subneighbor != remote_elem &&
965 2362 : subneighbor.active() &&
966 973 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
967 598 : if ((subneighbor.p_level() > my_p_level &&
968 34 : subneighbor.p_refinement_flag() != Elem::COARSEN)
969 1014 : || (subneighbor.p_level() == my_p_level &&
970 0 : subneighbor.p_refinement_flag() == Elem::REFINE))
971 : {
972 194 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
973 12 : my_flag_changed = true;
974 12 : break;
975 : }
976 239 : if (my_flag_changed)
977 20 : break;
978 : }
979 : }
980 : }
981 : }
982 :
983 : // If the current element's flag changed, we hadn't
984 : // satisfied the level one rule.
985 20684848 : if (my_flag_changed)
986 10006 : level_one_satisfied = false;
987 :
988 : // Additionally, if it has non-local neighbors, and
989 : // we're not in serial, then we'll eventually have to
990 : // return compatible_with_refinement = false, because
991 : // our change has to propagate to neighboring
992 : // processors.
993 2132007 : if (my_flag_changed && !_mesh.is_serial())
994 15939 : for (auto n : elem->side_index_range())
995 : {
996 : Elem * neigh =
997 14876 : topological_neighbor(elem, point_locator.get(), n);
998 :
999 14876 : if (!neigh)
1000 2927 : continue;
1001 11949 : if (neigh == remote_elem ||
1002 10994 : neigh->processor_id() !=
1003 0 : this->processor_id())
1004 : {
1005 5345 : compatible_with_refinement = false;
1006 5345 : break;
1007 : }
1008 : // FIXME - for non-level one meshes we should
1009 : // test all descendants
1010 0 : if (neigh->has_children())
1011 2074 : for (auto & child : neigh->child_ref_range())
1012 1732 : if (&child == remote_elem ||
1013 1732 : child.processor_id() !=
1014 0 : this->processor_id())
1015 : {
1016 258 : compatible_with_refinement = false;
1017 258 : break;
1018 : }
1019 : }
1020 42216 : }
1021 : }
1022 46074 : while (!level_one_satisfied);
1023 :
1024 : } // end if (_face_level_mismatch_limit)
1025 :
1026 :
1027 : // Next we look at all of the ancestor cells.
1028 : // If there is a parent cell with all of its children
1029 : // wanting to be unrefined then the element is a candidate
1030 : // for unrefinement. If all the children don't
1031 : // all want to be unrefined then ALL of them need to have their
1032 : // unrefinement flags cleared.
1033 326922 : for (int level = max_level; level >= 0; level--)
1034 59530896 : for (auto & elem : as_range(_mesh.level_elements_begin(level), _mesh.level_elements_end(level)))
1035 31179738 : if (elem->ancestor())
1036 : {
1037 : // right now the element hasn't been disqualified
1038 : // as a candidate for unrefinement
1039 509812 : bool is_a_candidate = true;
1040 509812 : bool found_remote_child = false;
1041 :
1042 37613101 : for (auto & child : elem->child_ref_range())
1043 : {
1044 29888768 : if (&child == remote_elem)
1045 40 : found_remote_child = true;
1046 29468375 : else if ((child.refinement_flag() != Elem::COARSEN) ||
1047 129366 : !child.active() )
1048 1939286 : is_a_candidate = false;
1049 : }
1050 :
1051 7410315 : if (!is_a_candidate && !found_remote_child)
1052 : {
1053 6517415 : elem->set_refinement_flag(Elem::INACTIVE);
1054 :
1055 32685915 : for (auto & child : elem->child_ref_range())
1056 : {
1057 26168500 : if (&child == remote_elem)
1058 0 : continue;
1059 26168500 : if (child.refinement_flag() == Elem::COARSEN)
1060 : {
1061 27278 : level_one_satisfied = false;
1062 27278 : child.set_refinement_flag(Elem::DO_NOTHING);
1063 : }
1064 : }
1065 : }
1066 237590 : }
1067 :
1068 63350 : if (!level_one_satisfied && _face_level_mismatch_limit) goto repeat;
1069 :
1070 :
1071 : // If all the children of a parent are set to be coarsened
1072 : // then flag the parent so that they can kill their kids.
1073 :
1074 : // On a distributed mesh, we won't always be able to determine this
1075 : // on parent elements with remote children, even if we own the
1076 : // parent, without communication.
1077 : //
1078 : // We'll first communicate *to* parents' owners when we determine
1079 : // they cannot be coarsened, then we'll sync the final refinement
1080 : // flag *from* the parents.
1081 :
1082 : // uncoarsenable_parents[p] live on processor id p
1083 53943 : const processor_id_type n_proc = _mesh.n_processors();
1084 3214 : const processor_id_type my_proc_id = _mesh.processor_id();
1085 53943 : const bool distributed_mesh = !_mesh.is_replicated();
1086 :
1087 : std::vector<std::vector<dof_id_type>>
1088 53943 : uncoarsenable_parents(n_proc);
1089 :
1090 10376647 : for (auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1091 : {
1092 : // Presume all the children are flagged for coarsening and
1093 : // then look for a contradiction
1094 430656 : bool all_children_flagged_for_coarsening = true;
1095 :
1096 8686418 : for (auto & child : elem->child_ref_range())
1097 : {
1098 8232271 : if (&child != remote_elem &&
1099 470882 : child.refinement_flag() != Elem::COARSEN)
1100 : {
1101 417450 : all_children_flagged_for_coarsening = false;
1102 5844726 : if (!distributed_mesh)
1103 417298 : break;
1104 1187266 : if (child.processor_id() != elem->processor_id())
1105 : {
1106 95755 : uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
1107 0 : break;
1108 : }
1109 : }
1110 : }
1111 :
1112 5347627 : if (all_children_flagged_for_coarsening)
1113 365293 : elem->set_refinement_flag(Elem::COARSEN_INACTIVE);
1114 : else
1115 5078089 : elem->set_refinement_flag(Elem::INACTIVE);
1116 49025 : }
1117 :
1118 : // If we have a distributed mesh, we might need to sync up
1119 : // INACTIVE vs. COARSEN_INACTIVE flags.
1120 53943 : if (distributed_mesh)
1121 : {
1122 : // We'd better still be in sync here
1123 8 : parallel_object_only();
1124 :
1125 : Parallel::MessageTag
1126 29510 : uncoarsenable_tag = this->comm().get_unique_tag();
1127 29498 : std::vector<Parallel::Request> uncoarsenable_push_requests(n_proc-1);
1128 :
1129 343126 : for (processor_id_type p = 0; p != n_proc; ++p)
1130 : {
1131 313628 : if (p == my_proc_id)
1132 29494 : continue;
1133 :
1134 : Parallel::Request &request =
1135 284130 : uncoarsenable_push_requests[p - (p > my_proc_id)];
1136 :
1137 284130 : _mesh.comm().send
1138 284130 : (p, uncoarsenable_parents[p], request, uncoarsenable_tag);
1139 : }
1140 :
1141 313628 : for (processor_id_type p = 1; p != n_proc; ++p)
1142 : {
1143 16 : std::vector<dof_id_type> my_uncoarsenable_parents;
1144 284130 : _mesh.comm().receive
1145 284130 : (Parallel::any_source, my_uncoarsenable_parents,
1146 12 : uncoarsenable_tag);
1147 :
1148 355541 : for (const auto & id : my_uncoarsenable_parents)
1149 : {
1150 71411 : Elem & elem = _mesh.elem_ref(id);
1151 0 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1152 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1153 0 : elem.set_refinement_flag(Elem::INACTIVE);
1154 : }
1155 : }
1156 :
1157 29498 : Parallel::wait(uncoarsenable_push_requests);
1158 :
1159 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1160 29498 : &Elem::set_refinement_flag);
1161 : sync_dofobject_data_by_id
1162 58988 : (this->comm(), _mesh.not_local_elements_begin(),
1163 29510 : _mesh.not_local_elements_end(),
1164 : // We'd like a smaller sync, but this leads to bugs?
1165 : // SyncCoarsenInactive(),
1166 : hsync);
1167 29486 : }
1168 :
1169 : // If one processor finds an incompatibility, we're globally
1170 : // incompatible
1171 53943 : this->comm().min(compatible_with_refinement);
1172 :
1173 107886 : return compatible_with_refinement;
1174 49025 : }
1175 :
1176 :
1177 :
1178 :
1179 :
1180 :
1181 :
1182 :
1183 55006 : bool MeshRefinement::make_refinement_compatible()
1184 : {
1185 : // This function must be run on all processors at once
1186 3264 : parallel_object_only();
1187 :
1188 : // We may need a PointLocator for topological_neighbor() tests
1189 : // later, which we need to make sure gets constructed on all
1190 : // processors at once.
1191 53276 : std::unique_ptr<PointLocatorBase> point_locator;
1192 :
1193 : #ifdef LIBMESH_ENABLE_PERIODIC
1194 : bool has_periodic_boundaries =
1195 55006 : _periodic_boundaries && !_periodic_boundaries->empty();
1196 3264 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
1197 :
1198 3264 : if (has_periodic_boundaries)
1199 532 : point_locator = _mesh.sub_point_locator();
1200 : #endif
1201 :
1202 3264 : LOG_SCOPE ("make_refinement_compatible()", "MeshRefinement");
1203 :
1204 : // Unless we encounter a specific situation we will be compatible
1205 : // with any selected coarsening flags
1206 55006 : bool compatible_with_coarsening = true;
1207 :
1208 : // This loop enforces the level-1 rule. We should only
1209 : // execute it if the user indeed wants level-1 satisfied!
1210 55006 : if (_face_level_mismatch_limit)
1211 : {
1212 : // Unless we encounter a specific situation level-one
1213 : // will be satisfied after executing this loop just once
1214 1560 : bool level_one_satisfied = true;
1215 :
1216 218 : do
1217 : {
1218 1778 : level_one_satisfied = true;
1219 :
1220 25447222 : for (auto & elem : _mesh.active_element_ptr_range())
1221 : {
1222 13367727 : const unsigned short n_sides = elem->n_sides();
1223 :
1224 13900045 : if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1225 : // h refinement flag is set
1226 : {
1227 443819 : const unsigned int my_level = elem->level();
1228 :
1229 2159944 : for (unsigned short side = 0; side != n_sides;
1230 : ++side)
1231 : {
1232 : Elem * neighbor =
1233 1716125 : topological_neighbor(elem, point_locator.get(), side);
1234 :
1235 1721192 : if (neighbor != nullptr && // I have a
1236 1798469 : neighbor != remote_elem && // neighbor here
1237 164688 : neighbor->active()) // and it is active
1238 : {
1239 : // Case 1: The neighbor is at the same level I am.
1240 : // 1a: The neighbor will be refined -> NO PROBLEM
1241 : // 1b: The neighbor won't be refined -> NO PROBLEM
1242 : // 1c: The neighbor wants to be coarsened -> PROBLEM
1243 1223504 : if (neighbor->level() == my_level)
1244 : {
1245 1105677 : if (neighbor->refinement_flag() == Elem::COARSEN)
1246 : {
1247 2 : neighbor->set_refinement_flag(Elem::DO_NOTHING);
1248 79 : if (neighbor->parent())
1249 2 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1250 77 : compatible_with_coarsening = false;
1251 2 : level_one_satisfied = false;
1252 : }
1253 : }
1254 :
1255 :
1256 : // Case 2: The neighbor is one level lower than I am.
1257 : // The neighbor thus MUST be refined to satisfy
1258 : // the level-one rule, regardless of whether it
1259 : // was originally flagged for refinement. If it
1260 : // wasn't flagged already we need to repeat
1261 : // this process.
1262 117827 : else if ((neighbor->level()+1) == my_level)
1263 : {
1264 117827 : if (neighbor->refinement_flag() != Elem::REFINE)
1265 : {
1266 620 : neighbor->set_refinement_flag(Elem::REFINE);
1267 13349 : if (neighbor->parent())
1268 460 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1269 12729 : compatible_with_coarsening = false;
1270 620 : level_one_satisfied = false;
1271 : }
1272 : }
1273 : #ifdef DEBUG
1274 : // Note that the only other possibility is that the
1275 : // neighbor is already refined, in which case it isn't
1276 : // active and we should never get here.
1277 : else
1278 0 : libmesh_error_msg("ERROR: Neighbor level must be equal or 1 higher than mine.");
1279 : #endif
1280 : }
1281 : }
1282 : }
1283 13900045 : if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1284 : // p refinement flag is set
1285 : {
1286 2784 : const unsigned int my_p_level = elem->p_level();
1287 :
1288 129835 : for (unsigned int side=0; side != n_sides; side++)
1289 : {
1290 : Elem * neighbor =
1291 109604 : topological_neighbor(elem, point_locator.get(), side);
1292 :
1293 109604 : if (neighbor != nullptr && // I have a
1294 90721 : neighbor != remote_elem) // neighbor here
1295 : {
1296 11884 : if (neighbor->active()) // and it is active
1297 : {
1298 86179 : if (neighbor->p_level() < my_p_level &&
1299 18 : neighbor->p_refinement_flag() != Elem::REFINE)
1300 : {
1301 2 : neighbor->set_p_refinement_flag(Elem::REFINE);
1302 2 : level_one_satisfied = false;
1303 24 : compatible_with_coarsening = false;
1304 : }
1305 91963 : if (neighbor->p_level() == my_p_level &&
1306 5802 : neighbor->p_refinement_flag() == Elem::COARSEN)
1307 : {
1308 0 : neighbor->set_p_refinement_flag(Elem::DO_NOTHING);
1309 0 : level_one_satisfied = false;
1310 0 : compatible_with_coarsening = false;
1311 : }
1312 : }
1313 : else // I have an inactive neighbor
1314 : {
1315 50 : libmesh_assert(neighbor->has_children());
1316 1728 : for (auto & subneighbor : neighbor->child_ref_range())
1317 2099 : if (&subneighbor != remote_elem && subneighbor.active() &&
1318 947 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
1319 : {
1320 664 : if (subneighbor.p_level() < my_p_level &&
1321 38 : subneighbor.p_refinement_flag() != Elem::REFINE)
1322 : {
1323 : // We should already be level one
1324 : // compatible
1325 6 : libmesh_assert_greater (subneighbor.p_level() + 2u,
1326 : my_p_level);
1327 6 : subneighbor.set_p_refinement_flag(Elem::REFINE);
1328 6 : level_one_satisfied = false;
1329 68 : compatible_with_coarsening = false;
1330 : }
1331 638 : if (subneighbor.p_level() == my_p_level &&
1332 12 : subneighbor.p_refinement_flag() == Elem::COARSEN)
1333 : {
1334 0 : subneighbor.set_p_refinement_flag(Elem::DO_NOTHING);
1335 0 : level_one_satisfied = false;
1336 0 : compatible_with_coarsening = false;
1337 : }
1338 : }
1339 : }
1340 : }
1341 : }
1342 : }
1343 30740 : }
1344 : }
1345 :
1346 33614 : while (!level_one_satisfied);
1347 : } // end if (_face_level_mismatch_limit)
1348 :
1349 : // If we're not compatible on one processor, we're globally not
1350 : // compatible
1351 55006 : this->comm().min(compatible_with_coarsening);
1352 :
1353 60000 : return compatible_with_coarsening;
1354 50012 : }
1355 :
1356 :
1357 :
1358 :
1359 45303 : bool MeshRefinement::_coarsen_elements ()
1360 : {
1361 : // This function must be run on all processors at once
1362 1542 : parallel_object_only();
1363 :
1364 3084 : LOG_SCOPE ("_coarsen_elements()", "MeshRefinement");
1365 :
1366 : // Flags indicating if this call actually changes the mesh
1367 45303 : bool mesh_changed = false;
1368 45303 : bool mesh_p_changed = false;
1369 :
1370 : // Clear the unused_elements data structure.
1371 : // The elements have been packed since it was built,
1372 : // so there are _no_ unused elements. We cannot trust
1373 : // any iterators currently in this data structure.
1374 : // _unused_elements.clear();
1375 :
1376 : // Loop over the elements first to determine if the mesh will
1377 : // undergo h-coarsening. If it will, then we'll need to communicate
1378 : // more ghosted elements. We need to communicate them *before* we
1379 : // do the coarsening; otherwise it is possible to coarsen away a
1380 : // one-element-thick layer partition and leave the partitions on
1381 : // either side unable to figure out how to talk to each other.
1382 24182046 : for (auto & elem : _mesh.element_ptr_range())
1383 13356040 : if (elem->refinement_flag() == Elem::COARSEN)
1384 : {
1385 7028 : mesh_changed = true;
1386 7028 : break;
1387 42219 : }
1388 :
1389 : // If the mesh changed on any processor, it changed globally
1390 45303 : this->comm().max(mesh_changed);
1391 :
1392 : // And then we may need to widen the ghosting layers.
1393 45303 : if (mesh_changed)
1394 7121 : MeshCommunication().send_coarse_ghosts(_mesh);
1395 :
1396 35960810 : for (auto & elem : _mesh.element_ptr_range())
1397 : {
1398 : // Make sure we transfer the children's boundary id(s)
1399 : // up to its parent when necessary before coarsening.
1400 19779137 : _mesh.get_boundary_info().transfer_boundary_ids_from_children(elem);
1401 :
1402 : // active elements flagged for coarsening will
1403 : // no longer be deleted until MeshRefinement::contract()
1404 19779137 : if (elem->refinement_flag() == Elem::COARSEN)
1405 : {
1406 : // Huh? no level-0 element should be active
1407 : // and flagged for coarsening.
1408 0 : libmesh_assert_not_equal_to (elem->level(), 0);
1409 :
1410 : // Remove this element from any neighbor
1411 : // lists that point to it.
1412 0 : elem->nullify_neighbors();
1413 :
1414 : // Remove any boundary information associated
1415 : // with this element if we do not allow children to have boundary info.
1416 : // Otherwise, we will do the removal in `transfer_boundary_ids_from_children`
1417 : // to make sure we don't delete the information before it is transferred
1418 0 : if (!_mesh.get_boundary_info().is_children_on_boundary_side())
1419 0 : _mesh.get_boundary_info().remove (elem);
1420 :
1421 : // Add this iterator to the _unused_elements
1422 : // data structure so we might fill it.
1423 : // The _unused_elements optimization is currently off.
1424 : // _unused_elements.push_back (it);
1425 :
1426 : // Don't delete the element until
1427 : // MeshRefinement::contract()
1428 : // _mesh.delete_elem(elem);
1429 : }
1430 :
1431 : // inactive elements flagged for coarsening
1432 : // will become active
1433 18857505 : else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1434 : {
1435 395420 : elem->coarsen();
1436 18488 : libmesh_assert (elem->active());
1437 :
1438 : // the mesh has certainly changed
1439 395420 : mesh_changed = true;
1440 : }
1441 19779137 : if (elem->p_refinement_flag() == Elem::COARSEN)
1442 : {
1443 220 : if (elem->p_level() > 0)
1444 : {
1445 14 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1446 206 : elem->set_p_level(elem->p_level() - 1);
1447 206 : mesh_p_changed = true;
1448 : }
1449 : else
1450 : {
1451 0 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
1452 : }
1453 : }
1454 42219 : }
1455 :
1456 45303 : this->comm().max(mesh_p_changed);
1457 :
1458 : // And we may need to update DistributedMesh values reflecting the changes
1459 45303 : if (mesh_changed)
1460 7121 : _mesh.update_parallel_id_counts();
1461 :
1462 : // Node processor ids may need to change if an element of that id
1463 : // was coarsened away
1464 45303 : if (mesh_changed && !_mesh.is_serial())
1465 : {
1466 : // Update the _new_nodes_map so that processors can
1467 : // find requested nodes
1468 1228 : this->update_nodes_map ();
1469 :
1470 1228 : MeshCommunication().make_nodes_parallel_consistent (_mesh);
1471 :
1472 : // Clear the _new_nodes_map
1473 1228 : this->clear();
1474 :
1475 : #ifdef DEBUG
1476 4 : MeshTools::libmesh_assert_valid_procids<Node>(_mesh);
1477 : #endif
1478 : }
1479 :
1480 : // If p levels changed all we need to do is make sure that parent p
1481 : // levels changed in sync
1482 45303 : if (mesh_p_changed && !_mesh.is_serial())
1483 : {
1484 192 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1485 : }
1486 :
1487 48387 : return (mesh_changed || mesh_p_changed);
1488 : }
1489 :
1490 :
1491 :
1492 61823 : bool MeshRefinement::_refine_elements ()
1493 : {
1494 2096 : libmesh_assert(_mesh.is_prepared() || _mesh.is_replicated());
1495 :
1496 : // This function must be run on all processors at once
1497 2096 : parallel_object_only();
1498 :
1499 : // Update the _new_nodes_map so that elements can
1500 : // find nodes to connect to.
1501 61823 : this->update_nodes_map ();
1502 :
1503 6288 : LOG_SCOPE ("_refine_elements()", "MeshRefinement");
1504 :
1505 : // Iterate over the elements, counting the elements
1506 : // flagged for h refinement.
1507 2096 : dof_id_type n_elems_flagged = 0;
1508 :
1509 35186460 : for (auto & elem : _mesh.element_ptr_range())
1510 19342191 : if (elem->refinement_flag() == Elem::REFINE)
1511 1302041 : n_elems_flagged++;
1512 :
1513 : // Construct a local vector of Elem * which have been
1514 : // previously marked for refinement. We reserve enough
1515 : // space to allow for every element to be refined.
1516 2096 : std::vector<Elem *> local_copy_of_elements;
1517 61823 : local_copy_of_elements.reserve(n_elems_flagged);
1518 :
1519 : // If mesh p levels changed, we might need to synchronize parent p
1520 : // levels on a distributed mesh.
1521 61823 : bool mesh_p_changed = false;
1522 :
1523 : // Iterate over the elements, looking for elements flagged for
1524 : // refinement.
1525 :
1526 : // If we are on a ReplicatedMesh, then we just do the refinement in
1527 : // the same order on every processor and everything stays in sync.
1528 :
1529 : // If we are on a DistributedMesh, that's impossible.
1530 : //
1531 : // If the mesh is distributed, we need to make sure that if we end
1532 : // up as the owner of a new node, which might happen if that node is
1533 : // attached to one of our own elements, then we have given it a
1534 : // legitimate node id and our own processor id. We generate
1535 : // legitimate node ids and use our own processor id when we are
1536 : // refining our own elements but not when we refine others'
1537 : // elements. Therefore we want to refine our own elements *first*,
1538 : // thereby generating all nodes which might belong to us, and then
1539 : // refine others' elements *after*, thereby generating nodes with
1540 : // temporary ids which we know we will discard.
1541 : //
1542 : // Even if the DistributedMesh is serialized, we can't just treat it
1543 : // like a ReplicatedMesh, because DistributedMesh doesn't *trust*
1544 : // users to refine partitioned elements in a serialized way, so it
1545 : // assigns temporary ids, so we need to synchronize ids afterward to
1546 : // be safe anyway, so we might as well use the distributed mesh code
1547 : // path.
1548 23075692 : for (auto & elem : _mesh.is_replicated() ? _mesh.active_element_ptr_range() : _mesh.active_local_element_ptr_range())
1549 : {
1550 12780204 : if (elem->refinement_flag() == Elem::REFINE)
1551 832582 : local_copy_of_elements.push_back(elem);
1552 12781502 : if (elem->p_refinement_flag() == Elem::REFINE &&
1553 2596 : elem->active())
1554 : {
1555 11033 : elem->set_p_level(elem->p_level()+1);
1556 9735 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1557 9735 : mesh_p_changed = true;
1558 : }
1559 57631 : }
1560 :
1561 61823 : if (!_mesh.is_replicated())
1562 : {
1563 77424 : for (auto & elem : as_range(_mesh.active_not_local_elements_begin(),
1564 1625704 : _mesh.active_not_local_elements_end()))
1565 : {
1566 736058 : if (elem->refinement_flag() == Elem::REFINE)
1567 411828 : local_copy_of_elements.push_back(elem);
1568 736058 : if (elem->p_refinement_flag() == Elem::REFINE &&
1569 0 : elem->active())
1570 : {
1571 8163 : elem->set_p_level(elem->p_level()+1);
1572 8163 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1573 8163 : mesh_p_changed = true;
1574 : }
1575 38523 : }
1576 : }
1577 :
1578 : // Now iterate over the local copies and refine each one.
1579 : // This may resize the mesh's internal container and invalidate
1580 : // any existing iterators.
1581 1306233 : for (auto & elem : local_copy_of_elements)
1582 1244410 : elem->refine(*this);
1583 :
1584 : // The mesh changed if there were elements h refined
1585 61823 : bool mesh_changed = !local_copy_of_elements.empty();
1586 :
1587 : // If the mesh changed on any processor, it changed globally
1588 61823 : this->comm().max(mesh_changed);
1589 61823 : this->comm().max(mesh_p_changed);
1590 :
1591 : // And we may need to update DistributedMesh values reflecting the changes
1592 61823 : if (mesh_changed)
1593 35708 : _mesh.update_parallel_id_counts();
1594 :
1595 61823 : if (mesh_changed && !_mesh.is_replicated())
1596 : {
1597 22593 : MeshCommunication().make_elems_parallel_consistent (_mesh);
1598 22593 : MeshCommunication().make_new_nodes_parallel_consistent (_mesh);
1599 : #ifdef DEBUG
1600 38 : _mesh.libmesh_assert_valid_parallel_ids();
1601 : #endif
1602 : }
1603 :
1604 : // If we're refining a ReplicatedMesh, then we haven't yet assigned
1605 : // node processor ids. But if we're refining a partitioned
1606 : // ReplicatedMesh, then we *need* to assign node processor ids.
1607 64203 : if (mesh_changed && _mesh.is_replicated() &&
1608 74938 : (_mesh.unpartitioned_elements_begin() ==
1609 86863 : _mesh.unpartitioned_elements_end()))
1610 12724 : Partitioner::set_node_processor_ids(_mesh);
1611 :
1612 61823 : if (mesh_p_changed && !_mesh.is_replicated())
1613 : {
1614 2368 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1615 : }
1616 :
1617 : // Clear the _new_nodes_map and _unused_elements data structures.
1618 61823 : this->clear();
1619 :
1620 123646 : return (mesh_changed || mesh_p_changed);
1621 : }
1622 :
1623 :
1624 70845 : void MeshRefinement::_smooth_flags(bool refining, bool coarsening)
1625 : {
1626 : // Smoothing can break in weird ways on a mesh with broken topology
1627 : #ifdef DEBUG
1628 2390 : MeshTools::libmesh_assert_valid_neighbors(_mesh);
1629 : #endif
1630 :
1631 : // Repeat until flag changes match on every processor
1632 0 : do
1633 : {
1634 : // Repeat until coarsening & refinement flags jive
1635 2390 : bool satisfied = false;
1636 196 : do
1637 : {
1638 : // If we're refining or coarsening, hit the corresponding
1639 : // face level test code. Short-circuiting || is our friend
1640 : const bool coarsening_satisfied =
1641 131735 : !coarsening ||
1642 52433 : this->make_coarsening_compatible();
1643 :
1644 : const bool refinement_satisfied =
1645 132774 : !refining ||
1646 53472 : this->make_refinement_compatible();
1647 :
1648 : bool smoothing_satisfied =
1649 79302 : !this->eliminate_unrefined_patches();// &&
1650 :
1651 79302 : if (_edge_level_mismatch_limit)
1652 0 : smoothing_satisfied = smoothing_satisfied &&
1653 0 : !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);
1654 :
1655 79302 : if (_node_level_mismatch_limit)
1656 0 : smoothing_satisfied = smoothing_satisfied &&
1657 0 : !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);
1658 :
1659 79302 : if (_overrefined_boundary_limit>=0)
1660 54134 : smoothing_satisfied = smoothing_satisfied &&
1661 26137 : !this->limit_overrefined_boundary(_overrefined_boundary_limit);
1662 :
1663 79302 : if (_underrefined_boundary_limit>=0)
1664 53566 : smoothing_satisfied = smoothing_satisfied &&
1665 25569 : !this->limit_underrefined_boundary(_underrefined_boundary_limit);
1666 :
1667 79302 : satisfied = (coarsening_satisfied &&
1668 81888 : refinement_satisfied &&
1669 : smoothing_satisfied);
1670 :
1671 2586 : libmesh_assert(this->comm().verify(satisfied));
1672 : }
1673 76716 : while (!satisfied);
1674 : }
1675 137332 : while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
1676 70845 : }
1677 :
1678 :
1679 0 : void MeshRefinement::uniformly_p_refine (unsigned int n)
1680 : {
1681 : // Refine n times
1682 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1683 0 : for (auto & elem : _mesh.active_element_ptr_range())
1684 : {
1685 : // P refine all the active elements
1686 0 : elem->set_p_level(elem->p_level()+1);
1687 0 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1688 0 : }
1689 0 : }
1690 :
1691 :
1692 :
1693 0 : void MeshRefinement::uniformly_p_coarsen (unsigned int n)
1694 : {
1695 : // Coarsen p times
1696 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1697 0 : for (auto & elem : _mesh.active_element_ptr_range())
1698 0 : if (elem->p_level() > 0)
1699 : {
1700 : // P coarsen all the active elements
1701 0 : elem->set_p_level(elem->p_level()-1);
1702 0 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1703 0 : }
1704 0 : }
1705 :
1706 :
1707 :
1708 14558 : void MeshRefinement::uniformly_refine (unsigned int n)
1709 : {
1710 : // Refine n times
1711 : // FIXME - this won't work if n>1 and the mesh
1712 : // has already been attached to an equation system
1713 31295 : for (unsigned int rstep=0; rstep<n; rstep++)
1714 : {
1715 : // Clean up the refinement flags
1716 16737 : this->clean_refinement_flags();
1717 :
1718 : // Flag all the active elements for refinement.
1719 1938158 : for (auto & elem : _mesh.active_element_ptr_range())
1720 1029496 : elem->set_refinement_flag(Elem::REFINE);
1721 :
1722 : // Refine all the elements we just flagged.
1723 16737 : this->_refine_elements();
1724 : }
1725 :
1726 : // Finally, the new mesh probably needs to be prepared for use
1727 14558 : if (n > 0)
1728 13373 : _mesh.prepare_for_use ();
1729 14558 : }
1730 :
1731 :
1732 :
1733 1132 : void MeshRefinement::uniformly_coarsen (unsigned int n)
1734 : {
1735 : // Coarsen n times
1736 2264 : for (unsigned int rstep=0; rstep<n; rstep++)
1737 : {
1738 : // Clean up the refinement flags
1739 1132 : this->clean_refinement_flags();
1740 :
1741 : // Flag all the active elements for coarsening.
1742 471746 : for (auto & elem : _mesh.active_element_ptr_range())
1743 : {
1744 260509 : elem->set_refinement_flag(Elem::COARSEN);
1745 286261 : if (elem->parent())
1746 25752 : elem->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
1747 1068 : }
1748 :
1749 : // On a distributed mesh, we may have parent elements with
1750 : // remote active children. To keep flags consistent, we'll need
1751 : // a communication step.
1752 1132 : if (!_mesh.is_replicated())
1753 : {
1754 908 : const processor_id_type n_proc = _mesh.n_processors();
1755 4 : const processor_id_type my_proc_id = _mesh.processor_id();
1756 :
1757 : std::vector<std::vector<dof_id_type>>
1758 916 : parents_to_coarsen(n_proc);
1759 :
1760 117848 : for (const auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1761 57878 : if (elem->processor_id() != my_proc_id &&
1762 48 : elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1763 14280 : parents_to_coarsen[elem->processor_id()].push_back(elem->id());
1764 :
1765 : Parallel::MessageTag
1766 916 : coarsen_tag = this->comm().get_unique_tag();
1767 908 : std::vector<Parallel::Request> coarsen_push_requests(n_proc-1);
1768 :
1769 10452 : for (processor_id_type p = 0; p != n_proc; ++p)
1770 : {
1771 9544 : if (p == my_proc_id)
1772 904 : continue;
1773 :
1774 : Parallel::Request &request =
1775 8636 : coarsen_push_requests[p - (p > my_proc_id)];
1776 :
1777 8636 : _mesh.comm().send
1778 8636 : (p, parents_to_coarsen[p], request, coarsen_tag);
1779 : }
1780 :
1781 9544 : for (processor_id_type p = 1; p != n_proc; ++p)
1782 : {
1783 8 : std::vector<dof_id_type> my_parents_to_coarsen;
1784 8636 : _mesh.comm().receive
1785 8636 : (Parallel::any_source, my_parents_to_coarsen,
1786 8 : coarsen_tag);
1787 :
1788 21980 : for (const auto & id : my_parents_to_coarsen)
1789 : {
1790 13344 : Elem & elem = _mesh.elem_ref(id);
1791 36 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1792 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1793 36 : elem.set_refinement_flag(Elem::COARSEN_INACTIVE);
1794 : }
1795 : }
1796 :
1797 908 : Parallel::wait(coarsen_push_requests);
1798 :
1799 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1800 908 : &Elem::set_refinement_flag);
1801 : sync_dofobject_data_by_id
1802 1812 : (this->comm(), _mesh.not_local_elements_begin(),
1803 916 : _mesh.not_local_elements_end(),
1804 : // We'd like a smaller sync, but this leads to bugs?
1805 : // SyncCoarsenInactive(),
1806 : hsync);
1807 900 : }
1808 :
1809 : // Coarsen all the elements we just flagged.
1810 1132 : this->_coarsen_elements();
1811 : }
1812 :
1813 :
1814 : // Finally, the new mesh probably needs to be prepared for use
1815 1132 : if (n > 0)
1816 1132 : _mesh.prepare_for_use ();
1817 1132 : }
1818 :
1819 :
1820 :
1821 19127040 : Elem * MeshRefinement::topological_neighbor(Elem * elem,
1822 : const PointLocatorBase * point_locator,
1823 : const unsigned int side) const
1824 : {
1825 : #ifdef LIBMESH_ENABLE_PERIODIC
1826 19127040 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1827 : {
1828 488472 : libmesh_assert(point_locator);
1829 786435 : return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1830 : }
1831 : #endif
1832 18880191 : return elem->neighbor_ptr(side);
1833 : }
1834 :
1835 :
1836 :
1837 1920 : bool MeshRefinement::has_topological_neighbor(const Elem * elem,
1838 : const PointLocatorBase * point_locator,
1839 : const Elem * neighbor) const
1840 : {
1841 : #ifdef LIBMESH_ENABLE_PERIODIC
1842 1920 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1843 : {
1844 0 : libmesh_assert(point_locator);
1845 0 : return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1846 : }
1847 : #endif
1848 1772 : return elem->has_neighbor(neighbor);
1849 : }
1850 :
1851 :
1852 :
1853 : } // namespace libMesh
1854 :
1855 :
1856 : #endif
|
