From 34f9c07b2d0cad24d0d355d3a9ca3e1c8dbfaebe Mon Sep 17 00:00:00 2001 From: notxor Date: Wed, 27 Feb 2019 09:11:42 +0100 Subject: [PATCH] adapt pvtu output for tensors --- .../changes/minor/20190227DanielJodlbauer | 3 + source/base/data_out_base.cc | 41 ++-- tests/data_out/data_out_13.cc | 200 ++++++++++++++++++ tests/data_out/data_out_13.output | 101 +++++++++ 4 files changed, 330 insertions(+), 15 deletions(-) create mode 100644 doc/news/changes/minor/20190227DanielJodlbauer create mode 100644 tests/data_out/data_out_13.cc create mode 100644 tests/data_out/data_out_13.output diff --git a/doc/news/changes/minor/20190227DanielJodlbauer b/doc/news/changes/minor/20190227DanielJodlbauer new file mode 100644 index 0000000000..a7afb2bcae --- /dev/null +++ b/doc/news/changes/minor/20190227DanielJodlbauer @@ -0,0 +1,3 @@ +Fixed: Tensor-valued pvtu output. +
+(Daniel Jodlbauer, 2019/02/27) diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index 65336ef191..939c0c9050 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -5990,20 +5990,30 @@ namespace DataOutBase std::vector data_set_written(n_data_sets, false); for (const auto &nonscalar_data_range : nonscalar_data_ranges) { - AssertThrow(std::get<1>(nonscalar_data_range) >= - std::get<0>(nonscalar_data_range), - ExcLowerRange(std::get<1>(nonscalar_data_range), - std::get<0>(nonscalar_data_range))); - AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets, - ExcIndexRange(std::get<1>(nonscalar_data_range), - 0, - n_data_sets)); - AssertThrow(std::get<1>(nonscalar_data_range) + 1 - - std::get<0>(nonscalar_data_range) <= - 3, - ExcMessage( - "Can't declare a vector with more than 3 components " - "in VTK")); + const auto first_component = std::get<0>(nonscalar_data_range); + const auto last_component = std::get<1>(nonscalar_data_range); + const bool is_tensor = + (std::get<3>(nonscalar_data_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")); + } + else + { + Assert((last_component + 1 - first_component <= 3), + ExcMessage( + "Can't declare a vector with more than 3 components " + "in VTK")); + } // mark these components as already written: for (unsigned int i = std::get<0>(nonscalar_data_range); @@ -6026,7 +6036,8 @@ namespace DataOutBase out << data_names[std::get<1>(nonscalar_data_range)]; } - out << "\" NumberOfComponents=\"3\" format=\"ascii\"/>\n"; + out << "\" NumberOfComponents=\"" << n_components + << "\" format=\"ascii\"/>\n"; } for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set) diff --git a/tests/data_out/data_out_13.cc b/tests/data_out/data_out_13.cc new file mode 100644 index 0000000000..582285c0d2 --- /dev/null +++ b/tests/data_out/data_out_13.cc @@ -0,0 +1,200 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Output tensor-valued data in pvtu + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + + + +template +void +check() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, 0., 1.); + + FESystem fe(FE_Q(1), + dim, // vector valued + FE_Q(1), + dim * dim); // tensor valued + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector v(dof_handler.n_dofs()); + + + // 1 cell is enough to test: + const unsigned int dofs_per_cell = v.size(); + + // take 8 tensors from tensors.vtk file inside + // http://paraview.org/Wiki/images/8/81/TensorGlyph_v_1_0_3.zip or + // http://paraview.org/Wiki/images/3/39/SuperquadricTensorGlyph_v3.zip + const std::vector>> tensors_3d = { + {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, + {{-2, 0, 0}, {0, -2, 0}, {0, 0, 3}}, + {{1, 0, 0}, {0, 1, 0}, {0, 0, .001}}, + {{2, 1, 1}, {1, 2, 1}, {1, 1, 3}}, + {{0.247186415568, 0.490995206139, 0.131324884836}, + {0.490995206139, -0.371055707211, 0.719071682671}, + {0.131324884836, 0.719071682671, -0.156008182087}}, + {{0.280587657181, 0.467438945439, 0.934953136331}, + {0.467438945439, 0.0600321140579, -0.211376327727}, + {0.934953136331, -0.211376327727, 0.962975830149}}, + {{0.0628609549443, -0.0117908465212, -0.667617347633}, + {-0.0117908465212, -0.390601011028, -0.95780533241}, + {-0.667617347633, -0.95780533241, -0.193319773383}}, + {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}}; + + // coords: + const std::vector> coords_3d = {{0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {1, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {0, 1, 1}, + {1, 1, 1}}; + + const std::vector>> tensors_2d = { + {{1, 0}, {0, 1}}, {{3, 0}, {0, 1}}, {{1, 0}, {0, 3}}, {{3, 1}, {1, 2}}}; + + // coords: + const std::vector> coords_2d = {{0, 0}, + {1, 0}, + {0, 1}, + {1, 1}}; + + const std::vector>> &tensors = + (dim == 2 ? tensors_2d : tensors_3d); + const std::vector> &coords = + (dim == 2 ? coords_2d : coords_3d); + + const Quadrature support_quadrature(fe.get_unit_support_points()); + FEValues fe_values(fe, support_quadrature, update_quadrature_points); + for (const auto cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + const auto &qp = fe_values.get_quadrature_points(); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const auto & this_qp = qp[i]; + const unsigned int i_group = fe.system_to_base_index(i).first.first; + const unsigned int i_comp = fe.system_to_component_index(i).first; + if (i_group == 0) + // vector + { + if (i_comp == 0) + v(i) = 1; + else if (i_comp == 1) + v(i) = this_qp[1]; + else + v(i) = 0.; + } + else + // tensor + { + // find matching location + unsigned int p = 0; + for (const auto c : coords) + { + Point point; + for (unsigned int d = 0; d < dim; ++d) + point[d] = c[d]; + + const auto diff = point - this_qp; + if (diff.norm() < 1e-10) + break; + + ++p; + } + + const unsigned int tens_comp = i_comp - dim; + const unsigned int ii = + Tensor<2, dim>::unrolled_to_component_indices(tens_comp)[0]; + const unsigned int jj = + Tensor<2, dim>::unrolled_to_component_indices(tens_comp)[1]; + v(i) = tensors[p][ii][jj]; + } + } + } + + std::vector + data_component_interpretation( + dim + dim * dim, + DataComponentInterpretation::component_is_part_of_tensor); + std::fill(data_component_interpretation.begin(), + data_component_interpretation.begin() + dim, + DataComponentInterpretation::component_is_part_of_vector); + + std::vector component_name(dim + dim * dim, "tensor"); + std::fill(component_name.begin(), component_name.begin() + dim, "vector"); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(v, + component_name, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches(); + + std::vector filenames; + filenames.push_back("output_" + Utilities::int_to_string(dim) + "d.vtu"); + + /* + std::ofstream master_output("output_" + Utilities::int_to_string(dim) + + "d.pvtu"); data_out.write_pvtu_record(master_output, filenames); + + std::ofstream out(filenames[0]); + data_out.write_vtu(out); + */ + + data_out.write_vtu(deallog.get_file_stream()); + data_out.write_pvtu_record(deallog.get_file_stream(), filenames); +} + + + +int +main() +{ + initlog(); + check<2>(); + check<3>(); +} diff --git a/tests/data_out/data_out_13.output b/tests/data_out/data_out_13.output new file mode 100644 index 0000000000..f410ee99e0 --- /dev/null +++ b/tests/data_out/data_out_13.output @@ -0,0 +1,101 @@ + + + + + + + + +AQAAADAAAAAwAAAAFAAAAA==eNpjYEAGDfYMWPkgGsIGADHwAv0= + + + + + +AQAAABAAAAAQAAAAEwAAAA==eNpjYGBgYARiZiBmAmIAADwABw== + + +AQAAAAQAAAAEAAAADAAAAA==eNpjYWBgAAAAFAAF + + +AQAAAAEAAAABAAAACQAAAA==eNrjBAAACgAK + + + + +AQAAADAAAAAwAAAAFAAAAA==eNpjYGiwZ4ADdDaMj2ADAGQuBHs= + +AQAAAJAAAACQAAAAIgAAAA==eNpjYGiwZ0AB6HwQcHAgrAZdDF0PTAymDkxjqAEAycUGOw== + + + + + + + + + + + + + + + + + + + + + + + + + +AQAAAGAAAABgAAAAGgAAAA==eNpjYEAGDfYMWPkgGpscsjyyGAIDAC3jCPU= + + + + + +AQAAACAAAAAgAAAAHQAAAA==eNpjYGBgYARiZiBmAmIWIGYFYnYgZgNiAAF4AB0= + + +AQAAAAQAAAAEAAAADAAAAA==eNrjYGBgAAAAJAAJ + + +AQAAAAEAAAABAAAACQAAAA==eNrjAQAADQAN + + + + +AQAAAGAAAABgAAAAFgAAAA==eNpjYGiwZ4ADdDaMj84mXj0An0sI9Q== + +AQAAACABAAAgAQAAhAAAAA==eNpjYGiwZ0ABWPkHUMXQ+Q4ODATMyRdqtgIphIjDMDrfwaFArtZuZ/JvO44qNjDN/HvvPlFRC3sQH0TPOSC/T39lv51GyHs78ch8exBt86TUlrM4Yh+ID6KntZbZ2+9tsE3Vc9zz44XWfhB9+8/xfbuNSveD+CB63lfXfQxEAACmEDx/ + + + + + + + + + + + + + + + + + + -- 2.39.5