libMesh
mesh_refinement_smoothing.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 
20 // Local includes
21 #include "libmesh/libmesh_config.h"
22 
23 // only compile these functions if the user requests AMR support
24 #ifdef LIBMESH_ENABLE_AMR
25 
26 #include "libmesh/elem.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/mesh_refinement.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/remote_elem.h"
31 
32 namespace libMesh
33 {
34 
35 
36 
37 //-----------------------------------------------------------------
38 // Mesh refinement methods
39 bool MeshRefinement::limit_level_mismatch_at_node (const unsigned int max_mismatch)
40 {
41  // This function must be run on all processors at once
42  parallel_object_only();
43 
44  bool flags_changed = false;
45 
46 
47  // Vector holding the maximum element level that touches a node.
48  std::vector<unsigned char> max_level_at_node (_mesh.n_nodes(), 0);
49  std::vector<unsigned char> max_p_level_at_node (_mesh.n_nodes(), 0);
50 
51  // Loop over all the active elements & fill the vector
52  for (auto & elem : _mesh.active_element_ptr_range())
53  {
54  const unsigned char elem_level =
55  cast_int<unsigned char>(elem->level() +
56  ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
57  const unsigned char elem_p_level =
58  cast_int<unsigned char>(elem->p_level() +
59  ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
60 
61  // Set the max_level at each node
62  for (unsigned int n=0; n<elem->n_nodes(); n++)
63  {
64  const dof_id_type node_number = elem->node_id(n);
65 
66  libmesh_assert_less (node_number, max_level_at_node.size());
67 
68  max_level_at_node[node_number] =
69  std::max (max_level_at_node[node_number], elem_level);
70  max_p_level_at_node[node_number] =
71  std::max (max_p_level_at_node[node_number], elem_p_level);
72  }
73  }
74 
75 
76  // Now loop over the active elements and flag the elements
77  // who violate the requested level mismatch. Alternatively, if
78  // _enforce_mismatch_limit_prior_to_refinement is true, swap refinement flags
79  // accordingly.
80  for (auto & elem : _mesh.active_element_ptr_range())
81  {
82  const unsigned int elem_level = elem->level();
83  const unsigned int elem_p_level = elem->p_level();
84 
85  // Skip the element if it is already fully flagged
86  // unless we are enforcing mismatch prior to refinement and may need to
87  // remove the refinement flag(s)
88  if (elem->refinement_flag() == Elem::REFINE &&
89  elem->p_refinement_flag() == Elem::REFINE
91  continue;
92 
93  // Loop over the nodes, check for possible mismatch
94  for (unsigned int n=0; n<elem->n_nodes(); n++)
95  {
96  const dof_id_type node_number = elem->node_id(n);
97 
98  // Flag the element for refinement if it violates
99  // the requested level mismatch
100  if ((elem_level + max_mismatch) < max_level_at_node[node_number]
101  && elem->refinement_flag() != Elem::REFINE)
102  {
103  elem->set_refinement_flag (Elem::REFINE);
104  flags_changed = true;
105  }
106  if ((elem_p_level + max_mismatch) < max_p_level_at_node[node_number]
107  && elem->p_refinement_flag() != Elem::REFINE)
108  {
109  elem->set_p_refinement_flag (Elem::REFINE);
110  flags_changed = true;
111  }
112 
113  // Possibly enforce limit mismatch prior to refinement
114  flags_changed |= this->enforce_mismatch_limit_prior_to_refinement(elem, POINT, max_mismatch);
115  }
116  }
117 
118  // If flags changed on any processor then they changed globally
119  this->comm().max(flags_changed);
120 
121  return flags_changed;
122 }
123 
124 
125 
126 bool MeshRefinement::limit_level_mismatch_at_edge (const unsigned int max_mismatch)
127 {
128  // This function must be run on all processors at once
129  parallel_object_only();
130 
131  bool flags_changed = false;
132 
133 
134  // Maps holding the maximum element level that touches an edge
135  std::map<std::pair<unsigned int, unsigned int>, unsigned char>
136  max_level_at_edge;
137  std::map<std::pair<unsigned int, unsigned int>, unsigned char>
138  max_p_level_at_edge;
139 
140  // Loop over all the active elements & fill the maps
141  for (auto & elem : _mesh.active_element_ptr_range())
142  {
143  const unsigned char elem_level =
144  cast_int<unsigned char>(elem->level() +
145  ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
146  const unsigned char elem_p_level =
147  cast_int<unsigned char>(elem->p_level() +
148  ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
149 
150  // Set the max_level at each edge
151  for (unsigned int n=0; n<elem->n_edges(); n++)
152  {
153  UniquePtr<const Elem> edge = elem->build_edge_ptr(n);
154  dof_id_type childnode0 = edge->node_id(0);
155  dof_id_type childnode1 = edge->node_id(1);
156  if (childnode1 < childnode0)
157  std::swap(childnode0, childnode1);
158 
159  for (const Elem * p = elem; p != libmesh_nullptr; p = p->parent())
160  {
161  UniquePtr<const Elem> pedge = p->build_edge_ptr(n);
162  dof_id_type node0 = pedge->node_id(0);
163  dof_id_type node1 = pedge->node_id(1);
164 
165  if (node1 < node0)
166  std::swap(node0, node1);
167 
168  // If elem does not share this edge with its ancestor
169  // p, refinement levels of elements sharing p's edge
170  // are not restricted by refinement levels of elem.
171  // Furthermore, elem will not share this edge with any
172  // of p's ancestors, so we can safely break out of the
173  // for loop early.
174  if (node0 != childnode0 && node1 != childnode1)
175  break;
176 
177  childnode0 = node0;
178  childnode1 = node1;
179 
180  std::pair<unsigned int, unsigned int> edge_key =
181  std::make_pair(node0, node1);
182 
183  if (max_level_at_edge.find(edge_key) ==
184  max_level_at_edge.end())
185  {
186  max_level_at_edge[edge_key] = elem_level;
187  max_p_level_at_edge[edge_key] = elem_p_level;
188  }
189  else
190  {
191  max_level_at_edge[edge_key] =
192  std::max (max_level_at_edge[edge_key], elem_level);
193  max_p_level_at_edge[edge_key] =
194  std::max (max_p_level_at_edge[edge_key], elem_p_level);
195  }
196  }
197  }
198  }
199 
200 
201  // Now loop over the active elements and flag the elements
202  // who violate the requested level mismatch
203  for (auto & elem : _mesh.active_element_ptr_range())
204  {
205  const unsigned int elem_level = elem->level();
206  const unsigned int elem_p_level = elem->p_level();
207 
208  // Skip the element if it is already fully flagged
209  if (elem->refinement_flag() == Elem::REFINE &&
210  elem->p_refinement_flag() == Elem::REFINE
212  continue;
213 
214  // Loop over the nodes, check for possible mismatch
215  for (unsigned int n=0; n<elem->n_edges(); n++)
216  {
217  UniquePtr<Elem> edge = elem->build_edge_ptr(n);
218  dof_id_type node0 = edge->node_id(0);
219  dof_id_type node1 = edge->node_id(1);
220  if (node1 < node0)
221  std::swap(node0, node1);
222 
223  std::pair<dof_id_type, dof_id_type> edge_key =
224  std::make_pair(node0, node1);
225 
226  // Flag the element for refinement if it violates
227  // the requested level mismatch
228  if ((elem_level + max_mismatch) < max_level_at_edge[edge_key]
229  && elem->refinement_flag() != Elem::REFINE)
230  {
231  elem->set_refinement_flag (Elem::REFINE);
232  flags_changed = true;
233  }
234 
235  if ((elem_p_level + max_mismatch) < max_p_level_at_edge[edge_key]
236  && elem->p_refinement_flag() != Elem::REFINE)
237  {
238  elem->set_p_refinement_flag (Elem::REFINE);
239  flags_changed = true;
240  }
241 
242  // Possibly enforce limit mismatch prior to refinement
243  flags_changed |= this->enforce_mismatch_limit_prior_to_refinement(elem, EDGE, max_mismatch);
244  } // loop over edges
245  } // loop over active elements
246 
247  // If flags changed on any processor then they changed globally
248  this->comm().max(flags_changed);
249 
250  return flags_changed;
251 }
252 
253 
254 
255 bool MeshRefinement::limit_overrefined_boundary(const signed char max_mismatch)
256 {
257  // This function must be run on all processors at once
258  parallel_object_only();
259 
260  bool flags_changed = false;
261 
262  // Loop over all the active elements & look for mismatches to fix.
263  for (auto & elem : _mesh.active_element_ptr_range())
264  {
265  // If we don't have an interior_parent then there's nothing to
266  // be mismatched with.
267  if ((elem->dim() >= LIBMESH_DIM) ||
268  !elem->interior_parent())
269  continue;
270 
271  const unsigned char elem_level =
272  cast_int<unsigned char>(elem->level() +
273  ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0));
274  const unsigned char elem_p_level =
275  cast_int<unsigned char>(elem->p_level() +
276  ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
277 
278  // get all relevant interior elements
279  std::set<const Elem *> neighbor_set;
280  elem->find_interior_neighbors(neighbor_set);
281 
282  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
283  for (; n_it != neighbor_set.end(); ++n_it)
284  {
285  // FIXME - non-const versions of the Elem set methods
286  // would be nice
287  Elem * neighbor = const_cast<Elem *>(*n_it);
288 
289  if (max_mismatch >= 0)
290  {
291  if ((elem_level > neighbor->level() + max_mismatch) &&
292  (neighbor->refinement_flag() != Elem::REFINE))
293  {
295  flags_changed = true;
296  }
297 
298  if ((elem_p_level > neighbor->p_level() + max_mismatch) &&
299  (neighbor->p_refinement_flag() != Elem::REFINE))
300  {
302  flags_changed = true;
303  }
304  }
305  } // loop over interior neighbors
306  }
307 
308  // If flags changed on any processor then they changed globally
309  this->comm().max(flags_changed);
310 
311  return flags_changed;
312 }
313 
314 
315 
316 bool MeshRefinement::limit_underrefined_boundary(const signed char max_mismatch)
317 {
318  // This function must be run on all processors at once
319  parallel_object_only();
320 
321  bool flags_changed = false;
322 
323  // Loop over all the active elements & look for mismatches to fix.
324  for (auto & elem : _mesh.active_element_ptr_range())
325  {
326  // If we don't have an interior_parent then there's nothing to
327  // be mismatched with.
328  if ((elem->dim() >= LIBMESH_DIM) ||
329  !elem->interior_parent())
330  continue;
331 
332  // get all relevant interior elements
333  std::set<const Elem *> neighbor_set;
334  elem->find_interior_neighbors(neighbor_set);
335 
336  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
337  for (; n_it != neighbor_set.end(); ++n_it)
338  {
339  // FIXME - non-const versions of the Elem set methods
340  // would be nice
341  const Elem * neighbor = *n_it;
342 
343  const unsigned char neighbor_level =
344  cast_int<unsigned char>(neighbor->level() +
345  ((neighbor->refinement_flag() == Elem::REFINE) ? 1 : 0));
346 
347  const unsigned char neighbor_p_level =
348  cast_int<unsigned char>(neighbor->p_level() +
349  ((neighbor->p_refinement_flag() == Elem::REFINE) ? 1 : 0));
350 
351  if (max_mismatch >= 0)
352  {
353  if ((neighbor_level >
354  elem->level() + max_mismatch) &&
355  (elem->refinement_flag() != Elem::REFINE))
356  {
358  flags_changed = true;
359  }
360 
361  if ((neighbor_p_level >
362  elem->p_level() + max_mismatch) &&
363  (elem->p_refinement_flag() != Elem::REFINE))
364  {
365  elem->set_p_refinement_flag(Elem::REFINE);
366  flags_changed = true;
367  }
368  }
369  } // loop over interior neighbors
370  }
371 
372  // If flags changed on any processor then they changed globally
373  this->comm().max(flags_changed);
374 
375  return flags_changed;
376 }
377 
378 
379 
381 {
382  // This function must be run on all processors at once
383  parallel_object_only();
384 
385  bool flags_changed = false;
386 
387  // Note: we *cannot* use a reference to the real pointer here, since
388  // the pointer may be reseated below and we don't want to reseat
389  // pointers held by the Mesh.
390  for (Elem * elem : _mesh.active_element_ptr_range())
391  {
392  // First assume that we'll have to flag this element for both h
393  // and p refinement, then change our minds if we see any
394  // neighbors that are as coarse or coarser than us.
395  bool h_flag_me = true,
396  p_flag_me = true;
397 
398 
399  // Skip the element if it is already fully flagged for refinement
400  if (elem->p_refinement_flag() == Elem::REFINE)
401  p_flag_me = false;
402  if (elem->refinement_flag() == Elem::REFINE)
403  {
404  h_flag_me = false;
405  if (!p_flag_me)
406  continue;
407  }
408  // Test the parent if that is already flagged for coarsening
409  else if (elem->refinement_flag() == Elem::COARSEN)
410  {
411  libmesh_assert(elem->parent());
412  elem = elem->parent();
413  // FIXME - this doesn't seem right - RHS
414  if (elem->refinement_flag() != Elem::COARSEN_INACTIVE)
415  continue;
416  p_flag_me = false;
417  }
418 
419  const unsigned int my_level = elem->level();
420  int my_p_adjustment = 0;
421  if (elem->p_refinement_flag() == Elem::REFINE)
422  my_p_adjustment = 1;
423  else if (elem->p_refinement_flag() == Elem::COARSEN)
424  {
425  libmesh_assert_greater (elem->p_level(), 0);
426  my_p_adjustment = -1;
427  }
428  const unsigned int my_new_p_level = elem->p_level() +
429  my_p_adjustment;
430 
431  // Check all the element neighbors
432  for (auto neighbor : elem->neighbor_ptr_range())
433  {
434  // Quit if the element is on a local boundary
435  if (neighbor == libmesh_nullptr || neighbor == remote_elem)
436  {
437  h_flag_me = false;
438  p_flag_me = false;
439  break;
440  }
441  // if the neighbor will be equally or less refined than
442  // we are, then we will not become an unrefined island.
443  // So if we are still considering h refinement:
444  if (h_flag_me &&
445  // If our neighbor is already at a lower level,
446  // it can't end up at a higher level even if it
447  // is flagged for refinement once
448  ((neighbor->level() < my_level) ||
449  // If our neighbor is at the same level but isn't
450  // flagged for refinement, it won't end up at a
451  // higher level
452  ((neighbor->active()) &&
453  (neighbor->refinement_flag() != Elem::REFINE)) ||
454  // If our neighbor is currently more refined but is
455  // a parent flagged for coarsening, it will end up
456  // at the same level.
457  (neighbor->refinement_flag() == Elem::COARSEN_INACTIVE)))
458  {
459  // We've proven we won't become an unrefined island,
460  // so don't h refine to avoid that.
461  h_flag_me = false;
462 
463  // If we've also proven we don't need to p refine,
464  // we don't need to check more neighbors
465  if (!p_flag_me)
466  break;
467  }
468  if (p_flag_me)
469  {
470  // if active neighbors will have a p level
471  // equal to or lower than ours, then we do not need to p
472  // refine ourselves.
473  if (neighbor->active())
474  {
475  int p_adjustment = 0;
476  if (neighbor->p_refinement_flag() == Elem::REFINE)
477  p_adjustment = 1;
478  else if (neighbor->p_refinement_flag() == Elem::COARSEN)
479  {
480  libmesh_assert_greater (neighbor->p_level(), 0);
481  p_adjustment = -1;
482  }
483  if (my_new_p_level >= neighbor->p_level() + p_adjustment)
484  {
485  p_flag_me = false;
486  if (!h_flag_me)
487  break;
488  }
489  }
490  // If we have inactive neighbors, we need to
491  // test all their active descendants which neighbor us
492  else if (neighbor->ancestor())
493  {
494  if (neighbor->min_new_p_level_by_neighbor(elem,
495  my_new_p_level + 2) <= my_new_p_level)
496  {
497  p_flag_me = false;
498  if (!h_flag_me)
499  break;
500  }
501  }
502  }
503  }
504 
505  if (h_flag_me)
506  {
507  // Parents that would create islands should no longer
508  // coarsen
509  if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
510  {
511  for (auto & child : elem->child_ref_range())
512  {
513  libmesh_assert_equal_to (child.refinement_flag(),
514  Elem::COARSEN);
515  child.set_refinement_flag(Elem::DO_NOTHING);
516  }
517  elem->set_refinement_flag(Elem::INACTIVE);
518  }
519  else
520  elem->set_refinement_flag(Elem::REFINE);
521  flags_changed = true;
522  }
523  if (p_flag_me)
524  {
525  if (elem->p_refinement_flag() == Elem::COARSEN)
526  elem->set_p_refinement_flag(Elem::DO_NOTHING);
527  else
528  elem->set_p_refinement_flag(Elem::REFINE);
529  flags_changed = true;
530  }
531  }
532 
533  // If flags changed on any processor then they changed globally
534  this->comm().max(flags_changed);
535 
536  return flags_changed;
537 }
538 
539 
540 
542  NeighborType nt,
543  unsigned max_mismatch)
544 {
545  // Eventual return value
546  bool flags_changed = false;
547 
548  // If we are enforcing the limit prior to refinement then we
549  // need to remove flags from any elements marked for refinement that
550  // would cause a mismatch
552  && elem->refinement_flag() == Elem::REFINE)
553  {
554  // get all the relevant neighbors since we may have to refine
555  // elements off edges or corners as well
556  std::set<const Elem *> neighbor_set;
557 
558  if (nt == POINT)
559  elem->find_point_neighbors(neighbor_set);
560  else if (nt == EDGE)
561  elem->find_edge_neighbors(neighbor_set);
562  else
563  libmesh_error_msg("Unrecognized NeighborType: " << nt);
564 
565  // Loop over the neighbors of element e
566  std::set<const Elem *>::iterator n_it = neighbor_set.begin();
567  for (; n_it != neighbor_set.end(); ++n_it)
568  {
569  const Elem * neighbor = *n_it;
570 
571  if ((elem->level() + 1 - max_mismatch) > neighbor->level())
572  {
574  flags_changed = true;
575  }
576  if ((elem->p_level() + 1 - max_mismatch) > neighbor->p_level())
577  {
579  flags_changed = true;
580  }
581  } // loop over edge/point neighbors
582  } // if _enforce_mismatch_limit_prior_to_refinement
583 
584  return flags_changed;
585 }
586 
587 } // namespace libMesh
588 
589 
590 #endif
bool limit_level_mismatch_at_edge(const unsigned int max_mismatch)
bool & enforce_mismatch_limit_prior_to_refinement()
Get/set the _enforce_mismatch_limit_prior_to_refinement flag.
NeighborType
This helper function enforces the desired mismatch limits prior to refinement.
bool limit_level_mismatch_at_node(const unsigned int max_mismatch)
This algorithm restricts the maximum level mismatch at any node in the mesh.
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.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
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
bool limit_overrefined_boundary(const signed char max_mismatch)
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
The libMesh namespace provides an interface to certain functionality in the library.
void find_point_neighbors(const Point &p, std::set< const Elem * > &neighbor_set) const
This function finds all active elements (including this one) which are in the same manifold as this e...
Definition: elem.C:605
long double max(long double a, double b)
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
void find_edge_neighbors(const Point &p1, const Point &p2, std::set< const Elem * > &neighbor_set) const
This function finds all active elements in the same manifold as this element which touch the current ...
Definition: elem.C:770
RefinementState p_refinement_flag() const
Definition: elem.h:2521
bool _enforce_mismatch_limit_prior_to_refinement
This option enforces the mismatch level prior to refinement by checking if refining any element marke...
const Parallel::Communicator & comm() const
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
unsigned int level() const
Definition: elem.h:2388
RefinementState refinement_flag() const
Definition: elem.h:2505
void set_p_refinement_flag(const RefinementState pflag)
Sets the value of the p-refinement flag for the element.
Definition: elem.h:2529
bool eliminate_unrefined_patches()
This algorithm selects an element for refinement if all of its neighbors are (or will be) refined...
virtual dof_id_type n_nodes() const =0
uint8_t dof_id_type
Definition: id_types.h:64
const RemoteElem * remote_elem
Definition: remote_elem.C:57