]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Merge SolutionTransfer and p::d::SolutionTransfer 15706/head
authorPasquale Africa <pafrica@sissa.it>
Mon, 10 Jul 2023 18:55:31 +0000 (18:55 +0000)
committerPasquale Africa <pafrica@sissa.it>
Wed, 14 Aug 2024 08:59:11 +0000 (08:59 +0000)
63 files changed:
doc/news/changes/incompatibilities/20240103AfricaBlaisMunch [new file with mode: 0644]
examples/step-15/step-15.cc
examples/step-26/step-26.cc
examples/step-28/step-28.cc
examples/step-31/step-31.cc
examples/step-33/step-33.cc
examples/step-43/step-43.cc
examples/step-57/step-57.cc
examples/step-72/step-72.cc
examples/step-77/step-77.cc
include/deal.II/distributed/solution_transfer.h
include/deal.II/distributed/tria_base.h
include/deal.II/grid/tria.h
include/deal.II/numerics/solution_transfer.h
source/distributed/CMakeLists.txt
source/distributed/fully_distributed_tria.cc
source/distributed/solution_transfer.cc [deleted file]
source/distributed/solution_transfer.inst.in [deleted file]
source/grid/tria.cc
source/numerics/solution_transfer.cc
source/numerics/solution_transfer.inst.in
tests/aniso/solution_transfer.cc
tests/bits/solution_transfer.cc
tests/codim_one/solution_transfer_01.cc
tests/distributed_grids/solution_transfer_01.cc
tests/fullydistributed_grids/solution_transfer_01.cc
tests/grid/save_load_01.cc
tests/grid/solution_transfer_01.cc [new file with mode: 0644]
tests/grid/solution_transfer_01.output [new file with mode: 0644]
tests/grid/solution_transfer_02.cc [new file with mode: 0644]
tests/grid/solution_transfer_02.output [new file with mode: 0644]
tests/hp/do_function_derivatives_01.cc
tests/hp/do_function_hessians_01.cc
tests/hp/do_function_laplacians_01.cc
tests/hp/fe_nothing_17.cc
tests/hp/laplace.h
tests/hp/solution_transfer.cc
tests/hp/solution_transfer_02.cc
tests/hp/solution_transfer_03.cc
tests/hp/solution_transfer_04.cc
tests/hp/solution_transfer_05.cc
tests/hp/solution_transfer_06.cc
tests/hp/solution_transfer_07.cc
tests/hp/solution_transfer_08.cc
tests/hp/solution_transfer_09.cc
tests/hp/solution_transfer_10.cc
tests/hp/solution_transfer_11.cc
tests/hp/solution_transfer_12.cc
tests/hp/solution_transfer_13.cc
tests/hp/solution_transfer_14.cc
tests/hp/solution_transfer_15.cc
tests/hp/solution_transfer_16.cc
tests/mpi/solution_transfer_02.cc
tests/mpi/solution_transfer_03.cc
tests/mpi/solution_transfer_07.cc
tests/mpi/solution_transfer_08.cc
tests/mpi/solution_transfer_09.cc
tests/mpi/solution_transfer_10.cc
tests/petsc/solution_transfer_01.cc
tests/petsc/solution_transfer_02.cc
tests/simplex/step-31.cc
tests/sundials/step-77.cc
tests/trilinos/step-77-with-nox.cc

diff --git a/doc/news/changes/incompatibilities/20240103AfricaBlaisMunch b/doc/news/changes/incompatibilities/20240103AfricaBlaisMunch
new file mode 100644 (file)
index 0000000..a1c898a
--- /dev/null
@@ -0,0 +1,10 @@
+Deprecated: The classes SolutionTransfer and parallel::distributed::SolutionTransfer
+have been unified to SolutionTransfer. The class now supports both serial and parallel
+meshes. The class has lost some functions: the function SolutionTransfer::interpolate()
+that takes the input vector as well as the less frequently used functions
+SolutionTransfer::prepare_for_pure_refinement() and 
+SolutionTransfer::refine_interpolate(). Please use the other functions to accomplish the
+same functionality. For the time being, the old implementation has been moved to the
+`Legacy` namespace. The old parallel::distributed::SolutionTransfer has been now early deprecated.
+<br>
+(Pasquale Claudio Africa, Bruno Blais, Peter Munch, 2024/01/03)
index acb0a61febad59bfe3163c4164ace4c1f1adbf3a..991b42c3764f506a18d2d2ab1a1a20ece09240ef 100644 (file)
@@ -381,7 +381,7 @@ namespace Step15
 
     setup_system();
 
-    solution_transfer.interpolate(coarse_solution, current_solution);
+    solution_transfer.interpolate(current_solution);
 
     // On the new mesh, there are different hanging nodes, computed in
     // `setup_system()` above. To be on the safe side, we should  make sure that
index a1493d564829b8354ad79838dd7e0efec8a80833..77e0d1bec2bcae1d5ca03b210e596a1a4c34fb6d 100644 (file)
@@ -398,7 +398,7 @@ namespace Step26
     triangulation.execute_coarsening_and_refinement();
     setup_system();
 
-    solution_trans.interpolate(previous_solution, solution);
+    solution_trans.interpolate(solution);
     constraints.distribute(solution);
   }
 
index ad664acc22ae334317bc941c8e0931907b484996..4b3452850a41355799714b7b889ff029404b2f60 100644 (file)
@@ -1081,7 +1081,7 @@ namespace Step28
     setup_linear_system();
 
     solution.reinit(dof_handler.n_dofs());
-    soltrans.interpolate(solution_old, solution);
+    soltrans.interpolate(solution);
 
     // enforce constraints to make the interpolated solution conforming on
     // the new mesh:
index 08c3c88cd74703275c522e5ebcef745ea4dcedcc..6ee6303a89e4675d3b438725d7989ecbd4aee6ca 100644 (file)
@@ -1977,7 +1977,7 @@ namespace Step31
     std::vector<TrilinosWrappers::MPI::Vector> tmp = {
       TrilinosWrappers::MPI::Vector(temperature_solution),
       TrilinosWrappers::MPI::Vector(temperature_solution)};
-    temperature_trans.interpolate(x_temperature, tmp);
+    temperature_trans.interpolate(tmp);
 
     temperature_solution     = tmp[0];
     old_temperature_solution = tmp[1];
@@ -1991,7 +1991,7 @@ namespace Step31
     // we do not need another temporary vector since we just interpolate a
     // single vector. In the end, we have to tell the program that the matrices
     // and preconditioners need to be regenerated, since the mesh has changed.
-    stokes_trans.interpolate(x_stokes, stokes_solution);
+    stokes_trans.interpolate(stokes_solution);
 
     stokes_constraints.distribute(stokes_solution);
 
index 0610f0d6046f1c89014d77bf84c3ae2c0eb602eb..5ba1fa882e53e73ccc6fa048203cb3051b1c22a8 100644 (file)
@@ -2300,7 +2300,7 @@ namespace Step33
     std::vector<Vector<double>> transfer_out = {
       Vector<double>(dof_handler.n_dofs()),
       Vector<double>(dof_handler.n_dofs())};
-    soltrans.interpolate(transfer_in, transfer_out);
+    soltrans.interpolate(transfer_out);
 
     old_solution = std::move(transfer_out[0]);
     predictor    = std::move(transfer_out[1]);
index 36d22c13f89fb1f7495d93b255eb8d253dd3121a..2b82b6f765c958c47dcd635ca48efcd68af17a7d 100644 (file)
@@ -1688,7 +1688,7 @@ namespace Step43
       tmp_saturation[0].reinit(saturation_solution);
       tmp_saturation[1].reinit(saturation_solution);
       tmp_saturation[2].reinit(saturation_solution);
-      saturation_soltrans.interpolate(x_saturation, tmp_saturation);
+      saturation_soltrans.interpolate(tmp_saturation);
 
       saturation_solution                              = tmp_saturation[0];
       old_saturation_solution                          = tmp_saturation[1];
@@ -1702,7 +1702,7 @@ namespace Step43
       std::vector<TrilinosWrappers::MPI::BlockVector> tmp_darcy(2);
       tmp_darcy[0].reinit(darcy_solution);
       tmp_darcy[1].reinit(darcy_solution);
-      darcy_soltrans.interpolate(x_darcy, tmp_darcy);
+      darcy_soltrans.interpolate(tmp_darcy);
 
       last_computed_darcy_solution        = tmp_darcy[0];
       second_last_computed_darcy_solution = tmp_darcy[1];
index 29fe6ef6718a889adb8e89215e6d57c213115cfe..68f03e665de90bd37f22ba4379f2318edcf7887a 100644 (file)
@@ -605,7 +605,7 @@ namespace Step57
     // Transfer solution from coarse to fine mesh and apply boundary value
     // constraints to the new transferred solution. Note that present_solution
     // is still a vector corresponding to the old mesh.
-    solution_transfer.interpolate(present_solution, tmp);
+    solution_transfer.interpolate(tmp);
     nonzero_constraints.distribute(tmp);
 
     // Finally set up matrix and vectors and set the present_solution to the
index e7d42365ed1d3a3f26cc2118187ce0cde828c6f1..4eb70f2e8d4c2e7c8598429e27134e39738fd985 100644 (file)
@@ -811,7 +811,7 @@ namespace Step72
 
     setup_system();
 
-    solution_transfer.interpolate(coarse_solution, current_solution);
+    solution_transfer.interpolate(current_solution);
     nonzero_constraints.distribute(current_solution);
   }
 
index dc4d1436e90e7be7d42427b6ff20e2f0ee6d011a..50216ca92b981326331f5e2465fbceb3a8eb39de 100644 (file)
@@ -440,8 +440,7 @@ namespace Step77
 
     setup_system();
 
-    solution_transfer.interpolate(coarse_solution, current_solution);
-
+    solution_transfer.interpolate(current_solution);
     nonzero_constraints.distribute(current_solution);
   }
 
