]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix the order in which we output complex-valued vectors/tensors.
authorWolfgang Bangerth <bangerth@colostate.edu>
Sat, 4 Apr 2020 00:12:28 +0000 (18:12 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 6 Apr 2020 21:52:57 +0000 (15:52 -0600)
include/deal.II/numerics/data_out_dof_data.templates.h
source/numerics/data_out.cc

index 1c29e0a86a41c6880a2a119d22dbaa499e143a0c..174ab684b5f2d8250ff7fda23b2fa1f203f0d411 100644 (file)
@@ -1690,6 +1690,9 @@ DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::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<DoFHandlerType, patch_dim, patch_space_dim>::
               // 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<DoFHandlerType, patch_dim, patch_space_dim>::
               // 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))
                 {
index a21ef902b70d04b0715029660a53411cb670e56c..676c6f45c104e12112813c1c42dfeec0e4aacb25 100644 (file)
@@ -685,7 +685,9 @@ DataOut<dim, DoFHandlerType>::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<dim, DoFHandlerType>::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<dim, DoFHandlerType>::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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.