mesh_refinement_smoothing.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
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  std::unique_ptr<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 != nullptr; p = p->parent())
160  {
161  std::unique_ptr<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  std::unique_ptr<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, see if there's any possibility we might have to flag
393  // this element for h and p refinement - do we have any visible
394  // neighbors? Next we'll check to see if any of those neighbors
395  // are as coarse or coarser than us.
396  bool h_flag_me = false,
397  p_flag_me = false;
398  for (auto neighbor : elem->neighbor_ptr_range())
399  {
400  // Quit if the element is not a local boundary
401  if (neighbor != nullptr && neighbor != remote_elem)
402  {
403  h_flag_me = true;
404  p_flag_me = true;
405  break;
406  }
407  }
408 
409  // Skip the element if it is already fully flagged for refinement
410  if (elem->p_refinement_flag() == Elem::REFINE)
411  p_flag_me = false;
412  if (elem->refinement_flag() == Elem::REFINE)
413  {
414  h_flag_me = false;
415  if (!p_flag_me)
416  continue;
417  }
418  // Test the parent if that is already flagged for coarsening
419  else if (elem->refinement_flag() == Elem::COARSEN)
420  {
421  libmesh_assert(elem->parent());
422  elem = elem->parent();
423  // FIXME - this doesn't seem right - RHS
424  if (elem->refinement_flag() != Elem::COARSEN_INACTIVE)
425  continue;
426  p_flag_me = false;
427  }
428 
429  const unsigned int my_level = elem->level();
430  int my_p_adjustment = 0;
431  if (elem->p_refinement_flag() == Elem::REFINE)
432  my_p_adjustment = 1;
433  else if (elem->p_refinement_flag() == Elem::COARSEN)
434  {
435  libmesh_assert_greater (elem->p_level(), 0);
436  my_p_adjustment = -1;
437  }
438  const unsigned int my_new_p_level = elem->p_level() +
439  my_p_adjustment;
440 
441  // Check all the element neighbors
442  for (auto neighbor : elem->neighbor_ptr_range())
443  {
444  // Quit if the element is on a local boundary
445  if (neighbor == nullptr || neighbor == remote_elem)
446  {
447  h_flag_me = false;
448  p_flag_me = false;
449  break;
450  }
451  // if the neighbor will be equally or less refined than
452  // we are, then we will not become an unrefined island.
453  // So if we are still considering h refinement:
454  if (h_flag_me &&
455  // If our neighbor is already at a lower level,
456  // it can't end up at a higher level even if it
457  // is flagged for refinement once
458  ((neighbor->level() < my_level) ||
459  // If our neighbor is at the same level but isn't
460  // flagged for refinement, it won't end up at a
461  // higher level
462  ((neighbor->active()) &&
463  (neighbor->refinement_flag() != Elem::REFINE)) ||
464  // If our neighbor is currently more refined but is
465  // a parent flagged for coarsening, it will end up
466  // at the same level.
467  (neighbor->refinement_flag() == Elem::COARSEN_INACTIVE)))
468  {
469  // We've proven we won't become an unrefined island,
470  // so don't h refine to avoid that.
471  h_flag_me = false;
472 
473  // If we've also proven we don't need to p refine,
474  // we don't need to check more neighbors
475  if (!p_flag_me)
476  break;
477  }
478  if (p_flag_me)
479  {
480  // if active neighbors will have a p level
481  // equal to or lower than ours, then we do not need to p
482  // refine ourselves.
483  if (neighbor->active())
484  {
485  int p_adjustment = 0;
486  if (neighbor->p_refinement_flag() == Elem::REFINE)
487  p_adjustment = 1;
488  else if (neighbor->p_refinement_flag() == Elem::COARSEN)
489  {
490  libmesh_assert_greater (neighbor->p_level(), 0);
491  p_adjustment = -1;
492  }
493  if (my_new_p_level >= neighbor->p_level() + p_adjustment)
494  {
495  p_flag_me = false;
496  if (!h_flag_me)
497  break;
498  }
499  }
500  // If we have inactive neighbors, we need to
501  // test all their active descendants which neighbor us
502  else if (neighbor->ancestor())
503  {
504  if (neighbor->min_new_p_level_by_neighbor(elem,
505  my_new_p_level + 2) <= my_new_p_level)
506  {
507  p_flag_me = false;
508  if (!h_flag_me)
509  break;
510  }
511  }
512  }
513  }
514 
515  if (h_flag_me)
516  {
517  // Parents that would create islands should no longer
518  // coarsen
519  if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
520  {
521  for (auto & child : elem->child_ref_range())
522  {
523  libmesh_assert_equal_to (child.refinement_flag(),
524  Elem::COARSEN);
525  child.set_refinement_flag(Elem::DO_NOTHING);
526  }
527  elem->set_refinement_flag(Elem::INACTIVE);
528  }
529  else
530  elem->set_refinement_flag(Elem::REFINE);
531  flags_changed = true;
532  }
533  if (p_flag_me)
534  {
535  if (elem->p_refinement_flag() == Elem::COARSEN)
536  elem->set_p_refinement_flag(Elem::DO_NOTHING);
537  else
538  elem->set_p_refinement_flag(Elem::REFINE);
539  flags_changed = true;
540  }
541  }
542 
543  // If flags changed on any processor then they changed globally
544  this->comm().max(flags_changed);
545 
546  return flags_changed;
547 }
548 
549 
550 
552  NeighborType nt,
553  unsigned max_mismatch)
554 {
555  // Eventual return value
556  bool flags_changed = false;
557 
558  // If we are enforcing the limit prior to refinement then we
559  // need to remove flags from any elements marked for refinement that
560  // would cause a mismatch
562  && elem->refinement_flag() == Elem::REFINE)
563  {
564  // get all the relevant neighbors since we may have to refine
565  // elements off edges or corners as well
566  std::set<const Elem *> neighbor_set;
567 
568  if (nt == POINT)
569  elem->find_point_neighbors(neighbor_set);
570  else if (nt == EDGE)
571  elem->find_edge_neighbors(neighbor_set);
572  else
573  libmesh_error_msg("Unrecognized NeighborType: " << nt);
574 
575  // Loop over the neighbors of element e
576  for (const auto & neighbor : neighbor_set)
577  {
578  if ((elem->level() + 1 - max_mismatch) > neighbor->level())
579  {
581  flags_changed = true;
582  }
583  if ((elem->p_level() + 1 - max_mismatch) > neighbor->p_level())
584  {
586  flags_changed = true;
587  }
588  } // loop over edge/point neighbors
589  } // if _enforce_mismatch_limit_prior_to_refinement
590 
591  return flags_changed;
592 }
593 
594 } // namespace libMesh
595 
596 
597 #endif
bool limit_level_mismatch_at_edge(const unsigned int max_mismatch)
bool & enforce_mismatch_limit_prior_to_refinement()
RefinementState refinement_flag() const
Definition: elem.h:2638
const Elem * parent() const
Definition: elem.h:2479
bool limit_level_mismatch_at_node(const unsigned int max_mismatch)
bool limit_underrefined_boundary(const signed char max_mismatch)
RefinementState p_refinement_flag() const
Definition: elem.h:2654
The base class for all geometric element types.
Definition: elem.h:100
void set_refinement_flag(const RefinementState rflag)
Definition: elem.h:2646
const Parallel::Communicator & comm() const
bool limit_overrefined_boundary(const signed char max_mismatch)
unsigned int p_level() const
Definition: elem.h:2555
void find_edge_neighbors(const Point &p1, const Point &p2, std::set< const Elem *> &neighbor_set) const
Definition: elem.C:637
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
long double max(long double a, double b)
void find_point_neighbors(const Point &p, std::set< const Elem *> &neighbor_set) const
Definition: elem.C:498
bool _enforce_mismatch_limit_prior_to_refinement
unsigned int level() const
Definition: elem.h:2521
void swap(Iterator &lhs, Iterator &rhs)
void set_p_refinement_flag(const RefinementState pflag)
Definition: elem.h:2662
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