Merge pull request #13035 from bangerth/grid-out

Enable GridOut::write_gmsh() for simplex meshes

doc/news/changes/minor/20211206Bangerth

@@ -0,0 +1,3 @@
+New: GridOut::write_msh() now also works for triangular and tetrahedral meshes.
+<br>
+(Wolfgang Bangerth, 2021/12/06)

include/deal.II/grid/reference_cell.h

@@ -518,6 +518,12 @@ class ReferenceCell
   unsigned int
   vtk_lagrange_type() const;
 
+  /**
+   * Return the GMSH element type code that corresponds to the reference cell.
+   */
+  unsigned int
+  gmsh_element_type() const;
+
   /**
    * @}
    */

source/grid/grid_out.cc

@@ -1068,67 +1068,28 @@ GridOut::write_msh(const Triangulation<dim, spacedim> &tria,
             (msh_flags.write_lines ? n_boundary_lines(tria) : 0))
       << '\n';
 
-  /*
-    elm-type
-    defines the geometrical type of the n-th element:
-    1
-    Line (2 nodes).
-    2
-    Triangle (3 nodes).
-    3
-    Quadrangle (4 nodes).
-    4
-    Tetrahedron (4 nodes).
-    5
-    Hexahedron (8 nodes).
-    6
-    Prism (6 nodes).
-    7
-    Pyramid (5 nodes).
-    8
-    Second order line (3 nodes: 2 associated with the vertices and 1 with the
-    edge).
-    9
-    Second order triangle (6 nodes: 3 associated with the vertices and 3 with
-    the edges). 10 Second order quadrangle (9 nodes: 4 associated with the
-    vertices, 4 with the edges and 1 with the face). 11 Second order tetrahedron
-    (10 nodes: 4 associated with the vertices and 6 with the edges). 12 Second
-    order hexahedron (27 nodes: 8 associated with the vertices, 12 with the
-    edges, 6 with the faces and 1 with the volume). 13 Second order prism (18
-    nodes: 6 associated with the vertices, 9 with the edges and 3 with the
-    quadrangular faces). 14 Second order pyramid (14 nodes: 5 associated with
-    the vertices, 8 with the edges and 1 with the quadrangular face). 15 Point
-    (1 node).
-  */
-  unsigned int elm_type;
-  switch (dim)
-    {
-      case 1:
-        elm_type = 1;
-        break;
-      case 2:
-        elm_type = 3;
-        break;
-      case 3:
-        elm_type = 5;
-        break;
-      default:
-        Assert(false, ExcNotImplemented());
-    }
-
   // write cells. Enumerate cells
   // consecutively, starting with 1
   for (const auto &cell : tria.active_cell_iterators())
     {
-      out << cell->active_cell_index() + 1 << ' ' << elm_type << ' '
+      out << cell->active_cell_index() + 1 << ' '
+          << cell->reference_cell().gmsh_element_type() << ' '
           << cell->material_id() << ' ' << cell->subdomain_id() << ' '
           << cell->n_vertices() << ' ';
 
       // Vertex numbering follows UCD conventions.
 
-      for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
-        out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
-            << ' ';
+      for (const unsigned int vertex : cell->vertex_indices())
+        {
+          if (cell->reference_cell() == ReferenceCells::get_hypercube<dim>())
+            out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) +
+                     1
+                << ' ';
+          else if (cell->reference_cell() == ReferenceCells::get_simplex<dim>())
+            out << cell->vertex_index(vertex) + 1 << ' ';
+          else
+            Assert(false, ExcNotImplemented());
+        }
       out << '\n';
     }
 
