]> https://gitweb.dealii.org/ - dealii.git/commitdiff
adapt pvtu output for tensors 7755/head
authornotxor <pmnotxor@gmx.at>
Wed, 27 Feb 2019 08:11:42 +0000 (09:11 +0100)
committernotxor <pmnotxor@gmx.at>
Wed, 27 Feb 2019 08:22:46 +0000 (09:22 +0100)
doc/news/changes/minor/20190227DanielJodlbauer [new file with mode: 0644]
source/base/data_out_base.cc
tests/data_out/data_out_13.cc [new file with mode: 0644]
tests/data_out/data_out_13.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190227DanielJodlbauer b/doc/news/changes/minor/20190227DanielJodlbauer
new file mode 100644 (file)
index 0000000..a7afb2b
--- /dev/null
@@ -0,0 +1,3 @@
+Fixed: Tensor-valued pvtu output.
+<br>
+(Daniel Jodlbauer, 2019/02/27)
index 65336ef191547ad05e7a69c6ae7aa86ea6b9d434..939c0c905042b9e4f28f9fb55b7c27391210f5e3 100644 (file)
@@ -5990,20 +5990,30 @@ namespace DataOutBase
     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);
@@ -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 (file)
index 0000000..582285c
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/data_out/data_out_13.output b/tests/data_out/data_out_13.output
new file mode 100644 (file)
index 0000000..f410ee9
--- /dev/null
@@ -0,0 +1,101 @@
+
+<?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>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.