index caf9320da81a11d1100849c069fb967111685d13..857976b04df50243a581ed769e79e3d082dc867f 100644 (file)
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/distributed/tria.h>
-
-#include <deal.II/dofs/dof_handler.h>
-
-#include <vector>
+#include <deal.II/numerics/solution_transfer.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -31,348 +27,11 @@ namespace parallel
   namespace distributed
   {
     /**
-     * Transfer a discrete FE function (like a solution vector) by
-     * interpolation while refining and/or coarsening a distributed grid and
-     * handles the necessary communication.
-     *
-     * @note It is important to note, that if you use more than one
-     * SolutionTransfer object at the same time, that the calls to prepare_*()
-     * and interpolate()/deserialize() need to be in the same order.
-     *
-     * <h3>Note on ghost elements</h3> In a parallel computation PETSc or
-     * Trilinos vector may contain ghost elements or not. For reading in
-     * information with prepare_for_coarsening_and_refinement() or
-     * prepare_for_serialization() you need to supply vectors with ghost
-     * elements, so that all locally_active elements can be read. On the other
-     * hand, ghosted vectors are generally not writable, so for calls to
-     * interpolate() or deserialize() you need to supply distributed vectors
-     * without ghost elements. More precisely, during interpolation the
-     * current algorithm writes into all locally active degrees of freedom.
-     *
-     * <h3>Transferring a solution</h3> Here VectorType is your favorite
-     * vector type, e.g. PETScWrappers::MPI::Vector,
-     * TrilinosWrappers::MPI::Vector, or corresponding block vectors.
-     * @code
-     * parallel::distributed::SolutionTransfer<dim, VectorType>
-     *   soltrans(dof_handler);
-     * // flag some cells for refinement and coarsening, e.g.
-     * parallel::distributed::GridRefinement::
-     *   refine_and_coarsen_fixed_fraction(tria,
-     *                                     error_indicators,
-     *                                     0.3,
-     *                                     0.05);
-     *
-     * // prepare the triangulation,
-     * tria.prepare_coarsening_and_refinement();
-     *
-     * // prepare the SolutionTransfer object for coarsening and refinement
-     * // and give the solution vector that we intend to interpolate later,
-     * soltrans.prepare_for_coarsening_and_refinement(solution);
-     *
-     * // actually execute the refinement,
-     * tria.execute_coarsening_and_refinement();
-     *
-     * // redistribute dofs,
-     * dof_handler.distribute_dofs(fe);
-     *
-     * // and interpolate the solution
-     * VectorType interpolated_solution;
-     *
-     * //create VectorType in the right size here
-     * soltrans.interpolate(interpolated_solution);
-     * @endcode
-     *
-     * As the grid is distributed, it is important to note that the old
-     * solution(s) must be copied to one that also provides access to the
-     * locally relevant DoF values (these values required for the interpolation
-     * process):
-     * @code
-     * // Create initial indexsets pertaining to the grid before refinement
-     * const IndexSet &locally_owned_dofs    = dof_handler.locally_owned_dofs();
-     * const IndexSet  locally_relevant_dofs =
-     *   DoFTools::extract_locally_relevant_dofs(dof_handler);
-     *
-     * // The solution vector only knows about locally owned DoFs
-     * TrilinosWrappers::MPI::Vector solution;
-     * solution.reinit(locally_owned_dofs,
-     *                 mpi_communicator);
-     * ...
-     * // Transfer solution to vector that provides access to
-     * // locally relevant DoFs
-     * TrilinosWrappers::MPI::Vector old_solution;
-     * old_solution.reinit(locally_owned_dofs,
-     *                     locally_relevant_dofs,
-     *                     mpi_communicator);
-     * old_solution = solution;
-     *
-     * // Initialize SolutionTransfer object
-     * parallel::distributed::SolutionTransfer<dim, VectorType>
-     *   soltrans(dof_handler);
-     * soltrans.prepare_for_coarsening_and_refinement(old_solution);
-     * ...
-     * // Refine grid
-     * // Recreate locally_owned_dofs and locally_relevant_dofs index sets
-     * ...
-     * solution.reinit(locally_owned_dofs, mpi_communicator);
-     * soltrans.interpolate(solution);
-     * @endcode
-     *
-     * Different from PETSc and Trilinos vectors,
-     * LinearAlgebra::distributed::Vector allows writing into ghost elements.
-     * For a ghosted vector the interpolation step can be accomplished via
-     * @code
-     * interpolated_solution.zero_out_ghost_values();
-     * soltrans.interpolate(interpolated_solution);
-     * interpolated_solution.update_ghost_values();
-     * @endcode
-     *
-     * <h3>Use for Serialization</h3>
-     *
-     * This class can be used to serialize and later deserialize a distributed
-     * mesh with solution vectors to a file. If you use more than one
-     * DoFHandler and therefore more than one SolutionTransfer object, they
-     * need to be serialized and deserialized in the same order.
-     *
-     * If vector has the locally relevant DoFs, serialization works as
-     * follows:
-     * @code
-     * parallel::distributed::SolutionTransfer<dim, VectorType>
-     *   sol_trans(dof_handler);
-     * sol_trans.prepare_for_serialization(vector);
-     *
-     * triangulation.save(filename);
-     * @endcode
-     * For deserialization the vector needs to be a distributed vector
-     * (without ghost elements):
-     * @code
-     * //[create coarse mesh...]
-     * triangulation.load(filename);
-     *
-     * parallel::distributed::SolutionTransfer<dim, VectorType>
-     *   sol_trans(dof_handler);
-     * sol_trans.deserialize(distributed_vector);
-     * @endcode
-     *
-     *
-     * <h3>Note on usage with DoFHandler with hp-capabilities</h3>
-     *
-     * Since data on DoFHandler objects with hp-capabilities is associated with
-     * many different FiniteElement objects, each cell's data has to be
-     * processed with its corresponding `future_fe_index`. Further, if
-     * refinement is involved, data will be packed on the parent cell with its
-     * `future_fe_index` and unpacked later with the same index on its children.
-     * If cells get coarsened into one, data will be packed on the children with
-     * the least dominant finite element of their common subspace, and unpacked
-     * on the parent with this particular finite element (consult
-     * hp::FECollection::find_dominated_fe_extended() for more information).
-     *
-     * Transferring a solution across refinement works exactly like in the
-     * non-hp-case. However, when considering serialization, we also have to
-     * store the active FE indices in an additional step. A code snippet
-     * demonstrating serialization with the
-     * parallel::distributed::SolutionTransfer class with DoFHandler objects
-     * with hp-capabilities is provided in the following. Here VectorType is
-     * your favorite vector type, e.g. PETScWrappers::MPI::Vector,
-     * TrilinosWrappers::MPI::Vector, or corresponding block vectors.
-     *
-     * If vector has the locally relevant DoFs, serialization works as follows:
-     * @code
-     * parallel::distributed::
-     *   SolutionTransfer<dim, VectorType, DoFHandler<dim,spacedim>>
-     *     sol_trans(hp_dof_handler);
-     *
-     * hp_dof_handler.prepare_for_serialization_of_active_fe_indices();
-     * sol_trans.prepare_for_serialization(vector);
-     *
-     * triangulation.save(filename);
-     * @endcode
-     *
-     * For deserialization the vector needs to be a distributed vector
-     * (without ghost elements):
-     * @code
-     * //[create coarse mesh...]
-     * triangulation.load(filename);
-     *
-     * hp::FECollection<dim,spacedim> fe_collection;
-     * //[prepare identical fe_collection...]
-     *
-     * DoFHandler<dim,spacedim> hp_dof_handler(triangulation);
-     * // We need to introduce our dof_handler to the fe_collection
-     * // before setting all active FE indices.
-     * hp_dof_handler.deserialize_active_fe_indices();
-     * hp_dof_handler.distribute_dofs(fe_collection);
-     *
-     * parallel::distributed::
-     *   SolutionTransfer<dim,VectorType,DoFHandler<dim,spacedim>>
-     *     sol_trans(hp_dof_handler);
-     * sol_trans.deserialize(distributed_vector);
-     * @endcode
-     *
-     *
-     * <h3>Interaction with hanging nodes</h3>
-     *
-     * In essence, this class implements the same steps as does
-     * dealii::SolutionTransfer (though the implementation is entirely
-     * separate). Consequently, the same issue with hanging nodes and
-     * coarsening can happen with this class as happens with
-     * dealii::SolutionTransfer. See there for an extended discussion.
-     *
-     * @ingroup distributed
+     * @deprecated Use dealii::SolutionTransfer instead.
      */
     template <int dim, typename VectorType, int spacedim = dim>
-    class SolutionTransfer
-    {
-    public:
-      /**
-       * Constructor.
-       *
-       * @param[in] dof_handler 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.
-       * @param[in] average_values Average the contribututions to the same
-       *   DoF coming from different cells. Note: averaging requires an
-       * additional communication step, since the valence of the DoF has to be
-       * determined.
-       */
-      SolutionTransfer(const DoFHandler<dim, spacedim> &dof_handler,
-                       const bool                       average_values = false);
-
-      /**
-       * Destructor.
-       */
-      ~SolutionTransfer() = default;
-
-      /**
-       * Prepare the current object for coarsening and refinement. It
-       * stores the dof indices of each cell and stores the dof values of the
-       * vectors in @p all_in in each cell that'll be coarsened. @p all_in
-       * includes all vectors that are to be interpolated onto the new
-       * (refined and/or coarsened) grid.
-       */
-      void
-      prepare_for_coarsening_and_refinement(
-        const std::vector<const VectorType *> &all_in);
-
-      /**
-       * Same as the previous function but for only one discrete function to be
-       * interpolated.
-       */
-      void
-      prepare_for_coarsening_and_refinement(const VectorType &in);
-
-      /**
-       * 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
-      interpolate(std::vector<VectorType *> &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_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
-       *
-       * Multiple calling of this function is NOT allowed. Interpolating
-       * several functions can be performed in one step by using
-       * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> and
-       * <tt>interpolate(all_out)</tt>.
-       */
-      void
-      interpolate(VectorType &out);
-
-      /**
-       * Prepare the serialization of the given vector. The serialization is
-       * done by Triangulation::save(). The given vector needs all information
-       * on the locally active DoFs (it must be ghosted). See documentation of
-       * this class for more information.
-       */
-      void
-      prepare_for_serialization(const VectorType &in);
-
-      /**
-       * Same as the function above, only for a list of vectors.
-       */
-      void
-      prepare_for_serialization(const std::vector<const VectorType *> &all_in);
-
-      /**
-       * Execute the deserialization of the given vector. This needs to be
-       * done after calling Triangulation::load(). The given vector must be a
-       * fully distributed vector without ghost elements. See documentation of
-       * this class for more information.
-       */
-      void
-      deserialize(VectorType &in);
-
-
-      /**
-       * Same as the function above, only for a list of vectors.
-       */
-      void
-      deserialize(std::vector<VectorType *> &all_in);
-
-    private:
-      /**
-       * Pointer to the degree of freedom handler to work with.
-       */
-      SmartPointer<const DoFHandler<dim, spacedim>,
-                   SolutionTransfer<dim, VectorType, spacedim>>
-        dof_handler;
-
-      /**
-       * Flag indicating if averaging should be performed.
-       */
-      const bool average_values;
-
-      /**
-       * 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 VectorType *> input_vectors;
-
-      /**
-       * 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;
-
-      /**
-       * 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 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 CellStatus                                            status,
-        const boost::iterator_range<std::vector<char>::const_iterator>
-                                  &data_range,
-        std::vector<VectorType *> &all_out,
-        VectorType                &valence);
-
-
-      /**
-       * 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();
-    };
+    using SolutionTransfer DEAL_II_DEPRECATED_EARLY =
+      dealii::SolutionTransfer<dim, VectorType, spacedim>;
   } // namespace distributed
 } // namespace parallel
 
index c9f8fc6cba9eb401943238a9e6be3c580c4c7560..2ebec1f0777eaa4492075c91c4709f2aec094c39 100644 (file)
@@ -394,14 +394,14 @@ namespace parallel
     virtual void
     update_number_cache();
 
-    void
-    update_reference_cells() override;
-
     /**
      * Reset global active cell indices and global level cell indices.
      */
     void
     reset_global_cell_indices();
+
+    void
+    update_reference_cells() override;
   };
 
 
index e0a956b1c16482c2e47279ef97f8625c665c4ca3..bbd55ec65126511e068335b22b5ca48fe7ce0190 100644 (file)
@@ -3904,6 +3904,22 @@ protected:
   void
   update_cell_relations();
 
+  /**
+   * Function to pack data for
+   * SolutionTransfer::prepare_for_coarsening_and_refinement() in the case of a
+   * serial triangulation.
+   */
+  void
+  pack_data_serial();
+
+
+  /**
+   * Function to unpack data for SolutionTransfer::interpolate() in the case of
+   * a serial triangulation.
+   */
+  void
+  unpack_data_serial();
+
   /**
    * Vector of pairs, each containing a deal.II cell iterator and its
    * respective CellStatus. To update its contents, use the
index 6817b1d5fef969b818c4c5bc05a9fd552eedf497..0fb3a5c385bb3502ab02c4e4811bd6e2d273b6a7 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+
 /**
  * This class implements the transfer of a discrete FE function (e.g. a
  * solution vector) from one mesh to another that is obtained from the first
  * by a single refinement and/or coarsening step. During interpolation the
- * vector is reinitialized to the new size and filled with the interpolated
- * values. This class is used in the step-15, step-26, step-31, and step-33
- * tutorial programs. A version of this class that works on parallel
- * triangulations is available as parallel::distributed::SolutionTransfer.
+ * vector is filled with the interpolated values. This class is used in the
+ * step-15, step-26, step-31, and step-33 tutorial programs. This class
+ * works both for serial and distributed meshes.
  *
  * <h3>Usage</h3>
  *
- * This class implements the algorithms in two different ways:
- * <ul>
- * <li> If the grid will only be refined (i.e. no cells are coarsened) then
- * use @p SolutionTransfer as follows:
- * @code
- * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
- *
- * // flag some cells for refinement, e.g.
- * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
- *                                                   error_indicators,
- *                                                   0.3,
- *                                                   0);
- * // prepare the triangulation for refinement,
- * tria->prepare_coarsening_and_refinement();
- *
- * // tell the SolutionTransfer object that we intend to do pure refinement,
- * soltrans.prepare_for_pure_refinement();
- *
- * // actually execute the refinement,
- * tria->execute_coarsening_and_refinement();
- *
- * // and redistribute dofs.
- * dof_handler->distribute_dofs (fe);
- * @endcode
- *
- * Then to proceed do
- * @code
- * // take a copy of the solution vector
- * Vector<double> solution_old(solution);
- *
- * // resize solution vector to the correct size, as the @p refine_interpolate
- * // function requires the vectors to be of right sizes
- * solution.reinit(dof_handler->n_dofs());
- *
- * // and finally interpolate
- * soltrans.refine_interpolate(solution_old, solution);
- * @endcode
+ * @note It is important to note, that if you use more than one
+ * SolutionTransfer object at the same time, that the calls to prepare_*()
+ * and interpolate()/deserialize() need to be in the same order.
  *
- * Although the @p refine_interpolate functions are allowed to be called
- * multiple times, e.g. for interpolating several solution vectors, there is
- * the following possibility of interpolating several functions simultaneously.
+ * <h3>Transferring a solution</h3> Here VectorType is your favorite
+ * vector type, e.g. Vector, LinearAlgebra::distributed::Vector,
+ * PETScWrappers::MPI::Vector, TrilinosWrappers::MPI::Vector, or corresponding
+ * block vectors.
  * @code
- * std::vector<Vector<double> > solutions_old(n_vectors, Vector<double> (n));
- * ...
- * std::vector<Vector<double> > solutions(n_vectors, Vector<double> (n));
- * soltrans.refine_interpolate(solutions_old, solutions);
- * @endcode
- * This is used in several of the tutorial programs, for example step-31
- * and step-33.
- *
- * <li> If the grid has cells that will be coarsened, then use @p
- * SolutionTransfer as follows:
- * @code
- * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
- *
+ * SolutionTransfer<dim, VectorType> soltrans(dof_handler);
  * // flag some cells for refinement and coarsening, e.g.
- * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
- *                                                   error_indicators,
- *                                                   0.3,
- *                                                   0.05);
+ * GridRefinement::
+ *   refine_and_coarsen_fixed_fraction(tria,
+ *                                     error_indicators,
+ *                                     0.3,
+ *                                     0.05);
+ * // Use parallel::distributed::GridRefinement for distributed triangulations.
  *
  * // prepare the triangulation,
- * tria->prepare_coarsening_and_refinement();
+ * tria.prepare_coarsening_and_refinement();
  *
- * // prepare the SolutionTransfer object for coarsening and refinement and give
- * // the solution vector that we intend to interpolate later,
+ * // prepare the SolutionTransfer object for coarsening and refinement
+ * // and give the solution vector that we intend to interpolate later,
  * soltrans.prepare_for_coarsening_and_refinement(solution);
  *
  * // actually execute the refinement,
- * tria->execute_coarsening_and_refinement ();
+ * tria.execute_coarsening_and_refinement();
  *
  * // redistribute dofs,
- * dof_handler->distribute_dofs (fe);
+ * dof_handler.distribute_dofs(fe);
  *
  * // and interpolate the solution
- * Vector<double> interpolate_solution(dof_handler->n_dofs());
- * soltrans.interpolate(solution, interpolated_solution);
+ * VectorType interpolated_solution;
+ *
+ * //create VectorType in the right size here
+ * soltrans.interpolate(interpolated_solution);
  * @endcode
  *
- * If the grid is partitioned across several MPI processes, then it is
- * important to note that the old solution(s) must be copied to one that
- * also provides access to the locally relevant DoF values (these values
- * required for the interpolation process):
+ * <h3>Usage on distributed grids</h3>
+ * If the grid is distributed, it is important to note that the old
+ * solution(s) must be copied to one that also provides access to the
+ * locally relevant DoF values (these values required for the interpolation
+ * process):
  * @code
  * // Create initial indexsets pertaining to the grid before refinement
  * const IndexSet &locally_owned_dofs    = dof_handler.locally_owned_dofs();
@@ -137,116 +97,121 @@ DEAL_II_NAMESPACE_OPEN
  * solution.reinit(locally_owned_dofs,
  *                 mpi_communicator);
  * ...
- * // Transfer solution to vector that provides access to locally relevant DoFs
+ * // Transfer solution to vector that provides access to
+ * // locally relevant DoFs
  * TrilinosWrappers::MPI::Vector old_solution;
  * old_solution.reinit(locally_owned_dofs,
  *                     locally_relevant_dofs,
  *                     mpi_communicator);
  * old_solution = solution;
+ *
+ * // Initialize SolutionTransfer object
+ * SolutionTransfer<dim, VectorType> soltrans(dof_handler);
+ * soltrans.prepare_for_coarsening_and_refinement(old_solution);
  * ...
  * // Refine grid
  * // Recreate locally_owned_dofs and locally_relevant_dofs index sets
  * ...
  * solution.reinit(locally_owned_dofs, mpi_communicator);
- * soltrans.refine_interpolate(old_solution, solution);
+ * soltrans.interpolate(solution);
+ * @endcode
+ *
+ * <h4>Note on ghost elements</h4> In a parallel computation PETSc or
+ * Trilinos vector may contain ghost elements or not. For reading in
+ * information with prepare_for_coarsening_and_refinement() or
+ * prepare_for_serialization() you need to supply vectors with ghost
+ * elements, so that all locally_active elements can be read. On the other
+ * hand, ghosted vectors are generally not writable, so for calls to
+ * interpolate() or deserialize() you need to supply distributed vectors
+ * without ghost elements. More precisely, during interpolation the
+ * current algorithm writes into all locally active degrees of freedom.
+ *
+ * Different from PETSc and Trilinos vectors, LinearAlgebra::distributed::Vector
+ * allows writing into ghost elements. For a ghosted vector the interpolation
+ * step can be accomplished via
+ * @code
+ * interpolated_solution.zero_out_ghost_values();
+ * soltrans.interpolate(interpolated_solution);
+ * interpolated_solution.update_ghost_values();
  * @endcode
  *
- * Multiple calls to the function <code>interpolate (const VectorType &in,
- * VectorType &out)</code> are NOT allowed. Interpolating several
- * functions can be performed in one step by using <tt>void interpolate (const
- * vector<VectorType> &all_in, vector<VectorType> &all_out)
- * const</tt>, and using the respective @p
- * prepare_for_coarsening_and_refinement function taking several vectors as
- * input before actually refining and coarsening the triangulation (see
- * there).
- * </ul>
- *
- * For deleting all stored data in @p SolutionTransfer and reinitializing it
- * use the <tt>clear()</tt> function.
- *
- * The template argument @p VectorType denotes the type of data container you
- * want to transfer.
- *
- *
- * <h3>Interpolating in the presence of hanging nodes and boundary values</h3>
- *
- * The interpolation onto the new mesh is a local operation, i.e., it
- * interpolates onto the new mesh only. If that new mesh has hanging nodes,
- * you will therefore get a solution that does not satisfy hanging node
- * constraints. The same is true with boundary values: the interpolated
- * solution will just be the interpolation of the old solution at the
- * boundary, and this may or may not satisfy boundary values at newly
- * introduced boundary nodes.
- *
- * Consequently, you may have to apply hanging node or boundary value
- * constraints after interpolation. step-15 and step-26 have examples of
- * dealing with this.
- *
- *
- * <h3>Implementation</h3>
- *
- * <ul>
- * <li> Solution transfer with only refinement. Assume that we have got a
- * solution vector on the current (original) grid. Each entry of this vector
- * belongs to one of the DoFs of the discretization. If we now refine the grid
- * then the calling of DoFHandler::distribute_dofs() will change at least some
- * of the DoF indices. Hence we need to store the DoF indices of all active
- * cells before the refinement. A pointer for each active cell is used to
- * point to the vector of these DoF indices of that cell. This is done by
- * prepare_for_pure_refinement().
- *
- * In the function <tt>refine_interpolate(in,out)</tt> and on each cell where
- * the pointer is set (i.e. the cells that were active in the original grid)
- * we can now access the local values of the solution vector @p in on that
- * cell by using the stored DoF indices. These local values are interpolated
- * and set into the vector @p out that is at the end the discrete function @p
- * in interpolated on the refined mesh.
- *
- * The <tt>refine_interpolate(in,out)</tt> function can be called multiple
- * times for arbitrary many discrete functions (solution vectors) on the
- * original grid.
- *
- * <li> Solution transfer with coarsening and refinement. After calling
- * Triangulation::prepare_coarsening_and_refinement the coarsen flags of
- * either all or none of the children of a (father-)cell are set. While
- * coarsening (Triangulation::execute_coarsening_and_refinement) the cells
- * that are not needed any more will be deleted from the Triangulation.
- *
- * For the interpolation from the (to be coarsenend) children to their father
- * the children cells are needed. Hence this interpolation and the storing of
- * the interpolated values of each of the discrete functions that we want to
- * interpolate needs to take place before these children cells are coarsened
- * (and deleted!!). Again a pointer for each relevant cell is set to point to
- * these values (see below). Additionally the DoF indices of the cells that
- * will not be coarsened need to be stored according to the solution transfer
- * with pure refinement (cf there). All this is performed by
- * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
- * <tt>vector<VectorType> all_in</tt> includes all discrete
- * functions to be interpolated onto the new grid.
- *
- * As we need two different kinds of pointers (<tt>vector<unsigned int> *</tt>
- * for the Dof indices and <tt>vector<VectorType> *</tt> for the
- * interpolated DoF values) we use the @p Pointerstruct that includes both of
- * these pointers and the pointer for each cell points to these @p
- * Pointerstructs. On each cell only one of the two different pointers is used
- * at one time hence we could use a <tt>void * pointer</tt> as
- * <tt>vector<unsigned int> *</tt> at one time and as
- * <tt>vector<VectorType> *</tt> at the other but using this @p
- * Pointerstruct in between makes the use of these pointers more safe and
- * gives better possibility to expand their usage.
- *
- * In <tt>interpolate(all_in, all_out)</tt> the refined cells are treated
- * according to the solution transfer while pure refinement. Additionally, on
- * each cell that is coarsened (hence previously was a father cell), the
- * values of the discrete functions in @p all_out are set to the stored local
- * interpolated values that are accessible due to the 'vector<VectorType>
- * *' pointer in @p Pointerstruct that is pointed to by the pointer of that
- * cell. It is clear that <tt>interpolate(all_in, all_out)</tt> only can be
- * called with the <tt>vector<VectorType> all_in</tt> that previously was
- * the parameter of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt>
- * function. Hence <tt>interpolate(all_in, all_out)</tt> can (in contrast to
- * <tt>refine_interpolate(in, out)</tt>) only be called once.
- * </ul>
+ * <h3>Use for serialization</h3>
+ *
+ * This class can be used to serialize and later deserialize a
+ * mesh with solution vectors to a file. If you use more than one
+ * DoFHandler and therefore more than one SolutionTransfer object, they
+ * need to be serialized and deserialized in the same order.
+ *
+ * If vector has the locally relevant DoFs, serialization works as
+ * follows:
+ * @code
+ * SolutionTransfer<dim, VectorType> sol_trans(dof_handler);
+ * sol_trans.prepare_for_serialization(vector);
+ *
+ * triangulation.save(filename);
+ * @endcode
+ * For deserialization the vector needs to be a distributed vector
+ * (without ghost elements):
+ * @code
+ * //[create coarse mesh...]
+ * triangulation.load(filename);
+ *
+ * SolutionTransfer<dim, VectorType> sol_trans(dof_handler);
+ * sol_trans.deserialize(distributed_vector);
+ * @endcode
+ *
+ *
+ * <h3>Note on usage with DoFHandler with hp-capabilities</h3>
+ *
+ * Since data on DoFHandler objects with hp-capabilities is associated with
+ * many different FiniteElement objects, each cell's data has to be
+ * processed with its corresponding `future_fe_index`. Further, if
+ * refinement is involved, data will be packed on the parent cell with its
+ * `future_fe_index` and unpacked later with the same index on its children.
+ * If cells get coarsened into one, data will be packed on the children with
+ * the least dominant finite element of their common subspace, and unpacked
+ * on the parent with this particular finite element (consult
+ * hp::FECollection::find_dominated_fe_extended() for more information).
+ *
+ * Transferring a solution across refinement works exactly like in the
+ * non-hp-case. However, when considering serialization, we also have to
+ * store the active FE indices in an additional step. A code snippet
+ * demonstrating serialization with the SolutionTransfer class with DoFHandler
+ * objects with hp-capabilities is provided in the following. Here VectorType is
+ * your favorite vector type, e.g. Vector, LinearAlgebra::distributed::Vector,
+ * PETScWrappers::MPI::Vector, TrilinosWrappers::MPI::Vector, or corresponding
+ * block vectors.
+ *
+ * If vector has the locally relevant DoFs, serialization works as follows:
+ * @code
+ * SolutionTransfer<dim, VectorType, DoFHandler<dim,spacedim>>
+ *     sol_trans(hp_dof_handler);
+ *
+ * hp_dof_handler.prepare_for_serialization_of_active_fe_indices();
+ * sol_trans.prepare_for_serialization(vector);
+ *
+ * triangulation.save(filename);
+ * @endcode
+ *
+ * For deserialization the vector needs to be a distributed vector
+ * (without ghost elements):
+ * @code
+ * //[create coarse mesh...]
+ * triangulation.load(filename);
+ *
+ * hp::FECollection<dim,spacedim> fe_collection;
+ * //[prepare identical fe_collection...]
+ *
+ * DoFHandler<dim,spacedim> hp_dof_handler(triangulation);
+ * // We need to introduce our dof_handler to the fe_collection
+ * // before setting all active FE indices.
+ * hp_dof_handler.deserialize_active_fe_indices();
+ * hp_dof_handler.distribute_dofs(fe_collection);
+ *
+ * SolutionTransfer<dim,VectorType,DoFHandler<dim,spacedim>>
+ *     sol_trans(hp_dof_handler);
+ * sol_trans.deserialize(distributed_vector);
+ * @endcode
  *
  *
  * <h3>Interaction with hanging nodes</h3>
@@ -293,34 +258,34 @@ DEAL_II_NAMESPACE_OPEN
  *
  * <h3>Implementation in the context of hp-finite elements</h3>
  *
- * In the case of DoFHandlers with hp-capabilities, nothing defines which of the
- * finite elements that are part of the hp::FECollection associated with the
- * DoFHandler, should be considered on cells that are not active (i.e., that
- * have children). This is because degrees of freedom are only allocated for
- * active cells and, in fact, it is not allowed to set an active FE index on
- * non-active cells using DoFAccessor::set_active_fe_index().
+ * In the case of DoFHandlers with hp-capabilities, nothing defines which of
+ * the finite elements that are part of the hp::FECollection associated with
+ * the DoFHandler, should be considered on cells that are not active (i.e.,
+ * that have children). This is because degrees of freedom are only allocated
+ * for active cells and, in fact, it is not allowed to set an active FE index
+ * on non-active cells using DoFAccessor::set_active_fe_index().
  *
  * It is, thus, not entirely natural what should happen if, for example, a few
  * cells are coarsened away. This class then implements the following
  * algorithm:
  * - If a cell is refined, then the values of the solution vector(s) are
- *   interpolated before refinement on the to-be-refined cell from the space of
- *   the active finite element to the one of the future finite element. These
+ *   interpolated before refinement on the to-be-refined cell from the space
+ * of the active finite element to the one of the future finite element. These
  *   values are then distributed on the finite element spaces of the children
  *   post-refinement. This may lose information if, for example, the old cell
  *   used a Q2 space and the children use Q1 spaces, or the information may be
  *   prolonged if the mother cell used a Q1 space and the children are Q2s.
  * - If cells are to be coarsened, then the values from the child cells are
- *   interpolated to the mother cell using the largest of the child cell future
- *   finite element spaces, which will be identified as the least dominant
- *   element following the FiniteElementDomination logic (consult
+ *   interpolated to the mother cell using the largest of the child cell
+ * future finite element spaces, which will be identified as the least
+ * dominant element following the FiniteElementDomination logic (consult
  *   hp::FECollection::find_dominated_fe_extended() for more information). For
  *   example, if the children of a cell use Q1, Q2 and Q3 spaces, then the
  *   values from the children are interpolated into a Q3 space on the mother
  *   cell. After refinement, this Q3 function on the mother cell is then
- *   interpolated into the space the user has selected for this cell (which may
- *   be different from Q3, in this example, if the user has set the
- *   active FE index for a different space post-refinement and before calling
+ *   interpolated into the space the user has selected for this cell (which
+ * may be different from Q3, in this example, if the user has set the active
+ * FE index for a different space post-refinement and before calling
  *   DoFHandler::distribute_dofs()).
  *
  * @note In the context of hp-refinement, if cells are coarsened or the
@@ -337,132 +302,113 @@ class SolutionTransfer
 {
 public:
   /**
-   * Constructor, takes the current DoFHandler as argument.
+   * Constructor.
+   *
+   * @param[in] dof_handler 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.
+   * @param[in] average_values Average the contribututions to the same
+   *   DoF coming from different cells. Note: averaging requires an
+   * additional communication step, since the valence of the DoF has to be
+   * determined.
    */
-  SolutionTransfer(const DoFHandler<dim, spacedim> &dof);
+  SolutionTransfer(const DoFHandler<dim, spacedim> &dof_handler,
+                   const bool                       average_values = false);
 
   /**
-   * Destructor
+   * Destructor.
    */
-  ~SolutionTransfer();
+  ~SolutionTransfer() = default;
 
   /**
-   * Reinit this class to the state that it has directly after calling the
-   * Constructor
+   * Prepare the current object for coarsening and refinement. It
+   * stores pointers to the vectors in @p all_in to be interpolated onto the new
+   * (refined and/or coarsened) grid, and registers this object for data
+   * transfer on the grid.
    */
   void
-  clear();
+  prepare_for_coarsening_and_refinement(
+    const std::vector<const VectorType *> &all_in);
 
   /**
-   * Prepares the @p SolutionTransfer for pure refinement. It stores the dof
-   * indices of each cell. After calling this function only calling the @p
-   * refine_interpolate functions is allowed.
-   */
-  void
-  prepare_for_pure_refinement();
-
-  /**
-   * Prepares the @p SolutionTransfer for coarsening and refinement. It stores
-   * the dof indices of each cell and stores the dof values of the vectors in
-   * @p all_in in each cell that'll be coarsened. @p all_in includes all
-   * vectors that are to be interpolated onto the new (refined and/or
-   * coarsenend) grid.
+   * Same as above but without pointers.
    */
   void
   prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in);
 
   /**
-   * Same as previous function but for only one discrete function to be
+   * Same as the previous function but for only one discrete function to be
    * interpolated.
    */
   void
   prepare_for_coarsening_and_refinement(const VectorType &in);
 
   /**
-   * This function interpolates the discrete function @p in, which is a vector
-   * on the grid before the refinement, to the function @p out which then is a
-   * vector on the refined grid. It assumes the vectors having the right sizes
-   * (i.e. <tt>in.size()==n_dofs_old</tt>,
-   * <tt>out.size()==n_dofs_refined</tt>)
-   *
-   * Calling this function is allowed only if @p prepare_for_pure_refinement
-   * is called and the refinement is executed before. Multiple calling of this
-   * function is allowed. e.g. for interpolating several functions.
+   * 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
-  refine_interpolate(const VectorType &in, VectorType &out) const;
+  interpolate(std::vector<VectorType *> &all_out);
 
   /**
-   * This function interpolates the discrete functions that are stored in @p
-   * all_in onto the refined and/or coarsenend grid. It assumes the vectors in
-   * @p all_in denote the same vectors as in @p all_in as parameter of
-   * <tt>prepare_for_refinement_and_coarsening(all_in)</tt>. However, there is
-   * no way of verifying this internally, so be careful here.
-   *
-   * Calling this function is allowed only if first
-   * Triangulation::prepare_coarsening_and_refinement, second @p
-   * SolutionTransfer::prepare_for_coarsening_and_refinement, an then third
-   * Triangulation::execute_coarsening_and_refinement are called before.
-   * Multiple calling of this function is NOT allowed. Interpolating several
-   * functions can be performed in one step.
-   *
-   * The number of output vectors is assumed to be the same as the number of
-   * input vectors. Also, the sizes of the output vectors are assumed to be of
-   * the right size (@p n_dofs_refined). Otherwise an assertion will be
-   * thrown.
+   * Same as above but without pointers.
    */
   void
