GridIn::read_msh() now supports triangular and tetrahedral elements

diff --git a/doc/news/changes/minor/20200720DanielPaukner b/doc/news/changes/minor/20200720DanielPaukner
new file mode 100644 (file)
index 0000000..8a6e86c
--- /dev/null
+++ b/doc/news/changes/minor/20200720DanielPaukner
@@ -0,0 +1,3 @@
+New: GridIn::read_msh() now supports triangular and tetrahedral meshes.
+<br>
+(Daniel Paukner, 2020/07/20)

diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h
index b76d8a461e98a4c56c7eaaa302e7a164877a5402..74a6188e13f637625e51d6077edfa7140bb4c138 100644 (file)
--- a/include/deal.II/grid/grid_in.h
+++ b/include/deal.II/grid/grid_in.h
@@ -504,6 +504,8 @@ public:
    * @note The input function of deal.II does not distinguish between newline
    * and other whitespace. Therefore, deal.II will be able to read files in a
    * slightly more general format than %Gmsh.
+   *
+   * @ingroup simplex
    */
   void
   read_msh(std::istream &in);

diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index 121f18783f6ae8999517e08e1120c8d6a641cb2d..0fe7e290983e6b9ef00b736f63d28aed690f6089 100644 (file)
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -1731,6 +1731,8 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
   std::vector<CellData<dim>>                 cells;
   SubCellData                                subcelldata;
   std::map<unsigned int, types::boundary_id> boundary_ids_1d;
+  bool                                       is_hex_mesh = false;
+  bool                                       is_tet_mesh = false;
 
   {
     unsigned int global_cell = 0;
@@ -1818,14 +1820,33 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                 for (unsigned int i = 1; i < n_tags; ++i)
                   in >> dummy;
 
-                nod_num = GeometryInfo<dim>::vertices_per_cell;
+                if (cell_type == 1) // line
+                  nod_num = 2;
+                else if (cell_type == 2) // tri
+                  nod_num = 3;
+                else if (cell_type == 3) // quad
+                  nod_num = 4;
+                else if (cell_type == 4) // tet
+                  nod_num = 4;
+                else if (cell_type == 5) // hex
+                  nod_num = 8;
               }
             else
               {
                 // ignore tag
                 int tag;
                 in >> tag;
-                nod_num = GeometryInfo<dim>::vertices_per_cell;
+
+                if (cell_type == 1) // line
+                  nod_num = 2;
+                else if (cell_type == 2) // tri
+                  nod_num = 3;
+                else if (cell_type == 3) // quad
+                  nod_num = 4;
+                else if (cell_type == 4) // tet
+                  nod_num = 4;
+                else if (cell_type == 5) // hex
+                  nod_num = 8;
               }
 
 
