]> https://gitweb.dealii.org/ - dealii.git/commitdiff
GridIn::read_vtk() for simplex meshes 10748/head
authorPeter Munch <peterrmuench@gmail.com>
Thu, 23 Jul 2020 19:28:26 +0000 (21:28 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 25 Jul 2020 19:01:49 +0000 (21:01 +0200)
doc/news/changes/minor/20200723Munch [new file with mode: 0644]
include/deal.II/grid/grid_in.h
source/grid/grid_in.cc
source/grid/tria.cc
tests/simplex/grid_in_vtk.cc [new file with mode: 0644]
tests/simplex/grid_in_vtk.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/grid_in_vtk/tet.vtk [new file with mode: 0644]
tests/simplex/grid_in_vtk/tri.vtk [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20200723Munch b/doc/news/changes/minor/20200723Munch
new file mode 100644 (file)
index 0000000..801696c
--- /dev/null
@@ -0,0 +1,3 @@
+New: GridIn::read_vtk() now supports triangular and tetrahedral meshes.
+<br>
+(Peter Munch, 2020/07/23)
index b76d8a461e98a4c56c7eaaa302e7a164877a5402..7bb542fc09c0f5da3829e4e672abf2b4308c67ed 100644 (file)
@@ -364,20 +364,22 @@ public:
 
   /**
    * 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.
+   * the following VTK cell types: VTK_HEXAHEDRON (12), VTK_TETRA (10),
+   * VTK_QUAD (9), VTK_TRIANGLE (5), and VTK_LINE (3).
    *
    * 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_HEXAHEDRON/VTK_TETRA cell types
+   * - VTK_QUAD/VTK_TRIANGLE 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_QUAD/VTK_TRIANGLE cell types
    * - VTK_LINE cell types, to specify optional boundary or interior edges
    *
    * In one dimension
@@ -396,6 +398,8 @@ public:
    *
    * The companion GridOut::write_vtk function can be used to write VTK files
    * compatible with this method.
+   *
+   * @ingroup simplex
    */
   void
   read_vtk(std::istream &in);
index 121f18783f6ae8999517e08e1120c8d6a641cb2d..9f3db2828df8b0aa746b2a3aebc1760e9e0c2951 100644 (file)
@@ -176,8 +176,33 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
   unsigned int n_geometric_objects = 0;
   unsigned int n_ints;
 
+  bool is_quad_or_hex_mesh = false;
+  bool is_tria_or_tet_mesh = false;
+
   if (keyword == "CELLS")
     {
+      // jump to the `CELL_TYPES` section and read in cell types
+      std::vector<unsigned int> cell_types;
+      {
+        std::streampos oldpos = in.tellg();
+
+
+        while (in >> keyword)
+          if (keyword == "CELL_TYPES")
+            {
+              in >> n_ints;
+
+              cell_types.resize(n_ints);
+
+              for (unsigned int i = 0; i < n_ints; ++i)
+                in >> cell_types[i];
+
+              break;
+            }
+
+        in.seekg(oldpos);
+      }
+
       in >> n_geometric_objects;
       in >> n_ints; // Ignore this, since we don't need it.
 
@@ -188,15 +213,20 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               unsigned int type;
               in >> type;
 
-              if (type == 8)
+              if (cell_types[count] == 10 || cell_types[count] == 12)
                 {
+                  if (cell_types[count] == 10)
+                    is_tria_or_tet_mesh = true;
+                  if (cell_types[count] == 12)
+                    is_quad_or_hex_mesh = true;
+
                   // we assume that the file contains first all cells,
                   // and only then any faces or lines
                   AssertThrow(subcelldata.boundary_quads.size() == 0 &&
                                 subcelldata.boundary_lines.size() == 0,
                               ExcNotImplemented());
 
-                  cells.emplace_back();
+                  cells.emplace_back(type);
 
                   for (unsigned int j = 0; j < type; j++) // loop to feed data
                     in >> cells.back().vertices[j];
@@ -204,14 +234,19 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                   cells.back().material_id = 0;
                 }
 
-              else if (type == 4)
+              else if (cell_types[count] == 5 || cell_types[count] == 9)
                 {
+                  if (cell_types[count] == 5)
+                    is_tria_or_tet_mesh = true;
+                  if (cell_types[count] == 9)
+                    is_quad_or_hex_mesh = true;
+
                   // 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();
+                  subcelldata.boundary_quads.emplace_back(type);
 
                   for (unsigned int j = 0; j < type;
                        j++) // loop to feed the data to the boundary
@@ -219,9 +254,9 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 
                   subcelldata.boundary_quads.back().material_id = 0;
                 }
-              else if (type == 2)
+              else if (cell_types[count] == 3)
                 {
-                  subcelldata.boundary_lines.emplace_back();
+                  subcelldata.boundary_lines.emplace_back(type);
 
                   for (unsigned int j = 0; j < type;
                        j++) // loop to feed the data to the boundary
@@ -237,7 +272,6 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                     "While reading VTK file, unknown file type encountered"));
             }
         }
