Allow passing physical units into a VTU file.

diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h
index dab3f2df5aef1abced3e29a17fab234fb168a714..d4e3560ead160da056165498e10b364fddbdf670 100644 (file)
--- a/include/deal.II/base/data_out_base.h
+++ b/include/deal.II/base/data_out_base.h
@@ -1170,14 +1170,39 @@ namespace DataOutBase
     bool write_higher_order_cells;
 
     /**
-     * Constructor.
+     * A map that describes for (some or all) of the output quantities what
+     * the physical units are. This field is ignored for VTK file format, but
+     * used for VTU format where it is attached to the individual scalar,
+     * vector, or tensor fields for use by visualization or other postprocessing
+     * tools. The default is to not attach any physical units to fields at all,
+     * i.e., an empty map.
+     *
+     * If the map does not contain an entry for a specific output variable, then
+     * no unit will be written into the output file. In other words, it is not
+     * an error to provide units for only some variables.
+     *
+     * step-19, step-44 and step-69 all demonstrate how to use this variable.
+     *
+     * @note While the functions that make use of this information do not care
+     *   about how physical units are specified, downstream postprocessing tools
+     *   should and do. As a consequence, these units should be specified in a
+     *   format that is understandable to these postprocessing tools. As an
+     *   example, the [unyt project](https://unyt.readthedocs.io/en/stable/)
+     *   describes a standard for describing and converting units.
+     */
+    std::map<std::string, std::string> physical_units;
+
+    /**
+     * Constructor. Initializes the member variables with names corresponding
+     * to the argument names of this function.
      */
     VtkFlags(
       const double       time  = std::numeric_limits<double>::min(),
       const unsigned int cycle = std::numeric_limits<unsigned int>::min(),
       const bool         print_date_and_time              = true,
       const ZlibCompressionLevel compression_level        = best_compression,
-      const bool                 write_higher_order_cells = false);
+      const bool                 write_higher_order_cells = false,
+      const std::map<std::string, std::string> &physical_units = {});
   };
 
 
@@ -2077,7 +2102,8 @@ namespace DataOutBase
                  unsigned int,
                  std::string,
                  DataComponentInterpretation::DataComponentInterpretation>>
