From: Marc Fehling Date: Fri, 7 Jun 2019 10:01:46 +0000 (+0200) Subject: Introduced hp::Refinement namespace. X-Git-Tag: v9.2.0-rc1~1407^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c32bb659acd6474fc9e696301d0123824c9d969d;p=dealii.git Introduced hp::Refinement namespace. --- diff --git a/doc/news/changes/minor/20190607Fehling b/doc/news/changes/minor/20190607Fehling new file mode 100644 index 0000000000..24635b882e --- /dev/null +++ b/doc/news/changes/minor/20190607Fehling @@ -0,0 +1,3 @@ +New: Namespace hp::Refinement offering decision tools for p adaptivity. +
+(Marc Fehling, 2019/06/07) diff --git a/include/deal.II/hp/refinement.h b/include/deal.II/hp/refinement.h new file mode 100644 index 0000000000..7101ab19f2 --- /dev/null +++ b/include/deal.II/hp/refinement.h @@ -0,0 +1,355 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_hp_refinement_h +#define dealii_hp_refinement_h + + +#include + +DEAL_II_NAMESPACE_OPEN + +// forward declarations +template +class Vector; + +namespace hp +{ + template + class DoFHandler; +} + +namespace hp +{ + /** + * We supply adaptive methods to align computational ressources with the + * complexity of the numerical solution. Error estimates are an appropriate + * means of determining where adjustments need to be made. + * + * However with hp adaptivity, we have two ways to realise these adjustments: + * For irregular solutions, h adaptive methods which dynamically assign cell + * sizes tend to reduce the approximation error, while for smooth solutions p + * adaptive methods are better suited in which function spaces will be + * selected dynamically. This namespace offers all tools to decide which type + * of adaptive methods to apply. + * + *

Usage

+ * + * To successfully apply hp adaptive methods, we recommend the following + * workflow: + *
    + *
  1. A suitable error estimate is the basis for any kind of adaptive method. + * Similar to pure grid refinement, we will determine error estimates in the + * usual way (i.e. KellyErrorEstimator) and mark cells for refinement or + * coarsening (i.e. GridRefinement). + * + * Calling Triangulation::execute_coarsening_and_refinement() at this stage + * will perform pure grid refinement as expected. + * + *
  2. Once all refinement and coarsening flags have been distributed on the + * mesh, we may determine if those qualify for p adaptive methods. + * Corresponding functions will set @p future_fe_indices on top of the + * refinement and coarsening flags if they fulfil a certain criterion. + * + * In case of refinement, the superordinate element of the underlying + * hp::FECollection will be assigned as the future finite element. + * Correspondingly, the subordinate element will be selected for coarsening. + * + * Triangulation::execute_coarsening_and_refinement() will now supply both h + * and p adaptive methods independently. + * + *
  3. Right now, there may be cells scheduled for both h and p adaptation. + * If we do not want to impose both methods at once, we need to decide which + * one to pick for each cell individually and unambiguously. Since grid + * refinement will be imposed by default and we only determine qualification + * for p adaptivity on top, we will always decide in favour of p adaptive + * methods. + * + * Calling Triangulation::execute_coarsening_and_refinement() will now perform + * either h or p adaptive methods uniquely on each cell. + * + *
  4. Up to this point, each cell knows its destiny in terms of adaptivity + * We can now move on to prepare all data structures to be transferred across + * mesh changes. Previously set refinement and coarsening flags as well as + * @p future_fe_indices will be used to update the data accordingly. + *
