]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New GridIn::read_vtk version.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 7 Nov 2018 21:01:17 +0000 (22:01 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 15 Nov 2018 19:37:57 +0000 (20:37 +0100)
This version works also for dim == 1, and is more robust w.r.t. the VTK
file format.

include/deal.II/grid/grid_in.h
source/grid/grid_in.cc
tests/grid/grid_in_vtk_2d/mesh

index 58c3dab83f361ee297d9a76a522cecb686001061..d646ead8f41f7bc11b24bdfc49d0d222872d077b 100644 (file)
@@ -361,9 +361,42 @@ public:
   read(const std::string &in, Format format = Default);
 
   /**
-   * Read grid data from an vtk file. Numerical data is ignored.
+   * Read grid data from a unstructured vtk file. The vtk file may contain
+   * the following VTK cell types: VTK_HEXAHEDRON, VTK_QUAD, and VTK_LINE.
+   *
+   * Depending on the template dimension, only some of the above are accepted.
+   *
+   * In particular, in three dimensions, this function expects the file to
+   * contain
+   *
+   * - VTK_HEXAHEDRON cell types
+   * - VTK_QUAD cell types, to specify optional boundary or interior quad faces
+   * - VTK_LINE cell types, to specify optional boundary or interior edges
+   *
+   * In two dimensions:
+   *
+   * - VTK_QUAD cell types
+   * - VTK_LINE cell types, to specify optional boundary or interior edges
+   *
+   * In one dimension
+   *
+   * - VTK_LINE cell types
+   *
+   * The input file may specify boundary ids, material ids, and manifold ids
+   * using the CELL_DATA section of the
+   * [VTK file format](http://www.vtk.org/VTK/img/file-formats.pdf).
+   *
+   * This function interprets two types of CELL_DATA contained in the input
+   * file: `SCALARS MaterialID`, used to specify the material_id of the cells,
+   * or the boundary_id of the faces and edges, and `SCALARS ManifoldID`, that
+   * can be used to specify the manifold id of any Triangulation object (cell,
+   * face, or edge).
+   *
+   * The companion GridOut::write_vtk function can be used to write VTK files
+   * compatible with this method.
    *
    * @author Mayank Sabharwal, Andreas Putz, 2013
+   * @author Luca Heltai, 2018
    */
   void
   read_vtk(std::istream &in);
index 37ba68e4ac3d5ec1b9226208a8e761b34415b4de..e003ac7c087d8f25aa081b93762d52723d1d179a 100644 (file)
@@ -113,7 +113,6 @@ template <int dim, int spacedim>
 void
 GridIn<dim, spacedim>::read_vtk(std::istream &in)
 {
-  Assert((dim == 2) || (dim == 3), ExcNotImplemented());
   std::string line;
 
   // verify that the first, third and fourth lines match
@@ -157,8 +156,7 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
       unsigned int n_vertices;
       in >> n_vertices;
 
-      in.ignore(256,
-                '\n'); // ignoring the number beside the total no. of points.
+      in >> keyword; // float, double, int, char, etc.
 
       for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
         {
@@ -177,30 +175,15 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                 ExcMessage(
                   "While reading VTK file, failed to find POINTS section"));
 
-
-  //////////////////ignoring space between points and cells
-  /// sections////////////////////
-  std::string checkline;
-  int         no;
-  in.ignore(256,
-            '\n'); // this move pointer to the next line ignoring unwanted no.
-  no = in.tellg();
-  getline(in, checkline);
-  if (checkline.compare("") != 0)
-    {
-      in.seekg(no);
-    }
-
   in >> keyword;
 
-  ///////////////////Processing the CELLS section that contains cells(cells) and
-  /// bound_quads(subcelldata)///////////////////////
+  unsigned int n_geometric_objects = 0;
+  unsigned int n_ints;
 
   if (keyword == "CELLS")
     {
-      unsigned int n_geometric_objects;
       in >> n_geometric_objects;
-      in.ignore(256, '\n');
+      in >> n_ints; // Ignore this, since we don't need it.
 
       if (dim == 3)
         {
@@ -227,6 +210,11 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 
               else if (type == 4)
                 {
+                  // we assume that the file contains first all cells,
+                  // then all faces, and finally all lines
+                  AssertThrow(subcelldata.boundary_lines.size() == 0,
+                              ExcNotImplemented());
+
                   subcelldata.boundary_quads.emplace_back();
 
                   for (unsigned int j = 0; j < type;
@@ -235,6 +223,16 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 
                   subcelldata.boundary_quads.back().material_id = 0;
                 }
+              else if (type == 2)
+                {
+                  subcelldata.boundary_lines.emplace_back();
+
+                  for (unsigned int j = 0; j < type;
+                       j++) // loop to feed the data to the boundary
+                    in >> subcelldata.boundary_lines.back().vertices[j];
+
+                  subcelldata.boundary_lines.back().material_id = 0;
+                }
 
               else
                 AssertThrow(
@@ -285,7 +283,30 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                 AssertThrow(
                   false,
                   ExcMessage(
-                    "While reading VTK file, unknown file type encountered"));
+                    "While reading VTK file, unknown cell type encountered"));
+            }
+        }
+      else if (dim == 1)
+        {
+          for (unsigned int count = 0; count < n_geometric_objects; count++)
+            {
+              unsigned int type;
+              in >> type;
+
+              if (type == 2)
+                {
+                  cells.emplace_back();
+
+                  for (unsigned int j = 0; j < type; j++) // loop to feed data
+                    in >> cells.back().vertices[j];
+
+                  cells.back().material_id = 0;
+                }
+              else
+                AssertThrow(
+                  false,
+                  ExcMessage(
+                    "While reading VTK file, unknown cell type encountered"));
             }
         }
       else
@@ -301,98 +322,144 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
       if (keyword ==
           "CELL_TYPES") // Entering the cell_types section and ignoring data.
         {
-          in.ignore(256, '\n');
-
-          while (!in.fail() && !in.eof())
-            {
-              in >> keyword;
-              if (keyword != "12" && keyword != "9")
-                {
-                  break;
-                }
-            }
+          in >> n_ints;
+          int tmp_int;
+          for (unsigned int i = 0; i < n_ints; ++i)
+            in >> tmp_int;
         }
+      else
+        AssertThrow(
+          false,
+          ExcMessage(std::string(
+            "While reading VTK file, missing CELL_TYPES section. Found <" +
+            keyword + "> instead.")));
+
+      // Ignore everything up to CELL_DATA
+      while (in >> keyword)
+        if (keyword == "CELL_DATA")
+          {
+            unsigned int n_ids;
+            in >> n_ids;
+
+            AssertThrow(n_ids == n_geometric_objects,
+                        ExcMessage("The VTK reader found a CELL_DATA statement "
+                                   "that lists a total of " +
+                                   Utilities::int_to_string(n_ids) +
+                                   " cell data objects, but this needs to "
+                                   "equal the number of cells (which is " +
+                                   Utilities::int_to_string(cells.size()) +
+                                   ") plus the number of quads (" +
+                                   Utilities::int_to_string(
+                                     subcelldata.boundary_quads.size()) +
+                                   " in 3d or the number of lines (" +
+                                   Utilities::int_to_string(
+                                     subcelldata.boundary_lines.size()) +
+                                   ") in 2d."));
+
+            const std::vector<std::string> data_sets{"MaterialID",
+                                                     "ManifoldID"};
+
+            for (unsigned int i = 0; i < data_sets.size(); ++i)
+              {
+                // Ignore everything until we get to a SCALARS data set
+                while (in >> keyword)
+                  if (keyword == "SCALARS")
+                    {
+                      // Now see if we know about this type of data set,
+                      // if not, just ignore everything till the next SCALARS
+                      // keyword
+                      std::string set = "";
+                      in >> keyword;
+                      for (auto set_cmp : data_sets)
+                        if (keyword == set_cmp)
+                          set = keyword;
+                      if (set == "")
+                        // keep ignoring everything until the next SCALARS
+                        // keyword
+                        continue;
+
+                      // Now we got somewhere. Proceed from here.
+                      // Ignore everything till the end of the line.
+                      // SCALARS MaterialID 1 // this "1" is optional
+                      in.ignore(256, '\n');
+
+                      in >> keyword;
+                      AssertThrow(
+                        keyword == "LOOKUP_TABLE",
+                        ExcMessage(
+                          "While reading VTK file, missing keyword LOOKUP_TABLE"));
 
-      ////////////////////////Processing the CELL_DATA
-      /// section/////////////////////////////
-
-      if (keyword == "CELL_DATA")
-        {
-          unsigned int n_ids;
-          in >> n_ids;
-
-          AssertThrow(
-            n_ids == cells.size() +
-                       (dim == 3 ?
-                          subcelldata.boundary_quads.size() :
-                          (dim == 2 ? subcelldata.boundary_lines.size() : 0)),
-            ExcMessage(
-              "The VTK reader found a CELL_DATA statement "
-              "that lists a total of " +
-              Utilities::int_to_string(n_ids) +
-              " cell data objects, but this needs to "
-              "equal the number of cells (which is " +
-              Utilities::int_to_string(cells.size()) +
-              ") plus the number of quads (" +
-              Utilities::int_to_string(subcelldata.boundary_quads.size()) +
-              " in 3d or the number of lines (" +
-              Utilities::int_to_string(subcelldata.boundary_lines.size()) +
-              ") in 2d."));
-
-
-          std::string linenew;
-          std::string textnew[2];
-          textnew[0] = "SCALARS MaterialID double";
-          textnew[1] = "LOOKUP_TABLE default";
-
-          in.ignore(256, '\n');
-
-          for (unsigned int i = 0; i < 2; i++)
-            {
-              getline(in, linenew);
-              if (i == 0)
-                if (linenew.size() > textnew[0].size())
-                  linenew.resize(textnew[0].size());
-
-              AssertThrow(linenew.compare(textnew[i]) == 0,
-                          ExcMessage(
-                            std::string(
-                              "While reading VTK file, failed to find <") +
-                            textnew[i] + "> section"));
-            }
-
-          // read material ids first for all cells, then for all
-          // faces. the assumption that cells come before all faces
-          // has been verified above via an assertion, so the order
-          // used in the following blocks makes sense
-          for (unsigned int i = 0; i < cells.size(); i++)
-            {
-              double id;
-              in >> id;
-              cells[i].material_id = id;
-            }
-
-          if (dim == 3)
-            {
-              for (unsigned int i = 0; i < subcelldata.boundary_quads.size();
-                   i++)
-                {
-                  double id;
-                  in >> id;
-                  subcelldata.boundary_quads[i].material_id = id;
-                }
-            }
-          else if (dim == 2)
-            {
-              for (unsigned int i = 0; i < subcelldata.boundary_lines.size();
-                   i++)
-                {
-                  double id;
-                  in >> id;
-                  subcelldata.boundary_lines[i].material_id = id;
-                }
-            }
-        }
+                      in >> keyword;
+                      AssertThrow(
+                        keyword == "default",
+                        ExcMessage(
+                          "While reading VTK file, missing keyword default"));
+
+                      // read material or manifold ids first for all cells,
+                      // then for all faces, and finally for all lines. the
+                      // assumption that cells come before all faces and
+                      // lines has been verified above via an assertion, so
+                      // the order used in the following blocks makes sense
+                      for (unsigned int i = 0; i < cells.size(); i++)
+                        {
+                          double id;
+                          in >> id;
+                          if (set == "MaterialID")
+                            cells[i].material_id = id;
+                          else if (set == "ManifoldID")
+                            cells[i].manifold_id = id;
+                          else
+                            Assert(false, ExcInternalError());
+                        }
+
+                      if (dim == 3)
+                        {
+                          for (unsigned int i = 0;
+                               i < subcelldata.boundary_quads.size();
+                               i++)
+                            {
+                              double id;
+                              in >> id;
+                              if (set == "MaterialID")
+                                subcelldata.boundary_quads[i].material_id = id;
+                              else if (set == "ManifoldID")
+                                subcelldata.boundary_quads[i].manifold_id = id;
+                              else
+                                Assert(false, ExcInternalError());
+                            }
+                          for (unsigned int i = 0;
+                               i < subcelldata.boundary_lines.size();
+                               i++)
+                            {
+                              double id;
+                              in >> id;
+                              if (set == "MaterialID")
+                                subcelldata.boundary_lines[i].material_id = id;
+                              else if (set == "ManifoldID")
+                                subcelldata.boundary_lines[i].manifold_id = id;
+                              else
+                                Assert(false, ExcInternalError());
+                            }
+                        }
+                      else if (dim == 2)
+                        {
+                          for (unsigned int i = 0;
+                               i < subcelldata.boundary_lines.size();
+                               i++)
+                            {
+                              double id;
+                              in >> id;
+                              if (set == "MaterialID")
+                                subcelldata.boundary_lines[i].material_id = id;
+                              else if (set == "ManifoldID")
+                                subcelldata.boundary_lines[i].manifold_id = id;
+                              else
+                                Assert(false, ExcInternalError());
+                            }
+                        }
+                    }
+              }
+          }
 
       Assert(subcelldata.check_consistency(dim), ExcInternalError());
 
@@ -414,7 +481,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 }
 
 
-
 template <int dim, int spacedim>
 void
 GridIn<dim, spacedim>::read_unv(std::istream &in)
index 6e1d90b7ed2ed6c6da1c6320f733525f46178d54..71ed09807b86511e776c3d3e17c4e587c7833a0e 100644 (file)
@@ -7,6 +7,8 @@ POINTS 4 float
 0.5 0.5 0 
 CELLS 1 5
 4 0 1 3 2
+CELL_TYPES 1
+9
 CELL_DATA 1
 SCALARS MaterialID double
 LOOKUP_TABLE default

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.