24#include <deal.II/dofs/dof_accessor.templates.h>
46 template <
int dim,
int spacedim>
57 std::vector<bool> p_flags(
65 template <
int dim,
int spacedim>
68 const std::vector<bool> &p_flags)
80 if (cell->is_locally_owned() && p_flags[cell->active_cell_index()])
82 if (cell->refine_flag_set())
84 const unsigned int super_fe_index =
86 cell->active_fe_index());
89 if (super_fe_index != cell->active_fe_index())
90 cell->set_future_fe_index(super_fe_index);
92 else if (cell->coarsen_flag_set())
94 const unsigned int sub_fe_index =
96 cell->active_fe_index());
99 if (sub_fe_index != cell->active_fe_index())
100 cell->set_future_fe_index(sub_fe_index);
107 template <
int dim,
typename Number,
int spacedim>
112 const Number p_refine_threshold,
113 const Number p_coarsen_threshold,
128 std::vector<bool> p_flags(
132 if (cell->is_locally_owned() &&
133 ((cell->refine_flag_set() &&
134 compare_refine(criteria[cell->active_cell_index()],
135 p_refine_threshold)) ||
136 (cell->coarsen_flag_set() &&
137 compare_coarsen(criteria[cell->active_cell_index()],
138 p_coarsen_threshold))))
139 p_flags[cell->active_cell_index()] =
true;
146 template <
int dim,
typename Number,
int spacedim>
151 const double p_refine_fraction,
152 const double p_coarsen_fraction,
166 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
168 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
173 Number max_criterion_refine = std::numeric_limits<Number>::lowest(),
174 min_criterion_refine = std::numeric_limits<Number>::max();
175 Number max_criterion_coarsen = max_criterion_refine,
176 min_criterion_coarsen = min_criterion_refine;
179 if (cell->is_locally_owned())
181 if (cell->refine_flag_set())
183 max_criterion_refine =
185 criteria(cell->active_cell_index()));
186 min_criterion_refine =
188 criteria(cell->active_cell_index()));
190 else if (cell->coarsen_flag_set())
192 max_criterion_coarsen =
194 criteria(cell->active_cell_index()));
195 min_criterion_coarsen =
197 criteria(cell->active_cell_index()));
204 if (parallel_tria !=
nullptr &&
208 max_criterion_refine =
211 min_criterion_refine =
214 max_criterion_coarsen =
217 min_criterion_coarsen =
224 const Number threshold_refine =
225 min_criterion_refine +
227 (max_criterion_refine - min_criterion_refine),
229 min_criterion_coarsen +
231 (max_criterion_coarsen - min_criterion_coarsen);
243 template <
int dim,
typename Number,
int spacedim>
248 const double p_refine_fraction,
249 const double p_coarsen_fraction,
263 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
265 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
273 [](
const Number &,
const Number &) {
return false; };
275 [](
const Number &,
const Number &) {
return true; };
279 unsigned int n_flags_refinement = 0;
280 unsigned int n_flags_coarsening = 0;
285 for (
const auto &cell :
287 if (!cell->is_artificial() && cell->is_locally_owned())
289 if (cell->refine_flag_set())
290 indicators_refinement(n_flags_refinement++) =
291 criteria(cell->active_cell_index());
292 else if (cell->coarsen_flag_set())
293 indicators_coarsening(n_flags_coarsening++) =
294 criteria(cell->active_cell_index());
315 Number threshold_refinement = 0.;
316 Number threshold_coarsening = 0.;
317 auto reference_compare_refine = std::cref(compare_refine);
318 auto reference_compare_coarsen = std::cref(compare_coarsen);
323 if (parallel_tria !=
nullptr &&
327#ifndef DEAL_II_WITH_P4EST
338 const unsigned int n_global_flags_refinement =
340 const unsigned int n_global_flags_coarsening =
343 const unsigned int target_index_refinement =
344 static_cast<unsigned int>(
345 std::floor(p_refine_fraction * n_global_flags_refinement));
346 const unsigned int target_index_coarsening =
347 static_cast<unsigned int>(
348 std::ceil((1 - p_coarsen_fraction) * n_global_flags_coarsening));
352 const std::pair<Number, Number> global_min_max_refinement =
357 const std::pair<Number, Number> global_min_max_coarsening =
363 if (target_index_refinement == 0)
364 reference_compare_refine = std::cref(compare_false);
365 else if (target_index_refinement == n_global_flags_refinement)
366 reference_compare_refine = std::cref(compare_true);
370 indicators_refinement,
371 global_min_max_refinement,
372 target_index_refinement,
375 if (target_index_coarsening == n_global_flags_coarsening)
376 reference_compare_coarsen = std::cref(compare_false);
377 else if (target_index_coarsening == 0)
378 reference_compare_coarsen = std::cref(compare_true);
382 indicators_coarsening,
383 global_min_max_coarsening,
384 target_index_coarsening,
395 const unsigned int n_p_refine_cells =
static_cast<unsigned int>(
396 std::floor(p_refine_fraction * n_flags_refinement));
397 const unsigned int n_p_coarsen_cells =
static_cast<unsigned int>(
398 std::floor(p_coarsen_fraction * n_flags_coarsening));
401 if (n_p_refine_cells == 0)
402 reference_compare_refine = std::cref(compare_false);
403 else if (n_p_refine_cells == n_flags_refinement)
404 reference_compare_refine = std::cref(compare_true);
407 std::nth_element(indicators_refinement.
begin(),
408 indicators_refinement.
begin() +
409 n_p_refine_cells - 1,
410 indicators_refinement.
end(),
411 std::greater<Number>());
412 threshold_refinement =
413 *(indicators_refinement.
begin() + n_p_refine_cells - 1);
416 if (n_p_coarsen_cells == 0)
417 reference_compare_coarsen = std::cref(compare_false);
418 else if (n_p_coarsen_cells == n_flags_coarsening)
419 reference_compare_coarsen = std::cref(compare_true);
422 std::nth_element(indicators_coarsening.
begin(),
423 indicators_coarsening.
begin() +
424 n_p_coarsen_cells - 1,
425 indicators_coarsening.
end(),
426 std::less<Number>());
427 threshold_coarsening =
428 *(indicators_coarsening.
begin() + n_p_coarsen_cells - 1);
435 threshold_refinement,
436 threshold_coarsening,
437 std::cref(reference_compare_refine),
439 reference_compare_coarsen));
444 template <
int dim,
typename Number,
int spacedim>
456 sobolev_indices.
size());
459 if (cell->is_locally_owned())
461 if (cell->refine_flag_set())
463 const unsigned int super_fe_index =
465 cell->active_fe_index());
468 if (super_fe_index != cell->active_fe_index())
470 const unsigned int super_fe_degree =
473 if (sobolev_indices[cell->active_cell_index()] >
475 cell->set_future_fe_index(super_fe_index);
478 else if (cell->coarsen_flag_set())
480 const unsigned int sub_fe_index =
482 cell->active_fe_index());
485 if (sub_fe_index != cell->active_fe_index())
487 const unsigned int sub_fe_degree =
490 if (sobolev_indices[cell->active_cell_index()] <
492 cell->set_future_fe_index(sub_fe_index);
500 template <
int dim,
typename Number,
int spacedim>
522 std::vector<bool> p_flags(
526 if (cell->is_locally_owned() &&
527 ((cell->refine_flag_set() &&
528 compare_refine(criteria[cell->active_cell_index()],
529 references[cell->active_cell_index()])) ||
530 (cell->coarsen_flag_set() &&
531 compare_coarsen(criteria[cell->active_cell_index()],
532 references[cell->active_cell_index()]))))
533 p_flags[cell->active_cell_index()] =
true;
543 template <
int dim,
typename Number,
int spacedim>
548 const double gamma_p,
549 const double gamma_h,
550 const double gamma_n)
557 error_indicators.
size());
559 predicted_errors.
size());
560 Assert(0 < gamma_p && gamma_p < 1,
570 std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
unsigned int>
571 future_fe_indices_on_coarsened_cells;
574 predicted_errors = error_indicators;
580 if (!(cell->future_fe_index_set()) && !(cell->refine_flag_set()) &&
581 !(cell->coarsen_flag_set()))
583 predicted_errors[cell->active_cell_index()] *= gamma_n;
589 if (cell->coarsen_flag_set())
592 ExcMessage(
"A coarse cell is flagged for coarsening. "
593 "Please read the note in the documentation "
594 "of predict_error()."));
598 const auto &parent = cell->parent();
599 if (future_fe_indices_on_coarsened_cells.find(parent) ==
600 future_fe_indices_on_coarsened_cells.end())
603 for (
const auto &child : parent->child_iterators())
604 Assert(child->is_active() && child->coarsen_flag_set(),
606 dim>::ExcInconsistentCoarseningFlags());
609 parent_future_fe_index =
610 internal::hp::DoFHandlerImplementation::
611 dominated_future_fe_on_children<dim, spacedim>(parent);
613 future_fe_indices_on_coarsened_cells.insert(
614 {parent, parent_future_fe_index});
618 parent_future_fe_index =
619 future_fe_indices_on_coarsened_cells[parent];
633 if (cell->future_fe_index_set())
635 if (future_fe_degree > cell->get_fe().degree)
636 predicted_errors[cell->active_cell_index()] *=
638 future_fe_degree - cell->get_fe().degree);
639 else if (future_fe_degree < cell->get_fe().degree)
640 predicted_errors[cell->active_cell_index()] /=
642 cell->get_fe().degree - future_fe_degree);
651 if (cell->refine_flag_set())
653 predicted_errors[cell->active_cell_index()] *=
659 else if (cell->coarsen_flag_set())
661 predicted_errors[cell->active_cell_index()] /=
675 template <
int dim,
int spacedim>
687 if (cell->is_locally_owned() && cell->future_fe_index_set())
689 cell->clear_refine_flag();
690 cell->clear_coarsen_flag();
696 template <
int dim,
int spacedim>
716 if (cell->is_locally_owned() && cell->future_fe_index_set())
718 cell->clear_refine_flag();
724 if (cell->coarsen_flag_set())
726 if (cell->level() == 0)
731 cell->clear_coarsen_flag();
735 const auto &parent = cell->parent();
736 const unsigned int n_children = parent->n_children();
738 unsigned int h_flagged_children = 0, p_flagged_children = 0;
739 for (
const auto &child : parent->child_iterators())
741 if (child->is_active())
743 if (child->is_locally_owned())
745 if (child->coarsen_flag_set())
746 ++h_flagged_children;
747 if (child->future_fe_index_set())
748 ++p_flagged_children;
750 else if (child->is_ghost())
756 (
dynamic_cast<const parallel::shared::
757 Triangulation<dim, spacedim> *
>(
761 if (child->coarsen_flag_set())
762 ++h_flagged_children;
768 if (internal::DoFCellAccessorImplementation::
770 future_fe_index_set<dim, spacedim, false>(
772 ++p_flagged_children;
783 if (h_flagged_children == n_children &&
784 p_flagged_children != n_children)
788 for (
const auto &child : parent->child_iterators())
793 if (child->is_locally_owned())
794 child->clear_future_fe_index();
801 for (
const auto &child : parent->child_iterators())
803 if (child->is_active() && child->is_locally_owned())
804 child->clear_coarsen_flag();
816 template <
int dim,
int spacedim>
819 const unsigned int max_difference,
820 const unsigned int contains_fe_index)
831 "This function does not serve any purpose for max_difference = 0."));
843 const auto invalid_level =
static_cast<level_type
>(-1);
847 const std::vector<unsigned int> fe_index_for_hierarchy_level =
854 std::vector<unsigned int> hierarchy_level_for_fe_index(
856 for (
unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l)
857 hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l;
870 if (
const auto parallel_tria =
875 parallel_tria->global_active_cell_index_partitioner().lock());
885 future_levels[cell->global_active_cell_index()] =
886 hierarchy_level_for_fe_index[cell->future_fe_index()];
901 const auto update_neighbor_level =
902 [&future_levels, max_difference, invalid_level](
903 const auto &neighbor,
const level_type cell_level) ->
bool {
908 if (neighbor->is_locally_owned())
910 const level_type neighbor_level =
static_cast<level_type
>(
911 future_levels[neighbor->global_active_cell_index()]);
914 if (neighbor_level == invalid_level)
917 if ((cell_level - max_difference) > neighbor_level)
919 future_levels[neighbor->global_active_cell_index()] =
920 cell_level - max_difference;
944 const auto prepare_level_for_parent = [&](
const auto &neighbor) ->
bool {
946 if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
948 const auto parent = neighbor->parent();
950 std::vector<unsigned int> future_levels_children;
951 future_levels_children.reserve(parent->n_children());
952 for (
const auto &child : parent->child_iterators())
954 Assert(child->is_active() && child->coarsen_flag_set(),
956 ExcInconsistentCoarseningFlags()));
958 const level_type child_level =
static_cast<level_type
>(
959 future_levels[child->global_active_cell_index()]);
960 Assert(child_level != invalid_level,
962 "The FiniteElement on one of the siblings of "
963 "a cell you are trying to coarsen is not part "
964 "of the registered p-adaptation hierarchy."));
965 future_levels_children.push_back(child_level);
969 const unsigned int max_level_children =
970 *std::max_element(future_levels_children.begin(),
971 future_levels_children.end());
973 bool children_changed =
false;
974 for (
const auto &child : parent->child_iterators())
978 if (child->is_locally_owned() &&
979 future_levels[child->global_active_cell_index()] !=
982 future_levels[child->global_active_cell_index()] =
984 children_changed =
true;
986 return children_changed;
992 bool levels_changed =
false;
993 bool levels_changed_in_cycle;
996 levels_changed_in_cycle =
false;
1001 if (!cell->is_artificial())
1003 const level_type cell_level =
static_cast<level_type
>(
1004 future_levels[cell->global_active_cell_index()]);
1007 if (cell_level == invalid_level)
1012 if (cell_level <= max_difference)
1015 for (
unsigned int f = 0; f < cell->n_faces(); ++f)
1016 if (cell->face(f)->at_boundary() ==
false)
1018 if (cell->face(f)->has_children())
1020 for (
unsigned int sf = 0;
1021 sf < cell->face(f)->n_children();
1024 const auto neighbor =
1025 cell->neighbor_child_on_subface(f, sf);
1027 levels_changed_in_cycle |=
1028 update_neighbor_level(neighbor, cell_level);
1030 levels_changed_in_cycle |=
1031 prepare_level_for_parent(neighbor);
1036 const auto neighbor = cell->neighbor(f);
1038 levels_changed_in_cycle |=
1039 update_neighbor_level(neighbor, cell_level);
1041 levels_changed_in_cycle |=
1042 prepare_level_for_parent(neighbor);
1047 levels_changed_in_cycle =
1050 levels_changed |= levels_changed_in_cycle;
1052 while (levels_changed_in_cycle);
1058 const level_type cell_level =
static_cast<level_type
>(
1059 future_levels[cell->global_active_cell_index()]);
1061 if (cell_level != invalid_level)
1063 const unsigned int fe_index =
1064 fe_index_for_hierarchy_level[cell_level];
1066 if (fe_index != cell->active_fe_index())
1067 cell->set_future_fe_index(fe_index);
1069 cell->clear_future_fe_index();
1073 return levels_changed;
1080#include "refinement.inst"
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
MPI_Comm get_communicator() const
void update_ghost_values() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
unsigned int n_active_cells() const
virtual size_type size() const override
void grow_or_shrink(const size_type N)
unsigned int size() const
unsigned int previous_in_hierarchy(const unsigned int fe_index) const
std::vector< unsigned int > get_hierarchy_sequence(const unsigned int fe_index=0) const
unsigned int next_in_hierarchy(const unsigned int fe_index) const
virtual MPI_Comm get_communicator() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInvalidParameterValue()
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_ASSERT_UNREACHABLE()
T sum(const T &t, const MPI_Comm mpi_communicator)
T logical_or(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
void force_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_reference(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const Vector< Number > &references, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen)
void full_p_adaptivity(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_regularity(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &sobolev_indices)
bool limit_p_level_difference(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int max_difference=1, const unsigned int contains_fe_index=0)
void predict_error(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &error_indicators, Vector< Number > &predicted_errors, const double gamma_p=std::sqrt(0.4), const double gamma_h=2., const double gamma_n=1.)
void p_adaptivity_from_absolute_threshold(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const Number p_refine_threshold, const Number p_coarsen_threshold, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
void choose_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_fixed_number(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
std::function< bool(const Number &, const Number &)> ComparisonFunction
void p_adaptivity_from_relative_threshold(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
void p_adaptivity_from_flags(const DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const types::global_cell_index n_target_cells, const MPI_Comm mpi_communicator)
std::pair< number, number > compute_global_min_and_max_at_root(const ::Vector< number > &criteria, const MPI_Comm mpi_communicator)
static const unsigned int invalid_unsigned_int
typename type_identity< T >::type type_identity_t
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index