mesh_refinement_smoothing.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 
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 
52  // Loop over all the active elements & fill the vector
53  {
56 
57  for (; elem_it != elem_end; ++elem_it)
58  {
59  const Elem* elem = *elem_it;
60  const unsigned char elem_level =
61  elem->level() + ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0);
62  const unsigned char elem_p_level =
63  elem->p_level() + ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0);
64 
65  // Set the max_level at each node
66  for (unsigned int n=0; n<elem->n_nodes(); n++)
67  {
68  const dof_id_type node_number = elem->node(n);
69 
70  libmesh_assert_less (node_number, max_level_at_node.size());
71 
72  max_level_at_node[node_number] =
73  std::max (max_level_at_node[node_number], elem_level);
74  max_p_level_at_node[node_number] =
75  std::max (max_p_level_at_node[node_number], elem_p_level);
76  }
77  }
78  }
79 
80 
81  // Now loop over the active elements and flag the elements
82  // who violate the requested level mismatch
83  {
86 
87  for (; elem_it != elem_end; ++elem_it)
88  {
89  Elem* elem = *elem_it;
90  const unsigned int elem_level = elem->level();
91  const unsigned int elem_p_level = elem->p_level();
92 
93  // Skip the element if it is already fully flagged
94  if (elem->refinement_flag() == Elem::REFINE &&
96  continue;
97 
98  // Loop over the nodes, check for possible mismatch
99  for (unsigned int n=0; n<elem->n_nodes(); n++)
100  {
101  const dof_id_type node_number = elem->node(n);
102 
103  // Flag the element for refinement if it violates
104  // the requested level mismatch
105  if ( (elem_level + max_mismatch) < max_level_at_node[node_number]
106  && elem->refinement_flag() != Elem::REFINE)
107  {
109  flags_changed = true;
110  }
111  if ( (elem_p_level + max_mismatch) < max_p_level_at_node[node_number]
112  && elem->p_refinement_flag() != Elem::REFINE)
113  {
115  flags_changed = true;
116  }
117  }
118  }
119  }
120 
121  // If flags changed on any processor then they changed globally
122  this->comm().max(flags_changed);
123 
124  return flags_changed;
125 }
126 
127 
128 
129 //-----------------------------------------------------------------
130 // Mesh refinement methods
131 bool MeshRefinement::limit_level_mismatch_at_edge (const unsigned int max_mismatch)
132 {
133  // This function must be run on all processors at once
134  parallel_object_only();
135 
136  bool flags_changed = false;
137 
138 
139  // Maps holding the maximum element level that touches an edge
140  std::map<std::pair<unsigned int, unsigned int>, unsigned char>
141  max_level_at_edge;
142  std::map<std::pair<unsigned int, unsigned int>, unsigned char>
143  max_p_level_at_edge;
144 
145  // Loop over all the active elements & fill the maps
146  {
149 
150  for (; elem_it != elem_end; ++elem_it)
151  {
152  const Elem* elem = *elem_it;
153  const unsigned char elem_level =
154  elem->level() + ((elem->refinement_flag() == Elem::REFINE) ? 1 : 0);
155  const unsigned char elem_p_level =
156  elem->p_level() + ((elem->p_refinement_flag() == Elem::REFINE) ? 1 : 0);
157 
158  // Set the max_level at each edge
159  for (unsigned int n=0; n<elem->n_edges(); n++)
160  {
161  AutoPtr<Elem> edge = elem->build_edge(n);
162  dof_id_type childnode0 = edge->node(0);
163  dof_id_type childnode1 = edge->node(1);
164  if (childnode1 < childnode0)
165  std::swap(childnode0, childnode1);
166 
167  for (const Elem *p = elem; p != NULL; p = p->parent())
168  {
169  AutoPtr<Elem> pedge = p->build_edge(n);
170  dof_id_type node0 = pedge->node(0);
171  dof_id_type node1 = pedge->node(1);
172 
173  if (node1 < node0)
174  std::swap(node0, node1);
175 
176  // If elem does not share this edge with its ancestor
177  // p, refinement levels of elements sharing p's edge
178  // are not restricted by refinement levels of elem.
179  // Furthermore, elem will not share this edge with any
180  // of p's ancestors, so we can safely break out of the
181  // for loop early.
182  if (node0 != childnode0 && node1 != childnode1)
183  break;
184 
185  childnode0 = node0;
186  childnode1 = node1;
187 
188  std::pair<unsigned int, unsigned int> edge_key =
189  std::make_pair(node0, node1);
190 
191  if (max_level_at_edge.find(edge_key) ==
192  max_level_at_edge.end())
193  {
194  max_level_at_edge[edge_key] = elem_level;
195  max_p_level_at_edge[edge_key] = elem_p_level;
196  }
197  else
198  {
199  max_level_at_edge[edge_key] =
200  std::max (max_level_at_edge[edge_key], elem_level);
201  max_p_level_at_edge[edge_key] =
202  std::max (max_p_level_at_edge[edge_key], elem_p_level);
203  }
204  }
205  }
206  }
207  }
208 
209 
210  // Now loop over the active elements and flag the elements
211  // who violate the requested level mismatch
212  {
215 
216  for (; elem_it != elem_end; ++elem_it)
217  {
218  Elem* elem = *elem_it;
219  const unsigned int elem_level = elem->level();
220  const unsigned int elem_p_level = elem->p_level();
221 
222  // Skip the element if it is already fully flagged
223  if (elem->refinement_flag() == Elem::REFINE &&
224  elem->p_refinement_flag() == Elem::REFINE)
225  continue;
226 
227  // Loop over the nodes, check for possible mismatch
228  for (unsigned int n=0; n<elem->n_edges(); n++)
229  {
230  AutoPtr<Elem> edge = elem->build_edge(n);
231  dof_id_type node0 = edge->node(0);
232  dof_id_type node1 = edge->node(1);
233  if (node1 < node0)
234  std::swap(node0, node1);
235 
236  std::pair<dof_id_type, dof_id_type> edge_key =
237  std::make_pair(node0, node1);
238 
239  // Flag the element for refinement if it violates
240  // the requested level mismatch
241  if ( (elem_level + max_mismatch) < max_level_at_edge[edge_key]
242  && elem->refinement_flag() != Elem::REFINE)
243  {
245  flags_changed = true;
246  }
247  if ( (elem_p_level + max_mismatch) < max_p_level_at_edge[edge_key]
248  && elem->p_refinement_flag() != Elem::REFINE)
249  {
251  flags_changed = true;
252  }
253  }
254  }
255  }
256 
257  // If flags changed on any processor then they changed globally
258  this->comm().max(flags_changed);
259 
260  return flags_changed;
261 }
262 
263 
264 
265 
267 {
268  // This function must be run on all processors at once
269  parallel_object_only();
270 
271  bool flags_changed = false;
272 
275 
276  for (; elem_it != elem_end; ++elem_it)
277  {
278  Elem* elem = *elem_it;
279  // First assume that we'll have to flag this element for both h
280  // and p refinement, then change our minds if we see any
281  // neighbors that are as coarse or coarser than us.
282  bool h_flag_me = true,
283  p_flag_me = true;
284 
285 
286  // Skip the element if it is already fully flagged for refinement
287  if (elem->p_refinement_flag() == Elem::REFINE)
288  p_flag_me = false;
289  if (elem->refinement_flag() == Elem::REFINE)
290  {
291  h_flag_me = false;
292  if (!p_flag_me)
293  continue;
294  }
295  // Test the parent if that is already flagged for coarsening
296  else if (elem->refinement_flag() == Elem::COARSEN)
297  {
298  libmesh_assert(elem->parent());
299  elem = elem->parent();
300  // FIXME - this doesn't seem right - RHS
301  if (elem->refinement_flag() != Elem::COARSEN_INACTIVE)
302  continue;
303  p_flag_me = false;
304  }
305 
306  const unsigned int my_level = elem->level();
307  int my_p_adjustment = 0;
308  if (elem->p_refinement_flag() == Elem::REFINE)
309  my_p_adjustment = 1;
310  else if (elem->p_refinement_flag() == Elem::COARSEN)
311  {
312  libmesh_assert_greater (elem->p_level(), 0);
313  my_p_adjustment = -1;
314  }
315  const unsigned int my_new_p_level = elem->p_level() +
316  my_p_adjustment;
317 
318  // Check all the element neighbors
319  for (unsigned int n=0; n<elem->n_neighbors(); n++)
320  {
321  const Elem *neighbor = elem->neighbor(n);
322  // Quit if the element is on a local boundary
323  if (neighbor == NULL || neighbor == remote_elem)
324  {
325  h_flag_me = false;
326  p_flag_me = false;
327  break;
328  }
329  // if the neighbor will be equally or less refined than
330  // we are, then we will not become an unrefined island.
331  // So if we are still considering h refinement:
332  if (h_flag_me &&
333  // If our neighbor is already at a lower level,
334  // it can't end up at a higher level even if it
335  // is flagged for refinement once
336  ((neighbor->level() < my_level) ||
337  // If our neighbor is at the same level but isn't
338  // flagged for refinement, it won't end up at a
339  // higher level
340  ((neighbor->active()) &&
341  (neighbor->refinement_flag() != Elem::REFINE)) ||
342  // If our neighbor is currently more refined but is
343  // a parent flagged for coarsening, it will end up
344  // at the same level.
345  (neighbor->refinement_flag() == Elem::COARSEN_INACTIVE)))
346  {
347  // We've proven we won't become an unrefined island,
348  // so don't h refine to avoid that.
349  h_flag_me = false;
350 
351  // If we've also proven we don't need to p refine,
352  // we don't need to check more neighbors
353  if (!p_flag_me)
354  break;
355  }
356  if (p_flag_me)
357  {
358  // if active neighbors will have a p level
359  // equal to or lower than ours, then we do not need to p
360  // refine ourselves.
361  if (neighbor->active())
362  {
363  int p_adjustment = 0;
364  if (neighbor->p_refinement_flag() == Elem::REFINE)
365  p_adjustment = 1;
366  else if (neighbor->p_refinement_flag() == Elem::COARSEN)
367  {
368  libmesh_assert_greater (neighbor->p_level(), 0);
369  p_adjustment = -1;
370  }
371  if (my_new_p_level >= neighbor->p_level() + p_adjustment)
372  {
373  p_flag_me = false;
374  if (!h_flag_me)
375  break;
376  }
377  }
378  // If we have inactive neighbors, we need to
379  // test all their active descendants which neighbor us
380  else if (neighbor->ancestor())
381  {
382  if (neighbor->min_new_p_level_by_neighbor(elem,
383  my_new_p_level + 2) <= my_new_p_level)
384  {
385  p_flag_me = false;
386  if (!h_flag_me)
387  break;
388  }
389  }
390  }
391  }
392 
393  if (h_flag_me)
394  {
395  // Parents that would create islands should no longer
396  // coarsen
398  {
399  for (unsigned int c=0; c<elem->n_children(); c++)
400  {
401  libmesh_assert_equal_to (elem->child(c)->refinement_flag(),
402  Elem::COARSEN);
404  }
406  }
407  else
409  flags_changed = true;
410  }
411  if (p_flag_me)
412  {
413  if (elem->p_refinement_flag() == Elem::COARSEN)
415  else
417  flags_changed = true;
418  }
419  }
420 
421  // If flags changed on any processor then they changed globally
422  this->comm().max(flags_changed);
423 
424  return flags_changed;
425 }
426 
427 } // namespace libMesh
428 
429 
430 #endif

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

Hosted By:
SourceForge.net Logo