-  interpolate(const std::vector<VectorType> &all_in,
-              std::vector<VectorType>       &all_out) const;
+  interpolate(std::vector<VectorType> &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_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
    *
-   * Multiple calling of this function is NOT allowed. Interpolating several
-   * functions can be performed in one step by using <tt>interpolate (all_in,
-   * all_out)</tt>
+   * Multiple calling of this function is NOT allowed. Interpolating
+   * several functions can be performed in one step by using
+   * <tt>interpolate (all_in, all_out)</tt>
+   */
+  void
+  interpolate(VectorType &out);
+
+  /**
+   * Prepare the serialization of the given vector. The serialization is
+   * done by Triangulation::save(). The given vector needs all information
+   * on the locally active DoFs (it must be ghosted). See documentation of
+   * this class for more information.
    */
   void
-  interpolate(const VectorType &in, VectorType &out) const;
+  prepare_for_serialization(const VectorType &in);
 
   /**
-   * Determine an estimate for the memory consumption (in bytes) of this
-   * object.
+   * Same as the function above, only for a list of vectors.
    */
-  std::size_t
-  memory_consumption() const;
+  void
+  prepare_for_serialization(const std::vector<const VectorType *> &all_in);
 
   /**
-   * Exception
+   * Execute the deserialization of the given vector. This needs to be
+   * done after calling Triangulation::load(). The given vector must be a
+   * fully distributed vector without ghost elements. See documentation of
+   * this class for more information.
    */
-  DeclExceptionMsg(ExcNotPrepared,
-                   "You are attempting an operation for which this object is "
-                   "not prepared. This may be because you either did not call "
-                   "one of the prepare_*() functions at all, or because you "
-                   "called the wrong one for the operation you are currently "
-                   "attempting.");
+  void
+  deserialize(VectorType &in);
+
 
   /**
-   * Exception
+   * Same as the function above, only for a list of vectors.
    */
-  DeclExceptionMsg(
-    ExcAlreadyPrepForRef,
-    "You are attempting to call one of the prepare_*() functions "
-    "of this object to prepare it for an operation for which it "
-    "is already prepared. Specifically, the object was "
-    "previously prepared for pure refinement.");
+  void
+  deserialize(std::vector<VectorType *> &all_in);
 
   /**
-   * Exception
+   * Reinit this class to the state that it has directly after calling the
+   * constructor.
    */
-  DeclExceptionMsg(
-    ExcAlreadyPrepForCoarseAndRef,
-    "You are attempting to call one of the prepare_*() functions "
-    "of this object to prepare it for an operation for which it "
-    "is already prepared. Specifically, the object was "
-    "previously prepared for both coarsening and refinement.");
+  DEAL_II_DEPRECATED_EARLY void
+  clear();
 
 private:
   /**
@@ -473,99 +419,599 @@ private:
     dof_handler;
 
   /**
-   * Stores the number of DoFs before the refinement and/or coarsening.
+   * Flag indicating if averaging should be performed.
    */
-  types::global_dof_index n_dofs_old;
+  const bool average_values;
 
   /**
-   * Declaration of @p PreparationState that denotes the three possible states
-   * of the @p SolutionTransfer: being prepared for 'pure refinement',
-   * prepared for 'coarsening and refinement' or not prepared.
+   * A vector that stores pointers to all the vectors we are supposed to
+   * copy over from the old to the new mesh.
    */
-  enum PreparationState
-  {
-    /**
-     * The SolutionTransfer is not yet prepared.
-     */
-    none,
-    /**
-     * The SolutionTransfer is prepared for purely refinement.
-     */
-    pure_refinement,
-    /**
-     * The SolutionTransfer is prepared for coarsening and refinement.
-     */
-    coarsening_and_refinement
-  };
+  std::vector<const VectorType *> input_vectors;
 
   /**
-   * Definition of the respective variable.
+   * The handle that the Triangulation has assigned to this object
+   * with which we can access our memory offset and our pack function.
    */
-  PreparationState prepared_for;
-
+  unsigned int handle;
 
   /**
-   * Is used for @p prepare_for_refining (of course also for @p
-   * repare_for_refining_and_coarsening) and stores all dof indices of the
-   * cells that'll be refined
+   * 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<std::vector<types::global_dof_index>> indices_on_cell;
+  std::vector<char>
+  pack_callback(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const CellStatus                                            status);
 
   /**
-   * All cell data (the dof indices and the dof values) should be accessible
-   * from each cell. As each cell has got only one @p user_pointer, multiple
-   * pointers to the data need to be packetized in a structure. Note that in
-   * our case on each cell either the <tt>vector<unsigned int> indices</tt>
-   * (if the cell will be refined) or the <tt>vector<double> dof_values</tt>
-   * (if the children of this cell will be deleted) is needed, hence one @p
-   * user_pointer should be sufficient, but to allow some error checks and to
-   * preserve the user from making user errors the @p user_pointer will be
-   * 'multiplied' by this structure.
+   * 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.
    */
-  struct Pointerstruct
-  {
-    Pointerstruct()
-      : indices_ptr(nullptr)
-      , dof_values_ptr(nullptr)
-      , active_fe_index(0)
-    {}
-    Pointerstruct(std::vector<types::global_dof_index> *indices_ptr_in,
-                  const unsigned int                    active_fe_index_in = 0)
-      : indices_ptr(indices_ptr_in)
-      , dof_values_ptr(nullptr)
-      , active_fe_index(active_fe_index_in)
-    {}
-    Pointerstruct(
-      std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr_in,
-      const unsigned int active_fe_index_in = 0)
-      : indices_ptr(nullptr)
-      , dof_values_ptr(dof_values_ptr_in)
-      , active_fe_index(active_fe_index_in)
-    {}
-    std::size_t
-    memory_consumption() const;
+  void
+  unpack_callback(
+    const typename Triangulation<dim, spacedim>::cell_iterator     &cell,
+    const CellStatus                                                status,
+    const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+    std::vector<VectorType *>                                      &all_out,
+    VectorType                                                     &valence);
 
-    std::vector<types::global_dof_index>                 *indices_ptr;
-    std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr;
-    unsigned int                                          active_fe_index;
-  };
 
   /**
-   * Map mapping from level and index of cell to the @p Pointerstructs (cf.
-   * there). This map makes it possible to keep all the information needed to
-   * transfer the solution inside this object rather than using user pointers
-   * of the Triangulation for this purpose.
+   * Registers the pack_callback() function to the Triangulation that has been
+   * assigned to the DoFHandler class member and stores the returning handle.
    */
-  std::map<std::pair<unsigned int, unsigned int>, Pointerstruct> cell_map;
+  void
+  register_data_attach();
+};
 
+
+
+namespace Legacy
+{
   /**
-   * Is used for @p prepare_for_refining_and_coarsening The interpolated dof
-   * values of all cells that'll be coarsened will be stored in this vector.
+   * This class implements the transfer of a discrete FE function (e.g. a
+   * solution vector) from one mesh to another that is obtained from the first
+   * by a single refinement and/or coarsening step. During interpolation the
+   * vector is reinitialized to the new size and filled with the interpolated
+   * values. This class is used in the step-15, step-26, step-31, and step-33
+   * tutorial programs. A version of this class that works on all types of
+   * triangulations, including distributed ones, is available as
+   * dealii::SolutionTransfer.
+   *
+   * <h3>Usage</h3>
+   *
+   * This class implements the algorithms in two different ways:
+   * <ul>
+   * <li> If the grid will only be refined (i.e. no cells are coarsened) then
+   * use @p SolutionTransfer as follows:
+   * @code
+   * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
+   *
+   * // flag some cells for refinement, e.g.
+   * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
+   *                                                   error_indicators,
+   *                                                   0.3,
+   *                                                   0);
+   * // prepare the triangulation for refinement,
+   * tria->prepare_coarsening_and_refinement();
+   *
+   * // tell the SolutionTransfer object that we intend to do pure refinement,
+   * soltrans.prepare_for_pure_refinement();
+   *
+   * // actually execute the refinement,
+   * tria->execute_coarsening_and_refinement();
+   *
+   * // and redistribute dofs.
+   * dof_handler->distribute_dofs (fe);
+   * @endcode
+   *
+   * Then to proceed do
+   * @code
+   * // take a copy of the solution vector
+   * Vector<double> solution_old(solution);
+   *
+   * // resize solution vector to the correct size, as the @p refine_interpolate
+   * // function requires the vectors to be of right sizes
+   * solution.reinit(dof_handler->n_dofs());
+   *
+   * // and finally interpolate
+   * soltrans.refine_interpolate(solution_old, solution);
+   * @endcode
+   *
+   * Although the @p refine_interpolate functions are allowed to be called
+   * multiple times, e.g. for interpolating several solution vectors, there is
+   * the following possibility of interpolating several functions
+   * simultaneously.
+   * @code
+   * std::vector<Vector<double> > solutions_old(n_vectors, Vector<double> (n));
+   * ...
+   * std::vector<Vector<double> > solutions(n_vectors, Vector<double> (n));
+   * soltrans.refine_interpolate(solutions_old, solutions);
+   * @endcode
+   * This is used in several of the tutorial programs, for example step-31
+   * and step-33.
+   *
+   * <li> If the grid has cells that will be coarsened, then use @p
+   * SolutionTransfer as follows:
+   * @code
+   * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
+   *
+   * // flag some cells for refinement and coarsening, e.g.
+   * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
+   *                                                   error_indicators,
+   *                                                   0.3,
+   *                                                   0.05);
+   *
+   * // prepare the triangulation,
+   * tria->prepare_coarsening_and_refinement();
+   *
+   * // prepare the SolutionTransfer object for coarsening and refinement and
+   * give
+   * // the solution vector that we intend to interpolate later,
+   * soltrans.prepare_for_coarsening_and_refinement(solution);
+   *
+   * // actually execute the refinement,
+   * tria->execute_coarsening_and_refinement ();
+   *
+   * // redistribute dofs,
+   * dof_handler->distribute_dofs (fe);
+   *
+   * // and interpolate the solution
+   * Vector<double> interpolate_solution(dof_handler->n_dofs());
+   * soltrans.interpolate(solution, interpolated_solution);
+   * @endcode
+   *
+   * If the grid is partitioned across several MPI processes, then it is
+   * important to note that the old solution(s) must be copied to one that
+   * also provides access to the locally relevant DoF values (these values
+   * required for the interpolation process):
+   * @code
+   * // Create initial indexsets pertaining to the grid before refinement
+   * const IndexSet &locally_owned_dofs    = dof_handler.locally_owned_dofs();
+   * const IndexSet  locally_relevant_dofs =
+   *   DoFTools::extract_locally_relevant_dofs(dof_handler);
+   *
+   * // The solution vector only knows about locally owned DoFs
+   * TrilinosWrappers::MPI::Vector solution;
+   * solution.reinit(locally_owned_dofs,
+   *                 mpi_communicator);
+   * ...
+   * // Transfer solution to vector that provides access to locally relevant
+   * DoFs TrilinosWrappers::MPI::Vector old_solution;
+   * old_solution.reinit(locally_owned_dofs,
+   *                     locally_relevant_dofs,
+   *                     mpi_communicator);
+   * old_solution = solution;
+   * ...
+   * // Refine grid
+   * // Recreate locally_owned_dofs and locally_relevant_dofs index sets
+   * ...
+   * solution.reinit(locally_owned_dofs, mpi_communicator);
+   * soltrans.refine_interpolate(old_solution, solution);
+   * @endcode
+   *
+   * Multiple calls to the function <code>interpolate (const VectorType &in,
+   * VectorType &out)</code> are NOT allowed. Interpolating several
+   * functions can be performed in one step by using <tt>void interpolate (const
+   * vector<VectorType> &all_in, vector<VectorType> &all_out)
+   * const</tt>, and using the respective @p
+   * prepare_for_coarsening_and_refinement function taking several vectors as
+   * input before actually refining and coarsening the triangulation (see
+   * there).
+   * </ul>
+   *
+   * For deleting all stored data in @p SolutionTransfer and reinitializing it
+   * use the <tt>clear()</tt> function.
+   *
+   * The template argument @p VectorType denotes the type of data container you
+   * want to transfer.
+   *
+   *
+   * <h3>Interpolating in the presence of hanging nodes and boundary values</h3>
+   *
+   * The interpolation onto the new mesh is a local operation, i.e., it
+   * interpolates onto the new mesh only. If that new mesh has hanging nodes,
+   * you will therefore get a solution that does not satisfy hanging node
+   * constraints. The same is true with boundary values: the interpolated
+   * solution will just be the interpolation of the old solution at the
+   * boundary, and this may or may not satisfy boundary values at newly
+   * introduced boundary nodes.
+   *
+   * Consequently, you may have to apply hanging node or boundary value
+   * constraints after interpolation. step-15 and step-26 have examples of
+   * dealing with this.
+   *
+   *
+   * <h3>Implementation</h3>
+   *
+   * <ul>
+   * <li> Solution transfer with only refinement. Assume that we have got a
+   * solution vector on the current (original) grid. Each entry of this vector
+   * belongs to one of the DoFs of the discretization. If we now refine the grid
+   * then the calling of DoFHandler::distribute_dofs() will change at least some
+   * of the DoF indices. Hence we need to store the DoF indices of all active
+   * cells before the refinement. A pointer for each active cell is used to
+   * point to the vector of these DoF indices of that cell. This is done by
+   * prepare_for_pure_refinement().
+   *
+   * In the function <tt>refine_interpolate(in,out)</tt> and on each cell where
+   * the pointer is set (i.e. the cells that were active in the original grid)
+   * we can now access the local values of the solution vector @p in on that
+   * cell by using the stored DoF indices. These local values are interpolated
+   * and set into the vector @p out that is at the end the discrete function @p
+   * in interpolated on the refined mesh.
+   *
+   * The <tt>refine_interpolate(in,out)</tt> function can be called multiple
+   * times for arbitrary many discrete functions (solution vectors) on the
+   * original grid.
+   *
+   * <li> Solution transfer with coarsening and refinement. After calling
+   * Triangulation::prepare_coarsening_and_refinement the coarsen flags of
+   * either all or none of the children of a (father-)cell are set. While
+   * coarsening (Triangulation::execute_coarsening_and_refinement) the cells
+   * that are not needed any more will be deleted from the Triangulation.
+   *
+   * For the interpolation from the (to be coarsenend) children to their father
+   * the children cells are needed. Hence this interpolation and the storing of
+   * the interpolated values of each of the discrete functions that we want to
+   * interpolate needs to take place before these children cells are coarsened
+   * (and deleted!!). Again a pointer for each relevant cell is set to point to
+   * these values (see below). Additionally the DoF indices of the cells that
+   * will not be coarsened need to be stored according to the solution transfer
+   * with pure refinement (cf there). All this is performed by
+   * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
+   * <tt>vector<VectorType> all_in</tt> includes all discrete
+   * functions to be interpolated onto the new grid.
+   *
+   * As we need two different kinds of pointers (<tt>vector<unsigned int> *</tt>
+   * for the Dof indices and <tt>vector<VectorType> *</tt> for the
+   * interpolated DoF values) we use the @p Pointerstruct that includes both of
+   * these pointers and the pointer for each cell points to these @p
+   * Pointerstructs. On each cell only one of the two different pointers is used
+   * at one time hence we could use a <tt>void * pointer</tt> as
+   * <tt>vector<unsigned int> *</tt> at one time and as
+   * <tt>vector<VectorType> *</tt> at the other but using this @p
+   * Pointerstruct in between makes the use of these pointers more safe and
+   * gives better possibility to expand their usage.
+   *
+   * In <tt>interpolate(all_in, all_out)</tt> the refined cells are treated
+   * according to the solution transfer while pure refinement. Additionally, on
+   * each cell that is coarsened (hence previously was a father cell), the
+   * values of the discrete functions in @p all_out are set to the stored local
+   * interpolated values that are accessible due to the 'vector<VectorType>
+   * *' pointer in @p Pointerstruct that is pointed to by the pointer of that
+   * cell. It is clear that <tt>interpolate(all_in, all_out)</tt> only can be
+   * called with the <tt>vector<VectorType> all_in</tt> that previously was
+   * the parameter of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt>
+   * function. Hence <tt>interpolate(all_in, all_out)</tt> can (in contrast to
+   * <tt>refine_interpolate(in, out)</tt>) only be called once.
+   * </ul>
+   *
+   *
+   * <h3>Interaction with hanging nodes</h3>
+   *
+   * This class does its best to represent on the new mesh the finite element
+   * function that existed on the old mesh, but this may lead to situations
+   * where the function on the new mesh is no longer conforming at hanging
+   * nodes. To this end, consider a situation of a twice refined mesh that
+   * started with a single square cell (i.e., we now have 16 cells). Consider
+   * also that we coarsen 4 of the cells back to the first refinement level. In
+   * this case, we end up with a mesh that will look as follows if we were to
+   * use a $Q_1$ element:
+   *
+   * @image html hanging_nodes.png ""
+   *
+   * The process of interpolating from the old to the new mesh would imply that
+   * the values of the finite element function will not change on all of the
+   * cells that remained as they are (i.e., the fine cells) but that on the
+   * coarse cell at the top right, the four values at the vertices are obtained
+   * by interpolating down from its former children.  If the original function
+   * was not linear, this implies that the marked hanging nodes will retain
+   * their old values which, in general, will not lead to a continuous function
+   * along the corresponding edges. In other words, the solution vector obtained
+   * after SolutionTransfer::interpolate() does not satisfy hanging node
+   * constraints: it corresponds to the pointwise interpolation, but not to the
+   * interpolation <i>onto the new finite element space that contains
+   * constraints from hanging nodes</i>.
+   *
+   * Whether this is a problem you need to worry about or not depends on your
+   * application. The situation is easily corrected, of course, by applying
+   * AffineConstraints::distribute() to your solution vector after transfer,
+   * using a constraints object computed on the new DoFHandler object (you
+   * probably need to create this object anyway if you have hanging nodes). This
+   * is also what is done, for example, in step-15.
+   *
+   * @note This situation can only happen if you do coarsening. If all cells
+   * remain as they are or are refined, then SolutionTransfer::interpolate()
+   * computes a new vector of nodal values, but the function represented is of
+   * course exactly the same because the old finite element space is a subspace
+   * of the new one. Thus, if the old function was conforming (i.e., satisfied
+   * hanging node constraints), then so does the new one, and it is not
+   * necessary to call AffineConstraints::distribute().
+   *
+   *
+   * <h3>Implementation in the context of hp-finite elements</h3>
+   *
+   * In the case of DoFHandlers with hp-capabilities, nothing defines which of
+   * the finite elements that are part of the hp::FECollection associated with
+   * the DoFHandler, should be considered on cells that are not active (i.e.,
+   * that have children). This is because degrees of freedom are only allocated
+   * for active cells and, in fact, it is not allowed to set an active FE index
+   * on non-active cells using DoFAccessor::set_active_fe_index().
+   *
+   * It is, thus, not entirely natural what should happen if, for example, a few
+   * cells are coarsened away. This class then implements the following
+   * algorithm:
+   * - If a cell is refined, then the values of the solution vector(s) are
+   *   interpolated before refinement on the to-be-refined cell from the space
+   * of the active finite element to the one of the future finite element. These
+   *   values are then distributed on the finite element spaces of the children
+   *   post-refinement. This may lose information if, for example, the old cell
+   *   used a Q2 space and the children use Q1 spaces, or the information may be
+   *   prolonged if the mother cell used a Q1 space and the children are Q2s.
+   * - If cells are to be coarsened, then the values from the child cells are
+   *   interpolated to the mother cell using the largest of the child cell
+   * future finite element spaces, which will be identified as the least
+   * dominant element following the FiniteElementDomination logic (consult
+   *   hp::FECollection::find_dominated_fe_extended() for more information). For
+   *   example, if the children of a cell use Q1, Q2 and Q3 spaces, then the
+   *   values from the children are interpolated into a Q3 space on the mother
+   *   cell. After refinement, this Q3 function on the mother cell is then
+   *   interpolated into the space the user has selected for this cell (which
+   * may be different from Q3, in this example, if the user has set the active
+   * FE index for a different space post-refinement and before calling
+   *   DoFHandler::distribute_dofs()).
+   *
+   * @note In the context of hp-refinement, if cells are coarsened or the
+   * polynomial degree is lowered on some cells, then the old finite element
+   * space is not a subspace of the new space and you may run into the same
+   * situation as discussed above with hanging nodes. You may want to consider
+   * calling AffineConstraints::distribute() on the vector obtained by
+   * transferring the solution.
+   *
+   * @ingroup numerics
+   *
+   * @deprecated Use dealii::SolutionTransfer instead.
    */
-  std::vector<std::vector<Vector<typename VectorType::value_type>>>
-    dof_values_on_cell;
-};
+  template <int dim, typename VectorType = Vector<double>, int spacedim = dim>
+  class SolutionTransfer
+  {
+  public:
+    /**
+     * Constructor, takes the current DoFHandler as argument.
+     */
+    SolutionTransfer(const DoFHandler<dim, spacedim> &dof);
+
+    /**
+     * Destructor
+     */
+    ~SolutionTransfer();
+
+    /**
+     * Reinit this class to the state that it has directly after calling the
+     * Constructor
+     */
+    void
+    clear();
+
+    /**
+     * Prepares the @p SolutionTransfer for pure refinement. It stores the dof
+     * indices of each cell. After calling this function only calling the @p
+     * refine_interpolate functions is allowed.
+     */
+    void
+    prepare_for_pure_refinement();
+
+    /**
+     * Prepares the @p SolutionTransfer for coarsening and refinement. It stores
+     * the dof indices of each cell and stores the dof values of the vectors in
+     * @p all_in in each cell that'll be coarsened. @p all_in includes all
+     * vectors that are to be interpolated onto the new (refined and/or
+     * coarsenend) grid.
+     */
+    void
+    prepare_for_coarsening_and_refinement(
+      const std::vector<VectorType> &all_in);
 
+    /**
+     * Same as previous function but for only one discrete function to be
+     * interpolated.
+     */
+    void
+    prepare_for_coarsening_and_refinement(const VectorType &in);
+
+    /**
+     * This function interpolates the discrete function @p in, which is a vector
+     * on the grid before the refinement, to the function @p out which then is a
+     * vector on the refined grid. It assumes the vectors having the right sizes
+     * (i.e. <tt>in.size()==n_dofs_old</tt>,
+     * <tt>out.size()==n_dofs_refined</tt>)
+     *
+     * Calling this function is allowed only if @p prepare_for_pure_refinement
+     * is called and the refinement is executed before. Multiple calling of this
+     * function is allowed. e.g. for interpolating several functions.
+     */
+    void
+    refine_interpolate(const VectorType &in, VectorType &out) const;
+
+    /**
+     * This function interpolates the discrete functions that are stored in @p
+     * all_in onto the refined and/or coarsenend grid. It assumes the vectors in
+     * @p all_in denote the same vectors as in @p all_in as parameter of
+     * <tt>prepare_for_refinement_and_coarsening(all_in)</tt>. However, there is
+     * no way of verifying this internally, so be careful here.
+     *
+     * Calling this function is allowed only if first
+     * Triangulation::prepare_coarsening_and_refinement, second @p
+     * SolutionTransfer::prepare_for_coarsening_and_refinement, an then third
+     * Triangulation::execute_coarsening_and_refinement are called before.
+     * Multiple calling of this function is NOT allowed. Interpolating several
+     * functions can be performed in one step.
+     *
+     * The number of output vectors is assumed to be the same as the number of
+     * input vectors. Also, the sizes of the output vectors are assumed to be of
+     * the right size (@p n_dofs_refined). Otherwise an assertion will be
+     * thrown.
+     */
+    void
+    interpolate(const std::vector<VectorType> &all_in,
+                std::vector<VectorType>       &all_out) const;
+
+    /**
+     * Same as the previous function. It interpolates only one function. It
+     * assumes the vectors having the right sizes (i.e.
+     * <tt>in.size()==n_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
+     *
+     * Multiple calling of this function is NOT allowed. Interpolating several
+     * functions can be performed in one step by using <tt>interpolate (all_in,
+     * all_out)</tt>
+     */
+    void
+    interpolate(const VectorType &in, VectorType &out) const;
+
+    /**
+     * Determine an estimate for the memory consumption (in bytes) of this
+     * object.
+     */
+    std::size_t
+    memory_consumption() const;
+
+    /**
+     * Exception
+     */
+    DeclExceptionMsg(
+      ExcNotPrepared,
+      "You are attempting an operation for which this object is "
+      "not prepared. This may be because you either did not call "
+      "one of the prepare_*() functions at all, or because you "
+      "called the wrong one for the operation you are currently "
+      "attempting.");
+
+    /**
+     * Exception
+     */
+    DeclExceptionMsg(
+      ExcAlreadyPrepForRef,
+      "You are attempting to call one of the prepare_*() functions "
+      "of this object to prepare it for an operation for which it "
+      "is already prepared. Specifically, the object was "
+      "previously prepared for pure refinement.");
+
+    /**
+     * Exception
+     */
+    DeclExceptionMsg(
+      ExcAlreadyPrepForCoarseAndRef,
+      "You are attempting to call one of the prepare_*() functions "
+      "of this object to prepare it for an operation for which it "
+      "is already prepared. Specifically, the object was "
+      "previously prepared for both coarsening and refinement.");
+
+  private:
+    /**
+     * Pointer to the degree of freedom handler to work with.
+     */
+    SmartPointer<const DoFHandler<dim, spacedim>,
+                 SolutionTransfer<dim, VectorType, spacedim>>
+      dof_handler;
+
+    /**
+     * Stores the number of DoFs before the refinement and/or coarsening.
+     */
+    types::global_dof_index n_dofs_old;
+
+    /**
+     * Declaration of @p PreparationState that denotes the three possible states
+     * of the @p SolutionTransfer: being prepared for 'pure refinement',
+     * prepared for 'coarsening and refinement' or not prepared.
+     */
+    enum PreparationState
+    {
+      /**
+       * The SolutionTransfer is not yet prepared.
+       */
+      none,
+      /**
+       * The SolutionTransfer is prepared for purely refinement.
+       */
+      pure_refinement,
+      /**
+       * The SolutionTransfer is prepared for coarsening and refinement.
+       */
+      coarsening_and_refinement
+    };
+
+    /**
+     * Definition of the respective variable.
+     */
+    PreparationState prepared_for;
+
+
+    /**
+     * Is used for @p prepare_for_refining (of course also for @p
+     * repare_for_refining_and_coarsening) and stores all dof indices of the
+     * cells that'll be refined
+     */
+    std::vector<std::vector<types::global_dof_index>> indices_on_cell;
+
+    /**
+     * All cell data (the dof indices and the dof values) should be accessible
+     * from each cell. As each cell has got only one @p user_pointer, multiple
+     * pointers to the data need to be packetized in a structure. Note that in
+     * our case on each cell either the <tt>vector<unsigned int> indices</tt>
+     * (if the cell will be refined) or the <tt>vector<double> dof_values</tt>
+     * (if the children of this cell will be deleted) is needed, hence one @p
+     * user_pointer should be sufficient, but to allow some error checks and to
+     * preserve the user from making user errors the @p user_pointer will be
+     * 'multiplied' by this structure.
+     */
+    struct Pointerstruct
+    {
+      Pointerstruct()
+        : indices_ptr(nullptr)
+        , dof_values_ptr(nullptr)
+        , active_fe_index(0)
+      {}
+      Pointerstruct(std::vector<types::global_dof_index> *indices_ptr_in,
+                    const unsigned int active_fe_index_in = 0)
+        : indices_ptr(indices_ptr_in)
+        , dof_values_ptr(nullptr)
+        , active_fe_index(active_fe_index_in)
+      {}
+      Pointerstruct(
+        std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr_in,
+        const unsigned int active_fe_index_in = 0)
+        : indices_ptr(nullptr)
+        , dof_values_ptr(dof_values_ptr_in)
+        , active_fe_index(active_fe_index_in)
+      {}
+      std::size_t
+      memory_consumption() const;
+
+      std::vector<types::global_dof_index>                 *indices_ptr;
+      std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr;
+      unsigned int                                          active_fe_index;
+    };
+
+    /**
+     * Map mapping from level and index of cell to the @p Pointerstructs (cf.
+     * there). This map makes it possible to keep all the information needed to
+     * transfer the solution inside this object rather than using user pointers
+     * of the Triangulation for this purpose.
+     */
+    std::map<std::pair<unsigned int, unsigned int>, Pointerstruct> cell_map;
+
+    /**
+     * Is used for @p prepare_for_refining_and_coarsening The interpolated dof
+     * values of all cells that'll be coarsened will be stored in this vector.
+     */
+    std::vector<std::vector<Vector<typename VectorType::value_type>>>
+      dof_values_on_cell;
+  };
+} // namespace Legacy
 
 DEAL_II_NAMESPACE_CLOSE
 
