]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduced hp::Refinement namespace.
authorMarc Fehling <marc.fehling@gmx.net>
Fri, 7 Jun 2019 10:01:46 +0000 (12:01 +0200)
committermarcfehling <marc.fehling@gmx.net>
Thu, 4 Jul 2019 08:39:42 +0000 (10:39 +0200)
doc/news/changes/minor/20190607Fehling [new file with mode: 0644]
include/deal.II/hp/refinement.h [new file with mode: 0644]
source/hp/CMakeLists.txt
source/hp/refinement.cc [new file with mode: 0644]
source/hp/refinement.inst.in [new file with mode: 0644]
tests/hp/hp_coarsening.cc [new file with mode: 0644]
tests/hp/hp_coarsening.output [new file with mode: 0644]
tests/hp/p_adaptivity_flags.cc [new file with mode: 0644]
tests/hp/p_adaptivity_flags.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190607Fehling b/doc/news/changes/minor/20190607Fehling
new file mode 100644 (file)
index 0000000..24635b8
--- /dev/null
@@ -0,0 +1,3 @@
+New: Namespace hp::Refinement offering decision tools for p adaptivity.
+<br>
+(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 (file)
index 0000000..7101ab1
--- /dev/null
@@ -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/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
index 4da18191407f219107057738642ced384105f3a9..d1fa66f01862bfecac96be8564352ee480187f6c 100644 (file)
@@ -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 (file)
index 0000000..a8ae27c
--- /dev/null
@@ -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 <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
diff --git a/source/hp/refinement.inst.in b/source/hp/refinement.inst.in
new file mode 100644 (file)
index 0000000..480c873
--- /dev/null
@@ -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<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
+  }
diff --git a/tests/hp/hp_coarsening.cc b/tests/hp/hp_coarsening.cc
new file mode 100644 (file)
index 0000000..1f13187
--- /dev/null
@@ -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 <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();
+}
diff --git a/tests/hp/hp_coarsening.output b/tests/hp/hp_coarsening.output
new file mode 100644 (file)
index 0000000..1425a65
--- /dev/null
@@ -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 (file)
index 0000000..67a591b
--- /dev/null
@@ -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 <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();
+}
diff --git a/tests/hp/p_adaptivity_flags.output b/tests/hp/p_adaptivity_flags.output
new file mode 100644 (file)
index 0000000..6061a75
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.