Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
20 : #include "libmesh/libmesh_config.h"
21 :
22 : // only compile these functions if the user requests AMR support
23 : #ifdef LIBMESH_ENABLE_AMR
24 :
25 : #include "libmesh/boundary_info.h"
26 : #include "libmesh/elem_range.h"
27 : #include "libmesh/error_vector.h"
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/mesh_base.h"
30 : #include "libmesh/mesh_communication.h"
31 : #include "libmesh/mesh_refinement.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/parallel_ghost_sync.h"
34 : #include "libmesh/partitioner.h"
35 : #include "libmesh/remote_elem.h"
36 : #include "libmesh/sync_refinement_flags.h"
37 : #include "libmesh/int_range.h"
38 :
39 : #ifdef DEBUG
40 : // Some extra validation for DistributedMesh
41 : #include "libmesh/mesh_tools.h"
42 : #endif // DEBUG
43 :
44 : #ifdef LIBMESH_ENABLE_PERIODIC
45 : #include "libmesh/periodic_boundaries.h"
46 : #endif
47 :
48 : // C++ includes
49 : #include <atomic>
50 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
51 : #include <cmath> // for isnan(), when it's defined
52 : #include <limits>
53 :
54 :
55 : // Anonymous namespace for helper functions
56 : // namespace {
57 : //
58 : // using namespace libMesh;
59 : //
60 : // struct SyncCoarsenInactive
61 : // {
62 : // bool operator() (const Elem * elem) const {
63 : // // If we're not an ancestor, there's no chance our coarsening
64 : // // settings need to be changed.
65 : // if (!elem->ancestor())
66 : // return false;
67 : //
68 : // // If we don't have any remote children, we already know enough to
69 : // // determine the correct refinement flag ourselves.
70 : // //
71 : // // If we know we have a child that isn't being coarsened, that
72 : // // also forces a specific flag.
73 : // //
74 : // // Either way there's nothing we need to communicate.
75 : // bool found_remote_child = false;
76 : // for (auto & child : elem->child_ref_range())
77 : // {
78 : // if (child.refinement_flag() != Elem::COARSEN)
79 : // return false;
80 : // if (&child == remote_elem)
81 : // found_remote_child = true;
82 : // }
83 : // return found_remote_child;
84 : // }
85 : // };
86 : // }
87 :
88 :
89 :
90 : namespace libMesh
91 : {
92 :
93 : //-----------------------------------------------------------------
94 : // Mesh refinement methods
95 309181 : MeshRefinement::MeshRefinement (MeshBase & m) :
96 : ParallelObject(m),
97 291233 : _mesh(m),
98 291233 : _use_member_parameters(false),
99 291233 : _coarsen_by_parents(false),
100 291233 : _refine_fraction(0.3),
101 291233 : _coarsen_fraction(0.0),
102 291233 : _max_h_level(libMesh::invalid_uint),
103 291233 : _coarsen_threshold(10),
104 291233 : _nelem_target(0),
105 291233 : _absolute_global_tolerance(0.0),
106 291233 : _face_level_mismatch_limit(1),
107 291233 : _edge_level_mismatch_limit(0),
108 291233 : _node_level_mismatch_limit(0),
109 291233 : _overrefined_boundary_limit(0),
110 291233 : _underrefined_boundary_limit(0),
111 291233 : _allow_unrefined_patches(false),
112 291233 : _enforce_mismatch_limit_prior_to_refinement(false)
113 : #ifdef LIBMESH_ENABLE_PERIODIC
114 318155 : , _periodic_boundaries(nullptr)
115 : #endif
116 : {
117 309181 : }
118 :
119 :
120 :
121 : #ifdef LIBMESH_ENABLE_PERIODIC
122 10 : void MeshRefinement::set_periodic_boundaries_ptr(PeriodicBoundaries * pb_ptr)
123 : {
124 10 : _periodic_boundaries = pb_ptr;
125 10 : }
126 : #endif
127 :
128 :
129 :
130 292700 : MeshRefinement::~MeshRefinement () = default;
131 :
132 :
133 :
134 62845 : void MeshRefinement::clear ()
135 : {
136 2096 : _new_nodes_map.clear();
137 62845 : }
138 :
139 :
140 :
141 47447129 : Node * MeshRefinement::add_node(Elem & parent,
142 : unsigned int child,
143 : unsigned int node,
144 : processor_id_type proc_id)
145 : {
146 5342560 : LOG_SCOPE("add_node()", "MeshRefinement");
147 :
148 47447129 : unsigned int parent_n = parent.as_parent_node(child, node);
149 :
150 47447129 : if (parent_n != libMesh::invalid_uint)
151 14918018 : return parent.node_ptr(parent_n);
152 :
153 : const std::vector<std::pair<dof_id_type, dof_id_type>>
154 35200407 : bracketing_nodes = parent.bracketing_nodes(child, node);
155 :
156 : // If we're not a parent node, we *must* be bracketed by at least
157 : // one pair of parent nodes
158 1902508 : libmesh_assert(bracketing_nodes.size());
159 :
160 : // Return the node if it already exists.
161 : //
162 : // We'll leave the processor_id untouched in this case - if we're
163 : // repartitioning later or if this is a new unpartitioned node,
164 : // we'll update it then, and if not then we don't want to update it.
165 33297899 : if (const auto new_node_id = _new_nodes_map.find(bracketing_nodes);
166 : new_node_id != DofObject::invalid_id)
167 23762132 : return _mesh.node_ptr(new_node_id);
168 :
169 : // Otherwise we need to add a new node.
170 : //
171 : // Figure out where to add the point:
172 :
173 525165 : Point p; // defaults to 0,0,0
174 :
175 116803995 : for (auto n : parent.node_index_range())
176 : {
177 : // The value from the embedding matrix
178 107268228 : const Real em_val = parent.embedding_matrix(child,node,n);
179 :
180 107268228 : if (em_val != 0.)
181 : {
182 53606029 : p.add_scaled (parent.point(n), em_val);
183 :
184 : // If we'd already found the node we shouldn't be here
185 2739372 : libmesh_assert_not_equal_to (em_val, 1);
186 : }
187 : }
188 :
189 : // Although we're leaving new nodes unpartitioned at first, with a
190 : // DistributedMesh we would need a default id based on the numbering
191 : // scheme for the requested processor_id.
192 9535767 : Node * new_node = _mesh.add_point (p, DofObject::invalid_id, proc_id);
193 :
194 525165 : libmesh_assert(new_node);
195 :
196 : // But then we'll make sure this node is marked as unpartitioned.
197 9535767 : new_node->processor_id() = DofObject::invalid_processor_id;
198 :
199 : // Add the node to the map.
200 9535767 : _new_nodes_map.add_node(*new_node, bracketing_nodes);
201 :
202 : // Return the address of the new node
203 525165 : return new_node;
204 : }
205 :
206 :
207 :
208 0 : Elem * MeshRefinement::add_elem (Elem * elem)
209 : {
210 0 : libmesh_assert(elem);
211 0 : return _mesh.add_elem (elem);
212 : }
213 :
214 6534077 : Elem * MeshRefinement::add_elem (std::unique_ptr<Elem> elem)
215 : {
216 389786 : libmesh_assert(elem);
217 7313649 : return _mesh.add_elem(std::move(elem));
218 : }
219 :
220 :
221 1590 : void MeshRefinement::create_parent_error_vector(const ErrorVector & error_per_cell,
222 : ErrorVector & error_per_parent,
223 : Real & parent_error_min,
224 : Real & parent_error_max)
225 : {
226 : // This function must be run on all processors at once
227 82 : parallel_object_only();
228 :
229 : // Make sure the input error vector is valid
230 : #ifdef DEBUG
231 11662 : for (const auto & val : error_per_cell)
232 : {
233 11580 : libmesh_assert_greater_equal (val, 0);
234 : // isnan() isn't standard C++ yet
235 : #ifdef isnan
236 : libmesh_assert(!isnan(val));
237 : #endif
238 : }
239 :
240 : // Use a reference to std::vector to avoid confusing
241 : // this->comm().verify
242 82 : std::vector<ErrorVectorReal> & epc = error_per_parent;
243 82 : libmesh_assert(this->comm().verify(epc));
244 : #endif // #ifdef DEBUG
245 :
246 : // error values on uncoarsenable elements will be left at -1
247 1590 : error_per_parent.clear();
248 1672 : error_per_parent.resize(error_per_cell.size(), 0.0);
249 :
250 : {
251 : // Find which elements are uncoarsenable
252 79017 : for (auto & elem : _mesh.active_local_element_ptr_range())
253 : {
254 42475 : Elem * parent = elem->parent();
255 :
256 : // Active elements are uncoarsenable
257 46986 : error_per_parent[elem->id()] = -1.0;
258 :
259 : // Grandparents and up are uncoarsenable
260 121879 : while (parent)
261 : {
262 22546 : parent = parent->parent();
263 107676 : if (parent)
264 : {
265 14563 : const dof_id_type parentid = parent->id();
266 7292 : libmesh_assert_less (parentid, error_per_parent.size());
267 78692 : error_per_parent[parentid] = -1.0;
268 : }
269 : }
270 1426 : }
271 :
272 : // Sync between processors.
273 : // Use a reference to std::vector to avoid confusing
274 : // this->comm().min
275 82 : std::vector<ErrorVectorReal> & epp = error_per_parent;
276 1590 : this->comm().min(epp);
277 : }
278 :
279 : // The parent's error is defined as the square root of the
280 : // sum of the children's errors squared, so errors that are
281 : // Hilbert norms remain Hilbert norms.
282 : //
283 : // Because the children may be on different processors, we
284 : // calculate local contributions to the parents' errors squared
285 : // first, then sum across processors and take the square roots
286 : // second.
287 : Threads::parallel_for
288 1590 : (_mesh.active_local_element_stored_range(),
289 50935 : [&error_per_cell, &error_per_parent](const ConstElemRange & range)
290 : {
291 44287 : for (const Elem * elem : range)
292 : {
293 9031 : const Elem * parent = elem->parent();
294 :
295 : // Calculate each contribution to parent cells
296 42475 : if (parent)
297 : {
298 7983 : const dof_id_type parentid = parent->id();
299 3996 : libmesh_assert_less (parentid, error_per_parent.size());
300 :
301 : // If the parent has grandchildren we won't be able to
302 : // coarsen it, so forget it. Otherwise, add this child's
303 : // contribution to the sum of the squared child errors
304 40242 : if (error_per_parent[parentid] != -1.0)
305 38448 : error_per_parent[parentid] += (error_per_cell[elem->id()] *
306 3468 : error_per_cell[elem->id()]);
307 : }
308 : }
309 1812 : });
310 :
311 : // Sum the vector across all processors
312 1590 : this->comm().sum(static_cast<std::vector<ErrorVectorReal> &>(error_per_parent));
313 :
314 : // Calculate the min and max as we loop
315 1590 : parent_error_min = std::numeric_limits<double>::max();
316 1590 : parent_error_max = 0.;
317 :
318 302904 : for (auto i : index_range(error_per_parent))
319 : {
320 : // If this element isn't a coarsenable parent with error, we
321 : // have nothing to do. Just flag it as -1 and move on
322 : // Note that this->comm().sum might have left uncoarsenable
323 : // elements with error_per_parent=-n_proc, so reset it to
324 : // error_per_parent=-1
325 312870 : if (error_per_parent[i] < 0.)
326 : {
327 257326 : error_per_parent[i] = -1.;
328 257326 : continue;
329 : }
330 :
331 : // The error estimator might have already given us an
332 : // estimate on the coarsenable parent elements; if so then
333 : // we want to retain that estimate
334 45738 : if (error_per_cell[i])
335 : {
336 0 : error_per_parent[i] = error_per_cell[i];
337 0 : continue;
338 : }
339 : // if not, then e_parent = sqrt(sum(e_child^2))
340 : else
341 43988 : error_per_parent[i] = std::sqrt(error_per_parent[i]);
342 :
343 43988 : parent_error_min = std::min (parent_error_min,
344 1750 : error_per_parent[i]);
345 50695 : parent_error_max = std::max (parent_error_max,
346 1750 : error_per_parent[i]);
347 : }
348 1590 : }
349 :
350 :
351 :
352 62845 : void MeshRefinement::update_nodes_map ()
353 : {
354 62845 : this->_new_nodes_map.init(_mesh);
355 62845 : }
356 :
357 :
358 :
359 1396 : bool MeshRefinement::test_level_one (bool libmesh_dbg_var(libmesh_assert_pass)) const
360 : {
361 : // This function must be run on all processors at once
362 1396 : parallel_object_only();
363 :
364 : // We may need a PointLocator for topological_neighbor() tests
365 : // later, which we need to make sure gets constructed on all
366 : // processors at once.
367 1396 : std::unique_ptr<PointLocatorBase> point_locator;
368 :
369 : #ifdef LIBMESH_ENABLE_PERIODIC
370 : bool has_periodic_boundaries =
371 1396 : _periodic_boundaries && !_periodic_boundaries->empty();
372 1396 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
373 :
374 1396 : if (has_periodic_boundaries)
375 200 : point_locator = _mesh.sub_point_locator();
376 : #endif
377 :
378 1396 : bool failure = false;
379 :
380 : #ifndef NDEBUG
381 1396 : Elem * failed_elem = nullptr;
382 1396 : Elem * failed_neighbor = nullptr;
383 : #endif // !NDEBUG
384 :
385 281267 : for (auto & elem : _mesh.active_local_element_ptr_range())
386 1356551 : for (auto n : elem->side_index_range())
387 : {
388 : Elem * neighbor =
389 1076680 : topological_neighbor(elem, point_locator.get(), n);
390 :
391 2047895 : if (!neighbor || !neighbor->active() ||
392 971215 : neighbor == remote_elem)
393 105465 : continue;
394 :
395 971215 : if ((neighbor->level() + 1 < elem->level()) ||
396 1942430 : (neighbor->p_level() + 1 < elem->p_level()) ||
397 971215 : (neighbor->p_level() > elem->p_level() + 1))
398 : {
399 0 : failure = true;
400 : #ifndef NDEBUG
401 0 : failed_elem = elem;
402 0 : failed_neighbor = neighbor;
403 : #endif // !NDEBUG
404 0 : break;
405 : }
406 0 : }
407 :
408 : // If any processor failed, we failed globally
409 1396 : this->comm().max(failure);
410 :
411 1396 : if (failure)
412 : {
413 : // We didn't pass the level one test, so libmesh_assert that
414 : // we're allowed not to
415 : #ifndef NDEBUG
416 0 : if (libmesh_assert_pass)
417 : {
418 0 : libMesh::out << "MeshRefinement Level one failure, element: "
419 0 : << *failed_elem
420 0 : << std::endl;
421 0 : libMesh::out << "MeshRefinement Level one failure, neighbor: "
422 0 : << *failed_neighbor
423 0 : << std::endl;
424 : }
425 : #endif // !NDEBUG
426 0 : libmesh_assert(!libmesh_assert_pass);
427 0 : return false;
428 : }
429 1396 : return true;
430 0 : }
431 :
432 :
433 :
434 654 : bool MeshRefinement::test_unflagged (bool libmesh_dbg_var(libmesh_assert_pass)) const
435 : {
436 : // This function must be run on all processors at once
437 654 : parallel_object_only();
438 :
439 654 : bool found_flag = false;
440 :
441 : #ifndef NDEBUG
442 654 : Elem * failed_elem = nullptr;
443 : #endif
444 :
445 : // Search for local flags
446 133452 : for (auto & elem : _mesh.active_local_element_ptr_range())
447 132798 : if (elem->refinement_flag() == Elem::REFINE ||
448 132798 : elem->refinement_flag() == Elem::COARSEN ||
449 398394 : elem->p_refinement_flag() == Elem::REFINE ||
450 132798 : elem->p_refinement_flag() == Elem::COARSEN)
451 : {
452 0 : found_flag = true;
453 : #ifndef NDEBUG
454 0 : failed_elem = elem;
455 : #endif
456 0 : break;
457 0 : }
458 :
459 : // If we found a flag on any processor, it counts
460 654 : this->comm().max(found_flag);
461 :
462 654 : if (found_flag)
463 : {
464 : #ifndef NDEBUG
465 0 : if (libmesh_assert_pass)
466 : {
467 : libMesh::out <<
468 0 : "MeshRefinement test_unflagged failure, element: " <<
469 0 : *failed_elem << std::endl;
470 : }
471 : #endif
472 : // We didn't pass the "elements are unflagged" test,
473 : // so libmesh_assert that we're allowed not to
474 0 : libmesh_assert(!libmesh_assert_pass);
475 0 : return false;
476 : }
477 654 : return true;
478 : }
479 :
480 :
481 :
482 19406 : bool MeshRefinement::refine_and_coarsen_elements ()
483 : {
484 : // This function must be run on all processors at once
485 682 : parallel_object_only();
486 :
487 : // We can't yet turn a non-level-one mesh into a level-one mesh
488 682 : if (_face_level_mismatch_limit)
489 682 : libmesh_assert(test_level_one(true));
490 :
491 : // Possibly clean up the refinement flags from
492 : // a previous step. While we're at it, see if this method should be
493 : // a no-op.
494 19406 : std::atomic<bool> elements_flagged{false};
495 :
496 : Threads::parallel_for
497 19406 : (_mesh.element_stored_range(),
498 6356117 : [&elements_flagged](const ElemRange & range)
499 : {
500 8380006 : for (Elem * elem : range)
501 : {
502 : // This might be left over from the last step
503 344310 : const Elem::RefinementState flag = elem->refinement_flag();
504 :
505 : // Set refinement flag to INACTIVE if the
506 : // element isn't active
507 344310 : if (!elem->active())
508 : {
509 84336 : elem->set_refinement_flag(Elem::INACTIVE);
510 84336 : elem->set_p_refinement_flag(Elem::INACTIVE);
511 : }
512 6283757 : else if (flag == Elem::JUST_REFINED)
513 0 : elem->set_refinement_flag(Elem::DO_NOTHING);
514 6283757 : else if (!elements_flagged)
515 : {
516 45452 : if (flag == Elem::REFINE || flag == Elem::COARSEN)
517 606 : elements_flagged = true;
518 : else
519 : {
520 : const Elem::RefinementState pflag =
521 1681 : elem->p_refinement_flag();
522 30685 : if (pflag == Elem::REFINE || pflag == Elem::COARSEN)
523 51 : elements_flagged = true;
524 : }
525 : }
526 : }
527 20820 : });
528 :
529 : // Did *any* processor find elements flagged for AMR/C?
530 19406 : bool any_elements_flagged = elements_flagged;
531 19406 : _mesh.comm().max(any_elements_flagged);
532 :
533 : // If we have nothing to do, let's not bother verifying that nothing
534 : // is compatible with nothing.
535 19406 : if (!any_elements_flagged)
536 28 : return false;
537 :
538 : // Parallel consistency has to come first, or coarsening
539 : // along processor boundaries might occasionally be falsely
540 : // prevented
541 : #ifdef DEBUG
542 654 : bool flags_were_consistent = this->make_flags_parallel_consistent();
543 :
544 654 : libmesh_assert (flags_were_consistent);
545 : #endif
546 :
547 : // Smooth refinement and coarsening flags
548 18412 : _smooth_flags(true, true);
549 :
550 : // First coarsen the flagged elements.
551 : const bool coarsening_changed_mesh =
552 18412 : this->_coarsen_elements ();
553 :
554 : // First coarsen the flagged elements.
555 : // FIXME: test_level_one now tests consistency across periodic
556 : // boundaries, which requires a point_locator, which just got
557 : // invalidated by _coarsen_elements() and hasn't yet been cleared by
558 : // prepare_for_use().
559 :
560 : // libmesh_assert(this->make_coarsening_compatible());
561 : // libmesh_assert(this->make_refinement_compatible());
562 :
563 : // FIXME: This won't pass unless we add a redundant find_neighbors()
564 : // call or replace find_neighbors() with on-the-fly neighbor updating
565 : // libmesh_assert(!this->eliminate_unrefined_patches());
566 :
567 : // We can't contract the mesh ourselves anymore - a System might
568 : // need to restrict old coefficient vectors first
569 : // _mesh.contract();
570 :
571 : // First coarsen the flagged elements.
572 : // Now refine the flagged elements. This will
573 : // take up some space, maybe more than what was freed.
574 : const bool refining_changed_mesh =
575 18412 : this->_refine_elements();
576 :
577 : // First coarsen the flagged elements.
578 : // Finally, the new mesh needs to be prepared for use
579 18412 : if (coarsening_changed_mesh || refining_changed_mesh)
580 : {
581 : #ifdef DEBUG
582 646 : _mesh.libmesh_assert_valid_parallel_ids();
583 : #endif
584 :
585 18132 : _mesh.prepare_for_use ();
586 :
587 646 : if (_face_level_mismatch_limit)
588 646 : libmesh_assert(test_level_one(true));
589 646 : libmesh_assert(test_unflagged(true));
590 646 : libmesh_assert(this->make_coarsening_compatible());
591 646 : libmesh_assert(this->make_refinement_compatible());
592 : // FIXME: This won't pass unless we add a redundant find_neighbors()
593 : // call or replace find_neighbors() with on-the-fly neighbor updating
594 : // libmesh_assert(!this->eliminate_unrefined_patches());
595 :
596 18132 : return true;
597 : }
598 : else
599 : {
600 8 : if (_face_level_mismatch_limit)
601 8 : libmesh_assert(test_level_one(true));
602 8 : libmesh_assert(test_unflagged(true));
603 8 : libmesh_assert(this->make_coarsening_compatible());
604 8 : libmesh_assert(this->make_refinement_compatible());
605 : }
606 :
607 : // Otherwise there was no change in the mesh,
608 : // let the user know. Also, there is no need
609 : // to prepare the mesh for use since it did not change.
610 8 : return false;
611 :
612 : }
613 :
614 :
615 :
616 :
617 :
618 :
619 :
620 25688 : bool MeshRefinement::coarsen_elements ()
621 : {
622 : // This function must be run on all processors at once
623 854 : parallel_object_only();
624 :
625 : // We can't yet turn a non-level-one mesh into a level-one mesh
626 854 : if (_face_level_mismatch_limit)
627 4 : libmesh_assert(test_level_one(true));
628 :
629 : // Possibly clean up the refinement flags from
630 : // a previous step
631 : Threads::parallel_for
632 25688 : (_mesh.element_stored_range(),
633 26066 : [](const ElemRange & range)
634 : {
635 10492487 : for (Elem * elem : range)
636 : {
637 : // Set refinement flag to INACTIVE if the
638 : // element isn't active
639 1083556 : if (!elem->active())
640 : {
641 184516 : elem->set_refinement_flag(Elem::INACTIVE);
642 184516 : elem->set_p_refinement_flag(Elem::INACTIVE);
643 : }
644 :
645 : // This might be left over from the last step
646 10464885 : if (elem->refinement_flag() == Elem::JUST_REFINED)
647 135804 : elem->set_refinement_flag(Elem::DO_NOTHING);
648 : }
649 26066 : });
650 :
651 : // Parallel consistency has to come first, or coarsening
652 : // along processor boundaries might occasionally be falsely
653 : // prevented
654 25688 : bool flags_were_consistent = this->make_flags_parallel_consistent();
655 :
656 : // In theory, we should be able to remove the above call, which can
657 : // be expensive and should be unnecessary. In practice, doing
658 : // consistent flagging in parallel is hard, it's impossible to
659 : // verify at the library level if it's being done by user code, and
660 : // we don't want to abort large parallel runs in opt mode... but we
661 : // do want to warn that they should be fixed.
662 854 : libmesh_assert(flags_were_consistent);
663 :
664 854 : if (!flags_were_consistent)
665 : libmesh_warning("Warning: Refinement flags were not consistent between processors! "
666 : "Correcting and continuing.\n");
667 :
668 : // Smooth coarsening flags
669 25688 : _smooth_flags(false, true);
670 :
671 : // Coarsen the flagged elements.
672 : const bool mesh_changed =
673 25688 : this->_coarsen_elements ();
674 :
675 854 : if (_face_level_mismatch_limit)
676 4 : libmesh_assert(test_level_one(true));
677 854 : libmesh_assert(this->make_coarsening_compatible());
678 : // FIXME: This won't pass unless we add a redundant find_neighbors()
679 : // call or replace find_neighbors() with on-the-fly neighbor updating
680 : // libmesh_assert(!this->eliminate_unrefined_patches());
681 :
682 : // We can't contract the mesh ourselves anymore - a System might
683 : // need to restrict old coefficient vectors first
684 : // _mesh.contract();
685 :
686 : // Finally, the new mesh may need to be prepared for use
687 25688 : if (mesh_changed)
688 142 : _mesh.prepare_for_use ();
689 :
690 25688 : return mesh_changed;
691 : }
692 :
693 :
694 :
695 :
696 :
697 :
698 :
699 26603 : bool MeshRefinement::refine_elements ()
700 : {
701 : // This function must be run on all processors at once
702 878 : parallel_object_only();
703 :
704 : // Prevent refinement when the mesh has disjoint neighbor boundary pairs.
705 : //
706 : // Disjoint neighbor boundary interfaces violate the topological continuity
707 : // assumed by the refinement algorithm. Refinement on such meshes can
708 : // produce invalid neighbor relationships (for example reverse indices
709 : // equal to invalid_uint), which may trigger assertions like
710 : // `rev < neigh->n_neighbors()` or corrupt neighbor data.
711 : //
712 : // For safety, refinement is disabled while disjoint neighbor boundary pairs are
713 : // present.
714 26603 : if (_mesh.get_disjoint_neighbor_boundary_pairs())
715 : {
716 211 : libmesh_not_implemented_msg(
717 : "Mesh refinement is not yet implemented for meshes with disjoint neighbor boundary interfaces.\n");
718 : }
719 :
720 876 : if (_face_level_mismatch_limit)
721 26 : libmesh_assert(test_level_one(true));
722 :
723 : // Possibly clean up the refinement flags from
724 : // a previous step
725 : Threads::parallel_for
726 26532 : (_mesh.element_stored_range(),
727 26918 : [](const ElemRange & range)
728 : {
729 9096948 : for (Elem * elem : range)
730 : {
731 : // Set refinement flag to INACTIVE if the
732 : // element isn't active
733 961564 : if (!elem->active())
734 : {
735 113024 : elem->set_refinement_flag(Elem::INACTIVE);
736 113024 : elem->set_p_refinement_flag(Elem::INACTIVE);
737 : }
738 :
739 : // This might be left over from the last step
740 9068458 : if (elem->refinement_flag() == Elem::JUST_REFINED)
741 6 : elem->set_refinement_flag(Elem::DO_NOTHING);
742 : }
743 26918 : });
744 :
745 : // Parallel consistency has to come first, or coarsening
746 : // along processor boundaries might occasionally be falsely
747 : // prevented
748 26532 : bool flags_were_consistent = this->make_flags_parallel_consistent();
749 :
750 : // In theory, we should be able to remove the above call, which can
751 : // be expensive and should be unnecessary. In practice, doing
752 : // consistent flagging in parallel is hard, it's impossible to
753 : // verify at the library level if it's being done by user code, and
754 : // we don't want to abort large parallel runs in opt mode... but we
755 : // do want to warn that they should be fixed.
756 876 : libmesh_assert(flags_were_consistent);
757 :
758 876 : if (!flags_were_consistent)
759 : libmesh_warning("Warning: Refinement flags were not consistent between processors! "
760 : "Correcting and continuing.\n");
761 :
762 : // Smooth refinement flags
763 26532 : _smooth_flags(true, false);
764 :
765 : // Now refine the flagged elements. This will
766 : // take up some space, maybe more than what was freed.
767 : const bool mesh_changed =
768 26532 : this->_refine_elements();
769 :
770 876 : if (_face_level_mismatch_limit)
771 26 : libmesh_assert(test_level_one(true));
772 876 : libmesh_assert(this->make_refinement_compatible());
773 :
774 : // FIXME: This won't pass unless we add a redundant find_neighbors()
775 : // call or replace find_neighbors() with on-the-fly neighbor updating
776 : // libmesh_assert(!this->eliminate_unrefined_patches());
777 :
778 : // Finally, the new mesh needs to be prepared for use
779 26532 : if (mesh_changed)
780 2193 : _mesh.prepare_for_use ();
781 :
782 26532 : return mesh_changed;
783 : }
784 :
785 :
786 :
787 96707 : bool MeshRefinement::make_flags_parallel_consistent()
788 : {
789 : // This function must be run on all processors at once
790 2426 : parallel_object_only();
791 :
792 2426 : LOG_SCOPE ("make_flags_parallel_consistent()", "MeshRefinement");
793 :
794 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
795 96707 : &Elem::set_refinement_flag);
796 : Parallel::sync_dofobject_data_by_id
797 190988 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), hsync);
798 :
799 : SyncRefinementFlags psync(_mesh, &Elem::p_refinement_flag,
800 96707 : &Elem::set_p_refinement_flag);
801 : Parallel::sync_dofobject_data_by_id
802 190988 : (this->comm(), _mesh.elements_begin(), _mesh.elements_end(), psync);
803 :
804 : // If we weren't consistent in both h and p on every processor then
805 : // we weren't globally consistent
806 191402 : bool parallel_consistent = hsync.parallel_consistent &&
807 96456 : psync.parallel_consistent;
808 96707 : this->comm().min(parallel_consistent);
809 :
810 99133 : return parallel_consistent;
811 : }
812 :
813 :
814 :
815 53986 : bool MeshRefinement::make_coarsening_compatible()
816 : {
817 : // This function must be run on all processors at once
818 3210 : parallel_object_only();
819 :
820 : // We may need a PointLocator for topological_neighbor() tests
821 : // later, which we need to make sure gets constructed on all
822 : // processors at once.
823 52282 : std::unique_ptr<PointLocatorBase> point_locator;
824 :
825 : #ifdef LIBMESH_ENABLE_PERIODIC
826 : bool has_periodic_boundaries =
827 53986 : _periodic_boundaries && !_periodic_boundaries->empty();
828 3210 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
829 :
830 3210 : if (has_periodic_boundaries)
831 532 : point_locator = _mesh.sub_point_locator();
832 : #endif
833 :
834 6420 : LOG_SCOPE ("make_coarsening_compatible()", "MeshRefinement");
835 :
836 : // Unless we encounter a specific situation level-one
837 : // will be satisfied after executing this loop just once
838 3210 : bool level_one_satisfied = true;
839 :
840 :
841 : // Unless we encounter a specific situation we will be compatible
842 : // with any selected refinement flags
843 53986 : bool compatible_with_refinement = true;
844 :
845 :
846 : // find the maximum h and p levels in the mesh
847 53986 : unsigned int max_level = 0;
848 53986 : unsigned int max_p_level = 0;
849 :
850 : {
851 : // First we look at all the active level-0 elements. Since it doesn't make
852 : // sense to coarsen them we must un-set their coarsen flags if
853 : // they are set.
854 31750946 : for (auto & elem : _mesh.active_element_ptr_range())
855 : {
856 16867270 : max_level = std::max(max_level, elem->level());
857 16867270 : max_p_level =
858 16867270 : std::max(max_p_level,
859 17600108 : static_cast<unsigned int>(elem->p_level()));
860 :
861 16880814 : if ((elem->level() == 0) &&
862 123463 : (elem->refinement_flag() == Elem::COARSEN))
863 1040 : elem->set_refinement_flag(Elem::DO_NOTHING);
864 :
865 18217644 : if ((elem->p_level() == 0) &&
866 1350374 : (elem->p_refinement_flag() == Elem::COARSEN))
867 10 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
868 49072 : }
869 : }
870 :
871 : // Even if there are no refined elements on this processor then
872 : // there may still be work for us to do on e.g. ancestor elements.
873 : // At the very least we need to be in the loop if a distributed mesh
874 : // needs to synchronize data.
875 : #if 0
876 : if (max_level == 0 && max_p_level == 0)
877 : {
878 : // But we still have to check with other processors
879 : this->comm().min(compatible_with_refinement);
880 :
881 : return compatible_with_refinement;
882 : }
883 : #endif
884 :
885 : // Loop over all the active elements. If an element is marked
886 : // for coarsening we better check its neighbors. If ANY of these neighbors
887 : // are marked for refinement AND are at the same level then there is a
888 : // conflict. By convention refinement wins, so we un-mark the element for
889 : // coarsening. Level-one would be violated in this case so we need to re-run
890 : // the loop.
891 53986 : if (_face_level_mismatch_limit)
892 : {
893 :
894 27519 : repeat:
895 1882 : level_one_satisfied = true;
896 :
897 376 : do
898 : {
899 2258 : level_one_satisfied = true;
900 :
901 40688356 : for (auto & elem : _mesh.active_element_ptr_range())
902 : {
903 1095864 : bool my_flag_changed = false;
904 :
905 22092744 : if (elem->refinement_flag() == Elem::COARSEN) // If the element is active and
906 : // the coarsen flag is set
907 : {
908 5088325 : const unsigned int my_level = elem->level();
909 :
910 21745124 : for (auto n : elem->side_index_range())
911 : {
912 : const Elem * neighbor =
913 16703126 : topological_neighbor(elem, point_locator.get(), n);
914 :
915 16703126 : if (neighbor != nullptr && // I have a
916 16373408 : neighbor != remote_elem) // neighbor here
917 : {
918 1266264 : if (neighbor->active()) // and it is active
919 : {
920 16692478 : if ((neighbor->level() == my_level) &&
921 543304 : (neighbor->refinement_flag() == Elem::REFINE)) // the neighbor is at my level
922 : // and wants to be refined
923 : {
924 24433 : elem->set_refinement_flag(Elem::DO_NOTHING);
925 1104 : my_flag_changed = true;
926 24433 : break;
927 : }
928 : }
929 : else // I have a neighbor and it is not active. That means it has children.
930 : { // While it _may_ be possible to coarsen us if all the children of
931 : // that element want to be coarsened, it is impossible to know at this
932 : // stage. Forget about it for the moment... This can be handled in
933 : // two steps.
934 208836 : elem->set_refinement_flag(Elem::DO_NOTHING);
935 8900 : my_flag_changed = true;
936 208836 : break;
937 : }
938 : }
939 : }
940 : }
941 22092744 : if (elem->p_refinement_flag() == Elem::COARSEN) // If
942 : // the element is active and the order reduction flag is set
943 : {
944 164 : const unsigned int my_p_level = elem->p_level();
945 :
946 7129 : for (auto n : elem->side_index_range())
947 : {
948 : const Elem * neighbor =
949 5896 : topological_neighbor(elem, point_locator.get(), n);
950 :
951 5896 : if (neighbor != nullptr && // I have a
952 4637 : neighbor != remote_elem) // neighbor here
953 : {
954 388 : if (neighbor->active()) // and it is active
955 : {
956 348 : if ((neighbor->p_level() > my_p_level &&
957 0 : neighbor->p_refinement_flag() != Elem::COARSEN)
958 3772 : || (neighbor->p_level() == my_p_level &&
959 116 : neighbor->p_refinement_flag() == Elem::REFINE))
960 : {
961 217 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
962 14 : my_flag_changed = true;
963 14 : break;
964 : }
965 : }
966 : else // I have a neighbor and it is not active.
967 : { // We need to find which of its children
968 : // have me as a neighbor, and maintain
969 : // level one p compatibility with them.
970 : // Because we currently have level one h
971 : // compatibility, we don't need to check
972 : // grandchildren
973 :
974 20 : libmesh_assert(neighbor->has_children());
975 1616 : for (auto & subneighbor : neighbor->child_ref_range())
976 1513 : if (&subneighbor != remote_elem &&
977 2362 : subneighbor.active() &&
978 973 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
979 598 : if ((subneighbor.p_level() > my_p_level &&
980 34 : subneighbor.p_refinement_flag() != Elem::COARSEN)
981 1014 : || (subneighbor.p_level() == my_p_level &&
982 0 : subneighbor.p_refinement_flag() == Elem::REFINE))
983 : {
984 194 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
985 12 : my_flag_changed = true;
986 12 : break;
987 : }
988 239 : if (my_flag_changed)
989 20 : break;
990 : }
991 : }
992 : }
993 : }
994 :
995 : // If the current element's flag changed, we hadn't
996 : // satisfied the level one rule.
997 21262079 : if (my_flag_changed)
998 10012 : level_one_satisfied = false;
999 :
1000 : // Additionally, if it has non-local neighbors, and
1001 : // we're not in serial, then we'll eventually have to
1002 : // return compatible_with_refinement = false, because
1003 : // our change has to propagate to neighboring
1004 : // processors.
1005 2139682 : if (my_flag_changed && !_mesh.is_serial())
1006 15948 : for (auto n : elem->side_index_range())
1007 : {
1008 : Elem * neigh =
1009 14881 : topological_neighbor(elem, point_locator.get(), n);
1010 :
1011 14881 : if (!neigh)
1012 2930 : continue;
1013 11951 : if (neigh == remote_elem ||
1014 10995 : neigh->processor_id() !=
1015 0 : this->processor_id())
1016 : {
1017 5339 : compatible_with_refinement = false;
1018 5339 : break;
1019 : }
1020 : // FIXME - for non-level one meshes we should
1021 : // test all descendants
1022 0 : if (neigh->has_children())
1023 2074 : for (auto & child : neigh->child_ref_range())
1024 1732 : if (&child == remote_elem ||
1025 1732 : child.processor_id() !=
1026 0 : this->processor_id())
1027 : {
1028 258 : compatible_with_refinement = false;
1029 258 : break;
1030 : }
1031 : }
1032 42328 : }
1033 : }
1034 46188 : while (!level_one_satisfied);
1035 :
1036 : } // end if (_face_level_mismatch_limit)
1037 :
1038 :
1039 : // Next we look at all of the ancestor cells.
1040 : // If there is a parent cell with all of its children
1041 : // wanting to be unrefined then the element is a candidate
1042 : // for unrefinement. If all the children don't
1043 : // all want to be unrefined then ALL of them need to have their
1044 : // unrefinement flags cleared.
1045 328082 : for (int level = max_level; level >= 0; level--)
1046 61109765 : for (auto & elem : as_range(_mesh.level_elements_begin(level), _mesh.level_elements_end(level)))
1047 31967598 : if (elem->ancestor())
1048 : {
1049 : // right now the element hasn't been disqualified
1050 : // as a candidate for unrefinement
1051 509880 : bool is_a_candidate = true;
1052 509880 : bool found_remote_child = false;
1053 :
1054 38582515 : for (auto & child : elem->child_ref_range())
1055 : {
1056 30664344 : if (&child == remote_elem)
1057 40 : found_remote_child = true;
1058 30247083 : else if ((child.refinement_flag() != Elem::COARSEN) ||
1059 129338 : !child.active() )
1060 1939622 : is_a_candidate = false;
1061 : }
1062 :
1063 7604169 : if (!is_a_candidate && !found_remote_child)
1064 : {
1065 6703647 : elem->set_refinement_flag(Elem::INACTIVE);
1066 :
1067 33617355 : for (auto & child : elem->child_ref_range())
1068 : {
1069 26913708 : if (&child == remote_elem)
1070 0 : continue;
1071 26913708 : if (child.refinement_flag() == Elem::COARSEN)
1072 : {
1073 27338 : level_one_satisfied = false;
1074 27338 : child.set_refinement_flag(Elem::DO_NOTHING);
1075 : }
1076 : }
1077 : }
1078 238705 : }
1079 :
1080 63393 : if (!level_one_satisfied && _face_level_mismatch_limit) goto repeat;
1081 :
1082 :
1083 : // If all the children of a parent are set to be coarsened
1084 : // then flag the parent so that they can kill their kids.
1085 :
1086 : // On a distributed mesh, we won't always be able to determine this
1087 : // on parent elements with remote children, even if we own the
1088 : // parent, without communication.
1089 : //
1090 : // We'll first communicate *to* parents' owners when we determine
1091 : // they cannot be coarsened, then we'll sync the final refinement
1092 : // flag *from* the parents.
1093 :
1094 : // uncoarsenable_parents[p] live on processor id p
1095 53986 : const processor_id_type n_proc = _mesh.n_processors();
1096 3210 : const processor_id_type my_proc_id = _mesh.processor_id();
1097 53986 : const bool distributed_mesh = !_mesh.is_replicated();
1098 :
1099 : std::vector<std::vector<dof_id_type>>
1100 53986 : uncoarsenable_parents(n_proc);
1101 :
1102 10674576 : for (auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1103 : {
1104 : // Presume all the children are flagged for coarsening and
1105 : // then look for a contradiction
1106 430684 : bool all_children_flagged_for_coarsening = true;
1107 :
1108 8854969 : for (auto & child : elem->child_ref_range())
1109 : {
1110 8396606 : if (&child != remote_elem &&
1111 470856 : child.refinement_flag() != Elem::COARSEN)
1112 : {
1113 417496 : all_children_flagged_for_coarsening = false;
1114 5984831 : if (!distributed_mesh)
1115 417344 : break;
1116 1184150 : if (child.processor_id() != elem->processor_id())
1117 : {
1118 97170 : uncoarsenable_parents[elem->processor_id()].push_back(elem->id());
1119 0 : break;
1120 : }
1121 : }
1122 : }
1123 :
1124 5495094 : if (all_children_flagged_for_coarsening)
1125 372141 : elem->set_refinement_flag(Elem::COARSEN_INACTIVE);
1126 : else
1127 5220123 : elem->set_refinement_flag(Elem::INACTIVE);
1128 49072 : }
1129 :
1130 : // If we have a distributed mesh, we might need to sync up
1131 : // INACTIVE vs. COARSEN_INACTIVE flags.
1132 53986 : if (distributed_mesh)
1133 : {
1134 : // We'd better still be in sync here
1135 8 : parallel_object_only();
1136 :
1137 : Parallel::MessageTag
1138 29436 : uncoarsenable_tag = this->comm().get_unique_tag();
1139 29424 : std::vector<Parallel::Request> uncoarsenable_push_requests(n_proc-1);
1140 :
1141 342314 : for (processor_id_type p = 0; p != n_proc; ++p)
1142 : {
1143 312890 : if (p == my_proc_id)
1144 29420 : continue;
1145 :
1146 : Parallel::Request &request =
1147 283466 : uncoarsenable_push_requests[p - (p > my_proc_id)];
1148 :
1149 283466 : _mesh.comm().send
1150 283466 : (p, uncoarsenable_parents[p], request, uncoarsenable_tag);
1151 : }
1152 :
1153 312890 : for (processor_id_type p = 1; p != n_proc; ++p)
1154 : {
1155 16 : std::vector<dof_id_type> my_uncoarsenable_parents;
1156 283466 : _mesh.comm().receive
1157 283466 : (Parallel::any_source, my_uncoarsenable_parents,
1158 12 : uncoarsenable_tag);
1159 :
1160 355529 : for (const auto & id : my_uncoarsenable_parents)
1161 : {
1162 72063 : Elem & elem = _mesh.elem_ref(id);
1163 0 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1164 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1165 0 : elem.set_refinement_flag(Elem::INACTIVE);
1166 : }
1167 : }
1168 :
1169 29424 : Parallel::wait(uncoarsenable_push_requests);
1170 :
1171 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1172 29424 : &Elem::set_refinement_flag);
1173 : sync_dofobject_data_by_id
1174 58840 : (this->comm(), _mesh.not_local_elements_begin(),
1175 29436 : _mesh.not_local_elements_end(),
1176 : // We'd like a smaller sync, but this leads to bugs?
1177 : // SyncCoarsenInactive(),
1178 : hsync);
1179 29412 : }
1180 :
1181 : // If one processor finds an incompatibility, we're globally
1182 : // incompatible
1183 53986 : this->comm().min(compatible_with_refinement);
1184 :
1185 107972 : return compatible_with_refinement;
1186 49072 : }
1187 :
1188 :
1189 :
1190 :
1191 :
1192 :
1193 :
1194 :
1195 54976 : bool MeshRefinement::make_refinement_compatible()
1196 : {
1197 : // This function must be run on all processors at once
1198 3256 : parallel_object_only();
1199 :
1200 : // We may need a PointLocator for topological_neighbor() tests
1201 : // later, which we need to make sure gets constructed on all
1202 : // processors at once.
1203 53248 : std::unique_ptr<PointLocatorBase> point_locator;
1204 :
1205 : #ifdef LIBMESH_ENABLE_PERIODIC
1206 : bool has_periodic_boundaries =
1207 54976 : _periodic_boundaries && !_periodic_boundaries->empty();
1208 3256 : libmesh_assert(this->comm().verify(has_periodic_boundaries));
1209 :
1210 3256 : if (has_periodic_boundaries)
1211 532 : point_locator = _mesh.sub_point_locator();
1212 : #endif
1213 :
1214 3256 : LOG_SCOPE ("make_refinement_compatible()", "MeshRefinement");
1215 :
1216 : // Unless we encounter a specific situation we will be compatible
1217 : // with any selected coarsening flags
1218 54976 : bool compatible_with_coarsening = true;
1219 :
1220 : // This loop enforces the level-1 rule. We should only
1221 : // execute it if the user indeed wants level-1 satisfied!
1222 54976 : if (_face_level_mismatch_limit)
1223 : {
1224 : // Unless we encounter a specific situation level-one
1225 : // will be satisfied after executing this loop just once
1226 1556 : bool level_one_satisfied = true;
1227 :
1228 220 : do
1229 : {
1230 1776 : level_one_satisfied = true;
1231 :
1232 26213812 : for (auto & elem : _mesh.active_element_ptr_range())
1233 : {
1234 13750658 : const unsigned short n_sides = elem->n_sides();
1235 :
1236 14282842 : if (elem->refinement_flag() == Elem::REFINE) // If the element is active and the
1237 : // h refinement flag is set
1238 : {
1239 456932 : const unsigned int my_level = elem->level();
1240 :
1241 2217714 : for (unsigned short side = 0; side != n_sides;
1242 : ++side)
1243 : {
1244 : Elem * neighbor =
1245 1760782 : topological_neighbor(elem, point_locator.get(), side);
1246 :
1247 1765404 : if (neighbor != nullptr && // I have a
1248 1843062 : neighbor != remote_elem && // neighbor here
1249 164694 : neighbor->active()) // and it is active
1250 : {
1251 : // Case 1: The neighbor is at the same level I am.
1252 : // 1a: The neighbor will be refined -> NO PROBLEM
1253 : // 1b: The neighbor won't be refined -> NO PROBLEM
1254 : // 1c: The neighbor wants to be coarsened -> PROBLEM
1255 1256082 : if (neighbor->level() == my_level)
1256 : {
1257 1135716 : if (neighbor->refinement_flag() == Elem::COARSEN)
1258 : {
1259 2 : neighbor->set_refinement_flag(Elem::DO_NOTHING);
1260 79 : if (neighbor->parent())
1261 2 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1262 77 : compatible_with_coarsening = false;
1263 2 : level_one_satisfied = false;
1264 : }
1265 : }
1266 :
1267 :
1268 : // Case 2: The neighbor is one level lower than I am.
1269 : // The neighbor thus MUST be refined to satisfy
1270 : // the level-one rule, regardless of whether it
1271 : // was originally flagged for refinement. If it
1272 : // wasn't flagged already we need to repeat
1273 : // this process.
1274 120366 : else if ((neighbor->level()+1) == my_level)
1275 : {
1276 120366 : if (neighbor->refinement_flag() != Elem::REFINE)
1277 : {
1278 624 : neighbor->set_refinement_flag(Elem::REFINE);
1279 13492 : if (neighbor->parent())
1280 464 : neighbor->parent()->set_refinement_flag(Elem::INACTIVE);
1281 12868 : compatible_with_coarsening = false;
1282 624 : level_one_satisfied = false;
1283 : }
1284 : }
1285 : #ifdef DEBUG
1286 : // Note that the only other possibility is that the
1287 : // neighbor is already refined, in which case it isn't
1288 : // active and we should never get here.
1289 : else
1290 0 : libmesh_error_msg("ERROR: Neighbor level must be equal or 1 higher than mine.");
1291 : #endif
1292 : }
1293 : }
1294 : }
1295 14282842 : if (elem->p_refinement_flag() == Elem::REFINE) // If the element is active and the
1296 : // p refinement flag is set
1297 : {
1298 2798 : const unsigned int my_p_level = elem->p_level();
1299 :
1300 129811 : for (unsigned int side=0; side != n_sides; side++)
1301 : {
1302 : Elem * neighbor =
1303 109588 : topological_neighbor(elem, point_locator.get(), side);
1304 :
1305 109588 : if (neighbor != nullptr && // I have a
1306 90693 : neighbor != remote_elem) // neighbor here
1307 : {
1308 11906 : if (neighbor->active()) // and it is active
1309 : {
1310 86119 : if (neighbor->p_level() < my_p_level &&
1311 18 : neighbor->p_refinement_flag() != Elem::REFINE)
1312 : {
1313 0 : neighbor->set_p_refinement_flag(Elem::REFINE);
1314 0 : level_one_satisfied = false;
1315 18 : compatible_with_coarsening = false;
1316 : }
1317 91929 : if (neighbor->p_level() == my_p_level &&
1318 5828 : neighbor->p_refinement_flag() == Elem::COARSEN)
1319 : {
1320 0 : neighbor->set_p_refinement_flag(Elem::DO_NOTHING);
1321 0 : level_one_satisfied = false;
1322 0 : compatible_with_coarsening = false;
1323 : }
1324 : }
1325 : else // I have an inactive neighbor
1326 : {
1327 62 : libmesh_assert(neighbor->has_children());
1328 1788 : for (auto & subneighbor : neighbor->child_ref_range())
1329 2179 : if (&subneighbor != remote_elem && subneighbor.active() &&
1330 987 : has_topological_neighbor(&subneighbor, point_locator.get(), elem))
1331 : {
1332 690 : if (subneighbor.p_level() < my_p_level &&
1333 44 : subneighbor.p_refinement_flag() != Elem::REFINE)
1334 : {
1335 : // We should already be level one
1336 : // compatible
1337 6 : libmesh_assert_greater (subneighbor.p_level() + 2u,
1338 : my_p_level);
1339 6 : subneighbor.set_p_refinement_flag(Elem::REFINE);
1340 6 : level_one_satisfied = false;
1341 64 : compatible_with_coarsening = false;
1342 : }
1343 664 : if (subneighbor.p_level() == my_p_level &&
1344 18 : subneighbor.p_refinement_flag() == Elem::COARSEN)
1345 : {
1346 0 : subneighbor.set_p_refinement_flag(Elem::DO_NOTHING);
1347 0 : level_one_satisfied = false;
1348 0 : compatible_with_coarsening = false;
1349 : }
1350 : }
1351 : }
1352 : }
1353 : }
1354 : }
1355 30843 : }
1356 : }
1357 :
1358 33719 : while (!level_one_satisfied);
1359 : } // end if (_face_level_mismatch_limit)
1360 :
1361 : // If we're not compatible on one processor, we're globally not
1362 : // compatible
1363 54976 : this->comm().min(compatible_with_coarsening);
1364 :
1365 59960 : return compatible_with_coarsening;
1366 49992 : }
1367 :
1368 :
1369 :
1370 :
1371 45232 : bool MeshRefinement::_coarsen_elements ()
1372 : {
1373 : // This function must be run on all processors at once
1374 1540 : parallel_object_only();
1375 :
1376 1540 : LOG_SCOPE ("_coarsen_elements()", "MeshRefinement");
1377 :
1378 : // Flags indicating if this call actually changes the mesh
1379 45232 : bool mesh_changed = false;
1380 45232 : bool mesh_p_changed = false;
1381 :
1382 : // Clear the unused_elements data structure.
1383 : // The elements have been packed since it was built,
1384 : // so there are _no_ unused elements. We cannot trust
1385 : // any iterators currently in this data structure.
1386 : // _unused_elements.clear();
1387 :
1388 : // Loop over the elements first to determine if the mesh will
1389 : // undergo h-coarsening. If it will, then we'll need to communicate
1390 : // more ghosted elements. We need to communicate them *before* we
1391 : // do the coarsening; otherwise it is possible to coarsen away a
1392 : // one-element-thick layer partition and leave the partitions on
1393 : // either side unable to figure out how to talk to each other.
1394 24785600 : for (auto & elem : _mesh.element_ptr_range())
1395 13658699 : if (elem->refinement_flag() == Elem::COARSEN)
1396 : {
1397 6937 : mesh_changed = true;
1398 6937 : break;
1399 42152 : }
1400 :
1401 : // If the mesh changed on any processor, it changed globally
1402 45232 : this->comm().max(mesh_changed);
1403 :
1404 : // And then we may need to widen the ghosting layers.
1405 45232 : if (mesh_changed)
1406 7015 : MeshCommunication().send_coarse_ghosts(_mesh);
1407 :
1408 36713120 : for (auto & elem : _mesh.element_ptr_range())
1409 : {
1410 : // Make sure we transfer the children's boundary id(s)
1411 : // up to its parent when necessary before coarsening.
1412 20156170 : _mesh.get_boundary_info().transfer_boundary_ids_from_children(elem);
1413 :
1414 : // active elements flagged for coarsening will
1415 : // no longer be deleted until MeshRefinement::contract()
1416 20156170 : if (elem->refinement_flag() == Elem::COARSEN)
1417 : {
1418 : // Huh? no level-0 element should be active
1419 : // and flagged for coarsening.
1420 0 : libmesh_assert_not_equal_to (elem->level(), 0);
1421 :
1422 : // Remove this element from any neighbor
1423 : // lists that point to it.
1424 0 : elem->nullify_neighbors();
1425 :
1426 : // Remove any boundary information associated
1427 : // with this element if we do not allow children to have boundary info.
1428 : // Otherwise, we will do the removal in `transfer_boundary_ids_from_children`
1429 : // to make sure we don't delete the information before it is transferred
1430 0 : if (!_mesh.get_boundary_info().is_children_on_boundary_side())
1431 0 : _mesh.get_boundary_info().remove (elem);
1432 :
1433 : // Add this iterator to the _unused_elements
1434 : // data structure so we might fill it.
1435 : // The _unused_elements optimization is currently off.
1436 : // _unused_elements.push_back (it);
1437 :
1438 : // Don't delete the element until
1439 : // MeshRefinement::contract()
1440 : // _mesh.delete_elem(elem);
1441 : }
1442 :
1443 : // inactive elements flagged for coarsening
1444 : // will become active
1445 19234150 : else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1446 : {
1447 396682 : elem->coarsen();
1448 18476 : libmesh_assert (elem->active());
1449 :
1450 : // the mesh has certainly changed
1451 396682 : mesh_changed = true;
1452 : }
1453 20156170 : if (elem->p_refinement_flag() == Elem::COARSEN)
1454 : {
1455 220 : if (elem->p_level() > 0)
1456 : {
1457 14 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1458 206 : elem->set_p_level(elem->p_level() - 1);
1459 206 : mesh_p_changed = true;
1460 : }
1461 : else
1462 : {
1463 0 : elem->set_p_refinement_flag(Elem::DO_NOTHING);
1464 : }
1465 : }
1466 42152 : }
1467 :
1468 45232 : this->comm().max(mesh_p_changed);
1469 :
1470 : // And we may need to update DistributedMesh values reflecting the changes
1471 45232 : if (mesh_changed)
1472 7015 : _mesh.update_parallel_id_counts();
1473 :
1474 : // Node processor ids may need to change if an element of that id
1475 : // was coarsened away
1476 45232 : if (mesh_changed && !_mesh.is_serial())
1477 : {
1478 : // Update the _new_nodes_map so that processors can
1479 : // find requested nodes
1480 1164 : this->update_nodes_map ();
1481 :
1482 1164 : MeshCommunication().make_nodes_parallel_consistent (_mesh);
1483 :
1484 : // Clear the _new_nodes_map
1485 1164 : this->clear();
1486 :
1487 : #ifdef DEBUG
1488 4 : MeshTools::libmesh_assert_valid_procids<Node>(_mesh);
1489 : #endif
1490 : }
1491 :
1492 : // If p levels changed all we need to do is make sure that parent p
1493 : // levels changed in sync
1494 45232 : if (mesh_p_changed && !_mesh.is_serial())
1495 : {
1496 192 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1497 : }
1498 :
1499 48312 : return (mesh_changed || mesh_p_changed);
1500 : }
1501 :
1502 :
1503 :
1504 61681 : bool MeshRefinement::_refine_elements ()
1505 : {
1506 2092 : libmesh_assert(_mesh.is_prepared() || _mesh.is_replicated());
1507 :
1508 : // This function must be run on all processors at once
1509 2092 : parallel_object_only();
1510 :
1511 : // Update the _new_nodes_map so that elements can
1512 : // find nodes to connect to.
1513 61681 : this->update_nodes_map ();
1514 :
1515 4184 : LOG_SCOPE ("_refine_elements()", "MeshRefinement");
1516 :
1517 : // Iterate over the elements, counting the elements
1518 : // flagged for h refinement.
1519 2092 : dof_id_type n_elems_flagged = 0;
1520 :
1521 35927794 : for (auto & elem : _mesh.element_ptr_range())
1522 19713886 : if (elem->refinement_flag() == Elem::REFINE)
1523 1306217 : n_elems_flagged++;
1524 :
1525 : // Construct a local vector of Elem * which have been
1526 : // previously marked for refinement. We reserve enough
1527 : // space to allow for every element to be refined.
1528 2092 : std::vector<Elem *> local_copy_of_elements;
1529 61681 : local_copy_of_elements.reserve(n_elems_flagged);
1530 :
1531 : // If mesh p levels changed, we might need to synchronize parent p
1532 : // levels on a distributed mesh.
1533 61681 : bool mesh_p_changed = false;
1534 :
1535 : // Iterate over the elements, looking for elements flagged for
1536 : // refinement.
1537 :
1538 : // If we are on a ReplicatedMesh, then we just do the refinement in
1539 : // the same order on every processor and everything stays in sync.
1540 :
1541 : // If we are on a DistributedMesh, that's impossible.
1542 : //
1543 : // If the mesh is distributed, we need to make sure that if we end
1544 : // up as the owner of a new node, which might happen if that node is
1545 : // attached to one of our own elements, then we have given it a
1546 : // legitimate node id and our own processor id. We generate
1547 : // legitimate node ids and use our own processor id when we are
1548 : // refining our own elements but not when we refine others'
1549 : // elements. Therefore we want to refine our own elements *first*,
1550 : // thereby generating all nodes which might belong to us, and then
1551 : // refine others' elements *after*, thereby generating nodes with
1552 : // temporary ids which we know we will discard.
1553 : //
1554 : // Even if the DistributedMesh is serialized, we can't just treat it
1555 : // like a ReplicatedMesh, because DistributedMesh doesn't *trust*
1556 : // users to refine partitioned elements in a serialized way, so it
1557 : // assigns temporary ids, so we need to synchronize ids afterward to
1558 : // be safe anyway, so we might as well use the distributed mesh code
1559 : // path.
1560 23622774 : for (auto & elem : _mesh.is_replicated() ? _mesh.active_element_ptr_range() : _mesh.active_local_element_ptr_range())
1561 : {
1562 13054584 : if (elem->refinement_flag() == Elem::REFINE)
1563 836810 : local_copy_of_elements.push_back(elem);
1564 13055890 : if (elem->p_refinement_flag() == Elem::REFINE &&
1565 2602 : elem->active())
1566 : {
1567 11029 : elem->set_p_level(elem->p_level()+1);
1568 9733 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1569 9733 : mesh_p_changed = true;
1570 : }
1571 57497 : }
1572 :
1573 61681 : if (!_mesh.is_replicated())
1574 : {
1575 77168 : for (auto & elem : as_range(_mesh.active_not_local_elements_begin(),
1576 1626150 : _mesh.active_not_local_elements_end()))
1577 : {
1578 736537 : if (elem->refinement_flag() == Elem::REFINE)
1579 411910 : local_copy_of_elements.push_back(elem);
1580 736537 : if (elem->p_refinement_flag() == Elem::REFINE &&
1581 0 : elem->active())
1582 : {
1583 8163 : elem->set_p_level(elem->p_level()+1);
1584 8163 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1585 8163 : mesh_p_changed = true;
1586 : }
1587 38395 : }
1588 : }
1589 :
1590 : // Now iterate over the local copies and refine each one.
1591 : // This may resize the mesh's internal container and invalidate
1592 : // any existing iterators.
1593 1310401 : for (auto & elem : local_copy_of_elements)
1594 1248720 : elem->refine(*this);
1595 :
1596 : // The mesh changed if there were elements h refined
1597 61681 : bool mesh_changed = !local_copy_of_elements.empty();
1598 :
1599 : // If the mesh changed on any processor, it changed globally
1600 61681 : this->comm().max(mesh_changed);
1601 61681 : this->comm().max(mesh_p_changed);
1602 :
1603 : // And we may need to update DistributedMesh values reflecting the changes
1604 61681 : if (mesh_changed)
1605 35637 : _mesh.update_parallel_id_counts();
1606 :
1607 61681 : if (mesh_changed && !_mesh.is_replicated())
1608 : {
1609 22529 : MeshCommunication().make_elems_parallel_consistent (_mesh);
1610 22529 : MeshCommunication().make_new_nodes_parallel_consistent (_mesh);
1611 : #ifdef DEBUG
1612 38 : _mesh.libmesh_assert_valid_parallel_ids();
1613 : #endif
1614 : }
1615 :
1616 : // If we're refining a ReplicatedMesh, then we haven't yet assigned
1617 : // node processor ids. But if we're refining a partitioned
1618 : // ReplicatedMesh, then we *need* to assign node processor ids.
1619 64057 : if (mesh_changed && _mesh.is_replicated() &&
1620 74789 : (_mesh.unpartitioned_elements_begin() ==
1621 86709 : _mesh.unpartitioned_elements_end()))
1622 12717 : Partitioner::set_node_processor_ids(_mesh);
1623 :
1624 61681 : if (mesh_p_changed && !_mesh.is_replicated())
1625 : {
1626 2368 : MeshCommunication().make_p_levels_parallel_consistent (_mesh);
1627 : }
1628 :
1629 : // Clear the _new_nodes_map and _unused_elements data structures.
1630 61681 : this->clear();
1631 :
1632 123362 : return (mesh_changed || mesh_p_changed);
1633 : }
1634 :
1635 :
1636 70632 : void MeshRefinement::_smooth_flags(bool refining, bool coarsening)
1637 : {
1638 : // Smoothing can break in weird ways on a mesh with broken topology
1639 : #ifdef DEBUG
1640 2384 : MeshTools::libmesh_assert_valid_neighbors(_mesh);
1641 : #endif
1642 :
1643 : // Repeat until flag changes match on every processor
1644 0 : do
1645 : {
1646 : // Repeat until coarsening & refinement flags jive
1647 2384 : bool satisfied = false;
1648 196 : do
1649 : {
1650 : // If we're refining or coarsening, hit the corresponding
1651 : // face level test code. Short-circuiting || is our friend
1652 : const bool coarsening_satisfied =
1653 131683 : !coarsening ||
1654 52478 : this->make_coarsening_compatible();
1655 :
1656 : const bool refinement_satisfied =
1657 132651 : !refining ||
1658 53446 : this->make_refinement_compatible();
1659 :
1660 : bool smoothing_satisfied =
1661 79205 : !this->eliminate_unrefined_patches();// &&
1662 :
1663 79205 : if (_edge_level_mismatch_limit)
1664 0 : smoothing_satisfied = smoothing_satisfied &&
1665 0 : !this->limit_level_mismatch_at_edge (_edge_level_mismatch_limit);
1666 :
1667 79205 : if (_node_level_mismatch_limit)
1668 0 : smoothing_satisfied = smoothing_satisfied &&
1669 0 : !this->limit_level_mismatch_at_node (_node_level_mismatch_limit);
1670 :
1671 79205 : if (_overrefined_boundary_limit>=0)
1672 54126 : smoothing_satisfied = smoothing_satisfied &&
1673 26084 : !this->limit_overrefined_boundary(_overrefined_boundary_limit);
1674 :
1675 79205 : if (_underrefined_boundary_limit>=0)
1676 53558 : smoothing_satisfied = smoothing_satisfied &&
1677 25516 : !this->limit_underrefined_boundary(_underrefined_boundary_limit);
1678 :
1679 79205 : satisfied = (coarsening_satisfied &&
1680 81785 : refinement_satisfied &&
1681 : smoothing_satisfied);
1682 :
1683 2580 : libmesh_assert(this->comm().verify(satisfied));
1684 : }
1685 76623 : while (!satisfied);
1686 : }
1687 136918 : while (!_mesh.is_serial() && !this->make_flags_parallel_consistent());
1688 70632 : }
1689 :
1690 :
1691 0 : void MeshRefinement::uniformly_p_refine (unsigned int n)
1692 : {
1693 : // Refine n times
1694 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1695 0 : for (auto & elem : _mesh.active_element_ptr_range())
1696 : {
1697 : // P refine all the active elements
1698 0 : elem->set_p_level(elem->p_level()+1);
1699 0 : elem->set_p_refinement_flag(Elem::JUST_REFINED);
1700 0 : }
1701 0 : }
1702 :
1703 :
1704 :
1705 0 : void MeshRefinement::uniformly_p_coarsen (unsigned int n)
1706 : {
1707 : // Coarsen p times
1708 0 : for (unsigned int rstep=0; rstep<n; rstep++)
1709 0 : for (auto & elem : _mesh.active_element_ptr_range())
1710 0 : if (elem->p_level() > 0)
1711 : {
1712 : // P coarsen all the active elements
1713 0 : elem->set_p_level(elem->p_level()-1);
1714 0 : elem->set_p_refinement_flag(Elem::JUST_COARSENED);
1715 0 : }
1716 0 : }
1717 :
1718 :
1719 :
1720 14558 : void MeshRefinement::uniformly_refine (unsigned int n)
1721 : {
1722 : // Refine n times
1723 : // FIXME - this won't work if n>1 and the mesh
1724 : // has already been attached to an equation system
1725 31295 : for (unsigned int rstep=0; rstep<n; rstep++)
1726 : {
1727 : // Clean up the refinement flags
1728 16737 : this->clean_refinement_flags();
1729 :
1730 : // Flag all the active elements for refinement.
1731 1938158 : for (auto & elem : _mesh.active_element_ptr_range())
1732 1029496 : elem->set_refinement_flag(Elem::REFINE);
1733 :
1734 : // Refine all the elements we just flagged.
1735 16737 : this->_refine_elements();
1736 : }
1737 :
1738 : // Finally, the new mesh probably needs to be prepared for use
1739 14558 : if (n > 0)
1740 13373 : _mesh.prepare_for_use ();
1741 14558 : }
1742 :
1743 :
1744 :
1745 1132 : void MeshRefinement::uniformly_coarsen (unsigned int n)
1746 : {
1747 : // Coarsen n times
1748 2264 : for (unsigned int rstep=0; rstep<n; rstep++)
1749 : {
1750 : // Clean up the refinement flags
1751 1132 : this->clean_refinement_flags();
1752 :
1753 : // Flag all the active elements for coarsening.
1754 471746 : for (auto & elem : _mesh.active_element_ptr_range())
1755 : {
1756 260509 : elem->set_refinement_flag(Elem::COARSEN);
1757 286261 : if (elem->parent())
1758 25752 : elem->parent()->set_refinement_flag(Elem::COARSEN_INACTIVE);
1759 1068 : }
1760 :
1761 : // On a distributed mesh, we may have parent elements with
1762 : // remote active children. To keep flags consistent, we'll need
1763 : // a communication step.
1764 1132 : if (!_mesh.is_replicated())
1765 : {
1766 908 : const processor_id_type n_proc = _mesh.n_processors();
1767 4 : const processor_id_type my_proc_id = _mesh.processor_id();
1768 :
1769 : std::vector<std::vector<dof_id_type>>
1770 916 : parents_to_coarsen(n_proc);
1771 :
1772 117848 : for (const auto & elem : as_range(_mesh.ancestor_elements_begin(), _mesh.ancestor_elements_end()))
1773 57878 : if (elem->processor_id() != my_proc_id &&
1774 48 : elem->refinement_flag() == Elem::COARSEN_INACTIVE)
1775 14280 : parents_to_coarsen[elem->processor_id()].push_back(elem->id());
1776 :
1777 : Parallel::MessageTag
1778 916 : coarsen_tag = this->comm().get_unique_tag();
1779 908 : std::vector<Parallel::Request> coarsen_push_requests(n_proc-1);
1780 :
1781 10452 : for (processor_id_type p = 0; p != n_proc; ++p)
1782 : {
1783 9544 : if (p == my_proc_id)
1784 904 : continue;
1785 :
1786 : Parallel::Request &request =
1787 8636 : coarsen_push_requests[p - (p > my_proc_id)];
1788 :
1789 8636 : _mesh.comm().send
1790 8636 : (p, parents_to_coarsen[p], request, coarsen_tag);
1791 : }
1792 :
1793 9544 : for (processor_id_type p = 1; p != n_proc; ++p)
1794 : {
1795 8 : std::vector<dof_id_type> my_parents_to_coarsen;
1796 8636 : _mesh.comm().receive
1797 8636 : (Parallel::any_source, my_parents_to_coarsen,
1798 8 : coarsen_tag);
1799 :
1800 21980 : for (const auto & id : my_parents_to_coarsen)
1801 : {
1802 13344 : Elem & elem = _mesh.elem_ref(id);
1803 36 : libmesh_assert(elem.refinement_flag() == Elem::INACTIVE ||
1804 : elem.refinement_flag() == Elem::COARSEN_INACTIVE);
1805 36 : elem.set_refinement_flag(Elem::COARSEN_INACTIVE);
1806 : }
1807 : }
1808 :
1809 908 : Parallel::wait(coarsen_push_requests);
1810 :
1811 : SyncRefinementFlags hsync(_mesh, &Elem::refinement_flag,
1812 908 : &Elem::set_refinement_flag);
1813 : sync_dofobject_data_by_id
1814 1812 : (this->comm(), _mesh.not_local_elements_begin(),
1815 916 : _mesh.not_local_elements_end(),
1816 : // We'd like a smaller sync, but this leads to bugs?
1817 : // SyncCoarsenInactive(),
1818 : hsync);
1819 900 : }
1820 :
1821 : // Coarsen all the elements we just flagged.
1822 1132 : this->_coarsen_elements();
1823 : }
1824 :
1825 :
1826 : // Finally, the new mesh probably needs to be prepared for use
1827 1132 : if (n > 0)
1828 1132 : _mesh.prepare_for_use ();
1829 1132 : }
1830 :
1831 :
1832 :
1833 19670953 : Elem * MeshRefinement::topological_neighbor(Elem * elem,
1834 : const PointLocatorBase * point_locator,
1835 : const unsigned int side) const
1836 : {
1837 : #ifdef LIBMESH_ENABLE_PERIODIC
1838 19670953 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1839 : {
1840 488472 : libmesh_assert(point_locator);
1841 786435 : return elem->topological_neighbor(side, _mesh, *point_locator, _periodic_boundaries);
1842 : }
1843 : #endif
1844 19424292 : return elem->neighbor_ptr(side);
1845 : }
1846 :
1847 :
1848 :
1849 1960 : bool MeshRefinement::has_topological_neighbor(const Elem * elem,
1850 : const PointLocatorBase * point_locator,
1851 : const Elem * neighbor) const
1852 : {
1853 : #ifdef LIBMESH_ENABLE_PERIODIC
1854 1960 : if (_periodic_boundaries && !_periodic_boundaries->empty())
1855 : {
1856 0 : libmesh_assert(point_locator);
1857 0 : return elem->has_topological_neighbor(neighbor, _mesh, *point_locator, _periodic_boundaries);
1858 : }
1859 : #endif
1860 1812 : return elem->has_neighbor(neighbor);
1861 : }
1862 :
1863 :
1864 :
1865 : } // namespace libMesh
1866 :
1867 :
1868 : #endif
|
