libMesh
mesh_refinement.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/remote_elem.h"
39 #include "libmesh/sync_refinement_flags.h"
40 
41 #ifdef DEBUG
42 // Some extra validation for DistributedMesh
43 #include "libmesh/mesh_tools.h"
44 #endif // DEBUG
45 
46 #ifdef LIBMESH_ENABLE_PERIODIC
47 #include "libmesh/periodic_boundaries.h"
48 #endif
49 
50 
51 
52 // Anonymous namespace for helper functions
53 namespace {
54 
55 using namespace libMesh;
56 
57 struct SyncCoarsenInactive
58 {
59  bool operator() (const Elem * elem) const {
60  // If we're not an ancestor, there's no chance our coarsening
61  // settings need to be changed.
62  if (!elem->ancestor())
63  return false;
64 
65  // If we don't have any remote children, we already know enough to
66  // determine the correct refinement flag ourselves.
67  //
68  // If we know we have a child that isn't being coarsened, that
69  // also forces a specific flag.
70  //
71  // Either way there's nothing we need to communicate.
72  bool found_remote_child = false;
73  for (auto & child : elem->child_ref_range())
74  {
75  if (child.refinement_flag() != Elem::COARSEN)
76  return false;
77  if (&child == remote_elem)
78  found_remote_child = true;
79  }
80  return found_remote_child;
81  }
82 };
83 }
84 
85 
86 
87 namespace libMesh
88 {
89 
90 //-----------------------------------------------------------------
91 // Mesh refinement methods
93  ParallelObject(m),
94  _mesh(m),
95  _use_member_parameters(false),
96  _coarsen_by_parents(false),
97  _refine_fraction(0.3),
98  _coarsen_fraction(0.0),
99  _max_h_level(libMesh::invalid_uint),
100  _coarsen_threshold(10),
101  _nelem_target(0),
102  _absolute_global_tolerance(0.0),
103  _face_level_mismatch_limit(1),
104  _edge_level_mismatch_limit(0),
105  _node_level_mismatch_limit(0),
106  _overrefined_boundary_limit(0),
107  _underrefined_boundary_limit(0),
108  _enforce_mismatch_limit_prior_to_refinement(false)
109 #ifdef LIBMESH_ENABLE_PERIODIC
110  , _periodic_boundaries(libmesh_nullptr)
111 #endif
112 {
113 }
114 
115 
116 
117 #ifdef LIBMESH_ENABLE_PERIODIC
119 {
120  _periodic_boundaries = pb_ptr;
121 }
122 #endif
123 
124 
125 
127 {
128  this->clear();
129 }
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  const dof_id_type new_node_id =
160  _new_nodes_map.find(bracketing_nodes);
161 
162  // Return the node if it already exists, but first update the
163  // processor_id if the node is now going to be attached to the
164  // element of a processor which may take ownership of it.
165  if (new_node_id != DofObject::invalid_id)
166  {
167  Node * node = _mesh.node_ptr(new_node_id);
168  if (proc_id < node->processor_id())
169  node->processor_id() = proc_id;
170  return node;
171  }
172 
173  // Otherwise we need to add a new node, with a default id and the
174  // requested processor_id. Figure out where to add the point:
175 
176  Point p; // defaults to 0,0,0
177 
178  for (auto n : parent.node_index_range())
179  {
180  // The value from the embedding matrix
181  const float em_val = parent.embedding_matrix(child,node,n);
182 
183  if (em_val != 0.)
184  {
185  p.add_scaled (parent.point(n), em_val);
186 
187  // If we'd already found the node we shouldn't be here
188  libmesh_assert_not_equal_to (em_val, 1);
189  }
190  }
191 
192  Node * new_node = _mesh.add_point (p, DofObject::invalid_id, proc_id);
193 
194  libmesh_assert(new_node);
195 
196  // Add the node to the map.
197  _new_nodes_map.add_node(*new_node, bracketing_nodes);
198 
199  // Return the address of the new node
200  return new_node;
201 }
202 
203 
204 
206 {
207  libmesh_assert(elem);
208 
209 
210  // // If the unused_elements has any iterators from
211  // // old elements, take the first one
212  // if (!_unused_elements.empty())
213  // {
214  // std::vector<Elem *>::iterator it = _unused_elements.front();
215 
216  // *it = elem;
217 
218  // _unused_elements.pop_front();
219  // }
220 
221  // // Otherwise, use the conventional add method
222  // else
223  // {
224  // _mesh.add_elem (elem);
225  // }
226 
227  // The _unused_elements optimization has been turned off.
228  _mesh.add_elem (elem);
229 
230  return elem;
231 }
232 
233 
234 
236  ErrorVector & error_per_parent,
237  Real & parent_error_min,
238  Real & parent_error_max)
239 {
240  // This function must be run on all processors at once
241  parallel_object_only();
242 
243  // Make sure the input error vector is valid
244 #ifdef DEBUG
245  for (std::size_t i=0; i != error_per_cell.size(); ++i)
246  {
247  libmesh_assert_greater_equal (error_per_cell[i], 0);
248  // isnan() isn't standard C++ yet
249 #ifdef isnan
250  libmesh_assert(!isnan(error_per_cell[i]));
251 #endif
252  }
253 
254  // Use a reference to std::vector to avoid confusing
255  // this->comm().verify
256  std::vector<ErrorVectorReal> & epc = error_per_parent;
257  libmesh_assert(this->comm().verify(epc));
258 #endif // #ifdef DEBUG
259 
260  // error values on uncoarsenable elements will be left at -1
261  error_per_parent.clear();
262  error_per_parent.resize(error_per_cell.size(), 0.0);
263 
264  {
265  // Find which elements are uncoarsenable
266  for (auto & elem : _mesh.active_local_element_ptr_range())
267  {
268  Elem * parent = elem->parent();
269 
270  // Active elements are uncoarsenable
271  error_per_parent[elem->id()] = -1.0;
272 
273  // Grandparents and up are uncoarsenable
274  while (parent)
275  {
276  parent = parent->parent();
277  if (parent)
278  {
279  const dof_id_type parentid = parent->id();
280  libmesh_assert_less (parentid, error_per_parent.size());
281  error_per_parent[parentid] = -1.0;
282  }
283  }
284  }
285 
286  // Sync between processors.
287  // Use a reference to std::vector to avoid confusing
288  // this->comm().min
289  std::vector<ErrorVectorReal> & epp = error_per_parent;
290  this->comm().min(epp);
291  }
292 
293  // The parent's error is defined as the square root of the
294  // sum of the children's errors squared, so errors that are
295  // Hilbert norms remain Hilbert norms.
296  //
297  // Because the children may be on different processors, we
298  // calculate local contributions to the parents' errors squared
299  // first, then sum across processors and take the square roots
300  // second.
301  for (auto & elem : _mesh.active_local_element_ptr_range())
302  {
303  Elem * parent = elem->parent();
304 
305  // Calculate each contribution to parent cells
306  if (parent)
307  {
308  const dof_id_type parentid = parent->id();
309  libmesh_assert_less (parentid, error_per_parent.size());
310 
311  // If the parent has grandchildren we won't be able to
312  // coarsen it, so forget it. Otherwise, add this child's
313  // contribution to the sum of the squared child errors
314  if (error_per_parent[parentid] != -1.0)
315  error_per_parent[parentid] += (error_per_cell[elem->id()] *
316  error_per_cell[elem->id()]);
317  }
318  }
319 
320  // Sum the vector across all processors
321  this->comm().sum(static_cast<std::vector<ErrorVectorReal> &>(error_per_parent));
322 
323  // Calculate the min and max as we loop
324  parent_error_min = std::numeric_limits<double>::max();
325  parent_error_max = 0.;
326 
327  for (std::size_t i = 0; i != error_per_parent.size(); ++i)
328  {
329  // If this element isn't a coarsenable parent with error, we
330  // have nothing to do. Just flag it as -1 and move on
331  // Note that this->comm().sum might have left uncoarsenable
332  // elements with error_per_parent=-n_proc, so reset it to
333  // error_per_parent=-1
334  if (error_per_parent[i] < 0.)
335  {
336  error_per_parent[i] = -1.;
337  continue;
338  }
339 
340  // The error estimator might have already given us an
341  // estimate on the coarsenable parent elements; if so then
342  // we want to retain that estimate
343  if (error_per_cell[i])
344  {
345  error_per_parent[i] = error_per_cell[i];
346  continue;
347  }
348  // if not, then e_parent = sqrt(sum(e_child^2))
349  else
350  error_per_parent[i] = std::sqrt(error_per_parent[i]);
351 
352  parent_error_min = std::min (parent_error_min,
353  error_per_parent[i]);
354  parent_error_max = std::max (parent_error_max,
355  error_per_parent[i]);
356  }
357 }
358 
359 
360 
362 {
363  this->_new_nodes_map.init(_mesh);
364 }
365 
366 
367 
368 bool MeshRefinement::test_level_one (bool libmesh_dbg_var(libmesh_assert_pass))
369 {
370  // This function must be run on all processors at once
371  parallel_object_only();
372 
373  // We may need a PointLocator for topological_neighbor() tests
374  // later, which we need to make sure gets constructed on all
375  // processors at once.
376  UniquePtr<PointLocatorBase> point_locator;
377 
378 #ifdef LIBMESH_ENABLE_PERIODIC
379  bool has_periodic_boundaries =
381  libmesh_assert(this->comm().verify(has_periodic_boundaries));
382 
383  if (has_periodic_boundaries)
384  point_locator = _mesh.sub_point_locator();
385 #endif
386 
387  bool failure = false;
388 
389 #ifndef NDEBUG
390  Elem * failed_elem = libmesh_nullptr;
391  Elem * failed_neighbor = libmesh_nullptr;
392 #endif // !NDEBUG
393 
394  for (auto & elem : _mesh.active_local_element_ptr_range())
395  for (auto n : elem->side_index_range())
396  {
397  Elem * neighbor =
398  topological_neighbor(elem, point_locator.get(), n);
399 
400  if (!neighbor || !neighbor->active() ||
401  neighbor == remote_elem)
402  continue;
403 
404  if ((neighbor->level() + 1 < elem->level()) ||
405  (neighbor->p_level() + 1 < elem->p_level()) ||
406  (neighbor->p_level() > elem->p_level() + 1))
407  {
408  failure = true;
409 #ifndef NDEBUG
410  failed_elem = elem;
411  failed_neighbor = neighbor;
412 #endif // !NDEBUG
413  break;
414  }
415  }
416 
417  // If any processor failed, we failed globally
418  this->comm().max(failure);
419 
420  if (failure)
421  {
422  // We didn't pass the level one test, so libmesh_assert that
423  // we're allowed not to
424 #ifndef NDEBUG
425  if (libmesh_assert_pass)
426  {
427  libMesh::out << "MeshRefinement Level one failure, element: "
428  << *failed_elem
429  << std::endl;
430  libMesh::out << "MeshRefinement Level one failure, neighbor: "
431  << *failed_neighbor
432  << std::endl;
433  }
434 #endif // !NDEBUG
435  libmesh_assert(!libmesh_assert_pass);
436  return false;
437  }
438  return true;
439 }
440 
441 
442 
443 bool MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass))
444 {
445  // This function must be run on all processors at once
446  parallel_object_only();
447 
448  bool found_flag = false;
449 
450 #ifndef NDEBUG
451  Elem * failed_elem = libmesh_nullptr;
452 #endif
453 
454  // Search for local flags
455  for (auto & elem : _mesh.active_local_element_ptr_range())
456  if (elem->refinement_flag() == Elem::REFINE ||
457  elem->refinement_flag() == Elem::COARSEN ||
458  elem->p_refinement_flag() == Elem::REFINE ||
459  elem->p_refinement_flag() == Elem::COARSEN)
460  {
461  found_flag = true;
462 #ifndef NDEBUG
463  failed_elem = elem;
464 #endif
465  break;
466  }
467 
468  // If we found a flag on any processor, it counts
469  this->comm().max(found_flag);
470 
471  if (found_flag)
472  {
473 #ifndef NDEBUG
474  if (libmesh_assert_pass)
475  {
476  libMesh::out <<
477  "MeshRefinement test_unflagged failure, element: " <<
478  *failed_elem << std::endl;
479  }
480 #endif
481  // We didn't pass the "elements are unflagged" test,
482  // so libmesh_assert that we're allowed not to
483  libmesh_assert(!libmesh_assert_pass);
484  return false;
485  }
486  return true;
487 }
488 
489 
490 
492 {
493  // This function must be run on all processors at once
494  parallel_object_only();
495 
496  // We can't yet turn a non-level-one mesh into a level-one mesh
499 
500  // Possibly clean up the refinement flags from
501  // a previous step. While we're at it, see if this method should be
502  // a no-op.
503  bool elements_flagged = false;
504 
505  for (auto & elem : _mesh.element_ptr_range())
506  {
507  // This might be left over from the last step
508  const Elem::RefinementState flag = elem->refinement_flag();
509 
510  // Set refinement flag to INACTIVE if the
511  // element isn't active
512  if ( !elem->active())
513  {
516  }
517  else if (flag == Elem::JUST_REFINED)
519  else if (!elements_flagged)
520  {
521  if (flag == Elem::REFINE || flag == Elem::COARSEN)
522  elements_flagged = true;
523  else
524  {
525  const Elem::RefinementState pflag =
526  elem->p_refinement_flag();
527  if (pflag == Elem::REFINE || pflag == Elem::COARSEN)
528  elements_flagged = true;
529  }
530  }
531  }
532 
533  // Did *any* processor find elements flagged for AMR/C?
534  _mesh.comm().max(elements_flagged);
535 
536  // If we have nothing to do, let's not bother verifying that nothing
537  // is compatible with nothing.
538  if (!elements_flagged)
539  return false;
540 
541  // Parallel consistency has to come first, or coarsening
542  // along processor boundaries might occasionally be falsely
543  // prevented
544 #ifdef DEBUG
545  bool flags_were_consistent = this->make_flags_parallel_consistent();
546 
547  libmesh_assert (flags_were_consistent);
548 #endif
549 
550  // Smooth refinement and coarsening flags
551  _smooth_flags(true, true);
552 
553  // First coarsen the flagged elements.
554  const bool coarsening_changed_mesh =
555  this->_coarsen_elements ();
556 
557  // First coarsen the flagged elements.
558  // FIXME: test_level_one now tests consistency across periodic
559  // boundaries, which requires a point_locator, which just got
560  // invalidated by _coarsen_elements() and hasn't yet been cleared by
561  // prepare_for_use().
562 
563  // libmesh_assert(this->make_coarsening_compatible());
564  // libmesh_assert(this->make_refinement_compatible());
565 
566  // FIXME: This won't pass unless we add a redundant find_neighbors()
567  // call or replace find_neighbors() with on-the-fly neighbor updating
568  // libmesh_assert(!this->eliminate_unrefined_patches());
569 
570  // We can't contract the mesh ourselves anymore - a System might
571  // need to restrict old coefficient vectors first
572  // _mesh.contract();
573 
574  // First coarsen the flagged elements.
575  // Now refine the flagged elements. This will
576  // take up some space, maybe more than what was freed.
577  const bool refining_changed_mesh =
578  this->_refine_elements();
579 
580  // First coarsen the flagged elements.
581  // Finally, the new mesh needs to be prepared for use
582  if (coarsening_changed_mesh || refining_changed_mesh)
583  {
584 #ifdef DEBUG
586 #endif
587 
588  _mesh.prepare_for_use (/*skip_renumber =*/false);
589 
595  // FIXME: This won't pass unless we add a redundant find_neighbors()
596  // call or replace find_neighbors() with on-the-fly neighbor updating
597  // libmesh_assert(!this->eliminate_unrefined_patches());
598 
599  return true;
600  }
601  else
602  {
608  }
609 
610  // Otherwise there was no change in the mesh,
611  // let the user know. Also, there is no need
612  // to prepare the mesh for use since it did not change.
613  return false;
614 
615 }
616 
617 
618 
619 
620 
621 
622 
624 {
625  // This function must be run on all processors at once
626  parallel_object_only();
627 
628  // We can't yet turn a non-level-one mesh into a level-one mesh
631 
632  // Possibly clean up the refinement flags from
633  // a previous step
634  for (auto & elem : _mesh.element_ptr_range())
635  {
636  // Set refinement flag to INACTIVE if the
637  // element isn't active
638  if (!elem->active())
639  {
642  }
643 
644  // This might be left over from the last step
645  if (elem->refinement_flag() == Elem::JUST_REFINED)
647  }
648 
649  // Parallel consistency has to come first, or coarsening
650  // along processor boundaries might occasionally be falsely
651  // prevented
652  bool flags_were_consistent = this->make_flags_parallel_consistent();
653 
654  // In theory, we should be able to remove the above call, which can
655  // be expensive and should be unnecessary. In practice, doing
656  // consistent flagging in parallel is hard, it's impossible to
657  // verify at the library level if it's being done by user code, and
658  // we don't want to abort large parallel runs in opt mode... but we
659  // do want to warn that they should be fixed.
660  libmesh_assert(flags_were_consistent);
661  if (!flags_were_consistent)
662  {
663  libMesh::out << "Refinement flags were not consistent between processors!\n"
664  << "Correcting and continuing.";
665  }
666 
667  // Smooth coarsening flags
668  _smooth_flags(false, true);
669 
670  // Coarsen the flagged elements.
671  const bool mesh_changed =
672  this->_coarsen_elements ();
673 
677  // FIXME: This won't pass unless we add a redundant find_neighbors()
678  // call or replace find_neighbors() with on-the-fly neighbor updating
679  // libmesh_assert(!this->eliminate_unrefined_patches());
680 
681  // We can't contract the mesh ourselves anymore - a System might
682  // need to restrict old coefficient vectors first
683  // _mesh.contract();
684 
685  // Finally, the new mesh may need to be prepared for use
686  if (mesh_changed)
687  _mesh.prepare_for_use (/*skip_renumber =*/false);
688 
689  return mesh_changed;
690 }
691 
692 
693 
694 
695 
696 
697 
699 {
700  // This function must be run on all processors at once
701  parallel_object_only();
702 
705 
706  // Possibly clean up the refinement flags from
707  // a previous step
708  for (auto & elem : _mesh.element_ptr_range())
709  {
710  // Set refinement flag to INACTIVE if the
711  // element isn't active
712  if (!elem->active())
713  {
716  }
717 
718  // This might be left over from the last step
719  if (elem->refinement_flag() == Elem::JUST_REFINED)
721  }
722 
723 
724 
725  // Parallel consistency has to come first, or coarsening
726  // along processor boundaries might occasionally be falsely
727  // prevented
728  bool flags_were_consistent = this->make_flags_parallel_consistent();
729 
730  // In theory, we should be able to remove the above call, which can
731  // be expensive and should be unnecessary. In practice, doing
732  // consistent flagging in parallel is hard, it's impossible to
733  // verify at the library level if it's being done by user code, and
734  // we don't want to abort large parallel runs in opt mode... but we
735  // do want to warn that they should be fixed.
736  libmesh_assert(flags_were_consistent);
737  if (!flags_were_consistent)
738  {
739  libMesh::out << "Refinement flags were not consistent between processors!\n"
740  << "Correcting and continuing.";
741  }
742 
743  // Smooth refinement flags
744  _smooth_flags(true, false);
745 
746  // Now refine the flagged elements. This will
747  // take up some space, maybe more than what was freed.
748  const bool mesh_changed =
749  this->_refine_elements();
750 
754  // FIXME: This won't pass unless we add a redundant find_neighbors()
755  // call or replace find_neighbors() with on-the-fly neighbor updating
756  // libmesh_assert(!this->eliminate_unrefined_patches());
757 
758  // Finally, the new mesh needs to be prepared for use
759  if (mesh_changed)
760  _mesh.prepare_for_use (/*skip_renumber =*/false);
761 
762  return mesh_changed;
763 }
764 
765 
766 
768 {
769  // This function must be run on all processors at once
770  parallel_object_only();
771 
772  LOG_SCOPE ("make_flags_parallel_consistent()", "MeshRefinement");
773 
777  (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
778 
782  (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), psync);
783 
784  // If we weren't consistent in both h and p on every processor then
785  // we weren't globally consistent
786  bool parallel_consistent = hsync.parallel_consistent &&
787  psync.parallel_consistent;
788  this->comm().min(parallel_consistent);
789 
790  return parallel_consistent;
791 }
792 
793 
794 
796 {
797  // This function must be run on all processors at once
798  parallel_object_only();
799 
800  // We may need a PointLocator for topological_neighbor() tests
801  // later, which we need to make sure gets constructed on all
802  // processors at once.
803  UniquePtr<PointLocatorBase> point_locator;
804 
805 #ifdef LIBMESH_ENABLE_PERIODIC
806  bool has_periodic_boundaries =
808  libmesh_assert(this->comm().verify(has_periodic_boundaries));
809 
810  if (has_periodic_boundaries)
811  point_locator = _mesh.sub_point_locator();
812 #endif
813 
814  LOG_SCOPE ("make_coarsening_compatible()", "MeshRefinement");
815 
816  // Unless we encounter a specific situation level-one
817  // will be satisfied after executing this loop just once
818  bool level_one_satisfied = true;
819 
820 
821  // Unless we encounter a specific situation we will be compatible
822  // with any selected refinement flags
823  bool compatible_with_refinement = true;
824 
825 
826  // find the maximum h and p levels in the mesh
827  unsigned int max_level = 0;
828  unsigned int max_p_level = 0;
829 
830  {
831  // First we look at all the active level-0 elements. Since it doesn't make
832  // sense to coarsen them we must un-set their coarsen flags if
833  // they are set.
834  for (auto & elem : _mesh.active_element_ptr_range())
835  {
836  max_level = std::max(max_level, elem->level());
837  max_p_level =
838  std::max(max_p_level,
839  static_cast<unsigned int>(elem->p_level()));
840 
841  if ((elem->level() == 0) &&
842  (elem->refinement_flag() == Elem::COARSEN))
844 
845  if ((elem->p_level() == 0) &&
846  (elem->p_refinement_flag() == Elem::COARSEN))
848  }
849  }
850 
851  // Even if there are no refined elements on this processor then
852  // there may still be work for us to do on e.g. ancestor elements.
853  // At the very least we need to be in the loop if a distributed mesh
854  // needs to synchronize data.
855 #if 0
856  if (max_level == 0 && max_p_level == 0)
857  {
858  // But we still have to check with other processors
859  this->comm().min(compatible_with_refinement);
860 
861  return compatible_with_refinement;
862  }
863 #endif
864 
865  // Loop over all the active elements. If an element is marked
866  // for coarsening we better check its neighbors. If ANY of these neighbors
867  // are marked for refinement AND are at the same level then there is a
868  // conflict. By convention refinement wins, so we un-mark the element for
869  // coarsening. Level-one would be violated in this case so we need to re-run
870  // the loop.
872  {
873 
874  repeat:
875  level_one_satisfied = true;
876 
877  do
878  {
879  level_one_satisfied = true;
880 
881  for (auto & elem : _mesh.active_element_ptr_range())
882  {
883  bool my_flag_changed = false;
884 
885  if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
886  // the coarsen flag is set
887  {
888  const unsigned int my_level = elem->level();
889 
890  for (auto n : elem->side_index_range())
891  {
892  const Elem * neighbor =
893  topological_neighbor(elem, point_locator.get(), n);
894 
895  if (neighbor != libmesh_nullptr && // I have a
896  neighbor != remote_elem) // neighbor here
897  {
898  if (neighbor->active()) // and it is active
899  {
900  if ((neighbor->level() == my_level) &&
901  (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
902  // and wants to be refined
903  {
905  my_flag_changed = true;
906  break;
907  }
908  }
909  else // I have a neighbor and it is not active. That means it has children.
910  { // While it _may_ be possible to coarsen us if all the children of
911  // that element want to be coarsened, it is impossible to know at this
912  // stage. Forget about it for the moment... This can be handled in
913  // two steps.
915  my_flag_changed = true;
916  break;
917  }
918  }
919  }
920  }
921  if (elem->p_refinement_flag() == Elem::COARSEN) // If
922  // the element is active and the order reduction flag is set
923  {
924  const unsigned int my_p_level = elem->p_level();
925 
926  for (auto n : elem->side_index_range())
927  {
928  const Elem * neighbor =
929  topological_neighbor(elem, point_locator.get(), n);
930 
931  if (neighbor != libmesh_nullptr && // I have a
932  neighbor != remote_elem) // neighbor here
933  {
934  if (neighbor->active()) // and it is active
935  {
936  if ((neighbor->p_level() > my_p_level &&
937  neighbor->p_refinement_flag() != Elem::COARSEN)
938  || (neighbor->p_level() == my_p_level &&
939  neighbor->p_refinement_flag() == Elem::REFINE))
940  {
942  my_flag_changed = true;
943  break;
944  }
945  }
946  else // I have a neighbor and it is not active.
947  { // We need to find which of its children
948  // have me as a neighbor, and maintain
949  // level one p compatibility with them.
950  // Because we currently have level one h
951  // compatibility, we don't need to check
952  // grandchildren
953 
954  libmesh_assert(neighbor->has_children());
955  for (auto & subneighbor : neighbor->child_ref_range())
956  if (&subneighbor != remote_elem &&
957  subneighbor.active() &&
958  has_topological_neighbor(&subneighbor, point_locator.get(), elem))
959  if ((subneighbor.p_level() > my_p_level &&
960  subneighbor.p_refinement_flag() != Elem::COARSEN)
961  || (subneighbor.p_level() == my_p_level &&
962  subneighbor.p_refinement_flag() == Elem::REFINE))
963  {
965  my_flag_changed = true;
966  break;
967  }
968  if (my_flag_changed)
969  break;
970  }
971  }
972  }
973  }
974 
975  // If the current element's flag changed, we hadn't
976  // satisfied the level one rule.
977  if (my_flag_changed)
978  level_one_satisfied = false;
979 
980  // Additionally, if it has non-local neighbors, and
981  // we're not in serial, then we'll eventually have to
982  // return compatible_with_refinement = false, because
983  // our change has to propagate to neighboring
984  // processors.
985  if (my_flag_changed && !_mesh.is_serial())
986  for (auto n : elem->side_index_range())
987  {
988  Elem * neigh =
989  topological_neighbor(elem, point_locator.get(), n);
990 
991  if (!neigh)
992  continue;
993  if (neigh == remote_elem ||
994  neigh->processor_id() !=
995  this->processor_id())
996  {
997  compatible_with_refinement = false;
998  break;
999  }
1000  // FIXME - for non-level one meshes we should
1001  // test all descendants
1002  if (neigh->has_children())
1003  for (auto & child : neigh->child_ref_range())
1004  if (&child == remote_elem ||
1005  child.processor_id() !=
1006  this->processor_id())
1007  {
1008  compatible_with_refinement = false;
1009  break;
1010  }
1011  }
1012  }
1013  }
1014  while (!level_one_satisfied);
1015 
1016  } // end if (_face_level_mismatch_limit)
1017 
1018 
1019  // Next we look at all of the ancestor cells.
1020  // If there is a parent cell with all of its children
1021  // wanting to be unrefined then the element is a candidate
1022  // for unrefinement. If all the children don't
1023  // all want to be unrefined then ALL of them need to have their
1024  // unrefinement flags cleared.
1025  for (int level=(max_level); level >= 0; level--)
1026  {
1028  const MeshBase::element_iterator end_el = _mesh.level_elements_end(level);
1029  for (; el != end_el; ++el)
1030  {
1031  Elem * elem = *el;
1032  if (elem->ancestor())
1033  {
1034 
1035  // right now the element hasn't been disqualified
1036  // as a candidate for unrefinement
1037  bool is_a_candidate = true;
1038  bool found_remote_child = false;
1039 
1040  for (auto & child : elem->child_ref_range())
1041  {
1042  if (&child == remote_elem)
1043  found_remote_child = true;
1044  else if ((child.refinement_flag() != Elem::COARSEN) ||
1045  !child.active() )
1046  is_a_candidate = false;
1047  }
1048 
1049  if (!is_a_candidate && !found_remote_child)
1050  {
1052 
1053  for (auto & child : elem->child_ref_range())
1054  {
1055  if (&child == remote_elem)
1056  continue;
1057  if (child.refinement_flag() == Elem::COARSEN)
1058  {
1059  level_one_satisfied = false;
1060  child.set_refinement_flag(Elem::DO_NOTHING);
1061  }
1062  }
1063  }
1064  }
1065  }
1066  }
1067 
1068  if (!level_one_satisfied && _face_level_mismatch_limit) goto repeat;
1069 
1070 
1071  // If all the children of a parent are set to be coarsened
1072  // then flag the parent so that they can kill their kids.
1073 
1074  // On a distributed mesh, we won't always be able to determine this
1075  // on parent elements with remote children, even if we own the
1076  // parent, without communication.
1077  //
1078  // We'll first communicate *to* parents' owners when we determine
1079  // they cannot be coarsened, then we'll sync the final refinement
1080  // flag *from* the parents.
1081 
1082  // uncoarsenable_parents[p] live on processor id p
1083  const processor_id_type n_proc = _mesh.n_processors();
1084  const processor_id_type my_proc_id = _mesh.processor_id();
1085  const bool distributed_mesh = !_mesh.is_replicated();
1086 
1087  std::vector<std::vector<dof_id_type>>
1088  uncoarsenable_parents(n_proc);
1089 
1090  MeshBase::element_iterator ancestor_el =
1092  const MeshBase::element_iterator ancestor_el_end =
1094  for (; ancestor_el != ancestor_el_end; ++ancestor_el)
1095  {
1096  Elem * elem = *ancestor_el;
1097 
1098  // Presume all the children are flagged for coarsening and
1099  // then look for a contradiction
1100  bool all_children_flagged_for_coarsening = true;
1101 
1102  for (auto & child : elem->child_ref_range())
1103  {
1104  if (&child != remote_elem &&
1105  child.refinement_flag() != Elem::COARSEN)
1106  {
1107  all_children_flagged_for_coarsening = false;
1108  if (!distributed_mesh)
1109  break;
1110  if (child.processor_id() != elem->processor_id())
1111  {
1112  uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
1113  break;
1114  }
1115  }
1116  }
1117 
1118  if (all_children_flagged_for_coarsening)
1120  else
1122  }
1123 
1124  // If we have a distributed mesh, we might need to sync up
1125  // INACTIVE vs. COARSEN_INACTIVE flags.
1126  if (distributed_mesh)
1127  {
1128  // We'd better still be in sync here
1129  parallel_object_only();
1130 
1132  uncoarsenable_tag = this->comm().get_unique_tag(2718);
1133  std::vector<Parallel::Request> uncoarsenable_push_requests(n_proc-1);
1134 
1135  for (processor_id_type p = 0; p != n_proc; ++p)
1136  {
1137  if (p == my_proc_id)
1138  continue;
1139 
1141  uncoarsenable_push_requests[p - (p > my_proc_id)];
1142 
1143  _mesh.comm().send
1144  (p, uncoarsenable_parents[p], request, uncoarsenable_tag);
1145  }
1146 
1147  for (processor_id_type p = 1; p != n_proc; ++p)
1148  {
1149  std::vector<dof_id_type> my_uncoarsenable_parents;
1150  _mesh.comm().receive
1151  (Parallel::any_source, my_uncoarsenable_parents,
1152  uncoarsenable_tag);
1153 
1154  for (std::vector<dof_id_type>::const_iterator
1155  it = my_uncoarsenable_parents.begin(),
1156  end = my_uncoarsenable_parents.end(); it != end; ++it)
1157  {
1158  Elem & elem = _mesh.elem_ref(*it);
1162  }
1163  }
1164 
1165  Parallel::wait(uncoarsenable_push_requests);
1166 
1170  (this->comm(), _mesh.not_local_elements_begin(),
1172  // We'd like a smaller sync, but this leads to bugs?
1173  // SyncCoarsenInactive(),
1174  hsync);
1175  }
1176 
1177  // If one processor finds an incompatibility, we're globally
1178  // incompatible
1179  this->comm().min(compatible_with_refinement);
1180 
1181  return compatible_with_refinement;
1182 }
1183 
1184 
1185 
1186 
1187 
1188 
1189 
1190 
1192 {
1193  // This function must be run on all processors at once
1194  parallel_object_only();
1195 
1196  // We may need a PointLocator for topological_neighbor() tests
1197  // later, which we need to make sure gets constructed on all
1198  // processors at once.
1199  UniquePtr<PointLocatorBase> point_locator;
1200 
1201 #ifdef LIBMESH_ENABLE_PERIODIC
1202  bool has_periodic_boundaries =
1204  libmesh_assert(this->comm().verify(has_periodic_boundaries));
1205 
1206  if (has_periodic_boundaries)
1207  point_locator = _mesh.sub_point_locator();
1208 #endif
1209 
1210  LOG_SCOPE ("make_refinement_compatible()", "MeshRefinement");
1211 
1212  // Unless we encounter a specific situation we will be compatible
1213  // with any selected coarsening flags
1214  bool compatible_with_coarsening = true;
1215 
1216  // This loop enforces the level-1 rule. We should only
1217  // execute it if the user indeed wants level-1 satisfied!
1219  {
1220  // Unless we encounter a specific situation level-one
1221  // will be satisfied after executing this loop just once
1222  bool level_one_satisfied = true;
1223 
1224  do
1225  {
1226  level_one_satisfied = true;
1227 
1228  for (auto & elem : _mesh.active_element_ptr_range())
1229  {
1230  const unsigned short n_sides = elem->n_sides();
1231 
1232  if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1233  // h refinement flag is set
1234  {
1235  const unsigned int my_level = elem->level();
1236 
1237  for (unsigned short side = 0; side != n_sides;
1238  ++side)
1239  {
1240  Elem * neighbor =
1241  topological_neighbor(elem, point_locator.get(), side);
1242 
1243  if (neighbor != libmesh_nullptr && // I have a
1244  neighbor != remote_elem && // neighbor here
1245  neighbor->active()) // and it is active
1246  {
1247  // Case 1: The neighbor is at the same level I am.
1248  // 1a: The neighbor will be refined -> NO PROBLEM
1249  // 1b: The neighbor won't be refined -> NO PROBLEM
1250  // 1c: The neighbor wants to be coarsened -> PROBLEM
1251  if (neighbor->level() == my_level)
1252  {
1253  if (neighbor->refinement_flag() == Elem::COARSEN)
1254  {
1256  if (neighbor->parent())
1258  compatible_with_coarsening = false;
1259  level_one_satisfied = false;
1260  }
1261  }
1262 
1263 
1264  // Case 2: The neighbor is one level lower than I am.
1265  // The neighbor thus MUST be refined to satisfy
1266  // the level-one rule, regardless of whether it
1267  // was originally flagged for refinement. If it
1268  // wasn't flagged already we need to repeat
1269  // this process.
1270  else if ((neighbor->level()+1) == my_level)
1271  {
1272  if (neighbor->refinement_flag() != Elem::REFINE)
1273  {
1274  neighbor->set_refinement_flag(Elem::REFINE);
1275  if (neighbor->parent())
1277  compatible_with_coarsening = false;
1278  level_one_satisfied = false;
1279  }
1280  }
1281 #ifdef DEBUG
1282  // Note that the only other possibility is that the
1283  // neighbor is already refined, in which case it isn't
1284  // active and we should never get here.
1285  else
1286  libmesh_error_msg("ERROR: Neighbor level must be equal or 1 higher than mine.");
1287 #endif
1288  }
1289  }
1290  }
1291  if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1292  // p refinement flag is set
1293  {
1294  const unsigned int my_p_level = elem->p_level();
1295 
1296  for (unsigned int side=0; side != n_sides; side++)
1297  {
1298  Elem * neighbor =
1299  topological_neighbor(elem, point_locator.get(), side);
1300 
1301  if (neighbor != libmesh_nullptr && // I have a
1302  neighbor != remote_elem) // neighbor here
1303  {
1304  if (neighbor->active()) // and it is active
1305  {
1306  if (neighbor->p_level() < my_p_level &&
1307  neighbor->p_refinement_flag() != Elem::REFINE)
1308  {
1310  level_one_satisfied = false;
1311  compatible_with_coarsening = false;
1312  }
1313  if (neighbor->p_level() == my_p_level &&
1314  neighbor->p_refinement_flag() == Elem::COARSEN)
1315  {
1317  level_one_satisfied = false;
1318  compatible_with_coarsening = false;
1319  }
1320  }
1321  else // I have an inactive neighbor
1322  {
1323  libmesh_assert(neighbor->has_children());
1324  for (auto & subneighbor : neighbor->child_ref_range())
1325  if (&subneighbor != remote_elem && subneighbor.active() &&
1326  has_topological_neighbor(&subneighbor, point_locator.get(), elem))
1327  {
1328  if (subneighbor.p_level() < my_p_level &&
1329  subneighbor.p_refinement_flag() != Elem::REFINE)
1330  {
1331  // We should already be level one
1332  // compatible
1333  libmesh_assert_greater (subneighbor.p_level() + 2u,
1334  my_p_level);
1335  subneighbor.set_p_refinement_flag(Elem::REFINE);
1336  level_one_satisfied = false;
1337  compatible_with_coarsening = false;
1338  }
1339  if (subneighbor.p_level() == my_p_level &&
1340  subneighbor.p_refinement_flag() == Elem::COARSEN)
1341  {
1342  subneighbor.set_p_refinement_flag(Elem::DO_NOTHING);
1343  level_one_satisfied = false;
1344  compatible_with_coarsening = false;
1345  }
1346  }
1347  }
1348  }
1349  }
1350  }
1351  }
1352  }
1353 
1354  while (!level_one_satisfied);
1355  } // end if (_face_level_mismatch_limit)
1356 
1357  // If we're not compatible on one processor, we're globally not
1358  // compatible
1359  this->comm().min(compatible_with_coarsening);
1360 
1361  return compatible_with_coarsening;
1362 }
1363 
1364 
1365 
1366 
1368 {
1369  // This function must be run on all processors at once
1370  parallel_object_only();
1371 
1372  LOG_SCOPE ("_coarsen_elements()", "MeshRefinement");
1373 
1374  // Flags indicating if this call actually changes the mesh
1375  bool mesh_changed = false;
1376  bool mesh_p_changed = false;
1377 
1378  // Clear the unused_elements data structure.
1379  // The elements have been packed since it was built,
1380  // so there are _no_ unused elements. We cannot trust
1381  // any iterators currently in this data structure.
1382  // _unused_elements.clear();
1383 
1384  // Loop over the elements first to determine if the mesh will
1385  // undergo h-coarsening. If it will, then we'll need to communicate
1386  // more ghosted elements. We need to communicate them *before* we
1387  // do the coarsening; otherwise it is possible to coarsen away a
1388  // one-element-thick layer partition and leave the partitions on
1389  // either side unable to figure out how to talk to each other.
1390  for (auto & elem : _mesh.element_ptr_range())
1391  if (elem->refinement_flag() == Elem::COARSEN)
1392  {
1393  mesh_changed = true;
1394  break;
1395  }
1396 
1397  // If the mesh changed on any processor, it changed globally
1398  this->comm().max(mesh_changed);
1399 
1400  // And then we may need to widen the ghosting layers.
1401  if (mesh_changed)
1403 
1404  for (auto & elem : _mesh.element_ptr_range())
1405  {
1406  // active elements flagged for coarsening will
1407  // no longer be deleted until MeshRefinement::contract()
1408  if (elem->refinement_flag() == Elem::COARSEN)
1409  {
1410  // Huh? no level-0 element should be active
1411  // and flagged for coarsening.
1412  libmesh_assert_not_equal_to (elem->level(), 0);
1413 
1414  // Remove this element from any neighbor
1415  // lists that point to it.
1416  elem->nullify_neighbors();
1417 
1418  // Remove any boundary information associated
1419  // with this element
1420  _mesh.get_boundary_info().remove (elem);
1421 
1422  // Add this iterator to the _unused_elements
1423  // data structure so we might fill it.
1424  // The _unused_elements optimization is currently off.
1425  // _unused_elements.push_back (it);
1426 
1427  // Don't delete the element until
1428  // MeshRefinement::contract()
1429  // _mesh.delete_elem(elem);
1430  }
1431 
1432  // inactive elements flagged for coarsening
1433  // will become active
1434  else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1435  {
1436  elem->coarsen();
1437  libmesh_assert (elem->active());
1438 
1439  // the mesh has certainly changed
1440  mesh_changed = true;
1441  }
1442  if (elem->p_refinement_flag() == Elem::COARSEN)
1443  {
1444  if (elem->p_level() > 0)
1445  {
1447  elem->set_p_level(elem->p_level() - 1);
1448  mesh_p_changed = true;
1449  }
1450  else
1451  {
1453  }
1454  }
1455  }
1456 
1457  this->comm().max(mesh_p_changed);
1458 
1459  // And we may need to update DistributedMesh values reflecting the changes
1460  if (mesh_changed)
1462 
1463  // Node processor ids may need to change if an element of that id
1464  // was coarsened away
1465  if (mesh_changed && !_mesh.is_serial())
1466  {
1467  // Update the _new_nodes_map so that processors can
1468  // find requested nodes
1469  this->update_nodes_map ();
1470 
1472 
1473  // Clear the _new_nodes_map
1474  this->clear();
1475 
1476 #ifdef DEBUG
1477  MeshTools::libmesh_assert_valid_procids<Node>(_mesh);
1478 #endif
1479  }
1480 
1481  // If p levels changed all we need to do is make sure that parent p
1482  // levels changed in sync
1483  if (mesh_p_changed && !_mesh.is_serial())
1484  {
1486  }
1487 
1488  return (mesh_changed || mesh_p_changed);
1489 }
1490 
1491 
1492 
1494 {
1495  // This function must be run on all processors at once
1496  parallel_object_only();
1497 
1498  // Update the _new_nodes_map so that elements can
1499  // find nodes to connect to.
1500  this->update_nodes_map ();
1501 
1502  LOG_SCOPE ("_refine_elements()", "MeshRefinement");
1503 
1504  // Iterate over the elements, counting the elements
1505  // flagged for h refinement.
1506  dof_id_type n_elems_flagged = 0;
1507 
1508  for (auto & elem : _mesh.element_ptr_range())
1509  if (elem->refinement_flag() == Elem::REFINE)
1510  n_elems_flagged++;
1511 
1512  // Construct a local vector of Elem * which have been
1513  // previously marked for refinement. We reserve enough
1514  // space to allow for every element to be refined.
1515  std::vector<Elem *> local_copy_of_elements;
1516  local_copy_of_elements.reserve(n_elems_flagged);
1517 
1518  // If mesh p levels changed, we might need to synchronize parent p
1519  // levels on a distributed mesh.
1520  bool mesh_p_changed = false;
1521 
1522  // Iterate over the elements, looking for elements flagged for
1523  // refinement.
1524 
1525  // If we are on a ReplicatedMesh, then we just do the refinement in
1526  // the same order on every processor and everything stays in sync.
1527 
1528  // If we are on a DistributedMesh, that's impossible.
1529  //
1530  // If the mesh is distributed, we need to make sure that if we end
1531  // up as the owner of a new node, which might happen if that node is
1532  // attached to one of our own elements, then we have given it a
1533  // legitimate node id and our own processor id. We generate
1534  // legitimate node ids and use our own processor id when we are
1535  // refining our own elements but not when we refine others'
1536  // elements. Therefore we want to refine our own elements *first*,
1537  // thereby generating all nodes which might belong to us, and then
1538  // refine others' elements *after*, thereby generating nodes with
1539  // temporary ids which we know we will discard.
1540  //
1541  // Even if the DistributedMesh is serialized, we can't just treat it
1542  // like a ReplicatedMesh, because DistributedMesh doesn't *trust*
1543  // users to refine partitioned elements in a serialized way, so it
1544  // assigns temporary ids, so we need to synchronize ids afterward to
1545  // be safe anyway, so we might as well use the distributed mesh code
1546  // path.
1548  {
1549  if (elem->refinement_flag() == Elem::REFINE)
1550  local_copy_of_elements.push_back(elem);
1551  if (elem->p_refinement_flag() == Elem::REFINE &&
1552  elem->active())
1553  {
1554  elem->set_p_level(elem->p_level()+1);
1556  mesh_p_changed = true;
1557  }
1558  }
1559 
1560  if (!_mesh.is_replicated())
1561  {
1564  elem_end = _mesh.active_not_local_elements_end();
1565  elem_it != elem_end; ++elem_it)
1566  {
1567  Elem * elem = *elem_it;
1568  if (elem->refinement_flag() == Elem::REFINE)
1569  local_copy_of_elements.push_back(elem);
1570  if (elem->p_refinement_flag() == Elem::REFINE &&
1571  elem->active())
1572  {
1573  elem->set_p_level(elem->p_level()+1);
1575  mesh_p_changed = true;
1576  }
1577  }
1578  }
1579 
1580  // Now iterate over the local copies and refine each one.
1581  // This may resize the mesh's internal container and invalidate
1582  // any existing iterators.
1583 
1584  for (std::size_t e = 0; e != local_copy_of_elements.size(); ++e)
1585  local_copy_of_elements[e]->refine(*this);
1586 
1587  // The mesh changed if there were elements h refined
1588  bool mesh_changed = !local_copy_of_elements.empty();
1589 
1590  // If the mesh changed on any processor, it changed globally
1591  this->comm().max(mesh_changed);
1592  this->comm().max(mesh_p_changed);
1593 
1594  // And we may need to update DistributedMesh values reflecting the changes
1595  if (mesh_changed)
1597 
1598  if (mesh_changed && !_mesh.is_replicated())
1599  {
1602 #ifdef DEBUG
1604 #endif
1605  }
1606 
1607  if (mesh_p_changed && !_mesh.is_replicated())
1608  {
1610  }
1611 
1612  // Clear the _new_nodes_map and _unused_elements data structures.
1613  this->clear();
1614 
1615  return (mesh_changed || mesh_p_changed);
1616 }
1617 
1618 
1619 void MeshRefinement::_smooth_flags(bool refining, bool coarsening)
1620 {
1621  // Smoothing can break in weird ways on a mesh with broken topology
1622 #ifdef DEBUG
1624 #endif
1625 
1626  // Repeat until flag changes match on every processor
1627  do
1628  {
1629  // Repeat until coarsening & refinement flags jive
1630  bool satisfied = false;
1631  do
1632  {
1633  // If we're refining or coarsening, hit the corresponding
1634  // face level test code. Short-circuiting || is our friend
1635  const bool coarsening_satisfied =
1636  !coarsening ||
1638 
1639  const bool refinement_satisfied =
1640  !refining ||
1642 
1643  bool smoothing_satisfied =
1644  !this->eliminate_unrefined_patches();// &&
1645 
1647  smoothing_satisfied = smoothing_satisfied &&
1649 
1651  smoothing_satisfied = smoothing_satisfied &&
1653 
1655  smoothing_satisfied = smoothing_satisfied &&
1657 
1659  smoothing_satisfied = smoothing_satisfied &&
1661 
1662  satisfied = (coarsening_satisfied &&
1663  refinement_satisfied &&
1664  smoothing_satisfied);
1665 
1666  libmesh_assert(this->comm().verify(satisfied));
1667  }
1668  while (!satisfied);
1669  }
1670  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
1671 }
1672 
1673 
1675 {
1676  // Refine n times
1677  for (unsigned int rstep=0; rstep<n; rstep++)
1678  for (auto & elem : _mesh.active_element_ptr_range())
1679  {
1680  // P refine all the active elements
1681  elem->set_p_level(elem->p_level()+1);
1683  }
1684 }
1685 
1686 
1687 
1689 {
1690  // Coarsen p times
1691  for (unsigned int rstep=0; rstep<n; rstep++)
1692  for (auto & elem : _mesh.active_element_ptr_range())
1693  if (elem->p_level() > 0)
1694  {
1695  // P coarsen all the active elements
1696  elem->set_p_level(elem->p_level()-1);
1698  }
1699 }
1700 
1701 
1702 
1704 {
1705  // Refine n times
1706  // FIXME - this won't work if n>1 and the mesh
1707  // has already been attached to an equation system
1708  for (unsigned int rstep=0; rstep<n; rstep++)
1709  {
1710  // Clean up the refinement flags
1711  this->clean_refinement_flags();
1712 
1713  // Flag all the active elements for refinement.
1714  for (auto & elem : _mesh.active_element_ptr_range())
1716 
1717  // Refine all the elements we just flagged.
1718  this->_refine_elements();
1719  }
1720 
1721  // Finally, the new mesh probably needs to be prepared for use
1722  if (n > 0)
1723  _mesh.prepare_for_use (/*skip_renumber =*/false);
1724 }
1725 
1726 
1727 
1729 {
1730  // Coarsen n times
1731  for (unsigned int rstep=0; rstep<n; rstep++)
1732  {
1733  // Clean up the refinement flags
1734  this->clean_refinement_flags();
1735 
1736  // Flag all the active elements for coarsening.
1737  for (auto & elem : _mesh.active_element_ptr_range())
1738  {
1740  if (elem->parent())
1742  }
1743 
1744  // On a distributed mesh, we may have parent elements with
1745  // remote active children. To keep flags consistent, we'll need
1746  // a communication step.
1747  if (!_mesh.is_replicated())
1748  {
1749  const processor_id_type n_proc = _mesh.n_processors();
1750  const processor_id_type my_proc_id = _mesh.processor_id();
1751 
1752  std::vector<std::vector<dof_id_type>>
1753  parents_to_coarsen(n_proc);
1754 
1757 
1758  for ( ; elem_it != elem_end; ++elem_it)
1759  {
1760  const Elem & elem = **elem_it;
1761  if (elem.processor_id() != my_proc_id &&
1763  parents_to_coarsen[elem.processor_id()].push_back(elem.id());
1764  }
1765 
1767  coarsen_tag = this->comm().get_unique_tag(271);
1768  std::vector<Parallel::Request> coarsen_push_requests(n_proc-1);
1769 
1770  for (processor_id_type p = 0; p != n_proc; ++p)
1771  {
1772  if (p == my_proc_id)
1773  continue;
1774 
1776  coarsen_push_requests[p - (p > my_proc_id)];
1777 
1778  _mesh.comm().send
1779  (p, parents_to_coarsen[p], request, coarsen_tag);
1780  }
1781 
1782  for (processor_id_type p = 1; p != n_proc; ++p)
1783  {
1784  std::vector<dof_id_type> my_parents_to_coarsen;
1785  _mesh.comm().receive
1786  (Parallel::any_source, my_parents_to_coarsen,
1787  coarsen_tag);
1788 
1789  for (std::vector<dof_id_type>::const_iterator
1790  it = my_parents_to_coarsen.begin(),
1791  end = my_parents_to_coarsen.end(); it != end; ++it)
1792  {
1793  Elem & elem = _mesh.elem_ref(*it);
1797  }
1798  }
1799 
1800  Parallel::wait(coarsen_push_requests);
1801 
1805  (this->comm(), _mesh.not_local_elements_begin(),
1807  // We'd like a smaller sync, but this leads to bugs?
1808  // SyncCoarsenInactive(),
1809  hsync);
1810  }
1811 
1812  // Coarsen all the elements we just flagged.
1813  this->_coarsen_elements();
1814  }
1815 
1816 
1817  // Finally, the new mesh probably needs to be prepared for use
1818  if (n > 0)
1819  _mesh.prepare_for_use (/*skip_renumber =*/false);
1820 }
1821 
1822 
1823 
1825  const PointLocatorBase * point_locator,
1826  const unsigned int side)
1827 {
1828 #ifdef LIBMESH_ENABLE_PERIODIC
1829  if (_periodic_boundaries && !_periodic_boundaries->empty())
1830  {
1831  libmesh_assert(point_locator);
1832  return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1833  }
1834 #endif
1835  return elem->neighbor_ptr(side);
1836 }
1837 
1838 
1839 
1841  const PointLocatorBase * point_locator,
1842  const Elem * neighbor)
1843 {
1844 #ifdef LIBMESH_ENABLE_PERIODIC
1845  if (_periodic_boundaries && !_periodic_boundaries->empty())
1846  {
1847  libmesh_assert(point_locator);
1848  return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1849  }
1850 #endif
1851  return elem->has_neighbor(neighbor);
1852 }
1853 
1854 
1855 
1856 } // namespace libMesh
1857 
1858 
1859 #endif
bool limit_level_mismatch_at_edge(const unsigned int max_mismatch)
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
bool ancestor() const
Definition: elem.C:1574
~MeshRefinement()
Destructor.
void set_p_level(const unsigned int p)
Sets the value of the p-refinement level for the element.
Definition: elem.h:2559
virtual element_iterator ancestor_elements_begin()=0
Iterate over elements for which elem->ancestor() is true.
bool has_children() const
Definition: elem.h:2295
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
Status wait(Request &r)
Wait for a non-blocking send or receive to finish.
Definition: parallel.h:565
bool _refine_elements()
Refines user-requested elements.
virtual void libmesh_assert_valid_parallel_ids() const
Verify id and processor_id consistency of our elements and nodes containers.
Definition: mesh_base.h:965
A Node is like a Point, but with more information.
Definition: node.h:52
virtual bool is_serial() const
Definition: mesh_base.h:140
bool active() const
Definition: elem.h:2257
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
virtual element_iterator not_local_elements_end()=0
bool limit_level_mismatch_at_node(const unsigned int max_mismatch)
This algorithm restricts the maximum level mismatch at any node in the mesh.
dof_id_type find(dof_id_type bracket_node1, dof_id_type bracket_node2) const
Definition: topology_map.C:118
virtual element_iterator level_elements_begin(unsigned int level)=0
Iterate over elements of a given level.
void update_nodes_map()
Updates the _new_nodes_map.
const unsigned int any_source
Accept from any source.
Definition: parallel.h:204
TopologyMap _new_nodes_map
Data structure that holds the new nodes information.
bool has_topological_neighbor(const Elem *elem, const MeshBase &mesh, const PointLocatorBase &point_locator, const PeriodicBoundaries *pb) const
Definition: elem.C:1104
bool limit_underrefined_boundary(const signed char max_mismatch)
unsigned int p_level() const
Definition: elem.h:2422
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
MeshBase & _mesh
Reference to the mesh.
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
void min(T &r) const
Take a local variable and replace it with the minimum of it&#39;s values on all processors.
virtual bool is_replicated() const
Definition: mesh_base.h:162
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...
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
bool refine_elements()
Only refines the user-requested elements.
void make_elems_parallel_consistent(MeshBase &)
Copy ids of ghost elements from their local processors.
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2065
const Elem * parent() const
Definition: elem.h:2346
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:52
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
unsigned short int side
Definition: xdr_io.C:49
processor_id_type n_processors() const
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
virtual element_iterator ancestor_elements_end()=0
uint8_t processor_id_type
Definition: id_types.h:99
MPI_Request request
Request object for non-blocking I/O.
Definition: parallel.h:171
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:1125
const class libmesh_nullptr_t libmesh_nullptr
void set_refinement_flag(const RefinementState rflag)
Sets the value of the refinement flag for the element.
Definition: elem.h:2513
void init(MeshBase &)
Definition: topology_map.C:35
bool limit_overrefined_boundary(const signed char max_mismatch)
virtual const Node * node_ptr(const dof_id_type i) const =0
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.
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
The libMesh namespace provides an interface to certain functionality in the library.
void send_coarse_ghosts(MeshBase &) const
Examine a just-coarsened mesh, and for any newly-coarsened elements, send the associated ghosted elem...
long double max(long double a, double b)
PeriodicBoundaries * _periodic_boundaries
unsigned int max_level(const MeshBase &mesh)
Find the maximum h-refinement level in a mesh.
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.
This is the MeshBase class.
Definition: mesh_base.h:68
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:1699
bool make_refinement_compatible()
Take user-specified refinement flags and augment them so that level-one dependency is satisfied...
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
bool coarsen_elements()
Only coarsens the user-requested elements.
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
virtual element_iterator elements_begin()=0
Iterate over all the elements in the Mesh.
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
virtual element_iterator level_elements_end(unsigned int level)=0
const Elem * topological_neighbor(const unsigned int i, const MeshBase &mesh, const PointLocatorBase &point_locator, const PeriodicBoundaries *pb) const
Definition: elem.C:1067
bool test_level_one(bool libmesh_assert_yes=false)
This is the MeshCommunication class.
virtual SimpleRange< element_iterator > element_ptr_range()=0
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
unsigned char _face_level_mismatch_limit
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
Encapsulates the MPI tag integers.
Definition: parallel.h:227
virtual void update_parallel_id_counts()=0
Updates parallel caches so that methods like n_elem() accurately reflect changes on other processors...
virtual element_iterator elements_end()=0
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.
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:324
RefinementState p_refinement_flag() const
Definition: elem.h:2521
bool make_flags_parallel_consistent()
Copy refinement flags on ghost elements from their local processors.
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
virtual const std::vector< std::pair< dof_id_type, dof_id_type > > bracketing_nodes(unsigned int c, unsigned int n) const
Definition: elem.C:2353
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
Blocking-send to one processor with data-defined type.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
This is the base class for point locators.
virtual unsigned int n_sides() const =0
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.
This class forms the base class for all other classes that are expected to be implemented in parallel...
bool has_topological_neighbor(const Elem *elem, const PointLocatorBase *point_locator, const Elem *neighbor)
Local dispatch function for checking the correct has_neighbor function from the Elem class...
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.
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:490
unsigned char _node_level_mismatch_limit
virtual element_iterator active_not_local_elements_end()=0
bool test_unflagged(bool libmesh_assert_yes=false)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
bool verify(const T &r, const Communicator &comm=Communicator_World)
Elem * add_elem(Elem *elem)
Adds the element elem to the mesh.
Encapsulates the MPI_Request.
Definition: parallel.h:517
void coarsen()
Coarsen the element.
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
const Parallel::Communicator & comm() const
OStreamProxy out
bool refine_and_coarsen_elements()
Refines and coarsens user-requested elements.
bool has_neighbor(const Elem *elem) const
Definition: elem.h:2010
unsigned int level() const
Definition: elem.h:2388
RefinementState refinement_flag() const
Definition: elem.h:2505
void clear()
Deletes all the data that are currently stored.
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
Blocking-receive from one processor with data-defined type.
void _smooth_flags(bool refining, bool coarsening)
Smooths refinement flags according to current settings.
virtual element_iterator active_not_local_elements_begin()=0
unsigned char _edge_level_mismatch_limit
void nullify_neighbors()
Replaces this element with NULL for all of its neighbors.
Definition: elem.C:2672
void set_p_refinement_flag(const RefinementState pflag)
Sets the value of the p-refinement flag for the element.
Definition: elem.h:2529
virtual unsigned int as_parent_node(unsigned int c, unsigned int n) const
Definition: elem.C:2080
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:1713
virtual element_iterator not_local_elements_begin()=0
bool _coarsen_elements()
Coarsens user-requested elements.
dof_id_type id() const
Definition: dof_object.h:632
bool eliminate_unrefined_patches()
This algorithm selects an element for refinement if all of its neighbors are (or will be) refined...
void make_new_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on new nodes from their local processors.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
void set_periodic_boundaries_ptr(PeriodicBoundaries *pb_ptr)
Sets the PeriodicBoundaries pointer.
Elem * topological_neighbor(Elem *elem, const PointLocatorBase *point_locator, const unsigned int side)
Local dispatch function for getting the correct topological neighbor from the Elem class...
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
UniquePtr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:534
MessageTag get_unique_tag(int tagvalue) const
Get a tag that is unique to this Communicator.
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
virtual float embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
const RemoteElem * remote_elem
Definition: remote_elem.C:57