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 299660 : MeshRefinement::MeshRefinement (MeshBase & m) :
95 : ParallelObject(m),
96 282192 : _mesh(m),
97 282192 : _use_member_parameters(false),
98 282192 : _coarsen_by_parents(false),
99 282192 : _refine_fraction(0.3),
100 282192 : _coarsen_fraction(0.0),
101 282192 : _max_h_level(libMesh::invalid_uint),
102 282192 : _coarsen_threshold(10),
103 282192 : _nelem_target(0),
104 282192 : _absolute_global_tolerance(0.0),
105 282192 : _face_level_mismatch_limit(1),
106 282192 : _edge_level_mismatch_limit(0),
107 282192 : _node_level_mismatch_limit(0),
108 282192 : _overrefined_boundary_limit(0),
109 282192 : _underrefined_boundary_limit(0),
110 282192 : _allow_unrefined_patches(false),
111 282192 : _enforce_mismatch_limit_prior_to_refinement(false)
112 : #ifdef LIBMESH_ENABLE_PERIODIC
113 308394 : , _periodic_boundaries(nullptr)
114 : #endif
115 : {
116 299660 : }
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 283659 : MeshRefinement::~MeshRefinement () = default;
130 :
131 :
132 :
133 54499 : void MeshRefinement::clear ()
134 : {
135 1858 : _new_nodes_map.clear();
136 54499 : }
137 :
138 :
139 :
140 45223509 : Node * MeshRefinement::add_node(Elem & parent,
141 : unsigned int child,
142 : unsigned int node,
143 : processor_id_type proc_id)
144 : {
145 5208048 : LOG_SCOPE("add_node()", "MeshRefinement");
146 :
147 45223509 : unsigned int parent_n = parent.as_parent_node(child, node);
148 :
149 45223509 : if (parent_n != libMesh::invalid_uint)
150 13997414 : return parent.node_ptr(parent_n);
151 :
152 : const std::vector<std::pair<dof_id_type, dof_id_type>>
153 33830127 : 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 1860164 : 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 31969963 : if (const auto new_node_id = _new_nodes_map.find(bracketing_nodes);
165 : new_node_id != DofObject::invalid_id)
166 23050857 : 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 507877 : Point p; // defaults to 0,0,0
173 :
174 111475166 : for (auto n : parent.node_index_range())
175 : {
176 : // The value from the embedding matrix
177 102556060 : const Real em_val = parent.embedding_matrix(child,node,n);
178 :
179 102556060 : 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 2665706 : 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 8919106 : Node * new_node = _mesh.add_point (p, DofObject::invalid_id, proc_id);
192 :
193 507877 : libmesh_assert(new_node);
194 :
195 : // But then we'll make sure this node is marked as unpartitioned.
196 8919106 : new_node->processor_id() = DofObject::invalid_processor_id;
197 :
198 : // Add the node to the map.
199 8919106 : _new_nodes_map.add_node(*new_node, bracketing_nodes);
200 :
201 : // Return the address of the new node
202 507877 : 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 6223269 : Elem * MeshRefinement::add_elem (std::unique_ptr<Elem> elem)
214 : {
215 380306 : libmesh_assert(elem);
216 6983889 : 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 121903 : while (parent)
260 : {
261 22576 : parent = parent->parent();
262 107691 : if (parent)
263 : {
264 14584 : const dof_id_type parentid = parent->id();
265 7292 : libmesh_assert_less (parentid, error_per_parent.size());
266 78719 : 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 38400 : error_per_parent[parentid] += (error_per_cell[elem->id()] *
301 3452 : 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 257358 : error_per_parent[i] = -1.;
322 257358 : 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 45722 : 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 43980 : error_per_parent[i] = std::sqrt(error_per_parent[i]);
336 :
337 43980 : parent_error_min = std::min (parent_error_min,
338 1742 : error_per_parent[i]);
339 50691 : parent_error_max = std::max (parent_error_max,
340 1742 : 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 251928 : for (auto & elem : _mesh.active_local_element_ptr_range())
380 1222904 : for (auto n : elem->side_index_range())
381 : {
382 : Elem * neighbor =
383 972172 : topological_neighbor(elem, point_locator.get(), n);
384 :
385 1849199 : if (!neighbor || !neighbor->active() ||
386 877027 : neighbor == remote_elem)
387 95145 : continue;
388 :
389 877027 : if ((neighbor->level() + 1 < elem->level()) ||
390 1754054 : (neighbor->p_level() + 1 < elem->p_level()) ||
391 877027 : (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 118930 : for (auto & elem : _mesh.active_local_element_ptr_range())
441 118376 : if (elem->refinement_flag() == Elem::REFINE ||
442 118376 : elem->refinement_flag() == Elem::COARSEN ||
443 355128 : elem->p_refinement_flag() == Elem::REFINE ||
444 118376 : 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 42880 : if (flag == Elem::REFINE || flag == Elem::COARSEN)
507 11275 : elements_flagged = true;
508 : else
509 : {
510 : const Elem::RefinementState pflag =
511 2228 : elem->p_refinement_flag();
512 31605 : if (pflag == Elem::REFINE || pflag == Elem::COARSEN)
513 635 : 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 16522594 : for (auto & elem : _mesh.element_ptr_range())
620 : {
621 : // Set refinement flag to INACTIVE if the
622 : // element isn't active
623 8737570 : if (!elem->active())
624 : {
625 169648 : elem->set_refinement_flag(Elem::INACTIVE);
626 169648 : elem->set_p_refinement_flag(Elem::INACTIVE);
627 : }
628 :
629 : // This might be left over from the last step
630 8737570 : if (elem->refinement_flag() == Elem::JUST_REFINED)
631 131004 : 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 23032 : bool MeshRefinement::refine_elements ()
683 : {
684 : // This function must be run on all processors at once
685 776 : parallel_object_only();
686 :
687 776 : if (_face_level_mismatch_limit)
688 26 : libmesh_assert(test_level_one(true));
689 :
690 : // Possibly clean up the refinement flags from
691 : // a previous step
692 14248430 : for (auto & elem : _mesh.element_ptr_range())
693 : {
694 : // Set refinement flag to INACTIVE if the
695 : // element isn't active
696 7544251 : if (!elem->active())
697 : {
698 103748 : elem->set_refinement_flag(Elem::INACTIVE);
699 103748 : elem->set_p_refinement_flag(Elem::INACTIVE);
700 : }
701 :
702 : // This might be left over from the last step
703 7544251 : if (elem->refinement_flag() == Elem::JUST_REFINED)
704 6 : elem->set_refinement_flag(Elem::DO_NOTHING);
705 21480 : }
706 :
707 : // Parallel consistency has to come first, or coarsening
708 : // along processor boundaries might occasionally be falsely
709 : // prevented
710 23032 : bool flags_were_consistent = this->make_flags_parallel_consistent();
711 :
712 : // In theory, we should be able to remove the above call, which can
713 : // be expensive and should be unnecessary. In practice, doing
714 : // consistent flagging in parallel is hard, it's impossible to
715 : // verify at the library level if it's being done by user code, and
716 : // we don't want to abort large parallel runs in opt mode... but we
717 : // do want to warn that they should be fixed.
718 776 : libmesh_assert(flags_were_consistent);
719 :
720 776 : if (!flags_were_consistent)
721 : libmesh_warning("Warning: Refinement flags were not consistent between processors! "
722 : "Correcting and continuing.\n");
723 :
724 : // Smooth refinement flags
725 23032 : _smooth_flags(true, false);
726 :
727 : // Now refine the flagged elements. This will
728 : // take up some space, maybe more than what was freed.
729 : const bool mesh_changed =
730 23032 : this->_refine_elements();
731 :
732 776 : if (_face_level_mismatch_limit)
733 26 : libmesh_assert(test_level_one(true));
734 776 : libmesh_assert(this->make_refinement_compatible());
735 :
736 : // FIXME: This won't pass unless we add a redundant find_neighbors()
737 : // call or replace find_neighbors() with on-the-fly neighbor updating
738 : // libmesh_assert(!this->eliminate_unrefined_patches());
739 :
740 : // Finally, the new mesh needs to be prepared for use
741 23032 : if (mesh_changed)
742 2193 : _mesh.prepare_for_use ();
743 :
744 23032 : return mesh_changed;
745 : }
746 :
747 :
748 :
749 89607 : bool MeshRefinement::make_flags_parallel_consistent()
750 : {
751 : // This function must be run on all processors at once
752 2126 : parallel_object_only();
753 :
754 2126 : LOG_SCOPE ("make_flags_parallel_consistent()", "MeshRefinement");
755 :
756 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
757 89607 : &Elem::set_refinement_flag);
758 : Parallel::sync_dofobject_data_by_id
759 177088 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
760 :
761 : SyncRefinementFlags psync(_mesh, &Elem::p_refinement_flag,
762 89607 : &Elem::set_p_refinement_flag);
763 : Parallel::sync_dofobject_data_by_id
764 177088 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), psync);
765 :
766 : // If we weren't consistent in both h and p on every processor then
767 : // we weren't globally consistent
768 177404 : bool parallel_consistent = hsync.parallel_consistent &&
769 89358 : psync.parallel_consistent;
770 89607 : this->comm().min(parallel_consistent);
771 :
772 91733 : return parallel_consistent;
773 : }
774 :
775 :
776 :
777 45499 : bool MeshRefinement::make_coarsening_compatible()
778 : {
779 : // This function must be run on all processors at once
780 2778 : parallel_object_only();
781 :
782 : // We may need a PointLocator for topological_neighbor() tests
783 : // later, which we need to make sure gets constructed on all
784 : // processors at once.
785 44029 : std::unique_ptr<PointLocatorBase> point_locator;
786 :
787 : #ifdef LIBMESH_ENABLE_PERIODIC
788 : bool has_periodic_boundaries =
789 45499 : _periodic_boundaries && !_periodic_boundaries->empty();
790 2778 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
791 :
792 2778 : if (has_periodic_boundaries)
793 532 : point_locator = _mesh.sub_point_locator();
794 : #endif
795 :
796 5556 : LOG_SCOPE ("make_coarsening_compatible()", "MeshRefinement");
797 :
798 : // Unless we encounter a specific situation level-one
799 : // will be satisfied after executing this loop just once
800 2778 : bool level_one_satisfied = true;
801 :
802 :
803 : // Unless we encounter a specific situation we will be compatible
804 : // with any selected refinement flags
805 45499 : bool compatible_with_refinement = true;
806 :
807 :
808 : // find the maximum h and p levels in the mesh
809 45499 : unsigned int max_level = 0;
810 45499 : unsigned int max_p_level = 0;
811 :
812 : {
813 : // First we look at all the active level-0 elements. Since it doesn't make
814 : // sense to coarsen them we must un-set their coarsen flags if
815 : // they are set.
816 26162560 : for (auto & elem : _mesh.active_element_ptr_range())
817 : {
818 13984984 : max_level = std::max(max_level, elem->level());
819 13984984 : max_p_level =
820 13984984 : std::max(max_p_level,
821 14650212 : static_cast<unsigned int>(elem->p_level()));
822 :
823 13997692 : if ((elem->level() == 0) &&
824 107667 : (elem->refinement_flag() == Elem::COARSEN))
825 850 : elem->set_refinement_flag(Elem::DO_NOTHING);
826 :
827 15210216 : if ((elem->p_level() == 0) &&
828 1225232 : (elem->p_refinement_flag() == Elem::COARSEN))
829 10 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
830 41251 : }
831 : }
832 :
833 : // Even if there are no refined elements on this processor then
834 : // there may still be work for us to do on e.g. ancestor elements.
835 : // At the very least we need to be in the loop if a distributed mesh
836 : // needs to synchronize data.
837 : #if 0
838 : if (max_level == 0 && max_p_level == 0)
839 : {
840 : // But we still have to check with other processors
841 : this->comm().min(compatible_with_refinement);
842 :
843 : return compatible_with_refinement;
844 : }
845 : #endif
846 :
847 : // Loop over all the active elements. If an element is marked
848 : // for coarsening we better check its neighbors. If ANY of these neighbors
849 : // are marked for refinement AND are at the same level then there is a
850 : // conflict. By convention refinement wins, so we un-mark the element for
851 : // coarsening. Level-one would be violated in this case so we need to re-run
852 : // the loop.
853 45499 : if (_face_level_mismatch_limit)
854 : {
855 :
856 22632 : repeat:
857 1548 : level_one_satisfied = true;
858 :
859 276 : do
860 : {
861 1824 : level_one_satisfied = true;
862 :
863 32803542 : for (auto & elem : _mesh.active_element_ptr_range())
864 : {
865 969192 : bool my_flag_changed = false;
866 :
867 17951694 : if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
868 : // the coarsen flag is set
869 : {
870 3719058 : const unsigned int my_level = elem->level();
871 :
872 15922958 : for (auto n : elem->side_index_range())
873 : {
874 : const Elem * neighbor =
875 12182963 : topological_neighbor(elem, point_locator.get(), n);
876 :
877 12182963 : if (neighbor != nullptr && // I have a
878 11960630 : neighbor != remote_elem) // neighbor here
879 : {
880 1043756 : if (neighbor->active()) // and it is active
881 : {
882 12280828 : if ((neighbor->level() == my_level) &&
883 455326 : (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
884 : // and wants to be refined
885 : {
886 11657 : elem->set_refinement_flag(Elem::DO_NOTHING);
887 740 : my_flag_changed = true;
888 11657 : break;
889 : }
890 : }
891 : else // I have a neighbor and it is not active. That means it has children.
892 : { // While it _may_ be possible to coarsen us if all the children of
893 : // that element want to be coarsened, it is impossible to know at this
894 : // stage. Forget about it for the moment... This can be handled in
895 : // two steps.
896 119730 : elem->set_refinement_flag(Elem::DO_NOTHING);
897 6592 : my_flag_changed = true;
898 119730 : break;
899 : }
900 : }
901 : }
902 : }
903 17951694 : if (elem->p_refinement_flag() == Elem::COARSEN) // If
904 : // the element is active and the order reduction flag is set
905 : {
906 164 : const unsigned int my_p_level = elem->p_level();
907 :
908 7129 : for (auto n : elem->side_index_range())
909 : {
910 : const Elem * neighbor =
911 5896 : topological_neighbor(elem, point_locator.get(), n);
912 :
913 5896 : if (neighbor != nullptr && // I have a
914 4637 : neighbor != remote_elem) // neighbor here
915 : {
916 388 : if (neighbor->active()) // and it is active
917 : {
918 348 : if ((neighbor->p_level() > my_p_level &&
919 0 : neighbor->p_refinement_flag() != Elem::COARSEN)
920 3772 : || (neighbor->p_level() == my_p_level &&
921 116 : neighbor->p_refinement_flag() == Elem::REFINE))
922 : {
923 217 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
924 14 : my_flag_changed = true;
925 14 : break;
926 : }
927 : }
928 : else // I have a neighbor and it is not active.
929 : { // We need to find which of its children
930 : // have me as a neighbor, and maintain
931 : // level one p compatibility with them.
932 : // Because we currently have level one h
933 : // compatibility, we don't need to check
934 : // grandchildren
935 :
936 20 : libmesh_assert(neighbor->has_children());
937 1616 : for (auto & subneighbor : neighbor->child_ref_range())
938 1513 : if (&subneighbor != remote_elem &&
939 2362 : subneighbor.active() &&
940 973 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
941 598 : if ((subneighbor.p_level() > my_p_level &&
942 34 : subneighbor.p_refinement_flag() != Elem::COARSEN)
943 1014 : || (subneighbor.p_level() == my_p_level &&
944 0 : subneighbor.p_refinement_flag() == Elem::REFINE))
945 : {
946 194 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
947 12 : my_flag_changed = true;
948 12 : break;
949 : }
950 239 : if (my_flag_changed)
951 20 : break;
952 : }
953 : }
954 : }
955 : }
956 :
957 : // If the current element's flag changed, we hadn't
958 : // satisfied the level one rule.
959 17219059 : if (my_flag_changed)
960 7340 : level_one_satisfied = false;
961 :
962 : // Additionally, if it has non-local neighbors, and
963 : // we're not in serial, then we'll eventually have to
964 : // return compatible_with_refinement = false, because
965 : // our change has to propagate to neighboring
966 : // processors.
967 1818442 : if (my_flag_changed && !_mesh.is_serial())
968 15959 : for (auto n : elem->side_index_range())
969 : {
970 : Elem * neigh =
971 14890 : topological_neighbor(elem, point_locator.get(), n);
972 :
973 14890 : if (!neigh)
974 2936 : continue;
975 11954 : if (neigh == remote_elem ||
976 10998 : neigh->processor_id() !=
977 0 : this->processor_id())
978 : {
979 5340 : compatible_with_refinement = false;
980 5340 : break;
981 : }
982 : // FIXME - for non-level one meshes we should
983 : // test all descendants
984 0 : if (neigh->has_children())
985 2074 : for (auto & child : neigh->child_ref_range())
986 1732 : if (&child == remote_elem ||
987 1732 : child.processor_id() !=
988 0 : this->processor_id())
989 : {
990 258 : compatible_with_refinement = false;
991 258 : break;
992 : }
993 : }
994 31143 : }
995 : }
996 34233 : while (!level_one_satisfied);
997 :
998 : } // end if (_face_level_mismatch_limit)
999 :
1000 :
1001 : // Next we look at all of the ancestor cells.
1002 : // If there is a parent cell with all of its children
1003 : // wanting to be unrefined then the element is a candidate
1004 : // for unrefinement. If all the children don't
1005 : // all want to be unrefined then ALL of them need to have their
1006 : // unrefinement flags cleared.
1007 246463 : for (int level = max_level; level >= 0; level--)
1008 50097680 : for (auto & elem : as_range(_mesh.level_elements_begin(level), _mesh.level_elements_end(level)))
1009 26385164 : if (elem->ancestor())
1010 : {
1011 : // right now the element hasn't been disqualified
1012 : // as a candidate for unrefinement
1013 460034 : bool is_a_candidate = true;
1014 460034 : bool found_remote_child = false;
1015 :
1016 31986335 : for (auto & child : elem->child_ref_range())
1017 : {
1018 25412440 : if (&child == remote_elem)
1019 40 : found_remote_child = true;
1020 24975490 : else if ((child.refinement_flag() != Elem::COARSEN) ||
1021 109666 : !child.active() )
1022 1759862 : is_a_candidate = false;
1023 : }
1024 :
1025 6291233 : if (!is_a_candidate && !found_remote_child)
1026 : {
1027 5516178 : elem->set_refinement_flag(Elem::INACTIVE);
1028 :
1029 27679850 : for (auto & child : elem->child_ref_range())
1030 : {
1031 22163672 : if (&child == remote_elem)
1032 0 : continue;
1033 22163672 : if (child.refinement_flag() == Elem::COARSEN)
1034 : {
1035 20714 : level_one_satisfied = false;
1036 20714 : child.set_refinement_flag(Elem::DO_NOTHING);
1037 : }
1038 : }
1039 : }
1040 174144 : }
1041 :
1042 51335 : if (!level_one_satisfied && _face_level_mismatch_limit) goto repeat;
1043 :
1044 :
1045 : // If all the children of a parent are set to be coarsened
1046 : // then flag the parent so that they can kill their kids.
1047 :
1048 : // On a distributed mesh, we won't always be able to determine this
1049 : // on parent elements with remote children, even if we own the
1050 : // parent, without communication.
1051 : //
1052 : // We'll first communicate *to* parents' owners when we determine
1053 : // they cannot be coarsened, then we'll sync the final refinement
1054 : // flag *from* the parents.
1055 :
1056 : // uncoarsenable_parents[p] live on processor id p
1057 45499 : const processor_id_type n_proc = _mesh.n_processors();
1058 2778 : const processor_id_type my_proc_id = _mesh.processor_id();
1059 45499 : const bool distributed_mesh = !_mesh.is_replicated();
1060 :
1061 : std::vector<std::vector<dof_id_type>>
1062 45499 : uncoarsenable_parents(n_proc);
1063 :
1064 8847179 : for (auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1065 : {
1066 : // Presume all the children are flagged for coarsening and
1067 : // then look for a contradiction
1068 390434 : bool all_children_flagged_for_coarsening = true;
1069 :
1070 7618968 : for (auto & child : elem->child_ref_range())
1071 : {
1072 7207257 : if (&child != remote_elem &&
1073 425278 : child.refinement_flag() != Elem::COARSEN)
1074 : {
1075 379022 : all_children_flagged_for_coarsening = false;
1076 5122993 : if (!distributed_mesh)
1077 378870 : break;
1078 1184248 : if (child.processor_id() != elem->processor_id())
1079 : {
1080 97195 : uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
1081 0 : break;
1082 : }
1083 : }
1084 : }
1085 :
1086 4562672 : if (all_children_flagged_for_coarsening)
1087 301651 : elem->set_refinement_flag(Elem::COARSEN_INACTIVE);
1088 : else
1089 4358216 : elem->set_refinement_flag(Elem::INACTIVE);
1090 41251 : }
1091 :
1092 : // If we have a distributed mesh, we might need to sync up
1093 : // INACTIVE vs. COARSEN_INACTIVE flags.
1094 45499 : if (distributed_mesh)
1095 : {
1096 : // We'd better still be in sync here
1097 8 : parallel_object_only();
1098 :
1099 : Parallel::MessageTag
1100 29443 : uncoarsenable_tag = this->comm().get_unique_tag();
1101 29431 : std::vector<Parallel::Request> uncoarsenable_push_requests(n_proc-1);
1102 :
1103 342370 : for (processor_id_type p = 0; p != n_proc; ++p)
1104 : {
1105 312939 : if (p == my_proc_id)
1106 29427 : continue;
1107 :
1108 : Parallel::Request &request =
1109 283508 : uncoarsenable_push_requests[p - (p > my_proc_id)];
1110 :
1111 283508 : _mesh.comm().send
1112 283508 : (p, uncoarsenable_parents[p], request, uncoarsenable_tag);
1113 : }
1114 :
1115 312939 : for (processor_id_type p = 1; p != n_proc; ++p)
1116 : {
1117 16 : std::vector<dof_id_type> my_uncoarsenable_parents;
1118 283508 : _mesh.comm().receive
1119 283508 : (Parallel::any_source, my_uncoarsenable_parents,
1120 12 : uncoarsenable_tag);
1121 :
1122 355589 : for (const auto & id : my_uncoarsenable_parents)
1123 : {
1124 72081 : Elem & elem = _mesh.elem_ref(id);
1125 0 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1126 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1127 0 : elem.set_refinement_flag(Elem::INACTIVE);
1128 : }
1129 : }
1130 :
1131 29431 : Parallel::wait(uncoarsenable_push_requests);
1132 :
1133 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1134 29431 : &Elem::set_refinement_flag);
1135 : sync_dofobject_data_by_id
1136 58854 : (this->comm(), _mesh.not_local_elements_begin(),
1137 29443 : _mesh.not_local_elements_end(),
1138 : // We'd like a smaller sync, but this leads to bugs?
1139 : // SyncCoarsenInactive(),
1140 : hsync);
1141 29419 : }
1142 :
1143 : // If one processor finds an incompatibility, we're globally
1144 : // incompatible
1145 45499 : this->comm().min(compatible_with_refinement);
1146 :
1147 90998 : return compatible_with_refinement;
1148 41251 : }
1149 :
1150 :
1151 :
1152 :
1153 :
1154 :
1155 :
1156 :
1157 46489 : bool MeshRefinement::make_refinement_compatible()
1158 : {
1159 : // This function must be run on all processors at once
1160 2824 : parallel_object_only();
1161 :
1162 : // We may need a PointLocator for topological_neighbor() tests
1163 : // later, which we need to make sure gets constructed on all
1164 : // processors at once.
1165 44995 : std::unique_ptr<PointLocatorBase> point_locator;
1166 :
1167 : #ifdef LIBMESH_ENABLE_PERIODIC
1168 : bool has_periodic_boundaries =
1169 46489 : _periodic_boundaries && !_periodic_boundaries->empty();
1170 2824 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
1171 :
1172 2824 : if (has_periodic_boundaries)
1173 532 : point_locator = _mesh.sub_point_locator();
1174 : #endif
1175 :
1176 2824 : LOG_SCOPE ("make_refinement_compatible()", "MeshRefinement");
1177 :
1178 : // Unless we encounter a specific situation we will be compatible
1179 : // with any selected coarsening flags
1180 46489 : bool compatible_with_coarsening = true;
1181 :
1182 : // This loop enforces the level-1 rule. We should only
1183 : // execute it if the user indeed wants level-1 satisfied!
1184 46489 : if (_face_level_mismatch_limit)
1185 : {
1186 : // Unless we encounter a specific situation level-one
1187 : // will be satisfied after executing this loop just once
1188 1324 : bool level_one_satisfied = true;
1189 :
1190 188 : do
1191 : {
1192 1512 : level_one_satisfied = true;
1193 :
1194 22170122 : for (auto & elem : _mesh.active_element_ptr_range())
1195 : {
1196 11673270 : const unsigned short n_sides = elem->n_sides();
1197 :
1198 12158160 : if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1199 : // h refinement flag is set
1200 : {
1201 362485 : const unsigned int my_level = elem->level();
1202 :
1203 1799040 : for (unsigned short side = 0; side != n_sides;
1204 : ++side)
1205 : {
1206 : Elem * neighbor =
1207 1436555 : topological_neighbor(elem, point_locator.get(), side);
1208 :
1209 1436881 : if (neighbor != nullptr && // I have a
1210 1511103 : neighbor != remote_elem && // neighbor here
1211 149110 : neighbor->active()) // and it is active
1212 : {
1213 : // Case 1: The neighbor is at the same level I am.
1214 : // 1a: The neighbor will be refined -> NO PROBLEM
1215 : // 1b: The neighbor won't be refined -> NO PROBLEM
1216 : // 1c: The neighbor wants to be coarsened -> PROBLEM
1217 1020523 : if (neighbor->level() == my_level)
1218 : {
1219 927130 : if (neighbor->refinement_flag() == Elem::COARSEN)
1220 : {
1221 0 : neighbor->set_refinement_flag(Elem::DO_NOTHING);
1222 7 : if (neighbor->parent())
1223 0 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1224 7 : compatible_with_coarsening = false;
1225 0 : level_one_satisfied = false;
1226 : }
1227 : }
1228 :
1229 :
1230 : // Case 2: The neighbor is one level lower than I am.
1231 : // The neighbor thus MUST be refined to satisfy
1232 : // the level-one rule, regardless of whether it
1233 : // was originally flagged for refinement. If it
1234 : // wasn't flagged already we need to repeat
1235 : // this process.
1236 93393 : else if ((neighbor->level()+1) == my_level)
1237 : {
1238 93393 : if (neighbor->refinement_flag() != Elem::REFINE)
1239 : {
1240 564 : neighbor->set_refinement_flag(Elem::REFINE);
1241 11170 : if (neighbor->parent())
1242 418 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1243 10602 : compatible_with_coarsening = false;
1244 564 : level_one_satisfied = false;
1245 : }
1246 : }
1247 : #ifdef DEBUG
1248 : // Note that the only other possibility is that the
1249 : // neighbor is already refined, in which case it isn't
1250 : // active and we should never get here.
1251 : else
1252 0 : libmesh_error_msg("ERROR: Neighbor level must be equal or 1 higher than mine.");
1253 : #endif
1254 : }
1255 : }
1256 : }
1257 12158160 : if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1258 : // p refinement flag is set
1259 : {
1260 2802 : const unsigned int my_p_level = elem->p_level();
1261 :
1262 129865 : for (unsigned int side=0; side != n_sides; side++)
1263 : {
1264 : Elem * neighbor =
1265 109624 : topological_neighbor(elem, point_locator.get(), side);
1266 :
1267 109624 : if (neighbor != nullptr && // I have a
1268 90735 : neighbor != remote_elem) // neighbor here
1269 : {
1270 11914 : if (neighbor->active()) // and it is active
1271 : {
1272 86177 : if (neighbor->p_level() < my_p_level &&
1273 18 : neighbor->p_refinement_flag() != Elem::REFINE)
1274 : {
1275 0 : neighbor->set_p_refinement_flag(Elem::REFINE);
1276 0 : level_one_satisfied = false;
1277 22 : compatible_with_coarsening = false;
1278 : }
1279 91967 : if (neighbor->p_level() == my_p_level &&
1280 5808 : neighbor->p_refinement_flag() == Elem::COARSEN)
1281 : {
1282 0 : neighbor->set_p_refinement_flag(Elem::DO_NOTHING);
1283 0 : level_one_satisfied = false;
1284 0 : compatible_with_coarsening = false;
1285 : }
1286 : }
1287 : else // I have an inactive neighbor
1288 : {
1289 56 : libmesh_assert(neighbor->has_children());
1290 1806 : for (auto & subneighbor : neighbor->child_ref_range())
1291 2203 : if (&subneighbor != remote_elem && subneighbor.active() &&
1292 999 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
1293 : {
1294 708 : if (subneighbor.p_level() < my_p_level &&
1295 44 : subneighbor.p_refinement_flag() != Elem::REFINE)
1296 : {
1297 : // We should already be level one
1298 : // compatible
1299 6 : libmesh_assert_greater (subneighbor.p_level() + 2u,
1300 : my_p_level);
1301 6 : subneighbor.set_p_refinement_flag(Elem::REFINE);
1302 6 : level_one_satisfied = false;
1303 70 : compatible_with_coarsening = false;
1304 : }
1305 676 : if (subneighbor.p_level() == my_p_level &&
1306 12 : subneighbor.p_refinement_flag() == Elem::COARSEN)
1307 : {
1308 0 : subneighbor.set_p_refinement_flag(Elem::DO_NOTHING);
1309 0 : level_one_satisfied = false;
1310 0 : compatible_with_coarsening = false;
1311 : }
1312 : }
1313 : }
1314 : }
1315 : }
1316 : }
1317 25183 : }
1318 : }
1319 :
1320 27627 : while (!level_one_satisfied);
1321 : } // end if (_face_level_mismatch_limit)
1322 :
1323 : // If we're not compatible on one processor, we're globally not
1324 : // compatible
1325 46489 : this->comm().min(compatible_with_coarsening);
1326 :
1327 50807 : return compatible_with_coarsening;
1328 42171 : }
1329 :
1330 :
1331 :
1332 :
1333 38232 : bool MeshRefinement::_coarsen_elements ()
1334 : {
1335 : // This function must be run on all processors at once
1336 1340 : parallel_object_only();
1337 :
1338 1340 : LOG_SCOPE ("_coarsen_elements()", "MeshRefinement");
1339 :
1340 : // Flags indicating if this call actually changes the mesh
1341 38232 : bool mesh_changed = false;
1342 38232 : bool mesh_p_changed = false;
1343 :
1344 : // Clear the unused_elements data structure.
1345 : // The elements have been packed since it was built,
1346 : // so there are _no_ unused elements. We cannot trust
1347 : // any iterators currently in this data structure.
1348 : // _unused_elements.clear();
1349 :
1350 : // Loop over the elements first to determine if the mesh will
1351 : // undergo h-coarsening. If it will, then we'll need to communicate
1352 : // more ghosted elements. We need to communicate them *before* we
1353 : // do the coarsening; otherwise it is possible to coarsen away a
1354 : // one-element-thick layer partition and leave the partitions on
1355 : // either side unable to figure out how to talk to each other.
1356 19605884 : for (auto & elem : _mesh.element_ptr_range())
1357 10938959 : if (elem->refinement_flag() == Elem::COARSEN)
1358 : {
1359 4914 : mesh_changed = true;
1360 4914 : break;
1361 35552 : }
1362 :
1363 : // If the mesh changed on any processor, it changed globally
1364 38232 : this->comm().max(mesh_changed);
1365 :
1366 : // And then we may need to widen the ghosting layers.
1367 38232 : if (mesh_changed)
1368 4992 : MeshCommunication().send_coarse_ghosts(_mesh);
1369 :
1370 30328710 : for (auto & elem : _mesh.element_ptr_range())
1371 : {
1372 : // Make sure we transfer the children's boundary id(s)
1373 : // up to its parent when necessary before coarsening.
1374 16805719 : _mesh.get_boundary_info().transfer_boundary_ids_from_children(elem);
1375 :
1376 : // active elements flagged for coarsening will
1377 : // no longer be deleted until MeshRefinement::contract()
1378 16805719 : if (elem->refinement_flag() == Elem::COARSEN)
1379 : {
1380 : // Huh? no level-0 element should be active
1381 : // and flagged for coarsening.
1382 0 : libmesh_assert_not_equal_to (elem->level(), 0);
1383 :
1384 : // Remove this element from any neighbor
1385 : // lists that point to it.
1386 0 : elem->nullify_neighbors();
1387 :
1388 : // Remove any boundary information associated
1389 : // with this element if we do not allow children to have boundary info.
1390 : // Otherwise, we will do the removal in `transfer_boundary_ids_from_children`
1391 : // to make sure we don't delete the information before it is transferred
1392 0 : if (!_mesh.get_boundary_info().is_children_on_boundary_side())
1393 0 : _mesh.get_boundary_info().remove (elem);
1394 :
1395 : // Add this iterator to the _unused_elements
1396 : // data structure so we might fill it.
1397 : // The _unused_elements optimization is currently off.
1398 : // _unused_elements.push_back (it);
1399 :
1400 : // Don't delete the element until
1401 : // MeshRefinement::contract()
1402 : // _mesh.delete_elem(elem);
1403 : }
1404 :
1405 : // inactive elements flagged for coarsening
1406 : // will become active
1407 15966255 : else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1408 : {
1409 345905 : elem->coarsen();
1410 17078 : libmesh_assert (elem->active());
1411 :
1412 : // the mesh has certainly changed
1413 345905 : mesh_changed = true;
1414 : }
1415 16805719 : if (elem->p_refinement_flag() == Elem::COARSEN)
1416 : {
1417 220 : if (elem->p_level() > 0)
1418 : {
1419 14 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1420 206 : elem->set_p_level(elem->p_level() - 1);
1421 206 : mesh_p_changed = true;
1422 : }
1423 : else
1424 : {
1425 0 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
1426 : }
1427 : }
1428 35552 : }
1429 :
1430 38232 : this->comm().max(mesh_p_changed);
1431 :
1432 : // And we may need to update DistributedMesh values reflecting the changes
1433 38232 : if (mesh_changed)
1434 4992 : _mesh.update_parallel_id_counts();
1435 :
1436 : // Node processor ids may need to change if an element of that id
1437 : // was coarsened away
1438 38232 : if (mesh_changed && !_mesh.is_serial())
1439 : {
1440 : // Update the _new_nodes_map so that processors can
1441 : // find requested nodes
1442 1164 : this->update_nodes_map ();
1443 :
1444 1164 : MeshCommunication().make_nodes_parallel_consistent (_mesh);
1445 :
1446 : // Clear the _new_nodes_map
1447 1164 : this->clear();
1448 :
1449 : #ifdef DEBUG
1450 4 : MeshTools::libmesh_assert_valid_procids<Node>(_mesh);
1451 : #endif
1452 : }
1453 :
1454 : // If p levels changed all we need to do is make sure that parent p
1455 : // levels changed in sync
1456 38232 : if (mesh_p_changed && !_mesh.is_serial())
1457 : {
1458 192 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1459 : }
1460 :
1461 40912 : return (mesh_changed || mesh_p_changed);
1462 : }
1463 :
1464 :
1465 :
1466 53335 : bool MeshRefinement::_refine_elements ()
1467 : {
1468 1854 : libmesh_assert(_mesh.is_prepared() || _mesh.is_replicated());
1469 :
1470 : // This function must be run on all processors at once
1471 1854 : parallel_object_only();
1472 :
1473 : // Update the _new_nodes_map so that elements can
1474 : // find nodes to connect to.
1475 53335 : this->update_nodes_map ();
1476 :
1477 3708 : LOG_SCOPE ("_refine_elements()", "MeshRefinement");
1478 :
1479 : // Iterate over the elements, counting the elements
1480 : // flagged for h refinement.
1481 1854 : dof_id_type n_elems_flagged = 0;
1482 :
1483 29877806 : for (auto & elem : _mesh.element_ptr_range())
1484 16541021 : if (elem->refinement_flag() == Elem::REFINE)
1485 1227156 : n_elems_flagged++;
1486 :
1487 : // Construct a local vector of Elem * which have been
1488 : // previously marked for refinement. We reserve enough
1489 : // space to allow for every element to be refined.
1490 1854 : std::vector<Elem *> local_copy_of_elements;
1491 53335 : local_copy_of_elements.reserve(n_elems_flagged);
1492 :
1493 : // If mesh p levels changed, we might need to synchronize parent p
1494 : // levels on a distributed mesh.
1495 53335 : bool mesh_p_changed = false;
1496 :
1497 : // Iterate over the elements, looking for elements flagged for
1498 : // refinement.
1499 :
1500 : // If we are on a ReplicatedMesh, then we just do the refinement in
1501 : // the same order on every processor and everything stays in sync.
1502 :
1503 : // If we are on a DistributedMesh, that's impossible.
1504 : //
1505 : // If the mesh is distributed, we need to make sure that if we end
1506 : // up as the owner of a new node, which might happen if that node is
1507 : // attached to one of our own elements, then we have given it a
1508 : // legitimate node id and our own processor id. We generate
1509 : // legitimate node ids and use our own processor id when we are
1510 : // refining our own elements but not when we refine others'
1511 : // elements. Therefore we want to refine our own elements *first*,
1512 : // thereby generating all nodes which might belong to us, and then
1513 : // refine others' elements *after*, thereby generating nodes with
1514 : // temporary ids which we know we will discard.
1515 : //
1516 : // Even if the DistributedMesh is serialized, we can't just treat it
1517 : // like a ReplicatedMesh, because DistributedMesh doesn't *trust*
1518 : // users to refine partitioned elements in a serialized way, so it
1519 : // assigns temporary ids, so we need to synchronize ids afterward to
1520 : // be safe anyway, so we might as well use the distributed mesh code
1521 : // path.
1522 19351956 : for (auto & elem : _mesh.is_replicated() ? _mesh.active_element_ptr_range() : _mesh.active_local_element_ptr_range())
1523 : {
1524 10817639 : if (elem->refinement_flag() == Elem::REFINE)
1525 771917 : local_copy_of_elements.push_back(elem);
1526 10818941 : if (elem->p_refinement_flag() == Elem::REFINE &&
1527 2606 : elem->active())
1528 : {
1529 11047 : elem->set_p_level(elem->p_level()+1);
1530 9743 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1531 9743 : mesh_p_changed = true;
1532 : }
1533 49627 : }
1534 :
1535 53335 : if (!_mesh.is_replicated())
1536 : {
1537 75120 : for (auto & elem : as_range(_mesh.active_not_local_elements_begin(),
1538 1609454 : _mesh.active_not_local_elements_end()))
1539 : {
1540 730237 : if (elem->refinement_flag() == Elem::REFINE)
1541 405612 : local_copy_of_elements.push_back(elem);
1542 730237 : if (elem->p_refinement_flag() == Elem::REFINE &&
1543 0 : elem->active())
1544 : {
1545 8163 : elem->set_p_level(elem->p_level()+1);
1546 8163 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1547 8163 : mesh_p_changed = true;
1548 : }
1549 37371 : }
1550 : }
1551 :
1552 : // Now iterate over the local copies and refine each one.
1553 : // This may resize the mesh's internal container and invalidate
1554 : // any existing iterators.
1555 1230864 : for (auto & elem : local_copy_of_elements)
1556 1177529 : elem->refine(*this);
1557 :
1558 : // The mesh changed if there were elements h refined
1559 53335 : bool mesh_changed = !local_copy_of_elements.empty();
1560 :
1561 : // If the mesh changed on any processor, it changed globally
1562 53335 : this->comm().max(mesh_changed);
1563 53335 : this->comm().max(mesh_p_changed);
1564 :
1565 : // And we may need to update DistributedMesh values reflecting the changes
1566 53335 : if (mesh_changed)
1567 31071 : _mesh.update_parallel_id_counts();
1568 :
1569 53335 : if (mesh_changed && !_mesh.is_replicated())
1570 : {
1571 21505 : MeshCommunication().make_elems_parallel_consistent (_mesh);
1572 21505 : MeshCommunication().make_new_nodes_parallel_consistent (_mesh);
1573 : #ifdef DEBUG
1574 38 : _mesh.libmesh_assert_valid_parallel_ids();
1575 : #endif
1576 : }
1577 :
1578 : // If we're refining a ReplicatedMesh, then we haven't yet assigned
1579 : // node processor ids. But if we're refining a partitioned
1580 : // ReplicatedMesh, then we *need* to assign node processor ids.
1581 55451 : if (mesh_changed && _mesh.is_replicated() &&
1582 62901 : (_mesh.unpartitioned_elements_begin() ==
1583 71409 : _mesh.unpartitioned_elements_end()))
1584 9252 : Partitioner::set_node_processor_ids(_mesh);
1585 :
1586 53335 : if (mesh_p_changed && !_mesh.is_replicated())
1587 : {
1588 2368 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1589 : }
1590 :
1591 : // Clear the _new_nodes_map and _unused_elements data structures.
1592 53335 : this->clear();
1593 :
1594 106670 : return (mesh_changed || mesh_p_changed);
1595 : }
1596 :
1597 :
1598 60132 : void MeshRefinement::_smooth_flags(bool refining, bool coarsening)
1599 : {
1600 : // Smoothing can break in weird ways on a mesh with broken topology
1601 : #ifdef DEBUG
1602 2084 : MeshTools::libmesh_assert_valid_neighbors(_mesh);
1603 : #endif
1604 :
1605 : // Repeat until flag changes match on every processor
1606 0 : do
1607 : {
1608 : // Repeat until coarsening & refinement flags jive
1609 2084 : bool satisfied = false;
1610 164 : do
1611 : {
1612 : // If we're refining or coarsening, hit the corresponding
1613 : // face level test code. Short-circuiting || is our friend
1614 : const bool coarsening_satisfied =
1615 111609 : !coarsening ||
1616 44191 : this->make_coarsening_compatible();
1617 :
1618 : const bool refinement_satisfied =
1619 112577 : !refining ||
1620 45159 : this->make_refinement_compatible();
1621 :
1622 : bool smoothing_satisfied =
1623 67418 : !this->eliminate_unrefined_patches();// &&
1624 :
1625 67418 : if (_edge_level_mismatch_limit)
1626 0 : smoothing_satisfied = smoothing_satisfied &&
1627 0 : !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);
1628 :
1629 67418 : if (_node_level_mismatch_limit)
1630 0 : smoothing_satisfied = smoothing_satisfied &&
1631 0 : !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);
1632 :
1633 67418 : if (_overrefined_boundary_limit>=0)
1634 45192 : smoothing_satisfied = smoothing_satisfied &&
1635 21937 : !this->limit_overrefined_boundary(_overrefined_boundary_limit);
1636 :
1637 67418 : if (_underrefined_boundary_limit>=0)
1638 44624 : smoothing_satisfied = smoothing_satisfied &&
1639 21369 : !this->limit_underrefined_boundary(_underrefined_boundary_limit);
1640 :
1641 67418 : satisfied = (coarsening_satisfied &&
1642 69666 : refinement_satisfied &&
1643 : smoothing_satisfied);
1644 :
1645 2248 : libmesh_assert(this->comm().verify(satisfied));
1646 : }
1647 65170 : while (!satisfied);
1648 : }
1649 116518 : while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
1650 60132 : }
1651 :
1652 :
1653 0 : void MeshRefinement::uniformly_p_refine (unsigned int n)
1654 : {
1655 : // Refine n times
1656 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1657 0 : for (auto & elem : _mesh.active_element_ptr_range())
1658 : {
1659 : // P refine all the active elements
1660 0 : elem->set_p_level(elem->p_level()+1);
1661 0 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1662 0 : }
1663 0 : }
1664 :
1665 :
1666 :
1667 0 : void MeshRefinement::uniformly_p_coarsen (unsigned int n)
1668 : {
1669 : // Coarsen p times
1670 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1671 0 : for (auto & elem : _mesh.active_element_ptr_range())
1672 0 : if (elem->p_level() > 0)
1673 : {
1674 : // P coarsen all the active elements
1675 0 : elem->set_p_level(elem->p_level()-1);
1676 0 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1677 0 : }
1678 0 : }
1679 :
1680 :
1681 :
1682 13372 : void MeshRefinement::uniformly_refine (unsigned int n)
1683 : {
1684 : // Refine n times
1685 : // FIXME - this won't work if n>1 and the mesh
1686 : // has already been attached to an equation system
1687 28763 : for (unsigned int rstep=0; rstep<n; rstep++)
1688 : {
1689 : // Clean up the refinement flags
1690 15391 : this->clean_refinement_flags();
1691 :
1692 : // Flag all the active elements for refinement.
1693 1888224 : for (auto & elem : _mesh.active_element_ptr_range())
1694 1003696 : elem->set_refinement_flag(Elem::REFINE);
1695 :
1696 : // Refine all the elements we just flagged.
1697 15391 : this->_refine_elements();
1698 : }
1699 :
1700 : // Finally, the new mesh probably needs to be prepared for use
1701 13372 : if (n > 0)
1702 12167 : _mesh.prepare_for_use ();
1703 13372 : }
1704 :
1705 :
1706 :
1707 1132 : void MeshRefinement::uniformly_coarsen (unsigned int n)
1708 : {
1709 : // Coarsen n times
1710 2264 : for (unsigned int rstep=0; rstep<n; rstep++)
1711 : {
1712 : // Clean up the refinement flags
1713 1132 : this->clean_refinement_flags();
1714 :
1715 : // Flag all the active elements for coarsening.
1716 471746 : for (auto & elem : _mesh.active_element_ptr_range())
1717 : {
1718 260509 : elem->set_refinement_flag(Elem::COARSEN);
1719 286261 : if (elem->parent())
1720 25752 : elem->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
1721 1068 : }
1722 :
1723 : // On a distributed mesh, we may have parent elements with
1724 : // remote active children. To keep flags consistent, we'll need
1725 : // a communication step.
1726 1132 : if (!_mesh.is_replicated())
1727 : {
1728 908 : const processor_id_type n_proc = _mesh.n_processors();
1729 4 : const processor_id_type my_proc_id = _mesh.processor_id();
1730 :
1731 : std::vector<std::vector<dof_id_type>>
1732 916 : parents_to_coarsen(n_proc);
1733 :
1734 117848 : for (const auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1735 57878 : if (elem->processor_id() != my_proc_id &&
1736 48 : elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1737 14280 : parents_to_coarsen[elem->processor_id()].push_back(elem->id());
1738 :
1739 : Parallel::MessageTag
1740 916 : coarsen_tag = this->comm().get_unique_tag();
1741 908 : std::vector<Parallel::Request> coarsen_push_requests(n_proc-1);
1742 :
1743 10452 : for (processor_id_type p = 0; p != n_proc; ++p)
1744 : {
1745 9544 : if (p == my_proc_id)
1746 904 : continue;
1747 :
1748 : Parallel::Request &request =
1749 8636 : coarsen_push_requests[p - (p > my_proc_id)];
1750 :
1751 8636 : _mesh.comm().send
1752 8636 : (p, parents_to_coarsen[p], request, coarsen_tag);
1753 : }
1754 :
1755 9544 : for (processor_id_type p = 1; p != n_proc; ++p)
1756 : {
1757 8 : std::vector<dof_id_type> my_parents_to_coarsen;
1758 8636 : _mesh.comm().receive
1759 8636 : (Parallel::any_source, my_parents_to_coarsen,
1760 8 : coarsen_tag);
1761 :
1762 21980 : for (const auto & id : my_parents_to_coarsen)
1763 : {
1764 13344 : Elem & elem = _mesh.elem_ref(id);
1765 36 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1766 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1767 36 : elem.set_refinement_flag(Elem::COARSEN_INACTIVE);
1768 : }
1769 : }
1770 :
1771 908 : Parallel::wait(coarsen_push_requests);
1772 :
1773 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1774 908 : &Elem::set_refinement_flag);
1775 : sync_dofobject_data_by_id
1776 1812 : (this->comm(), _mesh.not_local_elements_begin(),
1777 916 : _mesh.not_local_elements_end(),
1778 : // We'd like a smaller sync, but this leads to bugs?
1779 : // SyncCoarsenInactive(),
1780 : hsync);
1781 900 : }
1782 :
1783 : // Coarsen all the elements we just flagged.
1784 1132 : this->_coarsen_elements();
1785 : }
1786 :
1787 :
1788 : // Finally, the new mesh probably needs to be prepared for use
1789 1132 : if (n > 0)
1790 1132 : _mesh.prepare_for_use ();
1791 1132 : }
1792 :
1793 :
1794 :
1795 14722100 : Elem * MeshRefinement::topological_neighbor(Elem * elem,
1796 : const PointLocatorBase * point_locator,
1797 : const unsigned int side) const
1798 : {
1799 : #ifdef LIBMESH_ENABLE_PERIODIC
1800 14722100 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1801 : {
1802 488472 : libmesh_assert(point_locator);
1803 786435 : return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1804 : }
1805 : #endif
1806 14353293 : return elem->neighbor_ptr(side);
1807 : }
1808 :
1809 :
1810 :
1811 1972 : bool MeshRefinement::has_topological_neighbor(const Elem * elem,
1812 : const PointLocatorBase * point_locator,
1813 : const Elem * neighbor) const
1814 : {
1815 : #ifdef LIBMESH_ENABLE_PERIODIC
1816 1972 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1817 : {
1818 0 : libmesh_assert(point_locator);
1819 0 : return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1820 : }
1821 : #endif
1822 1800 : return elem->has_neighbor(neighbor);
1823 : }
1824 :
1825 :
1826 :
1827 : } // namespace libMesh
1828 :
1829 :
1830 : #endif
|
