DataOut: update ghost values

Use new LinearAlgebra::ReadWriteVector::import

Work on block vector

diff --git a/doc/news/changes/major/20210609Munch b/doc/news/changes/major/20210609Munch
new file mode 100644 (file)
index 0000000..050a38b
--- /dev/null
+++ b/doc/news/changes/major/20210609Munch
@@ -0,0 +1,6 @@
+Improved: DataOut_DoFData and the derived classes create internally
+a copy of vectors attached via add_data_vector() and are now 
+responsible for updating ghost values so that users do not
+need to call update_ghost_values() and potentially zero_out_ghost_values().
+<br>
+(Peter Munch, Magdalena Schreter, 2021/05/27)

diff --git a/include/deal.II/numerics/data_out_dof_data.h b/include/deal.II/numerics/data_out_dof_data.h
index 264294e24dbc0dd16bb8dd3cefbdfd6c17751e3c..d01f76b976b28254bc078878546a863b966037c7 100644 (file)
--- a/include/deal.II/numerics/data_out_dof_data.h
+++ b/include/deal.II/numerics/data_out_dof_data.h
@@ -522,11 +522,6 @@ namespace internal
  * @ref step_22 "step-22"
  * tutorial program).
  *
- * This class does not copy the vector given to it through the
- * add_data_vector() functions, for memory consumption reasons. It only stores
- * a reference to it, so it is in your responsibility to make sure that the
- * data vectors exist long enough.
- *
  * After adding all data vectors, you need to call a function which generates
  * the patches (i.e., some intermediate data representation) for output from
  * the stored data. Derived classes name this function build_patches().
@@ -726,11 +721,6 @@ public:
    * @note The actual type for the vector argument may be any vector type from
    * which FEValues can extract values on a cell using the
    * FEValuesBase::get_function_values() function.
-   *
-   * @note When working in parallel, the vector to be written needs to be ghosted
-   * with read access to all degrees of freedom on the locally owned cells, see
-   * the step-40 or step-37 tutorial programs for details, i.e., it might be
-   * necessary to call data.update_ghost_values().
    */
   template <class VectorType>
   void
@@ -897,13 +887,13 @@ public:
   clear_data_vectors();
 
   /**
-   * Release pointers to all input data elements, i.e. pointers to data
-   * vectors and to the DoF handler object. This function may be useful when
+   * Release pointers to all input data elements, i.e. pointers to
+   * to the DoF handler object. This function may be useful when
    * you have called the @p build_patches function of derived class, since
    * then the patches are built and the input data is no more needed, nor is
    * there a need to reference it. You can then output the patches detached
    * from the main thread and need not make sure anymore that the DoF handler
-   * object and vectors must not be deleted before the output thread is
+   * object must not be deleted before the output thread is
    * finished.
    */
   void

diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h
index d7db51a91937704ad633963509f3188f95e44357..28406b9d957a807a4e3554156d629d8abc8a8f4e 100644 (file)
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -28,6 +28,7 @@
 
 #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_dgq.h>
 #include <deal.II/fe/fe_pyramid_p.h>
@@ -44,6 +45,8 @@
 #include <deal.II/hp/fe_values.h>
 #include <deal.II/hp/q_collection.h>
 
+#include <deal.II/lac/block_vector_base.h>
+#include <deal.II/lac/read_write_vector.h>
 #include <deal.II/lac/vector.h>
 
 #include <deal.II/numerics/data_out.h>
@@ -765,6 +768,152 @@ namespace internal
 
 
 