index 252838c5a820c69d0f3cedcfc71e52f90903ea16..5eb431c6d38bc13544b92f418f310cca38d5d2ef 100644 (file)
@@ -18,7 +18,6 @@ set(_unity_include_src
   cell_data_transfer.cc
   fully_distributed_tria.cc
   repartitioning_policy_tools.cc
-  solution_transfer.cc
   tria.cc
   tria_base.cc
   shared_tria.cc
@@ -45,7 +44,6 @@ set(_inst
   field_transfer.inst.in
   fully_distributed_tria.inst.in
   repartitioning_policy_tools.inst.in
-  solution_transfer.inst.in
   tria.inst.in
   shared_tria.inst.in
   tria_base.inst.in
index a7f3fae58ff3fbed6a34a5f2285ec43c3a204269..0d96cf982a2d3b8ed24626ab3d433289b7c6860d 100644 (file)
@@ -491,11 +491,12 @@ namespace parallel
           std::ofstream f(fname);
           f << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_global_active_cells"
             << std::endl
-            << 4 << " "
-            << Utilities::MPI::n_mpi_processes(this->mpi_communicator) << " "
-            << this->cell_attached_data.pack_callbacks_fixed.size() << " "
-            << this->cell_attached_data.pack_callbacks_variable.size() << " "
-            << this->n_global_active_cells() << std::endl;
+            << ::dealii::internal::CellAttachedDataSerializer<dim, spacedim>::
+                 version_number
+            << " " << Utilities::MPI::n_mpi_processes(this->mpi_communicator)
+            << " " << this->cell_attached_data.pack_callbacks_fixed.size()
+            << " " << this->cell_attached_data.pack_callbacks_variable.size()
+            << " " << this->n_global_active_cells() << std::endl;
         }
 
       // Save cell attached data.
@@ -606,12 +607,15 @@ namespace parallel
         std::ifstream f(fname);
         AssertThrow(f.fail() == false, ExcIO());
         std::string firstline;
-        getline(f, firstline); // skip first line
+        getline(f, firstline);
         f >> version >> numcpus >> attached_count_fixed >>
           attached_count_variable >> n_global_active_cells;
       }
 