@@ -1834,9 +1855,15 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                      `1'
                      Line (2 nodes, 1 edge).
 
+                     `2'
+                     Triangle (3 nodes, 3 edges).
+
                      `3'
                      Quadrangle (4 nodes, 4 edges).
 
+                     `4'
+                     Tetrahedron (4 nodes, 6 edges, 6 faces).
+
                      `5'
                      Hexahedron (8 nodes, 12 edges, 6 faces).
 
@@ -1845,18 +1872,45 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
             */
 
             if (((cell_type == 1) && (dim == 1)) ||
+                ((cell_type == 2) && (dim == 2)) ||
                 ((cell_type == 3) && (dim == 2)) ||
+                ((cell_type == 4) && (dim == 3)) ||
                 ((cell_type == 5) && (dim == 3)))
               // found a cell
               {
-                AssertThrow(nod_num == GeometryInfo<dim>::vertices_per_cell,
+                unsigned int vertices_per_cell = 0;
+                if (cell_type == 1) // line
+                  vertices_per_cell = 2;
+                else if (cell_type == 2) // tri
+                  {
+                    vertices_per_cell = 3;
+                    is_tet_mesh       = true;
+                  }
+                else if (cell_type == 3) // quad
+                  {
+                    vertices_per_cell = 4;
+                    is_hex_mesh       = true;
+                  }
+                else if (cell_type == 4) // tet
+                  {
+                    vertices_per_cell = 4;
+                    is_tet_mesh       = true;
+                  }
+                else if (cell_type == 5) // hex
+                  {
+                    vertices_per_cell = 8;
+                    is_hex_mesh       = true;
+                  }
+
+                AssertThrow(nod_num == vertices_per_cell,
                             ExcMessage(
                               "Number of nodes does not coincide with the "
                               "number required for this object"));
 
                 // allocate and read indices
                 cells.emplace_back();
-                for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+                cells.back().vertices.resize(vertices_per_cell);
+                for (unsigned int i = 0; i < vertices_per_cell; ++i)
                   in >> cells.back().vertices[i];
 
                 // to make sure that the cast won't fail
@@ -1874,7 +1928,7 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
 
                 // transform from ucd to
                 // consecutive numbering
-                for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+                for (unsigned int i = 0; i < vertices_per_cell; ++i)
                   {
                     AssertThrow(
                       vertex_indices.find(cells.back().vertices[i]) !=
@@ -1926,14 +1980,30 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                       vertex = numbers::invalid_unsigned_int;
                     }
               }
-            else if ((cell_type == 3) && (dim == 3))
+            else if ((cell_type == 2 || cell_type == 3) && (dim == 3))
               // boundary info
               {
+                unsigned int vertices_per_cell = 0;
+                // check cell type
+                if (cell_type == 2) // tri
+                  {
+                    vertices_per_cell = 3;
+                    is_tet_mesh       = true;
+                  }
+                else if (cell_type == 3) // quad
+                  {
+                    vertices_per_cell = 4;
+                    is_hex_mesh       = true;
+                  }
+
                 subcelldata.boundary_quads.emplace_back();
-                in >> subcelldata.boundary_quads.back().vertices[0] >>
-                  subcelldata.boundary_quads.back().vertices[1] >>
-                  subcelldata.boundary_quads.back().vertices[2] >>
-                  subcelldata.boundary_quads.back().vertices[3];
+
+                // resize vertices
+                subcelldata.boundary_quads.back().vertices.resize(
+                  vertices_per_cell);
+                // for loop
+                for (unsigned int i = 0; i < vertices_per_cell; ++i)
+                  in >> subcelldata.boundary_quads.back().vertices[i];
 
                 // to make sure that the cast won't fail
                 Assert(material_id <=
@@ -1987,22 +2057,7 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
                   boundary_ids_1d[vertex_indices[node_index]] = material_id;
               }
             else
-              // cannot read this, so throw
-              // an exception. treat
-              // triangles and tetrahedra
-              // specially since this
-              // deserves a more explicit
-              // error message
               {
-                AssertThrow(cell_type != 2,
-                            ExcMessage("Found triangles while reading a file "
-                                       "in gmsh format. deal.II does not "
-                                       "support triangles"));
-                AssertThrow(cell_type != 11,
-                            ExcMessage("Found tetrahedra while reading a file "
-                                       "in gmsh format. deal.II does not "
-                                       "support tetrahedra"));
-
                 AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type));
               }
           }
@@ -2020,19 +2075,33 @@ GridIn<dim, spacedim>::read_msh(std::istream &in)
 
   AssertThrow(in, ExcIO());
 
-  // check that we actually read some
-  // cells.
+  // check that we actually read some cells.
   AssertThrow(cells.size() > 0, ExcGmshNoCellInformation());
 
-  // do some clean-up on
-  // vertices...
-  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
-  // ... and cells
-  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);
+  // 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_tet_mesh ^ is_hex_mesh), ExcNotImplemented());
+
+  if (dim == 1 || is_hex_mesh)
+    {
+      // do some clean-up on vertices...
+      GridTools::delete_unused_vertices(vertices, cells, subcelldata);
+      // ... and cells
+      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);
+    }
+  else
+    {
+      tria->create_triangulation(vertices, cells, subcelldata);
+    }
 
   // in 1d, we also have to attach boundary ids to vertices, which does not
   // currently work through the call above

diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index eb0b44e86cc67b03dac000c298549cb8164ba6ed..1134124599c030e20ad0db88b4fafab7c2f33b55 100644 (file)
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -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_msh.cc b/tests/simplex/grid_in_msh.cc
new file mode 100644 (file)
index 0000000..397ddab
--- /dev/null
+++ b/tests/simplex/grid_in_msh.cc
@@ -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, const unsigned int vertices_per_cell)
+{
+  Triangulation<dim, spacedim> tria;
+
+  GridIn<dim, spacedim> grid_in;
+  grid_in.attach_triangulation(tria);
+  std::ifstream input_file(file_name);
+  grid_in.read_msh(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_msh/tri.msh"), 3);
+  deallog.pop();
+
+  // dim = 2, spacedim = 3
+  deallog.push("triangluar_elements_dim2_spacedim3: ");
+  check_file<2, 3>(std::string(SOURCE_DIR "/grid_in_msh/tri.msh"), 3);
+  deallog.pop();
+
+
+
+  // TETRAHEDRAL ELEMENTS
+  // dim = spacedim = 3
+  deallog.push("tetrahedral_elements_dim3_spacedim3: ");
+  check_file<3, 3>(std::string(SOURCE_DIR "/grid_in_msh/tet.msh"), 4);
+  deallog.pop();
+}

diff --git a/tests/simplex/grid_in_msh.with_simplex_support=on.output b/tests/simplex/grid_in_msh.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..b1fa345
--- /dev/null
+++ b/tests/simplex/grid_in_msh.with_simplex_support=on.output
@@ -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
+1 1 1 1 
+
+
+
+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
+1 1 1 1 
+
+
+
+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
+1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
+
+
+
+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_msh/tet.msh b/tests/simplex/grid_in_msh/tet.msh
new file mode 100644 (file)
index 0000000..16e73d3
--- /dev/null
+++ b/tests/simplex/grid_in_msh/tet.msh
@@ -0,0 +1,107 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 1 0 
+2 0 0 0 0 
+3 0 1 1 0 
+4 0 1 0 0 
+5 1 0 1 0 
+6 1 0 0 0 
+7 1 1 1 0 
+8 1 1 0 0 
+1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1 
+2 -1e-07 -9.999999994736442e-08 0.9999999000000001 1e-07 1.0000001 1.0000001 0 2 1 -3 
+3 -1e-07 0.9999999000000001 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 2 4 -3 
+4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 2 -4 
+5 0.9999999000000001 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 2 6 -5 
+6 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 5 -7 
+7 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 2 8 -7 
+8 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 6 -8 
+9 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 2 -6 
+10 -9.999999994736442e-08 -1e-07 0.9999999000000001 1.0000001 1e-07 1.0000001 0 2 1 -5 
+11 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 4 -8 
+12 -9.999999994736442e-08 0.9999999000000001 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 3 -7 
+1 -1e-07 -9.999999994736442e-08 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 4 1 2 -3 -4 
+2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 5 6 -7 -8 
+3 -9.999999994736442e-08 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 4 9 5 -10 -1 
+4 -9.999999994736442e-08 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 11 7 -12 -3 
+5 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 4 11 -8 -9 
+6 -9.999999994736442e-08 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 4 2 12 -6 -10 
+1 -9.999999994736442e-08 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 1 1 6 1 2 3 4 5 6 
+$EndEntities
+$Nodes
+15 14 1 14
+0 1 0 1
+1
+0 0 1
+0 2 0 1
+2
+0 0 0
+0 3 0 1
+3
+0 1 1
+0 4 0 1
+4
+0 1 0
+0 5 0 1
+5
+1 0 1
+0 6 0 1
+6
+1 0 0
+0 7 0 1
+7
+1 1 1
+0 8 0 1
+8
+1 1 0
+2 1 0 1
+9
+0 0.5 0.5
+2 2 0 1
+10
+1 0.5 0.5
+2 3 0 1
+11
+0.5 0 0.5
+2 4 0 1
+12
+0.5 1 0.5
+2 5 0 1
+13
+0.5 0.5 0
+2 6 0 1
+14
+0.5 0.5 1
+3 1 0 0
+$EndNodes
+$Elements
+1 24 1 24
+3 1 4 24
+1 14 11 12 10 
+2 10 11 12 13 
+3 12 9 11 14 
+4 9 12 11 13 
+5 5 14 7 10 
+6 4 9 12 3 
+7 2 9 1 11 
+8 1 9 3 14 
+9 5 11 1 14 
+10 10 5 11 6 
+11 12 3 14 7 
+12 10 12 7 8 
+13 2 4 9 13 
+14 12 4 8 13 
+15 13 8 10 6 
+16 13 11 2 6 
+17 14 5 11 10 
+18 9 12 3 14 
+19 7 14 12 10 
+20 9 1 11 14 
+21 11 2 9 13 
+22 12 9 4 13 
+23 8 10 12 13 
+24 13 10 11 6 
+$EndElements

diff --git a/tests/simplex/grid_in_msh/tri.msh b/tests/simplex/grid_in_msh/tri.msh
new file mode 100644 (file)
index 0000000..822926e
--- /dev/null
+++ b/tests/simplex/grid_in_msh/tri.msh
@@ -0,0 +1,41 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+4 4 1 0
+1 0 0 0 0 
+2 1 0 0 0 
+3 1 1 0 0 
+4 0 1 0 0 
+1 0 0 0 1 0 0 0 2 1 -2 
+2 1 0 0 1 1 0 0 2 2 -3 
+3 0 1 0 1 1 0 0 2 3 -4 
+4 0 0 0 0 1 0 0 2 4 -1 
+1 0 0 0 1 1 0 1 1 4 3 4 1 2 
+$EndEntities
+$Nodes
+5 5 1 5
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+1 0 0
+0 3 0 1
+3
+1 1 0
+0 4 0 1
+4
+0 1 0
+2 1 0 1
+5
+0.5 0.5 0
+$EndNodes
+$Elements
+1 4 1 4
+2 1 2 4
+1 1 2 5 
+2 4 1 5 
+3 2 3 5 
+4 3 4 5 
+$EndElements