From c46002170930aafb3d78a5364c83a53ae5785e99 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 12 Mar 2021 21:41:30 +0100 Subject: [PATCH] DataOut: update ghost values Use new LinearAlgebra::ReadWriteVector::import Work on block vector --- doc/news/changes/major/20210609Munch | 6 + include/deal.II/numerics/data_out_dof_data.h | 16 +- .../numerics/data_out_dof_data.templates.h | 211 ++++++++++++++++-- 3 files changed, 197 insertions(+), 36 deletions(-) create mode 100644 doc/news/changes/major/20210609Munch diff --git a/doc/news/changes/major/20210609Munch b/doc/news/changes/major/20210609Munch new file mode 100644 index 0000000000..050a38bcb6 --- /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(). +
+(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 264294e24d..d01f76b976 100644 --- 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 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 d7db51a919..28406b9d95 100644 --- 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 #include +#include #include #include @@ -44,6 +45,8 @@ #include #include +#include +#include #include #include @@ -765,6 +768,152 @@ namespace internal + namespace + { + /** + * Copy the data from an arbitrary non-block vector to a + * LinearAlgebra::distributed::Vector. + */ + template + void + copy_locally_owned_data_from( + const VectorType &src, + LinearAlgebra::distributed::Vector + &dst) + { + LinearAlgebra::ReadWriteVector 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 ::value, + VectorType>::type * = nullptr> + void + create_dof_vector( + const DoFHandler &dof_handler, + const VectorType & src, + LinearAlgebra::distributed::BlockVector + &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 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 ::value, + VectorType>::type * = nullptr> + void + create_dof_vector( + const DoFHandler &dof_handler, + const VectorType & src, + LinearAlgebra::distributed::BlockVector + &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 ::value, + VectorType>::type * = nullptr> + void + create_cell_vector( + const VectorType &src, + LinearAlgebra::distributed::BlockVector + &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 ::value, + VectorType>::type * = nullptr> + void + create_cell_vector( + const VectorType &src, + LinearAlgebra::distributed::BlockVector + &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 DataEntry(const DoFHandler *dofs, const VectorType * data, const std::vector & 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 + vector; }; template + template DataEntry::DataEntry( const DoFHandler *dofs, const VectorType * data, const std::vector & names, const std::vector< DataComponentInterpretation::DataComponentInterpretation> - &data_component_interpretation) + & data_component_interpretation, + const DataVectorType actual_type) : DataEntryBase(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 *data_postprocessor) : DataEntryBase(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::get(*vector, cell_number), + internal::ElementAccess>::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 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> 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> 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::memory_consumption() const { - return (sizeof(vector) + + return (vector.memory_consumption() + MemoryConsumption::memory_consumption(this->names)); } @@ -1237,7 +1399,6 @@ namespace internal void DataEntry::clear() { - vector = nullptr; this->dof_handler = nullptr; } @@ -1694,7 +1855,11 @@ DataOut_DoFData:: // finally, add the data vector: auto new_entry = std::make_unique< internal::DataOutImplementation::DataEntry>( - 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)); -- 2.39.5