mesh_refinement.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
33 #include "libmesh/mesh_base.h"
36 #include "libmesh/parallel.h"
38 #include "libmesh/remote_elem.h"
39 
40 #ifdef DEBUG
41 // Some extra validation for ParallelMesh
42 #include "libmesh/mesh_tools.h"
43 #include "libmesh/parallel_mesh.h"
44 #endif // DEBUG
45 
46 #ifdef LIBMESH_ENABLE_PERIODIC
48 #endif
49 
50 namespace libMesh
51 {
52 
53 //-----------------------------------------------------------------
54 // Mesh refinement methods
56  ParallelObject(m),
57  _mesh(m),
58  _use_member_parameters(false),
59  _coarsen_by_parents(false),
60  _refine_fraction(0.3),
61  _coarsen_fraction(0.0),
62  _max_h_level(libMesh::invalid_uint),
63  _coarsen_threshold(10),
64  _nelem_target(0),
65  _absolute_global_tolerance(0.0),
66  _face_level_mismatch_limit(1),
67  _edge_level_mismatch_limit(0),
68  _node_level_mismatch_limit(0)
69 #ifdef LIBMESH_ENABLE_PERIODIC
70  , _periodic_boundaries(NULL)
71 #endif
72 {
73 }
74 
75 
76 
77 #ifdef LIBMESH_ENABLE_PERIODIC
79 {
80  _periodic_boundaries = pb_ptr;
81 }
82 #endif
83 
84 
85 
87 {
88  this->clear();
89 }
90 
91 
92 
94 {
95  _new_nodes_map.clear();
96 }
97 
98 
99 
102  const Real tol)
103 {
104  START_LOG("add_point()", "MeshRefinement");
105 
106  // Return the node if it already exists
107  Node *node = _new_nodes_map.find(p, tol);
108  if (node)
109  {
110  STOP_LOG("add_point()", "MeshRefinement");
111  return node;
112  }
113 
114  // Add the node, with a default id and the requested
115  // processor_id
116  node = _mesh.add_point (p, DofObject::invalid_id, processor_id);
117 
118  libmesh_assert(node);
119 
120  // Add the node to the map.
121  _new_nodes_map.insert(*node);
122 
123  // Return the address of the new node
124  STOP_LOG("add_point()", "MeshRefinement");
125  return node;
126 }
127 
128 
129 
131 {
132  libmesh_assert(elem);
133 
134 
135 // // If the unused_elements has any iterators from
136 // // old elements, take the first one
137 // if (!_unused_elements.empty())
138 // {
139 // std::vector<Elem*>::iterator it = _unused_elements.front();
140 
141 // *it = elem;
142 
143 // _unused_elements.pop_front();
144 // }
145 
146 // // Otherwise, use the conventional add method
147 // else
148 // {
149 // _mesh.add_elem (elem);
150 // }
151 
152  // The _unused_elements optimization has been turned off.
153  _mesh.add_elem (elem);
154 
155  return elem;
156 }
157 
158 
159 
161  (const ErrorVector& error_per_cell,
162  ErrorVector& error_per_parent,
163  Real& parent_error_min,
164  Real& parent_error_max)
165 {
166  // This function must be run on all processors at once
167  parallel_object_only();
168 
169  // Make sure the input error vector is valid
170 #ifdef DEBUG
171  for (std::size_t i=0; i != error_per_cell.size(); ++i)
172  {
173  libmesh_assert_greater_equal (error_per_cell[i], 0);
174  // isnan() isn't standard C++ yet
175  #ifdef isnan
176  libmesh_assert(!isnan(error_per_cell[i]));
177  #endif
178  }
179 
180  // Use a reference to std::vector to avoid confusing
181  // this->comm().verify
182  std::vector<ErrorVectorReal> &epc = error_per_parent;
183  libmesh_assert(this->comm().verify(epc));
184 #endif // #ifdef DEBUG
185 
186  // error values on uncoarsenable elements will be left at -1
187  error_per_parent.clear();
188  error_per_parent.resize(error_per_cell.size(), 0.0);
189 
190  {
191  // Find which elements are uncoarsenable
192  MeshBase::element_iterator elem_it = _mesh.active_local_elements_begin();
193  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();
194  for (; elem_it != elem_end; ++elem_it)
195  {
196  Elem* elem = *elem_it;
197  Elem* parent = elem->parent();
198 
199  // Active elements are uncoarsenable
200  error_per_parent[elem->id()] = -1.0;
201 
202  // Grandparents and up are uncoarsenable
203  while (parent)
204  {
205  parent = parent->parent();
206  if (parent)
207  {
208  const dof_id_type parentid = parent->id();
209  libmesh_assert_less (parentid, error_per_parent.size());
210  error_per_parent[parentid] = -1.0;
211  }
212  }
213  }
214 
215  // Sync between processors.
216  // Use a reference to std::vector to avoid confusing
217  // this->comm().min
218  std::vector<ErrorVectorReal> &epp = error_per_parent;
219  this->comm().min(epp);
220  }
221 
222  // The parent's error is defined as the square root of the
223  // sum of the children's errors squared, so errors that are
224  // Hilbert norms remain Hilbert norms.
225  //
226  // Because the children may be on different processors, we
227  // calculate local contributions to the parents' errors squared
228  // first, then sum across processors and take the square roots
229  // second.
230  {
231  MeshBase::element_iterator elem_it = _mesh.active_local_elements_begin();
232  const MeshBase::element_iterator elem_end = _mesh.active_local_elements_end();
233 
234  for (; elem_it != elem_end; ++elem_it)
235  {
236  Elem* elem = *elem_it;
237  Elem* parent = elem->parent();
238 
239  // Calculate each contribution to parent cells
240  if (parent)
241  {
242  const dof_id_type parentid = parent->id();
243  libmesh_assert_less (parentid, error_per_parent.size());
244 
245  // If the parent has grandchildren we won't be able to
246  // coarsen it, so forget it. Otherwise, add this child's
247  // contribution to the sum of the squared child errors
248  if (error_per_parent[parentid] != -1.0)
249  error_per_parent[parentid] += (error_per_cell[elem->id()] *
250  error_per_cell[elem->id()]);
251  }
252  }
253  }
254 
255  // Sum the vector across all processors
256  this->comm().sum(static_cast<std::vector<ErrorVectorReal>&>(error_per_parent));
257 
258  // Calculate the min and max as we loop
259  parent_error_min = std::numeric_limits<double>::max();
260  parent_error_max = 0.;
261 
262  for (std::size_t i = 0; i != error_per_parent.size(); ++i)
263  {
264  // If this element isn't a coarsenable parent with error, we
265  // have nothing to do. Just flag it as -1 and move on
266  // Note that this->comm().sum might have left uncoarsenable
267  // elements with error_per_parent=-n_proc, so reset it to
268  // error_per_parent=-1
269  if (error_per_parent[i] < 0.)
270  {
271  error_per_parent[i] = -1.;
272  continue;
273  }
274 
275  // The error estimator might have already given us an
276  // estimate on the coarsenable parent elements; if so then
277  // we want to retain that estimate
278  if (error_per_cell[i])
279  {
280  error_per_parent[i] = error_per_cell[i];
281  continue;
282  }
283  // if not, then e_parent = sqrt(sum(e_child^2))
284  else
285  error_per_parent[i] = std::sqrt(error_per_parent[i]);
286 
287  parent_error_min = std::min (parent_error_min,
288  error_per_parent[i]);
289  parent_error_max = std::max (parent_error_max,
290  error_per_parent[i]);
291  }
292 }
293 
294 
295 
297 {
298  this->_new_nodes_map.init(_mesh);
299 }
300 
301 
302 
303 bool MeshRefinement::test_level_one (bool libmesh_dbg_var(libmesh_assert_pass))
304 {
305  // This function must be run on all processors at once
306  parallel_object_only();
307 
308  // We may need a PointLocator for topological_neighbor() tests
309  // later, which we need to make sure gets constructed on all
310  // processors at once.
311  AutoPtr<PointLocatorBase> point_locator;
312 
313 #ifdef LIBMESH_ENABLE_PERIODIC
314  bool has_periodic_boundaries =
316  libmesh_assert(this->comm().verify(has_periodic_boundaries));
317 
318  if (has_periodic_boundaries)
319  point_locator = _mesh.sub_point_locator();
320 #endif
321 
324 
325  bool failure = false;
326 
327 #ifndef NDEBUG
328  Elem *failed_elem = NULL;
329  Elem *failed_neighbor = NULL;
330 #endif // !NDEBUG
331 
332  for ( ; elem_it != elem_end && !failure; ++elem_it)
333  {
334  // Pointer to the element
335  Elem *elem = *elem_it;
336 
337  for (unsigned int n=0; n<elem->n_neighbors(); n++)
338  {
339  Elem *neighbor =
340  topological_neighbor(elem, point_locator.get(), n);
341 
342  if (!neighbor || !neighbor->active() ||
343  neighbor == remote_elem)
344  continue;
345 
346  if ((neighbor->level() + 1 < elem->level()) ||
347  (neighbor->p_level() + 1 < elem->p_level()) ||
348  (neighbor->p_level() > elem->p_level() + 1))
349  {
350  failure = true;
351 #ifndef NDEBUG
352  failed_elem = elem;
353  failed_neighbor = neighbor;
354 #endif // !NDEBUG
355  break;
356  }
357  }
358  }
359 
360  // If any processor failed, we failed globally
361  this->comm().max(failure);
362 
363  if (failure)
364  {
365  // We didn't pass the level one test, so libmesh_assert that
366  // we're allowed not to
367 #ifndef NDEBUG
368  if (libmesh_assert_pass)
369  {
370  libMesh::out <<
371  "MeshRefinement Level one failure, element: " <<
372  *failed_elem << std::endl;
373  libMesh::out <<
374  "MeshRefinement Level one failure, neighbor: " <<
375  *failed_neighbor << std::endl;
376  }
377 #endif // !NDEBUG
378  libmesh_assert(!libmesh_assert_pass);
379  return false;
380  }
381  return true;
382 }
383 
384 
385 
386 bool MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass))
387 {
388  // This function must be run on all processors at once
389  parallel_object_only();
390 
391  bool found_flag = false;
392 
393  // Search for local flags
396 
397 #ifndef NDEBUG
398  Elem *failed_elem = NULL;
399 #endif
400 
401  for ( ; elem_it != elem_end; ++elem_it)
402  {
403  // Pointer to the element
404  Elem *elem = *elem_it;
405 
406  if (elem->refinement_flag() == Elem::REFINE ||
407  elem->refinement_flag() == Elem::COARSEN ||
408  elem->p_refinement_flag() == Elem::REFINE ||
409  elem->p_refinement_flag() == Elem::COARSEN)
410  {
411  found_flag = true;
412 #ifndef NDEBUG
413  failed_elem = elem;
414 #endif
415  break;
416  }
417  }
418 
419  // If we found a flag on any processor, it counts
420  this->comm().max(found_flag);
421 
422  if (found_flag)
423  {
424 #ifndef NDEBUG
425  if (libmesh_assert_pass)
426  {
427  libMesh::out <<
428  "MeshRefinement test_unflagged failure, element: " <<
429  *failed_elem << std::endl;
430  }
431 #endif
432  // We didn't pass the "elements are unflagged" test,
433  // so libmesh_assert that we're allowed not to
434  libmesh_assert(!libmesh_assert_pass);
435  return false;
436  }
437  return true;
438 }
439 
440 
441 
442 bool MeshRefinement::refine_and_coarsen_elements (const bool maintain_level_one)
443 {
444  // This function must be run on all processors at once
445  parallel_object_only();
446 
447  bool _maintain_level_one = maintain_level_one;
448 
449  // If the user used non-default parameters, let's warn that they're
450  // deprecated
451  if (!maintain_level_one)
452  {
453  libmesh_deprecated();
454  }
455  else
456  _maintain_level_one = _face_level_mismatch_limit;
457 
458  // We can't yet turn a non-level-one mesh into a level-one mesh
459  if (_maintain_level_one)
461 
462  // Possibly clean up the refinement flags from
463  // a previous step
465  const MeshBase::element_iterator elem_end = _mesh.elements_end();
466 
467  for ( ; elem_it != elem_end; ++elem_it)
468  {
469  // Pointer to the element
470  Elem *elem = *elem_it;
471 
472  // Set refinement flag to INACTIVE if the
473  // element isn't active
474  if ( !elem->active())
475  {
478  }
479 
480  // This might be left over from the last step
481  if (elem->refinement_flag() == Elem::JUST_REFINED)
483  }
484 
485  // Parallel consistency has to come first, or coarsening
486  // along processor boundaries might occasionally be falsely
487  // prevented
488  bool flags_were_consistent = this->make_flags_parallel_consistent();
489 
490  // In theory, we should be able to remove the above call, which can
491  // be expensive and should be unnecessary. In practice, doing
492  // consistent flagging in parallel is hard, it's impossible to
493  // verify at the library level if it's being done by user code, and
494  // we don't want to abort large parallel runs in opt mode... but we
495  // do want to warn that they should be fixed.
496  if (!flags_were_consistent)
497  {
498  libMesh::out << "Refinement flags were not consistent between processors!\n"
499  << "Correcting and continuing.";
500  }
501 
502  // Repeat until flag changes match on every processor
503  do
504  {
505  // Repeat until coarsening & refinement flags jive
506  bool satisfied = false;
507  do
508  {
509  const bool coarsening_satisfied =
510  this->make_coarsening_compatible(maintain_level_one);
511 
512  const bool refinement_satisfied =
513  this->make_refinement_compatible(maintain_level_one);
514 
515  bool smoothing_satisfied =
517 
519  smoothing_satisfied = smoothing_satisfied &&
521 
523  smoothing_satisfied = smoothing_satisfied &&
525 
526  satisfied = (coarsening_satisfied &&
527  refinement_satisfied &&
528  smoothing_satisfied);
529 #ifdef DEBUG
530  bool max_satisfied = satisfied,
531  min_satisfied = satisfied;
532  this->comm().max(max_satisfied);
533  this->comm().min(min_satisfied);
534  libmesh_assert_equal_to (satisfied, max_satisfied);
535  libmesh_assert_equal_to (satisfied, min_satisfied);
536 #endif
537  }
538  while (!satisfied);
539  }
540  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
541 
542  // First coarsen the flagged elements.
543  const bool coarsening_changed_mesh =
544  this->_coarsen_elements ();
545 
546 // FIXME: test_level_one now tests consistency across periodic
547 // boundaries, which requires a point_locator, which just got
548 // invalidated by _coarsen_elements() and hasn't yet been cleared by
549 // prepare_for_use().
550 
551 // if (_maintain_level_one)
552 // libmesh_assert(test_level_one(true));
553 // libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
554 // libmesh_assert(this->make_refinement_compatible(maintain_level_one));
555 
556 // FIXME: This won't pass unless we add a redundant find_neighbors()
557 // call or replace find_neighbors() with on-the-fly neighbor updating
558 // libmesh_assert(!this->eliminate_unrefined_patches());
559 
560  // We can't contract the mesh ourselves anymore - a System might
561  // need to restrict old coefficient vectors first
562  // _mesh.contract();
563 
564  // Now refine the flagged elements. This will
565  // take up some space, maybe more than what was freed.
566  const bool refining_changed_mesh =
567  this->_refine_elements();
568 
569  // Finally, the new mesh needs to be prepared for use
570  if (coarsening_changed_mesh || refining_changed_mesh)
571  {
572 #ifdef DEBUG
574 #endif
575 
576  _mesh.prepare_for_use (/*skip_renumber =*/false);
577 
578  if (_maintain_level_one)
581  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
582  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
583 // FIXME: This won't pass unless we add a redundant find_neighbors()
584 // call or replace find_neighbors() with on-the-fly neighbor updating
585 // libmesh_assert(!this->eliminate_unrefined_patches());
586 
587  return true;
588  }
589  else
590  {
591  if (_maintain_level_one)
594  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
595  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
596  }
597 
598  // Otherwise there was no change in the mesh,
599  // let the user know. Also, there is no need
600  // to prepare the mesh for use since it did not change.
601  return false;
602 
603 }
604 
605 
606 
607 
608 
609 
610 
611 bool MeshRefinement::coarsen_elements (const bool maintain_level_one)
612 {
613  // This function must be run on all processors at once
614  parallel_object_only();
615 
616  bool _maintain_level_one = maintain_level_one;
617 
618  // If the user used non-default parameters, let's warn that they're
619  // deprecated
620  if (!maintain_level_one)
621  {
622  libmesh_deprecated();
623  }
624  else
625  _maintain_level_one = _face_level_mismatch_limit;
626 
627  // We can't yet turn a non-level-one mesh into a level-one mesh
628  if (_maintain_level_one)
630 
631  // Possibly clean up the refinement flags from
632  // a previous step
634  const MeshBase::element_iterator elem_end = _mesh.elements_end();
635 
636  for ( ; elem_it != elem_end; ++elem_it)
637  {
638  // Pointer to the element
639  Elem* elem = *elem_it;
640 
641  // Set refinement flag to INACTIVE if the
642  // element isn't active
643  if ( !elem->active())
644  {
647  }
648 
649  // This might be left over from the last step
650  if (elem->refinement_flag() == Elem::JUST_REFINED)
652  }
653 
654  // Parallel consistency has to come first, or coarsening
655  // along processor boundaries might occasionally be falsely
656  // prevented
657  bool flags_were_consistent = this->make_flags_parallel_consistent();
658 
659  // In theory, we should be able to remove the above call, which can
660  // be expensive and should be unnecessary. In practice, doing
661  // consistent flagging in parallel is hard, it's impossible to
662  // verify at the library level if it's being done by user code, and
663  // we don't want to abort large parallel runs in opt mode... but we
664  // do want to warn that they should be fixed.
665  libmesh_assert(flags_were_consistent);
666  if (!flags_were_consistent)
667  {
668  libMesh::out << "Refinement flags were not consistent between processors!\n"
669  << "Correcting and continuing.";
670  }
671 
672  // Repeat until flag changes match on every processor
673  do
674  {
675  // Repeat until the flags form a conforming mesh.
676  bool satisfied = false;
677  do
678  {
679  const bool coarsening_satisfied =
680  this->make_coarsening_compatible(maintain_level_one);
681 
682  bool smoothing_satisfied =
683  !this->eliminate_unrefined_patches();// &&
684 
686  smoothing_satisfied = smoothing_satisfied &&
688 
690  smoothing_satisfied = smoothing_satisfied &&
692 
693  satisfied = (coarsening_satisfied &&
694  smoothing_satisfied);
695 #ifdef DEBUG
696  bool max_satisfied = satisfied,
697  min_satisfied = satisfied;
698  this->comm().max(max_satisfied);
699  this->comm().min(min_satisfied);
700  libmesh_assert_equal_to (satisfied, max_satisfied);
701  libmesh_assert_equal_to (satisfied, min_satisfied);
702 #endif
703  }
704  while (!satisfied);
705  }
706  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
707 
708  // Coarsen the flagged elements.
709  const bool mesh_changed =
710  this->_coarsen_elements ();
711 
712  if (_maintain_level_one)
714  libmesh_assert(this->make_coarsening_compatible(maintain_level_one));
715 // FIXME: This won't pass unless we add a redundant find_neighbors()
716 // call or replace find_neighbors() with on-the-fly neighbor updating
717 // libmesh_assert(!this->eliminate_unrefined_patches());
718 
719  // We can't contract the mesh ourselves anymore - a System might
720  // need to restrict old coefficient vectors first
721  // _mesh.contract();
722 
723  // Finally, the new mesh may need to be prepared for use
724  if (mesh_changed)
725  _mesh.prepare_for_use (/*skip_renumber =*/false);
726 
727  return mesh_changed;
728 }
729 
730 
731 
732 
733 
734 
735 
736 bool MeshRefinement::refine_elements (const bool maintain_level_one)
737 {
738  // This function must be run on all processors at once
739  parallel_object_only();
740 
741  bool _maintain_level_one = maintain_level_one;
742 
743  // If the user used non-default parameters, let's warn that they're
744  // deprecated
745  if (!maintain_level_one)
746  {
747  libmesh_deprecated();
748  }
749  else
750  _maintain_level_one = _face_level_mismatch_limit;
751 
752  if (_maintain_level_one)
754 
755  // Possibly clean up the refinement flags from
756  // a previous step
758  const MeshBase::element_iterator elem_end = _mesh.elements_end();
759 
760  for ( ; elem_it != elem_end; ++elem_it)
761  {
762  // Pointer to the element
763  Elem *elem = *elem_it;
764 
765  // Set refinement flag to INACTIVE if the
766  // element isn't active
767  if ( !elem->active())
768  {
771  }
772 
773  // This might be left over from the last step
774  if (elem->refinement_flag() == Elem::JUST_REFINED)
776  }
777 
778 
779 
780  // Parallel consistency has to come first, or coarsening
781  // along processor boundaries might occasionally be falsely
782  // prevented
783  bool flags_were_consistent = this->make_flags_parallel_consistent();
784 
785  // In theory, we should be able to remove the above call, which can
786  // be expensive and should be unnecessary. In practice, doing
787  // consistent flagging in parallel is hard, it's impossible to
788  // verify at the library level if it's being done by user code, and
789  // we don't want to abort large parallel runs in opt mode... but we
790  // do want to warn that they should be fixed.
791  libmesh_assert(flags_were_consistent);
792  if (!flags_were_consistent)
793  {
794  libMesh::out << "Refinement flags were not consistent between processors!\n"
795  << "Correcting and continuing.";
796  }
797 
798  // Repeat until flag changes match on every processor
799  do
800  {
801  // Repeat until coarsening & refinement flags jive
802  bool satisfied = false;
803  do
804  {
805  const bool refinement_satisfied =
806  this->make_refinement_compatible(maintain_level_one);
807 
808  bool smoothing_satisfied =
809  !this->eliminate_unrefined_patches();// &&
810 
812  smoothing_satisfied = smoothing_satisfied &&
814 
816  smoothing_satisfied = smoothing_satisfied &&
818 
819  satisfied = (refinement_satisfied &&
820  smoothing_satisfied);
821 #ifdef DEBUG
822  bool max_satisfied = satisfied,
823  min_satisfied = satisfied;
824  this->comm().max(max_satisfied);
825  this->comm().min(min_satisfied);
826  libmesh_assert_equal_to (satisfied, max_satisfied);
827  libmesh_assert_equal_to (satisfied, min_satisfied);
828 #endif
829  }
830  while (!satisfied);
831  }
832  while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
833 
834  // Now refine the flagged elements. This will
835  // take up some space, maybe more than what was freed.
836  const bool mesh_changed =
837  this->_refine_elements();
838 
839  if (_maintain_level_one)
841  libmesh_assert(this->make_refinement_compatible(maintain_level_one));
842 // FIXME: This won't pass unless we add a redundant find_neighbors()
843 // call or replace find_neighbors() with on-the-fly neighbor updating
844 // libmesh_assert(!this->eliminate_unrefined_patches());
845 
846  // Finally, the new mesh needs to be prepared for use
847  if (mesh_changed)
848  _mesh.prepare_for_use (/*skip_renumber =*/false);
849 
850  return mesh_changed;
851 }
852 
853 
854 // Functor for make_flags_parallel_consistent
855 namespace {
856 
857 struct SyncRefinementFlags
858 {
859 typedef unsigned char datum;
860 typedef Elem::RefinementState (Elem::*get_a_flag)() const;
861 typedef void (Elem::*set_a_flag)(const Elem::RefinementState);
862 
863 SyncRefinementFlags(MeshBase &_mesh,
864  get_a_flag _getter,
865  set_a_flag _setter) :
866  mesh(_mesh), parallel_consistent(true),
867  get_flag(_getter), set_flag(_setter) {}
868 
869 MeshBase &mesh;
871 get_a_flag get_flag;
872 set_a_flag set_flag;
873 // References to pointers to member functions segfault?
874 // get_a_flag& get_flag;
875 // set_a_flag& set_flag;
876 
877 // Find the refinement flag on each requested element
878 void gather_data (const std::vector<dof_id_type>& ids,
879  std::vector<datum>& flags)
880 {
881  flags.resize(ids.size());
882 
883  for (std::size_t i=0; i != ids.size(); ++i)
884  {
885  // Look for this element in the mesh
886  // We'd better find every element we're asked for
887  Elem *elem = mesh.elem(ids[i]);
888 
889  // Return the element's refinement flag
890  flags[i] = (elem->*get_flag)();
891  }
892 }
893 
894 void act_on_data (const std::vector<dof_id_type>& ids,
895  std::vector<datum>& flags)
896 {
897  for (std::size_t i=0; i != ids.size(); ++i)
898  {
899  Elem *elem = mesh.elem(ids[i]);
900 
901  datum old_flag = (elem->*get_flag)();
902  datum &new_flag = flags[i];
903 
904  if (old_flag != new_flag)
905  {
906  // It's possible for foreign flags to be (temporarily) more
907  // conservative than our own, such as when a refinement in
908  // one of the foreign processor's elements is mandated by a
909  // refinement in one of our neighboring elements it can see
910  // which was mandated by a refinement in one of our
911  // neighboring elements it can't see
912  // libmesh_assert (!(new_flag != Elem::REFINE &&
913  // old_flag == Elem::REFINE));
914  //
915  (elem->*set_flag)
916  (static_cast<Elem::RefinementState>(new_flag));
917  parallel_consistent = false;
918  }
919  }
920 }
921 };
922 }
923 
924 
925 
927 {
928  // This function must be run on all processors at once
929  parallel_object_only();
930 
931  START_LOG ("make_flags_parallel_consistent()", "MeshRefinement");
932 
933  SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
936  (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
937 
938  SyncRefinementFlags psync(_mesh, &Elem::p_refinement_flag,
941  (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), psync);
942 
943  // If we weren't consistent in both h and p on every processor then
944  // we weren't globally consistent
945  bool parallel_consistent = hsync.parallel_consistent &&
946  psync.parallel_consistent;
947  this->comm().min(parallel_consistent);
948 
949  STOP_LOG ("make_flags_parallel_consistent()", "MeshRefinement");
950 
951  return parallel_consistent;
952 }
953 
954 
955 
956 bool MeshRefinement::make_coarsening_compatible(const bool maintain_level_one)
957 {
958  // This function must be run on all processors at once
959  parallel_object_only();
960 
961  // We may need a PointLocator for topological_neighbor() tests
962  // later, which we need to make sure gets constructed on all
963  // processors at once.
964  AutoPtr<PointLocatorBase> point_locator;
965 
966 #ifdef LIBMESH_ENABLE_PERIODIC
967  bool has_periodic_boundaries =
969  libmesh_assert(this->comm().verify(has_periodic_boundaries));
970 
971  if (has_periodic_boundaries)
972  point_locator = _mesh.sub_point_locator();
973 #endif
974 
975  START_LOG ("make_coarsening_compatible()", "MeshRefinement");
976 
977  bool _maintain_level_one = maintain_level_one;
978 
979  // If the user used non-default parameters, let's warn that they're
980  // deprecated
981  if (!maintain_level_one)
982  {
983  libmesh_deprecated();
984  }
985  else
986  _maintain_level_one = _face_level_mismatch_limit;
987 
988 
989  // Unless we encounter a specific situation level-one
990  // will be satisfied after executing this loop just once
991  bool level_one_satisfied = true;
992 
993 
994  // Unless we encounter a specific situation we will be compatible
995  // with any selected refinement flags
996  bool compatible_with_refinement = true;
997 
998 
999  // find the maximum h and p levels in the mesh
1000  unsigned int max_level = 0;
1001  unsigned int max_p_level = 0;
1002 
1003  {
1004  // First we look at all the active level-0 elements. Since it doesn't make
1005  // sense to coarsen them we must un-set their coarsen flags if
1006  // they are set.
1009 
1010  for (; el != end_el; ++el)
1011  {
1012  Elem *elem = *el;
1013  max_level = std::max(max_level, elem->level());
1014  max_p_level =
1015  std::max(max_p_level,
1016  static_cast<unsigned int>(elem->p_level()));
1017 
1018  if ((elem->level() == 0) &&
1019  (elem->refinement_flag() == Elem::COARSEN))
1021 
1022  if ((elem->p_level() == 0) &&
1023  (elem->p_refinement_flag() == Elem::COARSEN))
1025  }
1026  }
1027 
1028  // if there are no refined elements on this processor then
1029  // there is no work for us to do
1030  if (max_level == 0 && max_p_level == 0)
1031  {
1032  STOP_LOG ("make_coarsening_compatible()", "MeshRefinement");
1033 
1034  // But we still have to check with other processors
1035  this->comm().min(compatible_with_refinement);
1036 
1037  return compatible_with_refinement;
1038  }
1039 
1040 
1041 
1042  // Loop over all the active elements. If an element is marked
1043  // for coarsening we better check its neighbors. If ANY of these neighbors
1044  // are marked for refinement AND are at the same level then there is a
1045  // conflict. By convention refinement wins, so we un-mark the element for
1046  // coarsening. Level-one would be violated in this case so we need to re-run
1047  // the loop.
1048  if (_maintain_level_one)
1049  {
1050 
1051  repeat:
1052  level_one_satisfied = true;
1053 
1054  do
1055  {
1056  level_one_satisfied = true;
1057 
1060 
1061  for (; el != end_el; ++el)
1062  {
1063  Elem* elem = *el;
1064  bool my_flag_changed = false;
1065 
1066  if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
1067  // the coarsen flag is set
1068  {
1069  const unsigned int my_level = elem->level();
1070 
1071  for (unsigned int n=0; n<elem->n_neighbors(); n++)
1072  {
1073  const Elem* neighbor =
1074  topological_neighbor(elem, point_locator.get(), n);
1075 
1076  if (neighbor != NULL && // I have a
1077  neighbor != remote_elem) // neighbor here
1078  {
1079  if (neighbor->active()) // and it is active
1080  {
1081  if ((neighbor->level() == my_level) &&
1082  (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
1083  // and wants to be refined
1084  {
1086  my_flag_changed = true;
1087  break;
1088  }
1089  }
1090  else // I have a neighbor and it is not active. That means it has children.
1091  { // While it _may_ be possible to coarsen us if all the children of
1092  // that element want to be coarsened, it is impossible to know at this
1093  // stage. Forget about it for the moment... This can be handled in
1094  // two steps.
1096  my_flag_changed = true;
1097  break;
1098  }
1099  }
1100  }
1101  }
1102  if (elem->p_refinement_flag() == Elem::COARSEN) // If
1103  // the element is active and the order reduction flag is set
1104  {
1105  const unsigned int my_p_level = elem->p_level();
1106 
1107  for (unsigned int n=0; n<elem->n_neighbors(); n++)
1108  {
1109  const Elem* neighbor =
1110  topological_neighbor(elem, point_locator.get(), n);
1111 
1112  if (neighbor != NULL && // I have a
1113  neighbor != remote_elem) // neighbor here
1114  {
1115  if (neighbor->active()) // and it is active
1116  {
1117  if ((neighbor->p_level() > my_p_level &&
1118  neighbor->p_refinement_flag() != Elem::COARSEN)
1119  || (neighbor->p_level() == my_p_level &&
1120  neighbor->p_refinement_flag() == Elem::REFINE))
1121  {
1123  my_flag_changed = true;
1124  break;
1125  }
1126  }
1127  else // I have a neighbor and it is not active.
1128  { // We need to find which of its children
1129  // have me as a neighbor, and maintain
1130  // level one p compatibility with them.
1131  // Because we currently have level one h
1132  // compatibility, we don't need to check
1133  // grandchildren
1134 
1135  libmesh_assert(neighbor->has_children());
1136  for (unsigned int c=0; c!=neighbor->n_children(); c++)
1137  {
1138  Elem *subneighbor = neighbor->child(c);
1139  if (subneighbor != remote_elem &&
1140  subneighbor->active() &&
1141  has_topological_neighbor(subneighbor, point_locator.get(), elem))
1142  if ((subneighbor->p_level() > my_p_level &&
1143  subneighbor->p_refinement_flag() != Elem::COARSEN)
1144  || (subneighbor->p_level() == my_p_level &&
1145  subneighbor->p_refinement_flag() == Elem::REFINE))
1146  {
1148  my_flag_changed = true;
1149  break;
1150  }
1151  }
1152  if (my_flag_changed)
1153  break;
1154  }
1155  }
1156  }
1157  }
1158 
1159  // If the current element's flag changed, we hadn't
1160  // satisfied the level one rule.
1161  if (my_flag_changed)
1162  level_one_satisfied = false;
1163 
1164  // Additionally, if it has non-local neighbors, and
1165  // we're not in serial, then we'll eventually have to
1166  // return compatible_with_refinement = false, because
1167  // our change has to propagate to neighboring
1168  // processors.
1169  if (my_flag_changed && !_mesh.is_serial())
1170  for (unsigned int n=0; n != elem->n_neighbors(); ++n)
1171  {
1172  Elem* neigh =
1173  topological_neighbor(elem, point_locator.get(), n);
1174 
1175  if (!neigh)
1176  continue;
1177  if (neigh == remote_elem ||
1178  neigh->processor_id() !=
1179  this->processor_id())
1180  {
1181  compatible_with_refinement = false;
1182  break;
1183  }
1184  // FIXME - for non-level one meshes we should
1185  // test all descendants
1186  if (neigh->has_children())
1187  for (unsigned int c=0; c != neigh->n_children(); ++c)
1188  if (neigh->child(c) == remote_elem ||
1189  neigh->child(c)->processor_id() !=
1190  this->processor_id())
1191  {
1192  compatible_with_refinement = false;
1193  break;
1194  }
1195  }
1196  }
1197  }
1198  while (!level_one_satisfied);
1199 
1200  } // end if (_maintain_level_one)
1201 
1202 
1203  // Next we look at all of the ancestor cells.
1204  // If there is a parent cell with all of its children
1205  // wanting to be unrefined then the element is a candidate
1206  // for unrefinement. If all the children don't
1207  // all want to be unrefined then ALL of them need to have their
1208  // unrefinement flags cleared.
1209  for (int level=(max_level); level >= 0; level--)
1210  {
1212  const MeshBase::element_iterator end_el = _mesh.level_elements_end(level);
1213  for (; el != end_el; ++el)
1214  {
1215  Elem *elem = *el;
1216  if (elem->ancestor())
1217  {
1218 
1219  // right now the element hasn't been disqualified
1220  // as a candidate for unrefinement
1221  bool is_a_candidate = true;
1222  bool found_remote_child = false;
1223 
1224  for (unsigned int c=0; c<elem->n_children(); c++)
1225  {
1226  Elem *child = elem->child(c);
1227  if (child == remote_elem)
1228  found_remote_child = true;
1229  else if ((child->refinement_flag() != Elem::COARSEN) ||
1230  !child->active() )
1231  is_a_candidate = false;
1232  }
1233 
1234  if (!is_a_candidate && !found_remote_child)
1235  {
1237 
1238  for (unsigned int c=0; c<elem->n_children(); c++)
1239  {
1240  Elem *child = elem->child(c);
1241  if (child == remote_elem)
1242  continue;
1243  if (child->refinement_flag() == Elem::COARSEN)
1244  {
1245  level_one_satisfied = false;
1247  }
1248  }
1249  }
1250  }
1251  }
1252  }
1253 
1254  if (!level_one_satisfied && _maintain_level_one) goto repeat;
1255 
1256 
1257  // If all the children of a parent are set to be coarsened
1258  // then flag the parent so that they can kill thier kids...
1260  const MeshBase::element_iterator all_el_end = _mesh.elements_end();
1261  for (; all_el != all_el_end; ++all_el)
1262  {
1263  Elem *elem = *all_el;
1264  if (elem->ancestor())
1265  {
1266 
1267  // Presume all the children are local and flagged for
1268  // coarsening and then look for a contradiction
1269  bool all_children_flagged_for_coarsening = true;
1270  bool found_remote_child = false;
1271 
1272  for (unsigned int c=0; c<elem->n_children(); c++)
1273  {
1274  Elem *child = elem->child(c);
1275  if (child == remote_elem)
1276  found_remote_child = true;
1277  else if (child->refinement_flag() != Elem::COARSEN)
1278  all_children_flagged_for_coarsening = false;
1279  }
1280 
1281  if (!found_remote_child &&
1282  all_children_flagged_for_coarsening)
1284  else if (!found_remote_child)
1286  }
1287  }
1288 
1289  STOP_LOG ("make_coarsening_compatible()", "MeshRefinement");
1290 
1291  // If one processor finds an incompatibility, we're globally
1292  // incompatible
1293  this->comm().min(compatible_with_refinement);
1294 
1295  return compatible_with_refinement;
1296 }
1297 
1298 
1299 
1300 
1301 
1302 
1303 
1304 
1305 bool MeshRefinement::make_refinement_compatible(const bool maintain_level_one)
1306 {
1307  // This function must be run on all processors at once
1308  parallel_object_only();
1309 
1310  // We may need a PointLocator for topological_neighbor() tests
1311  // later, which we need to make sure gets constructed on all
1312  // processors at once.
1313  AutoPtr<PointLocatorBase> point_locator;
1314 
1315 #ifdef LIBMESH_ENABLE_PERIODIC
1316  bool has_periodic_boundaries =
1318  libmesh_assert(this->comm().verify(has_periodic_boundaries));
1319 
1320  if (has_periodic_boundaries)
1321  point_locator = _mesh.sub_point_locator();
1322 #endif
1323 
1324  START_LOG ("make_refinement_compatible()", "MeshRefinement");
1325 
1326  bool _maintain_level_one = maintain_level_one;
1327 
1328  // If the user used non-default parameters, let's warn that they're
1329  // deprecated
1330  if (!maintain_level_one)
1331  {
1332  libmesh_deprecated();
1333  }
1334  else
1335  _maintain_level_one = _face_level_mismatch_limit;
1336 
1337  // Unless we encounter a specific situation level-one
1338  // will be satisfied after executing this loop just once
1339  bool level_one_satisfied = true;
1340 
1341  // Unless we encounter a specific situation we will be compatible
1342  // with any selected coarsening flags
1343  bool compatible_with_coarsening = true;
1344 
1345  // This loop enforces the level-1 rule. We should only
1346  // execute it if the user indeed wants level-1 satisfied!
1347  if (_maintain_level_one)
1348  {
1349  do
1350  {
1351  level_one_satisfied = true;
1352 
1355 
1356  for (; el != end_el; ++el)
1357  {
1358  Elem *elem = *el;
1359  if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1360  // h refinement flag is set
1361  {
1362  const unsigned int my_level = elem->level();
1363 
1364  for (unsigned int side=0; side != elem->n_sides(); side++)
1365  {
1366  Elem* neighbor =
1367  topological_neighbor(elem, point_locator.get(), side);
1368 
1369  if (neighbor != NULL && // I have a
1370  neighbor != remote_elem && // neighbor here
1371  neighbor->active()) // and it is active
1372  {
1373  // Case 1: The neighbor is at the same level I am.
1374  // 1a: The neighbor will be refined -> NO PROBLEM
1375  // 1b: The neighbor won't be refined -> NO PROBLEM
1376  // 1c: The neighbor wants to be coarsened -> PROBLEM
1377  if (neighbor->level() == my_level)
1378  {
1379  if (neighbor->refinement_flag() == Elem::COARSEN)
1380  {
1382  if (neighbor->parent())
1384  compatible_with_coarsening = false;
1385  level_one_satisfied = false;
1386  }
1387  }
1388 
1389 
1390  // Case 2: The neighbor is one level lower than I am.
1391  // The neighbor thus MUST be refined to satisfy
1392  // the level-one rule, regardless of whether it
1393  // was originally flagged for refinement. If it
1394  // wasn't flagged already we need to repeat
1395  // this process.
1396  else if ((neighbor->level()+1) == my_level)
1397  {
1398  if (neighbor->refinement_flag() != Elem::REFINE)
1399  {
1400  neighbor->set_refinement_flag(Elem::REFINE);
1401  if (neighbor->parent())
1403  compatible_with_coarsening = false;
1404  level_one_satisfied = false;
1405  }
1406  }
1407 #ifdef DEBUG
1408 
1409  // Sanity check. We should never get into a
1410  // case when our neighbot is more than one
1411  // level away.
1412  else if ((neighbor->level()+1) < my_level)
1413  {
1414  libmesh_error();
1415  }
1416 
1417 
1418  // Note that the only other possibility is that the
1419  // neighbor is already refined, in which case it isn't
1420  // active and we should never get here.
1421  else
1422  {
1423  libmesh_error();
1424  }
1425 #endif
1426  }
1427  }
1428  }
1429  if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1430  // p refinement flag is set
1431  {
1432  const unsigned int my_p_level = elem->p_level();
1433 
1434  for (unsigned int side=0; side != elem->n_sides(); side++)
1435  {
1436  Elem* neighbor =
1437  topological_neighbor(elem, point_locator.get(), side);
1438 
1439  if (neighbor != NULL && // I have a
1440  neighbor != remote_elem) // neighbor here
1441  {
1442  if (neighbor->active()) // and it is active
1443  {
1444  if (neighbor->p_level() < my_p_level &&
1445  neighbor->p_refinement_flag() != Elem::REFINE)
1446  {
1448  level_one_satisfied = false;
1449  compatible_with_coarsening = false;
1450  }
1451  if (neighbor->p_level() == my_p_level &&
1452  neighbor->p_refinement_flag() == Elem::COARSEN)
1453  {
1455  level_one_satisfied = false;
1456  compatible_with_coarsening = false;
1457  }
1458  }
1459  else // I have an inactive neighbor
1460  {
1461  libmesh_assert(neighbor->has_children());
1462  for (unsigned int c=0; c!=neighbor->n_children(); c++)
1463  {
1464  Elem *subneighbor = neighbor->child(c);
1465  if (subneighbor == remote_elem)
1466  continue;
1467  if (subneighbor->active() &&
1468  has_topological_neighbor(subneighbor, point_locator.get(), elem))
1469  {
1470  if (subneighbor->p_level() < my_p_level &&
1471  subneighbor->p_refinement_flag() != Elem::REFINE)
1472  {
1473  // We should already be level one
1474  // compatible
1475  libmesh_assert_greater (subneighbor->p_level() + 2u,
1476  my_p_level);
1477  subneighbor->set_p_refinement_flag(Elem::REFINE);
1478  level_one_satisfied = false;
1479  compatible_with_coarsening = false;
1480  }
1481  if (subneighbor->p_level() == my_p_level &&
1482  subneighbor->p_refinement_flag() == Elem::COARSEN)
1483  {
1485  level_one_satisfied = false;
1486  compatible_with_coarsening = false;
1487  }
1488  }
1489  }
1490  }
1491  }
1492  }
1493  }
1494  }
1495  }
1496 
1497  while (!level_one_satisfied);
1498  } // end if (_maintain_level_one)
1499 
1500  // If we're not compatible on one processor, we're globally not
1501  // compatible
1502  this->comm().min(compatible_with_coarsening);
1503 
1504  STOP_LOG ("make_refinement_compatible()", "MeshRefinement");
1505 
1506  return compatible_with_coarsening;
1507 }
1508 
1509 
1510 
1511 
1513 {
1514  // This function must be run on all processors at once
1515  parallel_object_only();
1516 
1517  START_LOG ("_coarsen_elements()", "MeshRefinement");
1518 
1519  // Flag indicating if this call actually changes the mesh
1520  bool mesh_changed = false;
1521 
1522  // Clear the unused_elements data structure.
1523  // The elements have been packed since it was built,
1524  // so there are _no_ unused elements. We cannot trust
1525  // any iterators currently in this data structure.
1526  // _unused_elements.clear();
1527 
1530 
1531  // Loop over the elements.
1532  for ( ; it != end; ++it)
1533  {
1534  Elem* elem = *it;
1535 
1536  // Not necessary when using elem_iterator
1537  // libmesh_assert(elem);
1538 
1539  // active elements flagged for coarsening will
1540  // no longer be deleted until MeshRefinement::contract()
1541  if (elem->refinement_flag() == Elem::COARSEN)
1542  {
1543  // Huh? no level-0 element should be active
1544  // and flagged for coarsening.
1545  libmesh_assert_not_equal_to (elem->level(), 0);
1546 
1547  // Remove this element from any neighbor
1548  // lists that point to it.
1549  elem->nullify_neighbors();
1550 
1551  // Remove any boundary information associated
1552  // with this element
1553  _mesh.boundary_info->remove (elem);
1554 
1555  // Add this iterator to the _unused_elements
1556  // data structure so we might fill it.
1557  // The _unused_elements optimization is currently off.
1558  // _unused_elements.push_back (it);
1559 
1560  // Don't delete the element until
1561  // MeshRefinement::contract()
1562  // _mesh.delete_elem(elem);
1563 
1564  // the mesh has certainly changed
1565  mesh_changed = true;
1566  }
1567 
1568  // inactive elements flagged for coarsening
1569  // will become active
1570  else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1571  {
1572  elem->coarsen();
1573  libmesh_assert (elem->active());
1574 
1575  // the mesh has certainly changed
1576  mesh_changed = true;
1577  }
1578  if (elem->p_refinement_flag() == Elem::COARSEN)
1579  {
1580  if (elem->p_level() > 0)
1581  {
1583  elem->set_p_level(elem->p_level() - 1);
1584  mesh_changed = true;
1585  }
1586  else
1587  {
1589  }
1590  }
1591  }
1592 
1593  // If the mesh changed on any processor, it changed globally
1594  this->comm().max(mesh_changed);
1595  // And we may need to update ParallelMesh values reflecting the changes
1596  if (mesh_changed)
1598 
1599  // Node processor ids may need to change if an element of that id
1600  // was coarsened away
1601  if (mesh_changed && !_mesh.is_serial())
1602  {
1603  // Update the _new_nodes_map so that processors can
1604  // find requested nodes
1605  this->update_nodes_map ();
1606 
1608  (_mesh, _new_nodes_map);
1609 
1610  // Clear the _new_nodes_map
1611  this->clear();
1612 
1613 #ifdef DEBUG
1615 #endif
1616  }
1617 
1618  STOP_LOG ("_coarsen_elements()", "MeshRefinement");
1619 
1620  return mesh_changed;
1621 }
1622 
1623 
1624 
1626 {
1627  // This function must be run on all processors at once
1628  parallel_object_only();
1629 
1630  // Update the _new_nodes_map so that elements can
1631  // find nodes to connect to.
1632  this->update_nodes_map ();
1633 
1634  START_LOG ("_refine_elements()", "MeshRefinement");
1635 
1636  // Iterate over the elements, counting the elements
1637  // flagged for h refinement.
1638  dof_id_type n_elems_flagged = 0;
1639 
1642 
1643  for (; it != end; ++it)
1644  {
1645  Elem* elem = *it;
1646  if (elem->refinement_flag() == Elem::REFINE)
1647  n_elems_flagged++;
1648  }
1649 
1650  // Construct a local vector of Elem* which have been
1651  // previously marked for refinement. We reserve enough
1652  // space to allow for every element to be refined.
1653  std::vector<Elem*> local_copy_of_elements;
1654  local_copy_of_elements.reserve(n_elems_flagged);
1655 
1656  // Iterate over the elements, looking for elements
1657  // flagged for refinement.
1658  for (it = _mesh.elements_begin(); it != end; ++it)
1659  {
1660  Elem* elem = *it;
1661  if (elem->refinement_flag() == Elem::REFINE)
1662  local_copy_of_elements.push_back(elem);
1663  if (elem->p_refinement_flag() == Elem::REFINE &&
1664  elem->active())
1665  {
1666  elem->set_p_level(elem->p_level()+1);
1668  }
1669  }
1670 
1671  // Now iterate over the local copies and refine each one.
1672  // This may resize the mesh's internal container and invalidate
1673  // any existing iterators.
1674 
1675  for (std::size_t e = 0; e != local_copy_of_elements.size(); ++e)
1676  local_copy_of_elements[e]->refine(*this);
1677 
1678  // The mesh changed if there were elements h refined
1679  bool mesh_changed = !local_copy_of_elements.empty();
1680 
1681  // If the mesh changed on any processor, it changed globally
1682  this->comm().max(mesh_changed);
1683 
1684  // And we may need to update ParallelMesh values reflecting the changes
1685  if (mesh_changed)
1687 
1688  if (mesh_changed && !_mesh.is_serial())
1689  {
1692  (_mesh, _new_nodes_map);
1693 #ifdef DEBUG
1695 #endif
1696  }
1697 
1698  // Clear the _new_nodes_map and _unused_elements data structures.
1699  this->clear();
1700 
1701  STOP_LOG ("_refine_elements()", "MeshRefinement");
1702 
1703  return mesh_changed;
1704 }
1705 
1706 
1707 
1709 {
1710  // Refine n times
1711  for (unsigned int rstep=0; rstep<n; rstep++)
1712  {
1713  // P refine all the active elements
1716 
1717  for ( ; elem_it != elem_end; ++elem_it)
1718  {
1719  (*elem_it)->set_p_level((*elem_it)->p_level()+1);
1720  (*elem_it)->set_p_refinement_flag(Elem::JUST_REFINED);
1721  }
1722  }
1723 }
1724 
1725 
1726 
1728 {
1729  // Coarsen p times
1730  for (unsigned int rstep=0; rstep<n; rstep++)
1731  {
1732  // P coarsen all the active elements
1735 
1736  for ( ; elem_it != elem_end; ++elem_it)
1737  {
1738  if ((*elem_it)->p_level() > 0)
1739  {
1740  (*elem_it)->set_p_level((*elem_it)->p_level()-1);
1741  (*elem_it)->set_p_refinement_flag(Elem::JUST_COARSENED);
1742  }
1743  }
1744  }
1745 }
1746 
1747 
1748 
1750 {
1751  // Refine n times
1752  // FIXME - this won't work if n>1 and the mesh
1753  // has already been attached to an equation system
1754  for (unsigned int rstep=0; rstep<n; rstep++)
1755  {
1756  // Clean up the refinement flags
1757  this->clean_refinement_flags();
1758 
1759  // Flag all the active elements for refinement.
1762 
1763  for ( ; elem_it != elem_end; ++elem_it)
1764  (*elem_it)->set_refinement_flag(Elem::REFINE);
1765 
1766  // Refine all the elements we just flagged.
1767  this->_refine_elements();
1768  }
1769 
1770  // Finally, the new mesh probably needs to be prepared for use
1771  if (n > 0)
1772  _mesh.prepare_for_use (/*skip_renumber =*/false);
1773 }
1774 
1775 
1776 
1778 {
1779  // Coarsen n times
1780  for (unsigned int rstep=0; rstep<n; rstep++)
1781  {
1782  // Clean up the refinement flags
1783  this->clean_refinement_flags();
1784 
1785  // Flag all the active elements for coarsening
1788 
1789  for ( ; elem_it != elem_end; ++elem_it)
1790  {
1791  (*elem_it)->set_refinement_flag(Elem::COARSEN);
1792  if ((*elem_it)->parent())
1793  (*elem_it)->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
1794  }
1795 
1796  // Coarsen all the elements we just flagged.
1797  this->_coarsen_elements();
1798  }
1799 
1800 
1801  // Finally, the new mesh probably needs to be prepared for use
1802  if (n > 0)
1803  _mesh.prepare_for_use (/*skip_renumber =*/false);
1804 }
1805 
1806 
1807 
1809  const PointLocatorBase* point_locator,
1810  const unsigned int side)
1811 {
1812 #ifdef LIBMESH_ENABLE_PERIODIC
1813  if (_periodic_boundaries && !_periodic_boundaries->empty())
1814  {
1815  libmesh_assert(point_locator);
1816  return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1817  }
1818 #endif
1819  return elem->neighbor(side);
1820 }
1821 
1822 
1823 
1825  const PointLocatorBase* point_locator,
1826  Elem* neighbor)
1827 {
1828 #ifdef LIBMESH_ENABLE_PERIODIC
1829  if (_periodic_boundaries && !_periodic_boundaries->empty())
1830  {
1831  libmesh_assert(point_locator);
1832  return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1833  }
1834 #endif
1835  return elem->has_neighbor(neighbor);
1836 }
1837 
1838 
1839 
1840 } // namespace libMesh
1841 
1842 
1843 #endif

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo