libMesh
mesh_refinement.C
Go to the documentation of this file.
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
95  ParallelObject(m),
96  _mesh(m),
97  _use_member_parameters(false),
98  _coarsen_by_parents(false),
99  _refine_fraction(0.3),
100  _coarsen_fraction(0.0),
101  _max_h_level(libMesh::invalid_uint),
102  _coarsen_threshold(10),
103  _nelem_target(0),
104  _absolute_global_tolerance(0.0),
105  _face_level_mismatch_limit(1),
106  _edge_level_mismatch_limit(0),
107  _node_level_mismatch_limit(0),
108  _overrefined_boundary_limit(0),
109  _underrefined_boundary_limit(0),
110  _allow_unrefined_patches(false),
111  _enforce_mismatch_limit_prior_to_refinement(false)
112 #ifdef LIBMESH_ENABLE_PERIODIC
113  , _periodic_boundaries(nullptr)
114 #endif
115 {
116 }
117 
118 
119 
120 #ifdef LIBMESH_ENABLE_PERIODIC
122 {
123  _periodic_boundaries = pb_ptr;
124 }
125 #endif
126 
127 
128 
130 
131 
132 
134 {
136 }
137 
138 
139 
141  unsigned int child,
142  unsigned int node,
143  processor_id_type proc_id)
144 {
145  LOG_SCOPE("add_node()", "MeshRefinement");
146 
147  unsigned int parent_n = parent.as_parent_node(child, node);
148 
149  if (parent_n != libMesh::invalid_uint)
150  return parent.node_ptr(parent_n);
151 
152  const std::vector<std::pair<dof_id_type, dof_id_type>>
153  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  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  if (const auto new_node_id = _new_nodes_map.find(bracketing_nodes);
165  new_node_id != DofObject::invalid_id)
166  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  Point p; // defaults to 0,0,0
173 
174  for (auto n : parent.node_index_range())
175  {
176  // The value from the embedding matrix
177  const Real em_val = parent.embedding_matrix(child,node,n);
178 
179  if (em_val != 0.)
180  {
181  p.add_scaled (parent.point(n), em_val);
182 
183  // If we'd already found the node we shouldn't be here
184  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  Node * new_node = _mesh.add_point (p, DofObject::invalid_id, proc_id);
192 
193  libmesh_assert(new_node);
194 
195  // But then we'll make sure this node is marked as unpartitioned.
197 
198  // Add the node to the map.
199  _new_nodes_map.add_node(*new_node, bracketing_nodes);
200 
201  // Return the address of the new node
202  return new_node;
203 }
204 
205 
206 
208 {
209  libmesh_assert(elem);
210  return _mesh.add_elem (elem);
211 }
212 
213 Elem * MeshRefinement::add_elem (std::unique_ptr<Elem> elem)
214 {
215  libmesh_assert(elem);
216  return _mesh.add_elem(std::move(elem));
217 }
218 
219 
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  parallel_object_only();
227 
228  // Make sure the input error vector is valid
229 #ifdef DEBUG
230  for (const auto & val : error_per_cell)
231  {
232  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  std::vector<ErrorVectorReal> & epc = error_per_parent;
242  libmesh_assert(this->comm().verify(epc));
243 #endif // #ifdef DEBUG
244 
245  // error values on uncoarsenable elements will be left at -1
246  error_per_parent.clear();
247  error_per_parent.resize(error_per_cell.size(), 0.0);
248 
249  {
250  // Find which elements are uncoarsenable
251  for (auto & elem : _mesh.active_local_element_ptr_range())
252  {
253  Elem * parent = elem->parent();
254 
255  // Active elements are uncoarsenable
256  error_per_parent[elem->id()] = -1.0;
257 
258  // Grandparents and up are uncoarsenable
259  while (parent)
260  {
261  parent = parent->parent();
262  if (parent)
263  {
264  const dof_id_type parentid = parent->id();
265  libmesh_assert_less (parentid, error_per_parent.size());
266  error_per_parent[parentid] = -1.0;
267  }
268  }
269  }
270 
271  // Sync between processors.
272  // Use a reference to std::vector to avoid confusing
273  // this->comm().min
274  std::vector<ErrorVectorReal> & epp = error_per_parent;
275  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  for (auto & elem : _mesh.active_local_element_ptr_range())
287  {
288  Elem * parent = elem->parent();
289 
290  // Calculate each contribution to parent cells
291  if (parent)
292  {
293  const dof_id_type parentid = parent->id();
294  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  if (error_per_parent[parentid] != -1.0)
300  error_per_parent[parentid] += (error_per_cell[elem->id()] *
301  error_per_cell[elem->id()]);
302  }
303  }
304 
305  // Sum the vector across all processors
306  this->comm().sum(static_cast<std::vector<ErrorVectorReal> &>(error_per_parent));
307 
308  // Calculate the min and max as we loop
309  parent_error_min = std::numeric_limits<double>::max();
310  parent_error_max = 0.;
311 
312  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  if (error_per_parent[i] < 0.)
320  {
321  error_per_parent[i] = -1.;
322  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  if (error_per_cell[i])
329  {
330  error_per_parent[i] = error_per_cell[i];
331  continue;
332  }
333  // if not, then e_parent = sqrt(sum(e_child^2))
334  else
335  error_per_parent[i] = std::sqrt(error_per_parent[i]);
336 
337  parent_error_min = std::min (parent_error_min,
338  error_per_parent[i]);
339  parent_error_max = std::max (parent_error_max,
340  error_per_parent[i]);
341  }
342 }
343 
344 
345 
347 {
348  this->_new_nodes_map.init(_mesh);
349 }
350 
351 
352 
353 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  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  std::unique_ptr<PointLocatorBase> point_locator;
362 
363 #ifdef LIBMESH_ENABLE_PERIODIC
364  bool has_periodic_boundaries =
366  libmesh_assert(this->comm().verify(has_periodic_boundaries));
367 
368  if (has_periodic_boundaries)
369  point_locator = _mesh.sub_point_locator();
370 #endif
371 
372  bool failure = false;
373 
374 #ifndef NDEBUG
375  Elem * failed_elem = nullptr;
376  Elem * failed_neighbor = nullptr;
377 #endif // !NDEBUG
378 
379  for (auto & elem : _mesh.active_local_element_ptr_range())
380  for (auto n : elem->side_index_range())
381  {
382  Elem * neighbor =
383  topological_neighbor(elem, point_locator.get(), n);
384 
385  if (!neighbor || !neighbor->active() ||
386  neighbor == remote_elem)
387  continue;
388 
389  if ((neighbor->level() + 1 < elem->level()) ||
390  (neighbor->p_level() + 1 < elem->p_level()) ||
391  (neighbor->p_level() > elem->p_level() + 1))
392  {
393  failure = true;
394 #ifndef NDEBUG
395  failed_elem = elem;
396  failed_neighbor = neighbor;
397 #endif // !NDEBUG
398  break;
399  }
400  }
401 
402  // If any processor failed, we failed globally
403  this->comm().max(failure);
404 
405  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  if (libmesh_assert_pass)
411  {
412  libMesh::out << "MeshRefinement Level one failure, element: "
413  << *failed_elem
414  << std::endl;
415  libMesh::out << "MeshRefinement Level one failure, neighbor: "
416  << *failed_neighbor
417  << std::endl;
418  }
419 #endif // !NDEBUG
420  libmesh_assert(!libmesh_assert_pass);
421  return false;
422  }
423  return true;
424 }
425 
426 
427 
428 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  parallel_object_only();
432 
433  bool found_flag = false;
434 
435 #ifndef NDEBUG
436  Elem * failed_elem = nullptr;
437 #endif
438 
439  // Search for local flags
440  for (auto & elem : _mesh.active_local_element_ptr_range())
441  if (elem->refinement_flag() == Elem::REFINE ||
442  elem->refinement_flag() == Elem::COARSEN ||
443  elem->p_refinement_flag() == Elem::REFINE ||
444  elem->p_refinement_flag() == Elem::COARSEN)
445  {
446  found_flag = true;
447 #ifndef NDEBUG
448  failed_elem = elem;
449 #endif
450  break;
451  }
452 
453  // If we found a flag on any processor, it counts
454  this->comm().max(found_flag);
455 
456  if (found_flag)
457  {
458 #ifndef NDEBUG
459  if (libmesh_assert_pass)
460  {
461  libMesh::out <<
462  "MeshRefinement test_unflagged failure, element: " <<
463  *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  libmesh_assert(!libmesh_assert_pass);
469  return false;
470  }
471  return true;
472 }
473 
474 
475 
477 {
478  // This function must be run on all processors at once
479  parallel_object_only();
480 
481  // We can't yet turn a non-level-one mesh into a level-one mesh
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  bool elements_flagged = false;
489 
490  for (auto & elem : _mesh.element_ptr_range())
491  {
492  // This might be left over from the last step
493  const Elem::RefinementState flag = elem->refinement_flag();
494 
495  // Set refinement flag to INACTIVE if the
496  // element isn't active
497  if ( !elem->active())
498  {
499  elem->set_refinement_flag(Elem::INACTIVE);
500  elem->set_p_refinement_flag(Elem::INACTIVE);
501  }
502  else if (flag == Elem::JUST_REFINED)
503  elem->set_refinement_flag(Elem::DO_NOTHING);
504  else if (!elements_flagged)
505  {
506  if (flag == Elem::REFINE || flag == Elem::COARSEN)
507  elements_flagged = true;
508  else
509  {
510  const Elem::RefinementState pflag =
511  elem->p_refinement_flag();
512  if (pflag == Elem::REFINE || pflag == Elem::COARSEN)
513  elements_flagged = true;
514  }
515  }
516  }
517 
518  // Did *any* processor find elements flagged for AMR/C?
519  _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  if (!elements_flagged)
524  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  bool flags_were_consistent = this->make_flags_parallel_consistent();
531 
532  libmesh_assert (flags_were_consistent);
533 #endif
534 
535  // Smooth refinement and coarsening flags
536  _smooth_flags(true, true);
537 
538  // First coarsen the flagged elements.
539  const bool coarsening_changed_mesh =
540  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  this->_refine_elements();
564 
565  // First coarsen the flagged elements.
566  // Finally, the new mesh needs to be prepared for use
567  if (coarsening_changed_mesh || refining_changed_mesh)
568  {
569 #ifdef DEBUG
571 #endif
572 
574 
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  return true;
585  }
586  else
587  {
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  return false;
599 
600 }
601 
602 
603 
604 
605 
606 
607 
609 {
610  // This function must be run on all processors at once
611  parallel_object_only();
612 
613  // We can't yet turn a non-level-one mesh into a level-one mesh
616 
617  // Possibly clean up the refinement flags from
618  // a previous step
619  for (auto & elem : _mesh.element_ptr_range())
620  {
621  // Set refinement flag to INACTIVE if the
622  // element isn't active
623  if (!elem->active())
624  {
625  elem->set_refinement_flag(Elem::INACTIVE);
626  elem->set_p_refinement_flag(Elem::INACTIVE);
627  }
628 
629  // This might be left over from the last step
630  if (elem->refinement_flag() == Elem::JUST_REFINED)
631  elem->set_refinement_flag(Elem::DO_NOTHING);
632  }
633 
634  // Parallel consistency has to come first, or coarsening
635  // along processor boundaries might occasionally be falsely
636  // prevented
637  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  libmesh_assert(flags_were_consistent);
646 
647  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  _smooth_flags(false, true);
653 
654  // Coarsen the flagged elements.
655  const bool mesh_changed =
656  this->_coarsen_elements ();
657 
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  if (mesh_changed)
672 
673  return mesh_changed;
674 }
675 
676 
677 
678 
679 
680 
681 
683 {
684  // This function must be run on all processors at once
685  parallel_object_only();
686 
689 
690  // Possibly clean up the refinement flags from
691  // a previous step
692  for (auto & elem : _mesh.element_ptr_range())
693  {
694  // Set refinement flag to INACTIVE if the
695  // element isn't active
696  if (!elem->active())
697  {
698  elem->set_refinement_flag(Elem::INACTIVE);
699  elem->set_p_refinement_flag(Elem::INACTIVE);
700  }
701 
702  // This might be left over from the last step
703  if (elem->refinement_flag() == Elem::JUST_REFINED)
704  elem->set_refinement_flag(Elem::DO_NOTHING);
705  }
706 
707  // Parallel consistency has to come first, or coarsening
708  // along processor boundaries might occasionally be falsely
709  // prevented
710  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  libmesh_assert(flags_were_consistent);
719 
720  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  _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  this->_refine_elements();
731 
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  if (mesh_changed)
743 
744  return mesh_changed;
745 }
746 
747 
748 
750 {
751  // This function must be run on all processors at once
752  parallel_object_only();
753 
754  LOG_SCOPE ("make_flags_parallel_consistent()", "MeshRefinement");
755 
759  (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
760 
764  (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  bool parallel_consistent = hsync.parallel_consistent &&
769  psync.parallel_consistent;
770  this->comm().min(parallel_consistent);
771 
772  return parallel_consistent;
773 }
774 
775 
776 
778 {
779  // This function must be run on all processors at once
780  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  std::unique_ptr<PointLocatorBase> point_locator;
786 
787 #ifdef LIBMESH_ENABLE_PERIODIC
788  bool has_periodic_boundaries =
790  libmesh_assert(this->comm().verify(has_periodic_boundaries));
791 
792  if (has_periodic_boundaries)
793  point_locator = _mesh.sub_point_locator();
794 #endif
795 
796  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  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  bool compatible_with_refinement = true;
806 
807 
808  // find the maximum h and p levels in the mesh
809  unsigned int max_level = 0;
810  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  for (auto & elem : _mesh.active_element_ptr_range())
817  {
818  max_level = std::max(max_level, elem->level());
819  max_p_level =
820  std::max(max_p_level,
821  static_cast<unsigned int>(elem->p_level()));
822 
823  if ((elem->level() == 0) &&
824  (elem->refinement_flag() == Elem::COARSEN))
825  elem->set_refinement_flag(Elem::DO_NOTHING);
826 
827  if ((elem->p_level() == 0) &&
828  (elem->p_refinement_flag() == Elem::COARSEN))
829  elem->set_p_refinement_flag(Elem::DO_NOTHING);
830  }
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.
854  {
855 
856  repeat:
857  level_one_satisfied = true;
858 
859  do
860  {
861  level_one_satisfied = true;
862 
863  for (auto & elem : _mesh.active_element_ptr_range())
864  {
865  bool my_flag_changed = false;
866 
867  if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
868  // the coarsen flag is set
869  {
870  const unsigned int my_level = elem->level();
871 
872  for (auto n : elem->side_index_range())
873  {
874  const Elem * neighbor =
875  topological_neighbor(elem, point_locator.get(), n);
876 
877  if (neighbor != nullptr && // I have a
878  neighbor != remote_elem) // neighbor here
879  {
880  if (neighbor->active()) // and it is active
881  {
882  if ((neighbor->level() == my_level) &&
883  (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
884  // and wants to be refined
885  {
887  my_flag_changed = true;
888  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  elem->set_refinement_flag(Elem::DO_NOTHING);
897  my_flag_changed = true;
898  break;
899  }
900  }
901  }
902  }
903  if (elem->p_refinement_flag() == Elem::COARSEN) // If
904  // the element is active and the order reduction flag is set
905  {
906  const unsigned int my_p_level = elem->p_level();
907 
908  for (auto n : elem->side_index_range())
909  {
910  const Elem * neighbor =
911  topological_neighbor(elem, point_locator.get(), n);
912 
913  if (neighbor != nullptr && // I have a
914  neighbor != remote_elem) // neighbor here
915  {
916  if (neighbor->active()) // and it is active
917  {
918  if ((neighbor->p_level() > my_p_level &&
919  neighbor->p_refinement_flag() != Elem::COARSEN)
920  || (neighbor->p_level() == my_p_level &&
921  neighbor->p_refinement_flag() == Elem::REFINE))
922  {
924  my_flag_changed = true;
925  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  libmesh_assert(neighbor->has_children());
937  for (auto & subneighbor : neighbor->child_ref_range())
938  if (&subneighbor != remote_elem &&
939  subneighbor.active() &&
940  has_topological_neighbor(&subneighbor, point_locator.get(), elem))
941  if ((subneighbor.p_level() > my_p_level &&
942  subneighbor.p_refinement_flag() != Elem::COARSEN)
943  || (subneighbor.p_level() == my_p_level &&
944  subneighbor.p_refinement_flag() == Elem::REFINE))
945  {
946  elem->set_p_refinement_flag(Elem::DO_NOTHING);
947  my_flag_changed = true;
948  break;
949  }
950  if (my_flag_changed)
951  break;
952  }
953  }
954  }
955  }
956 
957  // If the current element's flag changed, we hadn't
958  // satisfied the level one rule.
959  if (my_flag_changed)
960  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  if (my_flag_changed && !_mesh.is_serial())
968  for (auto n : elem->side_index_range())
969  {
970  Elem * neigh =
971  topological_neighbor(elem, point_locator.get(), n);
972 
973  if (!neigh)
974  continue;
975  if (neigh == remote_elem ||
976  neigh->processor_id() !=
977  this->processor_id())
978  {
979  compatible_with_refinement = false;
980  break;
981  }
982  // FIXME - for non-level one meshes we should
983  // test all descendants
984  if (neigh->has_children())
985  for (auto & child : neigh->child_ref_range())
986  if (&child == remote_elem ||
987  child.processor_id() !=
988  this->processor_id())
989  {
990  compatible_with_refinement = false;
991  break;
992  }
993  }
994  }
995  }
996  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  for (int level = max_level; level >= 0; level--)
1008  for (auto & elem : as_range(_mesh.level_elements_begin(level), _mesh.level_elements_end(level)))
1009  if (elem->ancestor())
1010  {
1011  // right now the element hasn't been disqualified
1012  // as a candidate for unrefinement
1013  bool is_a_candidate = true;
1014  bool found_remote_child = false;
1015 
1016  for (auto & child : elem->child_ref_range())
1017  {
1018  if (&child == remote_elem)
1019  found_remote_child = true;
1020  else if ((child.refinement_flag() != Elem::COARSEN) ||
1021  !child.active() )
1022  is_a_candidate = false;
1023  }
1024 
1025  if (!is_a_candidate && !found_remote_child)
1026  {
1028 
1029  for (auto & child : elem->child_ref_range())
1030  {
1031  if (&child == remote_elem)
1032  continue;
1033  if (child.refinement_flag() == Elem::COARSEN)
1034  {
1035  level_one_satisfied = false;
1036  child.set_refinement_flag(Elem::DO_NOTHING);
1037  }
1038  }
1039  }
1040  }
1041 
1042  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  const processor_id_type n_proc = _mesh.n_processors();
1058  const processor_id_type my_proc_id = _mesh.processor_id();
1059  const bool distributed_mesh = !_mesh.is_replicated();
1060 
1061  std::vector<std::vector<dof_id_type>>
1062  uncoarsenable_parents(n_proc);
1063 
1064  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  bool all_children_flagged_for_coarsening = true;
1069 
1070  for (auto & child : elem->child_ref_range())
1071  {
1072  if (&child != remote_elem &&
1073  child.refinement_flag() != Elem::COARSEN)
1074  {
1075  all_children_flagged_for_coarsening = false;
1076  if (!distributed_mesh)
1077  break;
1078  if (child.processor_id() != elem->processor_id())
1079  {
1080  uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
1081  break;
1082  }
1083  }
1084  }
1085 
1086  if (all_children_flagged_for_coarsening)
1087  elem->set_refinement_flag(Elem::COARSEN_INACTIVE);
1088  else
1089  elem->set_refinement_flag(Elem::INACTIVE);
1090  }
1091 
1092  // If we have a distributed mesh, we might need to sync up
1093  // INACTIVE vs. COARSEN_INACTIVE flags.
1094  if (distributed_mesh)
1095  {
1096  // We'd better still be in sync here
1097  parallel_object_only();
1098 
1100  uncoarsenable_tag = this->comm().get_unique_tag();
1101  std::vector<Parallel::Request> uncoarsenable_push_requests(n_proc-1);
1102 
1103  for (processor_id_type p = 0; p != n_proc; ++p)
1104  {
1105  if (p == my_proc_id)
1106  continue;
1107 
1109  uncoarsenable_push_requests[p - (p > my_proc_id)];
1110 
1111  _mesh.comm().send
1112  (p, uncoarsenable_parents[p], request, uncoarsenable_tag);
1113  }
1114 
1115  for (processor_id_type p = 1; p != n_proc; ++p)
1116  {
1117  std::vector<dof_id_type> my_uncoarsenable_parents;
1118  _mesh.comm().receive
1119  (Parallel::any_source, my_uncoarsenable_parents,
1120  uncoarsenable_tag);
1121 
1122  for (const auto & id : my_uncoarsenable_parents)
1123  {
1124  Elem & elem = _mesh.elem_ref(id);
1128  }
1129  }
1130 
1131  Parallel::wait(uncoarsenable_push_requests);
1132 
1136  (this->comm(), _mesh.not_local_elements_begin(),
1137  _mesh.not_local_elements_end(),
1138  // We'd like a smaller sync, but this leads to bugs?
1139  // SyncCoarsenInactive(),
1140  hsync);
1141  }
1142 
1143  // If one processor finds an incompatibility, we're globally
1144  // incompatible
1145  this->comm().min(compatible_with_refinement);
1146 
1147  return compatible_with_refinement;
1148 }
1149 
1150 
1151 
1152 
1153 
1154 
1155 
1156 
1158 {
1159  // This function must be run on all processors at once
1160  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  std::unique_ptr<PointLocatorBase> point_locator;
1166 
1167 #ifdef LIBMESH_ENABLE_PERIODIC
1168  bool has_periodic_boundaries =
1170  libmesh_assert(this->comm().verify(has_periodic_boundaries));
1171 
1172  if (has_periodic_boundaries)
1173  point_locator = _mesh.sub_point_locator();
1174 #endif
1175 
1176  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  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!
1185  {
1186  // Unless we encounter a specific situation level-one
1187  // will be satisfied after executing this loop just once
1188  bool level_one_satisfied = true;
1189 
1190  do
1191  {
1192  level_one_satisfied = true;
1193 
1194  for (auto & elem : _mesh.active_element_ptr_range())
1195  {
1196  const unsigned short n_sides = elem->n_sides();
1197 
1198  if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1199  // h refinement flag is set
1200  {
1201  const unsigned int my_level = elem->level();
1202 
1203  for (unsigned short side = 0; side != n_sides;
1204  ++side)
1205  {
1206  Elem * neighbor =
1207  topological_neighbor(elem, point_locator.get(), side);
1208 
1209  if (neighbor != nullptr && // I have a
1210  neighbor != remote_elem && // neighbor here
1211  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  if (neighbor->level() == my_level)
1218  {
1219  if (neighbor->refinement_flag() == Elem::COARSEN)
1220  {
1222  if (neighbor->parent())
1224  compatible_with_coarsening = false;
1225  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  else if ((neighbor->level()+1) == my_level)
1237  {
1238  if (neighbor->refinement_flag() != Elem::REFINE)
1239  {
1240  neighbor->set_refinement_flag(Elem::REFINE);
1241  if (neighbor->parent())
1243  compatible_with_coarsening = false;
1244  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  libmesh_error_msg("ERROR: Neighbor level must be equal or 1 higher than mine.");
1253 #endif
1254  }
1255  }
1256  }
1257  if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1258  // p refinement flag is set
1259  {
1260  const unsigned int my_p_level = elem->p_level();
1261 
1262  for (unsigned int side=0; side != n_sides; side++)
1263  {
1264  Elem * neighbor =
1265  topological_neighbor(elem, point_locator.get(), side);
1266 
1267  if (neighbor != nullptr && // I have a
1268  neighbor != remote_elem) // neighbor here
1269  {
1270  if (neighbor->active()) // and it is active
1271  {
1272  if (neighbor->p_level() < my_p_level &&
1273  neighbor->p_refinement_flag() != Elem::REFINE)
1274  {
1276  level_one_satisfied = false;
1277  compatible_with_coarsening = false;
1278  }
1279  if (neighbor->p_level() == my_p_level &&
1280  neighbor->p_refinement_flag() == Elem::COARSEN)
1281  {
1283  level_one_satisfied = false;
1284  compatible_with_coarsening = false;
1285  }
1286  }
1287  else // I have an inactive neighbor
1288  {
1289  libmesh_assert(neighbor->has_children());
1290  for (auto & subneighbor : neighbor->child_ref_range())
1291  if (&subneighbor != remote_elem && subneighbor.active() &&
1292  has_topological_neighbor(&subneighbor, point_locator.get(), elem))
1293  {
1294  if (subneighbor.p_level() < my_p_level &&
1295  subneighbor.p_refinement_flag() != Elem::REFINE)
1296  {
1297  // We should already be level one
1298  // compatible
1299  libmesh_assert_greater (subneighbor.p_level() + 2u,
1300  my_p_level);
1301  subneighbor.set_p_refinement_flag(Elem::REFINE);
1302  level_one_satisfied = false;
1303  compatible_with_coarsening = false;
1304  }
1305  if (subneighbor.p_level() == my_p_level &&
1306  subneighbor.p_refinement_flag() == Elem::COARSEN)
1307  {
1308  subneighbor.set_p_refinement_flag(Elem::DO_NOTHING);
1309  level_one_satisfied = false;
1310  compatible_with_coarsening = false;
1311  }
1312  }
1313  }
1314  }
1315  }
1316  }
1317  }
1318  }
1319 
1320  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  this->comm().min(compatible_with_coarsening);
1326 
1327  return compatible_with_coarsening;
1328 }
1329 
1330 
1331 
1332 
1334 {
1335  // This function must be run on all processors at once
1336  parallel_object_only();
1337 
1338  LOG_SCOPE ("_coarsen_elements()", "MeshRefinement");
1339 
1340  // Flags indicating if this call actually changes the mesh
1341  bool mesh_changed = false;
1342  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  for (auto & elem : _mesh.element_ptr_range())
1357  if (elem->refinement_flag() == Elem::COARSEN)
1358  {
1359  mesh_changed = true;
1360  break;
1361  }
1362 
1363  // If the mesh changed on any processor, it changed globally
1364  this->comm().max(mesh_changed);
1365 
1366  // And then we may need to widen the ghosting layers.
1367  if (mesh_changed)
1369 
1370  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.
1375 
1376  // active elements flagged for coarsening will
1377  // no longer be deleted until MeshRefinement::contract()
1378  if (elem->refinement_flag() == Elem::COARSEN)
1379  {
1380  // Huh? no level-0 element should be active
1381  // and flagged for coarsening.
1382  libmesh_assert_not_equal_to (elem->level(), 0);
1383 
1384  // Remove this element from any neighbor
1385  // lists that point to it.
1386  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
1393  _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  else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1408  {
1409  elem->coarsen();
1410  libmesh_assert (elem->active());
1411 
1412  // the mesh has certainly changed
1413  mesh_changed = true;
1414  }
1415  if (elem->p_refinement_flag() == Elem::COARSEN)
1416  {
1417  if (elem->p_level() > 0)
1418  {
1419  elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1420  elem->set_p_level(elem->p_level() - 1);
1421  mesh_p_changed = true;
1422  }
1423  else
1424  {
1425  elem->set_p_refinement_flag(Elem::DO_NOTHING);
1426  }
1427  }
1428  }
1429 
1430  this->comm().max(mesh_p_changed);
1431 
1432  // And we may need to update DistributedMesh values reflecting the changes
1433  if (mesh_changed)
1435 
1436  // Node processor ids may need to change if an element of that id
1437  // was coarsened away
1438  if (mesh_changed && !_mesh.is_serial())
1439  {
1440  // Update the _new_nodes_map so that processors can
1441  // find requested nodes
1442  this->update_nodes_map ();
1443 
1445 
1446  // Clear the _new_nodes_map
1447  this->clear();
1448 
1449 #ifdef DEBUG
1450  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  if (mesh_p_changed && !_mesh.is_serial())
1457  {
1459  }
1460 
1461  return (mesh_changed || mesh_p_changed);
1462 }
1463 
1464 
1465 
1467 {
1469 
1470  // This function must be run on all processors at once
1471  parallel_object_only();
1472 
1473  // Update the _new_nodes_map so that elements can
1474  // find nodes to connect to.
1475  this->update_nodes_map ();
1476 
1477  LOG_SCOPE ("_refine_elements()", "MeshRefinement");
1478 
1479  // Iterate over the elements, counting the elements
1480  // flagged for h refinement.
1481  dof_id_type n_elems_flagged = 0;
1482 
1483  for (auto & elem : _mesh.element_ptr_range())
1484  if (elem->refinement_flag() == Elem::REFINE)
1485  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  std::vector<Elem *> local_copy_of_elements;
1491  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  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  for (auto & elem : _mesh.is_replicated() ? _mesh.active_element_ptr_range() : _mesh.active_local_element_ptr_range())
1523  {
1524  if (elem->refinement_flag() == Elem::REFINE)
1525  local_copy_of_elements.push_back(elem);
1526  if (elem->p_refinement_flag() == Elem::REFINE &&
1527  elem->active())
1528  {
1529  elem->set_p_level(elem->p_level()+1);
1530  elem->set_p_refinement_flag(Elem::JUST_REFINED);
1531  mesh_p_changed = true;
1532  }
1533  }
1534 
1535  if (!_mesh.is_replicated())
1536  {
1537  for (auto & elem : as_range(_mesh.active_not_local_elements_begin(),
1538  _mesh.active_not_local_elements_end()))
1539  {
1540  if (elem->refinement_flag() == Elem::REFINE)
1541  local_copy_of_elements.push_back(elem);
1542  if (elem->p_refinement_flag() == Elem::REFINE &&
1543  elem->active())
1544  {
1545  elem->set_p_level(elem->p_level()+1);
1546  elem->set_p_refinement_flag(Elem::JUST_REFINED);
1547  mesh_p_changed = true;
1548  }
1549  }
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  for (auto & elem : local_copy_of_elements)
1556  elem->refine(*this);
1557 
1558  // The mesh changed if there were elements h refined
1559  bool mesh_changed = !local_copy_of_elements.empty();
1560 
1561  // If the mesh changed on any processor, it changed globally
1562  this->comm().max(mesh_changed);
1563  this->comm().max(mesh_p_changed);
1564 
1565  // And we may need to update DistributedMesh values reflecting the changes
1566  if (mesh_changed)
1568 
1569  if (mesh_changed && !_mesh.is_replicated())
1570  {
1573 #ifdef DEBUG
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  if (mesh_changed && _mesh.is_replicated() &&
1582  (_mesh.unpartitioned_elements_begin() ==
1583  _mesh.unpartitioned_elements_end()))
1585 
1586  if (mesh_p_changed && !_mesh.is_replicated())
1587  {
1589  }
1590 
1591  // Clear the _new_nodes_map and _unused_elements data structures.
1592  this->clear();
1593 
1594  return (mesh_changed || mesh_p_changed);
1595 }
1596 
1597 
1598 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
1603 #endif
1604 
1605  // Repeat until flag changes match on every processor
1606  do
1607  {
1608  // Repeat until coarsening & refinement flags jive
1609  bool satisfied = false;
1610  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  !coarsening ||
1617 
1618  const bool refinement_satisfied =
1619  !refining ||
1621 
1622  bool smoothing_satisfied =
1623  !this->eliminate_unrefined_patches();// &&
1624 
1626  smoothing_satisfied = smoothing_satisfied &&
1628 
1630  smoothing_satisfied = smoothing_satisfied &&
1632 
1634  smoothing_satisfied = smoothing_satisfied &&
1636 
1638  smoothing_satisfied = smoothing_satisfied &&
1640 
1641  satisfied = (coarsening_satisfied &&
1642  refinement_satisfied &&
1643  smoothing_satisfied);
1644 
1645  libmesh_assert(this->comm().verify(satisfied));
1646  }
1647  while (!satisfied);
1648  }
1649  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
1650 }
1651 
1652 
1654 {
1655  // Refine n times
1656  for (unsigned int rstep=0; rstep<n; rstep++)
1657  for (auto & elem : _mesh.active_element_ptr_range())
1658  {
1659  // P refine all the active elements
1660  elem->set_p_level(elem->p_level()+1);
1661  elem->set_p_refinement_flag(Elem::JUST_REFINED);
1662  }
1663 }
1664 
1665 
1666 
1668 {
1669  // Coarsen p times
1670  for (unsigned int rstep=0; rstep<n; rstep++)
1671  for (auto & elem : _mesh.active_element_ptr_range())
1672  if (elem->p_level() > 0)
1673  {
1674  // P coarsen all the active elements
1675  elem->set_p_level(elem->p_level()-1);
1676  elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1677  }
1678 }
1679 
1680 
1681 
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  for (unsigned int rstep=0; rstep<n; rstep++)
1688  {
1689  // Clean up the refinement flags
1690  this->clean_refinement_flags();
1691 
1692  // Flag all the active elements for refinement.
1693  for (auto & elem : _mesh.active_element_ptr_range())
1694  elem->set_refinement_flag(Elem::REFINE);
1695 
1696  // Refine all the elements we just flagged.
1697  this->_refine_elements();
1698  }
1699 
1700  // Finally, the new mesh probably needs to be prepared for use
1701  if (n > 0)
1703 }
1704 
1705 
1706 
1708 {
1709  // Coarsen n times
1710  for (unsigned int rstep=0; rstep<n; rstep++)
1711  {
1712  // Clean up the refinement flags
1713  this->clean_refinement_flags();
1714 
1715  // Flag all the active elements for coarsening.
1716  for (auto & elem : _mesh.active_element_ptr_range())
1717  {
1718  elem->set_refinement_flag(Elem::COARSEN);
1719  if (elem->parent())
1720  elem->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
1721  }
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  if (!_mesh.is_replicated())
1727  {
1728  const processor_id_type n_proc = _mesh.n_processors();
1729  const processor_id_type my_proc_id = _mesh.processor_id();
1730 
1731  std::vector<std::vector<dof_id_type>>
1732  parents_to_coarsen(n_proc);
1733 
1734  for (const auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1735  if (elem->processor_id() != my_proc_id &&
1736  elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1737  parents_to_coarsen[elem->processor_id()].push_back(elem->id());
1738 
1740  coarsen_tag = this->comm().get_unique_tag();
1741  std::vector<Parallel::Request> coarsen_push_requests(n_proc-1);
1742 
1743  for (processor_id_type p = 0; p != n_proc; ++p)
1744  {
1745  if (p == my_proc_id)
1746  continue;
1747 
1749  coarsen_push_requests[p - (p > my_proc_id)];
1750 
1751  _mesh.comm().send
1752  (p, parents_to_coarsen[p], request, coarsen_tag);
1753  }
1754 
1755  for (processor_id_type p = 1; p != n_proc; ++p)
1756  {
1757  std::vector<dof_id_type> my_parents_to_coarsen;
1758  _mesh.comm().receive
1759  (Parallel::any_source, my_parents_to_coarsen,
1760  coarsen_tag);
1761 
1762  for (const auto & id : my_parents_to_coarsen)
1763  {
1764  Elem & elem = _mesh.elem_ref(id);
1768  }
1769  }
1770 
1771  Parallel::wait(coarsen_push_requests);
1772 
1776  (this->comm(), _mesh.not_local_elements_begin(),
1777  _mesh.not_local_elements_end(),
1778  // We'd like a smaller sync, but this leads to bugs?
1779  // SyncCoarsenInactive(),
1780  hsync);
1781  }
1782 
1783  // Coarsen all the elements we just flagged.
1784  this->_coarsen_elements();
1785  }
1786 
1787 
1788  // Finally, the new mesh probably needs to be prepared for use
1789  if (n > 0)
1791 }
1792 
1793 
1794 
1796  const PointLocatorBase * point_locator,
1797  const unsigned int side) const
1798 {
1799 #ifdef LIBMESH_ENABLE_PERIODIC
1800  if (_periodic_boundaries && !_periodic_boundaries->empty())
1801  {
1802  libmesh_assert(point_locator);
1803  return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1804  }
1805 #endif
1806  return elem->neighbor_ptr(side);
1807 }
1808 
1809 
1810 
1812  const PointLocatorBase * point_locator,
1813  const Elem * neighbor) const
1814 {
1815 #ifdef LIBMESH_ENABLE_PERIODIC
1816  if (_periodic_boundaries && !_periodic_boundaries->empty())
1817  {
1818  libmesh_assert(point_locator);
1819  return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1820  }
1821 #endif
1822  return elem->has_neighbor(neighbor);
1823 }
1824 
1825 
1826 
1827 } // namespace libMesh
1828 
1829 
1830 #endif
bool limit_level_mismatch_at_edge(const unsigned int max_mismatch)
~MeshRefinement()
Destructor.
bool has_neighbor(const Elem *elem) const
Definition: elem.h:2628
RefinementState refinement_flag() const
Definition: elem.h:3210
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
MPI_Request request
const Elem * parent() const
Definition: elem.h:3030
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
bool is_prepared() const
Definition: mesh_base.h:198
virtual const std::vector< std::pair< dof_id_type, dof_id_type > > bracketing_nodes(unsigned int c, unsigned int n) const
Definition: elem.C:2657
bool _refine_elements()
Refines user-requested elements.
A Node is like a Point, but with more information.
Definition: node.h:52
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
bool test_level_one(bool libmesh_assert_yes=false) const
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1638
bool limit_level_mismatch_at_node(const unsigned int max_mismatch)
This algorithm restricts the maximum level mismatch at any node in the mesh.
void update_nodes_map()
Updates the _new_nodes_map.
const Elem * topological_neighbor(const unsigned int i, const MeshBase &mesh, const PointLocatorBase &point_locator, const PeriodicBoundaries *pb) const
Definition: elem.C:1313
TopologyMap _new_nodes_map
Data structure that holds the new nodes information.
bool limit_underrefined_boundary(const signed char max_mismatch)
bool has_topological_neighbor(const Elem *elem, const PointLocatorBase *point_locator, const Elem *neighbor) const
Local dispatch function for checking the correct has_neighbor function from the Elem class...
MeshBase & _mesh
Reference to the mesh.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
bool refine_elements()
Only refines the user-requested elements.
void make_elems_parallel_consistent(MeshBase &)
Copy ids of ghost elements from their local processors.
static void set_node_processor_ids(MeshBase &mesh)
This function is called after partitioning to set the processor IDs for the nodes.
Definition: partitioner.C:852
dof_id_type find(dof_id_type bracket_node1, dof_id_type bracket_node2) const
Definition: topology_map.C:115
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
RefinementState p_refinement_flag() const
Definition: elem.h:3226
void add_node(const Node &mid_node, const std::vector< std::pair< dof_id_type, dof_id_type >> &bracketing_nodes)
Add a node to the map, between each pair of specified bracketing nodes.
Definition: topology_map.C:53
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
void sum(T &r) const
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
RefinementState
Enumeration of possible element refinement states.
Definition: elem.h:1452
void set_refinement_flag(const RefinementState rflag)
Sets the value of the refinement flag for the element.
Definition: elem.h:3218
void init(MeshBase &)
Definition: topology_map.C:36
const Parallel::Communicator & comm() const
bool limit_overrefined_boundary(const signed char max_mismatch)
unsigned int p_level() const
Definition: elem.h:3108
void clean_refinement_flags()
Sets the refinement flag to Elem::DO_NOTHING for each element in the mesh.
void create_parent_error_vector(const ErrorVector &error_per_cell, ErrorVector &error_per_parent, Real &parent_error_min, Real &parent_error_max)
Calculates the error on all coarsenable parents.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
PeriodicBoundaries * _periodic_boundaries
signed char _underrefined_boundary_limit
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:2346
bool make_refinement_compatible()
Take user-specified refinement flags and augment them so that level-one dependency is satisfied...
bool coarsen_elements()
Only coarsens the user-requested elements.
bool has_topological_neighbor(const Elem *elem, const MeshBase &mesh, const PointLocatorBase &point_locator, const PeriodicBoundaries *pb) const
Definition: elem.C:1350
virtual Real embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
This is the MeshCommunication class.
dof_id_type id() const
Definition: dof_object.h:828
void min(const T &r, T &o, Request &req) const
unsigned char _face_level_mismatch_limit
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:493
virtual unsigned int as_parent_node(unsigned int c, unsigned int n) const
Definition: elem.C:2386
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual void update_parallel_id_counts()=0
Updates parallel caches so that methods like n_elem() accurately reflect changes on other processors...
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
void make_p_levels_parallel_consistent(MeshBase &)
Copy p levels of ghost elements from their local processors.
libmesh_assert(ctx)
bool test_unflagged(bool libmesh_assert_yes=false) const
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
bool make_flags_parallel_consistent()
Copy refinement flags on ghost elements from their local processors.
This is the base class for point locators.
signed char _overrefined_boundary_limit
bool make_coarsening_compatible()
Take user-specified coarsening flags and augment them so that level-one dependency is satisfied...
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
An object whose state is distributed along a set of processors.
MeshRefinement(MeshBase &mesh)
Constructor.
Node * add_node(Elem &parent, unsigned int child, unsigned int node, processor_id_type proc_id)
Add a node to the mesh.
void send_coarse_ghosts(MeshBase &) const
Examine a just-coarsened mesh, and for any newly-coarsened elements, send the associated ghosted elem...
unsigned char _node_level_mismatch_limit
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
unsigned int level() const
Definition: elem.h:3074
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void max(const T &r, T &o, Request &req) const
Elem * add_elem(Elem *elem)
Adds the element elem to the mesh.
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
OStreamProxy out
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
virtual bool is_replicated() const
Definition: mesh_base.h:233
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
void clear()
Deletes all the data that are currently stored.
Elem * topological_neighbor(Elem *elem, const PointLocatorBase *point_locator, const unsigned int side) const
Local dispatch function for getting the correct topological neighbor from the Elem class...
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
void _smooth_flags(bool refining, bool coarsening)
Smooths refinement flags according to current settings.
virtual void libmesh_assert_valid_parallel_ids() const
Verify id and processor_id consistency of our elements and nodes containers.
Definition: mesh_base.h:1516
bool is_children_on_boundary_side() const
unsigned char _edge_level_mismatch_limit
void set_p_refinement_flag(const RefinementState pflag)
Sets the value of the p-refinement flag for the element.
Definition: elem.h:3234
void libmesh_assert_valid_neighbors(const MeshBase &mesh, bool assert_valid_remote_elems=true)
A function for verifying that neighbor connectivity is correct (each element is a neighbor of or desc...
Definition: mesh_tools.C:2205
bool _coarsen_elements()
Coarsens user-requested elements.
virtual const Node * node_ptr(const dof_id_type i) const =0
bool eliminate_unrefined_patches()
This algorithm selects an element for refinement if all of its neighbors are (or will be) refined...
processor_id_type processor_id() const
void make_new_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on new nodes from their local processors.
bool active() const
Definition: elem.h:2941
processor_id_type processor_id() const
Definition: dof_object.h:905
void set_periodic_boundaries_ptr(PeriodicBoundaries *pb_ptr)
Sets the PeriodicBoundaries pointer.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2453
bool has_children() const
Definition: elem.h:2979
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
uint8_t dof_id_type
Definition: id_types.h:67
void transfer_boundary_ids_from_children(const Elem *const parent)
Update parent&#39;s boundary id list so that this information is consistent with its children.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
const RemoteElem * remote_elem
Definition: remote_elem.C:57