--- /dev/null
+New: Namespace hp::Refinement offering decision tools for p adaptivity.
+<br>
+(Marc Fehling, 2019/06/07)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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/base/config.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// forward declarations
+template <typename Number>
+class Vector;
+
+namespace hp
+{
+ template <int dim, int spacedim>
+ 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.
+ *
+ * <h3>Usage</h3>
+ *
+ * To successfully apply hp adaptive methods, we recommend the following
+ * workflow:
+ * <ol>
+ * <li> 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.
+ *
+ * <li> 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.
+ *
+ * <li> 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.
+ *
+ * <li> 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.
+ * </ol>
+ *
+ * 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<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);
+ *
+ * // 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 <int dim, int spacedim>
+ void
+ 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.
+ *
+ * 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 <int dim, int spacedim>
+ void
+ p_adaptivity_from_flags(const hp::DoFHandler<dim, spacedim> &dof_handler,
+ 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]$.
+ *
+ * @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_threshold(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_regularity(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ 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.
+ *
+ * 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_prediction(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & error_indicators,
+ const Vector<Number> & 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 <int dim, int spacedim>
+ void
+ force_p_over_h(const hp::DoFHandler<dim, spacedim> &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:
+ * <ol>
+ * <li> Not all siblings flagged for coarsening: p coarsening<br>
+ * We keep the @p future_fe_indices and clear the coarsen flags
+ * on all siblings.
+ * <li> All siblings flagged for coarsening, but not all for
+ * p adaptation: h coarsening<br>
+ * We keep the coarsen flags and clear all @p future_fe_indices
+ * on all siblings.
+ * <li> All siblings flagged for coarsening and p adaptation: p coarsening<br>
+ * We keep the @p future_fe_indices and clear the coarsen flags
+ * on all siblings.
+ * </ol>
+ *
+ * @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 <int dim, int spacedim>
+ void
+ choose_p_over_h(const hp::DoFHandler<dim, spacedim> &dof_handler);
+
+ /**
+ * @}
+ */
+ } // namespace Refinement
+} // namespace hp
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // dealii_hp_refinement_h
fe_collection.cc
fe_values.cc
mapping_collection.cc
+ refinement.cc
)
SET(_separate_src
fe_collection.inst.in
fe_values.inst.in
mapping_collection.inst.in
+ refinement.inst.in
)
FILE(GLOB _header
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/refinement.h>
+
+#include <deal.II/lac/vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace hp
+{
+ namespace Refinement
+ {
+ /**
+ * Setting p adaptivity flags
+ */
+ template <int dim, int spacedim>
+ void
+ full_p_adaptivity(const hp::DoFHandler<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ p_adaptivity_from_flags(const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const std::vector<bool> & 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_threshold(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_regularity(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & 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 <int dim, typename Number, int spacedim>
+ void
+ p_adaptivity_from_prediction(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & error_indicators,
+ const Vector<Number> & 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 <int dim, int spacedim>
+ void
+ force_p_over_h(const hp::DoFHandler<dim, spacedim> &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 <int dim, int spacedim>
+ void
+ choose_p_over_h(const hp::DoFHandler<dim, spacedim> &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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<deal_II_dimension, deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &);
+
+ template void
+ p_adaptivity_from_flags<deal_II_dimension, deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const std::vector<bool> &);
+
+ template void
+ force_p_over_h<deal_II_dimension, deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &);
+
+ template void
+ choose_p_over_h<deal_II_dimension, deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &);
+ \}
+ \}
+#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<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);
+
+ template void
+ p_adaptivity_from_regularity<deal_II_dimension,
+ S,
+ deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const Vector<S> &);
+
+ template void
+ p_adaptivity_from_prediction<deal_II_dimension,
+ S,
+ deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const Vector<S> &,
+ const Vector<S> &);
+ \}
+ \}
+#endif
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/refinement.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+validate(const Triangulation<dim> &tria, const hp::DoFHandler<dim> &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 <int dim>
+void
+setup(Triangulation<dim> & tria,
+ hp::DoFHandler<dim> & dh,
+ const hp::FECollection<dim> &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 <int dim>
+void
+test()
+{
+ hp::FECollection<dim> fes;
+ for (unsigned int d = 1; d <= 2; ++d)
+ fes.push_back(FE_Q<dim>(d));
+
+ deallog << "starting situation: ";
+ {
+ Triangulation<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+ validate(tria, dh);
+ }
+
+ deallog << "full h&p flags" << std::endl;
+ {
+ deallog << " default behaviour: ";
+ {
+ Triangulation<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ tria.execute_coarsening_and_refinement();
+
+ validate(tria, dh);
+ }
+
+ deallog << " force p over h : ";
+ {
+ Triangulation<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> 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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/geometry_info.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/refinement.h>
+
+#include <deal.II/lac/vector.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+validate(const Triangulation<dim> &tria, const hp::DoFHandler<dim> &dh)
+{
+ deallog << " fe_indices:";
+ for (const auto &cell : dh.active_cell_iterators())
+ deallog << " " << cell->future_fe_index();
+ deallog << std::endl;
+}
+
+
+
+template <int dim>
+void
+setup(Triangulation<dim> & tria,
+ hp::DoFHandler<dim> & dh,
+ const hp::FECollection<dim> &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<dim>::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<dim>::max_children_per_cell)
+ child->set_refine_flag();
+ else
+ child->set_coarsen_flag();
+ }
+ }
+}
+
+
+
+template <int dim>
+void
+test()
+{
+ hp::FECollection<dim> fes;
+ for (unsigned int d = 1; d <= 3; ++d)
+ fes.push_back(FE_Q<dim>(d));
+
+ deallog << "starting situation: ";
+ {
+ Triangulation<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ deallog << "ncells: " << tria.n_active_cells() << std::endl;
+ validate(tria, dh);
+ }
+
+ deallog << "full p adaptivity" << std::endl;
+ {
+ Triangulation<dim> tria;
+ hp::DoFHandler<dim> 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<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ unsigned int n_active = tria.n_active_cells();
+ std::vector<bool> 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<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ unsigned int n_active = tria.n_active_cells();
+ Vector<double> 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<dim> tria;
+ hp::DoFHandler<dim> dh;
+ setup(tria, dh, fes);
+
+ unsigned int n_active = tria.n_active_cells();
+ Vector<double> 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<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);
+ 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();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff