From cd14f7c0ad675ae198a8bfcda916fae81bed46bb Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 11 Nov 2022 13:13:19 -0700 Subject: [PATCH] Encapulate different parts of write_vtu_main() into lambda functions. --- source/base/data_out_base.cc | 921 +++++++++++++++++++---------------- 1 file changed, 498 insertions(+), 423 deletions(-) diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index b4c8b1aeab..f0c05d9484 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -6096,159 +6096,220 @@ namespace DataOutBase return create_global_data_table(patches); }); - //----------------------------- - // first make up a list of used vertices along with their coordinates - // - // note that according to the standard, we have to print d=1..3 dimensions, - // even if we are in reality in 2d, for example out << "\n"; - out << " \n"; - out << " \n"; + o << " \n"; - const std::vector> node_positions = - get_node_positions(patches); - std::vector node_coordinates_3d; - node_coordinates_3d.reserve(node_positions.size() * 3); - for (const auto &node_position : node_positions) - { - node_coordinates_3d.emplace_back(node_position[0]); + const std::vector> node_positions = + get_node_positions(patches); + + // note that according to the standard, we have to print d=1..3 + // dimensions, even if we are in reality in 2d, for example + std::vector node_coordinates_3d; + node_coordinates_3d.reserve(node_positions.size() * 3); + for (const auto &node_position : node_positions) + { + node_coordinates_3d.emplace_back(node_position[0]); - if (spacedim >= 2) - node_coordinates_3d.emplace_back(node_position[1]); - else - node_coordinates_3d.emplace_back(0.0f); + if (spacedim >= 2) + node_coordinates_3d.emplace_back(node_position[1]); + else + node_coordinates_3d.emplace_back(0.0f); - if (spacedim >= 3) - node_coordinates_3d.emplace_back(node_position[2]); - else - node_coordinates_3d.emplace_back(0.0f); - } - out << vtu_stringize_array(node_coordinates_3d, + if (spacedim >= 3) + node_coordinates_3d.emplace_back(node_position[2]); + else + node_coordinates_3d.emplace_back(0.0f); + } + o << vtu_stringize_array(node_coordinates_3d, flags.compression_level, - out.precision()) + output_precision) << '\n'; - out << " \n"; - out << " \n\n"; + o << " \n"; + o << " \n\n"; + + return o.str(); + }; + out << stringize_vertex_information(); + //------------------------------- - // now for the cells - out << " \n"; - out << " \n"; + o << " \n"; - if (flags.write_higher_order_cells) - { - std::ostringstream o; + if (flags.write_higher_order_cells) { - VtuStream vtu_out(o, flags); + std::ostringstream oo; + { + VtuStream vtu_out(oo, flags); - write_high_order_cells(patches, vtu_out, /* legacy_format = */ false); - vtu_out.flush_cells(); + write_high_order_cells(patches, + vtu_out, + /* legacy_format = */ false); + vtu_out.flush_cells(); + } + o << oo.str() << '\n'; } - out << o.str() << '\n'; - } - else - { - Assert(dim <= 3, ExcNotImplemented()); + else + { + Assert(dim <= 3, ExcNotImplemented()); - std::vector cells; - unsigned int first_vertex_of_patch = 0; + std::vector cells; + unsigned int first_vertex_of_patch = 0; - for (const auto &patch : patches) - { - // special treatment of simplices since they are not subdivided - if (patch.reference_cell != ReferenceCells::get_hypercube()) - { - const unsigned int n_points = patch.data.n_cols(); - static const std::array - pyramid_index_translation_table = {{0, 1, 3, 2, 4}}; + for (const auto &patch : patches) + { + // special treatment of simplices since they are not subdivided + if (patch.reference_cell != ReferenceCells::get_hypercube()) + { + const unsigned int n_points = patch.data.n_cols(); + static const std::array + pyramid_index_translation_table = {{0, 1, 3, 2, 4}}; - if (deal_ii_with_zlib && - (flags.compression_level != - DataOutBase::CompressionLevel::plain_text)) - { - for (unsigned int i = 0; i < n_points; ++i) - cells.push_back( - first_vertex_of_patch + - (patch.reference_cell == ReferenceCells::Pyramid ? - pyramid_index_translation_table[i] : - i)); - } - else - { - for (unsigned int i = 0; i < n_points; ++i) - out << '\t' + if (deal_ii_with_zlib && + (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + for (unsigned int i = 0; i < n_points; ++i) + cells.push_back( + first_vertex_of_patch + + (patch.reference_cell == ReferenceCells::Pyramid ? + pyramid_index_translation_table[i] : + i)); + } + else + { + for (unsigned int i = 0; i < n_points; ++i) + o << '\t' << first_vertex_of_patch + (patch.reference_cell == ReferenceCells::Pyramid ? pyramid_index_translation_table[i] : i); - out << '\n'; - } + o << '\n'; + } - first_vertex_of_patch += n_points; - } - else - { - const unsigned int n_subdivisions = patch.n_subdivisions; - const unsigned int n_points_per_direction = n_subdivisions + 1; + first_vertex_of_patch += n_points; + } + else + { + const unsigned int n_subdivisions = patch.n_subdivisions; + const unsigned int n_points_per_direction = + n_subdivisions + 1; - switch (dim) - { - case 0: - { - auto write_cell = - [&flags, &out, &cells](const unsigned int start) { - if (deal_ii_with_zlib && - (flags.compression_level != - DataOutBase::CompressionLevel::plain_text)) - { - cells.push_back(start); - } - else - { - out << start; - out << '\n'; - } - }; + switch (dim) + { + case 0: + { + auto write_cell = + [&flags, &o, &cells](const unsigned int start) { + if (deal_ii_with_zlib && + (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + cells.push_back(start); + } + else + { + o << start; + o << '\n'; + } + }; + + const unsigned int starting_offset = + first_vertex_of_patch; + write_cell(starting_offset); + break; + } - const unsigned int starting_offset = - first_vertex_of_patch; - write_cell(starting_offset); - break; - } + case 1: + { + auto write_cell = + [&flags, &o, &cells](const unsigned int start) { + if (deal_ii_with_zlib && + (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + cells.push_back(start); + cells.push_back(start + 1); + } + else + { + o << start << '\t' << start + 1; + o << '\n'; + } + }; - case 1: - { - auto write_cell = - [&flags, &out, &cells](const unsigned int start) { - if (deal_ii_with_zlib && - (flags.compression_level != - DataOutBase::CompressionLevel::plain_text)) - { - cells.push_back(start); - cells.push_back(start + 1); - } - else + for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1) + { + const unsigned int starting_offset = + first_vertex_of_patch + i1; + write_cell(starting_offset); + } + break; + } + + case 2: + { + auto write_cell = + [&flags, &o, &cells, n_points_per_direction]( + const unsigned int start) { + if (deal_ii_with_zlib && + (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + cells.push_back(start); + cells.push_back(start + 1); + cells.push_back(start + + n_points_per_direction + 1); + cells.push_back(start + + n_points_per_direction); + } + else + { + o << start << '\t' << start + 1 << '\t' + << start + n_points_per_direction + 1 + << '\t' << start + n_points_per_direction; + o << '\n'; + } + }; + + for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2) + for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1) { - out << start << '\t' << start + 1; - out << '\n'; + const unsigned int starting_offset = + first_vertex_of_patch + + i2 * n_points_per_direction + i1; + write_cell(starting_offset); } - }; - - for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1) - { - const unsigned int starting_offset = - first_vertex_of_patch + i1; - write_cell(starting_offset); - } - break; - } + break; + } - case 2: - { - auto write_cell = - [&flags, &out, &cells, n_points_per_direction]( - const unsigned int start) { + case 3: + { + auto write_cell = [&flags, + &o, + &cells, + n_points_per_direction]( + const unsigned int start) { if (deal_ii_with_zlib && (flags.compression_level != DataOutBase::CompressionLevel::plain_text)) @@ -6258,61 +6319,25 @@ namespace DataOutBase cells.push_back(start + n_points_per_direction + 1); cells.push_back(start + n_points_per_direction); + cells.push_back(start + + n_points_per_direction * + n_points_per_direction); + cells.push_back(start + + n_points_per_direction * + n_points_per_direction + + 1); + cells.push_back(start + + n_points_per_direction * + n_points_per_direction + + n_points_per_direction + 1); + cells.push_back(start + + n_points_per_direction * + n_points_per_direction + + n_points_per_direction); } else { - out << start << '\t' << start + 1 << '\t' - << start + n_points_per_direction + 1 - << '\t' << start + n_points_per_direction; - out << '\n'; - } - }; - - for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2) - for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1) - { - const unsigned int starting_offset = - first_vertex_of_patch + - i2 * n_points_per_direction + i1; - write_cell(starting_offset); - } - break; - } - - case 3: - { - auto write_cell = [&flags, - &out, - &cells, - n_points_per_direction]( - const unsigned int start) { - if (deal_ii_with_zlib && - (flags.compression_level != - DataOutBase::CompressionLevel::plain_text)) - { - cells.push_back(start); - cells.push_back(start + 1); - cells.push_back(start + n_points_per_direction + - 1); - cells.push_back(start + n_points_per_direction); - cells.push_back(start + n_points_per_direction * - n_points_per_direction); - cells.push_back(start + - n_points_per_direction * - n_points_per_direction + - 1); - cells.push_back(start + - n_points_per_direction * - n_points_per_direction + - n_points_per_direction + 1); - cells.push_back(start + - n_points_per_direction * - n_points_per_direction + - n_points_per_direction); - } - else - { - out << start << '\t' << start + 1 << '\t' + o << start << '\t' << start + 1 << '\t' << start + n_points_per_direction + 1 << '\t' << start + n_points_per_direction << '\t' << start + n_points_per_direction * @@ -6332,107 +6357,134 @@ namespace DataOutBase n_points_per_direction * n_points_per_direction + n_points_per_direction; - out << '\n'; - } - }; - - for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3) - for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2) - for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1) - { - const unsigned int starting_offset = - first_vertex_of_patch + - i3 * n_points_per_direction * - n_points_per_direction + - i2 * n_points_per_direction + i1; - write_cell(starting_offset); + o << '\n'; } - break; - } + }; - default: - Assert(false, ExcNotImplemented()); - } + for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3) + for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2) + for (unsigned int i1 = 0; i1 < n_subdivisions; + ++i1) + { + const unsigned int starting_offset = + first_vertex_of_patch + + i3 * n_points_per_direction * + n_points_per_direction + + i2 * n_points_per_direction + i1; + write_cell(starting_offset); + } + break; + } - // Finally update the number of the first vertex of this patch - first_vertex_of_patch += - Utilities::fixed_power(n_subdivisions + 1); - } - } + default: + Assert(false, ExcNotImplemented()); + } - // Flush the 'cells' object we created herein. - if (deal_ii_with_zlib && (flags.compression_level != - DataOutBase::CompressionLevel::plain_text)) - { - out << vtu_stringize_array(cells, + // Finally update the number of the first vertex of this patch + first_vertex_of_patch += + Utilities::fixed_power(n_subdivisions + 1); + } + } + + // Flush the 'cells' object we created herein. + if (deal_ii_with_zlib && (flags.compression_level != + DataOutBase::CompressionLevel::plain_text)) + { + o << vtu_stringize_array(cells, flags.compression_level, - out.precision()) + output_precision) << '\n'; - } - } - out << " \n"; + } + } + o << " \n"; - // XML VTU format uses offsets; this is different than the VTK format, which - // puts the number of nodes per cell in front of the connectivity list. - out << " \n"; + return o.str(); + }; + out << stringize_cell_to_vertex_information(); - std::vector offsets; - offsets.reserve(n_cells); - // std::uint8_t might be an alias to unsigned char which is then not printed - // as ascii integers - std::vector cell_types; - cell_types.reserve(n_cells); + //------------------------------- + // The second part of cell information is the offsets in + // the array built by the previous lambda function that indicate + // individual cells. + // + // Note that this separates XML VTU format from the VTK format; the latter + // puts the number of nodes per cell in front of the connectivity list for + // each cell, whereas the VTU format uses one large list of vertex indices + // and a separate array of offsets. + // + // The third piece to cell information is that we need to + // output the types of the cells. + // + // The following function does both of these pieces. + const auto stringize_cell_offset_and_type_information = + [&patches, + &flags, + ascii_or_binary, + n_cells, + output_precision = out.precision()]() { + std::ostringstream o; - unsigned int first_vertex_of_patch = 0; + o << " \n"; - for (const auto &patch : patches) - { - const auto vtk_cell_id = - extract_vtk_patch_info(patch, flags.write_higher_order_cells); + std::vector offsets; + offsets.reserve(n_cells); - for (unsigned int i = 0; i < vtk_cell_id[1]; ++i) + // std::uint8_t might be an alias to unsigned char which is then not + // printed as ascii integers + std::vector cell_types; + cell_types.reserve(n_cells); + + unsigned int first_vertex_of_patch = 0; + + for (const auto &patch : patches) { - cell_types.push_back(vtk_cell_id[0]); - first_vertex_of_patch += vtk_cell_id[2]; - offsets.push_back(first_vertex_of_patch); + const auto vtk_cell_id = + extract_vtk_patch_info(patch, flags.write_higher_order_cells); + + for (unsigned int i = 0; i < vtk_cell_id[1]; ++i) + { + cell_types.push_back(vtk_cell_id[0]); + first_vertex_of_patch += vtk_cell_id[2]; + offsets.push_back(first_vertex_of_patch); + } } - } - out << vtu_stringize_array(offsets, - flags.compression_level, - out.precision()); - out << '\n'; - out << " \n"; + o << vtu_stringize_array(offsets, + flags.compression_level, + output_precision); + o << '\n'; + o << " \n"; - // next output the types of the cells. since all cells are the same, this is - // simple - out << " \n"; + o << " \n"; - // this should compress well :-) - if (deal_ii_with_zlib && - (flags.compression_level != CompressionLevel::plain_text)) - { - std::vector cell_types_uint8_t(cell_types.size()); - for (unsigned int i = 0; i < cell_types.size(); ++i) - cell_types_uint8_t[i] = static_cast(cell_types[i]); + if (deal_ii_with_zlib && + (flags.compression_level != CompressionLevel::plain_text)) + { + std::vector cell_types_uint8_t(cell_types.size()); + for (unsigned int i = 0; i < cell_types.size(); ++i) + cell_types_uint8_t[i] = static_cast(cell_types[i]); - out << vtu_stringize_array(cell_types_uint8_t, - flags.compression_level, - out.precision()); - } - else - { - out << vtu_stringize_array(cell_types, - flags.compression_level, - out.precision()); - } + o << vtu_stringize_array(cell_types_uint8_t, + flags.compression_level, + output_precision); + } + else + { + o << vtu_stringize_array(cell_types, + flags.compression_level, + output_precision); + } - out << '\n'; - out << " \n"; - out << " \n"; + o << '\n'; + o << " \n"; + o << " \n"; + + return o.str(); + }; + out << stringize_cell_offset_and_type_information(); //------------------------------------- @@ -6442,190 +6494,213 @@ namespace DataOutBase // place const Table<2, float> data_vectors = std::move(*create_global_data_table_task.return_value()); + std::vector data_set_written(n_data_sets, false); // then write data. the 'POINT_DATA' means: node data (as opposed to cell // data, which we do not support explicitly here). all following data sets // are point data out << " \n"; - // when writing, first write out all vector data, then handle the scalar - // data sets that have been left over - std::vector data_set_written(n_data_sets, false); - for (const auto &range : nonscalar_data_ranges) - { - const auto first_component = std::get<0>(range); - const auto last_component = std::get<1>(range); - const auto &name = std::get<2>(range); - const bool is_tensor = - (std::get<3>(range) == - DataComponentInterpretation::component_is_part_of_tensor); - const unsigned int n_components = (is_tensor ? 9 : 3); - AssertThrow(last_component >= first_component, - ExcLowerRange(last_component, first_component)); - AssertThrow(last_component < n_data_sets, - ExcIndexRange(last_component, 0, n_data_sets)); - if (is_tensor) - { - AssertThrow((last_component + 1 - first_component <= 9), - ExcMessage( - "Can't declare a tensor with more than 9 components " - "in VTK/VTU format.")); - } - else - { - AssertThrow((last_component + 1 - first_component <= 3), - ExcMessage( - "Can't declare a vector with more than 3 components " - "in VTK/VTU format.")); - } + const auto stringize_nonscalar_data_range = [&flags, + &data_set_written, + &data_names, + &data_vectors, + ascii_or_binary, + n_data_sets, + n_nodes, + output_precision = + out.precision()]( + const auto &range) { + std::ostringstream o; + + const auto first_component = std::get<0>(range); + const auto last_component = std::get<1>(range); + const auto &name = std::get<2>(range); + const bool is_tensor = + (std::get<3>(range) == + DataComponentInterpretation::component_is_part_of_tensor); + const unsigned int n_components = (is_tensor ? 9 : 3); + AssertThrow(last_component >= first_component, + ExcLowerRange(last_component, first_component)); + AssertThrow(last_component < n_data_sets, + ExcIndexRange(last_component, 0, n_data_sets)); + if (is_tensor) + { + AssertThrow((last_component + 1 - first_component <= 9), + ExcMessage( + "Can't declare a tensor with more than 9 components " + "in VTK/VTU format.")); + } + else + { + AssertThrow((last_component + 1 - first_component <= 3), + ExcMessage( + "Can't declare a vector with more than 3 components " + "in VTK/VTU format.")); + } - // mark these components as already written: - for (unsigned int i = first_component; i <= last_component; ++i) - data_set_written[i] = true; + // mark these components as already written: + for (unsigned int i = first_component; i <= last_component; ++i) + data_set_written[i] = true; - // write the header. concatenate all the component names with double - // underscores unless a vector name has been specified - out << " \n"; + o << "\" NumberOfComponents=\"" << n_components << "\" format=\"" + << 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()) + o << " units=\"" << flags.physical_units.at(name) << "\""; + } + else + { + if (flags.physical_units.find(data_names[first_component]) != + flags.physical_units.end()) + o << " units=\"" + << flags.physical_units.at(data_names[first_component]) << "\""; + } + o << ">\n"; - // now write data. pad all vectors to have three components - std::vector data; - data.reserve(n_nodes * n_components); + // now write data. pad all vectors to have three components + std::vector data; + data.reserve(n_nodes * n_components); - for (unsigned int n = 0; n < n_nodes; ++n) - { - if (!is_tensor) - { - switch (last_component - first_component) - { - case 0: - data.push_back(data_vectors(first_component, n)); - data.push_back(0); - data.push_back(0); - break; - - case 1: - data.push_back(data_vectors(first_component, n)); - data.push_back(data_vectors(first_component + 1, n)); - data.push_back(0); - break; - - case 2: - data.push_back(data_vectors(first_component, n)); - data.push_back(data_vectors(first_component + 1, n)); - data.push_back(data_vectors(first_component + 2, n)); - break; - - default: - // Anything else is not yet implemented - Assert(false, ExcInternalError()); - } - } - else - { - Tensor<2, 3> vtk_data; - vtk_data = 0.; + for (unsigned int n = 0; n < n_nodes; ++n) + { + if (!is_tensor) + { + switch (last_component - first_component) + { + case 0: + data.push_back(data_vectors(first_component, n)); + data.push_back(0); + data.push_back(0); + break; - const unsigned int size = last_component - first_component + 1; - if (size == 1) - // 1D, 1 element - { - vtk_data[0][0] = data_vectors(first_component, n); - } - else if (size == 4) - // 2D, 4 elements - { - for (unsigned int c = 0; c < size; ++c) - { - const auto ind = - Tensor<2, 2>::unrolled_to_component_indices(c); - vtk_data[ind[0]][ind[1]] = - data_vectors(first_component + c, n); - } - } - else if (size == 9) - // 3D 9 elements - { - for (unsigned int c = 0; c < size; ++c) - { - const auto ind = - Tensor<2, 3>::unrolled_to_component_indices(c); - vtk_data[ind[0]][ind[1]] = - data_vectors(first_component + c, n); - } - } - else - { + case 1: + data.push_back(data_vectors(first_component, n)); + data.push_back(data_vectors(first_component + 1, n)); + data.push_back(0); + break; + + case 2: + data.push_back(data_vectors(first_component, n)); + data.push_back(data_vectors(first_component + 1, n)); + data.push_back(data_vectors(first_component + 2, n)); + break; + + default: + // Anything else is not yet implemented Assert(false, ExcInternalError()); - } + } + } + else + { + Tensor<2, 3> vtk_data; + vtk_data = 0.; - // now put the tensor into data - // note we padd with zeros because VTK format always wants to - // see a 3x3 tensor, regardless of dimension - for (unsigned int i = 0; i < 3; ++i) - for (unsigned int j = 0; j < 3; ++j) - data.push_back(vtk_data[i][j]); - } - } // loop over nodes + const unsigned int size = last_component - first_component + 1; + if (size == 1) + // 1D, 1 element + { + vtk_data[0][0] = data_vectors(first_component, n); + } + else if (size == 4) + // 2D, 4 elements + { + for (unsigned int c = 0; c < size; ++c) + { + const auto ind = + Tensor<2, 2>::unrolled_to_component_indices(c); + vtk_data[ind[0]][ind[1]] = + data_vectors(first_component + c, n); + } + } + else if (size == 9) + // 3D 9 elements + { + for (unsigned int c = 0; c < size; ++c) + { + const auto ind = + Tensor<2, 3>::unrolled_to_component_indices(c); + vtk_data[ind[0]][ind[1]] = + data_vectors(first_component + c, n); + } + } + else + { + Assert(false, ExcInternalError()); + } - out << vtu_stringize_array(data, - flags.compression_level, - out.precision()); - out << '\n'; - out << " \n"; + // now put the tensor into data + // note we padd with zeros because VTK format always wants to + // see a 3x3 tensor, regardless of dimension + for (unsigned int i = 0; i < 3; ++i) + for (unsigned int j = 0; j < 3; ++j) + data.push_back(vtk_data[i][j]); + } + } // loop over nodes - } // loop over ranges + o << vtu_stringize_array(data, flags.compression_level, output_precision); + o << '\n'; + o << " \n"; - // now do the left over scalar data sets + return o.str(); + }; + + const auto stringize_scalar_data_set = [&flags, + &data_names, + &data_vectors, + ascii_or_binary, + output_precision = out.precision()]( + const unsigned int data_set) { + std::ostringstream o; + + o << " \n"; + + const std::vector data(data_vectors[data_set].begin(), + data_vectors[data_set].end()); + o << vtu_stringize_array(data, flags.compression_level, output_precision); + o << '\n'; + o << " \n"; + + return o.str(); + }; + + + // When writing, first write out all vector and tensor data + for (const auto &range : nonscalar_data_ranges) + { + out << stringize_nonscalar_data_range(range); + } + + // Now do the left over scalar data sets for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set) if (data_set_written[data_set] == false) { - out << " \n"; - - const std::vector data(data_vectors[data_set].begin(), - data_vectors[data_set].end()); - out << vtu_stringize_array(data, - flags.compression_level, - out.precision()); - out << '\n'; - out << " \n"; + out << stringize_scalar_data_set(data_set); } out << " \n"; -- 2.39.5