+    namespace
+    {
+      /**
+       * Copy the data from an arbitrary non-block vector to a
+       * LinearAlgebra::distributed::Vector.
+       */
+      template <typename VectorType>
+      void
+      copy_locally_owned_data_from(
+        const VectorType &src,
+        LinearAlgebra::distributed::Vector<typename VectorType::value_type>
+          &dst)
+      {
+        LinearAlgebra::ReadWriteVector<typename VectorType::value_type> temp;
+        temp.reinit(src.locally_owned_elements());
+        temp.import(src, VectorOperation::insert);
+        dst.import(temp, VectorOperation::insert);
+      }
+
+      /**
+       * Create a ghosted-copy of a block dof vector.
+       */
+      template <int dim,
+                int spacedim,
+                typename VectorType,
+                typename std::enable_if<IsBlockVector<VectorType>::value,
+                                        VectorType>::type * = nullptr>
+      void
+      create_dof_vector(
+        const DoFHandler<dim, spacedim> &dof_handler,
+        const VectorType &               src,
+        LinearAlgebra::distributed::BlockVector<typename VectorType::value_type>
+          &dst)
+      {
+        IndexSet locally_relevant_dofs;
+        DoFTools::extract_locally_relevant_dofs(dof_handler,
+                                                locally_relevant_dofs);
+
+        const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
+
+        std::vector<types::global_dof_index> n_indices_per_block(
+          src.n_blocks());
+
+        for (unsigned int b = 0; b < src.n_blocks(); ++b)
+          n_indices_per_block[b] = src.get_block_indices().block_size(b);
+
+        const auto locally_owned_dofs_b =
+          locally_owned_dofs.split_by_block(n_indices_per_block);
+        const auto locally_relevant_dofs_b =
+          locally_relevant_dofs.split_by_block(n_indices_per_block);
+
+        dst.reinit(src.n_blocks());
+
+        for (unsigned int b = 0; b < src.n_blocks(); ++b)
+          {
+            dst.block(b).reinit(locally_owned_dofs_b[b],
+                                locally_relevant_dofs_b[b],
+                                dof_handler.get_communicator());
+            copy_locally_owned_data_from(src.block(b), dst.block(b));
+          }
+
+        dst.collect_sizes();
+
+        dst.update_ghost_values();
+      }
+
+      /**
+       * Create a ghosted-copy of a non-block dof vector.
+       */
+      template <int dim,
+                int spacedim,
+                typename VectorType,
+                typename std::enable_if<!IsBlockVector<VectorType>::value,
+                                        VectorType>::type * = nullptr>
+      void
+      create_dof_vector(
+        const DoFHandler<dim, spacedim> &dof_handler,
+        const VectorType &               src,
+        LinearAlgebra::distributed::BlockVector<typename VectorType::value_type>
+          &dst)
+      {
+        IndexSet locally_relevant_dofs;
+        DoFTools::extract_locally_relevant_dofs(dof_handler,
+                                                locally_relevant_dofs);
+
+        dst.reinit(1);
+
+        dst.block(0).reinit(dof_handler.locally_owned_dofs(),
+                            locally_relevant_dofs,
+                            dof_handler.get_communicator());
+        copy_locally_owned_data_from(src, dst.block(0));
+
+        dst.collect_sizes();
+
+        dst.update_ghost_values();
+      }
+
+      /**
+       * Create a ghosted-copy of a block cell vector.
+       */
+      template <typename VectorType,
+                typename std::enable_if<IsBlockVector<VectorType>::value,
+                                        VectorType>::type * = nullptr>
+      void
+      create_cell_vector(
+        const VectorType &src,
+        LinearAlgebra::distributed::BlockVector<typename VectorType::value_type>
+          &dst)
+      {
+        dst.reinit(src.n_blocks());
+
+        for (unsigned int b = 0; b < src.n_blocks(); ++b)
+          {
+            dst.block(b).reinit(src.get_block_indices().block_size(b));
+            copy_locally_owned_data_from(src.block(b), dst.block(b));
+          }
+
+        dst.collect_sizes();
+      }
+
+
+      /**
+       * Create a ghosted-copy of a non-block cell vector.
+       */
+      template <typename VectorType,
+                typename std::enable_if<!IsBlockVector<VectorType>::value,
+                                        VectorType>::type * = nullptr>
+      void
+      create_cell_vector(
+        const VectorType &src,
+        LinearAlgebra::distributed::BlockVector<typename VectorType::value_type>
+          &dst)
+      {
+        dst.reinit(1);
+
+        dst.block(0).reinit(src.size());
+        copy_locally_owned_data_from(src, dst.block(0));
+
+        dst.collect_sizes();
+
+        dst.update_ghost_values();
+      }
+    } // namespace
+
+
+
     /**
      * Class that stores a pointer to a vector of type equal to the template
      * argument, and provides the functions to extract data from it.
@@ -778,12 +927,14 @@ namespace internal
        * the vector and their interpretation as scalar or vector data. This
        * constructor assumes that no postprocessor is going to be used.
        */
+      template <typename DataVectorType>
       DataEntry(const DoFHandler<dim, spacedim> *dofs,
                 const VectorType *               data,
                 const std::vector<std::string> & names,
                 const std::vector<
                   DataComponentInterpretation::DataComponentInterpretation>
-                  &data_component_interpretation);
+                  &                  data_component_interpretation,
+                const DataVectorType actual_type);
 
       /**
        * Constructor when a data postprocessor is going to be used. In that
@@ -890,22 +1041,31 @@ namespace internal
        * Pointer to the data vector. Note that ownership of the vector pointed
        * to remains with the caller of this class.
        */
-      const VectorType *vector;
+      LinearAlgebra::distributed::BlockVector<typename VectorType::value_type>
+        vector;
     };
 
 
 
     template <int dim, int spacedim, typename VectorType>
+    template <typename DataVectorType>
     DataEntry<dim, spacedim, VectorType>::DataEntry(
       const DoFHandler<dim, spacedim> *dofs,
       const VectorType *               data,
       const std::vector<std::string> & names,
       const std::vector<
         DataComponentInterpretation::DataComponentInterpretation>
-        &data_component_interpretation)
+        &                  data_component_interpretation,
+      const DataVectorType actual_type)
       : DataEntryBase<dim, spacedim>(dofs, names, data_component_interpretation)
-      , vector(data)
-    {}
+    {
+      if (actual_type == DataVectorType::type_dof_data)
+        create_dof_vector(*dofs, *data, vector);
+      else if (actual_type == DataVectorType::type_cell_data)
+        create_cell_vector(*data, vector);
+      else
+        Assert(false, ExcInternalError());
+    }
 
 
 
