From 95fd595dddd23db5e95dc3b00065d6a5e1f63e7d Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 25 Feb 2022 19:18:04 +0100 Subject: [PATCH] Averaging in parallel::distributed::SolutionTransfer --- .../deal.II/distributed/solution_transfer.h | 21 ++++-- include/deal.II/dofs/dof_accessor.h | 17 +++++ source/distributed/solution_transfer.cc | 67 +++++++++++++++---- source/dofs/dof_accessor_set.cc | 27 ++++++++ source/dofs/dof_accessor_set.inst.in | 16 +++++ 5 files changed, 131 insertions(+), 17 deletions(-) diff --git a/include/deal.II/distributed/solution_transfer.h b/include/deal.II/distributed/solution_transfer.h index a2c0e21059..d717c90985 100644 --- a/include/deal.II/distributed/solution_transfer.h +++ b/include/deal.II/distributed/solution_transfer.h @@ -228,12 +228,17 @@ namespace parallel /** * Constructor. * - * @param[in] dof The DoFHandler on which all operations will happen. - * At the time when this constructor is called, the DoFHandler still - * points to the Triangulation before the refinement in question + * @param[in] dof_handler The DoFHandler on which all operations will + * happen. At the time when this constructor is called, the DoFHandler + * still points to the Triangulation before the refinement in question * happens. + * @param[in] average_values Average the contribututions to the same + * DoF coming from different cells. Note: averaging requires an + * additional communication step, since the valence of the DoF has to be + * determined. */ - SolutionTransfer(const DoFHandler &dof); + SolutionTransfer(const DoFHandler &dof_handler, + const bool average_values = false); /** * Destructor. @@ -319,6 +324,11 @@ namespace parallel SolutionTransfer> dof_handler; + /** + * Flag indicating if averaging should be performed. + */ + const bool average_values; + /** * A vector that stores pointers to all the vectors we are supposed to * copy over from the old to the new mesh. @@ -352,7 +362,8 @@ namespace parallel const typename Triangulation::CellStatus status, const boost::iterator_range::const_iterator> & data_range, - std::vector &all_out); + std::vector &all_out, + VectorType & valence); /** diff --git a/include/deal.II/dofs/dof_accessor.h b/include/deal.II/dofs/dof_accessor.h index 2f561c44ed..3f7285bbac 100644 --- a/include/deal.II/dofs/dof_accessor.h +++ b/include/deal.II/dofs/dof_accessor.h @@ -1766,6 +1766,23 @@ public: DoFHandler::invalid_fe_index, const bool perform_check = false) const; + /** + * Similar to set_dof_values_by_interpolation() with the difference that + * values are added into the vector. + * + * @note In parallel::distributed::SolutionTransfer, this function is used + * to accumulate the contributions of all cells to a DoF; with a + * subsequent multiplication with the inverse of the valence, finally, + * the average value is obtained. + */ + template + void + distribute_local_to_global_by_interpolation( + const Vector &local_values, + OutputVector & values, + const unsigned int fe_index = + DoFHandler::invalid_fe_index) const; + /** * Distribute a local (cell based) vector to a global one by mapping the * local numbering of the degrees of freedom to the global one and entering diff --git a/source/distributed/solution_transfer.cc b/source/distributed/solution_transfer.cc index 889658aad0..f9c476c530 100644 --- a/source/distributed/solution_transfer.cc +++ b/source/distributed/solution_transfer.cc @@ -116,8 +116,10 @@ namespace parallel { template SolutionTransfer::SolutionTransfer( - const DoFHandler &dof) + const DoFHandler &dof, + const bool average_values) : dof_handler(&dof, typeid(*this).name()) + , average_values(average_values) , handle(numbers::invalid_unsigned_int) { Assert( @@ -242,20 +244,45 @@ namespace parallel &dof_handler->get_triangulation()))); Assert(tria != nullptr, ExcInternalError()); + for (const auto 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]( + [this, &all_out, &valence]( const typename Triangulation::cell_iterator &cell_, const typename Triangulation::CellStatus status, const boost::iterator_range::const_iterator> &data_range) { - this->unpack_callback(cell_, status, data_range, all_out); + this->unpack_callback(cell_, status, data_range, all_out, valence); }); - for (typename std::vector::iterator it = all_out.begin(); - it != all_out.end(); - ++it) - (*it)->compress(::dealii::VectorOperation::insert); + if (average_values) + { + // finalize valence: compress and invert + valence.compress(::dealii::VectorOperation::add); + for (const auto i : valence.locally_owned_elements()) + valence[i] = valence[i] == 0.0 ? 0.0 : (1.0 / valence[i]); + valence.compress(::dealii::VectorOperation::insert); + + for (const auto vec : all_out) + { + // comress and weight with valence + vec->compress(::dealii::VectorOperation::add); + vec->scale(valence); + } + } + else + { + for (const auto vec : all_out) + vec->compress(::dealii::VectorOperation::insert); + } input_vectors.clear(); } @@ -350,7 +377,8 @@ namespace parallel const typename Triangulation::CellStatus status, const boost::iterator_range::const_iterator> & data_range, - std::vector &all_out) + std::vector &all_out, + VectorType & valence) { typename DoFHandler::cell_iterator cell(*cell_, dof_handler); @@ -422,10 +450,25 @@ namespace parallel auto it_input = dof_values.cbegin(); auto it_output = all_out.begin(); for (; it_input != dof_values.cend(); ++it_input, ++it_output) - cell->set_dof_values_by_interpolation(*it_input, - *(*it_output), - fe_index, - true); + if (average_values) + cell->distribute_local_to_global_by_interpolation(*it_input, + *(*it_output), + fe_index); + else + cell->set_dof_values_by_interpolation(*it_input, + *(*it_output), + fe_index, + true); + + if (average_values) + { + // compute valence vector if averaging should be performed + Vector ones(dofs_per_cell); + ones = 1.0; + cell->distribute_local_to_global_by_interpolation(ones, + valence, + fe_index); + } } } // namespace distributed } // namespace parallel diff --git a/source/dofs/dof_accessor_set.cc b/source/dofs/dof_accessor_set.cc index bd2f50e5c6..20dd8d9685 100644 --- a/source/dofs/dof_accessor_set.cc +++ b/source/dofs/dof_accessor_set.cc @@ -287,6 +287,33 @@ DoFCellAccessor::set_dof_values_by_interpolation( } +template +template +void +DoFCellAccessor:: + distribute_local_to_global_by_interpolation( + const Vector &local_values, + OutputVector & values, + const unsigned int fe_index_) const +{ + internal::process_by_interpolation( + *this, + local_values, + values, + fe_index_, + [](const DoFCellAccessor &cell, + const Vector & local_values, + OutputVector & values) { + std::vector dof_indices( + cell.get_fe().n_dofs_per_cell()); + cell.get_dof_indices(dof_indices); + AffineConstraints().distribute_local_to_global(local_values, + dof_indices, + values); + }); +} + + // -------------------------------------------------------------------------- // explicit instantiations #include "dof_accessor_set.inst" diff --git a/source/dofs/dof_accessor_set.inst.in b/source/dofs/dof_accessor_set.inst.in index 6bdc91315f..fabddab2a7 100644 --- a/source/dofs/dof_accessor_set.inst.in +++ b/source/dofs/dof_accessor_set.inst.in @@ -23,6 +23,11 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL) const unsigned int fe_index, const bool) const; + template void DoFCellAccessor:: + distribute_local_to_global_by_interpolation( + const Vector &, VEC &, const unsigned int fe_index) + const; + #if deal_II_dimension != 3 template void @@ -32,6 +37,12 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL) const unsigned int fe_index, const bool) const; + template void + DoFCellAccessor:: + distribute_local_to_global_by_interpolation( + const Vector &, VEC &, const unsigned int fe_index) + const; + #endif #if deal_II_dimension == 3 @@ -42,5 +53,10 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL) const unsigned int fe_index, const bool) const; + template void + DoFCellAccessor<1, 3, lda>::distribute_local_to_global_by_interpolation( + const Vector &, VEC &, const unsigned int fe_index) + const; + #endif } -- 2.39.5