Merge pull request #13457 from peterrum/distribute_local_to_global_by…

…_interpolation

Allow to average values in parallel::distributed::SolutionTransfer

doc/news/changes/minor/20220514Munch

@@ -0,0 +1,5 @@
+New: We added the possibility to average the contributions to a DoF from 
+different cells in parallel::distributed::SolutionTransfer. For this purpose,
+set the corresponding flag in the constructor. 
+<br>
+(Peter Munch, Magdalena Schreter, 2022/05/14)

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<dim, spacedim> &dof);
+      SolutionTransfer(const DoFHandler<dim, spacedim> &dof_handler,
+                       const bool                       average_values = false);
 
       /**
        * Destructor.
@@ -319,6 +324,11 @@ namespace parallel
                    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.
@@ -352,7 +362,8 @@ namespace parallel
         const typename Triangulation<dim, spacedim>::CellStatus     status,
         const boost::iterator_range<std::vector<char>::const_iterator>
           &                        data_range,
-        std::vector<VectorType *> &all_out);
+        std::vector<VectorType *> &all_out,
+        VectorType &               valence);
 
 
       /**

include/deal.II/dofs/dof_accessor.h

@@ -1767,6 +1767,23 @@ class DoFCellAccessor : public DoFAccessor<dimension_,
       DoFHandler<dimension_, space_dimension_>::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 <class OutputVector, typename number>
+  void
+  distribute_local_to_global_by_interpolation(
+    const Vector<number> &local_values,
+    OutputVector &        values,
+    const unsigned int    fe_index =
+      DoFHandler<dimension_, space_dimension_>::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

source/distributed/solution_transfer.cc

@@ -116,8 +116,10 @@ namespace parallel
   {
     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())
+      , average_values(average_values)
       , handle(numbers::invalid_unsigned_int)
     {
       Assert(
@@ -242,20 +244,46 @@ namespace parallel
             &dof_handler->get_triangulation())));
       Assert(tria != nullptr, ExcInternalError());
 
+      if (average_values)
+        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<dim, spacedim>::cell_iterator &cell_,
           const typename Triangulation<dim, spacedim>::CellStatus     status,
           const boost::iterator_range<std::vector<char>::const_iterator>
             &data_range) {
-          this->unpack_callback(cell_, status, data_range, all_out);
+          this->unpack_callback(cell_, status, data_range, all_out, valence);
         });
 
-      for (typename std::vector<VectorType *>::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)
+            {
+              // compress 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 +378,8 @@ namespace parallel
       const typename Triangulation<dim, spacedim>::CellStatus     status,
       const boost::iterator_range<std::vector<char>::const_iterator>
         &                        data_range,
-      std::vector<VectorType *> &all_out)
+      std::vector<VectorType *> &all_out,
+      VectorType &               valence)
     {
       typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_,
                                                              dof_handler);
@@ -422,10 +451,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<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

source/dofs/dof_accessor_set.cc

@@ -287,6 +287,33 @@ DoFCellAccessor<dim, spacedim, lda>::set_dof_values_by_interpolation(
 }
 
 
+template <int dim, int spacedim, bool lda>
+template <class OutputVector, typename number>
+void
+DoFCellAccessor<dim, spacedim, lda>::
+  distribute_local_to_global_by_interpolation(
+    const Vector<number> &local_values,
+    OutputVector &        values,
+    const unsigned int    fe_index_) const
+{
+  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
+    *this,
+    local_values,
+    values,
+    fe_index_,
+    [](const DoFCellAccessor<dim, spacedim, lda> &cell,
+       const Vector<number> &                     local_values,
+       OutputVector &                             values) {
+      std::vector<types::global_dof_index> dof_indices(
+        cell.get_fe().n_dofs_per_cell());
+      cell.get_dof_indices(dof_indices);
+      AffineConstraints<number>().distribute_local_to_global(local_values,
+                                                             dof_indices,
+                                                             values);
+    });
+}
+
+
 // --------------------------------------------------------------------------
 // explicit instantiations
 #include "dof_accessor_set.inst"

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<deal_II_dimension, deal_II_dimension, lda>::
+      distribute_local_to_global_by_interpolation(
+        const Vector<VEC::value_type> &, 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<deal_II_dimension, deal_II_dimension + 1, lda>::
+      distribute_local_to_global_by_interpolation(
+        const Vector<VEC::value_type> &, 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::value_type> &, VEC &, const unsigned int fe_index)
+      const;
+
 #endif
   }

tests/distributed_grids/solution_transfer_05.cc

@@ -0,0 +1,95 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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 distributed solution_transfer with averaging (without averaging the
+// program fails)
+
+#include <deal.II/base/mpi.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/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+int
+main(int argc, char **argv)
+{
+  constexpr unsigned int dim    = 2;
+  constexpr unsigned int degree = 1;
+  using Number                  = double;
+
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_rectangle(tria,
+                                            {2, 1},
+                                            {0., 0.0},
+                                            {2.0, 1.0});
+  if (tria.create_cell_iterator(CellId("0_0:"))->is_locally_owned())
+    tria.create_cell_iterator(CellId("0_0:"))->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FE_Q<dim>(degree));
+
+  LinearAlgebra::distributed::Vector<Number> v;
+  v.reinit(dof_handler.locally_owned_dofs(),
+           DoFTools::extract_locally_relevant_dofs(dof_handler),
+           MPI_COMM_WORLD);
+
+  std::map<CellId, Vector<double>> values;
+  values[CellId("1_0:")]  = Vector<double>{0.5, 0.5, 0.5, 0.5};
+  values[CellId("0_1:0")] = Vector<double>{0.5, 0.5, 0.5, 0.5};
+  values[CellId("0_1:1")] = Vector<double>{0.5, 0.5, 0.5, 0.1};
+  values[CellId("0_1:2")] = Vector<double>{0.5, 0.5, 0.5, 0.5};
+  values[CellId("0_1:3")] = Vector<double>{0.5, 0.1, 0.5, 0.5};
+
+  for (const auto cell : dof_handler.active_cell_iterators())
+    if (cell->is_artificial() == false)
+      cell->set_dof_values(values[cell->id()], v);
+
+  v.print(deallog.get_file_stream());
+
+  parallel::distributed::
+    SolutionTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
+      solution_trans(dof_handler, true /*enabling averaging*/);
+
+  v.update_ghost_values();
+  solution_trans.prepare_for_coarsening_and_refinement(v);
+
+  if (tria.create_cell_iterator(CellId("1_0:"))->is_locally_owned())
+    tria.create_cell_iterator(CellId("1_0:"))->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  dof_handler.distribute_dofs(FE_Q<dim>(degree));
+  v.reinit(dof_handler.locally_owned_dofs(),
+           DoFTools::extract_locally_relevant_dofs(dof_handler),
+           MPI_COMM_WORLD);
+  solution_trans.interpolate(v);
+
+  v.print(deallog.get_file_stream());
+}

tests/distributed_grids/solution_transfer_05.mpirun=1.with_p4est=true.output

@@ -0,0 +1,9 @@
+
+Process #0
+Local range: [0, 11), global size: 11
+Vector data:
+5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 1.000e-01 5.000e-01 5.000e-01 
+Process #0
+Local range: [0, 15), global size: 15
+Vector data:
+5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 3.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 

tests/distributed_grids/solution_transfer_05.mpirun=2.with_p4est=true.output

@@ -0,0 +1,19 @@
+
+Process #0
+Local range: [0, 9), global size: 11
+Vector data:
+5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 1.000e-01 5.000e-01 5.000e-01 5.000e-01 
+Process #0
+Local range: [0, 9), global size: 15
+Vector data:
+5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 3.000e-01 5.000e-01 5.000e-01 5.000e-01 
+
+Process #1
+Local range: [9, 11), global size: 11
+Vector data:
+5.000e-01 5.000e-01 
+Process #1
+Local range: [9, 15), global size: 15
+Vector data:
+5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 5.000e-01 
+