-      &nonscalar_data_ranges);
+      &             nonscalar_data_ranges,
+    const VtkFlags &flags);
 
   /**
    * In ParaView it is possible to visualize time-dependent data tagged with

diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index a4427526446f312a35e4cf63d9e1a11fae3ccaab..83fb060ed5c177c1937f542ff7e721b74220baef 100644 (file)
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -2579,12 +2579,14 @@ namespace DataOutBase
                      const unsigned int                   cycle,
                      const bool                           print_date_and_time,
                      const VtkFlags::ZlibCompressionLevel compression_level,
-                     const bool write_higher_order_cells)
+                     const bool write_higher_order_cells,
+                     const std::map<std::string, std::string> &physical_units)
     : time(time)
     , cycle(cycle)
     , print_date_and_time(print_date_and_time)
     , compression_level(compression_level)
     , write_higher_order_cells(write_higher_order_cells)
+    , physical_units(physical_units)
   {}
 
 
@@ -5735,7 +5737,24 @@ namespace DataOutBase
           }
 
         out << "\" NumberOfComponents=\"" << n_components << "\" format=\""
-            << ascii_or_binary << "\">\n";
+            << ascii_or_binary << "\"";
+        // If present, also list the physical units for this quantity. Look this
+        // up for either the name of the whole vector/tensor, or if that isn't
+        // listed, via its first component.
+        if (!name.empty())
+          {
+            if (flags.physical_units.find(name) != flags.physical_units.end())
+              out << " units=\"" << flags.physical_units.at(name) << "\"";
+          }
+        else
+          {
+            if (flags.physical_units.find(data_names[first_component]) !=
+                flags.physical_units.end())
+              out << " units=\""
+                  << flags.physical_units.at(data_names[first_component])
+                  << "\"";
+          }
+        out << ">\n";
 
         // now write data. pad all vectors to have three components
         std::vector<float> data;
@@ -5828,7 +5847,14 @@ namespace DataOutBase
         {
           out << "    <DataArray type=\"Float32\" Name=\""
               << data_names[data_set] << "\" format=\"" << ascii_or_binary
-              << "\">\n";
+              << "\"";
+          // If present, also list the physical units for this quantity.
+          if (flags.physical_units.find(data_names[data_set]) !=
+              flags.physical_units.end())
+            out << " units=\"" << flags.physical_units.at(data_names[data_set])
+                << "\"";
+
+          out << ">\n";
 
           std::vector<float> data(data_vectors[data_set].begin(),
                                   data_vectors[data_set].end());
@@ -5860,7 +5886,8 @@ namespace DataOutBase
                  unsigned int,
                  std::string,
                  DataComponentInterpretation::DataComponentInterpretation>>
-      &nonscalar_data_ranges)
+      &             nonscalar_data_ranges,
+    const VtkFlags &flags)
   {
     AssertThrow(out, ExcIO());
 
@@ -5917,8 +5944,9 @@ namespace DataOutBase
         // underscores unless a vector name has been specified
         out << "    <PDataArray type=\"Float32\" Name=\"";
 
-        if (!std::get<2>(nonscalar_data_range).empty())
-          out << std::get<2>(nonscalar_data_range);
+        const std::string name = std::get<2>(nonscalar_data_range);
+        if (!name.empty())
+          out << name;
         else
           {
             for (unsigned int i = std::get<0>(nonscalar_data_range);
@@ -5929,14 +5957,42 @@ namespace DataOutBase
           }
 
         out << "\" NumberOfComponents=\"" << n_components
-            << "\" format=\"ascii\"/>\n";
+            << "\" format=\"ascii\"";
+        // If present, also list the physical units for this quantity. Look this
+        // up for either the name of the whole vector/tensor, or if that isn't
+        // listed, via its first component.
+        if (!name.empty())
+          {
+            if (flags.physical_units.find(name) != flags.physical_units.end())
+              out << " units=\"" << flags.physical_units.at(name) << "\"";
+          }
+        else
+          {
+            if (flags.physical_units.find(
+                  data_names[std::get<1>(nonscalar_data_range)]) !=
+                flags.physical_units.end())
+              out << " units=\""
+                  << flags.physical_units.at(
+                       data_names[std::get<1>(nonscalar_data_range)])
+                  << "\"";
+          }
+
+        out << "/>\n";
       }
 
+    // Now for the scalar fields
     for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
       if (data_set_written[data_set] == false)
         {
           out << "    <PDataArray type=\"Float32\" Name=\""
-              << data_names[data_set] << "\" format=\"ascii\"/>\n";
+              << data_names[data_set] << "\" format=\"ascii\"";
+
+          if (flags.physical_units.find(data_names[data_set]) !=
+              flags.physical_units.end())
+            out << " units=\"" << flags.physical_units.at(data_names[data_set])
+                << "\"";
+
+          out << "/>\n";
         }
 
     out << "    </PPointData>\n";
@@ -7352,6 +7408,7 @@ DataOutInterface<dim, spacedim>::write_vtu_in_parallel(
 }
 
 
+
 template <int dim, int spacedim>
 void
 DataOutInterface<dim, spacedim>::write_pvtu_record(
@@ -7361,10 +7418,12 @@ DataOutInterface<dim, spacedim>::write_pvtu_record(
   DataOutBase::write_pvtu_record(out,
                                  piece_names,
                                  get_dataset_names(),
-                                 get_nonscalar_data_ranges());
+                                 get_nonscalar_data_ranges(),
+                                 vtk_flags);
 }
 
 
+
 template <int dim, int spacedim>
 std::string
 DataOutInterface<dim, spacedim>::write_vtu_with_pvtu_record(