From 443e9115c527e99d4f14e74b95178528c4f46174 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Sun, 8 Sep 2019 22:56:44 +0200 Subject: [PATCH] Major redesign of hp::Refinement. Generalizing function and parameter names. Introducing comparators. Distinguishing between absolute and relative thresholds. Elaborating on documentation. Bugfix: Clear coarsen flags only on active children in hp::Refinement::choose_p_over_h(). --- include/deal.II/hp/refinement.h | 296 +++++++++++++++++++---------- source/hp/refinement.cc | 181 ++++++++++-------- source/hp/refinement.inst.in | 31 ++- tests/hp/error_prediction.cc | 8 + tests/hp/p_adaptivity_flags.cc | 66 +++++-- tests/hp/p_adaptivity_flags.output | 36 ++-- tests/mpi/petsc_step-27.cc | 8 +- tests/mpi/trilinos_step-27.cc | 8 +- 8 files changed, 406 insertions(+), 228 deletions(-) diff --git a/include/deal.II/hp/refinement.h b/include/deal.II/hp/refinement.h index 70d9bc1128..c87db67890 100644 --- a/include/deal.II/hp/refinement.h +++ b/include/deal.II/hp/refinement.h @@ -19,10 +19,14 @@ #include -#include +#include + +#include + DEAL_II_NAMESPACE_OPEN + // forward declarations #ifndef DOXYGEN template @@ -35,6 +39,7 @@ namespace hp } #endif + namespace hp { /** @@ -95,21 +100,21 @@ namespace hp * @code * // step 1: flag cells for refinement or coarsening * Vector estimated_error_per_cell (triangulation.n_active_cells()); - * KellyErrorEstimator::estimate (hp_dof_handler, - * QGauss (quadrature_points), - * typename FunctionMap::type(), - * solution, - * estimated_error_per_cell); - * GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, - * estimated_error_per_cell, - * top_fraction, - * bottom_fraction); + * KellyErrorEstimator::estimate(hp_dof_handler, + * QGauss (quadrature_points), + * typename FunctionMap::type(), + * solution, + * estimated_error_per_cell); + * GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, + * estimated_error_per_cell, + * top_fraction, + * bottom_fraction); * * // step 2: set future finite element indices on flagged cells - * hp::Refinement::full_p_adaptivity (hp_dof_handler); + * hp::Refinement::full_p_adaptivity(hp_dof_handler); * * // step 3: decide whether h- or p-adaptive methods will be supplied - * hp::Refinement::force_p_over_h (hp_dof_handler); + * hp::Refinement::force_p_over_h(hp_dof_handler); * * // step 4: prepare solutions to be transferred * ... @@ -122,6 +127,20 @@ namespace hp */ namespace Refinement { + /** + * An alias that defines the characteristics of a function that can be used + * as a comparison criterion for deciding whether to perform h- or + * p-adaptation. + * + * Such functions take two numbers as arguments: The first one corresponds + * to the provided criterion, while the other one conforms to the reference. + * The result of the comparision will be returned as a boolean. + */ + template + using ComparisonFunction = + std::function::type &, + const typename identity::type &)>; + /** * @name Setting p-adaptivity flags * @{ @@ -140,10 +159,10 @@ namespace hp full_p_adaptivity(const hp::DoFHandler &dof_handler); /** - * Adapt the finite element on cells that have been specifically flagged for - * p-adaptation via the parameter @p p_flags. Future finite elements will - * only be assigned if cells have been flagged for refinement and coarsening - * beforehand. + * Adapt which finite element to use on cells that have been specifically + * flagged for p-adaptation via the parameter @p p_flags. Future finite + * elements will only be assigned if cells have been flagged for refinement + * and coarsening beforehand. * * Each entry of the parameter @p p_flags needs to correspond to an active * cell. @@ -158,25 +177,24 @@ namespace hp const std::vector & p_flags); /** - * Adapt the finite element on cells whose smoothness indicators meet a - * certain threshold. - * - * The threshold will be chosen for refined and coarsened cells - * individually. For each class of cells, we determine the maximal and - * minimal values of the smoothness indicators and determine the threshold - * by linear interpolation between these limits. Parameters - * @p p_refine_fraction and @p p_refine_coarsen are used as interpolation - * factors, where `0` corresponds to the minimal and `1` to the maximal - * value. By default, mean values are considered as thresholds. - * - * We consider a cell for p-refinement if it is flagged for refinement and - * its smoothness indicator is larger than the corresponding threshold. The - * same applies for p-coarsening, but the cell's indicator must be lower - * than the threshold. - * - * Each entry of the parameter @p smoothness_indicators needs to correspond - * to an active cell. Parameters @p p_refine_fraction and - * @p p_coarsen_fraction need to be in the interval $[0,1]$. + * Adapt which finite element to use on cells whose criteria meet a certain + * absolute threshold. + * + * For p-refinement and p-coarsening, two separate thresholds need to + * provided via parameters @p p_refine_threshold and @p p_coarsen_threshold. + * + * We consider a cell for p-adaptivity if it is currently flagged for + * refinement or coarsening and its criterion successfully compares to the + * corresponding threshold. Let us be more specific on the default case: We + * consider a cell for p-refinement if it is flagged for refinement and its + * criterion is larger than the corresponding threshold. The same applies + * for p-coarsening, but the cell's criterion must be lower than the + * threshold. However, different compare function objects can be supplied + * via the parameters @p compare_refine and @p compare_coarsen to impose + * different decision strategies. + * + * Each entry of the parameter @p criteria needs to correspond to an active + * cell. * * @note Preceeding calls of Triangulation::prepare_for_coarsening_and_refinement() * may change refine and coarsen flags, which will ultimately change the @@ -184,14 +202,56 @@ namespace hp */ template void - p_adaptivity_from_threshold( + p_adaptivity_from_absolute_threshold( const hp::DoFHandler &dof_handler, - const Vector & smoothness_indicators, + const Vector & criteria, + const Number p_refine_threshold, + const Number p_coarsen_threshold, + const ComparisonFunction &compare_refine = std::greater(), + const ComparisonFunction &compare_coarsen = std::less()); + + /** + * Adapt which finite element to use on cells whose criteria meet a certain + * threshold relative to the overall range of criterion values. + * + * The threshold will be determined for refined and coarsened cells + * separately based on the currently set refinement markers. For each class + * of cells, we determine the maximal and minimal values of all criteria and + * determine the threshold by linear interpolation between these limits. + * Parameters @p p_refine_fraction and @p p_refine_coarsen are used as + * interpolation factors, where `0` corresponds to the minimal and `1` to + * the maximal value. By default, mean values are considered as thresholds. + * + * We consider a cell for p-adaptivity if it is currently flagged for + * refinement or coarsening and its criterion successfully compares to the + * corresponding threshold. Let us be more specific on the default case: We + * consider a cell for p-refinement if it is flagged for refinement and its + * criterion is larger than the corresponding threshold. The same applies + * for p-coarsening, but the cell's criterion must be lower than the + * threshold. However, different compare function objects can be supplied + * via the parameters @p compare_refine and @p compare_coarsen to impose + * different decision strategies. + * + * Each entry of the parameter @p criteria needs to correspond to an active + * cell. Parameters @p p_refine_fraction and @p p_coarsen_fraction need to be + * in the interval $[0,1]$. + * + * @note Preceeding calls of Triangulation::prepare_for_coarsening_and_refinement() + * may change refine and coarsen flags, which will ultimately change the + * results of this function. + */ + template + void + p_adaptivity_from_relative_threshold( + const hp::DoFHandler &dof_handler, + const Vector & criteria, const double p_refine_fraction = 0.5, - const double p_coarsen_fraction = 0.5); + const double p_coarsen_fraction = 0.5, + const ComparisonFunction &compare_refine = std::greater(), + const ComparisonFunction &compare_coarsen = std::less()); /** - * Adapt the finite element on cells based on the regularity of the + * Adapt which finite element to use on cells based on the regularity of the * (unknown) analytical solution. * * With an approximation of the local Sobolev regularity index $k_K$, @@ -211,17 +271,17 @@ namespace hp * * For more theoretical details see * @code{.bib} - * @article{Houston2005, - * author = {Houston, Paul and S{\"u}li, Endre}, - * title = {A note on the design of hp-adaptive finite element - * methods for elliptic partial differential equations}, - * journal = {{Computer Methods in Applied Mechanics and Engineering}}, - * volume = {194}, - * number = {2}, - * pages = {229--243}, + * @article{Ainsworth1998, + * author = {Ainsworth, Mark and Senior, Bill}, + * title = {An adaptive refinement strategy for hp-finite element + * computations}, + * journal = {{Applied Numerical Mathematics}}, + * volume = {26}, + * number = {1--2}, + * pages = {165--178}, * publisher = {Elsevier}, - * year = {2005}, - * doi = {10.1016/j.cma.2004.04.009} + * year = {1998}, + * doi = {10.1016/S0168-9274(97)00083-4} * } * @endcode * @@ -236,39 +296,19 @@ namespace hp const Vector & sobolev_indices); /** - * Adapt the finite element on cells based on their refinement history - * or rather the predicted change of their error estimates. + * Adapt which finite element to use on each cell based on how its criterion + * relates to a reference. * - * If a cell is flagged for adaptation, we will perform p-adaptation once - * the associated error indicators $\eta_{K}$ on cell $K$ satisfy - * $\eta_{K} < \eta_{K,\text{pred}}$, where the subscript $\text{pred}$ - * denotes the predicted error. This corresponds to our assumption of - * smoothness being correct, else h-adaptation is supplied. - * - * For the very first adaptation step, the user needs to decide whether h- - * or p-adaptation is supposed to happen. An h-step will be applied with - * $\eta_{K,\text{pred}} = 0$, whereas $\eta_{K,\text{pred}} = \infty$ - * ensures a p-step. The latter may be realised with - * `std::numeric_limits::max()`. + * We consider a cell for p-adaptivity if it is currently flagged for + * refinement or coarsening and its criterion successfully compares to the + * corresponding reference. Other than functions + * p_adaptivity_from_absolute_threshold() and + * p_adaptivity_from_relative_threshold(), compare function objects have to + * be provided explicitly via the parameters @p compare_refine and + * @p compare_coarsen. * - * Each entry of the parameter @p error_indicators and @p predicted_errors - * needs to correspond to an active cell. - * - * For more theoretical details see - * @code{.bib} - * @article{Melenk2001, - * author = {Melenk, Jens Markus and Wohlmuth, Barbara I.}, - * title = {{On residual-based a posteriori error estimation - * in hp-FEM}}, - * journal = {{Advances in Computational Mathematics}}, - * volume = {15}, - * number = {1}, - * pages = {311--331}, - * publisher = {Springer US}, - * year = {2001}, - * doi = {10.1023/A:1014268310921} - * } - * @endcode + * Each entry of the parameters @p criteria and @p references needs to + * correspond to an active cell. * * @note Preceeding calls of Triangulation::prepare_for_coarsening_and_refinement() * may change refine and coarsen flags, which will ultimately change the @@ -276,11 +316,12 @@ namespace hp */ template void - p_adaptivity_from_prediction( + p_adaptivity_from_reference( const hp::DoFHandler &dof_handler, - const Vector & error_indicators, - const Vector & predicted_errors); - + const Vector & criteria, + const Vector & references, + const ComparisonFunction & compare_refine, + const ComparisonFunction & compare_coarsen); /** * @} */ @@ -319,21 +360,18 @@ namespace hp * confident to say that the error will not change by sole interpolation on * the larger finite element space. * - * Further, the function assumes that the local error on a cell - * that will be refined, will lead to errors on the $2^{dim}$ - * children that are all equal, whereas local errors on siblings - * will be summed up on the parent cell in case of - * coarsening. This assumption is often not satisfied in practice: - * For example, if a cell is at a corner singularity, then the one - * child cell that ends up closest to the singularity will inherit - * the majority of the remaining error -- but this function can - * not know where the singularity will be, and consequently - * assumes equal distribution. - * - * When transferring the predicted error to the coarsened mesh, - * make sure to configure your CellDataTransfer object with - * CoarseningStrategies::sum() as a coarsening - * strategy. + * Further, the function assumes that the local error on a cell that will be + * refined, will lead to errors on the $2^{dim}$ children that are all + * equal, whereas local errors on siblings will be summed up on the parent + * cell in case of coarsening. This assumption is often not satisfied in + * practice: For example, if a cell is at a corner singularity, then the one + * child cell that ends up closest to the singularity will inherit the + * majority of the remaining error -- but this function can not know where + * the singularity will be, and consequently assumes equal distribution. + * + * When transferring the predicted error to the coarsened mesh, make sure to + * configure your CellDataTransfer object with CoarseningStrategies::sum() + * as a coarsening strategy. * * For p-adaptation, the local error is expected to converge exponentially * with the polynomial degree of the assigned finite element. Each increase @@ -366,6 +404,66 @@ namespace hp * \quad \forall K_c \text{ children of } K$ * * + * With these predicted error estimates, we are capable of adapting the + * finite element on cells based on their refinement history or rather the + * predicted change of their error estimates. + * + * If a cell is flagged for adaptation, we want to perform p-adaptation once + * the associated error indicators $\eta_{K}$ on cell $K$ satisfy + * $\eta_{K} < \eta_{K,\text{pred}}$, where the subscript $\text{pred}$ + * denotes the predicted error. This corresponds to our assumption of + * smoothness being correct, else h-adaptation is applied. We achieve this + * with the function hp::Refinement::p_adaptivity_from_criteria() and a + * function object `std::less()` for both comparator parameters. + * + * For the very first adaptation step, the user needs to decide whether h- + * or p-adaptation is supposed to happen. An h-step will be applied with + * $\eta_{K,\text{pred}} = 0$, whereas $\eta_{K,\text{pred}} = \infty$ + * ensures a p-step. The latter may be realised with + * `std::numeric_limits::infinity()`. + * + * The following code snippet demonstrates how to impose hp-adaptivity based + * on refinement history in an application: + * @code + * // [initialisation...] + * Vector predicted_error_per_cell(triangulation.n_active_cells()); + * for(unsigned int i = 0; i < triangulation.n_active_cells(); ++i) + * predicted_error_per_cell[i] = std::numeric_limits::max(); + * + * // [during each refinement step...] + * // set h-adaptivity flags + * Vector estimated_error_per_cell(triangulation.n_active_cells()); + * KellyErrorEstimator::estimate(...); + * GridRefinemet::refine_and_coarsen_fixed_fraction(...); + * + * // set p-adaptivity flags + * hp::Refinement::p_adaptivity_from_reference( + * hp_dof_handler, + * estimated_error_per_cell, + * predicted_error_per_cell, + * std::less(), + * std::less()); + * hp::Refinement::{choose|force}_p_over_h(...); + * + * // predict error for the subsequent adaptation + * triangulation.prepare_coarsening_and_refinement(); + * hp::Refinement::predict_error( + * hp_dof_handler, + * estimated_error_per_cell, + * predicted_error_per_cell); + * + * // perform adaptation + * CellDataTransfer> + * cell_data_transfer(triangulation); + * cell_data_transfer.prepare_for_coarsening_and_refinement(); + * + * triangulation.execute_coarsening_and_refinement(); + * + * Vector transferred_errors(triangulation.n_active_cells()); + * cell_data_transfer.unpack(predicted_error_per_cell, transferred_errors); + * predicted_error_per_cell = std::move(transferred_errors); + * @endcode + * * For more theoretical details see * @code{.bib} * @article{Melenk2001, diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc index ff30ad67cb..2d60af474a 100644 --- a/source/hp/refinement.cc +++ b/source/hp/refinement.cc @@ -14,6 +14,8 @@ // --------------------------------------------------------------------- +#include + #include #include @@ -84,52 +86,83 @@ namespace hp template void - p_adaptivity_from_threshold( + p_adaptivity_from_absolute_threshold( const hp::DoFHandler &dof_handler, - const Vector & smoothness_indicators, + const Vector & criteria, + const Number p_refine_threshold, + const Number p_coarsen_threshold, + const ComparisonFunction & compare_refine, + const ComparisonFunction & compare_coarsen) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + criteria.size()); + + std::vector p_flags( + dof_handler.get_triangulation().n_active_cells(), false); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned() && + ((cell->refine_flag_set() && + compare_refine(criteria[cell->active_cell_index()], + p_refine_threshold)) || + (cell->coarsen_flag_set() && + compare_coarsen(criteria[cell->active_cell_index()], + p_coarsen_threshold)))) + p_flags[cell->active_cell_index()] = true; + + p_adaptivity_from_flags(dof_handler, p_flags); + } + + + + template + void + p_adaptivity_from_relative_threshold( + const hp::DoFHandler &dof_handler, + const Vector & criteria, const double p_refine_fraction, - const double p_coarsen_fraction) + const double p_coarsen_fraction, + const ComparisonFunction & compare_refine, + const ComparisonFunction & compare_coarsen) { AssertDimension(dof_handler.get_triangulation().n_active_cells(), - smoothness_indicators.size()); + criteria.size()); Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1), dealii::GridRefinement::ExcInvalidParameterValue()); Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1), dealii::GridRefinement::ExcInvalidParameterValue()); // We first have to determine the maximal and minimal values of the - // smoothness indicators of all flagged cells. We start with the minimal - // and maximal values of all cells, a range within which the minimal and - // maximal values on cells flagged for refinement must surely lie. - Number max_smoothness_refine = - *std::min_element(smoothness_indicators.begin(), - smoothness_indicators.end()), - min_smoothness_refine = - *std::max_element(smoothness_indicators.begin(), - smoothness_indicators.end()); - Number max_smoothness_coarsen = max_smoothness_refine, - min_smoothness_coarsen = min_smoothness_refine; + // criteria of all flagged cells. We start with the minimal and maximal + // values of all cells, a range within which the minimal and maximal + // values on cells flagged for refinement must surely lie. + Number max_criterion_refine = + *std::min_element(criteria.begin(), criteria.end()), + min_criterion_refine = + *std::max_element(criteria.begin(), criteria.end()); + Number max_criterion_coarsen = max_criterion_refine, + min_criterion_coarsen = min_criterion_refine; for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned()) { if (cell->refine_flag_set()) { - max_smoothness_refine = - std::max(max_smoothness_refine, - smoothness_indicators(cell->active_cell_index())); - min_smoothness_refine = - std::min(min_smoothness_refine, - smoothness_indicators(cell->active_cell_index())); + max_criterion_refine = + std::max(max_criterion_refine, + criteria(cell->active_cell_index())); + min_criterion_refine = + std::min(min_criterion_refine, + criteria(cell->active_cell_index())); } else if (cell->coarsen_flag_set()) { - max_smoothness_coarsen = - std::max(max_smoothness_coarsen, - smoothness_indicators(cell->active_cell_index())); - min_smoothness_coarsen = - std::min(min_smoothness_coarsen, - smoothness_indicators(cell->active_cell_index())); + max_criterion_coarsen = + std::max(max_criterion_coarsen, + criteria(cell->active_cell_index())); + min_criterion_coarsen = + std::min(min_criterion_coarsen, + criteria(cell->active_cell_index())); } } @@ -137,47 +170,37 @@ namespace hp dynamic_cast *>( &dof_handler.get_triangulation())) { - max_smoothness_refine = - Utilities::MPI::max(max_smoothness_refine, + max_criterion_refine = + Utilities::MPI::max(max_criterion_refine, parallel_tria->get_communicator()); - min_smoothness_refine = - Utilities::MPI::min(min_smoothness_refine, + min_criterion_refine = + Utilities::MPI::min(min_criterion_refine, parallel_tria->get_communicator()); - max_smoothness_coarsen = - Utilities::MPI::max(max_smoothness_coarsen, + max_criterion_coarsen = + Utilities::MPI::max(max_criterion_coarsen, parallel_tria->get_communicator()); - min_smoothness_coarsen = - Utilities::MPI::min(min_smoothness_coarsen, + min_criterion_coarsen = + Utilities::MPI::min(min_criterion_coarsen, parallel_tria->get_communicator()); } // Absent any better strategies, we will set the threshold by linear // interpolation for both classes of cells individually. - const Number threshold_smoothness_refine = - min_smoothness_refine + + const Number threshold_refine = + min_criterion_refine + p_refine_fraction * - (max_smoothness_refine - min_smoothness_refine), - threshold_smoothness_coarsen = - min_smoothness_coarsen + + (max_criterion_refine - min_criterion_refine), + threshold_coarsen = + min_criterion_coarsen + p_coarsen_fraction * - (max_smoothness_coarsen - min_smoothness_coarsen); - - // We then compare each cell's smoothness indicator with the corresponding - // threshold. - std::vector p_flags( - dof_handler.get_triangulation().n_active_cells(), false); - - for (const auto &cell : dof_handler.active_cell_iterators()) - if (cell->is_locally_owned() && - ((cell->refine_flag_set() && - (smoothness_indicators(cell->active_cell_index()) > - threshold_smoothness_refine)) || - (cell->coarsen_flag_set() && - (smoothness_indicators(cell->active_cell_index()) < - threshold_smoothness_coarsen)))) - p_flags[cell->active_cell_index()] = true; - - p_adaptivity_from_flags(dof_handler, p_flags); + (max_criterion_coarsen - min_criterion_coarsen); + + p_adaptivity_from_absolute_threshold(dof_handler, + criteria, + threshold_refine, + threshold_coarsen, + compare_refine, + compare_coarsen); } @@ -235,24 +258,29 @@ namespace hp template void - p_adaptivity_from_prediction( + p_adaptivity_from_reference( const hp::DoFHandler &dof_handler, - const Vector & error_indicators, - const Vector & predicted_errors) + const Vector & criteria, + const Vector & references, + const ComparisonFunction & compare_refine, + const ComparisonFunction & compare_coarsen) { AssertDimension(dof_handler.get_triangulation().n_active_cells(), - error_indicators.size()); + criteria.size()); AssertDimension(dof_handler.get_triangulation().n_active_cells(), - predicted_errors.size()); + references.size()); std::vector p_flags( dof_handler.get_triangulation().n_active_cells(), false); for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned() && - ((cell->refine_flag_set() || cell->coarsen_flag_set()) && - (error_indicators[cell->active_cell_index()] < - predicted_errors[cell->active_cell_index()]))) + ((cell->refine_flag_set() && + compare_refine(criteria[cell->active_cell_index()], + references[cell->active_cell_index()])) || + (cell->coarsen_flag_set() && + compare_coarsen(criteria[cell->active_cell_index()], + references[cell->active_cell_index()])))) p_flags[cell->active_cell_index()] = true; p_adaptivity_from_flags(dof_handler, p_flags); @@ -284,8 +312,6 @@ namespace hp for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned()) { - const unsigned int active_cell_index = cell->active_cell_index(); - if (cell->future_fe_index_set()) // p adaptation { Assert(!cell->refine_flag_set() && !cell->coarsen_flag_set(), @@ -297,8 +323,8 @@ namespace hp .degree - cell->get_fe().degree; - predicted_errors[active_cell_index] = - error_indicators[active_cell_index] * + predicted_errors[cell->active_cell_index()] = + error_indicators[cell->active_cell_index()] * std::pow(gamma_p, degree_difference); } else if (cell->refine_flag_set()) // h refinement @@ -309,20 +335,20 @@ namespace hp ExcMessage( "Error prediction is only valid for isotropic refinement!")); - predicted_errors[active_cell_index] = - error_indicators[active_cell_index] * + predicted_errors[cell->active_cell_index()] = + error_indicators[cell->active_cell_index()] * (gamma_h * std::pow(.5, dim + cell->get_fe().degree)); } else if (cell->coarsen_flag_set()) // h coarsening { - predicted_errors[active_cell_index] = - error_indicators[active_cell_index] / + predicted_errors[cell->active_cell_index()] = + error_indicators[cell->active_cell_index()] / (gamma_h * std::pow(.5, cell->get_fe().degree)); } else // no changes { - predicted_errors[active_cell_index] = - error_indicators[active_cell_index] * gamma_n; + predicted_errors[cell->active_cell_index()] = + error_indicators[cell->active_cell_index()] * gamma_n; } } } @@ -390,7 +416,8 @@ namespace hp // drop all h coarsening flags. for (unsigned int child_index = 0; child_index < n_children; ++child_index) - parent->child(child_index)->clear_coarsen_flag(); + if (parent->child(child_index)->active()) + parent->child(child_index)->clear_coarsen_flag(); } } } diff --git a/source/hp/refinement.inst.in b/source/hp/refinement.inst.in index 2980563f27..f0ffb98195 100644 --- a/source/hp/refinement.inst.in +++ b/source/hp/refinement.inst.in @@ -52,13 +52,26 @@ for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS; namespace Refinement \{ template void - p_adaptivity_from_threshold( + p_adaptivity_from_absolute_threshold( + const hp::DoFHandler &, + const Vector &, + const S, + const S, + const std::function &, + const std::function &); + + template void + p_adaptivity_from_relative_threshold( const hp::DoFHandler &, const Vector &, const double, - const double); + const double, + const std::function &, + const std::function &); template void p_adaptivity_from_regularity &); template void - p_adaptivity_from_prediction( + p_adaptivity_from_reference( const hp::DoFHandler &, const Vector &, - const Vector &); + const Vector &, + const std::function &, + const std::function &); template void predict_error( diff --git a/tests/hp/error_prediction.cc b/tests/hp/error_prediction.cc index f8837eecbc..8d76702f56 100644 --- a/tests/hp/error_prediction.cc +++ b/tests/hp/error_prediction.cc @@ -101,6 +101,14 @@ test() deallog << std::endl; } + // ----- check feature ----- + hp::Refinement::p_adaptivity_from_reference( + dh, + error_indicators, + predicted_errors, + /*compare_refine=*/std::less(), + /*compare_coarsen=*/std::less()); + deallog << "OK" << std::endl; } diff --git a/tests/hp/p_adaptivity_flags.cc b/tests/hp/p_adaptivity_flags.cc index 67a591b4af..ab65ac8d98 100644 --- a/tests/hp/p_adaptivity_flags.cc +++ b/tests/hp/p_adaptivity_flags.cc @@ -15,7 +15,7 @@ -// validate algorithms that will flag cells for p adaptivity +// validate algorithms that will flag cells for p-adaptivity #include @@ -101,7 +101,7 @@ test() validate(tria, dh); } - deallog << "full p adaptivity" << std::endl; + deallog << "full p-adaptivity" << std::endl; { Triangulation tria; hp::DoFHandler dh; @@ -117,7 +117,7 @@ test() // Ultimately, the first quarter of all cells will be flagged for // p refinement, and the last quarter for p coarsening. - deallog << "p adaptivity from flags" << std::endl; + deallog << "p-adaptivity from flags" << std::endl; { Triangulation tria; hp::DoFHandler dh; @@ -134,29 +134,54 @@ test() validate(tria, dh); } - deallog << "p adaptivity from threshold" << std::endl; + deallog << "p-adaptivity from absolute threshold" << std::endl; { Triangulation tria; hp::DoFHandler dh; setup(tria, dh, fes); unsigned int n_active = tria.n_active_cells(); - Vector smoothness_indicators(n_active); + Vector indicators(n_active); for (unsigned int i = 0; i < n_active; ++i) { if (i < .25 * n_active) - smoothness_indicators[i] = 2.; + indicators[i] = 2.; else if (i < .75 * n_active) - smoothness_indicators[i] = 1.; + indicators[i] = 1.; else - smoothness_indicators[i] = 0.; + indicators[i] = 0.; } - hp::Refinement::p_adaptivity_from_threshold(dh, smoothness_indicators); + hp::Refinement::p_adaptivity_from_absolute_threshold(dh, + indicators, + 1 + 1e-4, + 1 - 1e-4); validate(tria, dh); } - deallog << "p adaptivity from regularity" << std::endl; + deallog << "p-adaptivity from relative threshold" << std::endl; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + unsigned int n_active = tria.n_active_cells(); + Vector indicators(n_active); + for (unsigned int i = 0; i < n_active; ++i) + { + if (i < .25 * n_active) + indicators[i] = 2.; + else if (i < .75 * n_active) + indicators[i] = 1.; + else + indicators[i] = 0.; + } + hp::Refinement::p_adaptivity_from_relative_threshold(dh, indicators); + + validate(tria, dh); + } + + deallog << "p-adaptivity from regularity" << std::endl; { Triangulation tria; hp::DoFHandler dh; @@ -178,35 +203,34 @@ test() validate(tria, dh); } - deallog << "p adaptivity from prediction" << std::endl; + deallog << "p-adaptivity from reference" << std::endl; { Triangulation tria; hp::DoFHandler dh; setup(tria, dh, fes); unsigned int n_active = tria.n_active_cells(); - Vector predicted_errors(n_active), error_estimates(n_active); + Vector references(n_active), criteria(n_active); for (unsigned int i = 0; i < n_active; ++i) { if (i < .25 * n_active) { - predicted_errors[i] = 1. + 1e-4; - error_estimates[i] = 1.; + references[i] = 1. + 1e-4; + criteria[i] = 1.; } else if (i < .75 * n_active) { - predicted_errors[i] = 1.; - error_estimates[i] = 1.; + references[i] = 1.; + criteria[i] = 1.; } else { - predicted_errors[i] = 1. + 1e-4; - error_estimates[i] = 1.; + references[i] = 1. - 1e-4; + criteria[i] = 1.; } } - hp::Refinement::p_adaptivity_from_prediction(dh, - error_estimates, - predicted_errors); + hp::Refinement::p_adaptivity_from_reference( + dh, criteria, references, std::less(), std::greater()); validate(tria, dh); } diff --git a/tests/hp/p_adaptivity_flags.output b/tests/hp/p_adaptivity_flags.output index 6061a757bb..bb046b4c0e 100644 --- a/tests/hp/p_adaptivity_flags.output +++ b/tests/hp/p_adaptivity_flags.output @@ -1,40 +1,46 @@ DEAL:1d::starting situation: ncells: 4 DEAL:1d:: fe_indices: 1 1 1 1 -DEAL:1d::full p adaptivity +DEAL:1d::full p-adaptivity DEAL:1d:: fe_indices: 2 2 0 0 -DEAL:1d::p adaptivity from flags +DEAL:1d::p-adaptivity from flags DEAL:1d:: fe_indices: 2 1 1 0 -DEAL:1d::p adaptivity from threshold +DEAL:1d::p-adaptivity from absolute threshold DEAL:1d:: fe_indices: 2 1 1 0 -DEAL:1d::p adaptivity from regularity +DEAL:1d::p-adaptivity from relative threshold DEAL:1d:: fe_indices: 2 1 1 0 -DEAL:1d::p adaptivity from prediction +DEAL:1d::p-adaptivity from regularity +DEAL:1d:: fe_indices: 2 1 1 0 +DEAL:1d::p-adaptivity from reference DEAL:1d:: fe_indices: 2 1 1 0 DEAL:1d::OK DEAL:2d::starting situation: ncells: 16 DEAL:2d:: fe_indices: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -DEAL:2d::full p adaptivity +DEAL:2d::full p-adaptivity DEAL:2d:: fe_indices: 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 -DEAL:2d::p adaptivity from flags +DEAL:2d::p-adaptivity from flags +DEAL:2d:: fe_indices: 2 2 2 2 1 1 1 1 1 1 1 1 0 0 0 0 +DEAL:2d::p-adaptivity from absolute threshold DEAL:2d:: fe_indices: 2 2 2 2 1 1 1 1 1 1 1 1 0 0 0 0 -DEAL:2d::p adaptivity from threshold +DEAL:2d::p-adaptivity from relative threshold DEAL:2d:: fe_indices: 2 2 2 2 1 1 1 1 1 1 1 1 0 0 0 0 -DEAL:2d::p adaptivity from regularity +DEAL:2d::p-adaptivity from regularity DEAL:2d:: fe_indices: 2 2 2 2 1 1 1 1 1 1 1 1 0 0 0 0 -DEAL:2d::p adaptivity from prediction +DEAL:2d::p-adaptivity from reference DEAL:2d:: fe_indices: 2 2 2 2 1 1 1 1 1 1 1 1 0 0 0 0 DEAL:2d::OK DEAL:3d::starting situation: ncells: 64 DEAL:3d:: fe_indices: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -DEAL:3d::full p adaptivity +DEAL:3d::full p-adaptivity DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:3d::p adaptivity from flags +DEAL:3d::p-adaptivity from flags +DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:3d::p-adaptivity from absolute threshold DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:3d::p adaptivity from threshold +DEAL:3d::p-adaptivity from relative threshold DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:3d::p adaptivity from regularity +DEAL:3d::p-adaptivity from regularity DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:3d::p adaptivity from prediction +DEAL:3d::p-adaptivity from reference DEAL:3d:: fe_indices: 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:3d::OK diff --git a/tests/mpi/petsc_step-27.cc b/tests/mpi/petsc_step-27.cc index 1fba613fcf..78c3e8d12a 100644 --- a/tests/mpi/petsc_step-27.cc +++ b/tests/mpi/petsc_step-27.cc @@ -383,10 +383,10 @@ namespace Step27 parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( triangulation, estimated_error_per_cell, 0.3, 0.03); - hp::Refinement::p_adaptivity_from_threshold(dof_handler, - smoothness_indicators, - 0.5, - 0.); + hp::Refinement::p_adaptivity_from_relative_threshold(dof_handler, + smoothness_indicators, + 0.5, + 0.); hp::Refinement::choose_p_over_h(dof_handler); triangulation.execute_coarsening_and_refinement(); diff --git a/tests/mpi/trilinos_step-27.cc b/tests/mpi/trilinos_step-27.cc index ef0886d92b..7741928580 100644 --- a/tests/mpi/trilinos_step-27.cc +++ b/tests/mpi/trilinos_step-27.cc @@ -385,10 +385,10 @@ namespace Step27 parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( triangulation, estimated_error_per_cell, 0.3, 0.03); - hp::Refinement::p_adaptivity_from_threshold(dof_handler, - smoothness_indicators, - 0.5, - 0.); + hp::Refinement::p_adaptivity_from_relative_threshold(dof_handler, + smoothness_indicators, + 0.5, + 0.); hp::Refinement::choose_p_over_h(dof_handler); triangulation.execute_coarsening_and_refinement(); -- 2.39.5