Improve asserts in parallel::distributed::SolutionTransfer

doc/news/changes/minor/20220225Munch

@@ -0,0 +1,10 @@
+Improved: The asserts in parallel::distributed::SolutionTransfer have
+been improved. Before this change, an assert has been only called
+in the parallel case for distributed vectors for cells neighboring
+interprocess boundaries. Now, we check for all DoFs that the 
+contributions from all cells have the same values. This condition
+might be not given if the solution has high griadients at hanging
+nodes. For this reason, we furthermore added the possibility to 
+average the contributions to a DoF from different cells.
+<br>
+(Peter Munch, Magdalena Schreter, 2022/02/25)

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

@@ -1751,6 +1751,24 @@ class DoFCellAccessor : public DoFAccessor<dimension_,
   template <class OutputVector, typename number>
   void
   set_dof_values_by_interpolation(
+    const Vector<number> &local_values,
+    OutputVector &        values,
+    const unsigned int    fe_index =
+      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 =

include/deal.II/lac/la_parallel_block_vector.h

@@ -366,6 +366,12 @@ namespace LinearAlgebra
       bool
       has_ghost_elements() const;
 
+      /**
+       * Change the ghost state of this vector to @p ghosted.
+       */
+      void
+      set_ghost_state(const bool ghosted) const;
+
       /**
        * This method copies the data in the locally owned range from another
        * distributed vector @p src into the calling vector. As opposed to

include/deal.II/lac/la_parallel_block_vector.templates.h

@@ -421,6 +421,16 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number>
+    void
+    BlockVector<Number>::set_ghost_state(const bool ghosted) const
+    {
+      for (unsigned int block = 0; block < this->n_blocks(); ++block)
+        this->block(block).set_ghost_state(ghosted);
+    }
+
+
+
     template <typename Number>
     void
     BlockVector<Number>::reinit(const VectorSpaceVector<Number> &V,

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,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<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)
+            {
+              // 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<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,9 +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);
+        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