+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2019 - 2020 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_distributed_error_predictor_h
-#define dealii_distributed_error_predictor_h
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/smartpointer.h>
-
-#include <boost/range/iterator_range.hpp>
-
-#include <vector>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-
-// forward declarations
-#ifndef DOXYGEN
-template <int dim, int spacedim>
-class DoFHandler;
-
-template <typename Number>
-class Vector;
-
-namespace parallel
-{
- namespace distributed
- {
- template <int dim, int spacedim>
- class Triangulation;
- }
-} // namespace parallel
-#endif
-
-
-namespace parallel
-{
- namespace distributed
- {
- /**
- * Predict how current error indicators will change after refinement and
- * coarsening were to happen on the provided DoFHandler in context of a
- * parallel::distributed::Triangulation.
- *
- * This algorithm follows the same logic as the error prediction algorithm
- * presented by @cite melenk2001hp and has been implemented in deal.II via
- * the function hp::Refinement::predict_error(). See the documentation of
- * this particular function for details on the algorithm and all featured
- * formulas.
- *
- * In theory, the combination of hp::Refinement::predict_error() to predict
- * how the error will change with parallel::distributed::CellDataTransfer
- * to transfer data across meshes forms the intended way to apply the
- * prediction algorithm. However, p4est determines the details of grid
- * refinement in the parallel distributed case; consequently, it yields more
- * reliable and trustworthy results when we determine the predicted errors
- * during the adaptation process. This is achieved with this class, similar
- * to the parallel::distributed::SolutionTransfer and
- * parallel::distributed::CellDataTransfer classes.
- *
- * Before adaptation happens, one or multiple vectors containing error
- * estimates for every cell have to be registered with this class using the
- * prepare_for_coarsening_and_refinement() functions. Performing adaptation
- * on the triangulation initiates the calculation of the error prediction
- * based on how p4est performed grid adaptation. After the refinement
- * process, the predicted errors need to be unpacked on the updated mesh via
- * one of the unpack() functions.
- *
- * The following code snippet demonstrates how to impose hp-adaptivity based
- * on refinement history in a parallel distributed application:
- * @code
- * // [initialization...]
- * Vector<float> predicted_error_per_cell(triangulation.n_active_cells());
- * for(unsigned int i = 0; i < triangulation.n_active_cells(); ++i)
- * predicted_error_per_cell[i] = std::numeric_limits<float>::infinity();
- *
- * // [during each refinement step...]
- * // set h-adaptivity flags
- * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
- * KellyErrorEstimator<dim>::estimate(...);
- * parallel::distributed::GridRefinement::
- * refine_and_coarsen_fixed_{number|fraction}(...);
- *
- * // set p-adaptivity flags
- * hp::Refinement::p_adaptivity_from_reference(
- * hp_dof_handler,
- * estimated_error_per_cell,
- * predicted_error_per_cell
- * std::less<float>(),
- * std::less<float>());
- * hp::Refinement::{choose|force}_p_over_h(hp_dof_handler);
- *
- * // perform adaptation and predict error for the subsequent adaptation
- * parallel::distributed::ErrorPredictor<dim> predictor(hp_dof_handler);
- * predictor.prepare_for_coarsening_and_refinement(estimated_error_per_cell);
- * triangulation.execute_coarsening_and_refinement();
- *
- * // unpack predicted errors
- * predicted_error_per_cell.reinit(triangulation.n_active_cells());
- * predictor.unpack(predicted_error_per_cell);
- * @endcode
- *
- * For serialization of predicted errors, use the
- * parallel::distributed::CellDataTransfer class.
- *
- * @ingroup distributed
- */
- template <int dim, int spacedim = dim>
- class ErrorPredictor
- {
- public:
- /**
- * Constructor.
- *
- * @param[in] dof The DoFHandler on which all operations will
- * happen. At the time when this constructor is called, the
- * DoFHandler still points to the triangulation before the
- * refinement in question happens.
- */
- ErrorPredictor(const DoFHandler<dim, spacedim> &dof);
-
- /**
- * Prepare the current object for coarsening and refinement.
- * @p all_in includes all vectors that are to be interpolated onto the
- * new (refined and/or coarsened) grid.
- *
- * As specified in the error prediction algorithm in
- * hp::Refinement::predict_error(), the control parameters @p gamma_p,
- * @p gamma_h, @p gamma_n can be used to influence the predicted errors in
- * case of p-adaptation, h-adaptation, and no adapation, respectively.
- * Their default values comply to the studies of @cite melenk2001hp.
- */
- void
- prepare_for_coarsening_and_refinement(
- const std::vector<const Vector<float> *> &all_in,
- const double gamma_p = std::sqrt(0.4),
- const double gamma_h = 2.,
- const double gamma_n = 1.);
-
- /**
- * Same as the previous function but for only one discrete function to be
- * interpolated.
- */
- void
- prepare_for_coarsening_and_refinement(
- const Vector<float> &in,
- const double gamma_p = std::sqrt(0.4),
- const double gamma_h = 2.,
- const double gamma_n = 1.);
-
- /**
- * Interpolate the data previously stored in this object before the mesh
- * was refined or coarsened onto the current set of cells. Do so for
- * each of the vectors provided to
- * prepare_for_coarsening_and_refinement() and write the result into the
- * given set of vectors.
- */
- void
- unpack(std::vector<Vector<float> *> &all_out);
-
- /**
- * Same as the previous function. It interpolates only one function. It
- * assumes the vectors having the right sizes (i.e.
- * <tt>in.size()==n_cells_old</tt>, <tt>out.size()==n_cells_refined</tt>)
- *
- * Multiple calling of this function is NOT allowed. Interpolating
- * several functions can be performed in one step by using the above
- * function.
- */
- void
- unpack(Vector<float> &out);
-
- private:
- /**
- * Pointer to the degree of freedom handler to work with.
- */
- SmartPointer<const DoFHandler<dim, spacedim>,
- ErrorPredictor<dim, spacedim>>
- dof_handler;
-
- /**
- * A vector that stores pointers to all the vectors we are supposed to
- * copy over from the old to the new mesh.
- */
- std::vector<const Vector<float> *> error_indicators;
-
- /**
- * The handle that the Triangulation has assigned to this object
- * with which we can access our memory offset and our pack function.
- */
- unsigned int handle;
-
- /**
- * Control parameters for the prediction algorithm.
- */
- double gamma_p, gamma_h, gamma_n;
-
- /**
- * A callback function used to pack the data on the current mesh into
- * objects that can later be retrieved after refinement, coarsening and
- * repartitioning.
- */
- std::vector<char>
- pack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const typename Triangulation<dim, spacedim>::CellStatus status);
-
- /**
- * A callback function used to unpack the data on the current mesh that
- * has been packed up previously on the mesh before refinement,
- * coarsening and repartitioning.
- */
- void
- unpack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const typename Triangulation<dim, spacedim>::CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- & data_range,
- std::vector<Vector<float> *> &all_out);
-
-
- /**
- * Registers the pack_callback() function to the
- * parallel::distributed::Triangulation that has been assigned to the
- * DoFHandler class member and stores the returning handle.
- */
- void
- register_data_attach();
- };
- } // namespace distributed
-} // namespace parallel
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2019 - 2020 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/base/config.h>
-
-#ifdef DEAL_II_WITH_P4EST
-
-# include <deal.II/distributed/error_predictor.h>
-# include <deal.II/distributed/tria.h>
-
-# include <deal.II/dofs/dof_accessor.h>
-# include <deal.II/dofs/dof_handler.h>
-# include <deal.II/dofs/dof_tools.h>
-
-# include <deal.II/grid/tria_accessor.h>
-# include <deal.II/grid/tria_iterator.h>
-
-# include <deal.II/lac/block_vector.h>
-# include <deal.II/lac/la_parallel_block_vector.h>
-# include <deal.II/lac/la_parallel_vector.h>
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# include <deal.II/lac/trilinos_parallel_block_vector.h>
-# include <deal.II/lac/trilinos_vector.h>
-# include <deal.II/lac/vector.h>
-
-# include <functional>
-# include <numeric>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-
-namespace parallel
-{
- namespace distributed
- {
- template <int dim, int spacedim>
- ErrorPredictor<dim, spacedim>::ErrorPredictor(
- const DoFHandler<dim, spacedim> &dof)
- : dof_handler(&dof, typeid(*this).name())
- , handle(numbers::invalid_unsigned_int)
- , gamma_p(0.)
- , gamma_h(0.)
- , gamma_n(0.)
- {
- Assert((dynamic_cast<
- const parallel::distributed::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation()) != nullptr),
- ExcMessage("parallel::distributed::ErrorPredictor requires a "
- "parallel::distributed::Triangulation object."));
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::prepare_for_coarsening_and_refinement(
- const std::vector<const Vector<float> *> &all_in,
- const double gamma_p,
- const double gamma_h,
- const double gamma_n)
- {
- error_indicators = all_in;
- this->gamma_p = gamma_p;
- this->gamma_h = gamma_h;
- this->gamma_n = gamma_n;
- register_data_attach();
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::register_data_attach()
- {
- // TODO: casting away constness is bad
- parallel::distributed::Triangulation<dim, spacedim> *tria =
- (dynamic_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
- const_cast<dealii::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation())));
- Assert(tria != nullptr, ExcInternalError());
-
- handle = tria->register_data_attach(
- [this](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const typename Triangulation<dim, spacedim>::CellStatus status) {
- return this->pack_callback(cell, status);
- },
- /*returns_variable_size_data=*/false);
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::prepare_for_coarsening_and_refinement(
- const Vector<float> &in,
- const double gamma_p,
- const double gamma_h,
- const double gamma_n)
- {
- std::vector<const Vector<float> *> all_in(1, &in);
- prepare_for_coarsening_and_refinement(all_in, gamma_p, gamma_h, gamma_n);
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::unpack(std::vector<Vector<float> *> &all_out)
- {
- Assert(error_indicators.size() == all_out.size(),
- ExcDimensionMismatch(error_indicators.size(), all_out.size()));
-
- // TODO: casting away constness is bad
- parallel::distributed::Triangulation<dim, spacedim> *tria =
- (dynamic_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
- const_cast<dealii::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation())));
- Assert(tria != nullptr, ExcInternalError());
-
- tria->notify_ready_to_unpack(
- handle,
- [this, &all_out](
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const typename Triangulation<dim, spacedim>::CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- &data_range) {
- this->unpack_callback(cell, status, data_range, all_out);
- });
-
- for (const auto &out : all_out)
- out->compress(::dealii::VectorOperation::insert);
-
- error_indicators.clear();
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::unpack(Vector<float> &out)
- {
- std::vector<Vector<float> *> all_out(1, &out);
- unpack(all_out);
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<char>
- ErrorPredictor<dim, spacedim>::pack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
- const typename Triangulation<dim, spacedim>::CellStatus status)
- {
- typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
- dof_handler);
-
- // create buffer for each individual input vector
- std::vector<float> predicted_errors(error_indicators.size());
-
- auto predicted_error_it = predicted_errors.begin();
- auto estimated_error_it = error_indicators.cbegin();
- for (; estimated_error_it != error_indicators.cend();
- ++predicted_error_it, ++estimated_error_it)
- switch (status)
- {
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_PERSIST:
- {
- if (cell->future_fe_index_set())
- {
- const int degree_difference =
- dof_handler->get_fe_collection()[cell->future_fe_index()]
- .degree -
- cell->get_fe().degree;
-
- *predicted_error_it =
- (**estimated_error_it)[cell->active_cell_index()] *
- std::pow(gamma_p, degree_difference);
- }
- else
- {
- *predicted_error_it =
- (**estimated_error_it)[cell->active_cell_index()] *
- gamma_n;
- }
- break;
- }
-
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_REFINE:
- {
- // Determine the exponent by the finite element degree on the
- // adapted mesh.
- const unsigned int future_fe_degree =
- dof_handler->get_fe_collection()[cell->future_fe_index()]
- .degree;
-
- *predicted_error_it =
- (**estimated_error_it)[cell->active_cell_index()] *
- (gamma_h * std::pow(.5, future_fe_degree));
-
- // If the future FE index differs from the active one, also take
- // into account p-adaptation.
- if (cell->future_fe_index_set())
- *predicted_error_it *=
- std::pow(gamma_p,
- static_cast<int>(future_fe_degree -
- cell->get_fe().degree));
-
- break;
- }
-
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_COARSEN:
- {
- // First figure out which finite element will be assigned to the
- // parent cell after h-adaptation analogously to
- // dealii::internal::hp::DoFHandlerImplementation::Implementation::
- // collect_fe_indices_on_cells_to_be_refined()
-# ifdef DEBUG
- for (const auto &child : cell->child_iterators())
- Assert(child->is_active() && child->coarsen_flag_set(),
- typename dealii::Triangulation<
- dim>::ExcInconsistentCoarseningFlags());
-# endif
-
- const unsigned int future_fe_index =
- dealii::internal::hp::DoFHandlerImplementation::
- dominated_future_fe_on_children<dim, spacedim>(cell);
-
- const unsigned int future_fe_degree =
- dof_handler->get_fe_collection()[future_fe_index].degree;
-
- // Then determine the actually contirbution to the predicted
- // error of every single cells that is about to be coarsened.
- float sqrsum_of_predicted_errors = 0.;
- float predicted_error = 0.;
- int degree_difference = 0;
- for (const auto &child : cell->child_iterators())
- {
- predicted_error =
- (**estimated_error_it)[child->active_cell_index()] /
- (gamma_h * std::pow(.5, future_fe_degree));
-
- degree_difference =
- future_fe_degree - child->get_fe().degree;
-
- if (degree_difference != 0)
- predicted_error *= std::pow(gamma_p, degree_difference);
-
- sqrsum_of_predicted_errors +=
- predicted_error * predicted_error;
- }
- *predicted_error_it = std::sqrt(sqrsum_of_predicted_errors);
-
- break;
- }
-
- default:
- Assert(false, ExcInternalError());
- break;
- }
-
- // We don't have to pack the whole container if there is just one entry.
- if (error_indicators.size() == 1)
- return Utilities::pack(predicted_errors[0],
- /*allow_compression=*/false);
- else
- return Utilities::pack(predicted_errors, /*allow_compression=*/false);
- }
-
-
-
- template <int dim, int spacedim>
- void
- ErrorPredictor<dim, spacedim>::unpack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const typename Triangulation<dim, spacedim>::CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- & data_range,
- std::vector<Vector<float> *> &all_out)
- {
- std::vector<float> predicted_errors;
-
- if (all_out.size() == 1)
- predicted_errors.push_back(
- Utilities::unpack<float>(data_range.begin(),
- data_range.end(),
- /*allow_compression=*/false));
- else
- predicted_errors =
- Utilities::unpack<std::vector<float>>(data_range.begin(),
- data_range.end(),
- /*allow_compression=*/false);
-
- Assert(predicted_errors.size() == all_out.size(), ExcInternalError());
-
- auto it_input = predicted_errors.cbegin();
- auto it_output = all_out.begin();
- for (; it_input != predicted_errors.cend(); ++it_input, ++it_output)
- switch (status)
- {
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_PERSIST:
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_COARSEN:
- (**it_output)[cell->active_cell_index()] = *it_input;
- break;
-
-
- case parallel::distributed::Triangulation<dim,
- spacedim>::CELL_REFINE:
- for (const auto &child : cell->child_iterators())
- (**it_output)[child->active_cell_index()] =
- (*it_input) / std::sqrt(cell->n_children());
- break;
-
- default:
- Assert(false, ExcInternalError());
- break;
- }
- }
- } // namespace distributed
-} // namespace parallel
-
-
-// explicit instantiations
-# include "error_predictor.inst"
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif /* DEAL_II_WITH_P4EST */
projects / dealii.git / commitdiff