From a3f438a0489fd69f29acf3b482d7d8865162aa91 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 26 Oct 2017 17:51:09 -0600 Subject: [PATCH] Allow DataEntry functions to return either the real or imaginary part. --- include/deal.II/numerics/data_out_dof_data.h | 38 +++- .../numerics/data_out_dof_data.templates.h | 191 +++++++++++------- source/numerics/data_out.cc | 11 +- source/numerics/data_out_faces.cc | 11 +- source/numerics/data_out_rotation.cc | 11 +- 5 files changed, 181 insertions(+), 81 deletions(-) diff --git a/include/deal.II/numerics/data_out_dof_data.h b/include/deal.II/numerics/data_out_dof_data.h index 4de8a3a6a8..6b05d428ac 100644 --- a/include/deal.II/numerics/data_out_dof_data.h +++ b/include/deal.II/numerics/data_out_dof_data.h @@ -145,6 +145,35 @@ namespace internal { namespace DataOut { + /** + * The DataEntry classes abstract away the concrete data type of vectors + * users can attach to DataOut (and similar) objects and allow the underlying + * DataOut functions to query for individual elements of solution vectors + * without having to know the concrete vector type. This avoids that + * DataOut has to know what vectors are being used, but it has the downside + * that DataOut also doesn't know the underlying scalar type of these vectors. + * + * If the underlying scalar types all represent real numbers (in the + * mathematical sense -- i.e., the scalar type would be @p float, + * @p double, etc) then that is not a problem -- DataOut simply + * receives the values of individual vector components as @p double + * objects. On the other hand, if the vector type uses a std::complex + * scalar type, then DataEntry returning a @p double for a vector + * entry is not sufficient -- we need to provide DataOut with a way + * to query both the real and the imaginary part, so that they can + * be written into output files separately. + * + * This enum allows DataOut to tell a DataEntry function which component + * of a vector entry it wants to query, i.e., whether it wants the real + * or the imaginary part of a vector entry. + */ + enum class ComponentExtractor + { + real_part, + imaginary_part + }; + + /** * For each vector that has been added through the add_data_vector() * functions, we need to keep track of a pointer to it, and allow data @@ -190,7 +219,8 @@ namespace internal */ virtual double - get_cell_data_value (const unsigned int cell_number) const = 0; + get_cell_data_value (const unsigned int cell_number, + const ComponentExtractor extract_component) const = 0; /** * Given a FEValuesBase object, extract the values on the present cell @@ -200,6 +230,7 @@ namespace internal void get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector &patch_values) const = 0; /** @@ -211,6 +242,7 @@ namespace internal void get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_values_system) const = 0; /** @@ -221,6 +253,7 @@ namespace internal void get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_gradients) const = 0; /** @@ -232,6 +265,7 @@ namespace internal void get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_gradients_system) const = 0; /** @@ -242,6 +276,7 @@ namespace internal void get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_hessians) const = 0; /** @@ -253,6 +288,7 @@ namespace internal void get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_hessians_system) const = 0; /** 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 bf36a986ef..0777a9c96f 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -312,6 +313,66 @@ namespace internal { namespace DataOut { + /** + * Extract the specified component of a number. This template is used + * when the given value is assumed to be a real scalar, so asking + * for the real part is the only valid choice for the second argument. + */ + template + double get_component (const NumberType value, + const ComponentExtractor extract_component) + { + static_assert (numbers::NumberTraits::is_complex == false, + "This function must not be called for complex-valued data types."); + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + return value; + } + + + + /** + * Extract the specified component of a number. This template is used + * when the given value is a complex number + */ + template + double get_component (const std::complex &value, + const ComponentExtractor extract_component) + { + switch (extract_component) + { + case ComponentExtractor::real_part: + return value.real(); + + case ComponentExtractor::imaginary_part: + return value.imag(); + + default: + Assert (false, ExcInternalError()); + } + + return numbers::signaling_nan(); + } + + + + template + Tensor get_component (const Tensor &value, + const ComponentExtractor extract_component) + { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + + Tensor t; + for (unsigned int d=0; d DataEntryBase::DataEntryBase (const DoFHandlerType *dofs, @@ -407,7 +468,8 @@ namespace internal */ virtual double - get_cell_data_value (const unsigned int cell_number) const; + get_cell_data_value (const unsigned int cell_number, + const ComponentExtractor extract_component) const; /** * Given a FEValuesBase object, extract the values on the present cell @@ -417,6 +479,7 @@ namespace internal void get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector &patch_values) const; /** @@ -428,6 +491,7 @@ namespace internal void get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_values_system) const; /** @@ -438,6 +502,7 @@ namespace internal void get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_gradients) const; /** @@ -449,6 +514,7 @@ namespace internal void get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_gradients_system) const; /** @@ -459,6 +525,7 @@ namespace internal void get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_hessians) const; /** @@ -470,6 +537,7 @@ namespace internal void get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_hessians_system) const; /** @@ -526,7 +594,8 @@ namespace internal namespace { template - double + inline + typename VectorType::value_type get_vector_element (const VectorType &vector, const unsigned int cell_number) { @@ -534,7 +603,8 @@ namespace internal } - inline double + inline + double get_vector_element (const IndexSet &is, const unsigned int cell_number) { @@ -547,9 +617,11 @@ namespace internal template double DataEntry:: - get_cell_data_value (const unsigned int cell_number) const + get_cell_data_value (const unsigned int cell_number, + const ComponentExtractor extract_component) const { - return get_vector_element(*vector, cell_number); + return get_component (get_vector_element(*vector, cell_number), + extract_component); } @@ -558,21 +630,15 @@ namespace internal void DataEntry::get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_values_system) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_values (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -591,7 +657,8 @@ namespace internal fe_patch_values.get_function_values (*vector, tmp); for (unsigned int i = 0; i < patch_values_system.size(); i++) - patch_values_system[i] = tmp[i]; + for (unsigned int j=0; j::get_function_values (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector &patch_values) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_values (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -632,7 +693,7 @@ namespace internal fe_patch_values.get_function_values (*vector, tmp); for (unsigned int i = 0; i < tmp.size(); i++) - patch_values[i] = tmp[i]; + patch_values[i] = get_component (tmp[i], extract_component); } } @@ -642,21 +703,15 @@ namespace internal void DataEntry::get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_gradients_system) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_gradients (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -678,7 +733,7 @@ namespace internal for (unsigned int i = 0; i < tmp.size(); i++) for (unsigned int j = 0; j < tmp[i].size(); j++) - patch_gradients_system[i][j] = tmp[i][j]; + patch_gradients_system[i][j] = get_component (tmp[i][j], extract_component); } } @@ -688,21 +743,15 @@ namespace internal void DataEntry::get_function_gradients (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_gradients) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_gradients (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -720,7 +769,7 @@ namespace internal fe_patch_values.get_function_gradients (*vector, tmp); for (unsigned int i = 0; i < tmp.size(); i++) - patch_gradients[i] = tmp[i]; + patch_gradients[i] = get_component(tmp[i], extract_component); } } @@ -730,21 +779,15 @@ namespace internal void DataEntry::get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > > &patch_hessians_system) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_hessians (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -766,7 +809,7 @@ namespace internal for (unsigned int i = 0; i < tmp.size(); i++) for (unsigned int j = 0; j < tmp[i].size(); j++) - patch_hessians_system[i][j] = tmp[i][j]; + patch_hessians_system[i][j] = get_component (tmp[i][j], extract_component); } } @@ -776,21 +819,15 @@ namespace internal void DataEntry::get_function_hessians (const FEValuesBase &fe_patch_values, + const ComponentExtractor extract_component, std::vector > &patch_hessians) const { - // FIXME: FEValuesBase gives us data in types that match that of - // the solution vector. but this function needs to pass it back - // up as 'double' vectors. this requires the use of a temporary - // variable here if the data we get is not a 'double' vector. - // (of course, in reality, this also means that we may lose - // information to begin with.) - // - // the correct thing would be to also use the correct data type - // upstream somewhere, but this is complicated because we hide - // the actual data type from upstream. rather, we should at - // least make sure we can deal with complex numbers if (typeid(typename VectorType::value_type) == typeid(double)) { + Assert (extract_component == ComponentExtractor::real_part, + ExcMessage ("You cannot extract anything other than the real " + "part from a real number.")); + fe_patch_values.get_function_hessians (*vector, // reinterpret output argument type; because of // the 'if' statement above, this is the @@ -809,7 +846,7 @@ namespace internal fe_patch_values.get_function_hessians (*vector, tmp); for (unsigned int i = 0; i < tmp.size(); i++) - patch_hessians[i] = tmp[i]; + patch_hessians[i] = get_component(tmp[i], extract_component); } } diff --git a/source/numerics/data_out.cc b/source/numerics/data_out.cc index 04ff7b72cc..c1078d7c0f 100644 --- a/source/numerics/data_out.cc +++ b/source/numerics/data_out.cc @@ -156,12 +156,15 @@ build_one_patch // gradient etc. if (update_flags & update_values) this->dof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_scalar.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_scalar.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_scalar.solution_hessians); if (update_flags & update_quadrature_points) @@ -185,12 +188,15 @@ build_one_patch // derivative... if (update_flags & update_values) this->dof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_system.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_system.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_system.solution_hessians); if (update_flags & update_quadrature_points) @@ -221,6 +227,7 @@ build_one_patch if (n_components == 1) { this->dof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_scalar.solution_values); for (unsigned int q=0; qdof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, scratch_data.patch_values_system.solution_values); for (unsigned int component=0; componentcell_data.size(); ++dataset) { const double value - = this->cell_data[dataset]->get_cell_data_value (cell_and_index->second); + = this->cell_data[dataset]->get_cell_data_value (cell_and_index->second, + internal::DataOut::ComponentExtractor::real_part); for (unsigned int q=0; qdof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_hessians); if (update_flags & update_quadrature_points) @@ -186,12 +189,15 @@ build_one_patch (const FaceDescriptor *cell_and_face, data.resize_system_vectors(n_components); if (update_flags & update_values) this->dof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_hessians); if (update_flags & update_quadrature_points) @@ -224,6 +230,7 @@ build_one_patch (const FaceDescriptor *cell_and_face, if (n_components == 1) { this->dof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_values); for (unsigned int q=0; qdof_data[dataset]->get_function_values (this_fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_values); for (unsigned int component=0; component::active_cell_iterator(cell_and_face->first)); const double value - = this->cell_data[dataset]->get_cell_data_value (cell_number); + = this->cell_data[dataset]->get_cell_data_value (cell_number, + internal::DataOut::ComponentExtractor::real_part); for (unsigned int q=0; qdof_data[dataset]->get_function_values (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_hessians); if (update_flags & update_quadrature_points) @@ -233,12 +236,15 @@ build_one_patch (const cell_iterator // its derivative... if (update_flags & update_values) this->dof_data[dataset]->get_function_values (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_values); if (update_flags & update_gradients) this->dof_data[dataset]->get_function_gradients (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_gradients); if (update_flags & update_hessians) this->dof_data[dataset]->get_function_hessians (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_hessians); if (update_flags & update_quadrature_points) @@ -288,6 +294,7 @@ build_one_patch (const cell_iterator else if (n_components == 1) { this->dof_data[dataset]->get_function_values (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_scalar.solution_values); switch (dimension) @@ -320,6 +327,7 @@ build_one_patch (const cell_iterator { data.resize_system_vectors(n_components); this->dof_data[dataset]->get_function_values (fe_patch_values, + internal::DataOut::ComponentExtractor::real_part, data.patch_values_system.solution_values); for (unsigned int component=0; componenttriangulation->begin_active(), typename Triangulation::active_cell_iterator(*cell)); const double value - = this->cell_data[dataset]->get_cell_data_value (cell_number); + = this->cell_data[dataset]->get_cell_data_value (cell_number, + internal::DataOut::ComponentExtractor::real_part); switch (dimension) { case 1: -- 2.39.5