* @}
*/
+ /**
+ * @name Error prediction
+ * @{
+ */
+
+ /**
+ * Predict how the current @p error_indicators will adapt after refinement
+ * and coarsening has happened on the provided @p dof_handler, and write its
+ * results to @p predicted_errors. Each entry of @p error_indicators and
+ * @p predicted_errors corresponds to an active cell on the underlying
+ * Triangulation, thus each container has to be of size
+ * Triangulation::n_active_cells().
+ *
+ * For h adaptation, we expect the local error $\eta_K$ on cell $K$ to be
+ * proportional to $(h_K)^{p_K}$ in the energy norm, where $h_K$ denotes the
+ * cell diameter and $p_K$ the polynomial degree of the currently assigned
+ * finite element. Here, we assume that the finite element will not change
+ * in the adaptation process so that $p_K = \text{const}$. However during
+ * coarsening, the finite elements on siblings may be different, and their
+ * parent cell will be assigned to their least dominating finite element
+ * that belongs to its most general child. Thus, we will always interpolate
+ * on an enclosing finite element space. Additionaly assuming that the
+ * finite elements on the cells to be coarsened are sufficient to represent
+ * the solution correct (e.g. at least quadratic basis functions for a
+ * quadratic solution), we are confident to say that the error will not
+ * change by sole interpolation on the larger finite element space.
+ *
+ * Further, we expect that the local error will be divided equally on
+ * all $2^{dim}$ children during refinement, whereas local errors on
+ * siblings will be summed up on the parent cell in case of coarsening. When
+ * transferring the predicted error to the coarsened mesh, make sure to
+ * configure your CellDataTransfer object with
+ * GridTools::CoarseningStrategies::sum() as a coarsening strategy.
+ *
+ * For p adaptation, the local error is expected to converge exponentially
+ * with the polynomial degree of the assigned finite element. Each increase
+ * or decrease of the degree will thus change its value by a user-defined
+ * control parameter @p gamma_p. The assumption of exponential convergence
+ * is only valid if both h and p adaptive methods are combined. An exception
+ * is thrown if a cell is flagged for both h and p adaptation at once.
+ *
+ * The prediction algorithm is formulated as follows with control parameters
+ * @p gamma_p, @p gamma_h and @p gamma_n that may be used to influence
+ * prediction for each adaptation type individually.
+ * <table>
+ * <tr><th>Adaptation type <th colspan="2">Prediction formula
+ * <tr><td>no adaptation
+ * <td>$\eta_{K,\text{pred}} = \eta_{K} \, \gamma_\text{n}$
+ * <td>$\gamma_\text{n} \in (0,\infty)$
+ * <tr><td>p adaptation
+ * <td>$\eta_{K,\text{pred}} = \eta_{K} \,
+ * \gamma_\text{p}^{(p_{K,\text{future}} - p_K)}$
+ * <td>$\gamma_\text{p} \in (0,1)$
+ * <tr><td>h refinement
+ * <td>$\eta_{K_c,\text{pred}} = \eta_{K} \,
+ * \gamma_\text{h} \, 0.5^{p_K} \, 0.5^{\text{dim}}
+ * \quad \forall K_c \text{ children of } K$
+ * <td rowspan="2">$\gamma_\text{h} \in (0,\infty)$
+ * <tr><td>h coarsening
+ * <td>$\eta_{K,\text{pred}} = \sum\limits_{K_c} \eta_{K_c} /
+ * (\gamma_\text{h} \, 0.5^{p_{K_c}})
+ * \quad \forall K_c \text{ children of } K$
+ * </table>
+ *
+ * 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 This feature is currently only implemented for isotropic refinement.
+ *
+ * @note We want to predict the error by how adaptation will actually happen.
+ * Thus, this function needs to be called after
+ * Triangulation::prepare_for_coarsening_and_refinement().
+ */
+ template <int dim, typename Number, int spacedim>
+ void
+ predict_error(const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & error_indicators,
+ Vector<Number> & predicted_errors,
+ const double gamma_p = std::sqrt(0.1),
+ const double gamma_h = 1.,
+ const double gamma_n = 1.);
+
+ /**
+ * @}
+ */
+
/**
* @name Decide between h and p adaptivity
* @{
+ /**
+ * Error prediction
+ */
+ template <int dim, typename Number, int spacedim>
+ void
+ predict_error(const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const Vector<Number> & error_indicators,
+ Vector<Number> & predicted_errors,
+ const double gamma_p,
+ const double gamma_h,
+ const double gamma_n)
+ {
+ AssertDimension(dof_handler.get_triangulation().n_active_cells(),
+ error_indicators.size());
+ AssertDimension(dof_handler.get_triangulation().n_active_cells(),
+ predicted_errors.size());
+ Assert(0 < gamma_p && gamma_p < 1,
+ dealii::GridRefinement::ExcInvalidParameterValue());
+ Assert(0 < gamma_h, dealii::GridRefinement::ExcInvalidParameterValue());
+ Assert(0 < gamma_n, dealii::GridRefinement::ExcInvalidParameterValue());
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ const unsigned int active_cell_index = cell->active_cell_index();
+
+ if (cell->future_fe_index_set()) // p adaptation
+ {
+ Assert(!cell->refine_flag_set() && !cell->coarsen_flag_set(),
+ ExcMessage("Cell has to be either flagged for h or p "
+ "adaptation, and not for both!"));
+
+ const int degree_difference =
+ dof_handler.get_fe_collection()[cell->future_fe_index()]
+ .degree -
+ cell->get_fe().degree;
+
+ predicted_errors[active_cell_index] =
+ error_indicators[active_cell_index] *
+ std::pow(gamma_p, degree_difference);
+ }
+ else if (cell->refine_flag_set()) // h refinement
+ {
+ Assert(
+ cell->refine_flag_set() ==
+ RefinementCase<dim>::isotropic_refinement,
+ ExcMessage(
+ "Error prediction is only valid for isotropic refinement!"));
+
+ predicted_errors[active_cell_index] =
+ error_indicators[active_cell_index] *
+ (gamma_h * std::pow(.5, dim + cell->get_fe().degree));
+ }
+ else if (cell->coarsen_flag_set()) // h coarsening
+ {
+ predicted_errors[active_cell_index] =
+ error_indicators[active_cell_index] /
+ (gamma_h * std::pow(.5, cell->get_fe().degree));
+ }
+ else // no changes
+ {
+ predicted_errors[active_cell_index] =
+ error_indicators[active_cell_index] * gamma_n;
+ }
+ }
+ }
+
+
+
/**
* Decide between h and p adaptivity
*/
const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const Vector<S> &,
const Vector<S> &);
+
+ template void
+ predict_error<deal_II_dimension, S, deal_II_space_dimension>(
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const Vector<S> &,
+ Vector<S> &,
+ const double,
+ const double,
+ const double);
\}
\}
#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 error prediction algorithm for hp adaptive methods
+
+
+#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
+test()
+{
+ const unsigned int n_cells = 4;
+
+ // ----- setup -----
+ Triangulation<dim> tria;
+ std::vector<unsigned int> rep(dim, 1);
+ rep[0] = n_cells;
+ Point<dim> p1, p2;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ p1[d] = 0;
+ p2[d] = (d == 0) ? n_cells : 1;
+ }
+ GridGenerator::subdivided_hyper_rectangle(tria, rep, p1, p2);
+
+ hp::FECollection<dim> fes;
+ for (unsigned int d = 1; d <= 3; ++d)
+ fes.push_back(FE_Q<dim>(d));
+
+ hp::DoFHandler<dim> dh(tria);
+ dh.set_fe(fes);
+ for (const auto &cell : dh.active_cell_iterators())
+ {
+ if (cell->id().to_string() == "0_0:")
+ cell->set_refine_flag();
+ else if (cell->id().to_string() == "1_0:")
+ cell->set_coarsen_flag();
+ else if (cell->id().to_string() == "2_0:")
+ cell->set_future_fe_index(2);
+ }
+
+ // ----- predict -----
+ Vector<float> error_indicators, predicted_errors;
+ error_indicators.reinit(tria.n_active_cells());
+ predicted_errors.reinit(tria.n_active_cells());
+ for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
+ error_indicators[i] = 10;
+
+ hp::Refinement::predict_error(dh,
+ error_indicators,
+ predicted_errors,
+ /*gamma_p=*/0.5,
+ /*gamma_h=*/1.,
+ /*gamma_n=*/1.);
+
+ // ----- verify ------
+ deallog << "ncells:" << tria.n_active_cells() << std::endl;
+ for (const auto &cell : dh.active_cell_iterators())
+ {
+ deallog << " cell:" << cell->id().to_string()
+ << " fe_deg:" << cell->get_fe().degree
+ << " error:" << error_indicators[cell->active_cell_index()]
+ << " predict:" << predicted_errors[cell->active_cell_index()];
+
+ if (cell->refine_flag_set())
+ deallog << " refining";
+ else if (cell->coarsen_flag_set())
+ deallog << " coarsening";
+ else if (cell->future_fe_index_set())
+ deallog << " future_fe_deg:"
+ << dh.get_fe_collection()[cell->future_fe_index()].degree;
+
+ deallog << std::endl;
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+ deallog << std::setprecision(3);
+
+ 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::ncells:4
+DEAL:1d:: cell:0_0: fe_deg:1 error:10.0 predict:2.50 refining
+DEAL:1d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening
+DEAL:1d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3
+DEAL:1d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0
+DEAL:1d::OK
+DEAL:2d::ncells:4
+DEAL:2d:: cell:0_0: fe_deg:1 error:10.0 predict:1.25 refining
+DEAL:2d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening
+DEAL:2d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3
+DEAL:2d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0
+DEAL:2d::OK
+DEAL:3d::ncells:4
+DEAL:3d:: cell:0_0: fe_deg:1 error:10.0 predict:0.625 refining
+DEAL:3d:: cell:1_0: fe_deg:1 error:10.0 predict:20.0 coarsening
+DEAL:3d:: cell:2_0: fe_deg:1 error:10.0 predict:2.50 future_fe_deg:3
+DEAL:3d:: cell:3_0: fe_deg:1 error:10.0 predict:10.0
+DEAL:3d::OK
projects / dealii.git / commitdiff