--- /dev/null
+Deprecated: The classes SolutionTransfer and parallel::distributed::SolutionTransfer
+have been unified to SolutionTransfer. The class now supports both serial and parallel
+meshes. The class has lost some functions: the function SolutionTransfer::interpolate()
+that takes the input vector as well as the less frequently used functions
+SolutionTransfer::prepare_for_pure_refinement() and
+SolutionTransfer::refine_interpolate(). Please use the other functions to accomplish the
+same functionality. For the time being, the old implementation has been moved to the
+`Legacy` namespace. The old parallel::distributed::SolutionTransfer has been now early deprecated.
+<br>
+(Pasquale Claudio Africa, Bruno Blais, Peter Munch, 2024/01/03)
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
triangulation.execute_coarsening_and_refinement();
setup_system();
- solution_trans.interpolate(previous_solution, solution);
+ solution_trans.interpolate(solution);
constraints.distribute(solution);
}
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:
std::vector<TrilinosWrappers::MPI::Vector> tmp = {
TrilinosWrappers::MPI::Vector(temperature_solution),
TrilinosWrappers::MPI::Vector(temperature_solution)};
- temperature_trans.interpolate(x_temperature, tmp);
+ temperature_trans.interpolate(tmp);
temperature_solution = tmp[0];
old_temperature_solution = tmp[1];
// 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);
std::vector<Vector<double>> transfer_out = {
Vector<double>(dof_handler.n_dofs()),
Vector<double>(dof_handler.n_dofs())};
- soltrans.interpolate(transfer_in, transfer_out);
+ soltrans.interpolate(transfer_out);
old_solution = std::move(transfer_out[0]);
predictor = std::move(transfer_out[1]);
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];
std::vector<TrilinosWrappers::MPI::BlockVector> tmp_darcy(2);
tmp_darcy[0].reinit(darcy_solution);
tmp_darcy[1].reinit(darcy_solution);
- darcy_soltrans.interpolate(x_darcy, tmp_darcy);
+ darcy_soltrans.interpolate(tmp_darcy);
last_computed_darcy_solution = tmp_darcy[0];
second_last_computed_darcy_solution = tmp_darcy[1];
// 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
setup_system();
- solution_transfer.interpolate(coarse_solution, current_solution);
+ solution_transfer.interpolate(current_solution);
nonzero_constraints.distribute(current_solution);
}
setup_system();
- solution_transfer.interpolate(coarse_solution, current_solution);
-
+ solution_transfer.interpolate(current_solution);
nonzero_constraints.distribute(current_solution);
}
#include <deal.II/base/config.h>
-#include <deal.II/distributed/tria.h>
-
-#include <deal.II/dofs/dof_handler.h>
-
-#include <vector>
+#include <deal.II/numerics/solution_transfer.h>
DEAL_II_NAMESPACE_OPEN
namespace distributed
{
/**
- * Transfer a discrete FE function (like a solution vector) by
- * interpolation while refining and/or coarsening a distributed grid and
- * handles the necessary communication.
- *
- * @note It is important to note, that if you use more than one
- * SolutionTransfer object at the same time, that the calls to prepare_*()
- * and interpolate()/deserialize() need to be in the same order.
- *
- * <h3>Note on ghost elements</h3> In a parallel computation PETSc or
- * Trilinos vector may contain ghost elements or not. For reading in
- * information with prepare_for_coarsening_and_refinement() or
- * prepare_for_serialization() you need to supply vectors with ghost
- * elements, so that all locally_active elements can be read. On the other
- * hand, ghosted vectors are generally not writable, so for calls to
- * interpolate() or deserialize() you need to supply distributed vectors
- * without ghost elements. More precisely, during interpolation the
- * current algorithm writes into all locally active degrees of freedom.
- *
- * <h3>Transferring a solution</h3> Here VectorType is your favorite
- * vector type, e.g. PETScWrappers::MPI::Vector,
- * TrilinosWrappers::MPI::Vector, or corresponding block vectors.
- * @code
- * parallel::distributed::SolutionTransfer<dim, VectorType>
- * soltrans(dof_handler);
- * // flag some cells for refinement and coarsening, e.g.
- * parallel::distributed::GridRefinement::
- * refine_and_coarsen_fixed_fraction(tria,
- * error_indicators,
- * 0.3,
- * 0.05);
- *
- * // prepare the triangulation,
- * tria.prepare_coarsening_and_refinement();
- *
- * // prepare the SolutionTransfer object for coarsening and refinement
- * // and give the solution vector that we intend to interpolate later,
- * soltrans.prepare_for_coarsening_and_refinement(solution);
- *
- * // actually execute the refinement,
- * tria.execute_coarsening_and_refinement();
- *
- * // redistribute dofs,
- * dof_handler.distribute_dofs(fe);
- *
- * // and interpolate the solution
- * VectorType interpolated_solution;
- *
- * //create VectorType in the right size here
- * soltrans.interpolate(interpolated_solution);
- * @endcode
- *
- * As the grid is distributed, it is important to note that the old
- * solution(s) must be copied to one that also provides access to the
- * locally relevant DoF values (these values required for the interpolation
- * process):
- * @code
- * // Create initial indexsets pertaining to the grid before refinement
- * const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
- * const IndexSet locally_relevant_dofs =
- * DoFTools::extract_locally_relevant_dofs(dof_handler);
- *
- * // The solution vector only knows about locally owned DoFs
- * TrilinosWrappers::MPI::Vector solution;
- * solution.reinit(locally_owned_dofs,
- * mpi_communicator);
- * ...
- * // Transfer solution to vector that provides access to
- * // locally relevant DoFs
- * TrilinosWrappers::MPI::Vector old_solution;
- * old_solution.reinit(locally_owned_dofs,
- * locally_relevant_dofs,
- * mpi_communicator);
- * old_solution = solution;
- *
- * // Initialize SolutionTransfer object
- * parallel::distributed::SolutionTransfer<dim, VectorType>
- * soltrans(dof_handler);
- * soltrans.prepare_for_coarsening_and_refinement(old_solution);
- * ...
- * // Refine grid
- * // Recreate locally_owned_dofs and locally_relevant_dofs index sets
- * ...
- * solution.reinit(locally_owned_dofs, mpi_communicator);
- * soltrans.interpolate(solution);
- * @endcode
- *
- * Different from PETSc and Trilinos vectors,
- * LinearAlgebra::distributed::Vector allows writing into ghost elements.
- * For a ghosted vector the interpolation step can be accomplished via
- * @code
- * interpolated_solution.zero_out_ghost_values();
- * soltrans.interpolate(interpolated_solution);
- * interpolated_solution.update_ghost_values();
- * @endcode
- *
- * <h3>Use for Serialization</h3>
- *
- * This class can be used to serialize and later deserialize a distributed
- * mesh with solution vectors to a file. If you use more than one
- * DoFHandler and therefore more than one SolutionTransfer object, they
- * need to be serialized and deserialized in the same order.
- *
- * If vector has the locally relevant DoFs, serialization works as
- * follows:
- * @code
- * parallel::distributed::SolutionTransfer<dim, VectorType>
- * sol_trans(dof_handler);
- * sol_trans.prepare_for_serialization(vector);
- *
- * triangulation.save(filename);
- * @endcode
- * For deserialization the vector needs to be a distributed vector
- * (without ghost elements):
- * @code
- * //[create coarse mesh...]
- * triangulation.load(filename);
- *
- * parallel::distributed::SolutionTransfer<dim, VectorType>
- * sol_trans(dof_handler);
- * sol_trans.deserialize(distributed_vector);
- * @endcode
- *
- *
- * <h3>Note on usage with DoFHandler with hp-capabilities</h3>
- *
- * Since data on DoFHandler objects with hp-capabilities is associated with
- * many different FiniteElement objects, each cell's data has to be
- * processed with its corresponding `future_fe_index`. Further, if
- * refinement is involved, data will be packed on the parent cell with its
- * `future_fe_index` and unpacked later with the same index on its children.
- * If cells get coarsened into one, data will be packed on the children with
- * the least dominant finite element of their common subspace, and unpacked
- * on the parent with this particular finite element (consult
- * hp::FECollection::find_dominated_fe_extended() for more information).
- *
- * Transferring a solution across refinement works exactly like in the
- * non-hp-case. However, when considering serialization, we also have to
- * store the active FE indices in an additional step. A code snippet
- * demonstrating serialization with the
- * parallel::distributed::SolutionTransfer class with DoFHandler objects
- * with hp-capabilities is provided in the following. Here VectorType is
- * your favorite vector type, e.g. PETScWrappers::MPI::Vector,
- * TrilinosWrappers::MPI::Vector, or corresponding block vectors.
- *
- * If vector has the locally relevant DoFs, serialization works as follows:
- * @code
- * parallel::distributed::
- * SolutionTransfer<dim, VectorType, DoFHandler<dim,spacedim>>
- * sol_trans(hp_dof_handler);
- *
- * hp_dof_handler.prepare_for_serialization_of_active_fe_indices();
- * sol_trans.prepare_for_serialization(vector);
- *
- * triangulation.save(filename);
- * @endcode
- *
- * For deserialization the vector needs to be a distributed vector
- * (without ghost elements):
- * @code
- * //[create coarse mesh...]
- * triangulation.load(filename);
- *
- * hp::FECollection<dim,spacedim> fe_collection;
- * //[prepare identical fe_collection...]
- *
- * DoFHandler<dim,spacedim> hp_dof_handler(triangulation);
- * // We need to introduce our dof_handler to the fe_collection
- * // before setting all active FE indices.
- * hp_dof_handler.deserialize_active_fe_indices();
- * hp_dof_handler.distribute_dofs(fe_collection);
- *
- * parallel::distributed::
- * SolutionTransfer<dim,VectorType,DoFHandler<dim,spacedim>>
- * sol_trans(hp_dof_handler);
- * sol_trans.deserialize(distributed_vector);
- * @endcode
- *
- *
- * <h3>Interaction with hanging nodes</h3>
- *
- * In essence, this class implements the same steps as does
- * dealii::SolutionTransfer (though the implementation is entirely
- * separate). Consequently, the same issue with hanging nodes and
- * coarsening can happen with this class as happens with
- * dealii::SolutionTransfer. See there for an extended discussion.
- *
- * @ingroup distributed
+ * @deprecated Use dealii::SolutionTransfer instead.
*/
template <int dim, typename VectorType, int spacedim = dim>
- class SolutionTransfer
- {
- public:
- /**
- * Constructor.
- *
- * @param[in] dof_handler The DoFHandler on which all operations will
- * happen. At the time when this constructor is called, the DoFHandler
- * still points to the Triangulation before the refinement in question
- * happens.
- * @param[in] average_values Average the contribututions to the same
- * DoF coming from different cells. Note: averaging requires an
- * additional communication step, since the valence of the DoF has to be
- * determined.
- */
- SolutionTransfer(const DoFHandler<dim, spacedim> &dof_handler,
- const bool average_values = false);
-
- /**
- * Destructor.
- */
- ~SolutionTransfer() = default;
-
- /**
- * Prepare the current object for coarsening and refinement. It
- * stores the dof indices of each cell and stores the dof values of the
- * vectors in @p all_in in each cell that'll be coarsened. @p all_in
- * includes all vectors that are to be interpolated onto the new
- * (refined and/or coarsened) grid.
- */
- void
- prepare_for_coarsening_and_refinement(
- const std::vector<const VectorType *> &all_in);
-
- /**
- * Same as the previous function but for only one discrete function to be
- * interpolated.
- */
- void
- prepare_for_coarsening_and_refinement(const VectorType &in);
-
- /**
- * Interpolate the data previously stored in this object before the mesh
- * was refined or coarsened onto the current set of cells. Do so for
- * each of the vectors provided to
- * prepare_for_coarsening_and_refinement() and write the result into the
- * given set of vectors.
- */
- void
- interpolate(std::vector<VectorType *> &all_out);
-
- /**
- * Same as the previous function. It interpolates only one function. It
- * assumes the vectors having the right sizes (i.e.
- * <tt>in.size()==n_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
- *
- * Multiple calling of this function is NOT allowed. Interpolating
- * several functions can be performed in one step by using
- * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> and
- * <tt>interpolate(all_out)</tt>.
- */
- void
- interpolate(VectorType &out);
-
- /**
- * Prepare the serialization of the given vector. The serialization is
- * done by Triangulation::save(). The given vector needs all information
- * on the locally active DoFs (it must be ghosted). See documentation of
- * this class for more information.
- */
- void
- prepare_for_serialization(const VectorType &in);
-
- /**
- * Same as the function above, only for a list of vectors.
- */
- void
- prepare_for_serialization(const std::vector<const VectorType *> &all_in);
-
- /**
- * Execute the deserialization of the given vector. This needs to be
- * done after calling Triangulation::load(). The given vector must be a
- * fully distributed vector without ghost elements. See documentation of
- * this class for more information.
- */
- void
- deserialize(VectorType &in);
-
-
- /**
- * Same as the function above, only for a list of vectors.
- */
- void
- deserialize(std::vector<VectorType *> &all_in);
-
- private:
- /**
- * Pointer to the degree of freedom handler to work with.
- */
- SmartPointer<const DoFHandler<dim, spacedim>,
- SolutionTransfer<dim, VectorType, spacedim>>
- dof_handler;
-
- /**
- * Flag indicating if averaging should be performed.
- */
- const bool average_values;
-
- /**
- * A vector that stores pointers to all the vectors we are supposed to
- * copy over from the old to the new mesh.
- */
- std::vector<const VectorType *> input_vectors;
-
- /**
- * The handle that the Triangulation has assigned to this object
- * with which we can access our memory offset and our pack function.
- */
- unsigned int handle;
-
- /**
- * A callback function used to pack the data on the current mesh into
- * objects that can later be retrieved after refinement, coarsening and
- * repartitioning.
- */
- std::vector<char>
- pack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const CellStatus status);
-
- /**
- * A callback function used to unpack the data on the current mesh that
- * has been packed up previously on the mesh before refinement,
- * coarsening and repartitioning.
- */
- void
- unpack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- &data_range,
- std::vector<VectorType *> &all_out,
- VectorType &valence);
-
-
- /**
- * Registers the pack_callback() function to the
- * parallel::distributed::Triangulation that has been assigned to the
- * DoFHandler class member and stores the returning handle.
- */
- void
- register_data_attach();
- };
+ using SolutionTransfer DEAL_II_DEPRECATED_EARLY =
+ dealii::SolutionTransfer<dim, VectorType, spacedim>;
} // namespace distributed
} // namespace parallel
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;
};
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
DEAL_II_NAMESPACE_OPEN
+
/**
* This class implements the transfer of a discrete FE function (e.g. a
* solution vector) from one mesh to another that is obtained from the first
* by a single refinement and/or coarsening step. During interpolation the
- * vector is reinitialized to the new size and filled with the interpolated
- * values. This class is used in the step-15, step-26, step-31, and step-33
- * tutorial programs. A version of this class that works on parallel
- * triangulations is available as parallel::distributed::SolutionTransfer.
+ * vector is filled with the interpolated values. This class is used in the
+ * step-15, step-26, step-31, and step-33 tutorial programs. This class
+ * works both for serial and distributed meshes.
*
* <h3>Usage</h3>
*
- * This class implements the algorithms in two different ways:
- * <ul>
- * <li> If the grid will only be refined (i.e. no cells are coarsened) then
- * use @p SolutionTransfer as follows:
- * @code
- * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
- *
- * // flag some cells for refinement, e.g.
- * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
- * error_indicators,
- * 0.3,
- * 0);
- * // prepare the triangulation for refinement,
- * tria->prepare_coarsening_and_refinement();
- *
- * // tell the SolutionTransfer object that we intend to do pure refinement,
- * soltrans.prepare_for_pure_refinement();
- *
- * // actually execute the refinement,
- * tria->execute_coarsening_and_refinement();
- *
- * // and redistribute dofs.
- * dof_handler->distribute_dofs (fe);
- * @endcode
- *
- * Then to proceed do
- * @code
- * // take a copy of the solution vector
- * Vector<double> solution_old(solution);
- *
- * // resize solution vector to the correct size, as the @p refine_interpolate
- * // function requires the vectors to be of right sizes
- * solution.reinit(dof_handler->n_dofs());
- *
- * // and finally interpolate
- * soltrans.refine_interpolate(solution_old, solution);
- * @endcode
+ * @note It is important to note, that if you use more than one
+ * SolutionTransfer object at the same time, that the calls to prepare_*()
+ * and interpolate()/deserialize() need to be in the same order.
*
- * Although the @p refine_interpolate functions are allowed to be called
- * multiple times, e.g. for interpolating several solution vectors, there is
- * the following possibility of interpolating several functions simultaneously.
+ * <h3>Transferring a solution</h3> Here VectorType is your favorite
+ * vector type, e.g. Vector, LinearAlgebra::distributed::Vector,
+ * PETScWrappers::MPI::Vector, TrilinosWrappers::MPI::Vector, or corresponding
+ * block vectors.
* @code
- * std::vector<Vector<double> > solutions_old(n_vectors, Vector<double> (n));
- * ...
- * std::vector<Vector<double> > solutions(n_vectors, Vector<double> (n));
- * soltrans.refine_interpolate(solutions_old, solutions);
- * @endcode
- * This is used in several of the tutorial programs, for example step-31
- * and step-33.
- *
- * <li> If the grid has cells that will be coarsened, then use @p
- * SolutionTransfer as follows:
- * @code
- * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
- *
+ * SolutionTransfer<dim, VectorType> soltrans(dof_handler);
* // flag some cells for refinement and coarsening, e.g.
- * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
- * error_indicators,
- * 0.3,
- * 0.05);
+ * GridRefinement::
+ * refine_and_coarsen_fixed_fraction(tria,
+ * error_indicators,
+ * 0.3,
+ * 0.05);
+ * // Use parallel::distributed::GridRefinement for distributed triangulations.
*
* // prepare the triangulation,
- * tria->prepare_coarsening_and_refinement();
+ * tria.prepare_coarsening_and_refinement();
*
- * // prepare the SolutionTransfer object for coarsening and refinement and give
- * // the solution vector that we intend to interpolate later,
+ * // prepare the SolutionTransfer object for coarsening and refinement
+ * // and give the solution vector that we intend to interpolate later,
* soltrans.prepare_for_coarsening_and_refinement(solution);
*
* // actually execute the refinement,
- * tria->execute_coarsening_and_refinement ();
+ * tria.execute_coarsening_and_refinement();
*
* // redistribute dofs,
- * dof_handler->distribute_dofs (fe);
+ * dof_handler.distribute_dofs(fe);
*
* // and interpolate the solution
- * Vector<double> interpolate_solution(dof_handler->n_dofs());
- * soltrans.interpolate(solution, interpolated_solution);
+ * VectorType interpolated_solution;
+ *
+ * //create VectorType in the right size here
+ * soltrans.interpolate(interpolated_solution);
* @endcode
*
- * If the grid is partitioned across several MPI processes, then it is
- * important to note that the old solution(s) must be copied to one that
- * also provides access to the locally relevant DoF values (these values
- * required for the interpolation process):
+ * <h3>Usage on distributed grids</h3>
+ * If the grid is distributed, it is important to note that the old
+ * solution(s) must be copied to one that also provides access to the
+ * locally relevant DoF values (these values required for the interpolation
+ * process):
* @code
* // Create initial indexsets pertaining to the grid before refinement
* const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
* solution.reinit(locally_owned_dofs,
* mpi_communicator);
* ...
- * // Transfer solution to vector that provides access to locally relevant DoFs
+ * // Transfer solution to vector that provides access to
+ * // locally relevant DoFs
* TrilinosWrappers::MPI::Vector old_solution;
* old_solution.reinit(locally_owned_dofs,
* locally_relevant_dofs,
* mpi_communicator);
* old_solution = solution;
+ *
+ * // Initialize SolutionTransfer object
+ * SolutionTransfer<dim, VectorType> soltrans(dof_handler);
+ * soltrans.prepare_for_coarsening_and_refinement(old_solution);
* ...
* // Refine grid
* // Recreate locally_owned_dofs and locally_relevant_dofs index sets
* ...
* solution.reinit(locally_owned_dofs, mpi_communicator);
- * soltrans.refine_interpolate(old_solution, solution);
+ * soltrans.interpolate(solution);
+ * @endcode
+ *
+ * <h4>Note on ghost elements</h4> In a parallel computation PETSc or
+ * Trilinos vector may contain ghost elements or not. For reading in
+ * information with prepare_for_coarsening_and_refinement() or
+ * prepare_for_serialization() you need to supply vectors with ghost
+ * elements, so that all locally_active elements can be read. On the other
+ * hand, ghosted vectors are generally not writable, so for calls to
+ * interpolate() or deserialize() you need to supply distributed vectors
+ * without ghost elements. More precisely, during interpolation the
+ * current algorithm writes into all locally active degrees of freedom.
+ *
+ * Different from PETSc and Trilinos vectors, LinearAlgebra::distributed::Vector
+ * allows writing into ghost elements. For a ghosted vector the interpolation
+ * step can be accomplished via
+ * @code
+ * interpolated_solution.zero_out_ghost_values();
+ * soltrans.interpolate(interpolated_solution);
+ * interpolated_solution.update_ghost_values();
* @endcode
*
- * Multiple calls to the function <code>interpolate (const VectorType &in,
- * VectorType &out)</code> are NOT allowed. Interpolating several
- * functions can be performed in one step by using <tt>void interpolate (const
- * vector<VectorType> &all_in, vector<VectorType> &all_out)
- * const</tt>, and using the respective @p
- * prepare_for_coarsening_and_refinement function taking several vectors as
- * input before actually refining and coarsening the triangulation (see
- * there).
- * </ul>
- *
- * For deleting all stored data in @p SolutionTransfer and reinitializing it
- * use the <tt>clear()</tt> function.
- *
- * The template argument @p VectorType denotes the type of data container you
- * want to transfer.
- *
- *
- * <h3>Interpolating in the presence of hanging nodes and boundary values</h3>
- *
- * The interpolation onto the new mesh is a local operation, i.e., it
- * interpolates onto the new mesh only. If that new mesh has hanging nodes,
- * you will therefore get a solution that does not satisfy hanging node
- * constraints. The same is true with boundary values: the interpolated
- * solution will just be the interpolation of the old solution at the
- * boundary, and this may or may not satisfy boundary values at newly
- * introduced boundary nodes.
- *
- * Consequently, you may have to apply hanging node or boundary value
- * constraints after interpolation. step-15 and step-26 have examples of
- * dealing with this.
- *
- *
- * <h3>Implementation</h3>
- *
- * <ul>
- * <li> Solution transfer with only refinement. Assume that we have got a
- * solution vector on the current (original) grid. Each entry of this vector
- * belongs to one of the DoFs of the discretization. If we now refine the grid
- * then the calling of DoFHandler::distribute_dofs() will change at least some
- * of the DoF indices. Hence we need to store the DoF indices of all active
- * cells before the refinement. A pointer for each active cell is used to
- * point to the vector of these DoF indices of that cell. This is done by
- * prepare_for_pure_refinement().
- *
- * In the function <tt>refine_interpolate(in,out)</tt> and on each cell where
- * the pointer is set (i.e. the cells that were active in the original grid)
- * we can now access the local values of the solution vector @p in on that
- * cell by using the stored DoF indices. These local values are interpolated
- * and set into the vector @p out that is at the end the discrete function @p
- * in interpolated on the refined mesh.
- *
- * The <tt>refine_interpolate(in,out)</tt> function can be called multiple
- * times for arbitrary many discrete functions (solution vectors) on the
- * original grid.
- *
- * <li> Solution transfer with coarsening and refinement. After calling
- * Triangulation::prepare_coarsening_and_refinement the coarsen flags of
- * either all or none of the children of a (father-)cell are set. While
- * coarsening (Triangulation::execute_coarsening_and_refinement) the cells
- * that are not needed any more will be deleted from the Triangulation.
- *
- * For the interpolation from the (to be coarsenend) children to their father
- * the children cells are needed. Hence this interpolation and the storing of
- * the interpolated values of each of the discrete functions that we want to
- * interpolate needs to take place before these children cells are coarsened
- * (and deleted!!). Again a pointer for each relevant cell is set to point to
- * these values (see below). Additionally the DoF indices of the cells that
- * will not be coarsened need to be stored according to the solution transfer
- * with pure refinement (cf there). All this is performed by
- * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
- * <tt>vector<VectorType> all_in</tt> includes all discrete
- * functions to be interpolated onto the new grid.
- *
- * As we need two different kinds of pointers (<tt>vector<unsigned int> *</tt>
- * for the Dof indices and <tt>vector<VectorType> *</tt> for the
- * interpolated DoF values) we use the @p Pointerstruct that includes both of
- * these pointers and the pointer for each cell points to these @p
- * Pointerstructs. On each cell only one of the two different pointers is used
- * at one time hence we could use a <tt>void * pointer</tt> as
- * <tt>vector<unsigned int> *</tt> at one time and as
- * <tt>vector<VectorType> *</tt> at the other but using this @p
- * Pointerstruct in between makes the use of these pointers more safe and
- * gives better possibility to expand their usage.
- *
- * In <tt>interpolate(all_in, all_out)</tt> the refined cells are treated
- * according to the solution transfer while pure refinement. Additionally, on
- * each cell that is coarsened (hence previously was a father cell), the
- * values of the discrete functions in @p all_out are set to the stored local
- * interpolated values that are accessible due to the 'vector<VectorType>
- * *' pointer in @p Pointerstruct that is pointed to by the pointer of that
- * cell. It is clear that <tt>interpolate(all_in, all_out)</tt> only can be
- * called with the <tt>vector<VectorType> all_in</tt> that previously was
- * the parameter of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt>
- * function. Hence <tt>interpolate(all_in, all_out)</tt> can (in contrast to
- * <tt>refine_interpolate(in, out)</tt>) only be called once.
- * </ul>
+ * <h3>Use for serialization</h3>
+ *
+ * This class can be used to serialize and later deserialize a
+ * mesh with solution vectors to a file. If you use more than one
+ * DoFHandler and therefore more than one SolutionTransfer object, they
+ * need to be serialized and deserialized in the same order.
+ *
+ * If vector has the locally relevant DoFs, serialization works as
+ * follows:
+ * @code
+ * SolutionTransfer<dim, VectorType> sol_trans(dof_handler);
+ * sol_trans.prepare_for_serialization(vector);
+ *
+ * triangulation.save(filename);
+ * @endcode
+ * For deserialization the vector needs to be a distributed vector
+ * (without ghost elements):
+ * @code
+ * //[create coarse mesh...]
+ * triangulation.load(filename);
+ *
+ * SolutionTransfer<dim, VectorType> sol_trans(dof_handler);
+ * sol_trans.deserialize(distributed_vector);
+ * @endcode
+ *
+ *
+ * <h3>Note on usage with DoFHandler with hp-capabilities</h3>
+ *
+ * Since data on DoFHandler objects with hp-capabilities is associated with
+ * many different FiniteElement objects, each cell's data has to be
+ * processed with its corresponding `future_fe_index`. Further, if
+ * refinement is involved, data will be packed on the parent cell with its
+ * `future_fe_index` and unpacked later with the same index on its children.
+ * If cells get coarsened into one, data will be packed on the children with
+ * the least dominant finite element of their common subspace, and unpacked
+ * on the parent with this particular finite element (consult
+ * hp::FECollection::find_dominated_fe_extended() for more information).
+ *
+ * Transferring a solution across refinement works exactly like in the
+ * non-hp-case. However, when considering serialization, we also have to
+ * store the active FE indices in an additional step. A code snippet
+ * demonstrating serialization with the SolutionTransfer class with DoFHandler
+ * objects with hp-capabilities is provided in the following. Here VectorType is
+ * your favorite vector type, e.g. Vector, LinearAlgebra::distributed::Vector,
+ * PETScWrappers::MPI::Vector, TrilinosWrappers::MPI::Vector, or corresponding
+ * block vectors.
+ *
+ * If vector has the locally relevant DoFs, serialization works as follows:
+ * @code
+ * SolutionTransfer<dim, VectorType, DoFHandler<dim,spacedim>>
+ * sol_trans(hp_dof_handler);
+ *
+ * hp_dof_handler.prepare_for_serialization_of_active_fe_indices();
+ * sol_trans.prepare_for_serialization(vector);
+ *
+ * triangulation.save(filename);
+ * @endcode
+ *
+ * For deserialization the vector needs to be a distributed vector
+ * (without ghost elements):
+ * @code
+ * //[create coarse mesh...]
+ * triangulation.load(filename);
+ *
+ * hp::FECollection<dim,spacedim> fe_collection;
+ * //[prepare identical fe_collection...]
+ *
+ * DoFHandler<dim,spacedim> hp_dof_handler(triangulation);
+ * // We need to introduce our dof_handler to the fe_collection
+ * // before setting all active FE indices.
+ * hp_dof_handler.deserialize_active_fe_indices();
+ * hp_dof_handler.distribute_dofs(fe_collection);
+ *
+ * SolutionTransfer<dim,VectorType,DoFHandler<dim,spacedim>>
+ * sol_trans(hp_dof_handler);
+ * sol_trans.deserialize(distributed_vector);
+ * @endcode
*
*
* <h3>Interaction with hanging nodes</h3>
*
* <h3>Implementation in the context of hp-finite elements</h3>
*
- * In the case of DoFHandlers with hp-capabilities, nothing defines which of the
- * finite elements that are part of the hp::FECollection associated with the
- * DoFHandler, should be considered on cells that are not active (i.e., that
- * have children). This is because degrees of freedom are only allocated for
- * active cells and, in fact, it is not allowed to set an active FE index on
- * non-active cells using DoFAccessor::set_active_fe_index().
+ * In the case of DoFHandlers with hp-capabilities, nothing defines which of
+ * the finite elements that are part of the hp::FECollection associated with
+ * the DoFHandler, should be considered on cells that are not active (i.e.,
+ * that have children). This is because degrees of freedom are only allocated
+ * for active cells and, in fact, it is not allowed to set an active FE index
+ * on non-active cells using DoFAccessor::set_active_fe_index().
*
* It is, thus, not entirely natural what should happen if, for example, a few
* cells are coarsened away. This class then implements the following
* algorithm:
* - If a cell is refined, then the values of the solution vector(s) are
- * interpolated before refinement on the to-be-refined cell from the space of
- * the active finite element to the one of the future finite element. These
+ * interpolated before refinement on the to-be-refined cell from the space
+ * of the active finite element to the one of the future finite element. These
* values are then distributed on the finite element spaces of the children
* post-refinement. This may lose information if, for example, the old cell
* used a Q2 space and the children use Q1 spaces, or the information may be
* prolonged if the mother cell used a Q1 space and the children are Q2s.
* - If cells are to be coarsened, then the values from the child cells are
- * interpolated to the mother cell using the largest of the child cell future
- * finite element spaces, which will be identified as the least dominant
- * element following the FiniteElementDomination logic (consult
+ * interpolated to the mother cell using the largest of the child cell
+ * future finite element spaces, which will be identified as the least
+ * dominant element following the FiniteElementDomination logic (consult
* hp::FECollection::find_dominated_fe_extended() for more information). For
* example, if the children of a cell use Q1, Q2 and Q3 spaces, then the
* values from the children are interpolated into a Q3 space on the mother
* cell. After refinement, this Q3 function on the mother cell is then
- * interpolated into the space the user has selected for this cell (which may
- * be different from Q3, in this example, if the user has set the
- * active FE index for a different space post-refinement and before calling
+ * interpolated into the space the user has selected for this cell (which
+ * may be different from Q3, in this example, if the user has set the active
+ * FE index for a different space post-refinement and before calling
* DoFHandler::distribute_dofs()).
*
* @note In the context of hp-refinement, if cells are coarsened or the
{
public:
/**
- * Constructor, takes the current DoFHandler as argument.
+ * Constructor.
+ *
+ * @param[in] dof_handler The DoFHandler on which all operations will
+ * happen. At the time when this constructor is called, the DoFHandler
+ * still points to the Triangulation before the refinement in question
+ * happens.
+ * @param[in] average_values Average the contribututions to the same
+ * DoF coming from different cells. Note: averaging requires an
+ * additional communication step, since the valence of the DoF has to be
+ * determined.
*/
- SolutionTransfer(const DoFHandler<dim, spacedim> &dof);
+ SolutionTransfer(const DoFHandler<dim, spacedim> &dof_handler,
+ const bool average_values = false);
/**
- * Destructor
+ * Destructor.
*/
- ~SolutionTransfer();
+ ~SolutionTransfer() = default;
/**
- * Reinit this class to the state that it has directly after calling the
- * Constructor
+ * Prepare the current object for coarsening and refinement. It
+ * stores pointers to the vectors in @p all_in to be interpolated onto the new
+ * (refined and/or coarsened) grid, and registers this object for data
+ * transfer on the grid.
*/
void
- clear();
+ prepare_for_coarsening_and_refinement(
+ const std::vector<const VectorType *> &all_in);
/**
- * Prepares the @p SolutionTransfer for pure refinement. It stores the dof
- * indices of each cell. After calling this function only calling the @p
- * refine_interpolate functions is allowed.
- */
- void
- prepare_for_pure_refinement();
-
- /**
- * Prepares the @p SolutionTransfer for coarsening and refinement. It stores
- * the dof indices of each cell and stores the dof values of the vectors in
- * @p all_in in each cell that'll be coarsened. @p all_in includes all
- * vectors that are to be interpolated onto the new (refined and/or
- * coarsenend) grid.
+ * Same as above but without pointers.
*/
void
prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in);
/**
- * Same as previous function but for only one discrete function to be
+ * Same as the previous function but for only one discrete function to be
* interpolated.
*/
void
prepare_for_coarsening_and_refinement(const VectorType &in);
/**
- * This function interpolates the discrete function @p in, which is a vector
- * on the grid before the refinement, to the function @p out which then is a
- * vector on the refined grid. It assumes the vectors having the right sizes
- * (i.e. <tt>in.size()==n_dofs_old</tt>,
- * <tt>out.size()==n_dofs_refined</tt>)
- *
- * Calling this function is allowed only if @p prepare_for_pure_refinement
- * is called and the refinement is executed before. Multiple calling of this
- * function is allowed. e.g. for interpolating several functions.
+ * Interpolate the data previously stored in this object before the mesh
+ * was refined or coarsened onto the current set of cells. Do so for
+ * each of the vectors provided to
+ * prepare_for_coarsening_and_refinement() and write the result into the
+ * given set of vectors.
*/
void
- refine_interpolate(const VectorType &in, VectorType &out) const;
+ interpolate(std::vector<VectorType *> &all_out);
/**
- * This function interpolates the discrete functions that are stored in @p
- * all_in onto the refined and/or coarsenend grid. It assumes the vectors in
- * @p all_in denote the same vectors as in @p all_in as parameter of
- * <tt>prepare_for_refinement_and_coarsening(all_in)</tt>. However, there is
- * no way of verifying this internally, so be careful here.
- *
- * Calling this function is allowed only if first
- * Triangulation::prepare_coarsening_and_refinement, second @p
- * SolutionTransfer::prepare_for_coarsening_and_refinement, an then third
- * Triangulation::execute_coarsening_and_refinement are called before.
- * Multiple calling of this function is NOT allowed. Interpolating several
- * functions can be performed in one step.
- *
- * The number of output vectors is assumed to be the same as the number of
- * input vectors. Also, the sizes of the output vectors are assumed to be of
- * the right size (@p n_dofs_refined). Otherwise an assertion will be
- * thrown.
+ * Same as above but without pointers.
*/
void
- interpolate(const std::vector<VectorType> &all_in,
- std::vector<VectorType> &all_out) const;
+ interpolate(std::vector<VectorType> &all_out);
/**
* Same as the previous function. It interpolates only one function. It
* assumes the vectors having the right sizes (i.e.
* <tt>in.size()==n_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
*
- * Multiple calling of this function is NOT allowed. Interpolating several
- * functions can be performed in one step by using <tt>interpolate (all_in,
- * all_out)</tt>
+ * Multiple calling of this function is NOT allowed. Interpolating
+ * several functions can be performed in one step by using
+ * <tt>interpolate (all_in, all_out)</tt>
+ */
+ void
+ interpolate(VectorType &out);
+
+ /**
+ * Prepare the serialization of the given vector. The serialization is
+ * done by Triangulation::save(). The given vector needs all information
+ * on the locally active DoFs (it must be ghosted). See documentation of
+ * this class for more information.
*/
void
- interpolate(const VectorType &in, VectorType &out) const;
+ prepare_for_serialization(const VectorType &in);
/**
- * Determine an estimate for the memory consumption (in bytes) of this
- * object.
+ * Same as the function above, only for a list of vectors.
*/
- std::size_t
- memory_consumption() const;
+ void
+ prepare_for_serialization(const std::vector<const VectorType *> &all_in);
/**
- * Exception
+ * Execute the deserialization of the given vector. This needs to be
+ * done after calling Triangulation::load(). The given vector must be a
+ * fully distributed vector without ghost elements. See documentation of
+ * this class for more information.
*/
- DeclExceptionMsg(ExcNotPrepared,
- "You are attempting an operation for which this object is "
- "not prepared. This may be because you either did not call "
- "one of the prepare_*() functions at all, or because you "
- "called the wrong one for the operation you are currently "
- "attempting.");
+ void
+ deserialize(VectorType &in);
+
/**
- * Exception
+ * Same as the function above, only for a list of vectors.
*/
- DeclExceptionMsg(
- ExcAlreadyPrepForRef,
- "You are attempting to call one of the prepare_*() functions "
- "of this object to prepare it for an operation for which it "
- "is already prepared. Specifically, the object was "
- "previously prepared for pure refinement.");
+ void
+ deserialize(std::vector<VectorType *> &all_in);
/**
- * Exception
+ * Reinit this class to the state that it has directly after calling the
+ * constructor.
*/
- DeclExceptionMsg(
- ExcAlreadyPrepForCoarseAndRef,
- "You are attempting to call one of the prepare_*() functions "
- "of this object to prepare it for an operation for which it "
- "is already prepared. Specifically, the object was "
- "previously prepared for both coarsening and refinement.");
+ DEAL_II_DEPRECATED_EARLY void
+ clear();
private:
/**
dof_handler;
/**
- * Stores the number of DoFs before the refinement and/or coarsening.
+ * Flag indicating if averaging should be performed.
*/
- types::global_dof_index n_dofs_old;
+ const bool average_values;
/**
- * Declaration of @p PreparationState that denotes the three possible states
- * of the @p SolutionTransfer: being prepared for 'pure refinement',
- * prepared for 'coarsening and refinement' or not prepared.
+ * A vector that stores pointers to all the vectors we are supposed to
+ * copy over from the old to the new mesh.
*/
- enum PreparationState
- {
- /**
- * The SolutionTransfer is not yet prepared.
- */
- none,
- /**
- * The SolutionTransfer is prepared for purely refinement.
- */
- pure_refinement,
- /**
- * The SolutionTransfer is prepared for coarsening and refinement.
- */
- coarsening_and_refinement
- };
+ std::vector<const VectorType *> input_vectors;
/**
- * Definition of the respective variable.
+ * The handle that the Triangulation has assigned to this object
+ * with which we can access our memory offset and our pack function.
*/
- PreparationState prepared_for;
-
+ unsigned int handle;
/**
- * Is used for @p prepare_for_refining (of course also for @p
- * repare_for_refining_and_coarsening) and stores all dof indices of the
- * cells that'll be refined
+ * A callback function used to pack the data on the current mesh into
+ * objects that can later be retrieved after refinement, coarsening and
+ * repartitioning.
*/
- std::vector<std::vector<types::global_dof_index>> indices_on_cell;
+ std::vector<char>
+ pack_callback(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const CellStatus status);
/**
- * All cell data (the dof indices and the dof values) should be accessible
- * from each cell. As each cell has got only one @p user_pointer, multiple
- * pointers to the data need to be packetized in a structure. Note that in
- * our case on each cell either the <tt>vector<unsigned int> indices</tt>
- * (if the cell will be refined) or the <tt>vector<double> dof_values</tt>
- * (if the children of this cell will be deleted) is needed, hence one @p
- * user_pointer should be sufficient, but to allow some error checks and to
- * preserve the user from making user errors the @p user_pointer will be
- * 'multiplied' by this structure.
+ * A callback function used to unpack the data on the current mesh that
+ * has been packed up previously on the mesh before refinement,
+ * coarsening and repartitioning.
*/
- struct Pointerstruct
- {
- Pointerstruct()
- : indices_ptr(nullptr)
- , dof_values_ptr(nullptr)
- , active_fe_index(0)
- {}
- Pointerstruct(std::vector<types::global_dof_index> *indices_ptr_in,
- const unsigned int active_fe_index_in = 0)
- : indices_ptr(indices_ptr_in)
- , dof_values_ptr(nullptr)
- , active_fe_index(active_fe_index_in)
- {}
- Pointerstruct(
- std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr_in,
- const unsigned int active_fe_index_in = 0)
- : indices_ptr(nullptr)
- , dof_values_ptr(dof_values_ptr_in)
- , active_fe_index(active_fe_index_in)
- {}
- std::size_t
- memory_consumption() const;
+ void
+ unpack_callback(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const CellStatus status,
+ const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+ std::vector<VectorType *> &all_out,
+ VectorType &valence);
- std::vector<types::global_dof_index> *indices_ptr;
- std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr;
- unsigned int active_fe_index;
- };
/**
- * Map mapping from level and index of cell to the @p Pointerstructs (cf.
- * there). This map makes it possible to keep all the information needed to
- * transfer the solution inside this object rather than using user pointers
- * of the Triangulation for this purpose.
+ * Registers the pack_callback() function to the Triangulation that has been
+ * assigned to the DoFHandler class member and stores the returning handle.
*/
- std::map<std::pair<unsigned int, unsigned int>, Pointerstruct> cell_map;
+ void
+ register_data_attach();
+};
+
+
+namespace Legacy
+{
/**
- * Is used for @p prepare_for_refining_and_coarsening The interpolated dof
- * values of all cells that'll be coarsened will be stored in this vector.
+ * This class implements the transfer of a discrete FE function (e.g. a
+ * solution vector) from one mesh to another that is obtained from the first
+ * by a single refinement and/or coarsening step. During interpolation the
+ * vector is reinitialized to the new size and filled with the interpolated
+ * values. This class is used in the step-15, step-26, step-31, and step-33
+ * tutorial programs. A version of this class that works on all types of
+ * triangulations, including distributed ones, is available as
+ * dealii::SolutionTransfer.
+ *
+ * <h3>Usage</h3>
+ *
+ * This class implements the algorithms in two different ways:
+ * <ul>
+ * <li> If the grid will only be refined (i.e. no cells are coarsened) then
+ * use @p SolutionTransfer as follows:
+ * @code
+ * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
+ *
+ * // flag some cells for refinement, e.g.
+ * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
+ * error_indicators,
+ * 0.3,
+ * 0);
+ * // prepare the triangulation for refinement,
+ * tria->prepare_coarsening_and_refinement();
+ *
+ * // tell the SolutionTransfer object that we intend to do pure refinement,
+ * soltrans.prepare_for_pure_refinement();
+ *
+ * // actually execute the refinement,
+ * tria->execute_coarsening_and_refinement();
+ *
+ * // and redistribute dofs.
+ * dof_handler->distribute_dofs (fe);
+ * @endcode
+ *
+ * Then to proceed do
+ * @code
+ * // take a copy of the solution vector
+ * Vector<double> solution_old(solution);
+ *
+ * // resize solution vector to the correct size, as the @p refine_interpolate
+ * // function requires the vectors to be of right sizes
+ * solution.reinit(dof_handler->n_dofs());
+ *
+ * // and finally interpolate
+ * soltrans.refine_interpolate(solution_old, solution);
+ * @endcode
+ *
+ * Although the @p refine_interpolate functions are allowed to be called
+ * multiple times, e.g. for interpolating several solution vectors, there is
+ * the following possibility of interpolating several functions
+ * simultaneously.
+ * @code
+ * std::vector<Vector<double> > solutions_old(n_vectors, Vector<double> (n));
+ * ...
+ * std::vector<Vector<double> > solutions(n_vectors, Vector<double> (n));
+ * soltrans.refine_interpolate(solutions_old, solutions);
+ * @endcode
+ * This is used in several of the tutorial programs, for example step-31
+ * and step-33.
+ *
+ * <li> If the grid has cells that will be coarsened, then use @p
+ * SolutionTransfer as follows:
+ * @code
+ * SolutionTransfer<dim, Vector<double> > soltrans(*dof_handler);
+ *
+ * // flag some cells for refinement and coarsening, e.g.
+ * GridRefinement::refine_and_coarsen_fixed_fraction(*tria,
+ * error_indicators,
+ * 0.3,
+ * 0.05);
+ *
+ * // prepare the triangulation,
+ * tria->prepare_coarsening_and_refinement();
+ *
+ * // prepare the SolutionTransfer object for coarsening and refinement and
+ * give
+ * // the solution vector that we intend to interpolate later,
+ * soltrans.prepare_for_coarsening_and_refinement(solution);
+ *
+ * // actually execute the refinement,
+ * tria->execute_coarsening_and_refinement ();
+ *
+ * // redistribute dofs,
+ * dof_handler->distribute_dofs (fe);
+ *
+ * // and interpolate the solution
+ * Vector<double> interpolate_solution(dof_handler->n_dofs());
+ * soltrans.interpolate(solution, interpolated_solution);
+ * @endcode
+ *
+ * If the grid is partitioned across several MPI processes, then it is
+ * important to note that the old solution(s) must be copied to one that
+ * also provides access to the locally relevant DoF values (these values
+ * required for the interpolation process):
+ * @code
+ * // Create initial indexsets pertaining to the grid before refinement
+ * const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
+ * const IndexSet locally_relevant_dofs =
+ * DoFTools::extract_locally_relevant_dofs(dof_handler);
+ *
+ * // The solution vector only knows about locally owned DoFs
+ * TrilinosWrappers::MPI::Vector solution;
+ * solution.reinit(locally_owned_dofs,
+ * mpi_communicator);
+ * ...
+ * // Transfer solution to vector that provides access to locally relevant
+ * DoFs TrilinosWrappers::MPI::Vector old_solution;
+ * old_solution.reinit(locally_owned_dofs,
+ * locally_relevant_dofs,
+ * mpi_communicator);
+ * old_solution = solution;
+ * ...
+ * // Refine grid
+ * // Recreate locally_owned_dofs and locally_relevant_dofs index sets
+ * ...
+ * solution.reinit(locally_owned_dofs, mpi_communicator);
+ * soltrans.refine_interpolate(old_solution, solution);
+ * @endcode
+ *
+ * Multiple calls to the function <code>interpolate (const VectorType &in,
+ * VectorType &out)</code> are NOT allowed. Interpolating several
+ * functions can be performed in one step by using <tt>void interpolate (const
+ * vector<VectorType> &all_in, vector<VectorType> &all_out)
+ * const</tt>, and using the respective @p
+ * prepare_for_coarsening_and_refinement function taking several vectors as
+ * input before actually refining and coarsening the triangulation (see
+ * there).
+ * </ul>
+ *
+ * For deleting all stored data in @p SolutionTransfer and reinitializing it
+ * use the <tt>clear()</tt> function.
+ *
+ * The template argument @p VectorType denotes the type of data container you
+ * want to transfer.
+ *
+ *
+ * <h3>Interpolating in the presence of hanging nodes and boundary values</h3>
+ *
+ * The interpolation onto the new mesh is a local operation, i.e., it
+ * interpolates onto the new mesh only. If that new mesh has hanging nodes,
+ * you will therefore get a solution that does not satisfy hanging node
+ * constraints. The same is true with boundary values: the interpolated
+ * solution will just be the interpolation of the old solution at the
+ * boundary, and this may or may not satisfy boundary values at newly
+ * introduced boundary nodes.
+ *
+ * Consequently, you may have to apply hanging node or boundary value
+ * constraints after interpolation. step-15 and step-26 have examples of
+ * dealing with this.
+ *
+ *
+ * <h3>Implementation</h3>
+ *
+ * <ul>
+ * <li> Solution transfer with only refinement. Assume that we have got a
+ * solution vector on the current (original) grid. Each entry of this vector
+ * belongs to one of the DoFs of the discretization. If we now refine the grid
+ * then the calling of DoFHandler::distribute_dofs() will change at least some
+ * of the DoF indices. Hence we need to store the DoF indices of all active
+ * cells before the refinement. A pointer for each active cell is used to
+ * point to the vector of these DoF indices of that cell. This is done by
+ * prepare_for_pure_refinement().
+ *
+ * In the function <tt>refine_interpolate(in,out)</tt> and on each cell where
+ * the pointer is set (i.e. the cells that were active in the original grid)
+ * we can now access the local values of the solution vector @p in on that
+ * cell by using the stored DoF indices. These local values are interpolated
+ * and set into the vector @p out that is at the end the discrete function @p
+ * in interpolated on the refined mesh.
+ *
+ * The <tt>refine_interpolate(in,out)</tt> function can be called multiple
+ * times for arbitrary many discrete functions (solution vectors) on the
+ * original grid.
+ *
+ * <li> Solution transfer with coarsening and refinement. After calling
+ * Triangulation::prepare_coarsening_and_refinement the coarsen flags of
+ * either all or none of the children of a (father-)cell are set. While
+ * coarsening (Triangulation::execute_coarsening_and_refinement) the cells
+ * that are not needed any more will be deleted from the Triangulation.
+ *
+ * For the interpolation from the (to be coarsenend) children to their father
+ * the children cells are needed. Hence this interpolation and the storing of
+ * the interpolated values of each of the discrete functions that we want to
+ * interpolate needs to take place before these children cells are coarsened
+ * (and deleted!!). Again a pointer for each relevant cell is set to point to
+ * these values (see below). Additionally the DoF indices of the cells that
+ * will not be coarsened need to be stored according to the solution transfer
+ * with pure refinement (cf there). All this is performed by
+ * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
+ * <tt>vector<VectorType> all_in</tt> includes all discrete
+ * functions to be interpolated onto the new grid.
+ *
+ * As we need two different kinds of pointers (<tt>vector<unsigned int> *</tt>
+ * for the Dof indices and <tt>vector<VectorType> *</tt> for the
+ * interpolated DoF values) we use the @p Pointerstruct that includes both of
+ * these pointers and the pointer for each cell points to these @p
+ * Pointerstructs. On each cell only one of the two different pointers is used
+ * at one time hence we could use a <tt>void * pointer</tt> as
+ * <tt>vector<unsigned int> *</tt> at one time and as
+ * <tt>vector<VectorType> *</tt> at the other but using this @p
+ * Pointerstruct in between makes the use of these pointers more safe and
+ * gives better possibility to expand their usage.
+ *
+ * In <tt>interpolate(all_in, all_out)</tt> the refined cells are treated
+ * according to the solution transfer while pure refinement. Additionally, on
+ * each cell that is coarsened (hence previously was a father cell), the
+ * values of the discrete functions in @p all_out are set to the stored local
+ * interpolated values that are accessible due to the 'vector<VectorType>
+ * *' pointer in @p Pointerstruct that is pointed to by the pointer of that
+ * cell. It is clear that <tt>interpolate(all_in, all_out)</tt> only can be
+ * called with the <tt>vector<VectorType> all_in</tt> that previously was
+ * the parameter of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt>
+ * function. Hence <tt>interpolate(all_in, all_out)</tt> can (in contrast to
+ * <tt>refine_interpolate(in, out)</tt>) only be called once.
+ * </ul>
+ *
+ *
+ * <h3>Interaction with hanging nodes</h3>
+ *
+ * This class does its best to represent on the new mesh the finite element
+ * function that existed on the old mesh, but this may lead to situations
+ * where the function on the new mesh is no longer conforming at hanging
+ * nodes. To this end, consider a situation of a twice refined mesh that
+ * started with a single square cell (i.e., we now have 16 cells). Consider
+ * also that we coarsen 4 of the cells back to the first refinement level. In
+ * this case, we end up with a mesh that will look as follows if we were to
+ * use a $Q_1$ element:
+ *
+ * @image html hanging_nodes.png ""
+ *
+ * The process of interpolating from the old to the new mesh would imply that
+ * the values of the finite element function will not change on all of the
+ * cells that remained as they are (i.e., the fine cells) but that on the
+ * coarse cell at the top right, the four values at the vertices are obtained
+ * by interpolating down from its former children. If the original function
+ * was not linear, this implies that the marked hanging nodes will retain
+ * their old values which, in general, will not lead to a continuous function
+ * along the corresponding edges. In other words, the solution vector obtained
+ * after SolutionTransfer::interpolate() does not satisfy hanging node
+ * constraints: it corresponds to the pointwise interpolation, but not to the
+ * interpolation <i>onto the new finite element space that contains
+ * constraints from hanging nodes</i>.
+ *
+ * Whether this is a problem you need to worry about or not depends on your
+ * application. The situation is easily corrected, of course, by applying
+ * AffineConstraints::distribute() to your solution vector after transfer,
+ * using a constraints object computed on the new DoFHandler object (you
+ * probably need to create this object anyway if you have hanging nodes). This
+ * is also what is done, for example, in step-15.
+ *
+ * @note This situation can only happen if you do coarsening. If all cells
+ * remain as they are or are refined, then SolutionTransfer::interpolate()
+ * computes a new vector of nodal values, but the function represented is of
+ * course exactly the same because the old finite element space is a subspace
+ * of the new one. Thus, if the old function was conforming (i.e., satisfied
+ * hanging node constraints), then so does the new one, and it is not
+ * necessary to call AffineConstraints::distribute().
+ *
+ *
+ * <h3>Implementation in the context of hp-finite elements</h3>
+ *
+ * In the case of DoFHandlers with hp-capabilities, nothing defines which of
+ * the finite elements that are part of the hp::FECollection associated with
+ * the DoFHandler, should be considered on cells that are not active (i.e.,
+ * that have children). This is because degrees of freedom are only allocated
+ * for active cells and, in fact, it is not allowed to set an active FE index
+ * on non-active cells using DoFAccessor::set_active_fe_index().
+ *
+ * It is, thus, not entirely natural what should happen if, for example, a few
+ * cells are coarsened away. This class then implements the following
+ * algorithm:
+ * - If a cell is refined, then the values of the solution vector(s) are
+ * interpolated before refinement on the to-be-refined cell from the space
+ * of the active finite element to the one of the future finite element. These
+ * values are then distributed on the finite element spaces of the children
+ * post-refinement. This may lose information if, for example, the old cell
+ * used a Q2 space and the children use Q1 spaces, or the information may be
+ * prolonged if the mother cell used a Q1 space and the children are Q2s.
+ * - If cells are to be coarsened, then the values from the child cells are
+ * interpolated to the mother cell using the largest of the child cell
+ * future finite element spaces, which will be identified as the least
+ * dominant element following the FiniteElementDomination logic (consult
+ * hp::FECollection::find_dominated_fe_extended() for more information). For
+ * example, if the children of a cell use Q1, Q2 and Q3 spaces, then the
+ * values from the children are interpolated into a Q3 space on the mother
+ * cell. After refinement, this Q3 function on the mother cell is then
+ * interpolated into the space the user has selected for this cell (which
+ * may be different from Q3, in this example, if the user has set the active
+ * FE index for a different space post-refinement and before calling
+ * DoFHandler::distribute_dofs()).
+ *
+ * @note In the context of hp-refinement, if cells are coarsened or the
+ * polynomial degree is lowered on some cells, then the old finite element
+ * space is not a subspace of the new space and you may run into the same
+ * situation as discussed above with hanging nodes. You may want to consider
+ * calling AffineConstraints::distribute() on the vector obtained by
+ * transferring the solution.
+ *
+ * @ingroup numerics
+ *
+ * @deprecated Use dealii::SolutionTransfer instead.
*/
- std::vector<std::vector<Vector<typename VectorType::value_type>>>
- dof_values_on_cell;
-};
+ template <int dim, typename VectorType = Vector<double>, int spacedim = dim>
+ class SolutionTransfer
+ {
+ public:
+ /**
+ * Constructor, takes the current DoFHandler as argument.
+ */
+ SolutionTransfer(const DoFHandler<dim, spacedim> &dof);
+
+ /**
+ * Destructor
+ */
+ ~SolutionTransfer();
+
+ /**
+ * Reinit this class to the state that it has directly after calling the
+ * Constructor
+ */
+ void
+ clear();
+
+ /**
+ * Prepares the @p SolutionTransfer for pure refinement. It stores the dof
+ * indices of each cell. After calling this function only calling the @p
+ * refine_interpolate functions is allowed.
+ */
+ void
+ prepare_for_pure_refinement();
+
+ /**
+ * Prepares the @p SolutionTransfer for coarsening and refinement. It stores
+ * the dof indices of each cell and stores the dof values of the vectors in
+ * @p all_in in each cell that'll be coarsened. @p all_in includes all
+ * vectors that are to be interpolated onto the new (refined and/or
+ * coarsenend) grid.
+ */
+ void
+ prepare_for_coarsening_and_refinement(
+ const std::vector<VectorType> &all_in);
+ /**
+ * Same as previous function but for only one discrete function to be
+ * interpolated.
+ */
+ void
+ prepare_for_coarsening_and_refinement(const VectorType &in);
+
+ /**
+ * This function interpolates the discrete function @p in, which is a vector
+ * on the grid before the refinement, to the function @p out which then is a
+ * vector on the refined grid. It assumes the vectors having the right sizes
+ * (i.e. <tt>in.size()==n_dofs_old</tt>,
+ * <tt>out.size()==n_dofs_refined</tt>)
+ *
+ * Calling this function is allowed only if @p prepare_for_pure_refinement
+ * is called and the refinement is executed before. Multiple calling of this
+ * function is allowed. e.g. for interpolating several functions.
+ */
+ void
+ refine_interpolate(const VectorType &in, VectorType &out) const;
+
+ /**
+ * This function interpolates the discrete functions that are stored in @p
+ * all_in onto the refined and/or coarsenend grid. It assumes the vectors in
+ * @p all_in denote the same vectors as in @p all_in as parameter of
+ * <tt>prepare_for_refinement_and_coarsening(all_in)</tt>. However, there is
+ * no way of verifying this internally, so be careful here.
+ *
+ * Calling this function is allowed only if first
+ * Triangulation::prepare_coarsening_and_refinement, second @p
+ * SolutionTransfer::prepare_for_coarsening_and_refinement, an then third
+ * Triangulation::execute_coarsening_and_refinement are called before.
+ * Multiple calling of this function is NOT allowed. Interpolating several
+ * functions can be performed in one step.
+ *
+ * The number of output vectors is assumed to be the same as the number of
+ * input vectors. Also, the sizes of the output vectors are assumed to be of
+ * the right size (@p n_dofs_refined). Otherwise an assertion will be
+ * thrown.
+ */
+ void
+ interpolate(const std::vector<VectorType> &all_in,
+ std::vector<VectorType> &all_out) const;
+
+ /**
+ * Same as the previous function. It interpolates only one function. It
+ * assumes the vectors having the right sizes (i.e.
+ * <tt>in.size()==n_dofs_old</tt>, <tt>out.size()==n_dofs_refined</tt>)
+ *
+ * Multiple calling of this function is NOT allowed. Interpolating several
+ * functions can be performed in one step by using <tt>interpolate (all_in,
+ * all_out)</tt>
+ */
+ void
+ interpolate(const VectorType &in, VectorType &out) const;
+
+ /**
+ * Determine an estimate for the memory consumption (in bytes) of this
+ * object.
+ */
+ std::size_t
+ memory_consumption() const;
+
+ /**
+ * Exception
+ */
+ DeclExceptionMsg(
+ ExcNotPrepared,
+ "You are attempting an operation for which this object is "
+ "not prepared. This may be because you either did not call "
+ "one of the prepare_*() functions at all, or because you "
+ "called the wrong one for the operation you are currently "
+ "attempting.");
+
+ /**
+ * Exception
+ */
+ DeclExceptionMsg(
+ ExcAlreadyPrepForRef,
+ "You are attempting to call one of the prepare_*() functions "
+ "of this object to prepare it for an operation for which it "
+ "is already prepared. Specifically, the object was "
+ "previously prepared for pure refinement.");
+
+ /**
+ * Exception
+ */
+ DeclExceptionMsg(
+ ExcAlreadyPrepForCoarseAndRef,
+ "You are attempting to call one of the prepare_*() functions "
+ "of this object to prepare it for an operation for which it "
+ "is already prepared. Specifically, the object was "
+ "previously prepared for both coarsening and refinement.");
+
+ private:
+ /**
+ * Pointer to the degree of freedom handler to work with.
+ */
+ SmartPointer<const DoFHandler<dim, spacedim>,
+ SolutionTransfer<dim, VectorType, spacedim>>
+ dof_handler;
+
+ /**
+ * Stores the number of DoFs before the refinement and/or coarsening.
+ */
+ types::global_dof_index n_dofs_old;
+
+ /**
+ * Declaration of @p PreparationState that denotes the three possible states
+ * of the @p SolutionTransfer: being prepared for 'pure refinement',
+ * prepared for 'coarsening and refinement' or not prepared.
+ */
+ enum PreparationState
+ {
+ /**
+ * The SolutionTransfer is not yet prepared.
+ */
+ none,
+ /**
+ * The SolutionTransfer is prepared for purely refinement.
+ */
+ pure_refinement,
+ /**
+ * The SolutionTransfer is prepared for coarsening and refinement.
+ */
+ coarsening_and_refinement
+ };
+
+ /**
+ * Definition of the respective variable.
+ */
+ PreparationState prepared_for;
+
+
+ /**
+ * Is used for @p prepare_for_refining (of course also for @p
+ * repare_for_refining_and_coarsening) and stores all dof indices of the
+ * cells that'll be refined
+ */
+ std::vector<std::vector<types::global_dof_index>> indices_on_cell;
+
+ /**
+ * All cell data (the dof indices and the dof values) should be accessible
+ * from each cell. As each cell has got only one @p user_pointer, multiple
+ * pointers to the data need to be packetized in a structure. Note that in
+ * our case on each cell either the <tt>vector<unsigned int> indices</tt>
+ * (if the cell will be refined) or the <tt>vector<double> dof_values</tt>
+ * (if the children of this cell will be deleted) is needed, hence one @p
+ * user_pointer should be sufficient, but to allow some error checks and to
+ * preserve the user from making user errors the @p user_pointer will be
+ * 'multiplied' by this structure.
+ */
+ struct Pointerstruct
+ {
+ Pointerstruct()
+ : indices_ptr(nullptr)
+ , dof_values_ptr(nullptr)
+ , active_fe_index(0)
+ {}
+ Pointerstruct(std::vector<types::global_dof_index> *indices_ptr_in,
+ const unsigned int active_fe_index_in = 0)
+ : indices_ptr(indices_ptr_in)
+ , dof_values_ptr(nullptr)
+ , active_fe_index(active_fe_index_in)
+ {}
+ Pointerstruct(
+ std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr_in,
+ const unsigned int active_fe_index_in = 0)
+ : indices_ptr(nullptr)
+ , dof_values_ptr(dof_values_ptr_in)
+ , active_fe_index(active_fe_index_in)
+ {}
+ std::size_t
+ memory_consumption() const;
+
+ std::vector<types::global_dof_index> *indices_ptr;
+ std::vector<Vector<typename VectorType::value_type>> *dof_values_ptr;
+ unsigned int active_fe_index;
+ };
+
+ /**
+ * Map mapping from level and index of cell to the @p Pointerstructs (cf.
+ * there). This map makes it possible to keep all the information needed to
+ * transfer the solution inside this object rather than using user pointers
+ * of the Triangulation for this purpose.
+ */
+ std::map<std::pair<unsigned int, unsigned int>, Pointerstruct> cell_map;
+
+ /**
+ * Is used for @p prepare_for_refining_and_coarsening The interpolated dof
+ * values of all cells that'll be coarsened will be stored in this vector.
+ */
+ std::vector<std::vector<Vector<typename VectorType::value_type>>>
+ dof_values_on_cell;
+ };
+} // namespace Legacy
DEAL_II_NAMESPACE_CLOSE
cell_data_transfer.cc
fully_distributed_tria.cc
repartitioning_policy_tools.cc
- solution_transfer.cc
tria.cc
tria_base.cc
shared_tria.cc
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
std::ofstream f(fname);
f << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_global_active_cells"
<< std::endl
- << 4 << " "
- << Utilities::MPI::n_mpi_processes(this->mpi_communicator) << " "
- << this->cell_attached_data.pack_callbacks_fixed.size() << " "
- << this->cell_attached_data.pack_callbacks_variable.size() << " "
- << this->n_global_active_cells() << std::endl;
+ << ::dealii::internal::CellAttachedDataSerializer<dim, spacedim>::
+ version_number
+ << " " << Utilities::MPI::n_mpi_processes(this->mpi_communicator)
+ << " " << this->cell_attached_data.pack_callbacks_fixed.size()
+ << " " << this->cell_attached_data.pack_callbacks_variable.size()
+ << " " << this->n_global_active_cells() << std::endl;
}
// Save cell attached data.
std::ifstream f(fname);
AssertThrow(f.fail() == false, ExcIO());
std::string firstline;
- getline(f, firstline); // skip first line
+ getline(f, firstline);
f >> version >> numcpus >> attached_count_fixed >>
attached_count_variable >> n_global_active_cells;
}
- AssertThrow(version == 4,
+ const auto expected_version = ::dealii::internal::
+ CellAttachedDataSerializer<dim, spacedim>::version_number;
+
+ AssertThrow(version == expected_version,
ExcMessage("Incompatible version found in .info file."));
// Load description and construct the triangulation.
+++ /dev/null
-// ------------------------------------------------------------------------
-//
-// SPDX-License-Identifier: LGPL-2.1-or-later
-// Copyright (C) 2010 - 2024 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// Part of the source code is dual licensed under Apache-2.0 WITH
-// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
-// governing the source code and code contributions can be found in
-// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
-//
-// ------------------------------------------------------------------------
-
-
-#include <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_P4EST
-
-# include <deal.II/distributed/solution_transfer.h>
-# include <deal.II/distributed/tria.h>
-
-# include <deal.II/dofs/dof_accessor.h>
-# include <deal.II/dofs/dof_tools.h>
-
-# include <deal.II/grid/tria_accessor.h>
-# include <deal.II/grid/tria_iterator.h>
-
-# include <deal.II/lac/block_vector.h>
-# include <deal.II/lac/la_parallel_block_vector.h>
-# include <deal.II/lac/la_parallel_vector.h>
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# include <deal.II/lac/trilinos_parallel_block_vector.h>
-# include <deal.II/lac/trilinos_vector.h>
-# include <deal.II/lac/vector.h>
-
-# include <functional>
-# include <numeric>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-
-namespace
-{
- /**
- * Optimized pack function for values assigned on degrees of freedom.
- *
- * Given that the elements of @p dof_values are stored in consecutive
- * locations, we can just memcpy them. Since floating point values don't
- * compress well, we also waive the compression that the default
- * Utilities::pack() and Utilities::unpack() functions offer.
- */
- template <typename value_type>
- std::vector<char>
- pack_dof_values(std::vector<Vector<value_type>> &dof_values,
- const unsigned int dofs_per_cell)
- {
- for (const auto &values : dof_values)
- {
- AssertDimension(values.size(), dofs_per_cell);
- (void)values;
- }
-
- const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
-
- std::vector<char> buffer(dof_values.size() * bytes_per_entry);
- for (unsigned int i = 0; i < dof_values.size(); ++i)
- std::memcpy(&buffer[i * bytes_per_entry],
- &dof_values[i](0),
- bytes_per_entry);
-
- return buffer;
- }
-
-
-
- /**
- * Optimized unpack function for values assigned on degrees of freedom.
- */
- template <typename value_type>
- std::vector<Vector<value_type>>
- unpack_dof_values(
- const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
- const unsigned int dofs_per_cell)
- {
- const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
- const unsigned int n_elements = data_range.size() / bytes_per_entry;
-
- Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
-
- std::vector<Vector<value_type>> unpacked_data;
- unpacked_data.reserve(n_elements);
- for (unsigned int i = 0; i < n_elements; ++i)
- {
- Vector<value_type> dof_values(dofs_per_cell);
- std::memcpy(&dof_values(0),
- &(*std::next(data_range.begin(), i * bytes_per_entry)),
- bytes_per_entry);
- unpacked_data.emplace_back(std::move(dof_values));
- }
-
- return unpacked_data;
- }
-} // namespace
-
-
-
-namespace parallel
-{
- namespace distributed
- {
- template <int dim, typename VectorType, int spacedim>
- SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
- const DoFHandler<dim, spacedim> &dof,
- const bool average_values)
- : dof_handler(&dof, typeid(*this).name())
- , average_values(average_values)
- , handle(numbers::invalid_unsigned_int)
- {
- Assert(
- (dynamic_cast<
- const parallel::DistributedTriangulationBase<dim, spacedim> *>(
- &dof_handler->get_triangulation()) != nullptr),
- ExcMessage(
- "parallel::distributed::SolutionTransfer requires a parallel::distributed::Triangulation object."));
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::
- prepare_for_coarsening_and_refinement(
- const std::vector<const VectorType *> &all_in)
- {
- for (unsigned int i = 0; i < all_in.size(); ++i)
- Assert(all_in[i]->size() == dof_handler->n_dofs(),
- ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
-
- input_vectors = all_in;
- register_data_attach();
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::register_data_attach()
- {
- // TODO: casting away constness is bad
- parallel::DistributedTriangulationBase<dim, spacedim> *tria =
- (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
- const_cast<dealii::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation())));
- Assert(tria != nullptr, ExcInternalError());
-
- Assert(handle == numbers::invalid_unsigned_int,
- ExcMessage("You can only add one solution per "
- "SolutionTransfer object."));
-
- handle = tria->register_data_attach(
- [this](
- const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
- const CellStatus status) {
- return this->pack_callback(cell_, status);
- },
- /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::
- prepare_for_coarsening_and_refinement(const VectorType &in)
- {
- std::vector<const VectorType *> all_in(1, &in);
- prepare_for_coarsening_and_refinement(all_in);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
- const VectorType &in)
- {
- std::vector<const VectorType *> all_in(1, &in);
- prepare_for_serialization(all_in);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
- const std::vector<const VectorType *> &all_in)
- {
- prepare_for_coarsening_and_refinement(all_in);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::deserialize(VectorType &in)
- {
- std::vector<VectorType *> all_in(1, &in);
- deserialize(all_in);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::deserialize(
- std::vector<VectorType *> &all_in)
- {
- register_data_attach();
-
- // this makes interpolate() happy
- input_vectors.resize(all_in.size());
-
- interpolate(all_in);
- }
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::interpolate(
- std::vector<VectorType *> &all_out)
- {
- Assert(input_vectors.size() == all_out.size(),
- ExcDimensionMismatch(input_vectors.size(), all_out.size()));
- for (unsigned int i = 0; i < all_out.size(); ++i)
- Assert(all_out[i]->size() == dof_handler->n_dofs(),
- ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
-
- // TODO: casting away constness is bad
- parallel::DistributedTriangulationBase<dim, spacedim> *tria =
- (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
- const_cast<dealii::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation())));
- Assert(tria != nullptr, ExcInternalError());
-
- if (average_values)
- for (auto *const vec : all_out)
- *vec = 0.0;
-
- VectorType valence;
-
- // initialize valence vector only if we need to average
- if (average_values)
- valence.reinit(*all_out[0]);
-
- tria->notify_ready_to_unpack(
- handle,
- [this, &all_out, &valence](
- const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
- const CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- &data_range) {
- this->unpack_callback(cell_, status, data_range, all_out, valence);
- });
-
- if (average_values)
- {
- // finalize valence: compress and invert
- using Number = typename VectorType::value_type;
- valence.compress(VectorOperation::add);
- for (const auto i : valence.locally_owned_elements())
- valence[i] = (static_cast<Number>(valence[i]) == Number() ?
- Number() :
- (Number(1.0) / static_cast<Number>(valence[i])));
- valence.compress(VectorOperation::insert);
-
- for (auto *const vec : all_out)
- {
- // compress and weight with valence
- vec->compress(VectorOperation::add);
- vec->scale(valence);
- }
- }
- else
- {
- for (auto *const vec : all_out)
- vec->compress(VectorOperation::insert);
- }
-
- input_vectors.clear();
- handle = numbers::invalid_unsigned_int;
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::interpolate(VectorType &out)
- {
- std::vector<VectorType *> all_out(1, &out);
- interpolate(all_out);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- std::vector<char>
- SolutionTransfer<dim, VectorType, spacedim>::pack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
- const CellStatus status)
- {
- typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
- dof_handler);
-
- // create buffer for each individual object
- std::vector<::dealii::Vector<typename VectorType::value_type>> dof_values(
- input_vectors.size());
-
- unsigned int fe_index = 0;
- if (dof_handler->has_hp_capabilities())
- {
- switch (status)
- {
- case CellStatus::cell_will_persist:
- case CellStatus::cell_will_be_refined:
- {
- fe_index = cell->future_fe_index();
- break;
- }
-
- case CellStatus::children_will_be_coarsened:
- {
- // In case of coarsening, we need to find a suitable FE index
- // for the parent cell. We choose the 'least dominant fe'
- // on all children from the associated FECollection.
-# ifdef DEBUG
- for (const auto &child : cell->child_iterators())
- Assert(child->is_active() && child->coarsen_flag_set(),
- typename dealii::Triangulation<
- dim>::ExcInconsistentCoarseningFlags());
-# endif
-
- fe_index = dealii::internal::hp::DoFHandlerImplementation::
- dominated_future_fe_on_children<dim, spacedim>(cell);
- break;
- }
-
- default:
- DEAL_II_ASSERT_UNREACHABLE();
- break;
- }
- }
-
- const unsigned int dofs_per_cell =
- dof_handler->get_fe(fe_index).n_dofs_per_cell();
-
- if (dofs_per_cell == 0)
- return std::vector<char>(); // nothing to do for FE_Nothing
-
- auto it_input = input_vectors.cbegin();
- auto it_output = dof_values.begin();
- for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
- {
- it_output->reinit(dofs_per_cell);
- cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
- }
-
- return pack_dof_values<typename VectorType::value_type>(dof_values,
- dofs_per_cell);
- }
-
-
-
- template <int dim, typename VectorType, int spacedim>
- void
- SolutionTransfer<dim, VectorType, spacedim>::unpack_callback(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
- const CellStatus status,
- const boost::iterator_range<std::vector<char>::const_iterator>
- &data_range,
- std::vector<VectorType *> &all_out,
- VectorType &valence)
- {
- typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
- dof_handler);
-
- unsigned int fe_index = 0;
- if (dof_handler->has_hp_capabilities())
- {
- switch (status)
- {
- case CellStatus::cell_will_persist:
- case CellStatus::children_will_be_coarsened:
- {
- fe_index = cell->active_fe_index();
- break;
- }
-
- case CellStatus::cell_will_be_refined:
- {
- // After refinement, this particular cell is no longer active,
- // and its children have inherited its FE index. However, to
- // unpack the data on the old cell, we need to recover its FE
- // index from one of the children. Just to be sure, we also
- // check if all children have the same FE index.
- fe_index = cell->child(0)->active_fe_index();
- for (unsigned int child_index = 1;
- child_index < cell->n_children();
- ++child_index)
- Assert(cell->child(child_index)->active_fe_index() ==
- fe_index,
- ExcInternalError());
- break;
- }
-
- default:
- DEAL_II_ASSERT_UNREACHABLE();
- break;
- }
- }
-
- const unsigned int dofs_per_cell =
- dof_handler->get_fe(fe_index).n_dofs_per_cell();
-
- if (dofs_per_cell == 0)
- return; // nothing to do for FE_Nothing
-
- const std::vector<::dealii::Vector<typename VectorType::value_type>>
- dof_values =
- unpack_dof_values<typename VectorType::value_type>(data_range,
- dofs_per_cell);
-
- // check if sizes match
- Assert(dof_values.size() == all_out.size(), ExcInternalError());
-
- // check if we have enough dofs provided by the FE object
- // to interpolate the transferred data correctly
- for (auto it_dof_values = dof_values.begin();
- it_dof_values != dof_values.end();
- ++it_dof_values)
- Assert(
- dofs_per_cell == it_dof_values->size(),
- ExcMessage(
- "The transferred data was packed with a different number of dofs than the "
- "currently registered FE object assigned to the DoFHandler has."));
-
- // distribute data for each registered vector on mesh
- auto it_input = dof_values.cbegin();
- auto it_output = all_out.begin();
- for (; it_input != dof_values.cend(); ++it_input, ++it_output)
- if (average_values)
- cell->distribute_local_to_global_by_interpolation(*it_input,
- *(*it_output),
- fe_index);
- else
- cell->set_dof_values_by_interpolation(*it_input,
- *(*it_output),
- fe_index,
- true);
-
- if (average_values)
- {
- // compute valence vector if averaging should be performed
- Vector<typename VectorType::value_type> ones(dofs_per_cell);
- ones = 1.0;
- cell->distribute_local_to_global_by_interpolation(ones,
- valence,
- fe_index);
- }
- }
- } // namespace distributed
-} // namespace parallel
-
-
-// explicit instantiations
-# include "solution_transfer.inst"
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ------------------------------------------------------------------------
-//
-// SPDX-License-Identifier: LGPL-2.1-or-later
-// Copyright (C) 2010 - 2021 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// Part of the source code is dual licensed under Apache-2.0 WITH
-// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
-// governing the source code and code contributions can be found in
-// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
-//
-// ------------------------------------------------------------------------
-
-
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
- {
- namespace parallel
- \{
- namespace distributed
- \{
-#if deal_II_dimension <= deal_II_space_dimension
- template class SolutionTransfer<deal_II_dimension,
- ::dealii::Vector<double>,
- deal_II_space_dimension>;
- template class SolutionTransfer<
- deal_II_dimension,
- ::dealii::LinearAlgebra::distributed::Vector<double>,
- deal_II_space_dimension>;
- template class SolutionTransfer<
- deal_II_dimension,
- ::dealii::LinearAlgebra::distributed::Vector<float>,
- deal_II_space_dimension>;
- template class SolutionTransfer<
- deal_II_dimension,
- ::dealii::LinearAlgebra::distributed::BlockVector<double>,
- deal_II_space_dimension>;
- template class SolutionTransfer<
- deal_II_dimension,
- ::dealii::LinearAlgebra::distributed::BlockVector<float>,
- deal_II_space_dimension>;
-
-
-# ifdef DEAL_II_WITH_PETSC
- template class SolutionTransfer<deal_II_dimension,
- PETScWrappers::MPI::Vector,
- deal_II_space_dimension>;
- template class SolutionTransfer<deal_II_dimension,
- PETScWrappers::MPI::BlockVector,
- deal_II_space_dimension>;
-# endif
-
-# ifdef DEAL_II_WITH_TRILINOS
- template class SolutionTransfer<deal_II_dimension,
- TrilinosWrappers::MPI::Vector,
- deal_II_space_dimension>;
- template class SolutionTransfer<deal_II_dimension,
- TrilinosWrappers::MPI::BlockVector,
- deal_II_space_dimension>;
-# endif
-
-#endif
- \}
- \}
- }
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/tria_base.h>
#include <deal.II/grid/connectivity.h>
#include <deal.II/grid/grid_tools.h>
// can only tolerate one level of coarsening at a time, so
// check that the children are all active
Assert(dealii_cell->is_active() == false, ExcInternalError());
- for (unsigned int c = 0;
- c < GeometryInfo<dim>::max_children_per_cell;
- ++c)
+ for (unsigned int c = 0; c < dealii_cell->n_children(); ++c)
Assert(dealii_cell->child(c)->is_active(),
ExcInternalError());
break;
sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
(variable_size_data_stored ? 1 : 0));
-
// Read header data.
file.read(reinterpret_cast<char *>(sizes_fixed_cumulative.data()),
sizes_fixed_cumulative.size() * sizeof(unsigned int));
void Triangulation<dim, spacedim>::save(const std::string &file_basename) const
{
// Save triangulation information.
- std::ofstream ofs(file_basename + "_triangulation.data");
- boost::archive::text_oarchive oa(ofs, boost::archive::no_header);
- save(oa, 0);
+ {
+ std::ofstream ofs_tria(file_basename + "_triangulation.data");
+ AssertThrow(ofs_tria.fail() == false, ExcIO());
+
+ boost::archive::text_oarchive oa(ofs_tria, boost::archive::no_header);
+ save(oa,
+ internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
+ }
// Save attached data.
{
- std::ofstream ifs(file_basename + ".info");
- ifs
+ std::ofstream ofs_info(file_basename + ".info");
+ ofs_info
<< "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_active_cells"
<< std::endl
<< internal::CellAttachedDataSerializer<dim, spacedim>::version_number
void Triangulation<dim, spacedim>::load(const std::string &file_basename)
{
// Load triangulation information.
- std::ifstream ifs(file_basename + "_triangulation.data");
- AssertThrow(ifs.fail() == false, ExcIO());
+ {
+ std::ifstream ifs_tria(file_basename + "_triangulation.data");
+ AssertThrow(ifs_tria.fail() == false, ExcIO());
- boost::archive::text_iarchive ia(ifs, boost::archive::no_header);
- load(ia, 0);
+ boost::archive::text_iarchive ia(ifs_tria, boost::archive::no_header);
+ load(ia,
+ internal::CellAttachedDataSerializer<dim, spacedim>::version_number);
+ }
// Load attached data.
unsigned int version, numcpus, attached_count_fixed, attached_count_variable,
n_global_active_cells;
{
- std::ifstream ifs(std::string(file_basename) + ".info");
- AssertThrow(ifs.fail() == false, ExcIO());
+ std::ifstream ifs_info(std::string(file_basename) + ".info");
+ AssertThrow(ifs_info.fail() == false, ExcIO());
std::string firstline;
- getline(ifs, firstline);
- ifs >> version >> numcpus >> attached_count_fixed >>
+ std::getline(ifs_info, firstline);
+ ifs_info >> version >> numcpus >> attached_count_fixed >>
attached_count_variable >> n_global_active_cells;
}
+template <int dim, int spacedim>
+DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
+void Triangulation<dim, spacedim>::pack_data_serial()
+{
+ if (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
+ this))
+ return;
+
+ std::vector<CellId> active_cell_old;
+
+ // pack data before triangulation gets updated
+ if (this->cell_attached_data.n_attached_data_sets > 0)
+ {
+ // store old active cells to determine cell status after
+ // coarsening/refinement
+ active_cell_old.reserve(this->n_active_cells());
+
+ for (const auto &cell : this->active_cell_iterators())
+ {
+ const bool children_will_be_coarsened =
+ (cell->level() > 0) && (cell->coarsen_flag_set());
+
+ if (children_will_be_coarsened == false)
+ active_cell_old.emplace_back(cell->id());
+ else
+ {
+ if (cell->parent()->child(0) == cell)
+ active_cell_old.emplace_back(cell->parent()->id());
+ }
+ }
+
+ // update cell relations
+ this->local_cell_relations.clear();
+ this->local_cell_relations.reserve(this->n_global_active_cells());
+
+ std::vector<
+ std::pair<unsigned int,
+ typename internal::CellAttachedDataSerializer<dim, spacedim>::
+ cell_relation_t>>
+ local_cell_relations_tmp;
+
+ for (const auto &cell : this->active_cell_iterators())
+ {
+ if (std::find(active_cell_old.begin(),
+ active_cell_old.end(),
+ cell->id()) != active_cell_old.end())
+ {
+ const unsigned int index =
+ std::distance(active_cell_old.begin(),
+ std::find(active_cell_old.begin(),
+ active_cell_old.end(),
+ cell->id()));
+
+ ::dealii::CellStatus status =
+ cell->refine_flag_set() ?
+ ::dealii::CellStatus::cell_will_be_refined :
+ ::dealii::CellStatus::cell_will_persist;
+
+ local_cell_relations_tmp.emplace_back(
+ index,
+ typename internal::CellAttachedDataSerializer<dim, spacedim>::
+ cell_relation_t{cell, status});
+ }
+ else if (cell->level() > 0 &&
+ std::find(active_cell_old.begin(),
+ active_cell_old.end(),
+ cell->parent()->id()) != active_cell_old.end())
+ {
+ const unsigned int index =
+ std::distance(active_cell_old.begin(),
+ std::find(active_cell_old.begin(),
+ active_cell_old.end(),
+ cell->parent()->id()));
+
+ ::dealii::CellStatus status;
+
+ if (cell->parent()->child_iterator_to_index(cell) == 0)
+ status = ::dealii::CellStatus::children_will_be_coarsened;
+ else
+ status = ::dealii::CellStatus::cell_invalid;
+
+ local_cell_relations_tmp.emplace_back(
+ index,
+ typename internal::CellAttachedDataSerializer<dim, spacedim>::
+ cell_relation_t{cell->parent(), status});
+ }
+ else
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+ std::stable_sort(local_cell_relations_tmp.begin(),
+ local_cell_relations_tmp.end(),
+ [](const auto &a, const auto &b) {
+ return a.first < b.first;
+ });
+
+ for (const auto &tmp : local_cell_relations_tmp)
+ this->local_cell_relations.emplace_back(tmp.second);
+
+ // pack data
+ this->data_serializer.pack_data(
+ this->local_cell_relations,
+ this->cell_attached_data.pack_callbacks_fixed,
+ this->cell_attached_data.pack_callbacks_variable,
+ this->get_communicator());
+
+ // dummy copy of data
+ this->data_serializer.dest_data_fixed =
+ this->data_serializer.src_data_fixed;
+ this->data_serializer.dest_data_variable =
+ this->data_serializer.src_data_variable;
+ this->data_serializer.dest_sizes_variable =
+ this->data_serializer.src_sizes_variable;
+ }
+}
+
+
+
+template <int dim, int spacedim>
+DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
+void Triangulation<dim, spacedim>::unpack_data_serial()
+{
+ if (dynamic_cast<parallel::DistributedTriangulationBase<dim, spacedim> *>(
+ this))
+ return;
+
+ // transfer data after triangulation got updated
+ if (this->cell_attached_data.n_attached_data_sets > 0)
+ {
+ std::vector<typename internal::CellAttachedDataSerializer<dim, spacedim>::
+ cell_relation_t>
+ temp;
+
+ for (const auto &cell : local_cell_relations)
+ {
+ if (cell.first->has_children())
+ {
+ Assert(cell.second == ::dealii::CellStatus::cell_will_be_refined,
+ ExcInternalError());
+
+ temp.emplace_back(cell.first->child(0),
+ ::dealii::CellStatus::cell_will_be_refined);
+ }
+ else
+ temp.push_back(cell);
+ }
+
+ this->local_cell_relations = temp;
+ }
+}
+
+
+
template <int dim, int spacedim>
DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
void Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
// 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
// cast away constness
auto tria = const_cast<Triangulation<dim, spacedim> *>(this);
+ // each cell should have been flagged `CellStatus::cell_will_persist`
+ for (const auto &cell_rel : this->local_cell_relations)
+ {
+ (void)cell_rel;
+ Assert((cell_rel.second == // cell_status
+ dealii::CellStatus::cell_will_persist),
+ ExcInternalError());
+ }
+
if (this->cell_attached_data.n_attached_data_sets > 0)
{
// pack attached data first
//
// ------------------------------------------------------------------------
+#include <deal.II/base/config.h>
+
#include <deal.II/base/memory_consumption.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe.h>
#include <deal.II/numerics/solution_transfer.h>
+#include <functional>
+#include <numeric>
+
DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+ /**
+ * Optimized pack function for values assigned on degrees of freedom.
+ *
+ * Given that the elements of @p dof_values are stored in consecutive
+ * locations, we can just memcpy them. Since floating point values don't
+ * compress well, we also waive the compression that the default
+ * Utilities::pack() and Utilities::unpack() functions offer.
+ */
+ template <typename value_type>
+ std::vector<char>
+ pack_dof_values(std::vector<Vector<value_type>> &dof_values,
+ const unsigned int dofs_per_cell)
+ {
+ for (const auto &values : dof_values)
+ {
+ AssertDimension(values.size(), dofs_per_cell);
+ (void)values;
+ }
+
+ const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
+
+ std::vector<char> buffer(dof_values.size() * bytes_per_entry);
+ for (unsigned int i = 0; i < dof_values.size(); ++i)
+ std::memcpy(&buffer[i * bytes_per_entry],
+ &dof_values[i](0),
+ bytes_per_entry);
+
+ return buffer;
+ }
+
+
+
+ /**
+ * Optimized unpack function for values assigned on degrees of freedom.
+ */
+ template <typename value_type>
+ std::vector<Vector<value_type>>
+ unpack_dof_values(
+ const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+ const unsigned int dofs_per_cell)
+ {
+ const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
+ const unsigned int n_elements = data_range.size() / bytes_per_entry;
+
+ Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
+
+ std::vector<Vector<value_type>> unpacked_data;
+ unpacked_data.reserve(n_elements);
+ for (unsigned int i = 0; i < n_elements; ++i)
+ {
+ Vector<value_type> dof_values(dofs_per_cell);
+ std::memcpy(&dof_values(0),
+ &(*std::next(data_range.begin(), i * bytes_per_entry)),
+ bytes_per_entry);
+ unpacked_data.emplace_back(std::move(dof_values));
+ }
+
+ return unpacked_data;
+ }
+} // namespace
+
+
+
template <int dim, typename VectorType, int spacedim>
SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
- const DoFHandler<dim, spacedim> &dof)
+ const DoFHandler<dim, spacedim> &dof,
+ const bool average_values)
: dof_handler(&dof, typeid(*this).name())
- , n_dofs_old(0)
- , prepared_for(none)
+ , average_values(average_values)
+ , handle(numbers::invalid_unsigned_int)
+{}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+ prepare_for_coarsening_and_refinement(
+ const std::vector<const VectorType *> &all_in)
{
- Assert(
- (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> *>(
- &dof_handler->get_triangulation()) == nullptr),
- ExcMessage("You are calling the dealii::SolutionTransfer class "
- "with a DoFHandler that is built on a "
- "parallel::distributed::Triangulation. This will not "
- "work for parallel computations. You probably want to "
- "use the parallel::distributed::SolutionTransfer class."));
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
+
+ for (unsigned int i = 0; i < all_in.size(); ++i)
+ Assert(all_in[i]->size() == dof_handler->n_dofs(),
+ ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
+
+ input_vectors = all_in;
+ register_data_attach();
}
template <int dim, typename VectorType, int spacedim>
-SolutionTransfer<dim, VectorType, spacedim>::~SolutionTransfer()
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+ prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
{
- clear();
+ std::vector<const VectorType *> temp(all_in.size());
+
+ for (std::size_t i = 0; i < temp.size(); ++i)
+ temp[i] = &(all_in[i]);
+
+ this->prepare_for_coarsening_and_refinement(temp);
}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::clear()
+SolutionTransfer<dim, VectorType, spacedim>::register_data_attach()
{
- indices_on_cell.clear();
- dof_values_on_cell.clear();
- cell_map.clear();
+ // TODO: casting away constness is bad
+ auto tria = const_cast<dealii::Triangulation<dim, spacedim> *>(
+ &dof_handler->get_triangulation());
+ Assert(tria != nullptr, ExcInternalError());
+
+ Assert(handle == numbers::invalid_unsigned_int,
+ ExcMessage("You can only add one solution per "
+ "SolutionTransfer object."));
+
+ handle = tria->register_data_attach(
+ [this](const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+ const CellStatus status) {
+ return this->pack_callback(cell_, status);
+ },
+ /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
+}
- prepared_for = none;
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::
+ prepare_for_coarsening_and_refinement(const VectorType &in)
+{
+ std::vector<const VectorType *> all_in(1, &in);
+ prepare_for_coarsening_and_refinement(all_in);
}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::prepare_for_pure_refinement()
+SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
+ const VectorType &in)
{
- Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
- Assert(prepared_for != coarsening_and_refinement,
- ExcAlreadyPrepForCoarseAndRef());
-
- clear();
-
- // We need to access dof indices on the entire domain. For
- // parallel::shared::Triangulations, ownership of cells might change. If they
- // allow artificial cells, we need to restore the "true" cell owners
- // temporarily.
- // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
- // the current set of subdomain ids, set subdomain ids to the "true" owner of
- // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
- // and later restore these flags when it is destroyed.
- const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
- spacedim>
- subdomain_modifier(dof_handler->get_triangulation());
-
- const unsigned int n_active_cells =
- dof_handler->get_triangulation().n_active_cells();
- n_dofs_old = dof_handler->n_dofs();
-
- // efficient reallocation of indices_on_cell
- std::vector<std::vector<types::global_dof_index>>(n_active_cells)
- .swap(indices_on_cell);
-
- for (const auto &cell : dof_handler->active_cell_iterators())
- {
- const unsigned int i = cell->active_cell_index();
- indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
- // on each cell store the indices of the
- // dofs. after refining we get the values
- // on the children by taking these
- // indices, getting the respective values
- // out of the data vectors and prolonging
- // them to the children
- cell->get_dof_indices(indices_on_cell[i]);
- cell_map[std::make_pair(cell->level(), cell->index())] =
- Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
- }
- prepared_for = pure_refinement;
+ std::vector<const VectorType *> all_in(1, &in);
+ prepare_for_serialization(all_in);
}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
- const VectorType &in,
- VectorType &out) const
+SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
+ const std::vector<const VectorType *> &all_in)
{
- Assert(prepared_for == pure_refinement, ExcNotPrepared());
- Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
- Assert(out.size() == dof_handler->n_dofs(),
- ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
- Assert(&in != &out,
- ExcMessage("Vectors cannot be used as input and output"
- " at the same time!"));
-
- // We need to access dof indices on the entire domain. For
- // parallel::shared::Triangulations, ownership of cells might change. If they
- // allow artificial cells, we need to restore the "true" cell owners
- // temporarily.
- // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
- // the current set of subdomain ids, set subdomain ids to the "true" owner of
- // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
- // and later restore these flags when it is destroyed.
- const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
- spacedim>
- subdomain_modifier(dof_handler->get_triangulation());
-
- Vector<typename VectorType::value_type> local_values(0);
-
- typename std::map<std::pair<unsigned int, unsigned int>,
- Pointerstruct>::const_iterator pointerstruct,
- cell_map_end = cell_map.end();
-
- for (const auto &cell : dof_handler->cell_iterators())
- {
- pointerstruct =
- cell_map.find(std::make_pair(cell->level(), cell->index()));
-
- if (pointerstruct != cell_map_end)
- // this cell was refined or not
- // touched at all, so we can get
- // the new values by just setting
- // or interpolating to the children,
- // which is both done by one
- // function
- {
- const unsigned int this_fe_index =
- pointerstruct->second.active_fe_index;
- const unsigned int dofs_per_cell =
- cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
- local_values.reinit(dofs_per_cell, true);
-
- // make sure that the size of the stored indices is the same as
- // dofs_per_cell. since we store the desired fe_index, we know
- // what this size should be
- Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
- ExcInternalError());
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- local_values(i) = internal::ElementAccess<VectorType>::get(
- in, (*pointerstruct->second.indices_ptr)[i]);
- cell->set_dof_values_by_interpolation(local_values,
- out,
- this_fe_index,
- true);
- }
- }
+ prepare_for_coarsening_and_refinement(all_in);
}
-namespace internal
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::deserialize(VectorType &in)
{
- /**
- * Generate a table that contains
- * interpolation matrices between
- * each combination of finite
- * elements used in a DoFHandler of
- * some kind. Since not all
- * elements can be interpolated
- * onto each other, the table may
- * contain empty matrices for those
- * combinations of elements for
- * which no such interpolation is
- * implemented.
- */
- template <int dim, int spacedim>
- void
- extract_interpolation_matrices(const DoFHandler<dim, spacedim> &dof,
- dealii::Table<2, FullMatrix<double>> &matrices)
- {
- if (dof.has_hp_capabilities() == false)
- return;
-
- const dealii::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
- matrices.reinit(fe.size(), fe.size());
- for (unsigned int i = 0; i < fe.size(); ++i)
- for (unsigned int j = 0; j < fe.size(); ++j)
- if (i != j)
- {
- matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
- fe[j].n_dofs_per_cell());
-
- // see if we can get the interpolation matrices for this
- // combination of elements. if not, reset the matrix sizes to zero
- // to indicate that this particular combination isn't
- // supported. this isn't an outright error right away since we may
- // never need to actually interpolate between these two elements
- // on actual cells; we simply have to trigger an error if someone
- // actually tries
- try
- {
- fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
- }
- catch (const typename FiniteElement<dim, spacedim>::
- ExcInterpolationNotImplemented &)
- {
- matrices(i, j).reinit(0, 0);
- }
- }
- }
+ std::vector<VectorType *> all_in(1, &in);
+ deserialize(all_in);
+}
- template <int dim, int spacedim>
- void
- restriction_additive(const FiniteElement<dim, spacedim> &,
- std::vector<std::vector<bool>> &)
- {}
- template <int dim, int spacedim>
- void
- restriction_additive(const dealii::hp::FECollection<dim, spacedim> &fe,
- std::vector<std::vector<bool>> &restriction_is_additive)
- {
- restriction_is_additive.resize(fe.size());
- for (unsigned int f = 0; f < fe.size(); ++f)
- {
- restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
- for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
- restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
- }
- }
-} // namespace internal
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::deserialize(
+ std::vector<VectorType *> &all_in)
+{
+ register_data_attach();
+ // this makes interpolate() happy
+ input_vectors.resize(all_in.size());
+
+ interpolate(all_in);
+}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::
- prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+ std::vector<VectorType *> &all_out)
{
- Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
- Assert(prepared_for != coarsening_and_refinement,
- ExcAlreadyPrepForCoarseAndRef());
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
- clear();
- n_dofs_old = dof_handler->n_dofs();
- const unsigned int in_size = all_in.size();
+ Assert(input_vectors.size() == all_out.size(),
+ ExcDimensionMismatch(input_vectors.size(), all_out.size()));
+ for (unsigned int i = 0; i < all_out.size(); ++i)
+ Assert(all_out[i]->size() == dof_handler->n_dofs(),
+ ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
-#ifdef DEBUG
- Assert(in_size != 0,
- ExcMessage("The array of input vectors you pass to this "
- "function has no elements. This is not useful."));
- for (unsigned int i = 0; i < in_size; ++i)
+ // TODO: casting away constness is bad
+ auto tria = const_cast<dealii::Triangulation<dim, spacedim> *>(
+ &dof_handler->get_triangulation());
+ Assert(tria != nullptr, ExcInternalError());
+ Assert(
+ handle != numbers::invalid_unsigned_int,
+ ExcMessage(
+ "You can only call interpolate() once per SolutionTransfer object."));
+
+ if (average_values)
+ for (auto *const vec : all_out)
+ *vec = 0.0;
+
+ VectorType valence;
+
+ // initialize valence vector only if we need to average
+ if (average_values)
+ valence.reinit(*all_out[0]);
+
+ tria->notify_ready_to_unpack(
+ handle,
+ [this, &all_out, &valence](
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+ const CellStatus status,
+ const boost::iterator_range<std::vector<char>::const_iterator>
+ &data_range) {
+ this->unpack_callback(cell_, status, data_range, all_out, valence);
+ });
+
+ if (average_values)
{
- Assert(all_in[i].size() == n_dofs_old,
- ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+ // finalize valence: compress and invert
+ using Number = typename VectorType::value_type;
+ valence.compress(VectorOperation::add);
+ for (const auto i : valence.locally_owned_elements())
+ valence[i] = (static_cast<Number>(valence[i]) == Number() ?
+ Number() :
+ (Number(1.0) / static_cast<Number>(valence[i])));
+ valence.compress(VectorOperation::insert);
+
+ for (auto *const vec : all_out)
+ {
+ // compress and weight with valence
+ vec->compress(VectorOperation::add);
+ vec->scale(valence);
+ }
}
-#endif
-
- // We need to access dof indices on the entire domain. For
- // parallel::shared::Triangulations, ownership of cells might change. If they
- // allow artificial cells, we need to restore the "true" cell owners
- // temporarily.
- // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
- // the current set of subdomain ids, set subdomain ids to the "true" owner of
- // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
- // and later restore these flags when it is destroyed.
- const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
- spacedim>
- subdomain_modifier(dof_handler->get_triangulation());
-
- // first count the number
- // of cells that will be coarsened
- // and that'll stay or be refined
- unsigned int n_cells_to_coarsen = 0;
- unsigned int n_cells_to_stay_or_refine = 0;
- for (const auto &act_cell : dof_handler->active_cell_iterators())
+ else
{
- if (act_cell->coarsen_flag_set())
- ++n_cells_to_coarsen;
- else
- ++n_cells_to_stay_or_refine;
+ for (auto *const vec : all_out)
+ vec->compress(VectorOperation::insert);
}
- Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
- dof_handler->get_triangulation().n_active_cells(),
- ExcInternalError());
-
- unsigned int n_coarsen_fathers = 0;
- for (const auto &cell : dof_handler->cell_iterators())
- if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
- ++n_coarsen_fathers;
- Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
- (void)n_cells_to_coarsen;
-
- // allocate the needed memory. initialize
- // the following arrays in an efficient
- // way, without copying much
- std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
- .swap(indices_on_cell);
-
- std::vector<std::vector<Vector<typename VectorType::value_type>>>(
- n_coarsen_fathers,
- std::vector<Vector<typename VectorType::value_type>>(in_size))
- .swap(dof_values_on_cell);
-
- Table<2, FullMatrix<double>> interpolation_hp;
- std::vector<std::vector<bool>> restriction_is_additive;
-
- internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
- internal::restriction_additive(dof_handler->get_fe_collection(),
- restriction_is_additive);
-
- // we need counters for
- // the 'to_stay_or_refine' cells 'n_sr' and
- // the 'coarsen_fathers' cells 'n_cf',
- unsigned int n_sr = 0, n_cf = 0;
- for (const auto &cell : dof_handler->cell_iterators())
+
+ input_vectors.clear();
+ handle = numbers::invalid_unsigned_int;
+}
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+ std::vector<VectorType> &all_out)
+{
+ std::vector<VectorType *> temp(all_out.size());
+
+ for (std::size_t i = 0; i < temp.size(); ++i)
+ temp[i] = &(all_out[i]);
+
+ this->interpolate(temp);
+}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+void
+SolutionTransfer<dim, VectorType, spacedim>::interpolate(VectorType &out)
+{
+ std::vector<VectorType *> all_out(1, &out);
+ interpolate(all_out);
+}
+
+
+
+template <int dim, typename VectorType, int spacedim>
+std::vector<char>
+SolutionTransfer<dim, VectorType, spacedim>::pack_callback(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+ const CellStatus status)
+{
+ typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
+
+ // create buffer for each individual object
+ std::vector<::dealii::Vector<typename VectorType::value_type>> dof_values(
+ input_vectors.size());
+
+ unsigned int fe_index = 0;
+ if (dof_handler->has_hp_capabilities())
{
- // CASE 1: active cell that remains as it is
- if (cell->is_active() && !cell->coarsen_flag_set())
+ switch (status)
{
- const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
- indices_on_cell[n_sr].resize(dofs_per_cell);
- // cell will not be coarsened,
- // so we get away by storing the
- // dof indices and later
- // interpolating to the children
- cell->get_dof_indices(indices_on_cell[n_sr]);
- cell_map[std::make_pair(cell->level(), cell->index())] =
- Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
- ++n_sr;
- }
+ case CellStatus::cell_will_persist:
+ case CellStatus::cell_will_be_refined:
+ {
+ fe_index = cell->future_fe_index();
+ break;
+ }
- // CASE 2: cell is inactive but will become active
- else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
- {
- // we will need to interpolate from the children of this cell
- // to the current one. in the hp-context, this also means
- // we need to figure out which finite element space to interpolate
- // to since that is not implied by the global FE as in the non-hp-
- // case. we choose the 'least dominant fe' on all children from
- // the associated FECollection.
- std::set<unsigned int> fe_indices_children;
- for (const auto &child : cell->child_iterators())
+ case CellStatus::children_will_be_coarsened:
{
- Assert(child->is_active() && child->coarsen_flag_set(),
- typename dealii::Triangulation<
- dim>::ExcInconsistentCoarseningFlags());
+ // In case of coarsening, we need to find a suitable FE index
+ // for the parent cell. We choose the 'least dominant fe'
+ // on all children from the associated FECollection.
+#ifdef DEBUG
+ for (const auto &child : cell->child_iterators())
+ Assert(child->is_active() && child->coarsen_flag_set(),
+ typename dealii::Triangulation<
+ dim>::ExcInconsistentCoarseningFlags());
+#endif
- fe_indices_children.insert(child->active_fe_index());
+ fe_index = dealii::internal::hp::DoFHandlerImplementation::
+ dominated_future_fe_on_children<dim, spacedim>(cell);
+ break;
}
- Assert(!fe_indices_children.empty(), ExcInternalError());
-
- const unsigned int target_fe_index =
- dof_handler->get_fe_collection().find_dominated_fe_extended(
- fe_indices_children, /*codim=*/0);
-
- Assert(target_fe_index != numbers::invalid_unsigned_int,
- internal::hp::DoFHandlerImplementation::
- ExcNoDominatedFiniteElementOnChildren());
-
- const unsigned int dofs_per_cell =
- dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
-
- std::vector<Vector<typename VectorType::value_type>>(
- in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
- .swap(dof_values_on_cell[n_cf]);
-
-
- // store the data of each of the input vectors. get this data
- // as interpolated onto a finite element space that encompasses
- // that of all the children. note that
- // cell->get_interpolated_dof_values already does all of the
- // interpolations between spaces
- for (unsigned int j = 0; j < in_size; ++j)
- cell->get_interpolated_dof_values(all_in[j],
- dof_values_on_cell[n_cf][j],
- target_fe_index);
- cell_map[std::make_pair(cell->level(), cell->index())] =
- Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
- ++n_cf;
+
+ default:
+ Assert(false, ExcInternalError());
+ break;
}
}
- Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
- Assert(n_cf == n_coarsen_fathers, ExcInternalError());
- prepared_for = coarsening_and_refinement;
+ const unsigned int dofs_per_cell =
+ dof_handler->get_fe(fe_index).n_dofs_per_cell();
+
+ if (dofs_per_cell == 0)
+ return std::vector<char>(); // nothing to do for FE_Nothing
+
+ auto it_input = input_vectors.cbegin();
+ auto it_output = dof_values.begin();
+ for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
+ {
+ it_output->reinit(dofs_per_cell);
+ cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
+ }
+
+ return pack_dof_values<typename VectorType::value_type>(dof_values,
+ dofs_per_cell);
}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::
- prepare_for_coarsening_and_refinement(const VectorType &in)
+SolutionTransfer<dim, VectorType, spacedim>::unpack_callback(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
+ const CellStatus status,
+ const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
+ std::vector<VectorType *> &all_out,
+ VectorType &valence)
{
- std::vector<VectorType> all_in(1, in);
- prepare_for_coarsening_and_refinement(all_in);
+ typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
+
+ unsigned int fe_index = 0;
+ if (dof_handler->has_hp_capabilities())
+ {
+ switch (status)
+ {
+ case CellStatus::cell_will_persist:
+ case CellStatus::children_will_be_coarsened:
+ {
+ fe_index = cell->active_fe_index();
+ break;
+ }
+
+ case CellStatus::cell_will_be_refined:
+ {
+ // After refinement, this particular cell is no longer active,
+ // and its children have inherited its FE index. However, to
+ // unpack the data on the old cell, we need to recover its FE
+ // index from one of the children. Just to be sure, we also
+ // check if all children have the same FE index.
+ fe_index = cell->child(0)->active_fe_index();
+ for (unsigned int child_index = 1;
+ child_index < cell->n_children();
+ ++child_index)
+ Assert(cell->child(child_index)->active_fe_index() == fe_index,
+ ExcInternalError());
+ break;
+ }
+
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+ }
+
+ const unsigned int dofs_per_cell =
+ dof_handler->get_fe(fe_index).n_dofs_per_cell();
+
+ if (dofs_per_cell == 0)
+ return; // nothing to do for FE_Nothing
+
+ const std::vector<::dealii::Vector<typename VectorType::value_type>>
+ dof_values =
+ unpack_dof_values<typename VectorType::value_type>(data_range,
+ dofs_per_cell);
+
+ // check if sizes match
+ AssertDimension(dof_values.size(), all_out.size());
+
+ // check if we have enough dofs provided by the FE object
+ // to interpolate the transferred data correctly
+ for (auto it_dof_values = dof_values.begin();
+ it_dof_values != dof_values.end();
+ ++it_dof_values)
+ Assert(
+ dofs_per_cell == it_dof_values->size(),
+ ExcMessage(
+ "The transferred data was packed with a different number of dofs than the "
+ "currently registered FE object assigned to the DoFHandler has."));
+
+ // distribute data for each registered vector on mesh
+ auto it_input = dof_values.cbegin();
+ auto it_output = all_out.begin();
+ for (; it_input != dof_values.cend(); ++it_input, ++it_output)
+ if (average_values)
+ cell->distribute_local_to_global_by_interpolation(*it_input,
+ *(*it_output),
+ fe_index);
+ else
+ cell->set_dof_values_by_interpolation(*it_input,
+ *(*it_output),
+ fe_index,
+ true);
+
+ if (average_values)
+ {
+ // compute valence vector if averaging should be performed
+ Vector<typename VectorType::value_type> ones(dofs_per_cell);
+ ones = 1.0;
+ cell->distribute_local_to_global_by_interpolation(ones,
+ valence,
+ fe_index);
+ }
}
template <int dim, typename VectorType, int spacedim>
void
-SolutionTransfer<dim, VectorType, spacedim>::interpolate(
- const std::vector<VectorType> &all_in,
- std::vector<VectorType> &all_out) const
+SolutionTransfer<dim, VectorType, spacedim>::clear()
{
- const unsigned int size = all_in.size();
-#ifdef DEBUG
- Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
- Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
- for (unsigned int i = 0; i < size; ++i)
- Assert(all_in[i].size() == n_dofs_old,
- ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
- for (unsigned int i = 0; i < all_out.size(); ++i)
- Assert(all_out[i].size() == dof_handler->n_dofs(),
- ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
- for (unsigned int i = 0; i < size; ++i)
- for (unsigned int j = 0; j < size; ++j)
- Assert(&all_in[i] != &all_out[j],
- ExcMessage("Vectors cannot be used as input and output"
- " at the same time!"));
-#endif
+ // nothing to do
+}
- // We need to access dof indices on the entire domain. For
- // parallel::shared::Triangulations, ownership of cells might change. If they
- // allow artificial cells, we need to restore the "true" cell owners
- // temporarily.
- // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
- // the current set of subdomain ids, set subdomain ids to the "true" owner of
- // each cell upon construction of the TemporarilyRestoreSubdomainIds object,
- // and later restore these flags when it is destroyed.
- const internal::parallel::shared::TemporarilyRestoreSubdomainIds<dim,
- spacedim>
- subdomain_modifier(dof_handler->get_triangulation());
-
- Vector<typename VectorType::value_type> local_values;
- std::vector<types::global_dof_index> dofs;
-
- typename std::map<std::pair<unsigned int, unsigned int>,
- Pointerstruct>::const_iterator pointerstruct,
- cell_map_end = cell_map.end();
-
- Table<2, FullMatrix<double>> interpolation_hp;
- internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
- Vector<typename VectorType::value_type> tmp, tmp2;
-
- for (const auto &cell : dof_handler->cell_iterators())
- {
- pointerstruct =
- cell_map.find(std::make_pair(cell->level(), cell->index()));
- if (pointerstruct != cell_map_end)
- {
- const std::vector<types::global_dof_index> *const indexptr =
- pointerstruct->second.indices_ptr;
- const std::vector<Vector<typename VectorType::value_type>>
- *const valuesptr = pointerstruct->second.dof_values_ptr;
+namespace Legacy
+{
- // cell stayed as it was or was refined
- if (indexptr != nullptr)
- {
- Assert(valuesptr == nullptr, ExcInternalError());
+ template <int dim, typename VectorType, int spacedim>
+ SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
+ const DoFHandler<dim, spacedim> &dof)
+ : dof_handler(&dof, typeid(*this).name())
+ , n_dofs_old(0)
+ , prepared_for(none)
+ {
+ Assert((dynamic_cast<
+ const parallel::distributed::Triangulation<dim, spacedim> *>(
+ &dof_handler->get_triangulation()) == nullptr),
+ ExcMessage(
+ "You are calling the dealii::SolutionTransfer class "
+ "with a DoFHandler that is built on a "
+ "parallel::distributed::Triangulation. This will not "
+ "work for parallel computations. You probably want to "
+ "use the parallel::distributed::SolutionTransfer class."));
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ SolutionTransfer<dim, VectorType, spacedim>::~SolutionTransfer()
+ {
+ clear();
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::clear()
+ {
+ indices_on_cell.clear();
+ dof_values_on_cell.clear();
+ cell_map.clear();
+
+ prepared_for = none;
+ }
- const unsigned int old_fe_index =
- pointerstruct->second.active_fe_index;
- // get the values of each of the input data vectors on this cell
- // and prolong it to its children
- unsigned int in_size = indexptr->size();
- for (unsigned int j = 0; j < size; ++j)
+
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::prepare_for_pure_refinement()
+ {
+ Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
+ Assert(prepared_for != coarsening_and_refinement,
+ ExcAlreadyPrepForCoarseAndRef());
+
+ clear();
+
+ // We need to access dof indices on the entire domain. For
+ // parallel::shared::Triangulations, ownership of cells might change. If
+ // they allow artificial cells, we need to restore the "true" cell owners
+ // temporarily.
+ // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+ // the current set of subdomain ids, set subdomain ids to the "true" owner
+ // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+ // object, and later restore these flags when it is destroyed.
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
+
+ const unsigned int n_active_cells =
+ dof_handler->get_triangulation().n_active_cells();
+ n_dofs_old = dof_handler->n_dofs();
+
+ // efficient reallocation of indices_on_cell
+ std::vector<std::vector<types::global_dof_index>>(n_active_cells)
+ .swap(indices_on_cell);
+
+ for (const auto &cell : dof_handler->active_cell_iterators())
+ {
+ const unsigned int i = cell->active_cell_index();
+ indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
+ // on each cell store the indices of the
+ // dofs. after refining we get the values
+ // on the children by taking these
+ // indices, getting the respective values
+ // out of the data vectors and prolonging
+ // them to the children
+ cell->get_dof_indices(indices_on_cell[i]);
+ cell_map[std::make_pair(cell->level(), cell->index())] =
+ Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
+ }
+ prepared_for = pure_refinement;
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
+ const VectorType &in,
+ VectorType &out) const
+ {
+ Assert(prepared_for == pure_refinement, ExcNotPrepared());
+ Assert(in.size() == n_dofs_old,
+ ExcDimensionMismatch(in.size(), n_dofs_old));
+ Assert(out.size() == dof_handler->n_dofs(),
+ ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+ Assert(&in != &out,
+ ExcMessage("Vectors cannot be used as input and output"
+ " at the same time!"));
+
+ // We need to access dof indices on the entire domain. For
+ // parallel::shared::Triangulations, ownership of cells might change. If
+ // they allow artificial cells, we need to restore the "true" cell owners
+ // temporarily.
+ // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+ // the current set of subdomain ids, set subdomain ids to the "true" owner
+ // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+ // object, and later restore these flags when it is destroyed.
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
+
+ Vector<typename VectorType::value_type> local_values(0);
+
+ typename std::map<std::pair<unsigned int, unsigned int>,
+ Pointerstruct>::const_iterator pointerstruct,
+ cell_map_end = cell_map.end();
+
+ for (const auto &cell : dof_handler->cell_iterators())
+ {
+ pointerstruct =
+ cell_map.find(std::make_pair(cell->level(), cell->index()));
+
+ if (pointerstruct != cell_map_end)
+ // this cell was refined or not
+ // touched at all, so we can get
+ // the new values by just setting
+ // or interpolating to the children,
+ // which is both done by one
+ // function
+ {
+ const unsigned int this_fe_index =
+ pointerstruct->second.active_fe_index;
+ const unsigned int dofs_per_cell =
+ cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
+ local_values.reinit(dofs_per_cell, true);
+
+ // make sure that the size of the stored indices is the same as
+ // dofs_per_cell. since we store the desired fe_index, we know
+ // what this size should be
+ Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
+ ExcInternalError());
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ local_values(i) =
+ dealii::internal::ElementAccess<VectorType>::get(
+ in, (*pointerstruct->second.indices_ptr)[i]);
+ cell->set_dof_values_by_interpolation(local_values,
+ out,
+ this_fe_index,
+ true);
+ }
+ }
+ }
+
+
+
+ namespace internal
+ {
+ /**
+ * Generate a table that contains
+ * interpolation matrices between
+ * each combination of finite
+ * elements used in a DoFHandler of
+ * some kind. Since not all
+ * elements can be interpolated
+ * onto each other, the table may
+ * contain empty matrices for those
+ * combinations of elements for
+ * which no such interpolation is
+ * implemented.
+ */
+ template <int dim, int spacedim>
+ void
+ extract_interpolation_matrices(
+ const DoFHandler<dim, spacedim> &dof,
+ dealii::Table<2, FullMatrix<double>> &matrices)
+ {
+ if (dof.has_hp_capabilities() == false)
+ return;
+
+ const dealii::hp::FECollection<dim, spacedim> &fe =
+ dof.get_fe_collection();
+ matrices.reinit(fe.size(), fe.size());
+ for (unsigned int i = 0; i < fe.size(); ++i)
+ for (unsigned int j = 0; j < fe.size(); ++j)
+ if (i != j)
+ {
+ matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
+ fe[j].n_dofs_per_cell());
+
+ // see if we can get the interpolation matrices for this
+ // combination of elements. if not, reset the matrix sizes to zero
+ // to indicate that this particular combination isn't
+ // supported. this isn't an outright error right away since we may
+ // never need to actually interpolate between these two elements
+ // on actual cells; we simply have to trigger an error if someone
+ // actually tries
+ try
{
- tmp.reinit(in_size, true);
- for (unsigned int i = 0; i < in_size; ++i)
- tmp(i) =
- internal::ElementAccess<VectorType>::get(all_in[j],
- (*indexptr)[i]);
-
- cell->set_dof_values_by_interpolation(tmp,
- all_out[j],
- old_fe_index,
- true);
+ fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
}
- }
- else if (valuesptr)
- // the children of this cell were deleted
- {
- Assert(!cell->has_children(), ExcInternalError());
- Assert(indexptr == nullptr, ExcInternalError());
-
- const unsigned int dofs_per_cell =
- cell->get_fe().n_dofs_per_cell();
- dofs.resize(dofs_per_cell);
- // get the local
- // indices
- cell->get_dof_indices(dofs);
-
- // distribute the stored data to the new vectors
- for (unsigned int j = 0; j < size; ++j)
+ catch (const typename FiniteElement<dim, spacedim>::
+ ExcInterpolationNotImplemented &)
{
- // make sure that the size of the stored indices is the same
- // as dofs_per_cell. this is kind of a test if we use the same
- // FE in the hp-case. to really do that test we would have to
- // store the fe_index of all cells
- const Vector<typename VectorType::value_type> *data = nullptr;
- const unsigned int active_fe_index = cell->active_fe_index();
- if (active_fe_index != pointerstruct->second.active_fe_index)
- {
- const unsigned int old_index =
- pointerstruct->second.active_fe_index;
- const FullMatrix<double> &interpolation_matrix =
- interpolation_hp(active_fe_index, old_index);
- // The interpolation matrix might be empty when using
- // FE_Nothing.
- if (interpolation_matrix.empty())
- tmp.reinit(dofs_per_cell, false);
- else
- {
- tmp.reinit(dofs_per_cell, true);
- AssertDimension((*valuesptr)[j].size(),
- interpolation_matrix.n());
- AssertDimension(tmp.size(), interpolation_matrix.m());
- interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
- }
- data = &tmp;
- }
- else
- data = &(*valuesptr)[j];
-
-
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- internal::ElementAccess<VectorType>::set((*data)(i),
- dofs[i],
- all_out[j]);
+ matrices(i, j).reinit(0, 0);
}
}
- // undefined status
- else
- DEAL_II_ASSERT_UNREACHABLE();
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ restriction_additive(const FiniteElement<dim, spacedim> &,
+ std::vector<std::vector<bool>> &)
+ {}
+
+ template <int dim, int spacedim>
+ void
+ restriction_additive(
+ const dealii::hp::FECollection<dim, spacedim> &fe,
+ std::vector<std::vector<bool>> &restriction_is_additive)
+ {
+ restriction_is_additive.resize(fe.size());
+ for (unsigned int f = 0; f < fe.size(); ++f)
+ {
+ restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
+ for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
+ restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
}
}
+ } // namespace internal
- // We have written into the output vectors. If this was a PETSc vector, for
- // example, then we need to compress these to make future operations safe:
- for (auto &vec : all_out)
- vec.compress(VectorOperation::insert);
-}
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::
+ prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
+ {
+ Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
+ Assert(prepared_for != coarsening_and_refinement,
+ ExcAlreadyPrepForCoarseAndRef());
-template <int dim, typename VectorType, int spacedim>
-void
-SolutionTransfer<dim, VectorType, spacedim>::interpolate(const VectorType &in,
- VectorType &out) const
-{
- Assert(in.size() == n_dofs_old, ExcDimensionMismatch(in.size(), n_dofs_old));
- Assert(out.size() == dof_handler->n_dofs(),
- ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+ clear();
+ n_dofs_old = dof_handler->n_dofs();
+ const unsigned int in_size = all_in.size();
- std::vector<VectorType> all_in = {in};
- std::vector<VectorType> all_out = {out};
+#ifdef DEBUG
+ Assert(in_size != 0,
+ ExcMessage("The array of input vectors you pass to this "
+ "function has no elements. This is not useful."));
+ for (unsigned int i = 0; i < in_size; ++i)
+ {
+ Assert(all_in[i].size() == n_dofs_old,
+ ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+ }
+#endif
- interpolate(all_in, all_out);
+ // We need to access dof indices on the entire domain. For
+ // parallel::shared::Triangulations, ownership of cells might change. If
+ // they allow artificial cells, we need to restore the "true" cell owners
+ // temporarily.
+ // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+ // the current set of subdomain ids, set subdomain ids to the "true" owner
+ // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+ // object, and later restore these flags when it is destroyed.
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
+
+ // first count the number
+ // of cells that will be coarsened
+ // and that'll stay or be refined
+ unsigned int n_cells_to_coarsen = 0;
+ unsigned int n_cells_to_stay_or_refine = 0;
+ for (const auto &act_cell : dof_handler->active_cell_iterators())
+ {
+ if (act_cell->coarsen_flag_set())
+ ++n_cells_to_coarsen;
+ else
+ ++n_cells_to_stay_or_refine;
+ }
+ Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
+ dof_handler->get_triangulation().n_active_cells(),
+ ExcInternalError());
+
+ unsigned int n_coarsen_fathers = 0;
+ for (const auto &cell : dof_handler->cell_iterators())
+ if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
+ ++n_coarsen_fathers;
+ Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
+ (void)n_cells_to_coarsen;
+
+ // allocate the needed memory. initialize
+ // the following arrays in an efficient
+ // way, without copying much
+ std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
+ .swap(indices_on_cell);
+
+ std::vector<std::vector<Vector<typename VectorType::value_type>>>(
+ n_coarsen_fathers,
+ std::vector<Vector<typename VectorType::value_type>>(in_size))
+ .swap(dof_values_on_cell);
+
+ Table<2, FullMatrix<double>> interpolation_hp;
+ std::vector<std::vector<bool>> restriction_is_additive;
+
+ internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
+ internal::restriction_additive(dof_handler->get_fe_collection(),
+ restriction_is_additive);
+
+ // we need counters for
+ // the 'to_stay_or_refine' cells 'n_sr' and
+ // the 'coarsen_fathers' cells 'n_cf',
+ unsigned int n_sr = 0, n_cf = 0;
+ for (const auto &cell : dof_handler->cell_iterators())
+ {
+ // CASE 1: active cell that remains as it is
+ if (cell->is_active() && !cell->coarsen_flag_set())
+ {
+ const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
+ indices_on_cell[n_sr].resize(dofs_per_cell);
+ // cell will not be coarsened,
+ // so we get away by storing the
+ // dof indices and later
+ // interpolating to the children
+ cell->get_dof_indices(indices_on_cell[n_sr]);
+ cell_map[std::make_pair(cell->level(), cell->index())] =
+ Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
+ ++n_sr;
+ }
- out = all_out[0];
-}
+ // CASE 2: cell is inactive but will become active
+ else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
+ {
+ // we will need to interpolate from the children of this cell
+ // to the current one. in the hp-context, this also means
+ // we need to figure out which finite element space to interpolate
+ // to since that is not implied by the global FE as in the non-hp-
+ // case. we choose the 'least dominant fe' on all children from
+ // the associated FECollection.
+ std::set<unsigned int> fe_indices_children;
+ for (const auto &child : cell->child_iterators())
+ {
+ Assert(child->is_active() && child->coarsen_flag_set(),
+ typename dealii::Triangulation<
+ dim>::ExcInconsistentCoarseningFlags());
+ fe_indices_children.insert(child->active_fe_index());
+ }
+ Assert(!fe_indices_children.empty(), ExcInternalError());
+
+ const unsigned int target_fe_index =
+ dof_handler->get_fe_collection().find_dominated_fe_extended(
+ fe_indices_children, /*codim=*/0);
+
+ Assert(target_fe_index != numbers::invalid_unsigned_int,
+ dealii::internal::hp::DoFHandlerImplementation::
+ ExcNoDominatedFiniteElementOnChildren());
+
+ const unsigned int dofs_per_cell =
+ dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
+
+ std::vector<Vector<typename VectorType::value_type>>(
+ in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
+ .swap(dof_values_on_cell[n_cf]);
+
+
+ // store the data of each of the input vectors. get this data
+ // as interpolated onto a finite element space that encompasses
+ // that of all the children. note that
+ // cell->get_interpolated_dof_values already does all of the
+ // interpolations between spaces
+ for (unsigned int j = 0; j < in_size; ++j)
+ cell->get_interpolated_dof_values(all_in[j],
+ dof_values_on_cell[n_cf][j],
+ target_fe_index);
+ cell_map[std::make_pair(cell->level(), cell->index())] =
+ Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
+ ++n_cf;
+ }
+ }
+ Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
+ Assert(n_cf == n_coarsen_fathers, ExcInternalError());
+ prepared_for = coarsening_and_refinement;
+ }
-template <int dim, typename VectorType, int spacedim>
-std::size_t
-SolutionTransfer<dim, VectorType, spacedim>::memory_consumption() const
-{
- // at the moment we do not include the memory
- // consumption of the cell_map as we have no
- // real idea about memory consumption of a
- // std::map
- return (MemoryConsumption::memory_consumption(dof_handler) +
- MemoryConsumption::memory_consumption(n_dofs_old) +
- sizeof(prepared_for) +
- MemoryConsumption::memory_consumption(indices_on_cell) +
- MemoryConsumption::memory_consumption(dof_values_on_cell));
-}
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::
+ prepare_for_coarsening_and_refinement(const VectorType &in)
+ {
+ std::vector<VectorType> all_in(1, in);
+ prepare_for_coarsening_and_refinement(all_in);
+ }
-template <int dim, typename VectorType, int spacedim>
-std::size_t
-SolutionTransfer<dim, VectorType, spacedim>::Pointerstruct::memory_consumption()
- const
-{
- return sizeof(*this);
-}
+
+
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+ const std::vector<VectorType> &all_in,
+ std::vector<VectorType> &all_out) const
+ {
+ const unsigned int size = all_in.size();
+#ifdef DEBUG
+ Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
+ Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
+ for (unsigned int i = 0; i < size; ++i)
+ Assert(all_in[i].size() == n_dofs_old,
+ ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
+ for (unsigned int i = 0; i < all_out.size(); ++i)
+ Assert(all_out[i].size() == dof_handler->n_dofs(),
+ ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
+ for (unsigned int i = 0; i < size; ++i)
+ for (unsigned int j = 0; j < size; ++j)
+ Assert(&all_in[i] != &all_out[j],
+ ExcMessage("Vectors cannot be used as input and output"
+ " at the same time!"));
+#endif
+
+ // We need to access dof indices on the entire domain. For
+ // parallel::shared::Triangulations, ownership of cells might change. If
+ // they allow artificial cells, we need to restore the "true" cell owners
+ // temporarily.
+ // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
+ // the current set of subdomain ids, set subdomain ids to the "true" owner
+ // of each cell upon construction of the TemporarilyRestoreSubdomainIds
+ // object, and later restore these flags when it is destroyed.
+ const dealii::internal::parallel::shared::
+ TemporarilyRestoreSubdomainIds<dim, spacedim>
+ subdomain_modifier(dof_handler->get_triangulation());
+
+ Vector<typename VectorType::value_type> local_values;
+ std::vector<types::global_dof_index> dofs;
+
+ typename std::map<std::pair<unsigned int, unsigned int>,
+ Pointerstruct>::const_iterator pointerstruct,
+ cell_map_end = cell_map.end();
+
+ Table<2, FullMatrix<double>> interpolation_hp;
+ internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
+ Vector<typename VectorType::value_type> tmp, tmp2;
+
+ for (const auto &cell : dof_handler->cell_iterators())
+ {
+ pointerstruct =
+ cell_map.find(std::make_pair(cell->level(), cell->index()));
+
+ if (pointerstruct != cell_map_end)
+ {
+ const std::vector<types::global_dof_index> *const indexptr =
+ pointerstruct->second.indices_ptr;
+
+ const std::vector<Vector<typename VectorType::value_type>>
+ *const valuesptr = pointerstruct->second.dof_values_ptr;
+
+ // cell stayed as it was or was refined
+ if (indexptr != nullptr)
+ {
+ Assert(valuesptr == nullptr, ExcInternalError());
+
+ const unsigned int old_fe_index =
+ pointerstruct->second.active_fe_index;
+
+ // get the values of each of the input data vectors on this cell
+ // and prolong it to its children
+ unsigned int in_size = indexptr->size();
+ for (unsigned int j = 0; j < size; ++j)
+ {
+ tmp.reinit(in_size, true);
+ for (unsigned int i = 0; i < in_size; ++i)
+ tmp(i) = dealii::internal::ElementAccess<VectorType>::get(
+ all_in[j], (*indexptr)[i]);
+
+ cell->set_dof_values_by_interpolation(tmp,
+ all_out[j],
+ old_fe_index,
+ true);
+ }
+ }
+ else if (valuesptr)
+ // the children of this cell were deleted
+ {
+ Assert(!cell->has_children(), ExcInternalError());
+ Assert(indexptr == nullptr, ExcInternalError());
+
+ const unsigned int dofs_per_cell =
+ cell->get_fe().n_dofs_per_cell();
+ dofs.resize(dofs_per_cell);
+ // get the local
+ // indices
+ cell->get_dof_indices(dofs);
+
+ // distribute the stored data to the new vectors
+ for (unsigned int j = 0; j < size; ++j)
+ {
+ // make sure that the size of the stored indices is the same
+ // as dofs_per_cell. this is kind of a test if we use the
+ // same FE in the hp-case. to really do that test we would
+ // have to store the fe_index of all cells
+ const Vector<typename VectorType::value_type> *data =
+ nullptr;
+ const unsigned int active_fe_index =
+ cell->active_fe_index();
+ if (active_fe_index !=
+ pointerstruct->second.active_fe_index)
+ {
+ const unsigned int old_index =
+ pointerstruct->second.active_fe_index;
+ const FullMatrix<double> &interpolation_matrix =
+ interpolation_hp(active_fe_index, old_index);
+ // The interpolation matrix might be empty when using
+ // FE_Nothing.
+ if (interpolation_matrix.empty())
+ tmp.reinit(dofs_per_cell, false);
+ else
+ {
+ tmp.reinit(dofs_per_cell, true);
+ AssertDimension((*valuesptr)[j].size(),
+ interpolation_matrix.n());
+ AssertDimension(tmp.size(),
+ interpolation_matrix.m());
+ interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
+ }
+ data = &tmp;
+ }
+ else
+ data = &(*valuesptr)[j];
+
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ dealii::internal::ElementAccess<VectorType>::set(
+ (*data)(i), dofs[i], all_out[j]);
+ }
+ }
+ // undefined status
+ else
+ Assert(false, ExcInternalError());
+ }
+ }
+
+ // We have written into the output vectors. If this was a PETSc vector, for
+ // example, then we need to compress these to make future operations safe:
+ for (auto &vec : all_out)
+ vec.compress(VectorOperation::insert);
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ void
+ SolutionTransfer<dim, VectorType, spacedim>::interpolate(
+ const VectorType &in,
+ VectorType &out) const
+ {
+ Assert(in.size() == n_dofs_old,
+ ExcDimensionMismatch(in.size(), n_dofs_old));
+ Assert(out.size() == dof_handler->n_dofs(),
+ ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
+
+ std::vector<VectorType> all_in = {in};
+ std::vector<VectorType> all_out = {out};
+
+ interpolate(all_in, all_out);
+
+ out = all_out[0];
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ std::size_t
+ SolutionTransfer<dim, VectorType, spacedim>::memory_consumption() const
+ {
+ // at the moment we do not include the memory
+ // consumption of the cell_map as we have no
+ // real idea about memory consumption of a
+ // std::map
+ return (MemoryConsumption::memory_consumption(dof_handler) +
+ MemoryConsumption::memory_consumption(n_dofs_old) +
+ sizeof(prepared_for) +
+ MemoryConsumption::memory_consumption(indices_on_cell) +
+ MemoryConsumption::memory_consumption(dof_values_on_cell));
+ }
+
+
+
+ template <int dim, typename VectorType, int spacedim>
+ std::size_t
+ SolutionTransfer<dim, VectorType, spacedim>::Pointerstruct::
+ memory_consumption() const
+ {
+ return sizeof(*this);
+ }
+
+} // namespace Legacy
/*-------------- Explicit Instantiations -------------------------------*/
deal_II_space_dimension : SPACE_DIMENSIONS)
{
#if deal_II_dimension <= deal_II_space_dimension
+ template class Legacy::
+ SolutionTransfer<deal_II_dimension, VEC, deal_II_space_dimension>;
+
template class SolutionTransfer<deal_II_dimension,
VEC,
deal_II_space_dimension>;
cell->set_refine_flag(RefinementCase<dim>::cut_x);
tria.prepare_coarsening_and_refinement();
- soltrans.prepare_for_pure_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(solution);
tria.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
Vector<double> new_solution(dof_handler.n_dofs());
- soltrans.refine_interpolate(solution, new_solution);
+ soltrans.interpolate(new_solution);
solution.reinit(dof_handler.n_dofs());
solution = new_solution;
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");
cell->set_refine_flag();
tria.prepare_coarsening_and_refinement();
- q_soltrans.prepare_for_pure_refinement();
- dgq_soltrans.prepare_for_pure_refinement();
+ q_soltrans.prepare_for_coarsening_and_refinement(q_solution);
+ dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_solution);
tria.execute_coarsening_and_refinement();
q_dof_handler.distribute_dofs(fe_q);
dgq_dof_handler.distribute_dofs(fe_dgq);
- Vector<double> tmp_q(q_dof_handler.n_dofs());
- q_soltrans.refine_interpolate(q_solution, tmp_q);
q_solution.reinit(q_dof_handler.n_dofs());
- q_solution = tmp_q;
+ q_soltrans.interpolate(q_solution);
- Vector<double> tmp_dgq(dgq_dof_handler.n_dofs());
- dgq_soltrans.refine_interpolate(dgq_solution, tmp_dgq);
dgq_solution.reinit(dgq_dof_handler.n_dofs());
- dgq_solution = tmp_dgq;
-
+ dgq_soltrans.interpolate(dgq_solution);
q_data_out.clear_data_vectors();
q_data_out.add_data_vector(q_solution, "solution");
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");
// get the interpolated solution
// back
Vector<double> tmp(dh.n_dofs());
- soltrans.interpolate(solution, tmp);
+ soltrans.interpolate(tmp);
deallog << "New values:" << std::endl;
for (unsigned int i = 0; i < tmp.size(); ++i)
#include "../tests.h"
-#include "coarse_grid_common.h"
-
-
template <int dim>
void
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/solution_transfer.h>
+#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_q.h>
-
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_description.h>
-
-#include <deal.II/lac/vector.h>
#include "./tests.h"
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SolutionTransfer::load() and save() for serial iangulations.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "./tests.h"
+
+
+template <int dim>
+class InterpolationFunction : public Function<dim>
+{
+public:
+ InterpolationFunction()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const
+ {
+ return p.norm();
+ }
+};
+
+template <int dim, typename TriangulationType>
+void
+test(TriangulationType &triangulation)
+{
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(FE_Q<dim>(2));
+
+ using VectorType = Vector<double>;
+
+ VectorType vector(dof_handler.n_dofs());
+
+ VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), vector);
+
+ VectorType vector_loaded(dof_handler.n_dofs());
+
+ const std::string filename =
+ "solution_transfer_" + std::to_string(dim) + "d_out";
+
+ {
+ SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+ solution_transfer.prepare_for_serialization(vector);
+
+ triangulation.save(filename);
+ }
+
+ triangulation.clear();
+
+ {
+ triangulation.load(filename);
+ dof_handler.distribute_dofs(FE_Q<dim>(2));
+
+ SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+ solution_transfer.deserialize(vector_loaded);
+ }
+
+ // Verify that error is 0.
+ VectorType error(vector);
+ error.add(-1, vector_loaded);
+
+ deallog << (error.linfty_norm() < 1e-16 ? "PASSED" : "FAILED") << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ deallog.push("2d");
+ {
+ constexpr int dim = 2;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube(triangulation);
+ triangulation.refine_global(3);
+
+ test<dim>(triangulation);
+ }
+ deallog.pop();
+
+ deallog.push("3d");
+ {
+ constexpr int dim = 3;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube(triangulation);
+ triangulation.refine_global(3);
+
+ test<dim>(triangulation);
+ }
+ deallog.pop();
+}
--- /dev/null
+
+DEAL:2d::PASSED
+DEAL:3d::PASSED
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2008 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test SolutionTransfer::interpolate() for serial triangulations.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "./tests.h"
+
+
+template <int dim>
+class InterpolationFunction : public Function<dim>
+{
+public:
+ InterpolationFunction()
+ : Function<dim>(1)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const
+ {
+ return p[0];
+ }
+};
+
+template <int dim>
+void
+test(const unsigned int type)
+{
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube(triangulation);
+ triangulation.refine_global(3);
+
+ const FE_Q<dim> fe(2);
+
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ using VectorType = Vector<double>;
+
+ VectorType vector(dof_handler.n_dofs());
+ VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), vector);
+
+ SolutionTransfer<dim, VectorType> solution_transfer(dof_handler);
+
+ triangulation.prepare_coarsening_and_refinement();
+ solution_transfer.prepare_for_coarsening_and_refinement(vector);
+
+ if (type == 0)
+ {
+ triangulation.refine_global(1);
+ }
+ else if (type == 1)
+ {
+ for (const auto &cell : triangulation.active_cell_iterators())
+ if (cell->center()[0] < 0.5)
+ cell->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+ }
+ else if (type == 2)
+ {
+ for (const auto &cell : triangulation.active_cell_iterators())
+ if (cell->center()[0] < 0.5)
+ cell->set_coarsen_flag();
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+ dof_handler.distribute_dofs(fe);
+
+ vector.reinit(dof_handler.n_dofs());
+
+ solution_transfer.interpolate(vector);
+
+ VectorType error(dof_handler.n_dofs());
+ VectorTools::interpolate(dof_handler, InterpolationFunction<dim>(), error);
+
+ error -= vector;
+
+ deallog << (error.linfty_norm() < 1e-16 ? "PASSED" : "FAILED") << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ deallog.push("2d");
+ {
+ constexpr int dim = 2;
+
+ test<dim>(0);
+ test<dim>(1);
+ test<dim>(2);
+ }
+ deallog.pop();
+
+ deallog.push("3d");
+ {
+ constexpr int dim = 3;
+
+ test<dim>(0);
+ test<dim>(1);
+ test<dim>(2);
+ }
+ deallog.pop();
+}
--- /dev/null
+
+DEAL:2d::PASSED
+DEAL:2d::PASSED
+DEAL:2d::PASSED
+DEAL:3d::PASSED
+DEAL:3d::PASSED
+DEAL:3d::PASSED
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
hp::QCollection<2> q;
q.push_back(QMidpoint<2>());
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
hp::QCollection<2> q;
q.push_back(QMidpoint<2>());
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
hp::QCollection<2> q;
q.push_back(QMidpoint<2>());
cell->set_refine_flag();
tria.prepare_coarsening_and_refinement();
- soltrans.prepare_for_pure_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(solution);
tria.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
Vector<double> tmp_q(dof_handler.n_dofs());
- soltrans.refine_interpolate(solution, tmp_q);
+ soltrans.interpolate(tmp_q);
solution.reinit(dof_handler.n_dofs());
solution = tmp_q;
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;
}
setup_system();
// 3.7. Solution Transfer finish
- soltrans.interpolate(solution_coarse, solution);
+ soltrans.interpolate(solution);
}
cell->set_coarsen_flag();
}
- Vector<double> q_old_solution = q_solution, dgq_old_solution = dgq_solution;
- tria.prepare_coarsening_and_refinement();
- q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
- dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
- tria.execute_coarsening_and_refinement();
-
counter = 0;
{
typename DoFHandler<dim>::active_cell_iterator
for (; cell != endc; ++cell, ++celldg, ++counter)
{
if (counter > 20 && counter < 90)
- cell->set_active_fe_index(0);
+ cell->set_future_fe_index(0);
else
- cell->set_active_fe_index(Testing::rand() % max_degree);
+ cell->set_future_fe_index(Testing::rand() % max_degree);
if (counter > 20 && counter < 90)
- celldg->set_active_fe_index(0);
+ celldg->set_future_fe_index(0);
else
- celldg->set_active_fe_index(Testing::rand() % max_degree);
+ celldg->set_future_fe_index(Testing::rand() % max_degree);
}
}
+ Vector<double> q_old_solution = q_solution, dgq_old_solution = dgq_solution;
+ tria.prepare_coarsening_and_refinement();
+ q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+ dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
+ tria.execute_coarsening_and_refinement();
+
q_dof_handler.distribute_dofs(fe_q);
dgq_dof_handler.distribute_dofs(fe_dgq);
q_solution.reinit(q_dof_handler.n_dofs());
dgq_solution.reinit(dgq_dof_handler.n_dofs());
- q_soltrans.interpolate(q_old_solution, q_solution);
- dgq_soltrans.interpolate(dgq_old_solution, dgq_solution);
+ q_soltrans.interpolate(q_solution);
+ dgq_soltrans.interpolate(dgq_solution);
// check correctness by comparing the values
// on points of QGauss of order 2.
tria.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
- soltrans.interpolate(old_solution, solution);
+ soltrans.interpolate(solution);
}
// Interpolate solution
SolutionTransfer<2, Vector<double>> solultion_trans(dof_handler);
- solultion_trans.prepare_for_coarsening_and_refinement(solution);
- triangulation.execute_coarsening_and_refinement();
// Assign FE_Q_ to all cells
cell = dof_handler.begin_active();
for (; cell != endc; ++cell)
{
- cell->set_active_fe_index(0);
+ cell->set_future_fe_index(0);
}
+
+ solultion_trans.prepare_for_coarsening_and_refinement(solution);
+
+ triangulation.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
// Save output
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
// we are good if we made it to here
deallog << "OK" << std::endl;
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
// Define compression level for output data
DataOutBase::VtkFlags vtk_flags;
// set refine flag for the only cell we have, then do the refinement
SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
dof_handler.begin_active()->set_refine_flag();
- solution_trans.prepare_for_coarsening_and_refinement(solution);
- triangulation.execute_coarsening_and_refinement();
// now set the active_fe_index flags on the new set of fine level cells
for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
- dof_handler.begin(0)->child(c)->set_active_fe_index(1);
+ dof_handler.begin(0)->child(c)->set_future_fe_index(1);
+
+ solution_trans.prepare_for_coarsening_and_refinement(solution);
+ triangulation.execute_coarsening_and_refinement();
// distribute dofs and transfer solution there
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solution_trans.interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
// we should now have only 1s in the new_solution vector
for (unsigned int i = 0; i < new_solution.size(); ++i)
// set refine flag for the only cell we have, then do the refinement
SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
dof_handler.begin_active()->set_refine_flag();
- solution_trans.prepare_for_pure_refinement();
- triangulation.execute_coarsening_and_refinement();
// now set the active_fe_index flags on the new set of fine level cells
for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
- dof_handler.begin(0)->child(c)->set_active_fe_index(1);
+ dof_handler.begin(0)->child(c)->set_future_fe_index(1);
+
+ solution_trans.prepare_for_coarsening_and_refinement(solution);
+ triangulation.execute_coarsening_and_refinement();
// distribute dofs and transfer solution there
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solution_trans.refine_interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
// we should now have only 1s in the new_solution vector
for (unsigned int i = 0; i < new_solution.size(); ++i)
for (unsigned int c = 0; c < dof_handler.begin(0)->n_children(); ++c)
dof_handler.begin(0)->child(c)->set_coarsen_flag();
+ // now set the active_fe_index flags on the only remaining cell
+ dof_handler.begin_active()->set_future_fe_index(0);
+
solution_trans.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
- // now set the active_fe_index flags on the only remaining cell
- dof_handler.begin_active()->set_active_fe_index(0);
-
// distribute dofs and transfer solution there
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solution_trans.interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
// we should now have only 1s in the new_solution vector
for (unsigned int i = 0; i < new_solution.size(); ++i)
// set refine flag for the only cell we have, then do the refinement
SolutionTransfer<dim, Vector<double>> solution_trans(dof_handler);
dof_handler.begin_active()->set_refine_flag();
- solution_trans.prepare_for_pure_refinement();
+ solution_trans.prepare_for_coarsening_and_refinement(solution);
triangulation.execute_coarsening_and_refinement();
// now set the active_fe_index flags on the new set of fine level cells
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solution_trans.refine_interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
// we should now have only 1s in the new_solution vector
for (unsigned int i = 0; i < new_solution.size(); ++i)
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solultion_trans.interpolate(solution, new_solution);
+ solultion_trans.interpolate(new_solution);
deallog << "OK" << std::endl;
}
Vector<double> old_solution = solution;
tria.prepare_coarsening_and_refinement();
- soltrans.prepare_for_pure_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(solution);
tria.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
- soltrans.refine_interpolate(old_solution, solution);
+ soltrans.interpolate(solution);
}
cell->set_coarsen_flag();
}
- Vector<double> q_old_solution = q_solution;
- tria.prepare_coarsening_and_refinement();
- q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
- tria.execute_coarsening_and_refinement();
-
counter = 0;
{
typename DoFHandler<dim>::active_cell_iterator cell = q_dof_handler
for (; cell != endc; ++cell, ++counter)
{
if (counter > 20 && counter < 90)
- cell->set_active_fe_index(0);
+ cell->set_future_fe_index(0);
else
- cell->set_active_fe_index(Testing::rand() % max_degree);
+ cell->set_future_fe_index(Testing::rand() % max_degree);
}
}
+ Vector<double> q_old_solution = q_solution;
+ tria.prepare_coarsening_and_refinement();
+ q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
+ tria.execute_coarsening_and_refinement();
+
q_dof_handler.distribute_dofs(fe_q);
q_solution.reinit(q_dof_handler.n_dofs());
- q_soltrans.interpolate(q_old_solution, q_solution);
+ q_soltrans.interpolate(q_solution);
// check correctness by comparing the values
// on points of QGauss of order 2.
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution(dof_handler.n_dofs());
- solution_trans.interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
FE_Type.reinit(triangulation.n_active_cells());
cnt_cells = 0;
dof_handler.distribute_dofs(fe_collection);
Vector<double> new_solution2(dof_handler.n_dofs());
- solution_trans2.interpolate(solution, new_solution2);
+ solution_trans2.interpolate(new_solution2);
FE_Type.reinit(triangulation.n_active_cells());
cnt_cells = 0;
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);
}
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);
this_mpi_process);
solution.reinit(locally_owned_partitioning, mpi_communicator);
- soltrans.refine_interpolate(old_solution, solution);
+ soltrans.interpolate(solution);
}
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
triangulation.execute_coarsening_and_refinement();
- // Assign FE_Q(1) to all cells
- for (const auto &cell : dof_handler.active_cell_iterators())
- cell->set_active_fe_index(0);
-
dof_handler.distribute_dofs(fe_collection);
deallog << "Final number of dofs: " << dof_handler.n_dofs() << std::endl;
Vector<double> new_solution(dof_handler.n_dofs());
- new_solution = 1.;
- solution_trans.interpolate(solution, new_solution);
+ solution_trans.interpolate(new_solution);
deallog << "Vector after solution transfer:" << std::endl;
new_solution.print(deallog.get_file_stream());
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 =
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);
}
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);
this_mpi_process);
solution.reinit(locally_owned_partitioning, mpi_communicator);
- soltrans.refine_interpolate(old_solution, solution);
+ soltrans.interpolate(solution);
}
dh.distribute_dofs(fe);
Vector<double> sol_new(dh.n_dofs());
- soltrans.interpolate(sol_old, sol_new);
+ soltrans.interpolate(sol_new);
+
for (unsigned int i = 0; i < sol_new.size(); ++i)
AssertThrow(sol_new[i] == 1., ExcInternalError());
sol_old = 1.;
SolutionTransfer<dim> soltrans(dh);
- soltrans.prepare_for_pure_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(sol_old);
tria.execute_coarsening_and_refinement();
dh.distribute_dofs(fe);
Vector<double> sol_new(dh.n_dofs());
- soltrans.refine_interpolate(sol_old, sol_new);
+ soltrans.interpolate(sol_new);
for (unsigned int i = 0; i < sol_new.size(); ++i)
AssertThrow(sol_new[i] == 1., ExcInternalError());
dh.distribute_dofs(fe);
Vector<double> sol_new(dh.n_dofs());
- soltrans.interpolate(sol_old, sol_new);
+ soltrans.interpolate(sol_new);
return sol_new;
}
VectorTools::interpolate(MappingQ<dim>(1), dh, function, sol_old);
SolutionTransfer<dim> soltrans(dh);
- soltrans.prepare_for_pure_refinement();
+ soltrans.prepare_for_coarsening_and_refinement(sol_old);
tria.execute_coarsening_and_refinement();
dh.distribute_dofs(fe);
Vector<double> sol_new(dh.n_dofs());
- soltrans.refine_interpolate(sol_old, sol_new);
+ soltrans.interpolate(sol_new);
return sol_new;
}
triangulation.execute_coarsening_and_refinement();
- solution_trans.interpolate(previous_solution, solution);
+ solution_trans.interpolate(solution);
deallog << "OK" << std::endl;
}
triangulation.execute_coarsening_and_refinement();
- solution_trans.interpolate(previous_solution, solution);
+ solution_trans.interpolate(solution);
deallog << "OK" << std::endl;
}
std::vector<TrilinosWrappers::MPI::Vector> tmp(2);
tmp[0].reinit(temperature_solution);
tmp[1].reinit(temperature_solution);
- temperature_trans.interpolate(x_temperature, tmp);
+ temperature_trans.interpolate(tmp);
temperature_solution = tmp[0];
old_temperature_solution = tmp[1];
temperature_constraints.distribute(temperature_solution);
temperature_constraints.distribute(old_temperature_solution);
- stokes_trans.interpolate(x_stokes, stokes_solution);
+ stokes_trans.interpolate(stokes_solution);
stokes_constraints.distribute(stokes_solution);
rebuild_stokes_matrix = true;
rebuild_temperature_matrices = true;
dof_handler.distribute_dofs(fe);
Vector<double> tmp(dof_handler.n_dofs());
- solution_transfer.interpolate(current_solution, tmp);
+ solution_transfer.interpolate(tmp);
current_solution = std::move(tmp);
hanging_node_constraints.clear();
dof_handler.distribute_dofs(fe);
Vector<double> tmp(dof_handler.n_dofs());
- solution_transfer.interpolate(current_solution, tmp);
+ solution_transfer.interpolate(tmp);
current_solution = std::move(tmp);
hanging_node_constraints.clear();
projects / dealii.git / commitdiff