From b3b6590016b9b52f52b4a07036a4b3c3da1e6c40 Mon Sep 17 00:00:00 2001 From: Pasquale Africa Date: Mon, 10 Jul 2023 18:55:31 +0000 Subject: [PATCH] Merge SolutionTransfer and p::d::SolutionTransfer --- .../20240103AfricaBlaisMunch | 10 + examples/step-15/step-15.cc | 2 +- examples/step-26/step-26.cc | 2 +- examples/step-28/step-28.cc | 2 +- examples/step-31/step-31.cc | 4 +- examples/step-33/step-33.cc | 2 +- examples/step-43/step-43.cc | 4 +- examples/step-57/step-57.cc | 2 +- examples/step-72/step-72.cc | 2 +- examples/step-77/step-77.cc | 3 +- .../deal.II/distributed/solution_transfer.h | 349 +--- include/deal.II/distributed/tria_base.h | 6 +- include/deal.II/grid/tria.h | 16 + include/deal.II/numerics/solution_transfer.h | 1116 +++++++++---- source/distributed/CMakeLists.txt | 2 - source/distributed/fully_distributed_tria.cc | 18 +- source/distributed/solution_transfer.cc | 482 ------ source/distributed/solution_transfer.inst.in | 66 - source/grid/tria.cc | 211 ++- source/numerics/solution_transfer.cc | 1467 +++++++++++------ source/numerics/solution_transfer.inst.in | 3 + tests/aniso/solution_transfer.cc | 6 +- tests/bits/solution_transfer.cc | 17 +- tests/codim_one/solution_transfer_01.cc | 2 +- .../distributed_grids/solution_transfer_01.cc | 3 - .../solution_transfer_01.cc | 1 + tests/grid/save_load_01.cc | 5 - tests/grid/solution_transfer_01.cc | 123 ++ tests/grid/solution_transfer_01.output | 3 + tests/grid/solution_transfer_02.cc | 133 ++ tests/grid/solution_transfer_02.output | 7 + tests/hp/do_function_derivatives_01.cc | 2 +- tests/hp/do_function_hessians_01.cc | 2 +- tests/hp/do_function_laplacians_01.cc | 2 +- tests/hp/fe_nothing_17.cc | 6 +- tests/hp/laplace.h | 2 +- tests/hp/solution_transfer.cc | 24 +- tests/hp/solution_transfer_02.cc | 2 +- tests/hp/solution_transfer_03.cc | 10 +- tests/hp/solution_transfer_04.cc | 2 +- tests/hp/solution_transfer_05.cc | 2 +- tests/hp/solution_transfer_06.cc | 9 +- tests/hp/solution_transfer_07.cc | 9 +- tests/hp/solution_transfer_08.cc | 8 +- tests/hp/solution_transfer_09.cc | 4 +- tests/hp/solution_transfer_10.cc | 2 +- tests/hp/solution_transfer_11.cc | 4 +- tests/hp/solution_transfer_12.cc | 16 +- tests/hp/solution_transfer_13.cc | 4 +- tests/hp/solution_transfer_14.cc | 4 +- tests/hp/solution_transfer_15.cc | 4 +- tests/hp/solution_transfer_16.cc | 11 +- tests/mpi/solution_transfer_02.cc | 4 +- tests/mpi/solution_transfer_03.cc | 4 +- tests/mpi/solution_transfer_07.cc | 3 +- tests/mpi/solution_transfer_08.cc | 4 +- tests/mpi/solution_transfer_09.cc | 2 +- tests/mpi/solution_transfer_10.cc | 4 +- tests/petsc/solution_transfer_01.cc | 2 +- tests/petsc/solution_transfer_02.cc | 2 +- tests/simplex/step-31.cc | 4 +- tests/sundials/step-77.cc | 2 +- tests/trilinos/step-77-with-nox.cc | 2 +- 63 files changed, 2368 insertions(+), 1863 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20240103AfricaBlaisMunch delete mode 100644 source/distributed/solution_transfer.cc delete mode 100644 source/distributed/solution_transfer.inst.in create mode 100644 tests/grid/solution_transfer_01.cc create mode 100644 tests/grid/solution_transfer_01.output create mode 100644 tests/grid/solution_transfer_02.cc create mode 100644 tests/grid/solution_transfer_02.output diff --git a/doc/news/changes/incompatibilities/20240103AfricaBlaisMunch b/doc/news/changes/incompatibilities/20240103AfricaBlaisMunch new file mode 100644 index 0000000000..a1c898a17a --- /dev/null +++ b/doc/news/changes/incompatibilities/20240103AfricaBlaisMunch @@ -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. +
+(Pasquale Claudio Africa, Bruno Blais, Peter Munch, 2024/01/03) diff --git a/examples/step-15/step-15.cc b/examples/step-15/step-15.cc index acb0a61feb..991b42c376 100644 --- a/examples/step-15/step-15.cc +++ b/examples/step-15/step-15.cc @@ -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 diff --git a/examples/step-26/step-26.cc b/examples/step-26/step-26.cc index a1493d5648..77e0d1bec2 100644 --- a/examples/step-26/step-26.cc +++ b/examples/step-26/step-26.cc @@ -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); } diff --git a/examples/step-28/step-28.cc b/examples/step-28/step-28.cc index ad664acc22..4b3452850a 100644 --- a/examples/step-28/step-28.cc +++ b/examples/step-28/step-28.cc @@ -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: diff --git a/examples/step-31/step-31.cc b/examples/step-31/step-31.cc index 08c3c88cd7..6ee6303a89 100644 --- a/examples/step-31/step-31.cc +++ b/examples/step-31/step-31.cc @@ -1977,7 +1977,7 @@ namespace Step31 std::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); diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 0610f0d604..5ba1fa882e 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -2300,7 +2300,7 @@ namespace Step33 std::vector> transfer_out = { Vector(dof_handler.n_dofs()), Vector(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]); diff --git a/examples/step-43/step-43.cc b/examples/step-43/step-43.cc index 36d22c13f8..2b82b6f765 100644 --- a/examples/step-43/step-43.cc +++ b/examples/step-43/step-43.cc @@ -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 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]; diff --git a/examples/step-57/step-57.cc b/examples/step-57/step-57.cc index 29fe6ef671..68f03e665d 100644 --- a/examples/step-57/step-57.cc +++ b/examples/step-57/step-57.cc @@ -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 diff --git a/examples/step-72/step-72.cc b/examples/step-72/step-72.cc index e7d42365ed..4eb70f2e8d 100644 --- a/examples/step-72/step-72.cc +++ b/examples/step-72/step-72.cc @@ -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); } diff --git a/examples/step-77/step-77.cc b/examples/step-77/step-77.cc index dc4d1436e9..50216ca92b 100644 --- a/examples/step-77/step-77.cc +++ b/examples/step-77/step-77.cc @@ -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); } diff --git a/include/deal.II/distributed/solution_transfer.h b/include/deal.II/distributed/solution_transfer.h index caf9320da8..857976b04d 100644 --- a/include/deal.II/distributed/solution_transfer.h +++ b/include/deal.II/distributed/solution_transfer.h @@ -17,11 +17,7 @@ #include -#include - -#include - -#include +#include 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. - * - *

