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/base/config.h>
-#include <vector>
+#include <deal.II/base/template_constraints.h>
+
+#include <functional>
+
DEAL_II_NAMESPACE_OPEN
+
// forward declarations
#ifndef DOXYGEN
template <typename Number>
}
#endif
+
namespace hp
{
/**
* @code
* // step 1: flag cells for refinement or coarsening
* Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
- * KellyErrorEstimator<dim>::estimate (hp_dof_handler,
- * QGauss<dim-1> (quadrature_points),
- * typename FunctionMap<dim>::type(),
- * solution,
- * estimated_error_per_cell);
- * GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
- * estimated_error_per_cell,
- * top_fraction,
- * bottom_fraction);
+ * KellyErrorEstimator<dim>::estimate(hp_dof_handler,
+ * QGauss<dim-1> (quadrature_points),
+ * typename FunctionMap<dim>::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
* ...
*/
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 <typename Number>
+ using ComparisonFunction =
+ std::function<bool(const typename identity<Number>::type &,
+ const typename identity<Number>::type &)>;
+
/**
* @name Setting p-adaptivity flags
* @{
full_p_adaptivity(const hp::DoFHandler<dim, spacedim> &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.
const std::vector<bool> & 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
*/
template <int dim, typename Number, int spacedim>
void
- p_adaptivity_from_threshold(
+ p_adaptivity_from_absolute_threshold(
const hp::DoFHandler<dim, spacedim> &dof_handler,
- const Vector<Number> & smoothness_indicators,
+ const Vector<Number> & criteria,
+ const Number p_refine_threshold,
+ const Number p_coarsen_threshold,
+ const ComparisonFunction<Number> &compare_refine = std::greater<Number>(),
+ const ComparisonFunction<Number> &compare_coarsen = std::less<Number>());
+
+ /**
+ * 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_relative_threshold(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & criteria,
const double p_refine_fraction = 0.5,
- const double p_coarsen_fraction = 0.5);
+ const double p_coarsen_fraction = 0.5,
+ const ComparisonFunction<Number> &compare_refine = std::greater<Number>(),
+ const ComparisonFunction<Number> &compare_coarsen = std::less<Number>());
/**
- * 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$,
*
* 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
*
const Vector<Number> & 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
*/
template <int dim, typename Number, int spacedim>
void
- p_adaptivity_from_prediction(
+ p_adaptivity_from_reference(
const hp::DoFHandler<dim, spacedim> &dof_handler,
- const Vector<Number> & error_indicators,
- const Vector<Number> & predicted_errors);
-
+ const Vector<Number> & criteria,
+ const Vector<Number> & references,
+ const ComparisonFunction<Number> & compare_refine,
+ const ComparisonFunction<Number> & compare_coarsen);
/**
* @}
*/
* 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
* \quad \forall K_c \text{ children of } K$
* </table>
*
+ * 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<Number>()` 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<float> 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<float>::max();
+ *
+ * // [during each refinement step...]
+ * // set h-adaptivity flags
+ * Vector<float> 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<float>(),
+ * std::less<float>());
+ * 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<dim, spacedim, Vector<float>>
+ * cell_data_transfer(triangulation);
+ * cell_data_transfer.prepare_for_coarsening_and_refinement();
+ *
+ * triangulation.execute_coarsening_and_refinement();
+ *
+ * Vector<float> 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,
// ---------------------------------------------------------------------
+#include <deal.II/base/config.h>
+
#include <deal.II/base/mpi.h>
#include <deal.II/distributed/tria_base.h>
template <int dim, typename Number, int spacedim>
void
- p_adaptivity_from_threshold(
+ p_adaptivity_from_absolute_threshold(
const hp::DoFHandler<dim, spacedim> &dof_handler,
- const Vector<Number> & smoothness_indicators,
+ const Vector<Number> & criteria,
+ const Number p_refine_threshold,
+ const Number p_coarsen_threshold,
+ const ComparisonFunction<Number> & compare_refine,
+ const ComparisonFunction<Number> & compare_coarsen)
+ {
+ AssertDimension(dof_handler.get_triangulation().n_active_cells(),
+ criteria.size());
+
+ std::vector<bool> 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_relative_threshold(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & criteria,
const double p_refine_fraction,
- const double p_coarsen_fraction)
+ const double p_coarsen_fraction,
+ const ComparisonFunction<Number> & compare_refine,
+ const ComparisonFunction<Number> & 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()));
}
}
dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
&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<bool> 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);
}
template <int dim, typename Number, int spacedim>
void
- p_adaptivity_from_prediction(
+ p_adaptivity_from_reference(
const hp::DoFHandler<dim, spacedim> &dof_handler,
- const Vector<Number> & error_indicators,
- const Vector<Number> & predicted_errors)
+ const Vector<Number> & criteria,
+ const Vector<Number> & references,
+ const ComparisonFunction<Number> & compare_refine,
+ const ComparisonFunction<Number> & 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<bool> 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);
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(),
.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
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;
}
}
}
// 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();
}
}
}
namespace Refinement
\{
template void
- p_adaptivity_from_threshold<deal_II_dimension,
- S,
- deal_II_space_dimension>(
+ p_adaptivity_from_absolute_threshold<deal_II_dimension,
+ S,
+ deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const Vector<S> &,
+ const S,
+ const S,
+ const std::function<bool(const S &, const S &)> &,
+ const std::function<bool(const S &, const S &)> &);
+
+ template void
+ p_adaptivity_from_relative_threshold<deal_II_dimension,
+ S,
+ deal_II_space_dimension>(
const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const Vector<S> &,
const double,
- const double);
+ const double,
+ const std::function<bool(const S &, const S &)> &,
+ const std::function<bool(const S &, const S &)> &);
template void
p_adaptivity_from_regularity<deal_II_dimension,
const Vector<S> &);
template void
- p_adaptivity_from_prediction<deal_II_dimension,
- S,
- deal_II_space_dimension>(
+ p_adaptivity_from_reference<deal_II_dimension,
+ S,
+ deal_II_space_dimension>(
const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const Vector<S> &,
- const Vector<S> &);
+ const Vector<S> &,
+ const std::function<bool(const S &, const S &)> &,
+ const std::function<bool(const S &, const S &)> &);
template void
predict_error<deal_II_dimension, S, deal_II_space_dimension>(
deallog << std::endl;
}
+ // ----- check feature -----
+ hp::Refinement::p_adaptivity_from_reference(
+ dh,
+ error_indicators,
+ predicted_errors,
+ /*compare_refine=*/std::less<float>(),
+ /*compare_coarsen=*/std::less<float>());
+
deallog << "OK" << std::endl;
}
-// validate algorithms that will flag cells for p adaptivity
+// validate algorithms that will flag cells for p-adaptivity
#include <deal.II/base/geometry_info.h>
validate(tria, dh);
}
- deallog << "full p adaptivity" << std::endl;
+ deallog << "full p-adaptivity" << std::endl;
{
Triangulation<dim> tria;
hp::DoFHandler<dim> dh;
// 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<dim> tria;
hp::DoFHandler<dim> dh;
validate(tria, dh);
}
- deallog << "p adaptivity from threshold" << std::endl;
+ deallog << "p-adaptivity from absolute threshold" << std::endl;
{
Triangulation<dim> tria;
hp::DoFHandler<dim> dh;
setup(tria, dh, fes);
unsigned int n_active = tria.n_active_cells();
- Vector<double> smoothness_indicators(n_active);
+ Vector<double> 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<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ unsigned int n_active = tria.n_active_cells();
+ Vector<double> 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<dim> tria;
hp::DoFHandler<dim> dh;
validate(tria, dh);
}
- deallog << "p adaptivity from prediction" << std::endl;
+ deallog << "p-adaptivity from reference" << std::endl;
{
Triangulation<dim> tria;
hp::DoFHandler<dim> dh;
setup(tria, dh, fes);
unsigned int n_active = tria.n_active_cells();
- Vector<double> predicted_errors(n_active), error_estimates(n_active);
+ Vector<double> 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<double>(), std::greater<double>());
validate(tria, dh);
}
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
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();
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();
projects / dealii.git / commitdiff