std::vector<bool> 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);
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)
--- /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 pvtu
+
+#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::vector<std::string> 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>();
+}
--- /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"?>
+<!--
+#This file was generated
+-->
+<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <PUnstructuredGrid GhostLevel="0">
+ <PPointData Scalars="scalars">
+ <PDataArray type="Float32" Name="vector" NumberOfComponents="3" format="ascii"/>
+ <PDataArray type="Float32" Name="tensor" NumberOfComponents="9" format="ascii"/>
+ </PPointData>
+ <PPoints>
+ <PDataArray type="Float32" NumberOfComponents="3"/>
+ </PPoints>
+ <Piece Source="output_2d.vtu"/>
+ </PUnstructuredGrid>
+</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>
+<?xml version="1.0"?>
+<!--
+#This file was generated
+-->
+<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <PUnstructuredGrid GhostLevel="0">
+ <PPointData Scalars="scalars">
+ <PDataArray type="Float32" Name="vector" NumberOfComponents="3" format="ascii"/>
+ <PDataArray type="Float32" Name="tensor" NumberOfComponents="9" format="ascii"/>
+ </PPointData>
+ <PPoints>
+ <PDataArray type="Float32" NumberOfComponents="3"/>
+ </PPoints>
+ <Piece Source="output_3d.vtu"/>
+ </PUnstructuredGrid>
+</VTKFile>
projects / dealii.git / commitdiff