-      AssertThrow(version == 4,
+      const auto expected_version = ::dealii::internal::
+        CellAttachedDataSerializer<dim, spacedim>::version_number;
+
+      AssertThrow(version == expected_version,
                   ExcMessage("Incompatible version found in .info file."));
 
       // Load description and construct the triangulation.
diff --git a/source/distributed/solution_transfer.cc b/source/distributed/solution_transfer.cc
deleted file mode 100644 (file)
index 250c000..0000000
+++ /dev/null
@@ -1,482 +0,0 @@
-// ------------------------------------------------------------------------
-//
-// SPDX-License-Identifier: LGPL-2.1-or-later
-// Copyright (C) 2010 - 2024 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// Part of the source code is dual licensed under Apache-2.0 WITH
-// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
-// governing the source code and code contributions can be found in
-// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
-//
-// ------------------------------------------------------------------------
-
-
-#include <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_P4EST
-
-#  include <deal.II/distributed/solution_transfer.h>
-#  include <deal.II/distributed/tria.h>
-
-#  include <deal.II/dofs/dof_accessor.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
-{
-  /**
-   * Optimized pack function for values assigned on degrees of freedom.
-   *
-   * Given that the elements of @p dof_values are stored in consecutive
-   * locations, we can just memcpy them. Since floating point values don't
-   * compress well, we also waive the compression that the default
-   * Utilities::pack() and Utilities::unpack() functions offer.
-   */
-  template <typename value_type>
-  std::vector<char>
-  pack_dof_values(std::vector<Vector<value_type>> &dof_values,
-                  const unsigned int               dofs_per_cell)
-  {
-    for (const auto &values : dof_values)
-      {
-        AssertDimension(values.size(), dofs_per_cell);
-        (void)values;
-      }
-
-    const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
-
-    std::vector<char> buffer(dof_values.size() * bytes_per_entry);
-    for (unsigned int i = 0; i < dof_values.size(); ++i)
-      std::memcpy(&buffer[i * bytes_per_entry],
-                  &dof_values[i](0),
-                  bytes_per_entry);
-
-    return buffer;
-  }
-
-
-
-  /**
-   * Optimized unpack function for values assigned on degrees of freedom.
-   */
-  template <typename value_type>
-  std::vector<Vector<value_type>>
-  unpack_dof_values(
-    const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
-    const unsigned int dofs_per_cell)
-  {
-    const std::size_t  bytes_per_entry = sizeof(value_type) * dofs_per_cell;
-    const unsigned int n_elements      = data_range.size() / bytes_per_entry;
-
-    Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
-
-    std::vector<Vector<value_type>> unpacked_data;
-    unpacked_data.reserve(n_elements);
-    for (unsigned int i = 0; i < n_elements; ++i)
-      {
-        Vector<value_type> dof_values(dofs_per_cell);
-        std::memcpy(&dof_values(0),
-                    &(*std::next(data_range.begin(), i * bytes_per_entry)),
-                    bytes_per_entry);
-        unpacked_data.emplace_back(std::move(dof_values));
-      }
-
-    return unpacked_data;
-  }
-} // namespace
-
-
-
-namespace parallel
-{
-  namespace distributed
-  {
-    template <int dim, typename VectorType, int spacedim>
-    SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
-      const DoFHandler<dim, spacedim> &dof,
-      const bool                       average_values)
-      : dof_handler(&dof, typeid(*this).name())
-      , average_values(average_values)
-      , handle(numbers::invalid_unsigned_int)
-    {
-      Assert(
-        (dynamic_cast<
-           const parallel::DistributedTriangulationBase<dim, spacedim> *>(
-           &dof_handler->get_triangulation()) != nullptr),
-        ExcMessage(
-          "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::
-      prepare_for_coarsening_and_refinement(
-        const std::vector<const VectorType *> &all_in)
-    {
-      for (unsigned int i = 0; i < all_in.size(); ++i)
-        Assert(all_in[i]->size() == dof_handler->n_dofs(),
-               ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
-
-      input_vectors = all_in;
-      register_data_attach();
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::register_data_attach()
-    {
-      // TODO: casting away constness is bad
-      parallel::DistributedTriangulationBase<dim, spacedim> *tria =
-        (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
-          const_cast<dealii::Triangulation<dim, spacedim> *>(
-            &dof_handler->get_triangulation())));
-      Assert(tria != nullptr, ExcInternalError());
-
-      Assert(handle == numbers::invalid_unsigned_int,
-             ExcMessage("You can only add one solution per "
-                        "SolutionTransfer object."));
-
-      handle = tria->register_data_attach(
-        [this](
-          const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
-          const CellStatus                                            status) {
-          return this->pack_callback(cell_, status);
-        },
-        /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::
-      prepare_for_coarsening_and_refinement(const VectorType &in)
-    {
-      std::vector<const VectorType *> all_in(1, &in);
-      prepare_for_coarsening_and_refinement(all_in);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
-      const VectorType &in)
-    {
-      std::vector<const VectorType *> all_in(1, &in);
-      prepare_for_serialization(all_in);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
-      const std::vector<const VectorType *> &all_in)
-    {
-      prepare_for_coarsening_and_refinement(all_in);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::deserialize(VectorType &in)
-    {
-      std::vector<VectorType *> all_in(1, &in);
-      deserialize(all_in);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::deserialize(
-      std::vector<VectorType *> &all_in)
-    {
-      register_data_attach();
-
-      // this makes interpolate() happy
-      input_vectors.resize(all_in.size());
-
-      interpolate(all_in);
-    }
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::interpolate(
-      std::vector<VectorType *> &all_out)
-    {
-      Assert(input_vectors.size() == all_out.size(),
-             ExcDimensionMismatch(input_vectors.size(), all_out.size()));
-      for (unsigned int i = 0; i < all_out.size(); ++i)
-        Assert(all_out[i]->size() == dof_handler->n_dofs(),
-               ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
-
-      // TODO: casting away constness is bad
-      parallel::DistributedTriangulationBase<dim, spacedim> *tria =
-        (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
-          const_cast<dealii::Triangulation<dim, spacedim> *>(
-            &dof_handler->get_triangulation())));
-      Assert(tria != nullptr, ExcInternalError());
-
-      if (average_values)
-        for (auto *const vec : all_out)
-          *vec = 0.0;
-
-      VectorType valence;
-
-      // initialize valence vector only if we need to average
-      if (average_values)
-        valence.reinit(*all_out[0]);
-
-      tria->notify_ready_to_unpack(
-        handle,
-        [this, &all_out, &valence](
-          const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
-          const CellStatus                                            status,
-          const boost::iterator_range<std::vector<char>::const_iterator>
-            &data_range) {
-          this->unpack_callback(cell_, status, data_range, all_out, valence);
-        });
-
-      if (average_values)
-        {
-          // finalize valence: compress and invert
-          using Number = typename VectorType::value_type;
-          valence.compress(VectorOperation::add);
-          for (const auto i : valence.locally_owned_elements())
-            valence[i] = (static_cast<Number>(valence[i]) == Number() ?
-                            Number() :
-                            (Number(1.0) / static_cast<Number>(valence[i])));
-          valence.compress(VectorOperation::insert);
-
-          for (auto *const vec : all_out)
-            {
-              // compress and weight with valence
-              vec->compress(VectorOperation::add);
-              vec->scale(valence);
-            }
-        }
-      else
-        {
-          for (auto *const vec : all_out)
-            vec->compress(VectorOperation::insert);
-        }
-
-      input_vectors.clear();
-      handle = numbers::invalid_unsigned_int;
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::interpolate(VectorType &out)
-    {
-      std::vector<VectorType *> all_out(1, &out);
-      interpolate(all_out);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    std::vector<char>
-    SolutionTransfer<dim, VectorType, spacedim>::pack_callback(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
-      const CellStatus                                            status)
-    {
-      typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
-                                                             dof_handler);
-
-      // create buffer for each individual object
-      std::vector<::dealii::Vector<typename VectorType::value_type>> dof_values(
-        input_vectors.size());
-
-      unsigned int fe_index = 0;
-      if (dof_handler->has_hp_capabilities())
-        {
-          switch (status)
-            {
-              case CellStatus::cell_will_persist:
-              case CellStatus::cell_will_be_refined:
-                {
-                  fe_index = cell->future_fe_index();
-                  break;
-                }
-
-              case CellStatus::children_will_be_coarsened:
-                {
-                  // In case of coarsening, we need to find a suitable FE index
-                  // for the parent cell. We choose the 'least dominant fe'
-                  // on all children from the associated FECollection.
-#  ifdef DEBUG
-                  for (const auto &child : cell->child_iterators())
-                    Assert(child->is_active() && child->coarsen_flag_set(),
-                           typename dealii::Triangulation<
-                             dim>::ExcInconsistentCoarseningFlags());
-#  endif
-
-                  fe_index = dealii::internal::hp::DoFHandlerImplementation::
-                    dominated_future_fe_on_children<dim, spacedim>(cell);
-                  break;
-                }
-
-              default:
-                DEAL_II_ASSERT_UNREACHABLE();
-                break;
-            }
-        }
-
-      const unsigned int dofs_per_cell =
-        dof_handler->get_fe(fe_index).n_dofs_per_cell();
-
-      if (dofs_per_cell == 0)
-        return std::vector<char>(); // nothing to do for FE_Nothing
-
-      auto it_input  = input_vectors.cbegin();
-      auto it_output = dof_values.begin();
-      for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
-        {
-          it_output->reinit(dofs_per_cell);
-          cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
-        }
-
-      return pack_dof_values<typename VectorType::value_type>(dof_values,
-                                                              dofs_per_cell);
-    }
-
-
-
-    template <int dim, typename VectorType, int spacedim>
-    void
-    SolutionTransfer<dim, VectorType, spacedim>::unpack_callback(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
-      const CellStatus                                            status,
-      const boost::iterator_range<std::vector<char>::const_iterator>
-                                &data_range,
-      std::vector<VectorType *> &all_out,
-      VectorType                &valence)
-    {
-      typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
-                                                             dof_handler);
-
-      unsigned int fe_index = 0;
-      if (dof_handler->has_hp_capabilities())
-        {
-          switch (status)
-            {
-              case CellStatus::cell_will_persist:
-              case CellStatus::children_will_be_coarsened:
-                {
-                  fe_index = cell->active_fe_index();
-                  break;
-                }
-
-              case CellStatus::cell_will_be_refined:
-                {
-                  // After refinement, this particular cell is no longer active,
-                  // and its children have inherited its FE index. However, to
-                  // unpack the data on the old cell, we need to recover its FE
-                  // index from one of the children. Just to be sure, we also
-                  // check if all children have the same FE index.
-                  fe_index = cell->child(0)->active_fe_index();
-                  for (unsigned int child_index = 1;
-                       child_index < cell->n_children();
-                       ++child_index)
-                    Assert(cell->child(child_index)->active_fe_index() ==
-                             fe_index,
-                           ExcInternalError());
-                  break;
-                }
-
-              default:
-                DEAL_II_ASSERT_UNREACHABLE();
-                break;
-            }
-        }
-
-      const unsigned int dofs_per_cell =
-        dof_handler->get_fe(fe_index).n_dofs_per_cell();
-
-      if (dofs_per_cell == 0)
-        return; // nothing to do for FE_Nothing
-
-      const std::vector<::dealii::Vector<typename VectorType::value_type>>
-        dof_values =
-          unpack_dof_values<typename VectorType::value_type>(data_range,
-                                                             dofs_per_cell);
-
-      // check if sizes match
-      Assert(dof_values.size() == all_out.size(), ExcInternalError());
-
-      // check if we have enough dofs provided by the FE object
-      // to interpolate the transferred data correctly
-      for (auto it_dof_values = dof_values.begin();
-           it_dof_values != dof_values.end();
-           ++it_dof_values)
-        Assert(
-          dofs_per_cell == it_dof_values->size(),
-          ExcMessage(
-            "The transferred data was packed with a different number of dofs than the "
-            "currently registered FE object assigned to the DoFHandler has."));
-
-      // distribute data for each registered vector on mesh
-      auto it_input  = dof_values.cbegin();
-      auto it_output = all_out.begin();
-      for (; it_input != dof_values.cend(); ++it_input, ++it_output)
-        if (average_values)
-          cell->distribute_local_to_global_by_interpolation(*it_input,
-                                                            *(*it_output),
-                                                            fe_index);
-        else
-          cell->set_dof_values_by_interpolation(*it_input,
-                                                *(*it_output),
-                                                fe_index,
-                                                true);
-
-      if (average_values)
-        {
-          // compute valence vector if averaging should be performed
-          Vector<typename VectorType::value_type> ones(dofs_per_cell);
-          ones = 1.0;
-          cell->distribute_local_to_global_by_interpolation(ones,
-                                                            valence,
-                                                            fe_index);
-        }
-    }
-  } // namespace distributed
-} // namespace parallel
-
-
-// explicit instantiations
-#  include "solution_transfer.inst"
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
diff --git a/source/distributed/solution_transfer.inst.in b/source/distributed/solution_transfer.inst.in
deleted file mode 100644 (file)
index 931ddac..0000000
+++ /dev/null
@@ -1,66 +0,0 @@
-// ------------------------------------------------------------------------
-//
-// SPDX-License-Identifier: LGPL-2.1-or-later
-// Copyright (C) 2010 - 2021 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// Part of the source code is dual licensed under Apache-2.0 WITH
-// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
-// governing the source code and code contributions can be found in
-// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
-//
-// ------------------------------------------------------------------------
-
-
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
-  {
-    namespace parallel
-    \{
-      namespace distributed
-      \{
-#if deal_II_dimension <= deal_II_space_dimension
-        template class SolutionTransfer<deal_II_dimension,
-                                        ::dealii::Vector<double>,
-                                        deal_II_space_dimension>;
-        template class SolutionTransfer<
-          deal_II_dimension,
-          ::dealii::LinearAlgebra::distributed::Vector<double>,
-          deal_II_space_dimension>;
-        template class SolutionTransfer<
-          deal_II_dimension,
-          ::dealii::LinearAlgebra::distributed::Vector<float>,
-          deal_II_space_dimension>;
-        template class SolutionTransfer<
-          deal_II_dimension,
-          ::dealii::LinearAlgebra::distributed::BlockVector<double>,
-          deal_II_space_dimension>;
-        template class SolutionTransfer<
-          deal_II_dimension,
-          ::dealii::LinearAlgebra::distributed::BlockVector<float>,
-          deal_II_space_dimension>;
-
-
-#  ifdef DEAL_II_WITH_PETSC
-        template class SolutionTransfer<deal_II_dimension,
-                                        PETScWrappers::MPI::Vector,
-                                        deal_II_space_dimension>;
-        template class SolutionTransfer<deal_II_dimension,
-                                        PETScWrappers::MPI::BlockVector,
-                                        deal_II_space_dimension>;
-#  endif
-
-#  ifdef DEAL_II_WITH_TRILINOS
-        template class SolutionTransfer<deal_II_dimension,
-                                        TrilinosWrappers::MPI::Vector,
-                                        deal_II_space_dimension>;
-        template class SolutionTransfer<deal_II_dimension,
-                                        TrilinosWrappers::MPI::BlockVector,
-                                        deal_II_space_dimension>;
-#  endif
-
-#endif
-      \}
-    \}
-  }
index b1e5adec704daa9ed9ef9f8361234df038dd1915..57e891dbc525cc61690f0080a80b65d5e9878f02 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/tria_base.h>
 
 #include <deal.II/grid/connectivity.h>
 #include <deal.II/grid/grid_tools.h>
@@ -213,9 +214,7 @@ namespace internal
                 // can only tolerate one level of coarsening at a time, so
                 // check that the children are all active
                 Assert(dealii_cell->is_active() == false, ExcInternalError());
-                for (unsigned int c = 0;
-                     c < GeometryInfo<dim>::max_children_per_cell;
-                     ++c)
+                for (unsigned int c = 0; c < dealii_cell->n_children(); ++c)
                   Assert(dealii_cell->child(c)->is_active(),
                          ExcInternalError());
                 break;
@@ -1080,7 +1079,6 @@ namespace internal
 
           sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
                                         (variable_size_data_stored ? 1 : 0));
-
           // Read header data.
           file.read(reinterpret_cast<char *>(sizes_fixed_cumulative.data()),
                     sizes_fixed_cumulative.size() * sizeof(unsigned int));
@@ -13700,14 +13698,19 @@ DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
 void Triangulation<dim, spacedim>::save(const std::string &file_basename) const
 {
   // Save triangulation information.
-  std::ofstream                 ofs(file_basename + "_triangulation.data");
-  boost::archive::text_oarchive oa(ofs, boost::archive::no_header);
-  save(oa, 0);
+  {
+    std::ofstream ofs_tria(file_basename + "_triangulation.data");
+    AssertThrow(ofs_tria.fail() == false, ExcIO());
+
+    boost::archive::text_oarchive oa(ofs_tria, boost::archive::no_header);
+    save(oa,
+         internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
+  }
 
   // Save attached data.
   {
-    std::ofstream ifs(file_basename + ".info");
-    ifs
+    std::ofstream ofs_info(file_basename + ".info");
+    ofs_info
       << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_active_cells"
       << std::endl
       << internal::CellAttachedDataSerializer<dim, spacedim>::version_number
@@ -13726,21 +13729,24 @@ DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
 void Triangulation<dim, spacedim>::load(const std::string &file_basename)
 {
   // Load triangulation information.
-  std::ifstream ifs(file_basename + "_triangulation.data");
-  AssertThrow(ifs.fail() == false, ExcIO());
+  {
+    std::ifstream ifs_tria(file_basename + "_triangulation.data");
+    AssertThrow(ifs_tria.fail() == false, ExcIO());
 
-  boost::archive::text_iarchive ia(ifs, boost::archive::no_header);
-  load(ia, 0);
+    boost::archive::text_iarchive ia(ifs_tria, boost::archive::no_header);
+    load(ia,
+         internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
+  }
 
   // Load attached data.
   unsigned int version, numcpus, attached_count_fixed, attached_count_variable,
     n_global_active_cells;
   {
-    std::ifstream ifs(std::string(file_basename) + ".info");
-    AssertThrow(ifs.fail() == false, ExcIO());
+    std::ifstream ifs_info(std::string(file_basename) + ".info");
+    AssertThrow(ifs_info.fail() == false, ExcIO());
     std::string firstline;
-    getline(ifs, firstline);
-    ifs >> version >> numcpus >> attached_count_fixed >>
+    std::getline(ifs_info, firstline);
+    ifs_info >> version >> numcpus >> attached_count_fixed >>
       attached_count_variable >> n_global_active_cells;
   }
 
@@ -15652,6 +15658,161 @@ void Triangulation<dim, spacedim>::update_cell_relations()
 
 
 
+template <int dim, int spacedim>
+DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
+void Triangulation<dim, spacedim>::pack_data_serial()
+{
+  if (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
+        this))
+    return;
+
+  std::vector<CellId> active_cell_old;
+
+  // pack data before triangulation gets updated
+  if (this->cell_attached_data.n_attached_data_sets > 0)
+    {
+      // store old active cells to determine cell status after
+      // coarsening/refinement
+      active_cell_old.reserve(this->n_active_cells());
+
+      for (const auto &cell : this->active_cell_iterators())
+        {
+          const bool children_will_be_coarsened =
+            (cell->level() > 0) && (cell->coarsen_flag_set());
+
+          if (children_will_be_coarsened == false)
+            active_cell_old.emplace_back(cell->id());
+          else
+            {
+              if (cell->parent()->child(0) == cell)
+                active_cell_old.emplace_back(cell->parent()->id());
+            }
+        }
+
+      // update cell relations
+      this->local_cell_relations.clear();
+      this->local_cell_relations.reserve(this->n_global_active_cells());
+
+      std::vector<
+        std::pair<unsigned int,
+                  typename internal::CellAttachedDataSerializer<dim, spacedim>::
+                    cell_relation_t>>
+        local_cell_relations_tmp;
+
+      for (const auto &cell : this->active_cell_iterators())
+        {
+          if (std::find(active_cell_old.begin(),
+                        active_cell_old.end(),
+                        cell->id()) != active_cell_old.end())
+            {
+              const unsigned int index =
+                std::distance(active_cell_old.begin(),
+                              std::find(active_cell_old.begin(),
+                                        active_cell_old.end(),
+                                        cell->id()));
+
+              ::dealii::CellStatus status =
+                cell->refine_flag_set() ?
+                  ::dealii::CellStatus::cell_will_be_refined :
+                  ::dealii::CellStatus::cell_will_persist;
+
+              local_cell_relations_tmp.emplace_back(
+                index,
+                typename internal::CellAttachedDataSerializer<dim, spacedim>::
+                  cell_relation_t{cell, status});
+            }
+          else if (cell->level() > 0 &&
+                   std::find(active_cell_old.begin(),
+                             active_cell_old.end(),
+                             cell->parent()->id()) != active_cell_old.end())
+            {
+              const unsigned int index =
+                std::distance(active_cell_old.begin(),
+                              std::find(active_cell_old.begin(),
+                                        active_cell_old.end(),
+                                        cell->parent()->id()));
+
+              ::dealii::CellStatus status;
+
+              if (cell->parent()->child_iterator_to_index(cell) == 0)
+                status = ::dealii::CellStatus::children_will_be_coarsened;
+              else
+                status = ::dealii::CellStatus::cell_invalid;
+
+              local_cell_relations_tmp.emplace_back(
+                index,
+                typename internal::CellAttachedDataSerializer<dim, spacedim>::
+                  cell_relation_t{cell->parent(), status});
+            }
+          else
+            {
+              AssertThrow(false, ExcNotImplemented());
+            }
+        }
+
+      std::stable_sort(local_cell_relations_tmp.begin(),
+                       local_cell_relations_tmp.end(),
+                       [](const auto &a, const auto &b) {
+                         return a.first < b.first;
+                       });
+
+      for (const auto &tmp : local_cell_relations_tmp)
+        this->local_cell_relations.emplace_back(tmp.second);
+
+      // pack data
+      this->data_serializer.pack_data(
+        this->local_cell_relations,
+        this->cell_attached_data.pack_callbacks_fixed,
+        this->cell_attached_data.pack_callbacks_variable,
+        this->get_communicator());
+
+      // dummy copy of data
+      this->data_serializer.dest_data_fixed =
+        this->data_serializer.src_data_fixed;
+      this->data_serializer.dest_data_variable =
+        this->data_serializer.src_data_variable;
+      this->data_serializer.dest_sizes_variable =
+        this->data_serializer.src_sizes_variable;
+    }
+}
+
+
+
+template <int dim, int spacedim>
+DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
+void Triangulation<dim, spacedim>::unpack_data_serial()
+{
+  if (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
+        this))
+    return;
+
+  // transfer data after triangulation got updated
+  if (this->cell_attached_data.n_attached_data_sets > 0)
+    {
+      std::vector<typename internal::CellAttachedDataSerializer<dim, spacedim>::
+                    cell_relation_t>
+        temp;
+
+      for (const auto &cell : local_cell_relations)
+        {
+          if (cell.first->has_children())
+            {
+              Assert(cell.second == ::dealii::CellStatus::cell_will_be_refined,
+                     ExcInternalError());
+
+              temp.emplace_back(cell.first->child(0),
+                                ::dealii::CellStatus::cell_will_be_refined);
+            }
+          else
+            temp.push_back(cell);
+        }
+
+      this->local_cell_relations = temp;
+    }
+}
+
+
+
 template <int dim, int spacedim>
 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
 void Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
@@ -15672,10 +15833,17 @@ void Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
   // Inform all listeners about beginning of refinement.
   signals.pre_refinement();
 
+  this->pack_data_serial();
+
   execute_coarsening();
 
   const DistortedCellList cells_with_distorted_children = execute_refinement();
 
+  // We need to update the cell relations in order to be able to
+  // deserialize data. Later on, update_cell_relations is called to mark all
+  // active cells with the cell_will_persist status.
+  this->unpack_data_serial();
+
   reset_cell_vertex_indices_cache();
 
   // verify a case with which we have had
@@ -16072,6 +16240,15 @@ void Triangulation<dim, spacedim>::save_attached_data(
   // cast away constness
   auto tria = const_cast<Triangulation<dim, spacedim> *>(this);
 
+  // each cell should have been flagged `CellStatus::cell_will_persist`
+  for (const auto &cell_rel : this->local_cell_relations)
+    {
+      (void)cell_rel;
+      Assert((cell_rel.second == // cell_status
+              dealii::CellStatus::cell_will_persist),
+             ExcInternalError());
+    }
+
   if (this->cell_attached_data.n_attached_data_sets > 0)
     {
       // pack attached data first
index 81c87c503e26d43c361735091c74f57271bf8815..3f5e26856fe64936133d8fd839a96436bf1610a2 100644 (file)
@@ -12,6 +12,8 @@
 //
 // ------------------------------------------------------------------------
 
+#include <deal.II/base/config.h>
+
 #include <deal.II/base/memory_consumption.h>
 
 #include <deal.II/distributed/shared_tria.h>
@@ -19,6 +21,7 @@
 
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
 
 #include <deal.II/fe/fe.h>
 
 
 #include <deal.II/numerics/solution_transfer.h>
 
+#include <functional>
+#include <numeric>
+
 DEAL_II_NAMESPACE_OPEN
 
+
+namespace
+{
+  /**
+   * Optimized pack function for values assigned on degrees of freedom.
+   *
+   * Given that the elements of @p dof_values are stored in consecutive
+   * locations, we can just memcpy them. Since floating point values don't
+   * compress well, we also waive the compression that the default
+   * Utilities::pack() and Utilities::unpack() functions offer.
+   */
+  template <typename value_type>
+  std::vector<char>
+  pack_dof_values(std::vector<Vector<value_type>> &dof_values,
+                  const unsigned int               dofs_per_cell)
+  {
+    for (const auto &values : dof_values)
+      {
+        AssertDimension(values.size(), dofs_per_cell);
+        (void)values;
+      }
+
+    const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
+
+    std::vector<char> buffer(dof_values.size() * bytes_per_entry);
+    for (unsigned int i = 0; i < dof_values.size(); ++i)
+      std::memcpy(&buffer[i * bytes_per_entry],
+                  &dof_values[i](0),
+                  bytes_per_entry);
+
+    return buffer;
+  }
+
+
+
+  /**
+   * Optimized unpack function for values assigned on degrees of freedom.
+   */
+  template <typename value_type>
+  std::vector<Vector<value_type>>
+  unpack_dof_values(
+    const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+    const unsigned int dofs_per_cell)
+  {
+    const std::size_t  bytes_per_entry = sizeof(value_type) * dofs_per_cell;
+    const unsigned int n_elements      = data_range.size() / bytes_per_entry;
+
+    Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
+
+    std::vector<Vector<value_type>> unpacked_data;
+    unpacked_data.reserve(n_elements);
+    for (unsigned int i = 0; i < n_elements; ++i)
+      {
+        Vector<value_type> dof_values(dofs_per_cell);
+        std::memcpy(&dof_values(0),
+                    &(*std::next(data_range.begin(), i * bytes_per_entry)),
+                    bytes_per_entry);
+        unpacked_data.emplace_back(std::move(dof_values));
+      }
+
+    return unpacked_data;
+  }
+} // namespace
+
+
+
 template <int dim, typename VectorType, int spacedim>
 SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
-  const DoFHandler<dim, spacedim> &dof)
+  const DoFHandler<dim, spacedim> &dof,
+  const bool                       average_values)
   : dof_handler(&dof, typeid(*this).name())
-  , n_dofs_old(0)
-  , prepared_for(none)
+  , average_values(average_values)
+  , handle(numbers::invalid_unsigned_int)
+{}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+  prepare_for_coarsening_and_refinement(
+    const std::vector<const VectorType *> &all_in)
 {
-  Assert(
-    (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
-       &dof_handler->get_triangulation()) == nullptr),
-    ExcMessage("You are calling the dealii::SolutionTransfer class "
-               "with a DoFHandler that is built on a "
-               "parallel::distributed::Triangulation. This will not "
-               "work for parallel computations. You probably want to "
-               "use the parallel::distributed::SolutionTransfer class."));
+  const dealii::internal::parallel::shared::
+    TemporarilyRestoreSubdomainIds<dim, spacedim>
+      subdomain_modifier(dof_handler->get_triangulation());
+
+  for (unsigned int i = 0; i < all_in.size(); ++i)
+    Assert(all_in[i]->size() == dof_handler->n_dofs(),
+           ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
+
+  input_vectors = all_in;
+  register_data_attach();
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
-SolutionTransfer<dim, VectorType, spacedim>::~SolutionTransfer()
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
 {
-  clear();
+  std::vector<const VectorType *> temp(all_in.size());
+
+  for (std::size_t i = 0; i < temp.size(); ++i)
+    temp[i] = &(all_in[i]);
+
+  this->prepare_for_coarsening_and_refinement(temp);
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::clear()
+SolutionTransfer<dim, VectorType, spacedim>::register_data_attach()
 {
-  indices_on_cell.clear();
-  dof_values_on_cell.clear();
-  cell_map.clear();
+  // TODO: casting away constness is bad
+  auto tria = const_cast<dealii::Triangulation<dim, spacedim> *>(
+    &dof_handler->get_triangulation());
+  Assert(tria != nullptr, ExcInternalError());
+
+  Assert(handle == numbers::invalid_unsigned_int,
+         ExcMessage("You can only add one solution per "
+                    "SolutionTransfer object."));
+
+  handle = tria->register_data_attach(
+    [this](const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+           const CellStatus                                            status) {
+      return this->pack_callback(cell_, status);
+    },
+    /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
+}
 
-  prepared_for = none;
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+  prepare_for_coarsening_and_refinement(const VectorType &in)
+{
+  std::vector<const VectorType *> all_in(1, &in);
+  prepare_for_coarsening_and_refinement(all_in);
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::prepare_for_pure_refinement()
+SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
+  const VectorType &in)
 {
-  Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
-  Assert(prepared_for != coarsening_and_refinement,
-         ExcAlreadyPrepForCoarseAndRef());
-
-  clear();
-
-  // We need to access dof indices on the entire domain. For
-  // parallel::shared::Triangulations, ownership of cells might change. If they
-  // allow artificial cells, we need to restore the "true" cell owners
-  // temporarily.
-  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
-  // the current set of subdomain ids, set subdomain ids to the "true" owner of
-  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
-  // and later restore these flags when it is destroyed.
-  const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
-                                                                   spacedim>
-    subdomain_modifier(dof_handler->get_triangulation());
-
-  const unsigned int n_active_cells =
-    dof_handler->get_triangulation().n_active_cells();
-  n_dofs_old = dof_handler->n_dofs();
-
-  // efficient reallocation of indices_on_cell
-  std::vector<std::vector<types::global_dof_index>>(n_active_cells)
-    .swap(indices_on_cell);
-
-  for (const auto &cell : dof_handler->active_cell_iterators())
-    {
-      const unsigned int i = cell->active_cell_index();
-      indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
-      // on each cell store the indices of the
-      // dofs. after refining we get the values
-      // on the children by taking these
-      // indices, getting the respective values
-      // out of the data vectors and prolonging
-      // them to the children
-      cell->get_dof_indices(indices_on_cell[i]);
-      cell_map[std::make_pair(cell->level(), cell->index())] =
-        Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
-    }
-  prepared_for = pure_refinement;
+  std::vector<const VectorType *> all_in(1, &in);
+  prepare_for_serialization(all_in);
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
-  const VectorType &in,
-  VectorType       &out) const
+SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
+  const std::vector<const VectorType *> &all_in)
 {
-  Assert(prepared_for == pure_refinement, ExcNotPrepared());
-  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
-  Assert(out.size() == dof_handler->n_dofs(),
-         ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
-  Assert(&in != &out,
-         ExcMessage("Vectors cannot be used as input and output"
-                    " at the same time!"));
-
-  // We need to access dof indices on the entire domain. For
-  // parallel::shared::Triangulations, ownership of cells might change. If they
-  // allow artificial cells, we need to restore the "true" cell owners
-  // temporarily.
-  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
-  // the current set of subdomain ids, set subdomain ids to the "true" owner of
-  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
-  // and later restore these flags when it is destroyed.
-  const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
-                                                                   spacedim>
-    subdomain_modifier(dof_handler->get_triangulation());
-
-  Vector<typename VectorType::value_type> local_values(0);
-
-  typename std::map<std::pair<unsigned int, unsigned int>,
-                    Pointerstruct>::const_iterator pointerstruct,
-    cell_map_end = cell_map.end();
-
-  for (const auto &cell : dof_handler->cell_iterators())
-    {
-      pointerstruct =
-        cell_map.find(std::make_pair(cell->level(), cell->index()));
-
-      if (pointerstruct != cell_map_end)
-        // this cell was refined or not
-        // touched at all, so we can get
-        // the new values by just setting
-        // or interpolating to the children,
-        // which is both done by one
-        // function
-        {
-          const unsigned int this_fe_index =
-            pointerstruct->second.active_fe_index;
-          const unsigned int dofs_per_cell =
-            cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
-          local_values.reinit(dofs_per_cell, true);
-
-          // make sure that the size of the stored indices is the same as
-          // dofs_per_cell. since we store the desired fe_index, we know
-          // what this size should be
-          Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
-                 ExcInternalError());
-          for (unsigned int i = 0; i < dofs_per_cell; ++i)
-            local_values(i) = internal::ElementAccess<VectorType>::get(
-              in, (*pointerstruct->second.indices_ptr)[i]);
-          cell->set_dof_values_by_interpolation(local_values,
-                                                out,
-                                                this_fe_index,
-                                                true);
-        }
-    }
+  prepare_for_coarsening_and_refinement(all_in);
 }
 
 
 
-namespace internal
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::deserialize(VectorType &in)
 {
-  /**
-   * Generate a table that contains
-   * interpolation matrices between
-   * each combination of finite
-   * elements used in a DoFHandler of
-   * some kind. Since not all
-   * elements can be interpolated
-   * onto each other, the table may
-   * contain empty matrices for those
-   * combinations of elements for
-   * which no such interpolation is
-   * implemented.
-   */
-  template <int dim, int spacedim>
-  void
-  extract_interpolation_matrices(const DoFHandler<dim, spacedim>      &dof,
-                                 dealii::Table<2, FullMatrix<double>> &matrices)
-  {
-    if (dof.has_hp_capabilities() == false)
-      return;
-
-    const dealii::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
-    matrices.reinit(fe.size(), fe.size());
-    for (unsigned int i = 0; i < fe.size(); ++i)
-      for (unsigned int j = 0; j < fe.size(); ++j)
-        if (i != j)
-          {
-            matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
-                                  fe[j].n_dofs_per_cell());
-
-            // see if we can get the interpolation matrices for this
-            // combination of elements. if not, reset the matrix sizes to zero
-            // to indicate that this particular combination isn't
-            // supported. this isn't an outright error right away since we may
-            // never need to actually interpolate between these two elements
-            // on actual cells; we simply have to trigger an error if someone
-            // actually tries
-            try
-              {
-                fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
-              }
-            catch (const typename FiniteElement<dim, spacedim>::
-                     ExcInterpolationNotImplemented &)
-              {
-                matrices(i, j).reinit(0, 0);
-              }
-          }
-  }
+  std::vector<VectorType *> all_in(1, &in);
+  deserialize(all_in);
+}
 
 
-  template <int dim, int spacedim>
-  void
-  restriction_additive(const FiniteElement<dim, spacedim> &,
-                       std::vector<std::vector<bool>> &)
-  {}
 
-  template <int dim, int spacedim>
-  void
-  restriction_additive(const dealii::hp::FECollection<dim, spacedim> &fe,
-                       std::vector<std::vector<bool>> &restriction_is_additive)
-  {
-    restriction_is_additive.resize(fe.size());
-    for (unsigned int f = 0; f < fe.size(); ++f)
-      {
-        restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
-        for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
-          restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
-      }
-  }
-} // namespace internal
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::deserialize(
+  std::vector<VectorType *> &all_in)
+{
+  register_data_attach();
 
+  // this makes interpolate() happy
+  input_vectors.resize(all_in.size());
+
+  interpolate(all_in);
+}
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::
-  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+  std::vector<VectorType *> &all_out)
 {
-  Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
-  Assert(prepared_for != coarsening_and_refinement,
-         ExcAlreadyPrepForCoarseAndRef());
+  const dealii::internal::parallel::shared::
+    TemporarilyRestoreSubdomainIds<dim, spacedim>
+      subdomain_modifier(dof_handler->get_triangulation());
 
-  clear();
-  n_dofs_old                 = dof_handler->n_dofs();
-  const unsigned int in_size = all_in.size();
+  Assert(input_vectors.size() == all_out.size(),
+         ExcDimensionMismatch(input_vectors.size(), all_out.size()));
+  for (unsigned int i = 0; i < all_out.size(); ++i)
+    Assert(all_out[i]->size() == dof_handler->n_dofs(),
+           ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
 
-#ifdef DEBUG
-  Assert(in_size != 0,
-         ExcMessage("The array of input vectors you pass to this "
-                    "function has no elements. This is not useful."));
-  for (unsigned int i = 0; i < in_size; ++i)
+  // TODO: casting away constness is bad
+  auto tria = const_cast<dealii::Triangulation<dim, spacedim> *>(
+    &dof_handler->get_triangulation());
+  Assert(tria != nullptr, ExcInternalError());
+  Assert(
+    handle != numbers::invalid_unsigned_int,
+    ExcMessage(
+      "You can only call interpolate() once per SolutionTransfer object."));
+
+  if (average_values)
+    for (auto *const vec : all_out)
+      *vec = 0.0;
+
+  VectorType valence;
+
+  // initialize valence vector only if we need to average
+  if (average_values)
+    valence.reinit(*all_out[0]);
+
+  tria->notify_ready_to_unpack(
+    handle,
+    [this, &all_out, &valence](
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+      const CellStatus                                            status,
+      const boost::iterator_range<std::vector<char>::const_iterator>
+        &data_range) {
+      this->unpack_callback(cell_, status, data_range, all_out, valence);
+    });
+
+  if (average_values)
     {
-      Assert(all_in[i].size() == n_dofs_old,
-             ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+      // finalize valence: compress and invert
+      using Number = typename VectorType::value_type;
+      valence.compress(VectorOperation::add);
+      for (const auto i : valence.locally_owned_elements())
+        valence[i] = (static_cast<Number>(valence[i]) == Number() ?
+                        Number() :
+                        (Number(1.0) / static_cast<Number>(valence[i])));
+      valence.compress(VectorOperation::insert);
+
+      for (auto *const vec : all_out)
+        {
+          // compress and weight with valence
+          vec->compress(VectorOperation::add);
+          vec->scale(valence);
+        }
     }
-#endif
-
-  // We need to access dof indices on the entire domain. For
-  // parallel::shared::Triangulations, ownership of cells might change. If they
-  // allow artificial cells, we need to restore the "true" cell owners
-  // temporarily.
-  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
-  // the current set of subdomain ids, set subdomain ids to the "true" owner of
-  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
-  // and later restore these flags when it is destroyed.
-  const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
-                                                                   spacedim>
-    subdomain_modifier(dof_handler->get_triangulation());
-
-  // first count the number
-  // of cells that will be coarsened
-  // and that'll stay or be refined
-  unsigned int n_cells_to_coarsen        = 0;
-  unsigned int n_cells_to_stay_or_refine = 0;
-  for (const auto &act_cell : dof_handler->active_cell_iterators())
+  else
     {
-      if (act_cell->coarsen_flag_set())
-        ++n_cells_to_coarsen;
-      else
-        ++n_cells_to_stay_or_refine;
+      for (auto *const vec : all_out)
+        vec->compress(VectorOperation::insert);
     }
-  Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
-           dof_handler->get_triangulation().n_active_cells(),
-         ExcInternalError());
-
-  unsigned int n_coarsen_fathers = 0;
-  for (const auto &cell : dof_handler->cell_iterators())
-    if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
-      ++n_coarsen_fathers;
-  Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
-  (void)n_cells_to_coarsen;
-
-  // allocate the needed memory. initialize
-  // the following arrays in an efficient
-  // way, without copying much
-  std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
-    .swap(indices_on_cell);
-
-  std::vector<std::vector<Vector<typename VectorType::value_type>>>(
-    n_coarsen_fathers,
-    std::vector<Vector<typename VectorType::value_type>>(in_size))
-    .swap(dof_values_on_cell);
-
-  Table<2, FullMatrix<double>>   interpolation_hp;
-  std::vector<std::vector<bool>> restriction_is_additive;
-
-  internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
-  internal::restriction_additive(dof_handler->get_fe_collection(),
-                                 restriction_is_additive);
-
-  // we need counters for
-  // the 'to_stay_or_refine' cells 'n_sr' and
-  // the 'coarsen_fathers' cells 'n_cf',
-  unsigned int n_sr = 0, n_cf = 0;
-  for (const auto &cell : dof_handler->cell_iterators())
+
+  input_vectors.clear();
+  handle = numbers::invalid_unsigned_int;
+}
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+  std::vector<VectorType> &all_out)
+{
+  std::vector<VectorType *> temp(all_out.size());
+
+  for (std::size_t i = 0; i < temp.size(); ++i)
+    temp[i] = &(all_out[i]);
+
+  this->interpolate(temp);
+}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(VectorType &out)
+{
+  std::vector<VectorType *> all_out(1, &out);
+  interpolate(all_out);
+}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+std::vector<char>
+SolutionTransfer<dim, VectorType, spacedim>::pack_callback(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+  const CellStatus                                            status)
+{
+  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
+
+  // create buffer for each individual object
+  std::vector<::dealii::Vector<typename VectorType::value_type>> dof_values(
+    input_vectors.size());
+
+  unsigned int fe_index = 0;
+  if (dof_handler->has_hp_capabilities())
     {
-      // CASE 1: active cell that remains as it is
-      if (cell->is_active() && !cell->coarsen_flag_set())
+      switch (status)
         {
-          const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
-          indices_on_cell[n_sr].resize(dofs_per_cell);
-          // cell will not be coarsened,
-          // so we get away by storing the
-          // dof indices and later
-          // interpolating to the children
-          cell->get_dof_indices(indices_on_cell[n_sr]);
-          cell_map[std::make_pair(cell->level(), cell->index())] =
-            Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
-          ++n_sr;
-        }
+          case CellStatus::cell_will_persist:
+          case CellStatus::cell_will_be_refined:
+            {
+              fe_index = cell->future_fe_index();
+              break;
+            }
 
-      // CASE 2: cell is inactive but will become active
-      else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
-        {
-          // we will need to interpolate from the children of this cell
-          // to the current one. in the hp-context, this also means
-          // we need to figure out which finite element space to interpolate
-          // to since that is not implied by the global FE as in the non-hp-
-          // case. we choose the 'least dominant fe' on all children from
-          // the associated FECollection.
-          std::set<unsigned int> fe_indices_children;
-          for (const auto &child : cell->child_iterators())
+          case CellStatus::children_will_be_coarsened:
             {
-              Assert(child->is_active() && child->coarsen_flag_set(),
-                     typename dealii::Triangulation<
-                       dim>::ExcInconsistentCoarseningFlags());
+              // In case of coarsening, we need to find a suitable FE index
+              // for the parent cell. We choose the 'least dominant fe'
+              // on all children from the associated FECollection.
+#ifdef DEBUG
+              for (const auto &child : cell->child_iterators())
+                Assert(child->is_active() && child->coarsen_flag_set(),
+                       typename dealii::Triangulation<
+                         dim>::ExcInconsistentCoarseningFlags());
+#endif
 
-              fe_indices_children.insert(child->active_fe_index());
+              fe_index = dealii::internal::hp::DoFHandlerImplementation::
+                dominated_future_fe_on_children<dim, spacedim>(cell);
+              break;
             }
-          Assert(!fe_indices_children.empty(), ExcInternalError());
-
-          const unsigned int target_fe_index =
-            dof_handler->get_fe_collection().find_dominated_fe_extended(
-              fe_indices_children, /*codim=*/0);
-
-          Assert(target_fe_index != numbers::invalid_unsigned_int,
-                 internal::hp::DoFHandlerImplementation::
-                   ExcNoDominatedFiniteElementOnChildren());
-
-          const unsigned int dofs_per_cell =
-            dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
-
-          std::vector<Vector<typename VectorType::value_type>>(
-            in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
-            .swap(dof_values_on_cell[n_cf]);
-
-
-          // store the data of each of the input vectors. get this data
-          // as interpolated onto a finite element space that encompasses
-          // that of all the children. note that
-          // cell->get_interpolated_dof_values already does all of the
-          // interpolations between spaces
-          for (unsigned int j = 0; j < in_size; ++j)
-            cell->get_interpolated_dof_values(all_in[j],
-                                              dof_values_on_cell[n_cf][j],
-                                              target_fe_index);
-          cell_map[std::make_pair(cell->level(), cell->index())] =
-            Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
-          ++n_cf;
+
+          default:
+            Assert(false, ExcInternalError());
+            break;
         }
     }
-  Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
-  Assert(n_cf == n_coarsen_fathers, ExcInternalError());
 
-  prepared_for = coarsening_and_refinement;
+  const unsigned int dofs_per_cell =
+    dof_handler->get_fe(fe_index).n_dofs_per_cell();
+
+  if (dofs_per_cell == 0)
+    return std::vector<char>(); // nothing to do for FE_Nothing
+
+  auto it_input  = input_vectors.cbegin();
+  auto it_output = dof_values.begin();
+  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
+    {
+      it_output->reinit(dofs_per_cell);
+      cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
+    }
+
+  return pack_dof_values<typename VectorType::value_type>(dof_values,
+                                                          dofs_per_cell);
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::
-  prepare_for_coarsening_and_refinement(const VectorType &in)
+SolutionTransfer<dim, VectorType, spacedim>::unpack_callback(
+  const typename Triangulation<dim, spacedim>::cell_iterator     &cell_,
+  const CellStatus                                                status,
+  const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+  std::vector<VectorType *>                                      &all_out,
+  VectorType                                                     &valence)
 {
-  std::vector<VectorType> all_in(1, in);
-  prepare_for_coarsening_and_refinement(all_in);
+  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
+
+  unsigned int fe_index = 0;
+  if (dof_handler->has_hp_capabilities())
+    {
+      switch (status)
+        {
+          case CellStatus::cell_will_persist:
+          case CellStatus::children_will_be_coarsened:
+            {
+              fe_index = cell->active_fe_index();
+              break;
+            }
+
+          case CellStatus::cell_will_be_refined:
+            {
+              // After refinement, this particular cell is no longer active,
+              // and its children have inherited its FE index. However, to
+              // unpack the data on the old cell, we need to recover its FE
+              // index from one of the children. Just to be sure, we also
+              // check if all children have the same FE index.
+              fe_index = cell->child(0)->active_fe_index();
+              for (unsigned int child_index = 1;
+                   child_index < cell->n_children();
+                   ++child_index)
+                Assert(cell->child(child_index)->active_fe_index() == fe_index,
+                       ExcInternalError());
+              break;
+            }
+
+          default:
+            Assert(false, ExcInternalError());
+            break;
+        }
+    }
+
+  const unsigned int dofs_per_cell =
+    dof_handler->get_fe(fe_index).n_dofs_per_cell();
+
+  if (dofs_per_cell == 0)
+    return; // nothing to do for FE_Nothing
+
+  const std::vector<::dealii::Vector<typename VectorType::value_type>>
+    dof_values =
+      unpack_dof_values<typename VectorType::value_type>(data_range,
+                                                         dofs_per_cell);
+
+  // check if sizes match
+  AssertDimension(dof_values.size(), all_out.size());
+
+  // check if we have enough dofs provided by the FE object
+  // to interpolate the transferred data correctly
+  for (auto it_dof_values = dof_values.begin();
+       it_dof_values != dof_values.end();
+       ++it_dof_values)
+    Assert(
+      dofs_per_cell == it_dof_values->size(),
+      ExcMessage(
+        "The transferred data was packed with a different number of dofs than the "
+        "currently registered FE object assigned to the DoFHandler has."));
+
+  // distribute data for each registered vector on mesh
+  auto it_input  = dof_values.cbegin();
+  auto it_output = all_out.begin();
+  for (; it_input != dof_values.cend(); ++it_input, ++it_output)
+    if (average_values)
+      cell->distribute_local_to_global_by_interpolation(*it_input,
+                                                        *(*it_output),
+                                                        fe_index);
+    else
+      cell->set_dof_values_by_interpolation(*it_input,
+                                            *(*it_output),
+                                            fe_index,
+                                            true);
+
+  if (average_values)
+    {
+      // compute valence vector if averaging should be performed
+      Vector<typename VectorType::value_type> ones(dofs_per_cell);
+      ones = 1.0;
+      cell->distribute_local_to_global_by_interpolation(ones,
+                                                        valence,
+                                                        fe_index);
+    }
 }
 
 
 
 template <int dim, typename VectorType, int spacedim>
 void
-SolutionTransfer<dim, VectorType, spacedim>::interpolate(
-  const std::vector<VectorType> &all_in,
-  std::vector<VectorType>       &all_out) const
+SolutionTransfer<dim, VectorType, spacedim>::clear()
 {
-  const unsigned int size = all_in.size();
-#ifdef DEBUG
-  Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
-  Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
-  for (unsigned int i = 0; i < size; ++i)
-    Assert(all_in[i].size() == n_dofs_old,
-           ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
-  for (unsigned int i = 0; i < all_out.size(); ++i)
-    Assert(all_out[i].size() == dof_handler->n_dofs(),
-           ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
-  for (unsigned int i = 0; i < size; ++i)
-    for (unsigned int j = 0; j < size; ++j)
-      Assert(&all_in[i] != &all_out[j],
-             ExcMessage("Vectors cannot be used as input and output"
-                        " at the same time!"));
-#endif
+  // nothing to do
+}
 
-  // We need to access dof indices on the entire domain. For
-  // parallel::shared::Triangulations, ownership of cells might change. If they
-  // allow artificial cells, we need to restore the "true" cell owners
-  // temporarily.
-  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
-  // the current set of subdomain ids, set subdomain ids to the "true" owner of
-  // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
-  // and later restore these flags when it is destroyed.
-  const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
-                                                                   spacedim>
-    subdomain_modifier(dof_handler->get_triangulation());
-
-  Vector<typename VectorType::value_type> local_values;
-  std::vector<types::global_dof_index>    dofs;
-
-  typename std::map<std::pair<unsigned int, unsigned int>,
-                    Pointerstruct>::const_iterator pointerstruct,
-    cell_map_end = cell_map.end();
-
-  Table<2, FullMatrix<double>> interpolation_hp;
-  internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
-  Vector<typename VectorType::value_type> tmp, tmp2;
-
-  for (const auto &cell : dof_handler->cell_iterators())
-    {
-      pointerstruct =
-        cell_map.find(std::make_pair(cell->level(), cell->index()));
 
-      if (pointerstruct != cell_map_end)
-        {
-          const std::vector<types::global_dof_index> *const indexptr =
-            pointerstruct->second.indices_ptr;
 
-          const std::vector<Vector<typename VectorType::value_type>>
-            *const valuesptr = pointerstruct->second.dof_values_ptr;
+namespace Legacy
+{
 
-          // cell stayed as it was or was refined
-          if (indexptr != nullptr)
-            {
-              Assert(valuesptr == nullptr, ExcInternalError());
+  template <int dim, typename VectorType, int spacedim>
+  SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
+    const DoFHandler<dim, spacedim> &dof)
+    : dof_handler(&dof, typeid(*this).name())
+    , n_dofs_old(0)
+    , prepared_for(none)
+  {
+    Assert((dynamic_cast<
+              const parallel::distributed::Triangulation<dim, spacedim> *>(
+              &dof_handler->get_triangulation()) == nullptr),
+           ExcMessage(
+             "You are calling the dealii::SolutionTransfer class "
+             "with a DoFHandler that is built on a "
+             "parallel::distributed::Triangulation. This will not "
+             "work for parallel computations. You probably want to "
+             "use the parallel::distributed::SolutionTransfer class."));
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  SolutionTransfer<dim, VectorType, spacedim>::~SolutionTransfer()
+  {
+    clear();
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::clear()
+  {
+    indices_on_cell.clear();
+    dof_values_on_cell.clear();
+    cell_map.clear();
+
+    prepared_for = none;
+  }
 
-              const unsigned int old_fe_index =
-                pointerstruct->second.active_fe_index;
 
-              // get the values of each of the input data vectors on this cell
-              // and prolong it to its children
-              unsigned int in_size = indexptr->size();
-              for (unsigned int j = 0; j < size; ++j)
+
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::prepare_for_pure_refinement()
+  {
+    Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
+    Assert(prepared_for != coarsening_and_refinement,
+           ExcAlreadyPrepForCoarseAndRef());
+
+    clear();
+
+    // We need to access dof indices on the entire domain. For
+    // parallel::shared::Triangulations, ownership of cells might change. If
+    // they allow artificial cells, we need to restore the "true" cell owners
+    // temporarily.
+    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+    // the current set of subdomain ids, set subdomain ids to the "true" owner
+    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+    // object, and later restore these flags when it is destroyed.
+    const dealii::internal::parallel::shared::
+      TemporarilyRestoreSubdomainIds<dim, spacedim>
+        subdomain_modifier(dof_handler->get_triangulation());
+
+    const unsigned int n_active_cells =
+      dof_handler->get_triangulation().n_active_cells();
+    n_dofs_old = dof_handler->n_dofs();
+
+    // efficient reallocation of indices_on_cell
+    std::vector<std::vector<types::global_dof_index>>(n_active_cells)
+      .swap(indices_on_cell);
+
+    for (const auto &cell : dof_handler->active_cell_iterators())
+      {
+        const unsigned int i = cell->active_cell_index();
+        indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
+        // on each cell store the indices of the
+        // dofs. after refining we get the values
+        // on the children by taking these
+        // indices, getting the respective values
+        // out of the data vectors and prolonging
+        // them to the children
+        cell->get_dof_indices(indices_on_cell[i]);
+        cell_map[std::make_pair(cell->level(), cell->index())] =
+          Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
+      }
+    prepared_for = pure_refinement;
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
+    const VectorType &in,
+    VectorType       &out) const
+  {
+    Assert(prepared_for == pure_refinement, ExcNotPrepared());
+    Assert(in.size() == n_dofs_old,
+           ExcDimensionMismatch(in.size(), n_dofs_old));
+    Assert(out.size() == dof_handler->n_dofs(),
+           ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+    Assert(&in != &out,
+           ExcMessage("Vectors cannot be used as input and output"
+                      " at the same time!"));
+
+    // We need to access dof indices on the entire domain. For
+    // parallel::shared::Triangulations, ownership of cells might change. If
+    // they allow artificial cells, we need to restore the "true" cell owners
+    // temporarily.
+    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+    // the current set of subdomain ids, set subdomain ids to the "true" owner
+    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+    // object, and later restore these flags when it is destroyed.
+    const dealii::internal::parallel::shared::
+      TemporarilyRestoreSubdomainIds<dim, spacedim>
+        subdomain_modifier(dof_handler->get_triangulation());
+
+    Vector<typename VectorType::value_type> local_values(0);
+
+    typename std::map<std::pair<unsigned int, unsigned int>,
+                      Pointerstruct>::const_iterator pointerstruct,
+      cell_map_end = cell_map.end();
+
+    for (const auto &cell : dof_handler->cell_iterators())
+      {
+        pointerstruct =
+          cell_map.find(std::make_pair(cell->level(), cell->index()));
+
+        if (pointerstruct != cell_map_end)
+          // this cell was refined or not
+          // touched at all, so we can get
+          // the new values by just setting
+          // or interpolating to the children,
+          // which is both done by one
+          // function
+          {
+            const unsigned int this_fe_index =
+              pointerstruct->second.active_fe_index;
+            const unsigned int dofs_per_cell =
+              cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
+            local_values.reinit(dofs_per_cell, true);
+
+            // make sure that the size of the stored indices is the same as
+            // dofs_per_cell. since we store the desired fe_index, we know
+            // what this size should be
+            Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
+                   ExcInternalError());
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              local_values(i) =
+                dealii::internal::ElementAccess<VectorType>::get(
+                  in, (*pointerstruct->second.indices_ptr)[i]);
+            cell->set_dof_values_by_interpolation(local_values,
+                                                  out,
+                                                  this_fe_index,
+                                                  true);
+          }
+      }
+  }
+
+
+
+  namespace internal
+  {
+    /**
+     * Generate a table that contains
+     * interpolation matrices between
+     * each combination of finite
+     * elements used in a DoFHandler of
+     * some kind. Since not all
+     * elements can be interpolated
+     * onto each other, the table may
+     * contain empty matrices for those
+     * combinations of elements for
+     * which no such interpolation is
+     * implemented.
+     */
+    template <int dim, int spacedim>
+    void
+    extract_interpolation_matrices(
+      const DoFHandler<dim, spacedim>      &dof,
+      dealii::Table<2, FullMatrix<double>> &matrices)
+    {
+      if (dof.has_hp_capabilities() == false)
+        return;
+
+      const dealii::hp::FECollection<dim, spacedim> &fe =
+        dof.get_fe_collection();
+      matrices.reinit(fe.size(), fe.size());
+      for (unsigned int i = 0; i < fe.size(); ++i)
+        for (unsigned int j = 0; j < fe.size(); ++j)
+          if (i != j)
+            {
+              matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
+                                    fe[j].n_dofs_per_cell());
+
+              // see if we can get the interpolation matrices for this
+              // combination of elements. if not, reset the matrix sizes to zero
+              // to indicate that this particular combination isn't
+              // supported. this isn't an outright error right away since we may
+              // never need to actually interpolate between these two elements
+              // on actual cells; we simply have to trigger an error if someone
+              // actually tries
+              try
                 {
-                  tmp.reinit(in_size, true);
-                  for (unsigned int i = 0; i < in_size; ++i)
-                    tmp(i) =
-                      internal::ElementAccess<VectorType>::get(all_in[j],
-                                                               (*indexptr)[i]);
-
-                  cell->set_dof_values_by_interpolation(tmp,
-                                                        all_out[j],
-                                                        old_fe_index,
-                                                        true);
+                  fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
                 }
-            }
-          else if (valuesptr)
-            // the children of this cell were deleted
-            {
-              Assert(!cell->has_children(), ExcInternalError());
-              Assert(indexptr == nullptr, ExcInternalError());
-
-              const unsigned int dofs_per_cell =
-                cell->get_fe().n_dofs_per_cell();
-              dofs.resize(dofs_per_cell);
-              // get the local
-              // indices
-              cell->get_dof_indices(dofs);
-
-              // distribute the stored data to the new vectors
-              for (unsigned int j = 0; j < size; ++j)
+              catch (const typename FiniteElement<dim, spacedim>::
+                       ExcInterpolationNotImplemented &)
                 {
-                  // make sure that the size of the stored indices is the same
-                  // as dofs_per_cell. this is kind of a test if we use the same
-                  // FE in the hp-case. to really do that test we would have to
-                  // store the fe_index of all cells
-                  const Vector<typename VectorType::value_type> *data = nullptr;
-                  const unsigned int active_fe_index = cell->active_fe_index();
-                  if (active_fe_index != pointerstruct->second.active_fe_index)
-                    {
-                      const unsigned int old_index =
-                        pointerstruct->second.active_fe_index;
-                      const FullMatrix<double> &interpolation_matrix =
-                        interpolation_hp(active_fe_index, old_index);
-                      // The interpolation matrix might be empty when using
-                      // FE_Nothing.
-                      if (interpolation_matrix.empty())
-                        tmp.reinit(dofs_per_cell, false);
-                      else
-                        {
-                          tmp.reinit(dofs_per_cell, true);
-                          AssertDimension((*valuesptr)[j].size(),
-                                          interpolation_matrix.n());
-                          AssertDimension(tmp.size(), interpolation_matrix.m());
-                          interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
-                        }
-                      data = &tmp;
-                    }
-                  else
-                    data = &(*valuesptr)[j];
-
-
-                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
-                    internal::ElementAccess<VectorType>::set((*data)(i),
-                                                             dofs[i],
-                                                             all_out[j]);
+                  matrices(i, j).reinit(0, 0);
                 }
             }
-          // undefined status
-          else
-            DEAL_II_ASSERT_UNREACHABLE();
+    }
+
+
+    template <int dim, int spacedim>
+    void
+    restriction_additive(const FiniteElement<dim, spacedim> &,
+                         std::vector<std::vector<bool>> &)
+    {}
+
+    template <int dim, int spacedim>
+    void
+    restriction_additive(
+      const dealii::hp::FECollection<dim, spacedim> &fe,
+      std::vector<std::vector<bool>>                &restriction_is_additive)
+    {
+      restriction_is_additive.resize(fe.size());
+      for (unsigned int f = 0; f < fe.size(); ++f)
+        {
+          restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
+          for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
+            restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
         }
     }
+  } // namespace internal
 
-  // We have written into the output vectors. If this was a PETSc vector, for
-  // example, then we need to compress these to make future operations safe:
-  for (auto &vec : all_out)
-    vec.compress(VectorOperation::insert);
-}
 
 
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::
+    prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
+  {
+    Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
+    Assert(prepared_for != coarsening_and_refinement,
+           ExcAlreadyPrepForCoarseAndRef());
 
-template <int dim, typename VectorType, int spacedim>
-void
-SolutionTransfer<dim, VectorType, spacedim>::interpolate(const VectorType &in,
-                                                         VectorType &out) const
-{
-  Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
-  Assert(out.size() == dof_handler->n_dofs(),
-         ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+    clear();
+    n_dofs_old                 = dof_handler->n_dofs();
+    const unsigned int in_size = all_in.size();
 
-  std::vector<VectorType> all_in  = {in};
-  std::vector<VectorType> all_out = {out};
+#ifdef DEBUG
+    Assert(in_size != 0,
+           ExcMessage("The array of input vectors you pass to this "
+                      "function has no elements. This is not useful."));
+    for (unsigned int i = 0; i < in_size; ++i)
+      {
+        Assert(all_in[i].size() == n_dofs_old,
+               ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+      }
+#endif
 
-  interpolate(all_in, all_out);
+    // We need to access dof indices on the entire domain. For
+    // parallel::shared::Triangulations, ownership of cells might change. If
+    // they allow artificial cells, we need to restore the "true" cell owners
+    // temporarily.
+    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+    // the current set of subdomain ids, set subdomain ids to the "true" owner
+    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+    // object, and later restore these flags when it is destroyed.
+    const dealii::internal::parallel::shared::
+      TemporarilyRestoreSubdomainIds<dim, spacedim>
+        subdomain_modifier(dof_handler->get_triangulation());
+
+    // first count the number
+    // of cells that will be coarsened
+    // and that'll stay or be refined
+    unsigned int n_cells_to_coarsen        = 0;
+    unsigned int n_cells_to_stay_or_refine = 0;
+    for (const auto &act_cell : dof_handler->active_cell_iterators())
+      {
+        if (act_cell->coarsen_flag_set())
+          ++n_cells_to_coarsen;
+        else
+          ++n_cells_to_stay_or_refine;
+      }
+    Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
+             dof_handler->get_triangulation().n_active_cells(),
+           ExcInternalError());
+
+    unsigned int n_coarsen_fathers = 0;
+    for (const auto &cell : dof_handler->cell_iterators())
+      if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
+        ++n_coarsen_fathers;
+    Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
+    (void)n_cells_to_coarsen;
+
+    // allocate the needed memory. initialize
+    // the following arrays in an efficient
+    // way, without copying much
+    std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
+      .swap(indices_on_cell);
+
+    std::vector<std::vector<Vector<typename VectorType::value_type>>>(
+      n_coarsen_fathers,
+      std::vector<Vector<typename VectorType::value_type>>(in_size))
+      .swap(dof_values_on_cell);
+
+    Table<2, FullMatrix<double>>   interpolation_hp;
+    std::vector<std::vector<bool>> restriction_is_additive;
+
+    internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
+    internal::restriction_additive(dof_handler->get_fe_collection(),
+                                   restriction_is_additive);
+
+    // we need counters for
+    // the 'to_stay_or_refine' cells 'n_sr' and
+    // the 'coarsen_fathers' cells 'n_cf',
+    unsigned int n_sr = 0, n_cf = 0;
+    for (const auto &cell : dof_handler->cell_iterators())
+      {
+        // CASE 1: active cell that remains as it is
+        if (cell->is_active() && !cell->coarsen_flag_set())
+          {
+            const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
+            indices_on_cell[n_sr].resize(dofs_per_cell);
+            // cell will not be coarsened,
+            // so we get away by storing the
+            // dof indices and later
+            // interpolating to the children
+            cell->get_dof_indices(indices_on_cell[n_sr]);
+            cell_map[std::make_pair(cell->level(), cell->index())] =
+              Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
+            ++n_sr;
+          }
 
-  out = all_out[0];
-}
+        // CASE 2: cell is inactive but will become active
+        else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
+          {
+            // we will need to interpolate from the children of this cell
+            // to the current one. in the hp-context, this also means
+            // we need to figure out which finite element space to interpolate
+            // to since that is not implied by the global FE as in the non-hp-
+            // case. we choose the 'least dominant fe' on all children from
+            // the associated FECollection.
+            std::set<unsigned int> fe_indices_children;
+            for (const auto &child : cell->child_iterators())
+              {
+                Assert(child->is_active() && child->coarsen_flag_set(),
+                       typename dealii::Triangulation<
+                         dim>::ExcInconsistentCoarseningFlags());
 
+                fe_indices_children.insert(child->active_fe_index());
+              }
+            Assert(!fe_indices_children.empty(), ExcInternalError());
+
+            const unsigned int target_fe_index =
+              dof_handler->get_fe_collection().find_dominated_fe_extended(
+                fe_indices_children, /*codim=*/0);
+
+            Assert(target_fe_index != numbers::invalid_unsigned_int,
+                   dealii::internal::hp::DoFHandlerImplementation::
+                     ExcNoDominatedFiniteElementOnChildren());
+
+            const unsigned int dofs_per_cell =
+              dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
+
+            std::vector<Vector<typename VectorType::value_type>>(
+              in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
+              .swap(dof_values_on_cell[n_cf]);
+
+
+            // store the data of each of the input vectors. get this data
+            // as interpolated onto a finite element space that encompasses
+            // that of all the children. note that
+            // cell->get_interpolated_dof_values already does all of the
+            // interpolations between spaces
+            for (unsigned int j = 0; j < in_size; ++j)
+              cell->get_interpolated_dof_values(all_in[j],
+                                                dof_values_on_cell[n_cf][j],
+                                                target_fe_index);
+            cell_map[std::make_pair(cell->level(), cell->index())] =
+              Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
+            ++n_cf;
+          }
+      }
+    Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
+    Assert(n_cf == n_coarsen_fathers, ExcInternalError());
 
+    prepared_for = coarsening_and_refinement;
+  }
 
-template <int dim, typename VectorType, int spacedim>
-std::size_t
-SolutionTransfer<dim, VectorType, spacedim>::memory_consumption() const
-{
-  // at the moment we do not include the memory
-  // consumption of the cell_map as we have no
-  // real idea about memory consumption of a
-  // std::map
-  return (MemoryConsumption::memory_consumption(dof_handler) +
-          MemoryConsumption::memory_consumption(n_dofs_old) +
-          sizeof(prepared_for) +
-          MemoryConsumption::memory_consumption(indices_on_cell) +
-          MemoryConsumption::memory_consumption(dof_values_on_cell));
-}
 
 
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::
+    prepare_for_coarsening_and_refinement(const VectorType &in)
+  {
+    std::vector<VectorType> all_in(1, in);
+    prepare_for_coarsening_and_refinement(all_in);
+  }
 
-template <int dim, typename VectorType, int spacedim>
-std::size_t
-SolutionTransfer<dim, VectorType, spacedim>::Pointerstruct::memory_consumption()
-  const
-{
-  return sizeof(*this);
-}
+
+
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+    const std::vector<VectorType> &all_in,
+    std::vector<VectorType>       &all_out) const
+  {
+    const unsigned int size = all_in.size();
+#ifdef DEBUG
+    Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
+    Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
+    for (unsigned int i = 0; i < size; ++i)
+      Assert(all_in[i].size() == n_dofs_old,
+             ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+    for (unsigned int i = 0; i < all_out.size(); ++i)
+      Assert(all_out[i].size() == dof_handler->n_dofs(),
+             ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
+    for (unsigned int i = 0; i < size; ++i)
+      for (unsigned int j = 0; j < size; ++j)
+        Assert(&all_in[i] != &all_out[j],
+               ExcMessage("Vectors cannot be used as input and output"
+                          " at the same time!"));
+#endif
+
+    // We need to access dof indices on the entire domain. For
+    // parallel::shared::Triangulations, ownership of cells might change. If
+    // they allow artificial cells, we need to restore the "true" cell owners
+    // temporarily.
+    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+    // the current set of subdomain ids, set subdomain ids to the "true" owner
+    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+    // object, and later restore these flags when it is destroyed.
+    const dealii::internal::parallel::shared::
+      TemporarilyRestoreSubdomainIds<dim, spacedim>
+        subdomain_modifier(dof_handler->get_triangulation());
+
+    Vector<typename VectorType::value_type> local_values;
+    std::vector<types::global_dof_index>    dofs;
+
+    typename std::map<std::pair<unsigned int, unsigned int>,
+                      Pointerstruct>::const_iterator pointerstruct,
+      cell_map_end = cell_map.end();
+
+    Table<2, FullMatrix<double>> interpolation_hp;
+    internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
+    Vector<typename VectorType::value_type> tmp, tmp2;
+
+    for (const auto &cell : dof_handler->cell_iterators())
+      {
+        pointerstruct =
+          cell_map.find(std::make_pair(cell->level(), cell->index()));
+
+        if (pointerstruct != cell_map_end)
+          {
+            const std::vector<types::global_dof_index> *const indexptr =
+              pointerstruct->second.indices_ptr;
+
+            const std::vector<Vector<typename VectorType::value_type>>
+              *const valuesptr = pointerstruct->second.dof_values_ptr;
+
+            // cell stayed as it was or was refined
+            if (indexptr != nullptr)
+              {
+                Assert(valuesptr == nullptr, ExcInternalError());
+
+                const unsigned int old_fe_index =
+                  pointerstruct->second.active_fe_index;
+
+                // get the values of each of the input data vectors on this cell
+                // and prolong it to its children
+                unsigned int in_size = indexptr->size();
+                for (unsigned int j = 0; j < size; ++j)
+                  {
+                    tmp.reinit(in_size, true);
+                    for (unsigned int i = 0; i < in_size; ++i)
+                      tmp(i) = dealii::internal::ElementAccess<VectorType>::get(
+                        all_in[j], (*indexptr)[i]);
+
+                    cell->set_dof_values_by_interpolation(tmp,
+                                                          all_out[j],
+                                                          old_fe_index,
+                                                          true);
+                  }
+              }
+            else if (valuesptr)
+              // the children of this cell were deleted
+              {
+                Assert(!cell->has_children(), ExcInternalError());
+                Assert(indexptr == nullptr, ExcInternalError());
+
+                const unsigned int dofs_per_cell =
+                  cell->get_fe().n_dofs_per_cell();
+                dofs.resize(dofs_per_cell);
+                // get the local
+                // indices
+                cell->get_dof_indices(dofs);
+
+                // distribute the stored data to the new vectors
+                for (unsigned int j = 0; j < size; ++j)
+                  {
+                    // make sure that the size of the stored indices is the same
+                    // as dofs_per_cell. this is kind of a test if we use the
+                    // same FE in the hp-case. to really do that test we would
+                    // have to store the fe_index of all cells
+                    const Vector<typename VectorType::value_type> *data =
+                      nullptr;
+                    const unsigned int active_fe_index =
+                      cell->active_fe_index();
+                    if (active_fe_index !=
+                        pointerstruct->second.active_fe_index)
+                      {
+                        const unsigned int old_index =
+                          pointerstruct->second.active_fe_index;
+                        const FullMatrix<double> &interpolation_matrix =
+                          interpolation_hp(active_fe_index, old_index);
+                        // The interpolation matrix might be empty when using
+                        // FE_Nothing.
+                        if (interpolation_matrix.empty())
+                          tmp.reinit(dofs_per_cell, false);
+                        else
+                          {
+                            tmp.reinit(dofs_per_cell, true);
+                            AssertDimension((*valuesptr)[j].size(),
+                                            interpolation_matrix.n());
+                            AssertDimension(tmp.size(),
+                                            interpolation_matrix.m());
+                            interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
+                          }
+                        data = &tmp;
+                      }
+                    else
+                      data = &(*valuesptr)[j];
+
+
+                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                      dealii::internal::ElementAccess<VectorType>::set(
+                        (*data)(i), dofs[i], all_out[j]);
+                  }
+              }
+            // undefined status
+            else
+              Assert(false, ExcInternalError());
+          }
+      }
+
+    // We have written into the output vectors. If this was a PETSc vector, for
+    // example, then we need to compress these to make future operations safe:
+    for (auto &vec : all_out)
+      vec.compress(VectorOperation::insert);
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  void
+  SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+    const VectorType &in,
+    VectorType       &out) const
+  {
+    Assert(in.size() == n_dofs_old,
+           ExcDimensionMismatch(in.size(), n_dofs_old));
+    Assert(out.size() == dof_handler->n_dofs(),
+           ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+
+    std::vector<VectorType> all_in  = {in};
+    std::vector<VectorType> all_out = {out};
+
+    interpolate(all_in, all_out);
+
+    out = all_out[0];
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  std::size_t
+  SolutionTransfer<dim, VectorType, spacedim>::memory_consumption() const
+  {
+    // at the moment we do not include the memory
+    // consumption of the cell_map as we have no
+    // real idea about memory consumption of a
+    // std::map
+    return (MemoryConsumption::memory_consumption(dof_handler) +
+            MemoryConsumption::memory_consumption(n_dofs_old) +
+            sizeof(prepared_for) +
+            MemoryConsumption::memory_consumption(indices_on_cell) +
+            MemoryConsumption::memory_consumption(dof_values_on_cell));
+  }
+
+
+
+  template <int dim, typename VectorType, int spacedim>
+  std::size_t
+  SolutionTransfer<dim, VectorType, spacedim>::Pointerstruct::
+    memory_consumption() const
+  {
+    return sizeof(*this);
+  }
+
+} // namespace Legacy
 
 
 /*-------------- Explicit Instantiations -------------------------------*/
index a53d7569c808f914bbdea20b3477795126575414..ba92230013d9d91daec7487b3ce6de374e61af9c 100644 (file)
@@ -17,6 +17,9 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
      deal_II_space_dimension : SPACE_DIMENSIONS)
   {
 #if deal_II_dimension <= deal_II_space_dimension
+    template class Legacy::
+      SolutionTransfer<deal_II_dimension, VEC, deal_II_space_dimension>;
+
     template class SolutionTransfer<deal_II_dimension,
                                     VEC,
                                     deal_II_space_dimension>;
index b3e9c50c7970136345e8ed9f33fc6d2a52352bc6..1ad6174712dc76185fa37ee316ec96b9861740eb 100644 (file)
@@ -93,12 +93,12 @@ transfer(std::ostream &out)
     cell->set_refine_flag(RefinementCase<dim>::cut_x);
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(solution);
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  soltrans.refine_interpolate(solution, new_solution);
+  soltrans.interpolate(new_solution);
   solution.reinit(dof_handler.n_dofs());
   solution = new_solution;
 
@@ -122,7 +122,7 @@ transfer(std::ostream &out)
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   solution.reinit(dof_handler.n_dofs());
-  soltrans2.interpolate(old_solution, solution);
+  soltrans2.interpolate(solution);
 
   data_out.clear_data_vectors();
   data_out.add_data_vector(solution, "solution");
index 468b44c728cd66720baed81f4dc754a927433ea2..7b47a79f9fe907afb86a74649f21221f79825266 100644 (file)
@@ -115,22 +115,17 @@ transfer(std::ostream &out)
     cell->set_refine_flag();
 
   tria.prepare_coarsening_and_refinement();
-  q_soltrans.prepare_for_pure_refinement();
-  dgq_soltrans.prepare_for_pure_refinement();
+  q_soltrans.prepare_for_coarsening_and_refinement(q_solution);
+  dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_solution);
   tria.execute_coarsening_and_refinement();
   q_dof_handler.distribute_dofs(fe_q);
   dgq_dof_handler.distribute_dofs(fe_dgq);
 
-  Vector<double> tmp_q(q_dof_handler.n_dofs());
-  q_soltrans.refine_interpolate(q_solution, tmp_q);
   q_solution.reinit(q_dof_handler.n_dofs());
-  q_solution = tmp_q;
+  q_soltrans.interpolate(q_solution);
 
-  Vector<double> tmp_dgq(dgq_dof_handler.n_dofs());
-  dgq_soltrans.refine_interpolate(dgq_solution, tmp_dgq);
   dgq_solution.reinit(dgq_dof_handler.n_dofs());
-  dgq_solution = tmp_dgq;
-
+  dgq_soltrans.interpolate(dgq_solution);
 
   q_data_out.clear_data_vectors();
   q_data_out.add_data_vector(q_solution, "solution");
@@ -167,8 +162,8 @@ transfer(std::ostream &out)
   dgq_dof_handler.distribute_dofs(fe_dgq);
   q_solution.reinit(q_dof_handler.n_dofs());
   dgq_solution.reinit(dgq_dof_handler.n_dofs());
-  q_soltrans.interpolate(q_old_solution, q_solution);
-  dgq_soltrans.interpolate(dgq_old_solution, dgq_solution);
+  q_soltrans.interpolate(q_solution);
+  dgq_soltrans.interpolate(dgq_solution);
 
   q_data_out.clear_data_vectors();
   q_data_out.add_data_vector(q_solution, "solution");
index 45d1f2ff865c5cf6452e8ecad0181f4c3a95e6d9..6919afcf27cec2665a20222337f44c3119c7dbe7 100644 (file)
@@ -87,7 +87,7 @@ main()
   // get the interpolated solution
   // back
   Vector<double> tmp(dh.n_dofs());
-  soltrans.interpolate(solution, tmp);
+  soltrans.interpolate(tmp);
 
   deallog << "New values:" << std::endl;
   for (unsigned int i = 0; i < tmp.size(); ++i)
index 3dbcfa5c7d2bbfb50845db196038259efe5bc03e..237ea34e96b573323e9e7ecc8148a47a7adb03c9 100644 (file)
@@ -33,9 +33,6 @@
 
 #include "../tests.h"
 
-#include "coarse_grid_common.h"
-
-
 
 template <int dim>
 void
index 547fbea36e222f3b93d4423dba72fc1a91678e8e..f9027add3484d3efce0a9cfe1563e5317d6d19fa 100644 (file)
@@ -18,6 +18,7 @@
 
 #include <deal.II/distributed/fully_distributed_tria.h>
 #include <deal.II/distributed/solution_transfer.h>
+#include <deal.II/distributed/tria.h>
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
index fa0bb31eb36f013edc163560e984921e5bfb1f34..e3b6f2f2d8326479907ed32e904143d54428e422 100644 (file)
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
 
-#include <deal.II/fe/fe_q.h>
-
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_description.h>
-
-#include <deal.II/lac/vector.h>
 
 #include "./tests.h"
 
diff --git a/tests/grid/solution_transfer_01.cc b/tests/grid/solution_transfer_01.cc
new file mode 100644 (file)
index 0000000..c6b1c77
--- /dev/null
@@ -0,0 +1,123 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SolutionTransfer::load() and save() for serial iangulations.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "./tests.h"
+
+
+template <int dim>
+class InterpolationFunction : public Function<dim>
+{
+public:
+  InterpolationFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    return p.norm();
+  }
+};
+
+template <int dim, typename TriangulationType>
+void
+test(TriangulationType &triangulation)
+{
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(FE_Q<dim>(2));
+
+  using VectorType = Vector<double>;
+
+  VectorType vector(dof_handler.n_dofs());
+
+  VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), vector);
+
+  VectorType vector_loaded(dof_handler.n_dofs());
+
+  const std::string filename =
+    "solution_transfer_" + std::to_string(dim) + "d_out";
+
+  {
+    SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+    solution_transfer.prepare_for_serialization(vector);
+
+    triangulation.save(filename);
+  }
+
+  triangulation.clear();
+
+  {
+    triangulation.load(filename);
+    dof_handler.distribute_dofs(FE_Q<dim>(2));
+
+    SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+    solution_transfer.deserialize(vector_loaded);
+  }
+
+  // Verify that error is 0.
+  VectorType error(vector);
+  error.add(-1, vector_loaded);
+
+  deallog << (error.linfty_norm() < 1e-16 ? "PASSED" : "FAILED") << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    constexpr int dim = 2;
+
+    Triangulation<dim> triangulation;
+    GridGenerator::hyper_cube(triangulation);
+    triangulation.refine_global(3);
+
+    test<dim>(triangulation);
+  }
+  deallog.pop();
+
+  deallog.push("3d");
+  {
+    constexpr int dim = 3;
+
+    Triangulation<dim> triangulation;
+    GridGenerator::hyper_cube(triangulation);
+    triangulation.refine_global(3);
+
+    test<dim>(triangulation);
+  }
+  deallog.pop();
+}
diff --git a/tests/grid/solution_transfer_01.output b/tests/grid/solution_transfer_01.output
new file mode 100644 (file)
index 0000000..7e803b7
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:2d::PASSED
+DEAL:3d::PASSED
diff --git a/tests/grid/solution_transfer_02.cc b/tests/grid/solution_transfer_02.cc
new file mode 100644 (file)
index 0000000..eebb50c
--- /dev/null
@@ -0,0 +1,133 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SolutionTransfer::interpolate() for serial triangulations.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "./tests.h"
+
+
+template <int dim>
+class InterpolationFunction : public Function<dim>
+{
+public:
+  InterpolationFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    return p[0];
+  }
+};
+
+template <int dim>
+void
+test(const unsigned int type)
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+  triangulation.refine_global(3);
+
+  const FE_Q<dim> fe(2);
+
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  using VectorType = Vector<double>;
+
+  VectorType vector(dof_handler.n_dofs());
+  VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), vector);
+
+  SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+
+  triangulation.prepare_coarsening_and_refinement();
+  solution_transfer.prepare_for_coarsening_and_refinement(vector);
+
+  if (type == 0)
+    {
+      triangulation.refine_global(1);
+    }
+  else if (type == 1)
+    {
+      for (const auto &cell : triangulation.active_cell_iterators())
+        if (cell->center()[0] < 0.5)
+          cell->set_refine_flag();
+      triangulation.execute_coarsening_and_refinement();
+    }
+  else if (type == 2)
+    {
+      for (const auto &cell : triangulation.active_cell_iterators())
+        if (cell->center()[0] < 0.5)
+          cell->set_coarsen_flag();
+      triangulation.execute_coarsening_and_refinement();
+    }
+
+  dof_handler.distribute_dofs(fe);
+
+  vector.reinit(dof_handler.n_dofs());
+
+  solution_transfer.interpolate(vector);
+
+  VectorType error(dof_handler.n_dofs());
+  VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), error);
+
+  error -= vector;
+
+  deallog << (error.linfty_norm() < 1e-16 ? "PASSED" : "FAILED") << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.push("2d");
+  {
+    constexpr int dim = 2;
+
+    test<dim>(0);
+    test<dim>(1);
+    test<dim>(2);
+  }
+  deallog.pop();
+
+  deallog.push("3d");
+  {
+    constexpr int dim = 3;
+
+    test<dim>(0);
+    test<dim>(1);
+    test<dim>(2);
+  }
+  deallog.pop();
+}
diff --git a/tests/grid/solution_transfer_02.output b/tests/grid/solution_transfer_02.output
new file mode 100644 (file)
index 0000000..da1d891
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:2d::PASSED
+DEAL:2d::PASSED
+DEAL:2d::PASSED
+DEAL:3d::PASSED
+DEAL:3d::PASSED
+DEAL:3d::PASSED
index e748d4694de0925fcb7e328eabce5a6058ac0867..05c0c74688bb78885f283692be28e64ed8c7b3df 100644 (file)
@@ -86,7 +86,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   hp::QCollection<2> q;
   q.push_back(QMidpoint<2>());
index b27980d44b931f4dd87df40457e5fb9520e5454e..96b361bfcc31bc218562947f3f7b37efbce74222 100644 (file)
@@ -86,7 +86,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   hp::QCollection<2> q;
   q.push_back(QMidpoint<2>());
index a6d196df9f80705facb1da7cc4e3c67f11d334cd..800d217b72f90764016ca061436c44bb6a52d1bd 100644 (file)
@@ -86,7 +86,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   hp::QCollection<2> q;
   q.push_back(QMidpoint<2>());
index 7dbc05b7b4963b194884ce768b590916c52c420f..8bb35257bcfdb7c0a12ab64135d7184e18fc6041 100644 (file)
@@ -75,12 +75,12 @@ transfer(std::ostream &out)
     cell->set_refine_flag();
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(solution);
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
 
   Vector<double> tmp_q(dof_handler.n_dofs());
-  soltrans.refine_interpolate(solution, tmp_q);
+  soltrans.interpolate(tmp_q);
   solution.reinit(dof_handler.n_dofs());
   solution = tmp_q;
 
@@ -99,7 +99,7 @@ transfer(std::ostream &out)
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   solution.reinit(dof_handler.n_dofs());
-  soltrans.interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 
   deallog << "OK" << std::endl;
 }
index 65b7010e3e488f50a7c12bec4635793aa5cbedfe..d5a0bead530919ee5eeda1249f28d1a0dd515eef 100644 (file)
@@ -463,7 +463,7 @@ Laplace<dim>::refine_grid(const unsigned int cycle)
   setup_system();
 
   // 3.7. Solution Transfer finish
-  soltrans.interpolate(solution_coarse, solution);
+  soltrans.interpolate(solution);
 }
 
 
index ad13548bbc8eded0d74e12208d00a619ff64846c..e86cc1ef97a5cb87c12bbd639ba37629eb165fd1 100644 (file)
@@ -151,12 +151,6 @@ transfer(std::ostream &out)
         cell->set_coarsen_flag();
     }
 
-  Vector<double> q_old_solution = q_solution, dgq_old_solution = dgq_solution;
-  tria.prepare_coarsening_and_refinement();
-  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
-  dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
-  tria.execute_coarsening_and_refinement();
-
   counter = 0;
   {
     typename DoFHandler<dim>::active_cell_iterator
@@ -165,22 +159,28 @@ transfer(std::ostream &out)
     for (; cell != endc; ++cell, ++celldg, ++counter)
       {
         if (counter > 20 && counter < 90)
-          cell->set_active_fe_index(0);
+          cell->set_future_fe_index(0);
         else
-          cell->set_active_fe_index(Testing::rand() % max_degree);
+          cell->set_future_fe_index(Testing::rand() % max_degree);
         if (counter > 20 && counter < 90)
-          celldg->set_active_fe_index(0);
+          celldg->set_future_fe_index(0);
         else
-          celldg->set_active_fe_index(Testing::rand() % max_degree);
+          celldg->set_future_fe_index(Testing::rand() % max_degree);
       }
   }
 
+  Vector<double> q_old_solution = q_solution, dgq_old_solution = dgq_solution;
+  tria.prepare_coarsening_and_refinement();
+  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+  dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
+  tria.execute_coarsening_and_refinement();
+
   q_dof_handler.distribute_dofs(fe_q);
   dgq_dof_handler.distribute_dofs(fe_dgq);
   q_solution.reinit(q_dof_handler.n_dofs());
   dgq_solution.reinit(dgq_dof_handler.n_dofs());
-  q_soltrans.interpolate(q_old_solution, q_solution);
-  dgq_soltrans.interpolate(dgq_old_solution, dgq_solution);
+  q_soltrans.interpolate(q_solution);
+  dgq_soltrans.interpolate(dgq_solution);
 
   // check correctness by comparing the values
   // on points of QGauss of order 2.
index 622928532d76cc4d9810c99aabb6b4781fa68b1b..479898f4cd89010c5971d4586acacf15391911b3 100644 (file)
@@ -91,7 +91,7 @@ transfer(std::ostream &out)
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   solution.reinit(dof_handler.n_dofs());
-  soltrans.interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index 75edf39d6b346ebcf0bb97fbc90f6fdf4dc5b728..b59dd1ce3ee4adfab4b79439c279764c1f804d39 100644 (file)
@@ -92,19 +92,21 @@ main()
 
   // Interpolate solution
   SolutionTransfer<2, Vector<double>> solultion_trans(dof_handler);
-  solultion_trans.prepare_for_coarsening_and_refinement(solution);
 
-  triangulation.execute_coarsening_and_refinement();
   // Assign FE_Q_ to all cells
   cell = dof_handler.begin_active();
   for (; cell != endc; ++cell)
     {
-      cell->set_active_fe_index(0);
+      cell->set_future_fe_index(0);
     }
+
+  solultion_trans.prepare_for_coarsening_and_refinement(solution);
+
+  triangulation.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
 
   // Save output
index bb0c53d6de6238afbbb1e8d85354d503d7f5b600..e214319cb0b4d55058af95553beea75b81d1cc1b 100644 (file)
@@ -85,7 +85,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   // we are good if we made it to here
   deallog << "OK" << std::endl;
index a08d922c3c4309624a57fe5a56c146d3ee7a6f7f..2aca704e7c2ac0a65aa1b8a74676212abd233fc3 100644 (file)
@@ -85,7 +85,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   // Define compression level for output data
   DataOutBase::VtkFlags vtk_flags;
index 3350ff158e906ef1ba6d7fdf67af74673236b178..ee8647e5443cdf02d0e895b150bda6ca5cb9c382 100644 (file)
@@ -65,18 +65,19 @@ test()
   // set refine flag for the only cell we have, then do the refinement
   SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
   dof_handler.begin_active()->set_refine_flag();
-  solution_trans.prepare_for_coarsening_and_refinement(solution);
-  triangulation.execute_coarsening_and_refinement();
 
   // now set the active_fe_index flags on the new set of fine level cells
   for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
-    dof_handler.begin(0)->child(c)->set_active_fe_index(1);
+    dof_handler.begin(0)->child(c)->set_future_fe_index(1);
+
+  solution_trans.prepare_for_coarsening_and_refinement(solution);
+  triangulation.execute_coarsening_and_refinement();
 
   // distribute dofs and transfer solution there
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solution_trans.interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   // we should now have only 1s in the new_solution vector
   for (unsigned int i = 0; i < new_solution.size(); ++i)
index 08dbaa72586c1a879186ce4d2258391e520d25b9..ee8647e5443cdf02d0e895b150bda6ca5cb9c382 100644 (file)
@@ -65,18 +65,19 @@ test()
   // set refine flag for the only cell we have, then do the refinement
   SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
   dof_handler.begin_active()->set_refine_flag();
-  solution_trans.prepare_for_pure_refinement();
-  triangulation.execute_coarsening_and_refinement();
 
   // now set the active_fe_index flags on the new set of fine level cells
   for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
-    dof_handler.begin(0)->child(c)->set_active_fe_index(1);
+    dof_handler.begin(0)->child(c)->set_future_fe_index(1);
+
+  solution_trans.prepare_for_coarsening_and_refinement(solution);
+  triangulation.execute_coarsening_and_refinement();
 
   // distribute dofs and transfer solution there
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solution_trans.refine_interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   // we should now have only 1s in the new_solution vector
   for (unsigned int i = 0; i < new_solution.size(); ++i)
index 9321d9e9c81471a14e9be3a102a870a2de927f87..2f368b626201a65cd2406a6e663e3fbb958f7866 100644 (file)
@@ -66,17 +66,17 @@ test()
   for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
     dof_handler.begin(0)->child(c)->set_coarsen_flag();
 
+  // now set the active_fe_index flags on the only remaining cell
+  dof_handler.begin_active()->set_future_fe_index(0);
+
   solution_trans.prepare_for_coarsening_and_refinement(solution);
   triangulation.execute_coarsening_and_refinement();
 
-  // now set the active_fe_index flags on the only remaining cell
-  dof_handler.begin_active()->set_active_fe_index(0);
-
   // distribute dofs and transfer solution there
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solution_trans.interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   // we should now have only 1s in the new_solution vector
   for (unsigned int i = 0; i < new_solution.size(); ++i)
index 4695ed34546de0320b4687f6ca27aded08f74b61..9a21fa9209b1402b61fe94ce95962613954a9ab0 100644 (file)
@@ -63,7 +63,7 @@ test()
   // set refine flag for the only cell we have, then do the refinement
   SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
   dof_handler.begin_active()->set_refine_flag();
-  solution_trans.prepare_for_pure_refinement();
+  solution_trans.prepare_for_coarsening_and_refinement(solution);
   triangulation.execute_coarsening_and_refinement();
 
   // now set the active_fe_index flags on the new set of fine level cells
@@ -74,7 +74,7 @@ test()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solution_trans.refine_interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   // we should now have only 1s in the new_solution vector
   for (unsigned int i = 0; i < new_solution.size(); ++i)
index 7f0abd91dff5153a625f543a60dffdbe38c09be5..c9b84f0195fed979392f4d0187be79a7d3bcd53d 100644 (file)
@@ -105,7 +105,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solultion_trans.interpolate(solution, new_solution);
+  solultion_trans.interpolate(new_solution);
 
   deallog << "OK" << std::endl;
 }
index 678b0010181c3d010d3e2f8a989fc161d2a70d15..9dc5efc907ed50013e7a2281b6ef10260908ee6f 100644 (file)
@@ -74,11 +74,11 @@ transfer(std::ostream &out)
 
   Vector<double> old_solution = solution;
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(solution);
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   solution.reinit(dof_handler.n_dofs());
-  soltrans.refine_interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index a2ee19e2aaadf7d21051c9c007c81200f00146b3..9473c6e57be78dd5637633760ab3b152fe4c9ba1 100644 (file)
@@ -177,11 +177,6 @@ transfer(std::ostream &out)
         cell->set_coarsen_flag();
     }
 
-  Vector<double> q_old_solution = q_solution;
-  tria.prepare_coarsening_and_refinement();
-  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
-  tria.execute_coarsening_and_refinement();
-
   counter = 0;
   {
     typename DoFHandler<dim>::active_cell_iterator cell = q_dof_handler
@@ -191,15 +186,20 @@ transfer(std::ostream &out)
     for (; cell != endc; ++cell, ++counter)
       {
         if (counter > 20 && counter < 90)
-          cell->set_active_fe_index(0);
+          cell->set_future_fe_index(0);
         else
-          cell->set_active_fe_index(Testing::rand() % max_degree);
+          cell->set_future_fe_index(Testing::rand() % max_degree);
       }
   }
 
+  Vector<double> q_old_solution = q_solution;
+  tria.prepare_coarsening_and_refinement();
+  q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+  tria.execute_coarsening_and_refinement();
+
   q_dof_handler.distribute_dofs(fe_q);
   q_solution.reinit(q_dof_handler.n_dofs());
-  q_soltrans.interpolate(q_old_solution, q_solution);
+  q_soltrans.interpolate(q_solution);
 
   // check correctness by comparing the values
   // on points of QGauss of order 2.
index b0a79b637ae0af71055ad1e6952315ff35b55727..66f78351071319d6c7924045ff81f295c3d19751 100644 (file)
@@ -132,7 +132,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  solution_trans.interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   FE_Type.reinit(triangulation.n_active_cells());
   cnt_cells = 0;
@@ -182,7 +182,7 @@ main()
   dof_handler.distribute_dofs(fe_collection);
 
   Vector<double> new_solution2(dof_handler.n_dofs());
-  solution_trans2.interpolate(solution, new_solution2);
+  solution_trans2.interpolate(new_solution2);
 
   FE_Type.reinit(triangulation.n_active_cells());
   cnt_cells = 0;
index f197a50a064acae95366271f2aed02008fc46a8f..99afd752a786429aa6563257b61690e308aa22d5 100644 (file)
@@ -95,13 +95,13 @@ transfer(const MPI_Comm mpi_communicator)
   old_solution = solution;
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(old_solution);
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   locally_owned_dofs =
     DoFTools::locally_owned_dofs_per_subdomain(dof_handler)[this_mpi_process];
   solution.reinit(locally_owned_dofs, mpi_communicator);
-  soltrans.refine_interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index a220af6bff43d571a0ce261665ee8526b4a2c042..05e89b180ca640feb23e5110e8ff3edf9ffee4fe 100644 (file)
@@ -150,7 +150,7 @@ transfer(const MPI_Comm mpi_communicator)
   old_solution = solution;
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(old_solution);
   tria.execute_coarsening_and_refinement();
 
   dof_handler.distribute_dofs(fe);
@@ -165,7 +165,7 @@ transfer(const MPI_Comm mpi_communicator)
                        this_mpi_process);
 
   solution.reinit(locally_owned_partitioning, mpi_communicator);
-  soltrans.refine_interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index 73f4a7e4688917b9f829306b57431795f7e6fefc..d3c3968f54ecc5aebee54f509a9f42c1bcad3d70 100644 (file)
@@ -67,6 +67,10 @@ main()
   for (const auto &cell : dof_handler.active_cell_iterators())
     cell->set_coarsen_flag();
 
+  // Assign FE_Q(1) to all cells
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    cell->set_future_fe_index(0);
+
   triangulation.prepare_coarsening_and_refinement();
 
   // Interpolate solution
@@ -76,16 +80,11 @@ main()
 
   triangulation.execute_coarsening_and_refinement();
 
-  // Assign FE_Q(1) to all cells
-  for (const auto &cell : dof_handler.active_cell_iterators())
-    cell->set_active_fe_index(0);
-
   dof_handler.distribute_dofs(fe_collection);
   deallog << "Final number of dofs: " << dof_handler.n_dofs() << std::endl;
 
   Vector<double> new_solution(dof_handler.n_dofs());
-  new_solution = 1.;
-  solution_trans.interpolate(solution, new_solution);
+  solution_trans.interpolate(new_solution);
 
   deallog << "Vector after solution transfer:" << std::endl;
   new_solution.print(deallog.get_file_stream());
index 3da6d79d61db00fb93206a04de409aa1d383a562..ec31643421822fa3eb241d8fdaf52f78e6089241 100644 (file)
@@ -90,7 +90,7 @@ transfer(const MPI_Comm mpi_communicator)
   old_solution = solution;
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(old_solution);
   tria.execute_coarsening_and_refinement();
   dof_handler.distribute_dofs(fe);
   locally_owned_dofs =
@@ -98,7 +98,7 @@ transfer(const MPI_Comm mpi_communicator)
   locally_relevant_dofs = DoFTools::locally_relevant_dofs_per_subdomain(
     dof_handler)[this_mpi_process];
   solution.reinit(locally_owned_dofs, mpi_communicator);
-  soltrans.refine_interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index 82db17533f8acef1cf2bd314e55cbff9bc4dda4b..e2a3ba15c62afeed7a0db7b40f9b7c0a49d877dd 100644 (file)
@@ -144,7 +144,7 @@ transfer(const MPI_Comm mpi_communicator)
   old_solution = solution;
 
   tria.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(old_solution);
   tria.execute_coarsening_and_refinement();
 
   dof_handler.distribute_dofs(fe);
@@ -159,7 +159,7 @@ transfer(const MPI_Comm mpi_communicator)
                        this_mpi_process);
 
   solution.reinit(locally_owned_partitioning, mpi_communicator);
-  soltrans.refine_interpolate(old_solution, solution);
+  soltrans.interpolate(solution);
 }
 
 
index 9905d62448917d42e0f4f0275c83eb15c5504e0f..7c60b4c72ddee84cce973e66de6c68fcf79b8757 100644 (file)
@@ -62,7 +62,8 @@ test(const bool allow_artificial_cells)
   dh.distribute_dofs(fe);
 
   Vector<double> sol_new(dh.n_dofs());
-  soltrans.interpolate(sol_old, sol_new);
+  soltrans.interpolate(sol_new);
+
   for (unsigned int i = 0; i < sol_new.size(); ++i)
     AssertThrow(sol_new[i] == 1., ExcInternalError());
 
index 2f1dffded2bc4f5e782a9e86d1768086b3f13a1b..ca9d9578081d5bf9431745f3f6888133bbe9d290 100644 (file)
@@ -56,13 +56,13 @@ test(const bool allow_artificial_cells)
   sol_old = 1.;
 
   SolutionTransfer<dim> soltrans(dh);
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(sol_old);
 
   tria.execute_coarsening_and_refinement();
   dh.distribute_dofs(fe);
 
   Vector<double> sol_new(dh.n_dofs());
-  soltrans.refine_interpolate(sol_old, sol_new);
+  soltrans.interpolate(sol_new);
   for (unsigned int i = 0; i < sol_new.size(); ++i)
     AssertThrow(sol_new[i] == 1., ExcInternalError());
 
index 1fb9607b942a6589ba2153a41ab12c341f3cea0f..bbcad94f9bc9def143bfdd365fcf399c72e8744f 100644 (file)
@@ -61,7 +61,7 @@ refine_and_transfer(const Function<dim> &function, Triangulation<dim> &tria)
   dh.distribute_dofs(fe);
 
   Vector<double> sol_new(dh.n_dofs());
-  soltrans.interpolate(sol_old, sol_new);
+  soltrans.interpolate(sol_new);
 
   return sol_new;
 }