@@ -915,8 +1075,9 @@ namespace internal
       const VectorType *                 data,
       const DataPostprocessor<spacedim> *data_postprocessor)
       : DataEntryBase<dim, spacedim>(dofs, data_postprocessor)
-      , vector(data)
-    {}
+    {
+      create_dof_vector(*dofs, *data, vector);
+    }
 
 
 
@@ -927,7 +1088,8 @@ namespace internal
       const ComponentExtractor extract_component) const
     {
       return get_component(
-        internal::ElementAccess<VectorType>::get(*vector, cell_number),
+        internal::ElementAccess<LinearAlgebra::distributed::BlockVector<
+          typename VectorType::value_type>>::get(vector, cell_number),
         extract_component);
     }
 
@@ -947,7 +1109,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_values(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -973,7 +1135,7 @@ namespace internal
           for (unsigned int i = 0; i < n_eval_points; i++)
             tmp[i].reinit(n_components);
 
-          fe_patch_values.get_function_values(*vector, tmp);
+          fe_patch_values.get_function_values(vector, tmp);
 
           AssertDimension(patch_values_system.size(), n_eval_points);
           for (unsigned int i = 0; i < n_eval_points; i++)
@@ -1003,7 +1165,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_values(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -1015,7 +1177,7 @@ namespace internal
         {
           std::vector<typename VectorType::value_type> tmp(patch_values.size());
 
-          fe_patch_values.get_function_values(*vector, tmp);
+          fe_patch_values.get_function_values(vector, tmp);
 
           for (unsigned int i = 0; i < tmp.size(); i++)
             patch_values[i] = get_component(tmp[i], extract_component);
@@ -1039,7 +1201,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_gradients(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -1066,7 +1228,7 @@ namespace internal
           for (unsigned int i = 0; i < n_eval_points; i++)
             tmp[i].resize(n_components);
 
-          fe_patch_values.get_function_gradients(*vector, tmp);
+          fe_patch_values.get_function_gradients(vector, tmp);
 
           AssertDimension(patch_gradients_system.size(), n_eval_points);
           for (unsigned int i = 0; i < n_eval_points; i++)
@@ -1096,7 +1258,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_gradients(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -1110,7 +1272,7 @@ namespace internal
           std::vector<Tensor<1, spacedim, typename VectorType::value_type>> tmp;
           tmp.resize(patch_gradients.size());
 
-          fe_patch_values.get_function_gradients(*vector, tmp);
+          fe_patch_values.get_function_gradients(vector, tmp);
 
           for (unsigned int i = 0; i < tmp.size(); i++)
             patch_gradients[i] = get_component(tmp[i], extract_component);
@@ -1134,7 +1296,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_hessians(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -1161,7 +1323,7 @@ namespace internal
           for (unsigned int i = 0; i < n_eval_points; i++)
             tmp[i].resize(n_components);
 
-          fe_patch_values.get_function_hessians(*vector, tmp);
+          fe_patch_values.get_function_hessians(vector, tmp);
 
           AssertDimension(patch_hessians_system.size(), n_eval_points);
           for (unsigned int i = 0; i < n_eval_points; i++)
@@ -1191,7 +1353,7 @@ namespace internal
                             "part from a real number."));
 
           fe_patch_values.get_function_hessians(
-            *vector,
+            vector,
             // reinterpret output argument type; because of the 'if' statement
             // above, this is the identity cast whenever the code is executed,
             // but the cast is necessary to allow compilation even if we don't
@@ -1205,7 +1367,7 @@ namespace internal
           std::vector<Tensor<2, spacedim, typename VectorType::value_type>> tmp(
             patch_hessians.size());
 
-          fe_patch_values.get_function_hessians(*vector, tmp);
+          fe_patch_values.get_function_hessians(vector, tmp);
 
           for (unsigned int i = 0; i < tmp.size(); i++)
             patch_hessians[i] = get_component(tmp[i], extract_component);
@@ -1227,7 +1389,7 @@ namespace internal
     std::size_t
     DataEntry<dim, spacedim, VectorType>::memory_consumption() const
     {
-      return (sizeof(vector) +
+      return (vector.memory_consumption() +
               MemoryConsumption::memory_consumption(this->names));
     }
 
@@ -1237,7 +1399,6 @@ namespace internal
     void
     DataEntry<dim, spacedim, VectorType>::clear()
     {
-      vector            = nullptr;
       this->dof_handler = nullptr;
     }
 
@@ -1694,7 +1855,11 @@ DataOut_DoFData<dim, patch_dim, spacedim, patch_spacedim>::
   // finally, add the data vector:
   auto new_entry = std::make_unique<
     internal::DataOutImplementation::DataEntry<dim, spacedim, VectorType>>(
-    dof_handler, &data_vector, deduced_names, data_component_interpretation);
+    dof_handler,
+    &data_vector,
+    deduced_names,
+    data_component_interpretation,
+    actual_type);
 
   if (actual_type == type_dof_data)
     dof_data.emplace_back(std::move(new_entry));