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