From 6e54f4fc00713438d3af00d3dc5aed9ba62ff493 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 3 Apr 2020 18:12:28 -0600 Subject: [PATCH] Fix the order in which we output complex-valued vectors/tensors. --- .../numerics/data_out_dof_data.templates.h | 9 + source/numerics/data_out.cc | 277 +++++++++++++++--- 2 files changed, 252 insertions(+), 34 deletions(-) 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 1c29e0a86a..174ab684b5 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -1690,6 +1690,9 @@ DataOut_DoFData::get_dataset_names() // the 2nd, 4th, 6th component for the imaginary part of the vector field. // Rather, we need to put all real components of the same vector field into // consecutive components. + // + // This sort of logic is also explained in some detail in + // DataOut::build_one_patch(). for (const auto &input_data : dof_data) if (input_data->is_complex_valued() == false || (input_data->postprocessor != nullptr)) @@ -1853,6 +1856,9 @@ DataOut_DoFData:: // complex-valued vector and it is not postprocessed, then we need // to do it twice -- once for the real parts and once for the // imaginary parts + // + // This sort of logic is also explained in some detail in + // DataOut::build_one_patch(). if (input_data->is_complex_valued() == false || (input_data->postprocessor != nullptr)) { @@ -1923,6 +1929,9 @@ DataOut_DoFData:: // complex-valued tensor and it is not postprocessed, then we need // to do it twice -- once for the real parts and once for the // imaginary parts + // + // This sort of logic is also explained in some detail in + // DataOut::build_one_patch(). if (input_data->is_complex_valued() == false || (input_data->postprocessor != nullptr)) { diff --git a/source/numerics/data_out.cc b/source/numerics/data_out.cc index a21ef902b7..676c6f45c1 100644 --- a/source/numerics/data_out.cc +++ b/source/numerics/data_out.cc @@ -685,7 +685,9 @@ DataOut::build_one_patch( // efficiency reasons. if (n_components == 1) { - // first output the real part of the solution vector + Assert(dataset->n_output_variables == 1, ExcInternalError()); + + // First output the real part of the solution vector dataset->get_function_values( this_fe_patch_values, internal::DataOutImplementation::ComponentExtractor:: @@ -694,8 +696,13 @@ DataOut::build_one_patch( for (unsigned int q = 0; q < n_q_points; ++q) patch.data(offset, q) = scratch_data.patch_values_scalar.solution_values[q]; + offset += 1; - // and if there is one, also output the imaginary part + // And if there is one, also output the imaginary part. Note + // that the problem is scalar-valued, so we can freely add the + // imaginary part after the real part without having to worry + // that we are interleaving the real components of a vector + // with the imaginary components of the same vector. if (dataset->is_complex_valued() == true) { dataset->get_function_values( @@ -704,59 +711,261 @@ DataOut::build_one_patch( imaginary_part, scratch_data.patch_values_scalar.solution_values); for (unsigned int q = 0; q < n_q_points; ++q) - patch.data(offset + 1, q) = + patch.data(offset, q) = scratch_data.patch_values_scalar.solution_values[q]; + offset += 1; } } else { scratch_data.resize_system_vectors(n_components); - // same as above: first the real part - const unsigned int stride = - (dataset->is_complex_valued() ? 2 : 1); - dataset->get_function_values( - this_fe_patch_values, - internal::DataOutImplementation::ComponentExtractor:: - real_part, - scratch_data.patch_values_system.solution_values); - for (unsigned int component = 0; component < n_components; - ++component) - for (unsigned int q = 0; q < n_q_points; ++q) - patch.data(offset + component * stride, q) = - scratch_data.patch_values_system.solution_values[q]( - component); - - // and if there is one, also output the imaginary part - if (dataset->is_complex_valued() == true) + // So we have a multi-component DoFHandler here. That's more + // complicated. If the vector is real-valued, then we can just + // get everything at all quadrature points and copy them into + // the output array. In fact, we don't have to worry at all + // about the interpretation of the components. + if (dataset->is_complex_valued() == false) { dataset->get_function_values( this_fe_patch_values, internal::DataOutImplementation::ComponentExtractor:: - imaginary_part, + real_part, scratch_data.patch_values_system.solution_values); for (unsigned int component = 0; component < n_components; ++component) for (unsigned int q = 0; q < n_q_points; ++q) - patch.data(offset + component * stride + 1, q) = + patch.data(offset + component, q) = scratch_data.patch_values_system.solution_values[q]( component); + + // Increment the counter for the actual data record. + offset += dataset->n_output_variables; + } + else + // The situation is more complicated if the input vector is + // complex-valued. The easiest approach would have been to + // just have all real and then all imaginary components. + // This would have been conceptually easy, but it has the + // annoying downside that if you have a vector-valued + // problem (say, [u v]) then the output order would have + // been [u_re, v_re, u_im, v_im]. That's tolerable, but not + // quite so nice because one typically thinks of real and + // imaginary parts as belonging together. We would really + // like the output order to be [u_re, u_im, v_re, v_im]. + // That, too, would have been easy to implement because one + // just has to interleave real and imaginary parts. + // + // But that's also not what we want. That's because if one + // were, for example, to solve a complex-valued Stokes + // problem (e.g., computing eigenfunctions of the Stokes + // operator), then one has solution components + // [[u v] p] and the proper output order is + // [[u_re v_re] [u_im v_im] p_re p_im]. + // In other words, the order in which we want to output + // data depends on the *interpretation* of components. + // + // Doing this requires a bit more code, and also needs to + // be in sync with what we do in + // DataOut_DoFData::get_dataset_names() and + // DataOut_DoFData::get_nonscalar_data_ranges(). + { + // Given this description, first get the real parts of + // all components: + dataset->get_function_values( + this_fe_patch_values, + internal::DataOutImplementation::ComponentExtractor:: + real_part, + scratch_data.patch_values_system.solution_values); + + // Then we need to distribute them to the correct + // location. This requires knowledge of the interpretation + // of components as discussed above. + { + Assert(dataset->data_component_interpretation.size() == + n_components, + ExcInternalError()); + + unsigned int destination = offset; + for (unsigned int component = 0; + component < n_components; + /* component is updated below */) + { + switch ( + dataset->data_component_interpretation[component]) + { + case DataComponentInterpretation:: + component_is_scalar: + { + // OK, a scalar component. Put all of the + // values into the current row + // ('destination'); then move 'component' + // forward by one (so we treat the next + // component) and 'destination' forward by + // two (because we're going to put the + // imaginary part of the current component + // into the next slot). + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination, q) = + scratch_data.patch_values_system + .solution_values[q](component); + + ++component; + destination += 2; + + break; + } + + case DataComponentInterpretation:: + component_is_part_of_vector: + { + // A vector component. Put the + // DoFHandlerType::space_dimension + // components into the next set of + // contiguous rows + // ('destination+c'); then move 'component' + // forward by spacedim (so we get to the + // next component after the current vector) + // and 'destination' forward by two*spacedim + // (because we're going to put the imaginary + // part of the vector into the subsequent + // spacedim slots). + const unsigned int size = + DoFHandlerType::space_dimension; + for (unsigned int c = 0; c < size; ++c) + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination + c, q) = + scratch_data.patch_values_system + .solution_values[q](component + c); + + component += size; + destination += 2 * size; + + break; + } + + case DataComponentInterpretation:: + component_is_part_of_tensor: + { + // Same approach as for vectors above. + const unsigned int size = + DoFHandlerType::space_dimension * + DoFHandlerType::space_dimension; + for (unsigned int c = 0; c < size; ++c) + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination + c, q) = + scratch_data.patch_values_system + .solution_values[q](component + c); + + component += size; + destination += 2 * size; + + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + } + + // And now we need to do the same thing again for the + // imaginary parts, starting at the top of the list of + // components/destinations again. + dataset->get_function_values( + this_fe_patch_values, + internal::DataOutImplementation::ComponentExtractor:: + imaginary_part, + scratch_data.patch_values_system.solution_values); + { + unsigned int destination = offset; + for (unsigned int component = 0; + component < n_components; + /* component is updated below */) + { + switch ( + dataset->data_component_interpretation[component]) + { + case DataComponentInterpretation:: + component_is_scalar: + { + // OK, a scalar component. Put all of the + // values into the row past the current one + // ('destination+1') since 'destination' is + // occupied by the real part. + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination + 1, q) = + scratch_data.patch_values_system + .solution_values[q](component); + + ++component; + destination += 2; + + break; + } + + case DataComponentInterpretation:: + component_is_part_of_vector: + { + // A vector component. Put the + // DoFHandlerType::space_dimension + // components into the set of contiguous + // rows that follow the real parts + // ('destination+spacedim+c'). + const unsigned int size = + DoFHandlerType::space_dimension; + for (unsigned int c = 0; c < size; ++c) + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination + size + c, q) = + scratch_data.patch_values_system + .solution_values[q](component + c); + + component += size; + destination += 2 * size; + + break; + } + + case DataComponentInterpretation:: + component_is_part_of_tensor: + { + // Same as for vectors. + const unsigned int size = + DoFHandlerType::space_dimension * + DoFHandlerType::space_dimension; + for (unsigned int c = 0; c < size; ++c) + for (unsigned int q = 0; q < n_q_points; + ++q) + patch.data(destination + size + c, q) = + scratch_data.patch_values_system + .solution_values[q](component + c); + + component += size; + destination += 2 * size; + + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + } + + // Increment the counter for the actual data record. We + // need to move it forward a number of positions equal to + // the number of components of this data set, times two + // because we dealt with a complex-valued input vector + offset += dataset->n_output_variables * 2; } } } - // Increment the counter for the actual data record. We need to - // move it forward a number of positions equal to the number - // of components of this data set; if the input consisted - // of a complex-valued quantity and if it is not further - // processed by a postprocessor, then we need two output - // slots for each input variable. - offset += dataset->n_output_variables * - (dataset->is_complex_valued() && - (dataset->postprocessor == nullptr) ? - 2 : - 1); - // Also update the dataset_number index that we carry along with the // for-loop over all data sets. ++dataset_number; -- 2.39.5