+ * + * As an example, a realisation of pure p adaptive methods would look like the + * following: + * @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); + * + * // step 2: set future finite element indices on flagged cells + * 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); + * + * // step 4: prepare solutions to be transferred + * ... + * + * triangulation.execute_coarsening_and_refinement(); + * @endcode + * + * @ingroup hp + * @author Marc Fehling 2019 + */ + namespace Refinement + { + /** + * @name Setting p adaptivity flags + * @{ + */ + + /** + * Each cell flagged for h refinement will also be flagged for p refinement. + * The same applies to coarsening. + * + * @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 + 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. + * + * Each entry of the parameter @p p_flags 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 + * results of this function. + */ + template + void + p_adaptivity_from_flags(const hp::DoFHandler &dof_handler, + 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]$. + * + * @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_threshold( + const hp::DoFHandler &dof_handler, + const Vector & smoothness_indicators, + const double p_refine_fraction = 0.5, + const double p_coarsen_fraction = 0.5); + + /** + * Adapt the finite element on cells based on the regularity of the + * (unknown) analytical solution. + * + * With an approximation of the local Sobolev regularity index $k_K$, + * we may assess to which finite element space our local solution on cell + * $K$ belongs. Since the regularity index is only an etimate, we won't + * use it to assign the finite element space directly, but rather consider + * it as an indicator for adaptation. If a cell is flagged for refinement, + * we will perform p refinement once it satisfies + * $k_K > p_{K,\text{super}}$, where $p_{K,\text{super}}$ is + * the polynomial degree of the finite element superordinate to the + * currently active element on cell $K$. In case of coarsening, the + * criterion $k_K < p_{K,\text{sub}}$ has to be met, with + * $p_{K,\text{sub}}$ the degree of the subordinate element. + * + * Each entry of the parameter @p sobolev_indices needs to correspond + * to an active cell. + * + * 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}, + * publisher = {Elsevier}, + * year = {2005}, + * doi = {10.1016/j.cma.2004.04.009} + * } + * @endcode + * + * @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_regularity( + const hp::DoFHandler &dof_handler, + const Vector & sobolev_indices); + + /** + * Adapt 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 will perform p adaptation once + * the associated error indicators $\eta_{K}^2$ on cell $K$ satisfy + * $\eta_{K}^2 < \eta_{K,\text{pred}}^2$, 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 adapation step, the user needs to decide whether h or + * p adapatation 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()`. + * + * 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 + * + * @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_prediction( + const hp::DoFHandler &dof_handler, + const Vector & error_indicators, + const Vector & predicted_errors); + + /** + * @} + */ + + /** + * @name Decide between h and p adaptivity + * @{ + */ + + /** + * Choose p adaptivity over h adaptivity in any case. + * + * Removes all refine and coarsen flags on cells that have a + * @p future_fe_index assigned. + * + * @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 + force_p_over_h(const hp::DoFHandler &dof_handler); + + /** + * Choose p adaptivity over h adaptivity whenever it is invoked on all + * related cells. + * + * In case of refinement, information about finite elements will be + * inherited. Thus we will prefer p refinement over h refinement whenever + * desired, i.e. clear the refine flag and supply a corresponding + * @p future_fe_index. + * + * However for coarsening, we follow a different approach. Flagging a cell + * for h coarsening does not ultimately mean that it will be coarsened. Only + * if a cell and all of its siblings are flagged, they will be merged into + * their parent cell. If we consider p coarsening on top, we must decide for + * all siblings together how they will be coarsened. We distinguish between + * three different cases: + *
    + *
  1. Not all siblings flagged for coarsening: p coarsening
    + * We keep the @p future_fe_indices and clear the coarsen flags + * on all siblings. + *
  2. All siblings flagged for coarsening, but not all for + * p adaptation: h coarsening
    + * We keep the coarsen flags and clear all @p future_fe_indices + * on all siblings. + *
  3. All siblings flagged for coarsening and p adaptation: p coarsening
    + * We keep the @p future_fe_indices and clear the coarsen flags + * on all siblings. + *
