25 #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,
127 std::vector<bool> p_flags(
131 if (cell->is_locally_owned() &&
132 ((cell->refine_flag_set() &&
133 compare_refine(criteria[cell->active_cell_index()],
134 p_refine_threshold)) ||
135 (cell->coarsen_flag_set() &&
136 compare_coarsen(criteria[cell->active_cell_index()],
137 p_coarsen_threshold))))
138 p_flags[cell->active_cell_index()] =
true;
145 template <
int dim,
typename Number,
int spacedim>
150 const double p_refine_fraction,
151 const double p_coarsen_fraction,
164 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
166 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
171 Number max_criterion_refine = std::numeric_limits<Number>::lowest(),
173 Number max_criterion_coarsen = max_criterion_refine,
174 min_criterion_coarsen = min_criterion_refine;
177 if (cell->is_locally_owned())
179 if (cell->refine_flag_set())
181 max_criterion_refine =
183 criteria(cell->active_cell_index()));
184 min_criterion_refine =
186 criteria(cell->active_cell_index()));
188 else if (cell->coarsen_flag_set())
190 max_criterion_coarsen =
192 criteria(cell->active_cell_index()));
193 min_criterion_coarsen =
195 criteria(cell->active_cell_index()));
202 if (parallel_tria !=
nullptr &&
206 max_criterion_refine =
209 min_criterion_refine =
212 max_criterion_coarsen =
215 min_criterion_coarsen =
222 const Number threshold_refine =
223 min_criterion_refine +
225 (max_criterion_refine - min_criterion_refine),
227 min_criterion_coarsen +
229 (max_criterion_coarsen - min_criterion_coarsen);
241 template <
int dim,
typename Number,
int spacedim>
246 const double p_refine_fraction,
247 const double p_coarsen_fraction,
260 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
262 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
270 [](
const Number &,
const Number &) {
return false; };
272 [](
const Number &,
const Number &) {
return true; };
276 unsigned int n_flags_refinement = 0;
277 unsigned int n_flags_coarsening = 0;
282 for (
const auto &cell :
284 if (!cell->is_artificial() && cell->is_locally_owned())
286 if (cell->refine_flag_set())
287 indicators_refinement(n_flags_refinement++) =
288 criteria(cell->active_cell_index());
289 else if (cell->coarsen_flag_set())
290 indicators_coarsening(n_flags_coarsening++) =
291 criteria(cell->active_cell_index());
312 Number threshold_refinement = 0.;
313 Number threshold_coarsening = 0.;
314 auto reference_compare_refine = std::cref(compare_refine);
315 auto reference_compare_coarsen = std::cref(compare_coarsen);
320 if (parallel_tria !=
nullptr &&
324 #ifndef DEAL_II_WITH_P4EST
335 const unsigned int n_global_flags_refinement =
337 const unsigned int n_global_flags_coarsening =
340 const unsigned int target_index_refinement =
341 static_cast<unsigned int>(
342 std::floor(p_refine_fraction * n_global_flags_refinement));
343 const unsigned int target_index_coarsening =
344 static_cast<unsigned int>(
345 std::ceil((1 - p_coarsen_fraction) * n_global_flags_coarsening));
349 const std::pair<Number, Number> global_min_max_refinement =
354 const std::pair<Number, Number> global_min_max_coarsening =
360 if (target_index_refinement == 0)
361 reference_compare_refine = std::cref(compare_false);
362 else if (target_index_refinement == n_global_flags_refinement)
363 reference_compare_refine = std::cref(compare_true);
367 indicators_refinement,
368 global_min_max_refinement,
369 target_index_refinement,
372 if (target_index_coarsening == n_global_flags_coarsening)
373 reference_compare_coarsen = std::cref(compare_false);
374 else if (target_index_coarsening == 0)
375 reference_compare_coarsen = std::cref(compare_true);
379 indicators_coarsening,
380 global_min_max_coarsening,
381 target_index_coarsening,
392 const unsigned int n_p_refine_cells =
static_cast<unsigned int>(
393 std::floor(p_refine_fraction * n_flags_refinement));
394 const unsigned int n_p_coarsen_cells =
static_cast<unsigned int>(
395 std::floor(p_coarsen_fraction * n_flags_coarsening));
398 if (n_p_refine_cells == 0)
399 reference_compare_refine = std::cref(compare_false);
400 else if (n_p_refine_cells == n_flags_refinement)
401 reference_compare_refine = std::cref(compare_true);
404 std::nth_element(indicators_refinement.
begin(),
405 indicators_refinement.
begin() +
406 n_p_refine_cells - 1,
407 indicators_refinement.
end(),
408 std::greater<Number>());
409 threshold_refinement =
410 *(indicators_refinement.
begin() + n_p_refine_cells - 1);
413 if (n_p_coarsen_cells == 0)
414 reference_compare_coarsen = std::cref(compare_false);
415 else if (n_p_coarsen_cells == n_flags_coarsening)
416 reference_compare_coarsen = std::cref(compare_true);
419 std::nth_element(indicators_coarsening.
begin(),
420 indicators_coarsening.
begin() +
421 n_p_coarsen_cells - 1,
422 indicators_coarsening.
end(),
423 std::less<Number>());
424 threshold_coarsening =
425 *(indicators_coarsening.
begin() + n_p_coarsen_cells - 1);
432 threshold_refinement,
433 threshold_coarsening,
434 std::cref(reference_compare_refine),
436 reference_compare_coarsen));
441 template <
int dim,
typename Number,
int spacedim>
453 sobolev_indices.
size());
456 if (cell->is_locally_owned())
458 if (cell->refine_flag_set())
460 const unsigned int super_fe_index =
462 cell->active_fe_index());
465 if (super_fe_index != cell->active_fe_index())
467 const unsigned int super_fe_degree =
470 if (sobolev_indices[cell->active_cell_index()] >
472 cell->set_future_fe_index(super_fe_index);
475 else if (cell->coarsen_flag_set())
477 const unsigned int sub_fe_index =
479 cell->active_fe_index());
482 if (sub_fe_index != cell->active_fe_index())
484 const unsigned int sub_fe_degree =
487 if (sobolev_indices[cell->active_cell_index()] <
489 cell->set_future_fe_index(sub_fe_index);
497 template <
int dim,
typename Number,
int spacedim>
518 std::vector<bool> p_flags(
522 if (cell->is_locally_owned() &&
523 ((cell->refine_flag_set() &&
524 compare_refine(criteria[cell->active_cell_index()],
525 references[cell->active_cell_index()])) ||
526 (cell->coarsen_flag_set() &&
527 compare_coarsen(criteria[cell->active_cell_index()],
528 references[cell->active_cell_index()]))))
529 p_flags[cell->active_cell_index()] =
true;
539 template <
int dim,
typename Number,
int spacedim>
544 const double gamma_p,
545 const double gamma_h,
546 const double gamma_n)
553 error_indicators.
size());
555 predicted_errors.
size());
556 Assert(0 < gamma_p && gamma_p < 1,
566 std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
unsigned int>
567 future_fe_indices_on_coarsened_cells;
570 predicted_errors = error_indicators;
576 if (!(cell->future_fe_index_set()) && !(cell->refine_flag_set()) &&
577 !(cell->coarsen_flag_set()))
579 predicted_errors[cell->active_cell_index()] *= gamma_n;
585 if (cell->coarsen_flag_set())
589 const auto &parent = cell->parent();
590 if (future_fe_indices_on_coarsened_cells.find(parent) ==
591 future_fe_indices_on_coarsened_cells.end())
594 for (
const auto &child : parent->child_iterators())
595 Assert(child->is_active() && child->coarsen_flag_set(),
597 dim>::ExcInconsistentCoarseningFlags());
600 parent_future_fe_index =
601 internal::hp::DoFHandlerImplementation::
602 dominated_future_fe_on_children<dim, spacedim>(parent);
604 future_fe_indices_on_coarsened_cells.insert(
605 {parent, parent_future_fe_index});
609 parent_future_fe_index =
610 future_fe_indices_on_coarsened_cells[parent];
624 if (cell->future_fe_index_set())
626 predicted_errors[cell->active_cell_index()] *=
628 int(future_fe_degree) -
int(cell->get_fe().degree));
632 if (cell->refine_flag_set())
634 predicted_errors[cell->active_cell_index()] *=
635 (gamma_h * std::pow(.5, future_fe_degree));
640 else if (cell->coarsen_flag_set())
642 predicted_errors[cell->active_cell_index()] /=
643 (gamma_h * std::pow(.5, future_fe_degree));
656 template <
int dim,
int spacedim>
668 if (cell->is_locally_owned() && cell->future_fe_index_set())
670 cell->clear_refine_flag();
671 cell->clear_coarsen_flag();
677 template <
int dim,
int spacedim>
697 if (cell->is_locally_owned() && cell->future_fe_index_set())
699 cell->clear_refine_flag();
705 if (cell->coarsen_flag_set())
707 const auto & parent = cell->parent();
708 const unsigned int n_children = parent->n_children();
710 unsigned int h_flagged_children = 0, p_flagged_children = 0;
711 for (
const auto &child : parent->child_iterators())
713 if (child->is_active())
715 if (child->is_locally_owned())
717 if (child->coarsen_flag_set())
718 ++h_flagged_children;
719 if (child->future_fe_index_set())
720 ++p_flagged_children;
722 else if (child->is_ghost())
728 (
dynamic_cast<const parallel::shared::
729 Triangulation<dim, spacedim> *
>(
733 if (child->coarsen_flag_set())
734 ++h_flagged_children;
740 if (internal::DoFCellAccessorImplementation::
742 future_fe_index_set<dim, spacedim, false>(
744 ++p_flagged_children;
755 if (h_flagged_children == n_children &&
756 p_flagged_children != n_children)
760 for (
const auto &child : parent->child_iterators())
765 if (child->is_locally_owned())
766 child->clear_future_fe_index();
773 for (
const auto &child : parent->child_iterators())
775 if (child->is_active() && child->is_locally_owned())
776 child->clear_coarsen_flag();
788 template <
int dim,
int spacedim>
791 const unsigned int max_difference,
792 const unsigned int contains_fe_index)
803 "This function does not serve any purpose for max_difference = 0."));
816 const auto invalid_level =
static_cast<level_type
>(-1);
820 const std::vector<unsigned int> fe_index_for_hierarchy_level =
827 std::vector<unsigned int> hierarchy_level_for_fe_index(
829 for (
unsigned int l = 0;
l < fe_index_for_hierarchy_level.size(); ++
l)
830 hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[
l]] =
l;
843 if (
const auto parallel_tria =
848 parallel_tria->global_active_cell_index_partitioner().lock());
858 future_levels[cell->global_active_cell_index()] =
859 hierarchy_level_for_fe_index[cell->future_fe_index()];
874 const auto update_neighbor_level =
875 [&future_levels, max_difference, invalid_level](
876 const auto &neighbor,
const level_type cell_level) ->
bool {
881 if (neighbor->is_locally_owned())
883 const level_type neighbor_level =
static_cast<level_type
>(
884 future_levels[neighbor->global_active_cell_index()]);
887 if (neighbor_level == invalid_level)
890 if ((cell_level - max_difference) > neighbor_level)
892 future_levels[neighbor->global_active_cell_index()] =
893 cell_level - max_difference;
917 const auto prepare_level_for_parent = [&](
const auto &neighbor) ->
bool {
919 if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
921 const auto parent = neighbor->parent();
923 std::vector<unsigned int> future_levels_children;
924 future_levels_children.reserve(parent->n_children());
925 for (
const auto &child : parent->child_iterators())
927 Assert(child->is_active() && child->coarsen_flag_set(),
929 ExcInconsistentCoarseningFlags()));
931 const level_type child_level =
static_cast<level_type
>(
932 future_levels[child->global_active_cell_index()]);
933 Assert(child_level != invalid_level,
935 "The FiniteElement on one of the siblings of "
936 "a cell you are trying to coarsen is not part "
937 "of the registered p-adaptation hierarchy."));
938 future_levels_children.push_back(child_level);
942 const unsigned int max_level_children =
943 *std::max_element(future_levels_children.begin(),
944 future_levels_children.end());
946 bool children_changed =
false;
947 for (
const auto &child : parent->child_iterators())
951 if (child->is_locally_owned() &&
952 future_levels[child->global_active_cell_index()] !=
955 future_levels[child->global_active_cell_index()] =
957 children_changed =
true;
959 return children_changed;
965 bool levels_changed =
false;
966 bool levels_changed_in_cycle;
969 levels_changed_in_cycle =
false;
974 if (!cell->is_artificial())
976 const level_type cell_level =
static_cast<level_type
>(
977 future_levels[cell->global_active_cell_index()]);
980 if (cell_level == invalid_level)
985 if (cell_level <= max_difference)
988 for (
unsigned int f = 0; f < cell->n_faces(); ++f)
989 if (cell->face(f)->at_boundary() ==
false)
991 if (cell->face(f)->has_children())
993 for (
unsigned int sf = 0;
994 sf < cell->face(f)->n_children();
997 const auto neighbor =
998 cell->neighbor_child_on_subface(f, sf);
1000 levels_changed_in_cycle |=
1001 update_neighbor_level(neighbor, cell_level);
1003 levels_changed_in_cycle |=
1004 prepare_level_for_parent(neighbor);
1009 const auto neighbor = cell->neighbor(f);
1011 levels_changed_in_cycle |=
1012 update_neighbor_level(neighbor, cell_level);
1014 levels_changed_in_cycle |=
1015 prepare_level_for_parent(neighbor);
1020 levels_changed_in_cycle =
1023 levels_changed |= levels_changed_in_cycle;
1025 while (levels_changed_in_cycle);
1031 const level_type cell_level =
static_cast<level_type
>(
1032 future_levels[cell->global_active_cell_index()]);
1034 if (cell_level != invalid_level)
1036 const unsigned int fe_index =
1037 fe_index_for_hierarchy_level[cell_level];
1040 if (fe_index != cell->active_fe_index())
1041 cell->set_future_fe_index(fe_index);
1045 return levels_changed;
1052 #include "refinement.inst"
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
MPI_Comm get_communicator() const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void update_ghost_values() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
unsigned int n_active_cells() const
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
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidParameterValue()
Expression ceil(const Expression &x)
Expression floor(const Expression &x)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T logical_or(const T &t, const MPI_Comm &mpi_communicator)
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
T max(const T &t, const MPI_Comm &mpi_communicator)
void force_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
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_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< typename identity< Number >::type > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< typename identity< Number >::type > &compare_coarsen=std::less_equal< Number >())
void choose_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
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< typename identity< Number >::type > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< typename identity< Number >::type > &compare_coarsen=std::less_equal< Number >())
void p_adaptivity_from_reference(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const Vector< Number > &references, const ComparisonFunction< typename identity< Number >::type > &compare_refine, const ComparisonFunction< typename identity< Number >::type > &compare_coarsen)
std::function< bool(const Number &, const Number &)> ComparisonFunction
void p_adaptivity_from_flags(const DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
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< typename identity< Number >::type > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< typename identity< Number >::type > &compare_coarsen=std::less_equal< Number >())
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