From 8b9f8ff11cd529bcd173bf62298e8bce01ebe955 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 11 Jun 2021 08:07:51 +0200 Subject: [PATCH] Template DataEntry on value_type --- include/deal.II/base/numbers.h | 10 + .../deal.II/lac/read_write_vector.templates.h | 26 +-- .../numerics/data_out_dof_data.templates.h | 177 +++++++++--------- 3 files changed, 113 insertions(+), 100 deletions(-) diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h index 4e3432c7bc..32ccde7d29 100644 --- a/include/deal.II/base/numbers.h +++ b/include/deal.II/base/numbers.h @@ -428,6 +428,11 @@ namespace numbers */ using real_type = number; + /** + * For this data type, alias the corresponding double type. + */ + using double_type = double; + /** * Return the complex-conjugate of the given number. Since the general * template is selected if number is not a complex data type, this @@ -490,6 +495,11 @@ namespace numbers */ using real_type = number; + /** + * For this data type, alias the corresponding double type. + */ + using double_type = std::complex; + /** * Return the complex-conjugate of the given number. */ diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index ac4bc07f1a..e4ab936834 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -355,12 +355,13 @@ namespace LinearAlgebra if (n_elements() != in_vector.n_elements()) reinit(in_vector, true); - dealii::internal::VectorOperations::Vector_copy copier( - in_vector.values.get(), values.get()); - dealii::internal::VectorOperations::parallel_for(copier, - 0, - n_elements(), - thread_loop_partitioner); + if (n_elements() > 0) + { + dealii::internal::VectorOperations::Vector_copy copier( + in_vector.values.get(), values.get()); + dealii::internal::VectorOperations::parallel_for( + copier, 0, n_elements(), thread_loop_partitioner); + } return *this; } @@ -376,12 +377,13 @@ namespace LinearAlgebra if (n_elements() != in_vector.n_elements()) reinit(in_vector, true); - dealii::internal::VectorOperations::Vector_copy copier( - in_vector.values.get(), values.get()); - dealii::internal::VectorOperations::parallel_for(copier, - 0, - n_elements(), - thread_loop_partitioner); + if (n_elements() > 0) + { + dealii::internal::VectorOperations::Vector_copy copier( + in_vector.values.get(), values.get()); + dealii::internal::VectorOperations::parallel_for( + copier, 0, n_elements(), thread_loop_partitioner); + } return *this; } 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 0d386517c4..d2f5d89b49 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -774,17 +774,21 @@ namespace internal * Copy the data from an arbitrary non-block vector to a * LinearAlgebra::distributed::Vector. */ - template + template void copy_locally_owned_data_from( - const VectorType &src, - LinearAlgebra::distributed::Vector - &dst) + 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); + + LinearAlgebra::ReadWriteVector temp2; + temp2.reinit(temp, true); + temp2 = temp; + + dst.import(temp2, VectorOperation::insert); } /** @@ -793,14 +797,13 @@ namespace internal template ::value, VectorType>::type * = nullptr> void - create_dof_vector( - const DoFHandler &dof_handler, - const VectorType & src, - LinearAlgebra::distributed::BlockVector - &dst) + 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, @@ -840,14 +843,13 @@ namespace internal template ::value, VectorType>::type * = nullptr> void - create_dof_vector( - const DoFHandler &dof_handler, - const VectorType & src, - LinearAlgebra::distributed::BlockVector - &dst) + 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, @@ -869,13 +871,12 @@ namespace internal * 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) + create_cell_vector(const VectorType & src, + LinearAlgebra::distributed::BlockVector &dst) { dst.reinit(src.n_blocks()); @@ -893,13 +894,12 @@ namespace internal * 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) + create_cell_vector(const VectorType & src, + LinearAlgebra::distributed::BlockVector &dst) { dst.reinit(1); @@ -918,7 +918,7 @@ namespace internal * Class that stores a pointer to a vector of type equal to the template * argument, and provides the functions to extract data from it. */ - template + template class DataEntry : public DataEntryBase { public: @@ -927,7 +927,7 @@ namespace internal * the vector and their interpretation as scalar or vector data. This * constructor assumes that no postprocessor is going to be used. */ - template + template DataEntry(const DoFHandler *dofs, const VectorType * data, const std::vector & names, @@ -941,6 +941,7 @@ namespace internal * case, the names and vector declarations are going to be acquired from * the postprocessor. */ + template DataEntry(const DoFHandler * dofs, const VectorType * data, const DataPostprocessor *data_postprocessor); @@ -1040,15 +1041,14 @@ namespace internal * Pointer to the data vector. Note that ownership of the vector pointed * to remains with the caller of this class. */ - LinearAlgebra::distributed::BlockVector - vector; + LinearAlgebra::distributed::BlockVector vector; }; - template - template - DataEntry::DataEntry( + template + template + DataEntry::DataEntry( const DoFHandler *dofs, const VectorType * data, const std::vector & names, @@ -1068,8 +1068,9 @@ namespace internal - template - DataEntry::DataEntry( + template + template + DataEntry::DataEntry( const DoFHandler * dofs, const VectorType * data, const DataPostprocessor *data_postprocessor) @@ -1080,28 +1081,28 @@ namespace internal - template + template double - DataEntry::get_cell_data_value( + DataEntry::get_cell_data_value( const unsigned int cell_number, const ComponentExtractor extract_component) const { return get_component( internal::ElementAccess>::get(vector, cell_number), + ScalarType>>::get(vector, cell_number), extract_component); } - template + template void - DataEntry::get_function_values( + DataEntry::get_function_values( const FEValuesBase & fe_patch_values, const ComponentExtractor extract_component, std::vector> &patch_values_system) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1113,8 +1114,7 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast< - std::vector> &>( + reinterpret_cast> &>( patch_values_system)); } else @@ -1129,8 +1129,7 @@ namespace internal const unsigned int n_eval_points = fe_patch_values.n_quadrature_points; - std::vector> tmp( - n_eval_points); + std::vector> tmp(n_eval_points); for (unsigned int i = 0; i < n_eval_points; i++) tmp[i].reinit(n_components); @@ -1150,14 +1149,14 @@ namespace internal - template + template void - DataEntry::get_function_values( + DataEntry::get_function_values( const FEValuesBase &fe_patch_values, const ComponentExtractor extract_component, std::vector & patch_values) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1169,12 +1168,11 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast &>( - patch_values)); + reinterpret_cast &>(patch_values)); } else { - std::vector tmp(patch_values.size()); + std::vector tmp(patch_values.size()); fe_patch_values.get_function_values(vector, tmp); @@ -1185,15 +1183,15 @@ namespace internal - template + template void - DataEntry::get_function_gradients( + DataEntry::get_function_gradients( const FEValuesBase & fe_patch_values, const ComponentExtractor extract_component, std::vector>> &patch_gradients_system) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1205,8 +1203,8 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast>> &>( + reinterpret_cast< + std::vector>> &>( patch_gradients_system)); } else @@ -1221,9 +1219,8 @@ namespace internal const unsigned int n_eval_points = fe_patch_values.n_quadrature_points; - std::vector< - std::vector>> - tmp(n_eval_points); + std::vector>> tmp( + n_eval_points); for (unsigned int i = 0; i < n_eval_points; i++) tmp[i].resize(n_components); @@ -1243,14 +1240,14 @@ namespace internal - template + template void - DataEntry::get_function_gradients( + DataEntry::get_function_gradients( const FEValuesBase &fe_patch_values, const ComponentExtractor extract_component, std::vector> & patch_gradients) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1262,13 +1259,12 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast> &>( + reinterpret_cast> &>( patch_gradients)); } else { - std::vector> tmp; + std::vector> tmp; tmp.resize(patch_gradients.size()); fe_patch_values.get_function_gradients(vector, tmp); @@ -1280,15 +1276,15 @@ namespace internal - template + template void - DataEntry::get_function_hessians( + DataEntry::get_function_hessians( const FEValuesBase & fe_patch_values, const ComponentExtractor extract_component, std::vector>> &patch_hessians_system) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1300,8 +1296,8 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast>> &>( + reinterpret_cast< + std::vector>> &>( patch_hessians_system)); } else @@ -1316,9 +1312,8 @@ namespace internal const unsigned int n_eval_points = fe_patch_values.n_quadrature_points; - std::vector< - std::vector>> - tmp(n_eval_points); + std::vector>> tmp( + n_eval_points); for (unsigned int i = 0; i < n_eval_points; i++) tmp[i].resize(n_components); @@ -1338,14 +1333,14 @@ namespace internal - template + template void - DataEntry::get_function_hessians( + DataEntry::get_function_hessians( const FEValuesBase &fe_patch_values, const ComponentExtractor extract_component, std::vector> & patch_hessians) const { - if (typeid(typename VectorType::value_type) == typeid(double)) + if (typeid(ScalarType) == typeid(double)) { Assert(extract_component == ComponentExtractor::real_part, ExcMessage("You cannot extract anything other than the real " @@ -1357,13 +1352,12 @@ namespace internal // above, this is the identity cast whenever the code is executed, // but the cast is necessary to allow compilation even if we don't // get here - reinterpret_cast> &>( + reinterpret_cast> &>( patch_hessians)); } else { - std::vector> tmp( + std::vector> tmp( patch_hessians.size()); fe_patch_values.get_function_hessians(vector, tmp); @@ -1375,18 +1369,18 @@ namespace internal - template + template bool - DataEntry::is_complex_valued() const + DataEntry::is_complex_valued() const { - return numbers::NumberTraits::is_complex; + return numbers::NumberTraits::is_complex; } - template + template std::size_t - DataEntry::memory_consumption() const + DataEntry::memory_consumption() const { return (vector.memory_consumption() + MemoryConsumption::memory_consumption(this->names)); @@ -1394,9 +1388,9 @@ namespace internal - template + template void - DataEntry::clear() + DataEntry::clear() { this->dof_handler = nullptr; } @@ -1732,9 +1726,13 @@ DataOut_DoFData::add_data_vector( dof_handler.get_triangulation().n_active_cells())); - auto new_entry = std::make_unique< - internal::DataOutImplementation::DataEntry>( - &dof_handler, &vec, &data_postprocessor); + auto new_entry = std::make_unique::double_type>>(&dof_handler, + &vec, + &data_postprocessor); dof_data.emplace_back(std::move(new_entry)); } @@ -1852,8 +1850,11 @@ DataOut_DoFData:: DataComponentInterpretation::component_is_scalar)); // finally, add the data vector: - auto new_entry = std::make_unique< - internal::DataOutImplementation::DataEntry>( + auto new_entry = std::make_unique::double_type>>( dof_handler, &data_vector, deduced_names, -- 2.39.5