index b8a0c090fd8e42b5a5019b4c730da57015659aee..52e71b4f7a181d4786f8df035ddeccecd209187f 100644 (file)
@@ -55,13 +55,13 @@ refine_and_transfer(const Function<dim> &function, Triangulation<dim> &tria)
   VectorTools::interpolate(MappingQ<dim>(1), dh, function, sol_old);
 
   SolutionTransfer<dim> soltrans(dh);
-  soltrans.prepare_for_pure_refinement();
+  soltrans.prepare_for_coarsening_and_refinement(sol_old);
 
   tria.execute_coarsening_and_refinement();
   dh.distribute_dofs(fe);
 
   Vector<double> sol_new(dh.n_dofs());
-  soltrans.refine_interpolate(sol_old, sol_new);
+  soltrans.interpolate(sol_new);
 
   return sol_new;
 }
index 6c03f996806e66e5a9784f1fa3d64d39eb54d33c..495f020308b0c9492033883b66e0305f9f47f6f9 100644 (file)
@@ -69,7 +69,7 @@ test()
 
   triangulation.execute_coarsening_and_refinement();
 
-  solution_trans.interpolate(previous_solution, solution);
+  solution_trans.interpolate(solution);
 
   deallog << "OK" << std::endl;
 }
index 6c03f996806e66e5a9784f1fa3d64d39eb54d33c..495f020308b0c9492033883b66e0305f9f47f6f9 100644 (file)
@@ -69,7 +69,7 @@ test()
 
   triangulation.execute_coarsening_and_refinement();
 
