--- /dev/null
+New: Add support for output of tensor-valued data in
+DataOut::write_vtu().
+(Denis Davydov 2018/06/21)
}
+
+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 <int dim>
+ 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 <int dim>
+ 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 <int rank_, int dim, typename Number>
inline TableIndices<rank_>
Tensor<rank_, dim, Number>::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<dim>(remainder);
+ remainder = internal::div<dim>(remainder);
}
Assert(remainder == 0, ExcInternalError());
* 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
<< "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 <dim*dim> "
+ << "components. However, the tensor component at "
+ << "position " << arg1 << " with name <" << arg2
+ << "> does not satisfy these conditions.");
+
} // namespace DataOutImplementation
} // namespace Exceptions
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(),
"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
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
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
+#include <deal.II/numerics/data_component_interpretation.h>
+
#include <algorithm>
#include <cmath>
#include <cstring>
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)
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<float> 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 << " </DataArray>\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)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <fstream>
+#include <iomanip>
+#include <string>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+check()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria, 0., 1.);
+
+ FESystem<dim> fe(FE_Q<dim>(1),
+ dim, // vector valued
+ FE_Q<dim>(1),
+ dim * dim); // tensor valued
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> 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<std::vector<std::vector<double>>> 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<std::vector<double>> 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<std::vector<std::vector<double>>> tensors_2d = {
+ {{1, 0}, {0, 1}}, {{3, 0}, {0, 1}}, {{1, 0}, {0, 3}}, {{3, 1}, {1, 2}}};
+
+ // coords:
+ const std::vector<std::vector<double>> coords_2d = {{0, 0},
+ {1, 0},
+ {0, 1},
+ {1, 1}};
+
+ const std::vector<std::vector<std::vector<double>>> &tensors =
+ (dim == 2 ? tensors_2d : tensors_3d);
+ const std::vector<std::vector<double>> &coords =
+ (dim == 2 ? coords_2d : coords_3d);
+
+ const Quadrature<dim> support_quadrature(fe.get_unit_support_points());
+ FEValues<dim> 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<dim> 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<DataComponentInterpretation::DataComponentInterpretation>
+ 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<std::string> component_name(dim + dim * dim, "tensor");
+ std::fill(component_name.begin(), component_name.begin() + dim, "vector");
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(v,
+ component_name,
+ DataOut<dim>::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>();
+}
--- /dev/null
+
+<?xml version="1.0" ?>
+<!--
+# vtk DataFile Version 3.0
+#This file was generated
+-->
+<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="4" NumberOfCells="1" >
+ <Points>
+ <DataArray type="Float32" NumberOfComponents="3" format="binary">
+AQAAADAAAAAwAAAAFAAAAA==eNpjYEAGDfYMWPkgGsIGADHwAv0=
+ </DataArray>
+ </Points>
+
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" format="binary">
+AQAAABAAAAAQAAAAEwAAAA==eNpjYGBgYARiZiBmAmIAADwABw==
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" format="binary">
+AQAAAAQAAAAEAAAADAAAAA==eNpjYWBgAAAAFAAF
+ </DataArray>
+ <DataArray type="UInt8" Name="types" format="binary">
+AQAAAAEAAAABAAAACQAAAA==eNrjBAAACgAK
+ </DataArray>
+ </Cells>
+ <PointData Scalars="scalars">
+ <DataArray type="Float32" Name="vector" NumberOfComponents="3" format="binary">
+AQAAADAAAAAwAAAAFAAAAA==eNpjYGiwZ4ADdDaMj2ADAGQuBHs= </DataArray>
+ <DataArray type="Float32" Name="tensor" NumberOfComponents="9" format="binary">
+AQAAAJAAAACQAAAAIgAAAA==eNpjYGiwZ0AB6HwQcHAgrAZdDF0PTAymDkxjqAEAycUGOw== </DataArray>
+ </PointData>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
+<?xml version="1.0" ?>
+<!--
+# vtk DataFile Version 3.0
+#This file was generated
+-->
+<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
+<UnstructuredGrid>
+<Piece NumberOfPoints="8" NumberOfCells="1" >
+ <Points>
+ <DataArray type="Float32" NumberOfComponents="3" format="binary">
+AQAAAGAAAABgAAAAGgAAAA==eNpjYEAGDfYMWPkgGpscsjyyGAIDAC3jCPU=
+ </DataArray>
+ </Points>
+
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" format="binary">
+AQAAACAAAAAgAAAAHQAAAA==eNpjYGBgYARiZiBmAmIWIGYFYnYgZgNiAAF4AB0=
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" format="binary">
+AQAAAAQAAAAEAAAADAAAAA==eNrjYGBgAAAAJAAJ
+ </DataArray>
+ <DataArray type="UInt8" Name="types" format="binary">
+AQAAAAEAAAABAAAACQAAAA==eNrjAQAADQAN
+ </DataArray>
+ </Cells>
+ <PointData Scalars="scalars">
+ <DataArray type="Float32" Name="vector" NumberOfComponents="3" format="binary">
+AQAAAGAAAABgAAAAFgAAAA==eNpjYGiwZ4ADdDaMj84mXj0An0sI9Q== </DataArray>
+ <DataArray type="Float32" Name="tensor" NumberOfComponents="9" format="binary">
+AQAAACABAAAgAQAAhAAAAA==eNpjYGiwZ0ABWPkHUMXQ+Q4ODATMyRdqtgIphIjDMDrfwaFArtZuZ/JvO44qNjDN/HvvPlFRC3sQH0TPOSC/T39lv51GyHs78ch8exBt86TUlrM4Yh+ID6KntZbZ2+9tsE3Vc9zz44XWfhB9+8/xfbuNSveD+CB63lfXfQxEAACmEDx/ </DataArray>
+ </PointData>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
projects / dealii.git / commitdiff