From 86d5a15dc35437c62d72730a2f542be745ce133d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 9 Jul 2021 12:37:20 -0600 Subject: [PATCH] Allow passing physical units into a VTU file. --- include/deal.II/base/data_out_base.h | 32 ++++++++++-- source/base/data_out_base.cc | 77 ++++++++++++++++++++++++---- 2 files changed, 97 insertions(+), 12 deletions(-) diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h index dab3f2df5a..d4e3560ead 100644 --- 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 physical_units; + + /** + * Constructor. Initializes the member variables with names corresponding + * to the argument names of this function. */ VtkFlags( const double time = std::numeric_limits::min(), const unsigned int cycle = std::numeric_limits::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 &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 a442752644..83fb060ed5 100644 --- 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 &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 data; @@ -5828,7 +5847,14 @@ namespace DataOutBase { out << " \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 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 << " (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 << " \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 << " \n"; @@ -7352,6 +7408,7 @@ DataOutInterface::write_vtu_in_parallel( } + template void DataOutInterface::write_pvtu_record( @@ -7361,10 +7418,12 @@ DataOutInterface::write_pvtu_record( DataOutBase::write_pvtu_record(out, piece_names, get_dataset_names(), - get_nonscalar_data_ranges()); + get_nonscalar_data_ranges(), + vtk_flags); } + template std::string DataOutInterface::write_vtu_with_pvtu_record( -- 2.39.5