-  solution_trans.interpolate(previous_solution, solution);
+  solution_trans.interpolate(solution);
 
   deallog << "OK" << std::endl;
 }
index c3dfdd47bc61b89ad731a7515329d8526a05c08e..d732e30b9686438f306da87634b367f9047edcae 100644 (file)
@@ -1013,12 +1013,12 @@ namespace Step31
     std::vector<TrilinosWrappers::MPI::Vector> tmp(2);
     tmp[0].reinit(temperature_solution);
     tmp[1].reinit(temperature_solution);
-    temperature_trans.interpolate(x_temperature, tmp);
+    temperature_trans.interpolate(tmp);
     temperature_solution     = tmp[0];
     old_temperature_solution = tmp[1];
     temperature_constraints.distribute(temperature_solution);
     temperature_constraints.distribute(old_temperature_solution);
-    stokes_trans.interpolate(x_stokes, stokes_solution);
+    stokes_trans.interpolate(stokes_solution);
     stokes_constraints.distribute(stokes_solution);
     rebuild_stokes_matrix         = true;
     rebuild_temperature_matrices  = true;
index 859d9c6cce9ab23d83c8f86e1411731cb0241a0f..9b0b0dcec12465a54b0b06c21c256bfd3f1b6381 100644 (file)
@@ -339,7 +339,7 @@ namespace Step77
     dof_handler.distribute_dofs(fe);
 
     Vector<double> tmp(dof_handler.n_dofs());
-    solution_transfer.interpolate(current_solution, tmp);
+    solution_transfer.interpolate(tmp);
     current_solution = std::move(tmp);
 
     hanging_node_constraints.clear();
index bb2d05ee38568aff6fc853b75164fbb35d2e1763..737c298aec400a76338917ba38cd9ace8766cda6 100644 (file)
@@ -345,7 +345,7 @@ namespace Step77
     dof_handler.distribute_dofs(fe);
 
     Vector<double> tmp(dof_handler.n_dofs());
-    solution_transfer.interpolate(current_solution, tmp);
+    solution_transfer.interpolate(tmp);
     current_solution = std::move(tmp);
 
     hanging_node_constraints.clear();

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.