@@ -3951,29 +3912,25 @@ GridOut::write_msh_faces(const Triangulation<dim, spacedim> &tria,
   for (const auto &face : tria.active_face_iterators())
     if (face->at_boundary() && (face->boundary_id() != 0))
       {
-        out << current_element_index << ' ';
-        switch (dim)
-          {
-            case 2:
-              out << 1 << ' ';
-              break;
-            case 3:
-              out << 3 << ' ';
-              break;
-            default:
-              Assert(false, ExcNotImplemented());
-          }
+        out << current_element_index << ' '
+            << face->reference_cell().gmsh_element_type() << ' ';
         out << static_cast<unsigned int>(face->boundary_id()) << ' '
             << static_cast<unsigned int>(face->boundary_id()) << ' '
-            << GeometryInfo<dim>::vertices_per_face;
+            << face->n_vertices();
         // note: vertex numbers are 1-base
-        for (unsigned int vertex = 0;
-             vertex < GeometryInfo<dim>::vertices_per_face;
-             ++vertex)
-          out << ' '
-              << face->vertex_index(
-                   GeometryInfo<dim - 1>::ucd_to_deal[vertex]) +
-                   1;
+        for (unsigned int vertex : face->vertex_indices())
+          {
+            if (face->reference_cell() == ReferenceCells::Quadrilateral)
+              out << ' '
+                  << face->vertex_index(
+                       GeometryInfo<dim - 1>::ucd_to_deal[vertex]) +
+                       1;
+            else if ((face->reference_cell() == ReferenceCells::Triangle) ||
+                     (face->reference_cell() == ReferenceCells::Line))
+              out << ' ' << face->vertex_index(vertex) + 1;
+            else
+              Assert(false, ExcInternalError());
+          }
         out << '\n';
 
         ++current_element_index;
@@ -3982,6 +3939,7 @@ GridOut::write_msh_faces(const Triangulation<dim, spacedim> &tria,
 }
 
 
+
 template <int dim, int spacedim>
 unsigned int
 GridOut::write_msh_lines(const Triangulation<dim, spacedim> &tria,
@@ -4004,10 +3962,11 @@ GridOut::write_msh_lines(const Triangulation<dim, spacedim> &tria,
       if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
           (cell->line(l)->user_flag_set() == false))
         {
-          out << next_element_index << " 1 ";
+          out << next_element_index << ' '
+              << ReferenceCells::Line.gmsh_element_type() << ' ';
           out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
               << static_cast<unsigned int>(cell->line(l)->boundary_id())
-              << " 2 ";
+              << " 2 "; // two vertex indices to follow
           // note: vertex numbers are 1-base
           for (unsigned int vertex = 0; vertex < 2; ++vertex)
             out << ' '

source/grid/reference_cell.cc

@@ -458,6 +458,91 @@ ReferenceCell::vtk_lagrange_type() const
 
 
 
+unsigned int
+ReferenceCell::gmsh_element_type() const
+{
+  /*
+    From the GMSH documentation:
+
+    elm-type
+    defines the geometrical type of the n-th element:
+
+    1
+    Line (2 nodes).
+
+    2
+    Triangle (3 nodes).
+
+    3
+    Quadrangle (4 nodes).
+
+    4
+    Tetrahedron (4 nodes).
+
+    5
+    Hexahedron (8 nodes).
+
+    6
+    Prism (6 nodes).
+
+    7
+    Pyramid (5 nodes).
+
+    8
+    Second order line (3 nodes: 2 associated with the vertices and 1 with the
+    edge).
+
+    9
+    Second order triangle (6 nodes: 3 associated with the vertices and 3 with
+    the edges).
+
+    10 Second order quadrangle (9 nodes: 4 associated with the
+    vertices, 4 with the edges and 1 with the face).
+
+    11 Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
+    with the edges).
+
+    12 Second order hexahedron (27 nodes: 8 associated with the vertices, 12
+    with the edges, 6 with the faces and 1 with the volume).
+
+    13 Second order prism (18 nodes: 6 associated with the vertices, 9 with the
+    edges and 3 with the quadrangular faces).
+
+    14 Second order pyramid (14 nodes: 5 associated with the vertices, 8 with
+    the edges and 1 with the quadrangular face).
+
+    15 Point (1 node).
+  */
+
+  if (*this == ReferenceCells::Vertex)
+    return 15;
+  else if (*this == ReferenceCells::Line)
+    return 1;
+  else if (*this == ReferenceCells::Triangle)
+    return 2;
+  else if (*this == ReferenceCells::Quadrilateral)
+    return 3;
+  else if (*this == ReferenceCells::Tetrahedron)
+    return 4;
+  else if (*this == ReferenceCells::Pyramid)
+    return 7;
+  else if (*this == ReferenceCells::Wedge)
+    return 6;
+  else if (*this == ReferenceCells::Hexahedron)
+    return 5;
+  else if (*this == ReferenceCells::Invalid)
+    {
+      Assert(false, ExcNotImplemented());
+      return numbers::invalid_unsigned_int;
+    }
+
+  Assert(false, ExcNotImplemented());
+
+  return numbers::invalid_unsigned_int;
+}
+
+
+
 std::ostream &
 operator<<(std::ostream &out, const ReferenceCell &reference_cell)
 {

tests/simplex/grid_out_msh.cc

@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Create a triangulation with simplices and output it in gmsh format
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+check()
+{
+  Triangulation<dim> triangulation;
+
+  {
+    Triangulation<dim> x;
+    GridGenerator::subdivided_cylinder(x, 4, 1, 2);
+
+    // For the cylinder geometry
+    for (auto cell : x.active_cell_iterators())
+      for (const unsigned int face : cell->face_indices())
+        if (cell->at_boundary(face))
+          {
+            if (std::fabs(cell->face(face)->center()[0] - 0) <
+                0.1 * cell->diameter())
+              cell->face(face)->set_boundary_id(1);
+            else if (std::fabs(cell->face(face)->center()[0] - 4) <
+                     0.1 * cell->diameter())
+              cell->face(face)->set_boundary_id(2);
+          }
+    GridGenerator::convert_hypercube_to_simplex_mesh(x, triangulation);
+  }
+
+  {
+    GridOut go;
+    go.set_flags(GridOutFlags::Msh(true));
+
+    go.write_msh(triangulation, deallog.get_file_stream());
+  }
+
+  deallog << "OK!" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  check<3>();
+}