]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Averaging in parallel::distributed::SolutionTransfer
authorPeter Munch <peterrmuench@gmail.com>
Fri, 25 Feb 2022 18:18:04 +0000 (19:18 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 14 May 2022 11:38:08 +0000 (13:38 +0200)
include/deal.II/distributed/solution_transfer.h
include/deal.II/dofs/dof_accessor.h
source/distributed/solution_transfer.cc
source/dofs/dof_accessor_set.cc
source/dofs/dof_accessor_set.inst.in

index a2c0e21059c4c70bdcf010c44d5f55cfdfb7b25f..d717c90985d442ce65a852e9dbf5cc1d97c624e1 100644 (file)
@@ -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);
 
 
       /**
index 2f561c44ed455574c7987130a6769feb52e01721..3f7285bbac860a5dc16b94b8f47260a275fe23ac 100644 (file)
@@ -1766,6 +1766,23 @@ public:
       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
index 889658aad0ce08e4ebf412203d108f802d09dcf7..f9c476c53060cd017534457630292246b461c67b 100644 (file)
@@ -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,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<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
index bd2f50e5c68fe30648f8c424305a7dfbb0c48ec0..20dd8d96855622109817837f6d3ae1d38f53e47f 100644 (file)
@@ -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"
index 6bdc91315faf0d906f9f65b4f4d81231033caff2..fabddab2a7dbd7de1d5c7c3757a91b90a5c945cf 100644 (file)
@@ -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
   }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.