+ * + * @note The function Triangulation::prepare_coarsening_and_refinement() + * will clean up all h coarsening flags if they are not shared among + * all siblings. In the hp case, we need to bring forward this decision: + * If the cell will not be coarsened, but qualifies for p adaptivity, + * we have to set all flags accordingly. So this function anticipates + * the decision that Triangulation::prepare_coarsening_and_refinement() + * would have made later on. + * + * @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 + choose_p_over_h(const hp::DoFHandler &dof_handler); + + /** + * @} + */ + } // namespace Refinement +} // namespace hp + + +DEAL_II_NAMESPACE_CLOSE + +#endif // dealii_hp_refinement_h diff --git a/source/hp/CMakeLists.txt b/source/hp/CMakeLists.txt index 4da1819140..d1fa66f018 100644 --- a/source/hp/CMakeLists.txt +++ b/source/hp/CMakeLists.txt @@ -21,6 +21,7 @@ SET(_unity_include_src fe_collection.cc fe_values.cc mapping_collection.cc + refinement.cc ) SET(_separate_src @@ -41,6 +42,7 @@ SET(_inst fe_collection.inst.in fe_values.inst.in mapping_collection.inst.in + refinement.inst.in ) FILE(GLOB _header diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc new file mode 100644 index 0000000000..a8ae27ca22 --- /dev/null +++ b/source/hp/refinement.cc @@ -0,0 +1,367 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include + +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace hp +{ + namespace Refinement + { + /** + * Setting p adaptivity flags + */ + template + void + full_p_adaptivity(const hp::DoFHandler &dof_handler) + { + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->refine_flag_set()) + { + const unsigned int super_fe_index = + dof_handler.get_fe_collection().next_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already most superordinate element. + if (super_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(super_fe_index); + } + + if (cell->coarsen_flag_set()) + { + const unsigned int sub_fe_index = + dof_handler.get_fe_collection().previous_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already least subordinate element. + if (sub_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(sub_fe_index); + } + } + } + + + + template + void + p_adaptivity_from_flags(const hp::DoFHandler &dof_handler, + const std::vector & p_flags) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + p_flags.size()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->refine_flag_set() && p_flags[cell->active_cell_index()]) + { + const unsigned int super_fe_index = + dof_handler.get_fe_collection().next_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already most superordinate element. + if (super_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(super_fe_index); + } + + if (cell->coarsen_flag_set() && p_flags[cell->active_cell_index()]) + { + const unsigned int sub_fe_index = + dof_handler.get_fe_collection().previous_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already least subordinate element. + if (sub_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(sub_fe_index); + } + } + } + + + + template + void + p_adaptivity_from_threshold( + const hp::DoFHandler &dof_handler, + const Vector & smoothness_indicators, + const double p_refine_fraction, + const double p_coarsen_fraction) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + smoothness_indicators.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_coarsen; + + 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())); + } + 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())); + } + } + + // 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 + + p_refine_fraction * + (max_smoothness_refine - min_smoothness_refine), + threshold_smoothness_coarsen = + min_smoothness_coarsen + + p_coarsen_fraction * + (max_smoothness_coarsen - min_smoothness_coarsen); + + // We then compare each cell's smoothness indicator with the corresponding + // threshold. + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->refine_flag_set() && + (smoothness_indicators(cell->active_cell_index()) > + threshold_smoothness_refine)) + { + const unsigned int super_fe_index = + dof_handler.get_fe_collection().next_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already most superordinate element. + if (super_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(super_fe_index); + } + if (cell->coarsen_flag_set() && + (smoothness_indicators(cell->active_cell_index()) < + threshold_smoothness_coarsen)) + { + const unsigned int sub_fe_index = + dof_handler.get_fe_collection().previous_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already least subordinate element. + if (sub_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(sub_fe_index); + } + } + } + + + + template + void + p_adaptivity_from_regularity( + const hp::DoFHandler &dof_handler, + const Vector & sobolev_indices) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + sobolev_indices.size()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->refine_flag_set()) + { + const unsigned int super_fe_index = + dof_handler.get_fe_collection().next_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already most superordinate element. + if (super_fe_index != cell->active_fe_index()) + { + const unsigned int super_fe_degree = + dof_handler.get_fe_collection()[super_fe_index].degree; + + if (sobolev_indices[cell->active_cell_index()] > + super_fe_degree) + cell->set_future_fe_index(super_fe_index); + } + } + if (cell->coarsen_flag_set()) + { + const unsigned int sub_fe_index = + dof_handler.get_fe_collection().previous_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already least subordinate element. + if (sub_fe_index != cell->active_fe_index()) + { + const unsigned int sub_fe_degree = + dof_handler.get_fe_collection()[sub_fe_index].degree; + + if (sobolev_indices[cell->active_cell_index()] < + sub_fe_degree) + cell->set_future_fe_index(sub_fe_index); + } + } + } + } + + + + template + void + p_adaptivity_from_prediction( + const hp::DoFHandler &dof_handler, + const Vector & error_indicators, + const Vector & predicted_errors) + { + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + error_indicators.size()); + AssertDimension(dof_handler.get_triangulation().n_active_cells(), + predicted_errors.size()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + if (cell->refine_flag_set() && + (error_indicators[cell->active_cell_index()] < + predicted_errors[cell->active_cell_index()])) + { + const unsigned int super_fe_index = + dof_handler.get_fe_collection().next_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already most superordinate element. + if (super_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(super_fe_index); + } + + if (cell->coarsen_flag_set() && + (error_indicators[cell->active_cell_index()] < + predicted_errors[cell->active_cell_index()])) + { + const unsigned int sub_fe_index = + dof_handler.get_fe_collection().previous_in_hierarchy( + cell->active_fe_index()); + + // Reject update if already least subordinate element. + if (sub_fe_index != cell->active_fe_index()) + cell->set_future_fe_index(sub_fe_index); + } + } + } + + + + /** + * Decide between h and p adaptivity + */ + template + void + force_p_over_h(const hp::DoFHandler &dof_handler) + { + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned() && cell->future_fe_index_set()) + { + cell->clear_refine_flag(); + cell->clear_coarsen_flag(); + } + } + + + + template + void + choose_p_over_h(const hp::DoFHandler &dof_handler) + { + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned() && cell->future_fe_index_set()) + { + cell->clear_refine_flag(); + + // A cell will only be coarsened into its parent if all of its + // siblings are flagged for h coarsening as well. We must take this + // into account for our decision whether we would like to impose h + // or p adaptivity. + if (cell->coarsen_flag_set()) + { + const auto & parent = cell->parent(); + const unsigned int n_children = parent->n_children(); + + unsigned int h_flagged_children = 0, p_flagged_children = 0; + for (unsigned int child_index = 0; child_index < n_children; + ++child_index) + { + const auto &child = parent->child(child_index); + if (child->active()) + { + if (child->coarsen_flag_set()) + ++h_flagged_children; + if (child->future_fe_index_set()) + ++p_flagged_children; + } + } + + if (h_flagged_children == n_children && + p_flagged_children != n_children) + // Perform pure h coarsening and + // drop all p adaptation flags. + for (unsigned int child_index = 0; child_index < n_children; + ++child_index) + parent->child(child_index)->clear_future_fe_index(); + else + // Perform p adaptation on all children and + // drop all h coarsening flags. + for (unsigned int child_index = 0; child_index < n_children; + ++child_index) + parent->child(child_index)->clear_coarsen_flag(); + } + } + } + } // namespace Refinement +} // namespace hp + + +// explicit instantiations +#include "refinement.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/hp/refinement.inst.in b/source/hp/refinement.inst.in new file mode 100644 index 0000000000..480c873bd6 --- /dev/null +++ b/source/hp/refinement.inst.in @@ -0,0 +1,80 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace hp + \{ + namespace Refinement + \{ + template void + full_p_adaptivity( + const hp::DoFHandler &); + + template void + p_adaptivity_from_flags( + const hp::DoFHandler &, + const std::vector &); + + template void + force_p_over_h( + const hp::DoFHandler &); + + template void + choose_p_over_h( + const hp::DoFHandler &); + \} + \} +#endif + } + +for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS; + deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace hp + \{ + namespace Refinement + \{ + template void + p_adaptivity_from_threshold( + const hp::DoFHandler &, + const Vector &, + const double, + const double); + + template void + p_adaptivity_from_regularity( + const hp::DoFHandler &, + const Vector &); + + template void + p_adaptivity_from_prediction( + const hp::DoFHandler &, + const Vector &, + const Vector &); + \} + \} +#endif + } diff --git a/tests/hp/hp_coarsening.cc b/tests/hp/hp_coarsening.cc new file mode 100644 index 0000000000..1f13187e9d --- /dev/null +++ b/tests/hp/hp_coarsening.cc @@ -0,0 +1,225 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// validate hp decision algorithms on grid coarsening +// that depend on the composition of h and p adaptivity flags + + +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + + + +template +void +validate(const Triangulation &tria, const hp::DoFHandler &dh) +{ + deallog << "ncells: " << tria.n_global_active_cells() << " fe_indices:"; + for (const auto &cell : dh.active_cell_iterators()) + deallog << " " << cell->active_fe_index(); + deallog << std::endl; +} + + + +template +void +setup(Triangulation & tria, + hp::DoFHandler & dh, + const hp::FECollection &fes) +{ + // Initialize triangulation and dofhandler. + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + dh.initialize(tria, fes); + + // Set h and p flags on all cells. + for (const auto &cell : dh.active_cell_iterators()) + { + cell->set_coarsen_flag(); + cell->set_future_fe_index(1); + } +} + + + +template +void +test() +{ + hp::FECollection fes; + for (unsigned int d = 1; d <= 2; ++d) + fes.push_back(FE_Q(d)); + + deallog << "starting situation: "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + validate(tria, dh); + } + + deallog << "full h&p flags" << std::endl; + { + deallog << " default behaviour: "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " force p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + hp::Refinement::force_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " choose p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + hp::Refinement::choose_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + } + + + deallog << "full p flags" << std::endl; + { + deallog << " default behaviour: "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_coarsen_flag(); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " force p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_coarsen_flag(); + hp::Refinement::force_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " choose p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_coarsen_flag(); + hp::Refinement::choose_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + } + + + deallog << "full h flags" << std::endl; + { + deallog << " default behaviour: "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_future_fe_index(); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " force p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_future_fe_index(); + hp::Refinement::force_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + + deallog << " choose p over h : "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + dh.begin_active()->clear_future_fe_index(); + hp::Refinement::choose_p_over_h(dh); + tria.execute_coarsening_and_refinement(); + + validate(tria, dh); + } + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("1d"); + test<1>(); + deallog.pop(); + deallog.push("2d"); + test<2>(); + deallog.pop(); + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/hp/hp_coarsening.output b/tests/hp/hp_coarsening.output new file mode 100644 index 0000000000..1425a6504b --- /dev/null +++ b/tests/hp/hp_coarsening.output @@ -0,0 +1,43 @@ + +DEAL:1d::starting situation: ncells: 2 fe_indices: 0 0 +DEAL:1d::full h&p flags +DEAL:1d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:1d:: force p over h : ncells: 2 fe_indices: 1 1 +DEAL:1d:: choose p over h : ncells: 2 fe_indices: 1 1 +DEAL:1d::full p flags +DEAL:1d:: default behaviour: ncells: 2 fe_indices: 1 1 +DEAL:1d:: force p over h : ncells: 2 fe_indices: 1 1 +DEAL:1d:: choose p over h : ncells: 2 fe_indices: 1 1 +DEAL:1d::full h flags +DEAL:1d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:1d:: force p over h : ncells: 2 fe_indices: 0 1 +DEAL:1d:: choose p over h : ncells: 1 fe_indices: 0 +DEAL:1d::OK +DEAL:2d::starting situation: ncells: 4 fe_indices: 0 0 0 0 +DEAL:2d::full h&p flags +DEAL:2d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:2d:: force p over h : ncells: 4 fe_indices: 1 1 1 1 +DEAL:2d:: choose p over h : ncells: 4 fe_indices: 1 1 1 1 +DEAL:2d::full p flags +DEAL:2d:: default behaviour: ncells: 4 fe_indices: 1 1 1 1 +DEAL:2d:: force p over h : ncells: 4 fe_indices: 1 1 1 1 +DEAL:2d:: choose p over h : ncells: 4 fe_indices: 1 1 1 1 +DEAL:2d::full h flags +DEAL:2d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:2d:: force p over h : ncells: 4 fe_indices: 0 1 1 1 +DEAL:2d:: choose p over h : ncells: 1 fe_indices: 0 +DEAL:2d::OK +DEAL:3d::starting situation: ncells: 8 fe_indices: 0 0 0 0 0 0 0 0 +DEAL:3d::full h&p flags +DEAL:3d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:3d:: force p over h : ncells: 8 fe_indices: 1 1 1 1 1 1 1 1 +DEAL:3d:: choose p over h : ncells: 8 fe_indices: 1 1 1 1 1 1 1 1 +DEAL:3d::full p flags +DEAL:3d:: default behaviour: ncells: 8 fe_indices: 1 1 1 1 1 1 1 1 +DEAL:3d:: force p over h : ncells: 8 fe_indices: 1 1 1 1 1 1 1 1 +DEAL:3d:: choose p over h : ncells: 8 fe_indices: 1 1 1 1 1 1 1 1 +DEAL:3d::full h flags +DEAL:3d:: default behaviour: ncells: 1 fe_indices: 1 +DEAL:3d:: force p over h : ncells: 8 fe_indices: 0 1 1 1 1 1 1 1 +DEAL:3d:: choose p over h : ncells: 1 fe_indices: 0 +DEAL:3d::OK diff --git a/tests/hp/p_adaptivity_flags.cc b/tests/hp/p_adaptivity_flags.cc new file mode 100644 index 0000000000..67a591b4af --- /dev/null +++ b/tests/hp/p_adaptivity_flags.cc @@ -0,0 +1,233 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// validate algorithms that will flag cells for p adaptivity + + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + + +template +void +validate(const Triangulation &tria, const hp::DoFHandler &dh) +{ + deallog << " fe_indices:"; + for (const auto &cell : dh.active_cell_iterators()) + deallog << " " << cell->future_fe_index(); + deallog << std::endl; +} + + + +template +void +setup(Triangulation & tria, + hp::DoFHandler & dh, + const hp::FECollection &fes) +{ + // Initialize triangulation and dofhandler. + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + dh.initialize(tria, fes); + + // Set all active fe indices to 1. + // Flag first half of cells for refinement, and the other half for coarsening. + typename hp::DoFHandler::cell_iterator cell = dh.begin(1), + endc = dh.end(1); + for (unsigned int counter = 0; cell != endc; ++counter, ++cell) + { + Assert(!cell->active(), ExcInternalError()); + for (unsigned int child_index = 0; child_index < cell->n_children(); + ++child_index) + { + const auto &child = cell->child(child_index); + Assert(child->active(), ExcInternalError()); + + child->set_active_fe_index(1); + + if (counter < 0.5 * GeometryInfo::max_children_per_cell) + child->set_refine_flag(); + else + child->set_coarsen_flag(); + } + } +} + + + +template +void +test() +{ + hp::FECollection fes; + for (unsigned int d = 1; d <= 3; ++d) + fes.push_back(FE_Q(d)); + + deallog << "starting situation: "; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + deallog << "ncells: " << tria.n_active_cells() << std::endl; + validate(tria, dh); + } + + deallog << "full p adaptivity" << std::endl; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + hp::Refinement::full_p_adaptivity(dh); + + validate(tria, dh); + } + + // In the following cases, we flag the first half of all cells to be refined + // and the last half of all cells to be coarsened for p adapativity. + // 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; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + unsigned int n_active = tria.n_active_cells(); + std::vector p_flags(n_active, false); + + std::fill(p_flags.begin(), p_flags.begin() + .25 * n_active, true); + std::fill(p_flags.end() - .25 * n_active, p_flags.end(), true); + + hp::Refinement::p_adaptivity_from_flags(dh, p_flags); + + validate(tria, dh); + } + + deallog << "p adaptivity from 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); + for (unsigned int i = 0; i < n_active; ++i) + { + if (i < .25 * n_active) + smoothness_indicators[i] = 2.; + else if (i < .75 * n_active) + smoothness_indicators[i] = 1.; + else + smoothness_indicators[i] = 0.; + } + hp::Refinement::p_adaptivity_from_threshold(dh, smoothness_indicators); + + validate(tria, dh); + } + + deallog << "p adaptivity from regularity" << std::endl; + { + Triangulation tria; + hp::DoFHandler dh; + setup(tria, dh, fes); + + unsigned int n_active = tria.n_active_cells(); + Vector sobolev_indices(n_active); + for (unsigned int i = 0; i < n_active; ++i) + { + if (i < .25 * n_active) + sobolev_indices[i] = 3. + 1e-4; + else if (i < .75 * n_active) + sobolev_indices[i] = 2.; + else + sobolev_indices[i] = 1. - 1e-4; + } + hp::Refinement::p_adaptivity_from_regularity(dh, sobolev_indices); + + validate(tria, dh); + } + + deallog << "p adaptivity from prediction" << 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); + for (unsigned int i = 0; i < n_active; ++i) + { + if (i < .25 * n_active) + { + predicted_errors[i] = 1. + 1e-4; + error_estimates[i] = 1.; + } + else if (i < .75 * n_active) + { + predicted_errors[i] = 1.; + error_estimates[i] = 1.; + } + else + { + predicted_errors[i] = 1. + 1e-4; + error_estimates[i] = 1.; + } + } + hp::Refinement::p_adaptivity_from_prediction(dh, + error_estimates, + predicted_errors); + + validate(tria, dh); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("1d"); + test<1>(); + deallog.pop(); + deallog.push("2d"); + test<2>(); + deallog.pop(); + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/hp/p_adaptivity_flags.output b/tests/hp/p_adaptivity_flags.output new file mode 100644 index 0000000000..6061a757bb --- /dev/null +++ b/tests/hp/p_adaptivity_flags.output @@ -0,0 +1,40 @@ + +DEAL:1d::starting situation: ncells: 4 +DEAL:1d:: fe_indices: 1 1 1 1 +DEAL:1d::full p adaptivity +DEAL:1d:: fe_indices: 2 2 0 0 +DEAL:1d::p adaptivity from flags +DEAL:1d:: fe_indices: 2 1 1 0 +DEAL:1d::p adaptivity from threshold +DEAL:1d:: fe_indices: 2 1 1 0 +DEAL:1d::p adaptivity from regularity +DEAL:1d:: fe_indices: 2 1 1 0 +DEAL:1d::p adaptivity from prediction +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:: 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:: 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:: 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:: 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:: 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:: 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:: 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:: 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:: 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:: 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