Note on ghost elements

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. - * - *

Transferring a solution

Here VectorType is your favorite - * vector type, e.g. PETScWrappers::MPI::Vector, - * TrilinosWrappers::MPI::Vector, or corresponding block vectors. - * @code - * parallel::distributed::SolutionTransfer - * 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 - * 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 - * - *

Use for Serialization

- * - * 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 - * 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 - * sol_trans(dof_handler); - * sol_trans.deserialize(distributed_vector); - * @endcode - * - * - *

Note on usage with DoFHandler with hp-capabilities

- * - * 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> - * 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 fe_collection; - * //[prepare identical fe_collection...] - * - * DoFHandler 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> - * sol_trans(hp_dof_handler); - * sol_trans.deserialize(distributed_vector); - * @endcode - * - * - *

Interaction with hanging nodes

- * - * 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 - 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 &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 &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 &all_out); - - /** - * Same as the previous function. It interpolates only one function. It - * assumes the vectors having the right sizes (i.e. - * in.size()==n_dofs_old, out.size()==n_dofs_refined) - * - * Multiple calling of this function is NOT allowed. Interpolating - * several functions can be performed in one step by using - * prepare_for_coarsening_and_refinement(all_in) and - * interpolate(all_out). - */ - 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 &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 &all_in); - - private: - /** - * Pointer to the degree of freedom handler to work with. - */ - SmartPointer, - SolutionTransfer> - 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 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 - pack_callback( - const typename Triangulation::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::cell_iterator &cell, - const CellStatus status, - const boost::iterator_range::const_iterator> - &data_range, - std::vector &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; } // namespace distributed } // namespace parallel diff --git a/include/deal.II/distributed/tria_base.h b/include/deal.II/distributed/tria_base.h index c9f8fc6cba..2ebec1f077 100644 --- a/include/deal.II/distributed/tria_base.h +++ b/include/deal.II/distributed/tria_base.h @@ -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; }; diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index e0a956b1c1..bbd55ec651 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -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 diff --git a/include/deal.II/numerics/solution_transfer.h b/include/deal.II/numerics/solution_transfer.h index 6817b1d5fe..0fb3a5c385 100644 --- a/include/deal.II/numerics/solution_transfer.h +++ b/include/deal.II/numerics/solution_transfer.h @@ -32,100 +32,60 @@ 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. * *

Usage

* - * This class implements the algorithms in two different ways: - *
    - *
  • If the grid will only be refined (i.e. no cells are coarsened) then - * use @p SolutionTransfer as follows: - * @code - * SolutionTransfer > 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 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. + *

    Transferring a solution

    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 > solutions_old(n_vectors, Vector (n)); - * ... - * std::vector > solutions(n_vectors, Vector (n)); - * soltrans.refine_interpolate(solutions_old, solutions); - * @endcode - * This is used in several of the tutorial programs, for example step-31 - * and step-33. - * - *
  • If the grid has cells that will be coarsened, then use @p - * SolutionTransfer as follows: - * @code - * SolutionTransfer > soltrans(*dof_handler); - * + * SolutionTransfer 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 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): + *

    Usage on distributed grids

    + * 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 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 + * + *

    Note on ghost elements

    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 interpolate (const VectorType &in, - * VectorType &out) are NOT allowed. Interpolating several - * functions can be performed in one step by using void interpolate (const - * vector &all_in, vector &all_out) - * const, 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). - *
- * - * For deleting all stored data in @p SolutionTransfer and reinitializing it - * use the clear() function. - * - * The template argument @p VectorType denotes the type of data container you - * want to transfer. - * - * - *

Interpolating in the presence of hanging nodes and boundary values

- * - * 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. - * - * - *

Implementation