-
       else if (dim == 2)
         {
           for (unsigned int count = 0; count < n_geometric_objects; count++)
@@ -245,14 +279,19 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               unsigned int type;
               in >> type;
 
-              if (type == 4)
+              if (cell_types[count] == 5 || cell_types[count] == 9)
                 {
                   // we assume that the file contains first all cells,
                   // and only then any faces
                   AssertThrow(subcelldata.boundary_lines.size() == 0,
                               ExcNotImplemented());
 
-                  cells.emplace_back();
+                  if (cell_types[count] == 5)
+                    is_tria_or_tet_mesh = true;
+                  if (cell_types[count] == 9)
+                    is_quad_or_hex_mesh = true;
+
+                  cells.emplace_back(type);
 
                   for (unsigned int j = 0; j < type; j++) // loop to feed data
                     in >> cells.back().vertices[j];
@@ -260,11 +299,11 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
                   cells.back().material_id = 0;
                 }
 
-              else if (type == 2)
+              else if (cell_types[count] == 3)
                 {
                   // If this is encountered, the pointer comes out of the loop
                   // and starts processing boundaries.
-                  subcelldata.boundary_lines.emplace_back();
+                  subcelldata.boundary_lines.emplace_back(type);
 
                   for (unsigned int j = 0; j < type;
                        j++) // loop to feed the data to the boundary
@@ -290,10 +329,10 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
               in >> type;
 
               AssertThrow(
-                type == 2,
+                cell_types[count] == 3 && type == 2,
                 ExcMessage(
                   "While reading VTK file, unknown cell type encountered"));
-              cells.emplace_back();
+              cells.emplace_back(type);
 
               for (unsigned int j = 0; j < type; j++) // loop to feed data
                 in >> cells.back().vertices[j];
@@ -476,16 +515,36 @@ GridIn<dim, spacedim>::read_vtk(std::istream &in)
 
       Assert(subcelldata.check_consistency(dim), ExcInternalError());
 
-      GridTools::delete_unused_vertices(vertices, cells, subcelldata);
 
-      if (dim == spacedim)
-        GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(
-          vertices, cells);
+      // make sure that only either simplex or hypercube cells are available
+      //
+      // TODO: the functions below (GridTools::delete_unused_vertices(),
+      // GridTools::invert_all_cells_of_negative_grid(),
+      // GridReordering::reorder_cells(),
+      // Triangulation::create_triangulation_compatibility()) need to be
+      // revisited for simplex meshes
+      AssertThrow(dim == 1 || (is_tria_or_tet_mesh ^ is_quad_or_hex_mesh),
+                  ExcNotImplemented());
+
+      if (dim == 1 || is_quad_or_hex_mesh)
+        {
+          GridTools::delete_unused_vertices(vertices, cells, subcelldata);
+
+          if (dim == spacedim)
+            GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(
+              vertices, cells);
 
-      GridReordering<dim, spacedim>::reorder_cells(cells);
-      tria->create_triangulation_compatibility(vertices, cells, subcelldata);
+          GridReordering<dim, spacedim>::reorder_cells(cells);
+          tria->create_triangulation_compatibility(vertices,
+                                                   cells,
+                                                   subcelldata);
 
-      return;
+          return;
+        }
+      else
+        {
+          tria->create_triangulation(vertices, cells, subcelldata);
+        }
     }
   else
     AssertThrow(false,
index c4e207e428508f4411a944b6b2e90995cdf01869..3805f69cb1525afb0a4403c8b55e1c95c8ffbe11 100644 (file)
@@ -9743,7 +9743,7 @@ Triangulation<dim, spacedim>::create_triangulation(
                cell != this_round.end();
                ++cell)
             {
-              for (const unsigned int i : GeometryInfo<dim>::face_indices())
+              for (const unsigned int i : (*cell)->face_indices())
                 {
                   if (!((*cell)->face(i)->at_boundary()))
                     {
diff --git a/tests/simplex/grid_in_vtk.cc b/tests/simplex/grid_in_vtk.cc
new file mode 100644 (file)
index 0000000..4b32c6d
--- /dev/null
@@ -0,0 +1,73 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Read a file in the MSH format with (linear) triangular and tetrahedral
+// elements created by the GMSH program.
+
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+check_file(const std::string &file_name)
+{
+  Triangulation<dim, spacedim> tria;
+
+  GridIn<dim, spacedim> grid_in;
+  grid_in.attach_triangulation(tria);
+  std::ifstream input_file(file_name);
+  grid_in.read_vtk(input_file);
+
+  GridOut grid_out;
+#if false
+  std::ofstream out("mesh.out.vtk");
+  grid_out.write_vtk(tria, out);
+#else
+  grid_out.write_vtk(tria, deallog.get_file_stream());
+#endif
+
+  deallog << "OK!" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  // TRIANGULAR ELEMENTS
+  // dim = spacedim = 2
+  deallog.push("triangluar_elements_dim2_spacedim2: ");
+  check_file<2, 2>(std::string(SOURCE_DIR "/grid_in_vtk/tri.vtk"));
+  deallog.pop();
+
+  // dim = 2, spacedim = 3
+  deallog.push("triangluar_elements_dim2_spacedim3: ");
+  check_file<2, 3>(std::string(SOURCE_DIR "/grid_in_vtk/tri.vtk"));
+  deallog.pop();
+
+
+
+  // TETRAHEDRAL ELEMENTS
+  // dim = spacedim = 3
+  deallog.push("tetrahedral_elements_dim3_spacedim3: ");
+  check_file<3, 3>(std::string(SOURCE_DIR "/grid_in_vtk/tet.vtk"));
+  deallog.pop();
+}
diff --git a/tests/simplex/grid_in_vtk.with_simplex_support=on.output b/tests/simplex/grid_in_vtk.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..446eb9c
--- /dev/null
@@ -0,0 +1,135 @@
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 5 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+1.00000 1.00000 0
+0.00000 1.00000 0
+0.500000 0.500000 0
+
+CELLS 4 16
+3 0 1 4
+3 3 0 4
+3 1 2 4
+3 2 3 4
+
+CELL_TYPES 4
+5 5 5 5 
+
+
+
+CELL_DATA 4
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+
+
+DEAL:triangluar_elements_dim2_spacedim2: ::OK!
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 5 double
+0.00000 0.00000 0.00000
+1.00000 0.00000 0.00000
+1.00000 1.00000 0.00000
+0.00000 1.00000 0.00000
+0.500000 0.500000 0.00000
+
+CELLS 4 16
+3 0 1 4
+3 3 0 4
+3 1 2 4
+3 2 3 4
+
+CELL_TYPES 4
+5 5 5 5 
+
+
+
+CELL_DATA 4
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+
+
+DEAL:triangluar_elements_dim2_spacedim3: ::OK!
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 14 double
+0.00000 0.00000 1.00000
+0.00000 0.00000 0.00000
+0.00000 1.00000 1.00000
+0.00000 1.00000 0.00000
+1.00000 0.00000 1.00000
+1.00000 0.00000 0.00000
+1.00000 1.00000 1.00000
+1.00000 1.00000 0.00000
+0.00000 0.500000 0.500000
+1.00000 0.500000 0.500000
+0.500000 0.00000 0.500000
+0.500000 1.00000 0.500000
+0.500000 0.500000 0.00000
+0.500000 0.500000 1.00000
+
+CELLS 24 120
+4 13 10 11 9
+4 9 10 11 12
+4 11 8 10 13
+4 8 11 10 12
+4 4 13 6 9
+4 3 8 11 2
+4 1 8 0 10
+4 0 8 2 13
+4 4 10 0 13
+4 9 4 10 5
+4 11 2 13 6
+4 9 11 6 7
+4 1 3 8 12
+4 11 3 7 12
+4 12 7 9 5
+4 12 10 1 5
+4 13 4 10 9
+4 8 11 2 13
+4 6 13 11 9
+4 8 0 10 13
+4 10 1 8 12
+4 11 8 3 12
+4 7 9 11 12
+4 12 9 10 5
+
+CELL_TYPES 24
+10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 
+
+
+
+CELL_DATA 24
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+
+
+DEAL:tetrahedral_elements_dim3_spacedim3: ::OK!
diff --git a/tests/simplex/grid_in_vtk/tet.vtk b/tests/simplex/grid_in_vtk/tet.vtk
new file mode 100644 (file)
index 0000000..661934c
--- /dev/null
@@ -0,0 +1,99 @@
+# vtk DataFile Version 3.0
+tet, Created by Gmsh
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 14 double
+0 0 1
+0 0 0
+0 1 1
+0 1 0
+1 0 1
+1 0 0
+1 1 1
+1 1 0
+0 0.5 0.5
+1 0.5 0.5
+0.5 0 0.5
+0.5 1 0.5
+0.5 0.5 0
+0.5 0.5 1
+
+CELLS 24 120
+4 13 10 11 9
+4 9 10 11 12
+4 11 8 10 13
+4 8 11 10 12
+4 4 13 6 9
+4 3 8 11 2
+4 1 8 0 10
+4 0 8 2 13
+4 4 10 0 13
+4 9 4 10 5
+4 11 2 13 6
+4 9 11 6 7
+4 1 3 8 12
+4 11 3 7 12
+4 12 7 9 5
+4 12 10 1 5
+4 13 4 10 9
+4 8 11 2 13
+4 6 13 11 9
+4 8 0 10 13
+4 10 1 8 12
+4 11 8 3 12
+4 7 9 11 12
+4 12 9 10 5
+
+CELL_TYPES 24
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+10
+
+CELL_DATA 24
+SCALARS CellEntityIds int 1
+LOOKUP_TABLE default
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
diff --git a/tests/simplex/grid_in_vtk/tri.vtk b/tests/simplex/grid_in_vtk/tri.vtk
new file mode 100644 (file)
index 0000000..672448d
--- /dev/null
@@ -0,0 +1,30 @@
+# vtk DataFile Version 3.0
+tri, Created by Gmsh
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 5 double
+0 0 0
+1 0 0
+1 1 0
+0 1 0
+0.5 0.5 0
+
+CELLS 4 16
+3 0 1 4
+3 3 0 4
+3 1 2 4
+3 2 3 4
+
+CELL_TYPES 4
+5
+5
+5
+5
+
+CELL_DATA 4
+SCALARS CellEntityIds int 1
+LOOKUP_TABLE default
+1
+1
+1
+1

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.