From: Denis Davydov Date: Thu, 21 Jun 2018 10:41:16 +0000 (+0200) Subject: add component_is_part_of_tensor and its support to vtu X-Git-Tag: v9.1.0-rc1~950^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=58a148a855c93c5b09e4070df592f9ae18dd4a22;p=dealii.git add component_is_part_of_tensor and its support to vtu --- diff --git a/doc/news/changes/minor/20180621DenisDavydov_b b/doc/news/changes/minor/20180621DenisDavydov_b new file mode 100644 index 0000000000..cd9ef39ad6 --- /dev/null +++ b/doc/news/changes/minor/20180621DenisDavydov_b @@ -0,0 +1,3 @@ +New: Add support for output of tensor-valued data in +DataOut::write_vtu(). +(Denis Davydov 2018/06/21) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 1b3b06a06e..88919243c8 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -1357,6 +1357,46 @@ Tensor::component_to_unrolled_index( } + +namespace internal +{ + // unrolled_to_component_indices is instantiated from DataOut for dim==0 + // and rank=2. Make sure we don't have compiler warnings. + + template + inline unsigned int + mod(const unsigned int x) + { + return x % dim; + } + + template <> + inline unsigned int + mod<0>(const unsigned int x) + { + Assert(false, ExcInternalError()); + return x; + } + + template + inline unsigned int + div(const unsigned int x) + { + return x / dim; + } + + template <> + inline unsigned int + div<0>(const unsigned int x) + { + Assert(false, ExcInternalError()); + return x; + } + +} // namespace internal + + + template inline TableIndices Tensor::unrolled_to_component_indices(const unsigned int i) @@ -1369,8 +1409,8 @@ Tensor::unrolled_to_component_indices(const unsigned int i) unsigned int remainder = i; for (int r = rank_ - 1; r >= 0; --r) { - indices[r] = (remainder % dim); - remainder /= dim; + indices[r] = internal::mod(remainder); + remainder = internal::div(remainder); } Assert(remainder == 0, ExcInternalError()); diff --git a/include/deal.II/numerics/data_component_interpretation.h b/include/deal.II/numerics/data_component_interpretation.h index 91611e4ab0..cc722cd639 100644 --- a/include/deal.II/numerics/data_component_interpretation.h +++ b/include/deal.II/numerics/data_component_interpretation.h @@ -58,7 +58,13 @@ namespace DataComponentInterpretation * Indicates that a component of a data set is part of a vector-valued * quantity. */ - component_is_part_of_vector + component_is_part_of_vector, + + /** + * Indicates that a component of a data set is part of a tensor-valued + * (2nd order) quantity. + */ + component_is_part_of_tensor }; } // namespace DataComponentInterpretation diff --git a/include/deal.II/numerics/data_out_dof_data.h b/include/deal.II/numerics/data_out_dof_data.h index 88e0bca599..ea11144f7e 100644 --- a/include/deal.II/numerics/data_out_dof_data.h +++ b/include/deal.II/numerics/data_out_dof_data.h @@ -149,6 +149,19 @@ namespace Exceptions << "components. However, the vector component at " << "position " << arg1 << " with name <" << arg2 << "> does not satisfy these conditions."); + + DeclException2(ExcInvalidTensorDeclaration, + int, + std::string, + << "When declaring that a number of components in a data " + << "set to be output logically form a tensor instead of " + << "simply a set of scalar fields, you need to specify " + << "this for all relevant components. Furthermore, " + << "tensors must always consist of exactly " + << "components. However, the tensor component at " + << "position " << arg1 << " with name <" << arg2 + << "> does not satisfy these conditions."); + } // namespace DataOutImplementation } // namespace Exceptions diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h index e93ef62ac7..dd5f5343f0 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -1212,8 +1212,8 @@ DataOut_DoFData:: if (actual_type == type_dof_data) { // output vectors for which at least part of the data is to be interpreted - // as vector fields cannot be complex-valued, because we cannot visualize - // complex-valued vector fields + // as vector or tensor fields cannot be complex-valued, because we cannot + // visualize complex-valued vector fields Assert(!((std::find( data_component_interpretation.begin(), data_component_interpretation.end(), @@ -1228,6 +1228,20 @@ DataOut_DoFData:: "this vector as a collection of scalar fields that can then " "be visualized by their real and imaginary parts separately.")); + Assert(!((std::find( + data_component_interpretation.begin(), + data_component_interpretation.end(), + DataComponentInterpretation::component_is_part_of_tensor) != + data_component_interpretation.end()) && + new_entry->is_complex_valued()), + ExcMessage( + "Complex-valued vectors added to a DataOut-like object " + "cannot contain components that shall be interpreted as " + "tensor fields because one can not visualize complex-valued " + "tensor fields. However, you may want to try to output " + "this vector as a collection of scalar fields that can then " + "be visualized by their real and imaginary parts separately.")); + dof_data.emplace_back(std::move(new_entry)); } else @@ -1391,6 +1405,45 @@ DataOut_DoFData:: output_component += patch_space_dim; i += patch_space_dim; } + else if ((*d)->data_component_interpretation[i] == + DataComponentInterpretation::component_is_part_of_tensor) + { + const unsigned int size = patch_space_dim * patch_space_dim; + // ensure that there is a continuous number of next + // space_dim*space_dim components that all deal with tensors + Assert(i + size <= (*d)->n_output_variables, + Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration( + i, (*d)->names[i])); + for (unsigned int dd = 1; dd < size; ++dd) + Assert( + (*d)->data_component_interpretation[i + dd] == + DataComponentInterpretation::component_is_part_of_tensor, + Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration( + i, (*d)->names[i])); + + // all seems alright, so figure out whether there is a common name to + // these components. if not, leave the name empty and let the output + // format writer decide what to do here + std::string name = (*d)->names[i]; + for (unsigned int dd = 1; dd < size; ++dd) + if (name != (*d)->names[i + dd]) + { + name = ""; + break; + } + + // finally add a corresponding range + ranges.push_back(std::make_tuple( + output_component, + output_component + size - 1, + name, + DataComponentInterpretation::component_is_part_of_tensor)); + + // increase the 'component' counter by the appropriate amount, same + // for 'i', since we have already dealt with all these components + output_component += size; + i += size; + } else { // just move one component forward by one, or two if the vector diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index d064def0c3..a6f1f86197 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -38,6 +38,8 @@ #include #include +#include + #include #include #include @@ -5571,14 +5573,28 @@ namespace DataOutBase 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)); - AssertThrow(last_component + 1 - first_component <= 3, - ExcMessage( - "Can't declare a vector with more than 3 components " - "in VTK")); + if (is_tensor) + { + AssertThrow((last_component + 1 - first_component <= 9), + ExcMessage( + "Can't declare a tensor with more than 9 components " + "in VTK")); + } + else + { + AssertThrow((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 = first_component; i <= last_component; ++i) @@ -5597,43 +5613,82 @@ namespace DataOutBase out << data_names[last_component]; } - out << "\" NumberOfComponents=\"3\" format=\"" << ascii_or_binary - << "\">\n"; + out << "\" NumberOfComponents=\"" << n_components << "\" format=\"" + << ascii_or_binary << "\">\n"; // now write data. pad all vectors to have three components std::vector data; - data.reserve(n_nodes * dim); + data.reserve(n_nodes * n_components); for (unsigned int n = 0; n < n_nodes; ++n) { - switch (last_component - first_component) + if (!is_tensor) { - case 0: - data.push_back(data_vectors(first_component, n)); - data.push_back(0); - data.push_back(0); - break; + 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.; - 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; + 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) || (size == 9)) + // 2D, 4 elements or 3D 9 elements + { + for (unsigned int c = 0; c < size; ++c) + { + const auto ind = + Tensor<2, dim>::unrolled_to_component_indices(c); + vtk_data[ind[0]][ind[1]] = + data_vectors(first_component + c, n); + } + } + else + { + Assert(false, ExcInternalError()); + } - default: - // VTK doesn't support anything else than vectors with 1, 2, - // or 3 components - Assert(false, ExcInternalError()); + // 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 + vtu_out << data; out << " \n"; - } + + } // loop over ranges // now do the left over scalar data sets for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set) diff --git a/tests/data_out/data_out_12.cc b/tests/data_out/data_out_12.cc new file mode 100644 index 0000000000..d4a064a8bf --- /dev/null +++ b/tests/data_out/data_out_12.cc @@ -0,0 +1,193 @@ +// --------------------------------------------------------------------- +// +// 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 vtu + +#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::ofstream out("output_" + Utilities::int_to_string(dim) + "d.vtu"); + data_out.write_vtu(out); + */ + + data_out.write_vtu(deallog.get_file_stream()); +} + + + +int +main() +{ + initlog(); + check<2>(); + check<3>(); +} diff --git a/tests/data_out/data_out_12.output b/tests/data_out/data_out_12.output new file mode 100644 index 0000000000..e1298df815 --- /dev/null +++ b/tests/data_out/data_out_12.output @@ -0,0 +1,69 @@ + + + + + + + + +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/ + + + +