- * - *
    - *
  • 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 refine_interpolate(in,out) 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 refine_interpolate(in,out) function can be called multiple - * times for arbitrary many discrete functions (solution vectors) on the - * original grid. - * - *
  • 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 - * prepare_for_coarsening_and_refinement(all_in) where the - * vector all_in includes all discrete - * functions to be interpolated onto the new grid. - * - * As we need two different kinds of pointers (vector * - * for the Dof indices and vector * 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 void * pointer as - * vector * at one time and as - * vector * 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 interpolate(all_in, all_out) 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 - * *' pointer in @p Pointerstruct that is pointed to by the pointer of that - * cell. It is clear that interpolate(all_in, all_out) only can be - * called with the vector all_in that previously was - * the parameter of the prepare_for_coarsening_and_refinement(all_in) - * function. Hence interpolate(all_in, all_out) can (in contrast to - * refine_interpolate(in, out)) only be called once. - *
+ *

Use for serialization

+ * + * 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 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 sol_trans(dof_handler); + * sol_trans.deserialize(distributed_vector); + * @endcode + * + * + *

Note on usage with DoFHandler with hp-capabilities

+ * + * 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> + * 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 fe_collection; + * //[prepare identical fe_collection...] + * + * DoFHandler 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> + * sol_trans(hp_dof_handler); + * sol_trans.deserialize(distributed_vector); + * @endcode * * *

Interaction with hanging nodes

@@ -293,34 +258,34 @@ DEAL_II_NAMESPACE_OPEN * *

Implementation in the context of hp-finite elements

* - * 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 &dof); + SolutionTransfer(const DoFHandler &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 &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 &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. in.size()==n_dofs_old, - * out.size()==n_dofs_refined) - * - * 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 &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 - * prepare_for_refinement_and_coarsening(all_in). 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 &all_in, - std::vector &all_out) const; + interpolate(std::vector &all_out); /** * Same as the previous function. It interpolates only one function. It * assumes the vectors having the right sizes (i.e. * in.size()==n_dofs_old, out.size()==n_dofs_refined) * - * Multiple calling of this function is NOT allowed. Interpolating several - * functions can be performed in one step by using interpolate (all_in, - * all_out) + * Multiple calling of this function is NOT allowed. Interpolating + * several functions can be performed in one step by using + * interpolate (all_in, all_out) + */ + 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 &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 &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 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> indices_on_cell; + std::vector + pack_callback( + const typename Triangulation::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 vector indices - * (if the cell will be refined) or the vector dof_values - * (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 *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> *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::cell_iterator &cell, + const CellStatus status, + const boost::iterator_range::const_iterator> &data_range, + std::vector &all_out, + VectorType &valence); - std::vector *indices_ptr; - std::vector> *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, 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. + * + *

Usage

+ * + * This class implements the algorithms in two different ways: + *
    + *
  • If the grid will only be refined (i.e. no cells are coarsened) then + * use @p SolutionTransfer as follows: + * @code + * SolutionTransfer > 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 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 > solutions_old(n_vectors, Vector (n)); + * ... + * std::vector > solutions(n_vectors, Vector (n)); + * soltrans.refine_interpolate(solutions_old, solutions); + * @endcode + * This is used in several of the tutorial programs, for example step-31 + * and step-33. + * + *
  • If the grid has cells that will be coarsened, then use @p + * SolutionTransfer as follows: + * @code + * SolutionTransfer > 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 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 interpolate (const VectorType &in, + * VectorType &out) are NOT allowed. Interpolating several + * functions can be performed in one step by using void interpolate (const + * vector &all_in, vector &all_out) + * const, 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). + *
+ * + * For deleting all stored data in @p SolutionTransfer and reinitializing it + * use the clear() function. + * + * The template argument @p VectorType denotes the type of data container you + * want to transfer. + * + * + *

Interpolating in the presence of hanging nodes and boundary values

+ * + * 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. + * + * + *

Implementation

