]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add component_is_part_of_tensor and its support to vtu 6807/head
authorDenis Davydov <davydden@gmail.com>
Thu, 21 Jun 2018 10:41:16 +0000 (12:41 +0200)
committerDenis Davydov <davydden@gmail.com>
Thu, 5 Jul 2018 07:31:06 +0000 (09:31 +0200)
doc/news/changes/minor/20180621DenisDavydov_b [new file with mode: 0644]
include/deal.II/base/tensor.h
include/deal.II/numerics/data_component_interpretation.h
include/deal.II/numerics/data_out_dof_data.h
include/deal.II/numerics/data_out_dof_data.templates.h
source/base/data_out_base.cc
tests/data_out/data_out_12.cc [new file with mode: 0644]
tests/data_out/data_out_12.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180621DenisDavydov_b b/doc/news/changes/minor/20180621DenisDavydov_b
new file mode 100644 (file)
index 0000000..cd9ef39
--- /dev/null
@@ -0,0 +1,3 @@
+New: Add support for output of tensor-valued data in
+DataOut::write_vtu().
+(Denis Davydov 2018/06/21)
index 1b3b06a06ea2bf34dd64b0fdd33e38bfaa9b20f9..88919243c87b1597fa334d7acc8bfd3dd3d2abbf 100644 (file)
@@ -1357,6 +1357,46 @@ Tensor<rank_, dim, Number>::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 <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)
@@ -1369,8 +1409,8 @@ 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());
 
index 91611e4ab065b64ef7318f05f72548ce538cc69f..cc722cd639ba5ef4739861b62a75c48f0369440c 100644 (file)
@@ -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
 
index 88e0bca5990c5b35847a97f063e2e1108aee5576..ea11144f7e3f3d5ecdced89e2fc3ecf47133ba7f 100644 (file)
@@ -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 <dim*dim> "
+                   << "components. However, the tensor component at "
+                   << "position " << arg1 << " with name <" << arg2
+                   << "> does not satisfy these conditions.");
+
   } // namespace DataOutImplementation
 } // namespace Exceptions
 
index e93ef62ac7077b7beb628c1c975ec7581cae965b..dd5f5343f0ab441274326a9a478ac13e7dafb5e8 100644 (file)
@@ -1212,8 +1212,8 @@ DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
   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<DoFHandlerType, patch_dim, patch_space_dim>::
                "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<DoFHandlerType, patch_dim, patch_space_dim>::
           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
index d064def0c3ce2453d9e8f84385c0dcfa121142e5..a6f1f8619711994260932dd23ce9246b8a5aaf7e 100644 (file)
@@ -38,6 +38,8 @@
 #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>
@@ -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<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)
diff --git a/tests/data_out/data_out_12.cc b/tests/data_out/data_out_12.cc
new file mode 100644 (file)
index 0000000..d4a064a
--- /dev/null
@@ -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 <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>();
+}
diff --git a/tests/data_out/data_out_12.output b/tests/data_out/data_out_12.output
new file mode 100644 (file)
index 0000000..e1298df
--- /dev/null
@@ -0,0 +1,69 @@
+
+<?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>

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.