Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // C++ includes
27 : #include <algorithm> // for std::sort
28 :
29 : // Local includes
30 : #include "libmesh/elem.h"
31 : #include "libmesh/error_vector.h"
32 : #include "libmesh/mesh_refinement.h"
33 : #include "libmesh/mesh_base.h"
34 : #include "libmesh/parallel.h"
35 : #include "libmesh/remote_elem.h"
36 :
37 : namespace libMesh
38 : {
39 :
40 :
41 :
42 : //-----------------------------------------------------------------
43 : // Mesh refinement methods
44 10908 : void MeshRefinement::flag_elements_by_error_fraction (const ErrorVector & error_per_cell,
45 : const Real refine_frac,
46 : const Real coarsen_frac,
47 : const unsigned int max_l)
48 : {
49 404 : parallel_object_only();
50 :
51 : // Verify that our error vector is consistent, using std::vector to
52 : // avoid confusing this->comm().verify
53 404 : libmesh_assert(this->comm().verify(dynamic_cast<const std::vector<ErrorVectorReal> &>(error_per_cell)));
54 :
55 : // The function arguments are currently just there for
56 : // backwards_compatibility
57 10908 : if (!_use_member_parameters)
58 : {
59 : // If the user used non-default parameters, lets warn
60 : // that they're deprecated
61 : if (refine_frac != 0.3 ||
62 : coarsen_frac != 0.0 ||
63 : max_l != libMesh::invalid_uint)
64 : libmesh_deprecated();
65 :
66 0 : _refine_fraction = refine_frac;
67 0 : _coarsen_fraction = coarsen_frac;
68 0 : _max_h_level = max_l;
69 : }
70 :
71 : // Check for valid fractions..
72 : // The fraction values must be in [0,1]
73 404 : libmesh_assert_greater_equal (_refine_fraction, 0);
74 404 : libmesh_assert_less_equal (_refine_fraction, 1);
75 404 : libmesh_assert_greater_equal (_coarsen_fraction, 0);
76 404 : libmesh_assert_less_equal (_coarsen_fraction, 1);
77 :
78 : // Clean up the refinement flags. These could be left
79 : // over from previous refinement steps.
80 10908 : this->clean_refinement_flags();
81 :
82 : // We're getting the minimum and maximum error values
83 : // for the ACTIVE elements
84 10908 : Real error_min = 1.e30;
85 10908 : Real error_max = 0.;
86 :
87 : // And, if necessary, for their parents
88 10908 : Real parent_error_min = 1.e30;
89 10908 : Real parent_error_max = 0.;
90 :
91 : // Prepare another error vector if we need to sum parent errors
92 808 : ErrorVector error_per_parent;
93 10908 : if (_coarsen_by_parents)
94 : {
95 0 : create_parent_error_vector(error_per_cell,
96 : error_per_parent,
97 : parent_error_min,
98 : parent_error_max);
99 : }
100 :
101 : // We need to loop over all active elements to find the minimum
102 979722 : for (auto & elem : _mesh.active_local_element_ptr_range())
103 : {
104 958310 : const dof_id_type id = elem->id();
105 108750 : libmesh_assert_less (id, error_per_cell.size());
106 :
107 1067060 : error_max = std::max (error_max, error_per_cell[id]);
108 1020739 : error_min = std::min (error_min, error_per_cell[id]);
109 10100 : }
110 10908 : this->comm().max(error_max);
111 10908 : this->comm().min(error_min);
112 :
113 : // Compute the cutoff values for coarsening and refinement
114 10908 : const Real error_delta = (error_max - error_min);
115 10908 : const Real parent_error_delta = parent_error_max - parent_error_min;
116 :
117 10908 : const Real refine_cutoff = (1.- _refine_fraction)*error_max;
118 10908 : const Real coarsen_cutoff = _coarsen_fraction*error_delta + error_min;
119 10908 : const Real parent_cutoff = _coarsen_fraction*parent_error_delta + error_min;
120 :
121 : // // Print information about the error
122 : // libMesh::out << " Error Information:" << std::endl
123 : // << " ------------------" << std::endl
124 : // << " min: " << error_min << std::endl
125 : // << " max: " << error_max << std::endl
126 : // << " delta: " << error_delta << std::endl
127 : // << " refine_cutoff: " << refine_cutoff << std::endl
128 : // << " coarsen_cutoff: " << coarsen_cutoff << std::endl;
129 :
130 :
131 :
132 : // Loop over the elements and flag them for coarsening or
133 : // refinement based on the element error
134 9568674 : for (auto & elem : _mesh.active_element_ptr_range())
135 : {
136 4991131 : const dof_id_type id = elem->id();
137 :
138 217500 : libmesh_assert_less (id, error_per_cell.size());
139 :
140 4991131 : const ErrorVectorReal elem_error = error_per_cell[id];
141 :
142 4991131 : if (_coarsen_by_parents)
143 : {
144 0 : Elem * parent = elem->parent();
145 0 : if (parent)
146 : {
147 0 : const dof_id_type parentid = parent->id();
148 0 : if (error_per_parent[parentid] >= 0. &&
149 0 : error_per_parent[parentid] <= parent_cutoff)
150 0 : elem->set_refinement_flag(Elem::COARSEN);
151 : }
152 : }
153 : // Flag the element for coarsening if its error
154 : // is <= coarsen_fraction*delta + error_min
155 4991131 : else if (elem_error <= coarsen_cutoff)
156 : {
157 72200 : elem->set_refinement_flag(Elem::COARSEN);
158 : }
159 :
160 : // Flag the element for refinement if its error
161 : // is >= refinement_cutoff.
162 4991131 : if (elem_error >= refine_cutoff)
163 1000769 : if (elem->level() < _max_h_level)
164 135626 : elem->set_refinement_flag(Elem::REFINE);
165 10100 : }
166 10908 : }
167 :
168 :
169 :
170 0 : void MeshRefinement::flag_elements_by_error_tolerance (const ErrorVector & error_per_cell_in)
171 : {
172 0 : parallel_object_only();
173 :
174 : // Verify that our error vector is consistent, using std::vector to
175 : // avoid confusing this->comm().verify
176 0 : libmesh_assert(this->comm().verify(dynamic_cast<const std::vector<ErrorVectorReal> &>(error_per_cell_in)));
177 :
178 0 : libmesh_assert_greater (_coarsen_threshold, 0);
179 :
180 : // Check for valid fractions..
181 : // The fraction values must be in [0,1]
182 0 : libmesh_assert_greater_equal (_refine_fraction, 0);
183 0 : libmesh_assert_less_equal (_refine_fraction, 1);
184 0 : libmesh_assert_greater_equal (_coarsen_fraction, 0);
185 0 : libmesh_assert_less_equal (_coarsen_fraction, 1);
186 :
187 : // How much error per cell will we tolerate?
188 : const Real local_refinement_tolerance =
189 0 : _absolute_global_tolerance / std::sqrt(static_cast<Real>(_mesh.n_active_elem()));
190 0 : const Real local_coarsening_tolerance =
191 0 : local_refinement_tolerance * _coarsen_threshold;
192 :
193 : // Prepare another error vector if we need to sum parent errors
194 0 : ErrorVector error_per_parent;
195 0 : if (_coarsen_by_parents)
196 : {
197 : Real parent_error_min, parent_error_max;
198 :
199 0 : create_parent_error_vector(error_per_cell_in,
200 : error_per_parent,
201 : parent_error_min,
202 : parent_error_max);
203 : }
204 :
205 0 : for (auto & elem : _mesh.active_element_ptr_range())
206 : {
207 0 : Elem * parent = elem->parent();
208 0 : const dof_id_type elem_number = elem->id();
209 0 : const ErrorVectorReal elem_error = error_per_cell_in[elem_number];
210 :
211 0 : if (elem_error > local_refinement_tolerance &&
212 0 : elem->level() < _max_h_level)
213 0 : elem->set_refinement_flag(Elem::REFINE);
214 :
215 0 : if (!_coarsen_by_parents && elem_error <
216 : local_coarsening_tolerance)
217 0 : elem->set_refinement_flag(Elem::COARSEN);
218 :
219 0 : if (_coarsen_by_parents && parent)
220 : {
221 0 : ErrorVectorReal parent_error = error_per_parent[parent->id()];
222 0 : if (parent_error >= 0.)
223 : {
224 : const Real parent_coarsening_tolerance =
225 0 : std::sqrt(parent->n_children() *
226 : local_coarsening_tolerance *
227 : local_coarsening_tolerance);
228 0 : if (parent_error < parent_coarsening_tolerance)
229 0 : elem->set_refinement_flag(Elem::COARSEN);
230 : }
231 : }
232 0 : }
233 0 : }
234 :
235 :
236 :
237 1590 : bool MeshRefinement::flag_elements_by_nelem_target (const ErrorVector & error_per_cell)
238 : {
239 82 : parallel_object_only();
240 :
241 : // Verify that our error vector is consistent, using std::vector to
242 : // avoid confusing this->comm().verify
243 82 : libmesh_assert(this->comm().verify(dynamic_cast<const std::vector<ErrorVectorReal> &>(error_per_cell)));
244 :
245 : // Check for valid fractions..
246 : // The fraction values must be in [0,1]
247 82 : libmesh_assert_greater_equal (_refine_fraction, 0);
248 82 : libmesh_assert_less_equal (_refine_fraction, 1);
249 82 : libmesh_assert_greater_equal (_coarsen_fraction, 0);
250 82 : libmesh_assert_less_equal (_coarsen_fraction, 1);
251 :
252 : // This function is currently only coded to work when coarsening by
253 : // parents - it's too hard to guess how many coarsenings will be
254 : // performed otherwise.
255 82 : libmesh_assert (_coarsen_by_parents);
256 :
257 : // The number of active elements in the mesh - hopefully less than
258 : // 2 billion on 32 bit machines
259 1590 : const dof_id_type n_active_elem = _mesh.n_active_elem();
260 :
261 : // The maximum number of active elements to flag for coarsening
262 1590 : const dof_id_type max_elem_coarsen =
263 1590 : static_cast<dof_id_type>(_coarsen_fraction * n_active_elem) + 1;
264 :
265 : // The maximum number of elements to flag for refinement
266 1590 : const dof_id_type max_elem_refine =
267 1590 : static_cast<dof_id_type>(_refine_fraction * n_active_elem) + 1;
268 :
269 : // Clean up the refinement flags. These could be left
270 : // over from previous refinement steps.
271 1590 : this->clean_refinement_flags();
272 :
273 : // The target number of elements to add or remove
274 1590 : const std::ptrdiff_t n_elem_new =
275 1590 : std::ptrdiff_t(_nelem_target) - std::ptrdiff_t(n_active_elem);
276 :
277 : // Create an vector with active element errors and ids,
278 : // sorted by highest errors first
279 1590 : const dof_id_type max_elem_id = _mesh.max_elem_id();
280 164 : std::vector<std::pair<ErrorVectorReal, dof_id_type>> sorted_error;
281 :
282 1590 : sorted_error.reserve (n_active_elem);
283 :
284 : // On a DistributedMesh, we need to communicate to know which remote ids
285 : // correspond to active elements.
286 : {
287 1672 : std::vector<bool> is_active(max_elem_id, false);
288 :
289 79026 : for (auto & elem : _mesh.active_local_element_ptr_range())
290 : {
291 42484 : const dof_id_type eid = elem->id();
292 12646 : is_active[eid] = true;
293 4520 : libmesh_assert_less (eid, error_per_cell.size());
294 47004 : sorted_error.emplace_back(error_per_cell[eid], eid);
295 1426 : }
296 :
297 1590 : this->comm().max(is_active);
298 :
299 1590 : this->comm().allgather(sorted_error);
300 : }
301 :
302 : // Default sort works since pairs are sorted lexicographically
303 1590 : std::sort (sorted_error.begin(), sorted_error.end());
304 1508 : std::reverse (sorted_error.begin(), sorted_error.end());
305 :
306 : // Create a sorted error vector with coarsenable parent elements
307 : // only, sorted by lowest errors first
308 164 : ErrorVector error_per_parent;
309 164 : std::vector<std::pair<ErrorVectorReal, dof_id_type>> sorted_parent_error;
310 : Real parent_error_min, parent_error_max;
311 :
312 1590 : create_parent_error_vector(error_per_cell,
313 : error_per_parent,
314 : parent_error_min,
315 : parent_error_max);
316 :
317 : // create_parent_error_vector sets values for non-parents and
318 : // non-coarsenable parents to -1. Get rid of them.
319 302928 : for (auto i : index_range(error_per_parent))
320 312918 : if (error_per_parent[i] != -1)
321 43980 : sorted_parent_error.emplace_back(error_per_parent[i], i);
322 :
323 1590 : std::sort (sorted_parent_error.begin(), sorted_parent_error.end());
324 :
325 : // Keep track of how many elements we plan to coarsen & refine
326 82 : dof_id_type coarsen_count = 0;
327 82 : dof_id_type refine_count = 0;
328 :
329 1590 : const unsigned int dim = _mesh.mesh_dimension();
330 82 : unsigned int twotodim = 1;
331 4770 : for (unsigned int i=0; i!=dim; ++i)
332 3180 : twotodim *= 2;
333 :
334 : // First, let's try to get our element count to target_nelem
335 1590 : if (n_elem_new >= 0)
336 : {
337 : // Every element refinement creates at least
338 : // 2^dim-1 new elements
339 1590 : refine_count =
340 3262 : std::min(cast_int<dof_id_type>(n_elem_new / (twotodim-1)),
341 82 : max_elem_refine);
342 : }
343 : else
344 : {
345 : // Every successful element coarsening is likely to destroy
346 : // 2^dim-1 net elements.
347 0 : coarsen_count =
348 0 : std::min(cast_int<dof_id_type>(-n_elem_new / (twotodim-1)),
349 0 : max_elem_coarsen);
350 : }
351 :
352 : // Next, let's see if we can trade any refinement for coarsening
353 1672 : while (coarsen_count < max_elem_coarsen &&
354 82 : refine_count < max_elem_refine &&
355 0 : coarsen_count < sorted_parent_error.size() &&
356 1672 : refine_count < sorted_error.size() &&
357 0 : sorted_error[refine_count].first >
358 0 : sorted_parent_error[coarsen_count].first * _coarsen_threshold)
359 : {
360 0 : coarsen_count++;
361 0 : refine_count++;
362 : }
363 :
364 : // On a DistributedMesh, we need to communicate to know which remote ids
365 : // correspond to refinable elements
366 82 : dof_id_type successful_refine_count = 0;
367 : {
368 1672 : std::vector<bool> is_refinable(max_elem_id, false);
369 :
370 239189 : for (const auto & pr : sorted_error)
371 : {
372 237599 : dof_id_type eid = pr.second;
373 237599 : Elem * elem = _mesh.query_elem_ptr(eid);
374 237599 : if (elem && elem->level() < _max_h_level)
375 18080 : is_refinable[eid] = true;
376 : }
377 1590 : this->comm().max(is_refinable);
378 :
379 164 : if (refine_count > max_elem_refine)
380 0 : refine_count = max_elem_refine;
381 29362 : for (const auto & pr : sorted_error)
382 : {
383 29362 : if (successful_refine_count >= refine_count)
384 82 : break;
385 :
386 27772 : dof_id_type eid = pr.second;
387 27772 : Elem * elem = _mesh.query_elem_ptr(eid);
388 28690 : if (is_refinable[eid])
389 : {
390 27772 : if (elem)
391 918 : elem->set_refinement_flag(Elem::REFINE);
392 27772 : successful_refine_count++;
393 : }
394 : }
395 : }
396 :
397 : // If we couldn't refine enough elements, don't coarsen too many
398 : // either
399 1590 : if (coarsen_count < (refine_count - successful_refine_count))
400 0 : coarsen_count = 0;
401 : else
402 1590 : coarsen_count -= (refine_count - successful_refine_count);
403 :
404 164 : if (coarsen_count > max_elem_coarsen)
405 0 : coarsen_count = max_elem_coarsen;
406 :
407 82 : dof_id_type successful_coarsen_count = 0;
408 1590 : if (coarsen_count)
409 : {
410 0 : for (const auto & pr : sorted_parent_error)
411 : {
412 0 : if (successful_coarsen_count >= coarsen_count * twotodim)
413 0 : break;
414 :
415 0 : dof_id_type parent_id = pr.second;
416 0 : Elem * parent = _mesh.query_elem_ptr(parent_id);
417 :
418 : // On a DistributedMesh we skip remote elements
419 0 : if (!parent)
420 0 : continue;
421 :
422 0 : libmesh_assert(parent->has_children());
423 0 : for (auto & elem : parent->child_ref_range())
424 : {
425 0 : if (&elem != remote_elem)
426 : {
427 0 : libmesh_assert(elem.active());
428 0 : elem.set_refinement_flag(Elem::COARSEN);
429 0 : successful_coarsen_count++;
430 : }
431 : }
432 : }
433 : }
434 :
435 : // Return true if we've done all the AMR/C we can
436 1590 : if (!successful_coarsen_count &&
437 : !successful_refine_count)
438 0 : return true;
439 : // And false if there may still be more to do.
440 82 : return false;
441 : }
442 :
443 :
444 :
445 2982 : void MeshRefinement::flag_elements_by_elem_fraction (const ErrorVector & error_per_cell,
446 : const Real refine_frac,
447 : const Real coarsen_frac,
448 : const unsigned int max_l)
449 : {
450 84 : parallel_object_only();
451 :
452 : // Verify that our error vector is consistent, using std::vector to
453 : // avoid confusing this->comm().verify
454 84 : libmesh_assert(this->comm().verify(dynamic_cast<const std::vector<ErrorVectorReal> &>(error_per_cell)));
455 :
456 : // The function arguments are currently just there for
457 : // backwards_compatibility
458 2982 : if (!_use_member_parameters)
459 : {
460 : // If the user used non-default parameters, lets warn
461 : // that they're deprecated
462 : if (refine_frac != 0.3 ||
463 : coarsen_frac != 0.0 ||
464 : max_l != libMesh::invalid_uint)
465 : libmesh_deprecated();
466 :
467 0 : _refine_fraction = refine_frac;
468 0 : _coarsen_fraction = coarsen_frac;
469 0 : _max_h_level = max_l;
470 : }
471 :
472 : // Check for valid fractions..
473 : // The fraction values must be in [0,1]
474 84 : libmesh_assert_greater_equal (_refine_fraction, 0);
475 84 : libmesh_assert_less_equal (_refine_fraction, 1);
476 84 : libmesh_assert_greater_equal (_coarsen_fraction, 0);
477 84 : libmesh_assert_less_equal (_coarsen_fraction, 1);
478 :
479 : // The number of active elements in the mesh
480 2982 : const dof_id_type n_active_elem = _mesh.n_active_elem();
481 :
482 : // The number of elements to flag for coarsening
483 2982 : const dof_id_type n_elem_coarsen =
484 2982 : static_cast<dof_id_type>(_coarsen_fraction * n_active_elem);
485 :
486 : // The number of elements to flag for refinement
487 2982 : const dof_id_type n_elem_refine =
488 2982 : static_cast<dof_id_type>(_refine_fraction * n_active_elem);
489 :
490 :
491 :
492 : // Clean up the refinement flags. These could be left
493 : // over from previous refinement steps.
494 2982 : this->clean_refinement_flags();
495 :
496 :
497 : // This vector stores the error and element number for all the
498 : // active elements. It will be sorted and the top & bottom
499 : // elements will then be flagged for coarsening & refinement
500 168 : std::vector<ErrorVectorReal> sorted_error;
501 :
502 2982 : sorted_error.reserve (n_active_elem);
503 :
504 : // Loop over the active elements and create the entry
505 : // in the sorted_error vector
506 47980 : for (auto & elem : _mesh.active_local_element_ptr_range())
507 27676 : sorted_error.push_back (error_per_cell[elem->id()]);
508 :
509 2982 : this->comm().allgather(sorted_error);
510 :
511 : // Now sort the sorted_error vector
512 2982 : std::sort (sorted_error.begin(), sorted_error.end());
513 :
514 : // If we're coarsening by parents:
515 : // Create a sorted error vector with coarsenable parent elements
516 : // only, sorted by lowest errors first
517 168 : ErrorVector error_per_parent, sorted_parent_error;
518 2982 : if (_coarsen_by_parents)
519 : {
520 : Real parent_error_min, parent_error_max;
521 :
522 0 : create_parent_error_vector(error_per_cell,
523 : error_per_parent,
524 : parent_error_min,
525 : parent_error_max);
526 :
527 0 : sorted_parent_error = error_per_parent;
528 0 : std::sort (sorted_parent_error.begin(), sorted_parent_error.end());
529 :
530 : // All the other error values will be 0., so get rid of them.
531 0 : sorted_parent_error.erase (std::remove(sorted_parent_error.begin(),
532 0 : sorted_parent_error.end(), 0.),
533 0 : sorted_parent_error.end());
534 : }
535 :
536 :
537 84 : ErrorVectorReal top_error= 0., bottom_error = 0.;
538 :
539 : // Get the maximum error value corresponding to the
540 : // bottom n_elem_coarsen elements
541 2982 : if (_coarsen_by_parents && n_elem_coarsen)
542 : {
543 0 : const unsigned int dim = _mesh.mesh_dimension();
544 0 : unsigned int twotodim = 1;
545 0 : for (unsigned int i=0; i!=dim; ++i)
546 0 : twotodim *= 2;
547 :
548 0 : dof_id_type n_parent_coarsen = n_elem_coarsen / (twotodim - 1);
549 :
550 0 : if (n_parent_coarsen)
551 0 : bottom_error = sorted_parent_error[n_parent_coarsen - 1];
552 : }
553 2982 : else if (n_elem_coarsen)
554 : {
555 0 : bottom_error = sorted_error[n_elem_coarsen - 1];
556 : }
557 :
558 2982 : if (n_elem_refine)
559 2100 : top_error = sorted_error[sorted_error.size() - n_elem_refine];
560 :
561 : // Finally, let's do the element flagging
562 93738 : for (auto & elem : _mesh.active_element_ptr_range())
563 : {
564 47741 : Elem * parent = elem->parent();
565 :
566 47741 : if (_coarsen_by_parents && parent && n_elem_coarsen &&
567 0 : error_per_parent[parent->id()] <= bottom_error)
568 0 : elem->set_refinement_flag(Elem::COARSEN);
569 :
570 47741 : if (!_coarsen_by_parents && n_elem_coarsen &&
571 0 : error_per_cell[elem->id()] <= bottom_error)
572 0 : elem->set_refinement_flag(Elem::COARSEN);
573 :
574 50937 : if (n_elem_refine &&
575 51497 : elem->level() < _max_h_level &&
576 50881 : error_per_cell[elem->id()] >= top_error)
577 996 : elem->set_refinement_flag(Elem::REFINE);
578 2814 : }
579 2982 : }
580 :
581 :
582 :
583 0 : void MeshRefinement::flag_elements_by_mean_stddev (const ErrorVector & error_per_cell,
584 : const Real refine_frac,
585 : const Real coarsen_frac,
586 : const unsigned int max_l)
587 : {
588 : // Verify that our error vector is consistent, using std::vector to
589 : // avoid confusing this->comm().verify
590 0 : libmesh_assert(this->comm().verify(dynamic_cast<const std::vector<ErrorVectorReal> &>(error_per_cell)));
591 :
592 : // The function arguments are currently just there for
593 : // backwards_compatibility
594 0 : if (!_use_member_parameters)
595 : {
596 : // If the user used non-default parameters, lets warn
597 : // that they're deprecated
598 : if (refine_frac != 0.3 ||
599 : coarsen_frac != 0.0 ||
600 : max_l != libMesh::invalid_uint)
601 : libmesh_deprecated();
602 :
603 0 : _refine_fraction = refine_frac;
604 0 : _coarsen_fraction = coarsen_frac;
605 0 : _max_h_level = max_l;
606 : }
607 :
608 : // Get the mean value from the error vector
609 0 : const Real mean = error_per_cell.mean();
610 :
611 : // Get the standard deviation. This equals the
612 : // square-root of the variance
613 0 : const Real stddev = std::sqrt (error_per_cell.variance());
614 :
615 : // Check for valid fractions
616 0 : libmesh_assert_greater_equal (_refine_fraction, 0);
617 0 : libmesh_assert_less_equal (_refine_fraction, 1);
618 0 : libmesh_assert_greater_equal (_coarsen_fraction, 0);
619 0 : libmesh_assert_less_equal (_coarsen_fraction, 1);
620 :
621 : // The refine and coarsen cutoff
622 0 : const Real refine_cutoff = mean + _refine_fraction * stddev;
623 0 : const Real coarsen_cutoff = std::max(mean - _coarsen_fraction * stddev, 0.);
624 :
625 : // Loop over the elements and flag them for coarsening or
626 : // refinement based on the element error
627 0 : for (auto & elem : _mesh.active_element_ptr_range())
628 : {
629 0 : const dof_id_type id = elem->id();
630 :
631 0 : libmesh_assert_less (id, error_per_cell.size());
632 :
633 0 : const ErrorVectorReal elem_error = error_per_cell[id];
634 :
635 : // Possibly flag the element for coarsening ...
636 0 : if (elem_error <= coarsen_cutoff)
637 0 : elem->set_refinement_flag(Elem::COARSEN);
638 :
639 : // ... or refinement
640 0 : if ((elem_error >= refine_cutoff) && (elem->level() < _max_h_level))
641 0 : elem->set_refinement_flag(Elem::REFINE);
642 0 : }
643 0 : }
644 :
645 :
646 :
647 0 : void MeshRefinement::flag_elements_by (ElementFlagging & element_flagging)
648 : {
649 0 : element_flagging.flag_elements();
650 0 : }
651 :
652 :
653 :
654 1349 : void MeshRefinement::switch_h_to_p_refinement ()
655 : {
656 40330 : for (auto & elem : _mesh.element_ptr_range())
657 : {
658 20193 : if (elem->active())
659 : {
660 1196 : elem->set_p_refinement_flag(elem->refinement_flag());
661 2392 : elem->set_refinement_flag(Elem::DO_NOTHING);
662 : }
663 : else
664 : {
665 162 : elem->set_p_refinement_flag(elem->refinement_flag());
666 324 : elem->set_refinement_flag(Elem::INACTIVE);
667 : }
668 1273 : }
669 1349 : }
670 :
671 :
672 :
673 639 : void MeshRefinement::add_p_to_h_refinement ()
674 : {
675 5402 : for (auto & elem : _mesh.element_ptr_range())
676 2808 : elem->set_p_refinement_flag(elem->refinement_flag());
677 639 : }
678 :
679 :
680 :
681 276852 : void MeshRefinement::clean_refinement_flags ()
682 : {
683 : // Possibly clean up the refinement flags from
684 : // a previous step
685 36362504 : for (auto & elem : _mesh.element_ptr_range())
686 : {
687 19280618 : if (elem->active())
688 : {
689 1103786 : elem->set_refinement_flag(Elem::DO_NOTHING);
690 2207572 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
691 : }
692 : else
693 : {
694 268372 : elem->set_refinement_flag(Elem::INACTIVE);
695 536744 : elem->set_p_refinement_flag(Elem::INACTIVE);
696 : }
697 260612 : }
698 276852 : }
699 :
700 : } // namespace libMesh
701 :
702 : #endif
|