+ * + *
    + *
  • 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 refine_interpolate(in,out) 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 refine_interpolate(in,out) function can be called multiple + * times for arbitrary many discrete functions (solution vectors) on the + * original grid. + * + *
  • 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 + * prepare_for_coarsening_and_refinement(all_in) where the + * vector all_in includes all discrete + * functions to be interpolated onto the new grid. + * + * As we need two different kinds of pointers (vector * + * for the Dof indices and vector * 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 void * pointer as + * vector * at one time and as + * vector * 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 interpolate(all_in, all_out) 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 + * *' pointer in @p Pointerstruct that is pointed to by the pointer of that + * cell. It is clear that interpolate(all_in, all_out) only can be + * called with the vector all_in that previously was + * the parameter of the prepare_for_coarsening_and_refinement(all_in) + * function. Hence interpolate(all_in, all_out) can (in contrast to + * refine_interpolate(in, out)) only be called once. + *
+ * + * + *

Interaction with hanging nodes

+ * + * 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 onto the new finite element space that contains + * constraints from hanging nodes. + * + * 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(). + * + * + *

Implementation in the context of hp-finite elements

+ * + * 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>> - dof_values_on_cell; -}; + template , int spacedim = dim> + class SolutionTransfer + { + public: + /** + * Constructor, takes the current DoFHandler as argument. + */ + SolutionTransfer(const DoFHandler &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 &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. in.size()==n_dofs_old, + * out.size()==n_dofs_refined) + * + * 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 + * prepare_for_refinement_and_coarsening(all_in). 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 &all_in, + std::vector &all_out) const; + + /** + * Same as the previous function. It interpolates only one function. It + * assumes the vectors having the right sizes (i.e. + * in.size()==n_dofs_old, out.size()==n_dofs_refined) + * + * Multiple calling of this function is NOT allowed. Interpolating several + * functions can be performed in one step by using interpolate (all_in, + * all_out) + */ + 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, + SolutionTransfer> + 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> 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 vector indices + * (if the cell will be refined) or the vector dof_values + * (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 *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> *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 *indices_ptr; + std::vector> *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, 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>> + dof_values_on_cell; + }; +} // namespace Legacy DEAL_II_NAMESPACE_CLOSE diff --git a/source/distributed/CMakeLists.txt b/source/distributed/CMakeLists.txt index 252838c5a8..5eb431c6d3 100644 --- a/source/distributed/CMakeLists.txt +++ b/source/distributed/CMakeLists.txt @@ -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 diff --git a/source/distributed/fully_distributed_tria.cc b/source/distributed/fully_distributed_tria.cc index a7f3fae58f..0d96cf982a 100644 --- a/source/distributed/fully_distributed_tria.cc +++ b/source/distributed/fully_distributed_tria.cc @@ -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:: + 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::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 index 250c000b48..0000000000 --- a/source/distributed/solution_transfer.cc +++ /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 - -#ifdef DEAL_II_WITH_P4EST - -# include -# include - -# include -# include - -# include -# include - -# include -# include -# include -# include -# include -# include -# include -# include - -# include -# include - - -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 - std::vector - pack_dof_values(std::vector> &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 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 - std::vector> - unpack_dof_values( - const boost::iterator_range::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> unpacked_data; - unpacked_data.reserve(n_elements); - for (unsigned int i = 0; i < n_elements; ++i) - { - Vector 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 - SolutionTransfer::SolutionTransfer( - const DoFHandler &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 *>( - &dof_handler->get_triangulation()) != nullptr), - ExcMessage( - "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object.")); - } - - - - template - void - SolutionTransfer:: - prepare_for_coarsening_and_refinement( - const std::vector &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 - void - SolutionTransfer::register_data_attach() - { - // TODO: casting away constness is bad - parallel::DistributedTriangulationBase *tria = - (dynamic_cast *>( - const_cast *>( - &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::cell_iterator &cell_, - const CellStatus status) { - return this->pack_callback(cell_, status); - }, - /*returns_variable_size_data=*/dof_handler->has_hp_capabilities()); - } - - - - template - void - SolutionTransfer:: - prepare_for_coarsening_and_refinement(const VectorType &in) - { - std::vector all_in(1, &in); - prepare_for_coarsening_and_refinement(all_in); - } - - - - template - void - SolutionTransfer::prepare_for_serialization( - const VectorType &in) - { - std::vector all_in(1, &in); - prepare_for_serialization(all_in); - } - - - - template - void - SolutionTransfer::prepare_for_serialization( - const std::vector &all_in) - { - prepare_for_coarsening_and_refinement(all_in); - } - - - - template - void - SolutionTransfer::deserialize(VectorType &in) - { - std::vector all_in(1, &in); - deserialize(all_in); - } - - - - template - void - SolutionTransfer::deserialize( - std::vector &all_in) - { - register_data_attach(); - - // this makes interpolate() happy - input_vectors.resize(all_in.size()); - - interpolate(all_in); - } - - - template - void - SolutionTransfer::interpolate( - std::vector &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 *tria = - (dynamic_cast *>( - const_cast *>( - &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::cell_iterator &cell_, - const CellStatus status, - const boost::iterator_range::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(valence[i]) == Number() ? - Number() : - (Number(1.0) / static_cast(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 - void - SolutionTransfer::interpolate(VectorType &out) - { - std::vector all_out(1, &out); - interpolate(all_out); - } - - - - template - std::vector - SolutionTransfer::pack_callback( - const typename Triangulation::cell_iterator &cell_, - const CellStatus status) - { - typename DoFHandler::cell_iterator cell(*cell_, - dof_handler); - - // create buffer for each individual object - std::vector<::dealii::Vector> 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(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(); // 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(dof_values, - dofs_per_cell); - } - - - - template - void - SolutionTransfer::unpack_callback( - const typename Triangulation::cell_iterator &cell_, - const CellStatus status, - const boost::iterator_range::const_iterator> - &data_range, - std::vector &all_out, - VectorType &valence) - { - typename DoFHandler::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> - dof_values = - unpack_dof_values(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 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 index 931ddac307..0000000000 --- a/source/distributed/solution_transfer.inst.in +++ /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_space_dimension>; - template class SolutionTransfer< - deal_II_dimension, - ::dealii::LinearAlgebra::distributed::Vector, - deal_II_space_dimension>; - template class SolutionTransfer< - deal_II_dimension, - ::dealii::LinearAlgebra::distributed::Vector, - deal_II_space_dimension>; - template class SolutionTransfer< - deal_II_dimension, - ::dealii::LinearAlgebra::distributed::BlockVector, - deal_II_space_dimension>; - template class SolutionTransfer< - deal_II_dimension, - ::dealii::LinearAlgebra::distributed::BlockVector, - deal_II_space_dimension>; - - -# ifdef DEAL_II_WITH_PETSC - template class SolutionTransfer; - template class SolutionTransfer; -# endif - -# ifdef DEAL_II_WITH_TRILINOS - template class SolutionTransfer; - template class SolutionTransfer; -# endif - -#endif - \} - \} - } diff --git a/source/grid/tria.cc b/source/grid/tria.cc index b1e5adec70..57e891dbc5 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -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::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(sizes_fixed_cumulative.data()), sizes_fixed_cumulative.size() * sizeof(unsigned int)); @@ -13700,14 +13698,19 @@ DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim)) void Triangulation::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::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::version_number @@ -13726,21 +13729,24 @@ DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim)) void Triangulation::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::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::update_cell_relations() +template +DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim)) +void Triangulation::pack_data_serial() +{ + if (dynamic_cast *>( + this)) + return; + + std::vector 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:: + 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:: + 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:: + 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 +DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim)) +void Triangulation::unpack_data_serial() +{ + if (dynamic_cast *>( + this)) + return; + + // transfer data after triangulation got updated + if (this->cell_attached_data.n_attached_data_sets > 0) + { + std::vector:: + 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 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim)) void Triangulation::execute_coarsening_and_refinement() @@ -15672,10 +15833,17 @@ void Triangulation::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::save_attached_data( // cast away constness auto tria = const_cast *>(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 diff --git a/source/numerics/solution_transfer.cc b/source/numerics/solution_transfer.cc index 81c87c503e..3f5e26856f 100644 --- a/source/numerics/solution_transfer.cc +++ b/source/numerics/solution_transfer.cc @@ -12,6 +12,8 @@ // // ------------------------------------------------------------------------ +#include + #include #include @@ -19,6 +21,7 @@ #include #include +#include #include @@ -38,601 +41,1081 @@ #include +#include +#include + 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 + std::vector + pack_dof_values(std::vector> &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 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 + std::vector> + unpack_dof_values( + const boost::iterator_range::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> unpacked_data; + unpacked_data.reserve(n_elements); + for (unsigned int i = 0; i < n_elements; ++i) + { + Vector 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 SolutionTransfer::SolutionTransfer( - const DoFHandler &dof) + const DoFHandler &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 +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement( + const std::vector &all_in) { - Assert( - (dynamic_cast *>( - &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 + 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 -SolutionTransfer::~SolutionTransfer() +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement(const std::vector &all_in) { - clear(); + std::vector 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 void -SolutionTransfer::clear() +SolutionTransfer::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 *>( + &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::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 +void +SolutionTransfer:: + prepare_for_coarsening_and_refinement(const VectorType &in) +{ + std::vector all_in(1, &in); + prepare_for_coarsening_and_refinement(all_in); } template void -SolutionTransfer::prepare_for_pure_refinement() +SolutionTransfer::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 - 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>(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 all_in(1, &in); + prepare_for_serialization(all_in); } template void -SolutionTransfer::refine_interpolate( - const VectorType &in, - VectorType &out) const +SolutionTransfer::prepare_for_serialization( + const std::vector &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 - subdomain_modifier(dof_handler->get_triangulation()); - - Vector local_values(0); - - typename std::map, - 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::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 +void +SolutionTransfer::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 - void - extract_interpolation_matrices(const DoFHandler &dof, - dealii::Table<2, FullMatrix> &matrices) - { - if (dof.has_hp_capabilities() == false) - return; - - const dealii::hp::FECollection &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:: - ExcInterpolationNotImplemented &) - { - matrices(i, j).reinit(0, 0); - } - } - } + std::vector all_in(1, &in); + deserialize(all_in); +} - template - void - restriction_additive(const FiniteElement &, - std::vector> &) - {} - template - void - restriction_additive(const dealii::hp::FECollection &fe, - std::vector> &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 +void +SolutionTransfer::deserialize( + std::vector &all_in) +{ + register_data_attach(); + // this makes interpolate() happy + input_vectors.resize(all_in.size()); + + interpolate(all_in); +} template void -SolutionTransfer:: - prepare_for_coarsening_and_refinement(const std::vector &all_in) +SolutionTransfer::interpolate( + std::vector &all_out) { - Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); - Assert(prepared_for != coarsening_and_refinement, - ExcAlreadyPrepForCoarseAndRef()); + const dealii::internal::parallel::shared:: + TemporarilyRestoreSubdomainIds + 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 *>( + &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::cell_iterator &cell_, + const CellStatus status, + const boost::iterator_range::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(valence[i]) == Number() ? + Number() : + (Number(1.0) / static_cast(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 - 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>(n_cells_to_stay_or_refine) - .swap(indices_on_cell); - - std::vector>>( - n_coarsen_fathers, - std::vector>(in_size)) - .swap(dof_values_on_cell); - - Table<2, FullMatrix> interpolation_hp; - std::vector> 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 +void +SolutionTransfer::interpolate( + std::vector &all_out) +{ + std::vector temp(all_out.size()); + + for (std::size_t i = 0; i < temp.size(); ++i) + temp[i] = &(all_out[i]); + + this->interpolate(temp); +} + + + +template +void +SolutionTransfer::interpolate(VectorType &out) +{ + std::vector all_out(1, &out); + interpolate(all_out); +} + + + +template +std::vector +SolutionTransfer::pack_callback( + const typename Triangulation::cell_iterator &cell_, + const CellStatus status) +{ + typename DoFHandler::cell_iterator cell(*cell_, dof_handler); + + // create buffer for each individual object + std::vector<::dealii::Vector> 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 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(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>( - in_size, Vector(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(); // 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(dof_values, + dofs_per_cell); } template void -SolutionTransfer:: - prepare_for_coarsening_and_refinement(const VectorType &in) +SolutionTransfer::unpack_callback( + const typename Triangulation::cell_iterator &cell_, + const CellStatus status, + const boost::iterator_range::const_iterator> &data_range, + std::vector &all_out, + VectorType &valence) { - std::vector all_in(1, in); - prepare_for_coarsening_and_refinement(all_in); + typename DoFHandler::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> + dof_values = + unpack_dof_values(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 ones(dofs_per_cell); + ones = 1.0; + cell->distribute_local_to_global_by_interpolation(ones, + valence, + fe_index); + } } template void -SolutionTransfer::interpolate( - const std::vector &all_in, - std::vector &all_out) const +SolutionTransfer::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 - subdomain_modifier(dof_handler->get_triangulation()); - - Vector local_values; - std::vector dofs; - - typename std::map, - Pointerstruct>::const_iterator pointerstruct, - cell_map_end = cell_map.end(); - - Table<2, FullMatrix> interpolation_hp; - internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); - Vector 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 *const indexptr = - pointerstruct->second.indices_ptr; - const std::vector> - *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 + SolutionTransfer::SolutionTransfer( + const DoFHandler &dof) + : dof_handler(&dof, typeid(*this).name()) + , n_dofs_old(0) + , prepared_for(none) + { + Assert((dynamic_cast< + const parallel::distributed::Triangulation *>( + &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 + SolutionTransfer::~SolutionTransfer() + { + clear(); + } + + + + template + void + SolutionTransfer::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 + void + SolutionTransfer::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 + 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>(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 + void + SolutionTransfer::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 + subdomain_modifier(dof_handler->get_triangulation()); + + Vector local_values(0); + + typename std::map, + 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::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 + void + extract_interpolation_matrices( + const DoFHandler &dof, + dealii::Table<2, FullMatrix> &matrices) + { + if (dof.has_hp_capabilities() == false) + return; + + const dealii::hp::FECollection &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::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:: + 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 *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 &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::set((*data)(i), - dofs[i], - all_out[j]); + matrices(i, j).reinit(0, 0); } } - // undefined status - else - DEAL_II_ASSERT_UNREACHABLE(); + } + + + template + void + restriction_additive(const FiniteElement &, + std::vector> &) + {} + + template + void + restriction_additive( + const dealii::hp::FECollection &fe, + std::vector> &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 + void + SolutionTransfer:: + prepare_for_coarsening_and_refinement(const std::vector &all_in) + { + Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef()); + Assert(prepared_for != coarsening_and_refinement, + ExcAlreadyPrepForCoarseAndRef()); -template -void -SolutionTransfer::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 all_in = {in}; - std::vector 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 + 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>(n_cells_to_stay_or_refine) + .swap(indices_on_cell); + + std::vector>>( + n_coarsen_fathers, + std::vector>(in_size)) + .swap(dof_values_on_cell); + + Table<2, FullMatrix> interpolation_hp; + std::vector> 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 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>( + in_size, Vector(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 -std::size_t -SolutionTransfer::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 + void + SolutionTransfer:: + prepare_for_coarsening_and_refinement(const VectorType &in) + { + std::vector all_in(1, in); + prepare_for_coarsening_and_refinement(all_in); + } -template -std::size_t -SolutionTransfer::Pointerstruct::memory_consumption() - const -{ - return sizeof(*this); -} + + + template + void + SolutionTransfer::interpolate( + const std::vector &all_in, + std::vector &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 + subdomain_modifier(dof_handler->get_triangulation()); + + Vector local_values; + std::vector dofs; + + typename std::map, + Pointerstruct>::const_iterator pointerstruct, + cell_map_end = cell_map.end(); + + Table<2, FullMatrix> interpolation_hp; + internal::extract_interpolation_matrices(*dof_handler, interpolation_hp); + Vector 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 *const indexptr = + pointerstruct->second.indices_ptr; + + const std::vector> + *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::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 *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 &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::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 + void + SolutionTransfer::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 all_in = {in}; + std::vector all_out = {out}; + + interpolate(all_in, all_out); + + out = all_out[0]; + } + + + + template + std::size_t + SolutionTransfer::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 + std::size_t + SolutionTransfer::Pointerstruct:: + memory_consumption() const + { + return sizeof(*this); + } + +} // namespace Legacy /*-------------- Explicit Instantiations -------------------------------*/ diff --git a/source/numerics/solution_transfer.inst.in b/source/numerics/solution_transfer.inst.in index a53d7569c8..ba92230013 100644 --- a/source/numerics/solution_transfer.inst.in +++ b/source/numerics/solution_transfer.inst.in @@ -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; + template class SolutionTransfer; diff --git a/tests/aniso/solution_transfer.cc b/tests/aniso/solution_transfer.cc index b3e9c50c79..1ad6174712 100644 --- a/tests/aniso/solution_transfer.cc +++ b/tests/aniso/solution_transfer.cc @@ -93,12 +93,12 @@ transfer(std::ostream &out) cell->set_refine_flag(RefinementCase::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 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"); diff --git a/tests/bits/solution_transfer.cc b/tests/bits/solution_transfer.cc index 468b44c728..7b47a79f9f 100644 --- a/tests/bits/solution_transfer.cc +++ b/tests/bits/solution_transfer.cc @@ -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 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 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"); diff --git a/tests/codim_one/solution_transfer_01.cc b/tests/codim_one/solution_transfer_01.cc index 45d1f2ff86..6919afcf27 100644 --- a/tests/codim_one/solution_transfer_01.cc +++ b/tests/codim_one/solution_transfer_01.cc @@ -87,7 +87,7 @@ main() // get the interpolated solution // back Vector 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) diff --git a/tests/distributed_grids/solution_transfer_01.cc b/tests/distributed_grids/solution_transfer_01.cc index 3dbcfa5c7d..237ea34e96 100644 --- a/tests/distributed_grids/solution_transfer_01.cc +++ b/tests/distributed_grids/solution_transfer_01.cc @@ -33,9 +33,6 @@ #include "../tests.h" -#include "coarse_grid_common.h" - - template void diff --git a/tests/fullydistributed_grids/solution_transfer_01.cc b/tests/fullydistributed_grids/solution_transfer_01.cc index 547fbea36e..f9027add34 100644 --- a/tests/fullydistributed_grids/solution_transfer_01.cc +++ b/tests/fullydistributed_grids/solution_transfer_01.cc @@ -18,6 +18,7 @@ #include #include +#include #include #include diff --git a/tests/grid/save_load_01.cc b/tests/grid/save_load_01.cc index fa0bb31eb3..e3b6f2f2d8 100644 --- a/tests/grid/save_load_01.cc +++ b/tests/grid/save_load_01.cc @@ -21,13 +21,8 @@ #include #include -#include - #include #include -#include - -#include #include "./tests.h" diff --git a/tests/grid/solution_transfer_01.cc b/tests/grid/solution_transfer_01.cc new file mode 100644 index 0000000000..c6b1c77267 --- /dev/null +++ b/tests/grid/solution_transfer_01.cc @@ -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 +#include + +#include + +#include +#include +#include + +#include + +#include +#include + +#include "./tests.h" + + +template +class InterpolationFunction : public Function +{ +public: + InterpolationFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + return p.norm(); + } +}; + +template +void +test(TriangulationType &triangulation) +{ + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(FE_Q(2)); + + using VectorType = Vector; + + VectorType vector(dof_handler.n_dofs()); + + VectorTools::interpolate(dof_handler, InterpolationFunction(), vector); + + VectorType vector_loaded(dof_handler.n_dofs()); + + const std::string filename = + "solution_transfer_" + std::to_string(dim) + "d_out"; + + { + SolutionTransfer solution_transfer(dof_handler); + solution_transfer.prepare_for_serialization(vector); + + triangulation.save(filename); + } + + triangulation.clear(); + + { + triangulation.load(filename); + dof_handler.distribute_dofs(FE_Q(2)); + + SolutionTransfer 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 triangulation; + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + test(triangulation); + } + deallog.pop(); + + deallog.push("3d"); + { + constexpr int dim = 3; + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + test(triangulation); + } + deallog.pop(); +} diff --git a/tests/grid/solution_transfer_01.output b/tests/grid/solution_transfer_01.output new file mode 100644 index 0000000000..7e803b74e4 --- /dev/null +++ b/tests/grid/solution_transfer_01.output @@ -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 index 0000000000..eebb50ce0f --- /dev/null +++ b/tests/grid/solution_transfer_02.cc @@ -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 +#include + +#include + +#include +#include +#include + +#include + +#include +#include + +#include "./tests.h" + + +template +class InterpolationFunction : public Function +{ +public: + InterpolationFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + return p[0]; + } +}; + +template +void +test(const unsigned int type) +{ + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + + const FE_Q fe(2); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + using VectorType = Vector; + + VectorType vector(dof_handler.n_dofs()); + VectorTools::interpolate(dof_handler, InterpolationFunction(), vector); + + SolutionTransfer 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(), 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(0); + test(1); + test(2); + } + deallog.pop(); + + deallog.push("3d"); + { + constexpr int dim = 3; + + test(0); + test(1); + test(2); + } + deallog.pop(); +} diff --git a/tests/grid/solution_transfer_02.output b/tests/grid/solution_transfer_02.output new file mode 100644 index 0000000000..da1d891e2a --- /dev/null +++ b/tests/grid/solution_transfer_02.output @@ -0,0 +1,7 @@ + +DEAL:2d::PASSED +DEAL:2d::PASSED +DEAL:2d::PASSED +DEAL:3d::PASSED +DEAL:3d::PASSED +DEAL:3d::PASSED diff --git a/tests/hp/do_function_derivatives_01.cc b/tests/hp/do_function_derivatives_01.cc index e748d4694d..05c0c74688 100644 --- a/tests/hp/do_function_derivatives_01.cc +++ b/tests/hp/do_function_derivatives_01.cc @@ -86,7 +86,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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>()); diff --git a/tests/hp/do_function_hessians_01.cc b/tests/hp/do_function_hessians_01.cc index b27980d44b..96b361bfcc 100644 --- a/tests/hp/do_function_hessians_01.cc +++ b/tests/hp/do_function_hessians_01.cc @@ -86,7 +86,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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>()); diff --git a/tests/hp/do_function_laplacians_01.cc b/tests/hp/do_function_laplacians_01.cc index a6d196df9f..800d217b72 100644 --- a/tests/hp/do_function_laplacians_01.cc +++ b/tests/hp/do_function_laplacians_01.cc @@ -86,7 +86,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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>()); diff --git a/tests/hp/fe_nothing_17.cc b/tests/hp/fe_nothing_17.cc index 7dbc05b7b4..8bb35257bc 100644 --- a/tests/hp/fe_nothing_17.cc +++ b/tests/hp/fe_nothing_17.cc @@ -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 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; } diff --git a/tests/hp/laplace.h b/tests/hp/laplace.h index 65b7010e3e..d5a0bead53 100644 --- a/tests/hp/laplace.h +++ b/tests/hp/laplace.h @@ -463,7 +463,7 @@ Laplace::refine_grid(const unsigned int cycle) setup_system(); // 3.7. Solution Transfer finish - soltrans.interpolate(solution_coarse, solution); + soltrans.interpolate(solution); } diff --git a/tests/hp/solution_transfer.cc b/tests/hp/solution_transfer.cc index ad13548bbc..e86cc1ef97 100644 --- a/tests/hp/solution_transfer.cc +++ b/tests/hp/solution_transfer.cc @@ -151,12 +151,6 @@ transfer(std::ostream &out) cell->set_coarsen_flag(); } - Vector 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::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 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. diff --git a/tests/hp/solution_transfer_02.cc b/tests/hp/solution_transfer_02.cc index 622928532d..479898f4cd 100644 --- a/tests/hp/solution_transfer_02.cc +++ b/tests/hp/solution_transfer_02.cc @@ -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); } diff --git a/tests/hp/solution_transfer_03.cc b/tests/hp/solution_transfer_03.cc index 75edf39d6b..b59dd1ce3e 100644 --- a/tests/hp/solution_transfer_03.cc +++ b/tests/hp/solution_transfer_03.cc @@ -92,19 +92,21 @@ main() // Interpolate solution SolutionTransfer<2, Vector> 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 new_solution(dof_handler.n_dofs()); - solultion_trans.interpolate(solution, new_solution); + solultion_trans.interpolate(new_solution); // Save output diff --git a/tests/hp/solution_transfer_04.cc b/tests/hp/solution_transfer_04.cc index bb0c53d6de..e214319cb0 100644 --- a/tests/hp/solution_transfer_04.cc +++ b/tests/hp/solution_transfer_04.cc @@ -85,7 +85,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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; diff --git a/tests/hp/solution_transfer_05.cc b/tests/hp/solution_transfer_05.cc index a08d922c3c..2aca704e7c 100644 --- a/tests/hp/solution_transfer_05.cc +++ b/tests/hp/solution_transfer_05.cc @@ -85,7 +85,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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; diff --git a/tests/hp/solution_transfer_06.cc b/tests/hp/solution_transfer_06.cc index 3350ff158e..ee8647e544 100644 --- a/tests/hp/solution_transfer_06.cc +++ b/tests/hp/solution_transfer_06.cc @@ -65,18 +65,19 @@ test() // set refine flag for the only cell we have, then do the refinement SolutionTransfer> 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 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) diff --git a/tests/hp/solution_transfer_07.cc b/tests/hp/solution_transfer_07.cc index 08dbaa7258..ee8647e544 100644 --- a/tests/hp/solution_transfer_07.cc +++ b/tests/hp/solution_transfer_07.cc @@ -65,18 +65,19 @@ test() // set refine flag for the only cell we have, then do the refinement SolutionTransfer> 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 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) diff --git a/tests/hp/solution_transfer_08.cc b/tests/hp/solution_transfer_08.cc index 9321d9e9c8..2f368b6262 100644 --- a/tests/hp/solution_transfer_08.cc +++ b/tests/hp/solution_transfer_08.cc @@ -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 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) diff --git a/tests/hp/solution_transfer_09.cc b/tests/hp/solution_transfer_09.cc index 4695ed3454..9a21fa9209 100644 --- a/tests/hp/solution_transfer_09.cc +++ b/tests/hp/solution_transfer_09.cc @@ -63,7 +63,7 @@ test() // set refine flag for the only cell we have, then do the refinement SolutionTransfer> 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 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) diff --git a/tests/hp/solution_transfer_10.cc b/tests/hp/solution_transfer_10.cc index 7f0abd91df..c9b84f0195 100644 --- a/tests/hp/solution_transfer_10.cc +++ b/tests/hp/solution_transfer_10.cc @@ -105,7 +105,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector new_solution(dof_handler.n_dofs()); - solultion_trans.interpolate(solution, new_solution); + solultion_trans.interpolate(new_solution); deallog << "OK" << std::endl; } diff --git a/tests/hp/solution_transfer_11.cc b/tests/hp/solution_transfer_11.cc index 678b001018..9dc5efc907 100644 --- a/tests/hp/solution_transfer_11.cc +++ b/tests/hp/solution_transfer_11.cc @@ -74,11 +74,11 @@ transfer(std::ostream &out) Vector 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); } diff --git a/tests/hp/solution_transfer_12.cc b/tests/hp/solution_transfer_12.cc index a2ee19e2aa..9473c6e57b 100644 --- a/tests/hp/solution_transfer_12.cc +++ b/tests/hp/solution_transfer_12.cc @@ -177,11 +177,6 @@ transfer(std::ostream &out) cell->set_coarsen_flag(); } - Vector 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::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 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. diff --git a/tests/hp/solution_transfer_13.cc b/tests/hp/solution_transfer_13.cc index b0a79b637a..66f7835107 100644 --- a/tests/hp/solution_transfer_13.cc +++ b/tests/hp/solution_transfer_13.cc @@ -132,7 +132,7 @@ main() dof_handler.distribute_dofs(fe_collection); Vector 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 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; diff --git a/tests/hp/solution_transfer_14.cc b/tests/hp/solution_transfer_14.cc index f197a50a06..99afd752a7 100644 --- a/tests/hp/solution_transfer_14.cc +++ b/tests/hp/solution_transfer_14.cc @@ -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); } diff --git a/tests/hp/solution_transfer_15.cc b/tests/hp/solution_transfer_15.cc index a220af6bff..05e89b180c 100644 --- a/tests/hp/solution_transfer_15.cc +++ b/tests/hp/solution_transfer_15.cc @@ -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); } diff --git a/tests/hp/solution_transfer_16.cc b/tests/hp/solution_transfer_16.cc index 73f4a7e468..d3c3968f54 100644 --- a/tests/hp/solution_transfer_16.cc +++ b/tests/hp/solution_transfer_16.cc @@ -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 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()); diff --git a/tests/mpi/solution_transfer_02.cc b/tests/mpi/solution_transfer_02.cc index 3da6d79d61..ec31643421 100644 --- a/tests/mpi/solution_transfer_02.cc +++ b/tests/mpi/solution_transfer_02.cc @@ -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); } diff --git a/tests/mpi/solution_transfer_03.cc b/tests/mpi/solution_transfer_03.cc index 82db17533f..e2a3ba15c6 100644 --- a/tests/mpi/solution_transfer_03.cc +++ b/tests/mpi/solution_transfer_03.cc @@ -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); } diff --git a/tests/mpi/solution_transfer_07.cc b/tests/mpi/solution_transfer_07.cc index 9905d62448..7c60b4c72d 100644 --- a/tests/mpi/solution_transfer_07.cc +++ b/tests/mpi/solution_transfer_07.cc @@ -62,7 +62,8 @@ test(const bool allow_artificial_cells) dh.distribute_dofs(fe); Vector 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()); diff --git a/tests/mpi/solution_transfer_08.cc b/tests/mpi/solution_transfer_08.cc index 2f1dffded2..ca9d957808 100644 --- a/tests/mpi/solution_transfer_08.cc +++ b/tests/mpi/solution_transfer_08.cc @@ -56,13 +56,13 @@ test(const bool allow_artificial_cells) sol_old = 1.; SolutionTransfer 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 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()); diff --git a/tests/mpi/solution_transfer_09.cc b/tests/mpi/solution_transfer_09.cc index 1fb9607b94..bbcad94f9b 100644 --- a/tests/mpi/solution_transfer_09.cc +++ b/tests/mpi/solution_transfer_09.cc @@ -61,7 +61,7 @@ refine_and_transfer(const Function &function, Triangulation &tria) dh.distribute_dofs(fe); Vector sol_new(dh.n_dofs()); - soltrans.interpolate(sol_old, sol_new); + soltrans.interpolate(sol_new); return sol_new; } diff --git a/tests/mpi/solution_transfer_10.cc b/tests/mpi/solution_transfer_10.cc index b8a0c090fd..52e71b4f7a 100644 --- a/tests/mpi/solution_transfer_10.cc +++ b/tests/mpi/solution_transfer_10.cc @@ -55,13 +55,13 @@ refine_and_transfer(const Function &function, Triangulation &tria) VectorTools::interpolate(MappingQ(1), dh, function, sol_old); SolutionTransfer 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 sol_new(dh.n_dofs()); - soltrans.refine_interpolate(sol_old, sol_new); + soltrans.interpolate(sol_new); return sol_new; } diff --git a/tests/petsc/solution_transfer_01.cc b/tests/petsc/solution_transfer_01.cc index 6c03f99680..495f020308 100644 --- a/tests/petsc/solution_transfer_01.cc +++ b/tests/petsc/solution_transfer_01.cc @@ -69,7 +69,7 @@ test() triangulation.execute_coarsening_and_refinement(); - solution_trans.interpolate(previous_solution, solution); + solution_trans.interpolate(solution); deallog << "OK" << std::endl; } diff --git a/tests/petsc/solution_transfer_02.cc b/tests/petsc/solution_transfer_02.cc index 6c03f99680..495f020308 100644 --- a/tests/petsc/solution_transfer_02.cc +++ b/tests/petsc/solution_transfer_02.cc @@ -69,7 +69,7 @@ test() triangulation.execute_coarsening_and_refinement(); - solution_trans.interpolate(previous_solution, solution); + solution_trans.interpolate(solution); deallog << "OK" << std::endl; } diff --git a/tests/simplex/step-31.cc b/tests/simplex/step-31.cc index c3dfdd47bc..d732e30b96 100644 --- a/tests/simplex/step-31.cc +++ b/tests/simplex/step-31.cc @@ -1013,12 +1013,12 @@ namespace Step31 std::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; diff --git a/tests/sundials/step-77.cc b/tests/sundials/step-77.cc index 859d9c6cce..9b0b0dcec1 100644 --- a/tests/sundials/step-77.cc +++ b/tests/sundials/step-77.cc @@ -339,7 +339,7 @@ namespace Step77 dof_handler.distribute_dofs(fe); Vector tmp(dof_handler.n_dofs()); - solution_transfer.interpolate(current_solution, tmp); + solution_transfer.interpolate(tmp); current_solution = std::move(tmp); hanging_node_constraints.clear(); diff --git a/tests/trilinos/step-77-with-nox.cc b/tests/trilinos/step-77-with-nox.cc index bb2d05ee38..737c298aec 100644 --- a/tests/trilinos/step-77-with-nox.cc +++ b/tests/trilinos/step-77-with-nox.cc @@ -345,7 +345,7 @@ namespace Step77 dof_handler.distribute_dofs(fe); Vector tmp(dof_handler.n_dofs()); - solution_transfer.interpolate(current_solution, tmp); + solution_transfer.interpolate(tmp); current_solution = std::move(tmp); hanging_node_constraints.clear(); -- 2.39.5