mesh_refinement_flagging.C
Go to the documentation of this file.
1 
2 // The libMesh Finite Element Library.
3 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
4 
5 // This library is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU Lesser General Public
7 // License as published by the Free Software Foundation; either
8 // version 2.1 of the License, or (at your option) any later version.
9 
10 // This library is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
14 
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with this library; if not, write to the Free Software
17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 
19 
20 
21 // Local includes
22 #include "libmesh/libmesh_config.h"
23 
24 // only compile these functions if the user requests AMR support
25 #ifdef LIBMESH_ENABLE_AMR
26 
27 // C++ includes
28 #include <algorithm> // for std::sort
29 
30 // Local includes
31 #include "libmesh/elem.h"
32 #include "libmesh/error_vector.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/parallel.h"
36 #include "libmesh/remote_elem.h"
37 
38 namespace libMesh
39 {
40 
41 
42 
43 //-----------------------------------------------------------------
44 // Mesh refinement methods
46  const Real refine_frac,
47  const Real coarsen_frac,
48  const unsigned int max_l)
49 {
50  parallel_object_only();
51 
52  // The function arguments are currently just there for
53  // backwards_compatibility
55  {
56  // If the user used non-default parameters, lets warn
57  // that they're deprecated
58  if (refine_frac != 0.3 ||
59  coarsen_frac != 0.0 ||
60  max_l != libMesh::invalid_uint)
61  libmesh_deprecated();
62 
63  _refine_fraction = refine_frac;
64  _coarsen_fraction = coarsen_frac;
65  _max_h_level = max_l;
66  }
67 
68  // Check for valid fractions..
69  // The fraction values must be in [0,1]
70  libmesh_assert_greater_equal (_refine_fraction, 0);
71  libmesh_assert_less_equal (_refine_fraction, 1);
72  libmesh_assert_greater_equal (_coarsen_fraction, 0);
73  libmesh_assert_less_equal (_coarsen_fraction, 1);
74 
75  // Clean up the refinement flags. These could be left
76  // over from previous refinement steps.
77  this->clean_refinement_flags();
78 
79  // We're getting the minimum and maximum error values
80  // for the ACTIVE elements
81  Real error_min = 1.e30;
82  Real error_max = 0.;
83 
84  // And, if necessary, for their parents
85  Real parent_error_min = 1.e30;
86  Real parent_error_max = 0.;
87 
88  // Prepare another error vector if we need to sum parent errors
89  ErrorVector error_per_parent;
91  {
92  create_parent_error_vector(error_per_cell,
93  error_per_parent,
94  parent_error_min,
95  parent_error_max);
96  }
97 
98  // We need to loop over all active elements to find the minimum
101  const MeshBase::element_iterator el_end =
103 
104  for (; el_it != el_end; ++el_it)
105  {
106  const dof_id_type id = (*el_it)->id();
107  libmesh_assert_less (id, error_per_cell.size());
108 
109  error_max = std::max (error_max, error_per_cell[id]);
110  error_min = std::min (error_min, error_per_cell[id]);
111  }
112  this->comm().max(error_max);
113  this->comm().min(error_min);
114 
115  // Compute the cutoff values for coarsening and refinement
116  const Real error_delta = (error_max - error_min);
117  const Real parent_error_delta = parent_error_max - parent_error_min;
118 
119  const Real refine_cutoff = (1.- _refine_fraction)*error_max;
120  const Real coarsen_cutoff = _coarsen_fraction*error_delta + error_min;
121  const Real parent_cutoff = _coarsen_fraction*parent_error_delta + error_min;
122 
123 // // Print information about the error
124 // libMesh::out << " Error Information:" << std::endl
125 // << " ------------------" << std::endl
126 // << " min: " << error_min << std::endl
127 // << " max: " << error_max << std::endl
128 // << " delta: " << error_delta << std::endl
129 // << " refine_cutoff: " << refine_cutoff << std::endl
130 // << " coarsen_cutoff: " << coarsen_cutoff << std::endl;
131 
132 
133 
134  // Loop over the elements and flag them for coarsening or
135  // refinement based on the element error
136 
139  const MeshBase::element_iterator e_end =
141  for (; e_it != e_end; ++e_it)
142  {
143  Elem* elem = *e_it;
144  const dof_id_type id = elem->id();
145 
146  libmesh_assert_less (id, error_per_cell.size());
147 
148  const ErrorVectorReal elem_error = error_per_cell[id];
149 
151  {
152  Elem* parent = elem->parent();
153  if (parent)
154  {
155  const dof_id_type parentid = parent->id();
156  if (error_per_parent[parentid] >= 0. &&
157  error_per_parent[parentid] <= parent_cutoff)
159  }
160  }
161  // Flag the element for coarsening if its error
162  // is <= coarsen_fraction*delta + error_min
163  else if (elem_error <= coarsen_cutoff)
164  {
166  }
167 
168  // Flag the element for refinement if its error
169  // is >= refinement_cutoff.
170  if (elem_error >= refine_cutoff)
171  if (elem->level() < _max_h_level)
173  }
174 }
175 
176 
177 
179 {
180  parallel_object_only();
181 
183 
184  // Check for valid fractions..
185  // The fraction values must be in [0,1]
186  libmesh_assert_greater_equal (_refine_fraction, 0);
187  libmesh_assert_less_equal (_refine_fraction, 1);
188  libmesh_assert_greater_equal (_coarsen_fraction, 0);
189  libmesh_assert_less_equal (_coarsen_fraction, 1);
190 
191  // How much error per cell will we tolerate?
192  const Real local_refinement_tolerance =
193  _absolute_global_tolerance / std::sqrt(static_cast<Real>(_mesh.n_active_elem()));
194  const Real local_coarsening_tolerance =
195  local_refinement_tolerance * _coarsen_threshold;
196 
197  // Prepare another error vector if we need to sum parent errors
198  ErrorVector error_per_parent;
200  {
201  Real parent_error_min, parent_error_max;
202 
203  create_parent_error_vector(error_per_cell_in,
204  error_per_parent,
205  parent_error_min,
206  parent_error_max);
207  }
208 
211 
212  for (; elem_it != elem_end; ++elem_it)
213  {
214  Elem* elem = *elem_it;
215  Elem* parent = elem->parent();
216  const dof_id_type elem_number = elem->id();
217  const ErrorVectorReal elem_error = error_per_cell_in[elem_number];
218 
219  if (elem_error > local_refinement_tolerance &&
220  elem->level() < _max_h_level)
222 
223  if (!_coarsen_by_parents && elem_error <
224  local_coarsening_tolerance)
226 
227  if (_coarsen_by_parents && parent)
228  {
229  ErrorVectorReal parent_error = error_per_parent[parent->id()];
230  if (parent_error >= 0.)
231  {
232  const Real parent_coarsening_tolerance =
233  std::sqrt(parent->n_children() *
234  local_coarsening_tolerance *
235  local_coarsening_tolerance);
236  if (parent_error < parent_coarsening_tolerance)
238  }
239  }
240  }
241 }
242 
243 
244 
246 {
247  parallel_object_only();
248 
249  // Check for valid fractions..
250  // The fraction values must be in [0,1]
251  libmesh_assert_greater_equal (_refine_fraction, 0);
252  libmesh_assert_less_equal (_refine_fraction, 1);
253  libmesh_assert_greater_equal (_coarsen_fraction, 0);
254  libmesh_assert_less_equal (_coarsen_fraction, 1);
255 
256  // This function is currently only coded to work when coarsening by
257  // parents - it's too hard to guess how many coarsenings will be
258  // performed otherwise.
260 
261  // The number of active elements in the mesh - hopefully less than
262  // 2 billion on 32 bit machines
263  const dof_id_type n_active_elem = _mesh.n_active_elem();
264 
265  // The maximum number of active elements to flag for coarsening
266  const dof_id_type max_elem_coarsen =
267  static_cast<dof_id_type>(_coarsen_fraction * n_active_elem) + 1;
268 
269  // The maximum number of elements to flag for refinement
270  const dof_id_type max_elem_refine =
271  static_cast<dof_id_type>(_refine_fraction * n_active_elem) + 1;
272 
273  // Clean up the refinement flags. These could be left
274  // over from previous refinement steps.
275  this->clean_refinement_flags();
276 
277  // The target number of elements to add or remove
278  const std::ptrdiff_t n_elem_new = _nelem_target - n_active_elem;
279 
280  // Create an vector with active element errors and ids,
281  // sorted by highest errors first
282  const dof_id_type max_elem_id = _mesh.max_elem_id();
283  std::vector<std::pair<ErrorVectorReal, dof_id_type> > sorted_error;
284 
285  sorted_error.reserve (n_active_elem);
286 
287  // On a ParallelMesh, we need to communicate to know which remote ids
288  // correspond to active elements.
289  {
290  std::vector<bool> is_active(max_elem_id, false);
291 
294  for (; elem_it != elem_end; ++elem_it)
295  {
296  const dof_id_type eid = (*elem_it)->id();
297  is_active[eid] = true;
298  libmesh_assert_less (eid, error_per_cell.size());
299  sorted_error.push_back
300  (std::make_pair(error_per_cell[eid], eid));
301  }
302 
303  this->comm().max(is_active);
304 
305  this->comm().allgather(sorted_error);
306  }
307 
308  // Default sort works since pairs are sorted lexicographically
309  std::sort (sorted_error.begin(), sorted_error.end());
310  std::reverse (sorted_error.begin(), sorted_error.end());
311 
312  // Create a sorted error vector with coarsenable parent elements
313  // only, sorted by lowest errors first
314  ErrorVector error_per_parent;
315  std::vector<std::pair<ErrorVectorReal, dof_id_type> > sorted_parent_error;
316  Real parent_error_min, parent_error_max;
317 
318  create_parent_error_vector(error_per_cell,
319  error_per_parent,
320  parent_error_min,
321  parent_error_max);
322 
323  // create_parent_error_vector sets values for non-parents and
324  // non-coarsenable parents to -1. Get rid of them.
325  for (dof_id_type i=0; i != error_per_parent.size(); ++i)
326  if (error_per_parent[i] != -1)
327  sorted_parent_error.push_back(std::make_pair(error_per_parent[i], i));
328 
329  std::sort (sorted_parent_error.begin(), sorted_parent_error.end());
330 
331  // Keep track of how many elements we plan to coarsen & refine
332  dof_id_type coarsen_count = 0;
333  dof_id_type refine_count = 0;
334 
335  const unsigned int dim = _mesh.mesh_dimension();
336  unsigned int twotodim = 1;
337  for (unsigned int i=0; i!=dim; ++i)
338  twotodim *= 2;
339 
340  // First, let's try to get our element count to target_nelem
341  if (n_elem_new >= 0)
342  {
343  // Every element refinement creates at least
344  // 2^dim-1 new elements
345  refine_count =
346  std::min(libmesh_cast_int<dof_id_type>(n_elem_new / (twotodim-1)),
347  max_elem_refine);
348  }
349  else
350  {
351  // Every successful element coarsening is likely to destroy
352  // 2^dim-1 net elements.
353  coarsen_count =
354  std::min(libmesh_cast_int<dof_id_type>(-n_elem_new / (twotodim-1)),
355  max_elem_coarsen);
356  }
357 
358  // Next, let's see if we can trade any refinement for coarsening
359  while (coarsen_count < max_elem_coarsen &&
360  refine_count < max_elem_refine &&
361  coarsen_count < sorted_parent_error.size() &&
362  refine_count < sorted_error.size() &&
363  sorted_error[refine_count].first >
364  sorted_parent_error[coarsen_count].first * _coarsen_threshold)
365  {
366  coarsen_count++;
367  refine_count++;
368  }
369 
370  // On a ParallelMesh, we need to communicate to know which remote ids
371  // correspond to refinable elements
372  dof_id_type successful_refine_count = 0;
373  {
374  std::vector<bool> is_refinable(max_elem_id, false);
375 
376  for (dof_id_type i=0; i != sorted_error.size(); ++i)
377  {
378  dof_id_type eid = sorted_error[i].second;
379  Elem *elem = _mesh.query_elem(eid);
380  if (elem && elem->level() < _max_h_level)
381  is_refinable[eid] = true;
382  }
383  this->comm().max(is_refinable);
384 
385  if (refine_count > max_elem_refine)
386  refine_count = max_elem_refine;
387  for (dof_id_type i=0; i != sorted_error.size(); ++i)
388  {
389  if (successful_refine_count >= refine_count)
390  break;
391 
392  dof_id_type eid = sorted_error[i].second;
393  Elem *elem = _mesh.query_elem(eid);
394  if (is_refinable[eid])
395  {
396  if (elem)
398  successful_refine_count++;
399  }
400  }
401  }
402 
403  // If we couldn't refine enough elements, don't coarsen too many
404  // either
405  if (coarsen_count < (refine_count - successful_refine_count))
406  coarsen_count = 0;
407  else
408  coarsen_count -= (refine_count - successful_refine_count);
409 
410  if (coarsen_count > max_elem_coarsen)
411  coarsen_count = max_elem_coarsen;
412 
413  dof_id_type successful_coarsen_count = 0;
414  if (coarsen_count)
415  {
416  for (dof_id_type i=0; i != sorted_parent_error.size(); ++i)
417  {
418  if (successful_coarsen_count >= coarsen_count * twotodim)
419  break;
420 
421  dof_id_type parent_id = sorted_parent_error[i].second;
422  Elem *parent = _mesh.query_elem(parent_id);
423 
424  // On a ParallelMesh we skip remote elements
425  if (!parent)
426  continue;
427 
428  libmesh_assert(parent->has_children());
429  for (unsigned int c=0; c != parent->n_children(); ++c)
430  {
431  Elem *elem = parent->child(c);
432  if (elem && elem != remote_elem)
433  {
434  libmesh_assert(elem->active());
436  successful_coarsen_count++;
437  }
438  }
439  }
440  }
441 
442  // Return true if we've done all the AMR/C we can
443  if (!successful_coarsen_count &&
444  !successful_refine_count)
445  return true;
446  // And false if there may still be more to do.
447  return false;
448 }
449 
450 
451 
453  const Real refine_frac,
454  const Real coarsen_frac,
455  const unsigned int max_l)
456 {
457  parallel_object_only();
458 
459  // The function arguments are currently just there for
460  // backwards_compatibility
462  {
463  // If the user used non-default parameters, lets warn
464  // that they're deprecated
465  if (refine_frac != 0.3 ||
466  coarsen_frac != 0.0 ||
467  max_l != libMesh::invalid_uint)
468  libmesh_deprecated();
469 
470  _refine_fraction = refine_frac;
471  _coarsen_fraction = coarsen_frac;
472  _max_h_level = max_l;
473  }
474 
475  // Check for valid fractions..
476  // The fraction values must be in [0,1]
477  libmesh_assert_greater_equal (_refine_fraction, 0);
478  libmesh_assert_less_equal (_refine_fraction, 1);
479  libmesh_assert_greater_equal (_coarsen_fraction, 0);
480  libmesh_assert_less_equal (_coarsen_fraction, 1);
481 
482  // The number of active elements in the mesh
483  const dof_id_type n_active_elem = _mesh.n_elem();
484 
485  // The number of elements to flag for coarsening
486  const dof_id_type n_elem_coarsen =
487  static_cast<dof_id_type>(_coarsen_fraction * n_active_elem);
488 
489  // The number of elements to flag for refinement
490  const dof_id_type n_elem_refine =
491  static_cast<dof_id_type>(_refine_fraction * n_active_elem);
492 
493 
494 
495  // Clean up the refinement flags. These could be left
496  // over from previous refinement steps.
497  this->clean_refinement_flags();
498 
499 
500  // This vector stores the error and element number for all the
501  // active elements. It will be sorted and the top & bottom
502  // elements will then be flagged for coarsening & refinement
503  std::vector<ErrorVectorReal> sorted_error;
504 
505  sorted_error.reserve (n_active_elem);
506 
507  // Loop over the active elements and create the entry
508  // in the sorted_error vector
511 
512  for (; elem_it != elem_end; ++elem_it)
513  sorted_error.push_back (error_per_cell[(*elem_it)->id()]);
514 
515  this->comm().allgather(sorted_error);
516 
517  // Now sort the sorted_error vector
518  std::sort (sorted_error.begin(), sorted_error.end());
519 
520  // If we're coarsening by parents:
521  // Create a sorted error vector with coarsenable parent elements
522  // only, sorted by lowest errors first
523  ErrorVector error_per_parent, sorted_parent_error;
525  {
526  Real parent_error_min, parent_error_max;
527 
528  create_parent_error_vector(error_per_cell,
529  error_per_parent,
530  parent_error_min,
531  parent_error_max);
532 
533  sorted_parent_error = error_per_parent;
534  std::sort (sorted_parent_error.begin(), sorted_parent_error.end());
535 
536  // All the other error values will be 0., so get rid of them.
537  sorted_parent_error.erase (std::remove(sorted_parent_error.begin(),
538  sorted_parent_error.end(), 0.),
539  sorted_parent_error.end());
540  }
541 
542 
543  ErrorVectorReal top_error= 0., bottom_error = 0.;
544 
545  // Get the maximum error value corresponding to the
546  // bottom n_elem_coarsen elements
547  if (_coarsen_by_parents && n_elem_coarsen)
548  {
549  const unsigned int dim = _mesh.mesh_dimension();
550  unsigned int twotodim = 1;
551  for (unsigned int i=0; i!=dim; ++i)
552  twotodim *= 2;
553 
554  dof_id_type n_parent_coarsen = n_elem_coarsen / (twotodim - 1);
555 
556  if (n_parent_coarsen)
557  bottom_error = sorted_parent_error[n_parent_coarsen - 1];
558  }
559  else if (n_elem_coarsen)
560  {
561  bottom_error = sorted_error[n_elem_coarsen - 1];
562  }
563 
564  if (n_elem_refine)
565  top_error = sorted_error[sorted_error.size() - n_elem_refine];
566 
567  // Finally, let's do the element flagging
568  elem_it = _mesh.active_elements_begin();
569  for (; elem_it != elem_end; ++elem_it)
570  {
571  Elem* elem = *elem_it;
572  Elem* parent = elem->parent();
573 
574  if (_coarsen_by_parents && parent && n_elem_coarsen &&
575  error_per_parent[parent->id()] <= bottom_error)
577 
578  if (!_coarsen_by_parents && n_elem_coarsen &&
579  error_per_cell[elem->id()] <= bottom_error)
581 
582  if (n_elem_refine &&
583  elem->level() < _max_h_level &&
584  error_per_cell[elem->id()] >= top_error)
586  }
587 }
588 
589 
590 
592  const Real refine_frac,
593  const Real coarsen_frac,
594  const unsigned int max_l)
595 {
596  // The function arguments are currently just there for
597  // backwards_compatibility
599  {
600  // If the user used non-default parameters, lets warn
601  // that they're deprecated
602  if (refine_frac != 0.3 ||
603  coarsen_frac != 0.0 ||
604  max_l != libMesh::invalid_uint)
605  libmesh_deprecated();
606 
607  _refine_fraction = refine_frac;
608  _coarsen_fraction = coarsen_frac;
609  _max_h_level = max_l;
610  }
611 
612  // Get the mean value from the error vector
613  const Real mean = error_per_cell.mean();
614 
615  // Get the standard deviation. This equals the
616  // square-root of the variance
617  const Real stddev = std::sqrt (error_per_cell.variance());
618 
619  // Check for valid fractions
620  libmesh_assert_greater_equal (_refine_fraction, 0);
621  libmesh_assert_less_equal (_refine_fraction, 1);
622  libmesh_assert_greater_equal (_coarsen_fraction, 0);
623  libmesh_assert_less_equal (_coarsen_fraction, 1);
624 
625  // The refine and coarsen cutoff
626  const Real refine_cutoff = mean + _refine_fraction * stddev;
627  const Real coarsen_cutoff = std::max(mean - _coarsen_fraction * stddev, 0.);
628 
629  // Loop over the elements and flag them for coarsening or
630  // refinement based on the element error
633 
634  for (; elem_it != elem_end; ++elem_it)
635  {
636  Elem* elem = *elem_it;
637  const dof_id_type id = elem->id();
638 
639  libmesh_assert_less (id, error_per_cell.size());
640 
641  const ErrorVectorReal elem_error = error_per_cell[id];
642 
643  // Possibly flag the element for coarsening ...
644  if (elem_error <= coarsen_cutoff)
646 
647  // ... or refinement
648  if ((elem_error >= refine_cutoff) && (elem->level() < _max_h_level))
650  }
651 }
652 
653 
654 
656 {
657  element_flagging.flag_elements();
658 }
659 
660 
661 
663 {
665  const MeshBase::element_iterator elem_end = _mesh.elements_end();
666 
667  for ( ; elem_it != elem_end; ++elem_it)
668  {
669  if ((*elem_it)->active())
670  {
671  (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
672  (*elem_it)->set_refinement_flag(Elem::DO_NOTHING);
673  }
674  else
675  {
676  (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
677  (*elem_it)->set_refinement_flag(Elem::INACTIVE);
678  }
679  }
680 }
681 
682 
683 
685 {
687  const MeshBase::element_iterator elem_end = _mesh.elements_end();
688 
689  for ( ; elem_it != elem_end; ++elem_it)
690  (*elem_it)->set_p_refinement_flag((*elem_it)->refinement_flag());
691 }
692 
693 
694 
696 {
697  // Possibly clean up the refinement flags from
698  // a previous step
699 // elem_iterator elem_it (_mesh.elements_begin());
700 // const elem_iterator elem_end(_mesh.elements_end());
701 
703  const MeshBase::element_iterator elem_end = _mesh.elements_end();
704 
705  for ( ; elem_it != elem_end; ++elem_it)
706  {
707  if ((*elem_it)->active())
708  {
709  (*elem_it)->set_refinement_flag(Elem::DO_NOTHING);
710  (*elem_it)->set_p_refinement_flag(Elem::DO_NOTHING);
711  }
712  else
713  {
714  (*elem_it)->set_refinement_flag(Elem::INACTIVE);
715  (*elem_it)->set_p_refinement_flag(Elem::INACTIVE);
716  }
717  }
718 }
719 
720 } // namespace libMesh
721 
722 #endif

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

Hosted By:
SourceForge.net Logo