]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add Luca Heltai's reader for gmsh format.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 21 Apr 2004 21:28:06 +0000 (21:28 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 21 Apr 2004 21:28:06 +0000 (21:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@9080 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_in.h
deal.II/deal.II/source/grid/grid_in.cc
tests/bits/Makefile
tests/bits/grid_in_msh_01.2d.msh [new file with mode: 0644]
tests/bits/grid_in_msh_01.3d.msh [new file with mode: 0644]
tests/bits/grid_in_msh_01.cc [new file with mode: 0644]

index 35199ae262884c587856a14295dcc8970b6c7eb9..b13588b715e8ff9fbd0ab4425227e7c7ca8ccc26 100644 (file)
@@ -92,6 +92,11 @@ class SubCellData;
  * code. We don't have an exact specification of the format, but the reader
  * can read in several example files. If the reader does not grok your files,
  * it should be fairly simple to extend it.
+ *
+ * <li> <tt>Gmsh mesh</tt> format: this format is used by the @p GMSH
+ * mesh generator (see http://www.geuz.org/gmsh/ ). The documentation
+ * in the @p GMSH manual explains how to generate meshes compatible with
+ * the deal.II library (i.e. quads rather than triangles).
  * </ul>
  *
  *
@@ -148,7 +153,7 @@ class SubCellData;
  * that class if you experience unexpected problems when reading grids
  * through this class.
  *
- * @author Wolfgang Bangerth, 1998, 2000
+ * @author Wolfgang Bangerth, 1998, 2000, Luca Heltai, 2004
  */
 template <int dim>
 class GridIn
@@ -185,6 +190,11 @@ class GridIn
                                      */
     void read_xda (std::istream &);
     
+                                    /**
+                                     * Read grid data from an msh file.
+                                     */
+    void read_msh (std::istream &);
+
                                     /**
                                      * Exception
                                      */
@@ -223,7 +233,22 @@ class GridIn
                    int,
                    << "The specified dimension " << arg1
                    << " is not the same as that of the triangulation to be created.");
-    
+                                   
+    DeclException1 (ExcInvalidGMSHInput,
+                   std::string,
+                   << "The string <" << arg1 << "> is not recognized at the present"
+                   << " position of a Gmsh Mesh file.");
+
+    DeclException1 (ExcGmshUnsupportedGeometry,
+                   int,
+                   << "The Element Identifier <" << arg1 << "> is not "
+                   << "supported in the Deal.II Library.\n"
+                   << "Supported elements are: \n"
+                   << "ELM-TYPE\n"
+                   << "1 Line (2 nodes, 1 edge).\n"
+                   << "3 Quadrangle (4 nodes, 4 edges).\n"
+                   << "5 Hexahedron (8 nodes, 12 edges, 6 faces).\n"
+                   << "15 Point (1 node).");
   private:
                                     /**
                                      * Store address of the triangulation to
index f65b079cf5192b00ec53690e356449d7bf5ebe9f..4af96c5da7dd5c9a2162c80cdb33a5ae7d5e2ba1 100644 (file)
@@ -527,6 +527,207 @@ void GridIn<3>::read_xda (std::istream &in)
 
 
 
+template <int dim>
+void GridIn<dim>::read_msh (std::istream &in)
+{
+  Assert (tria != 0, ExcNoTriangulationSelected());
+  AssertThrow (in, ExcIO());
+
+  unsigned int n_vertices;
+  unsigned int n_cells;
+  unsigned int dummy;
+  std::string line;
+  
+  getline(in, line);
+
+  AssertThrow (line=="$NOD",
+              ExcInvalidGMSHInput(line));
+
+  in >> n_vertices;
+
+  std::vector<Point<dim> >     vertices (n_vertices);
+                                  // set up mapping between numbering
+                                  // in msh-file (nod) and in the
+                                  // vertices vector
+  std::map<int,int> vertex_indices;
+  
+  for (unsigned int vertex=0; vertex<n_vertices; ++vertex) 
+    {
+      int vertex_number;
+      double x[3];
+
+                                      // read vertex
+      in >> vertex_number
+        >> x[0] >> x[1] >> x[2];
+     
+      for (unsigned int d=0; d<dim; ++d)
+       vertices[vertex](d) = x[d];
+                                      // store mapping;
+      vertex_indices[vertex_number] = vertex;
+    };
+  
+                                  // This is needed to flush the last
+                                  // new line
+  getline (in, line);
+  
+                                  // Now read in next bit
+  getline (in, line);
+  AssertThrow (line=="$ENDNOD",
+              ExcInvalidGMSHInput(line));
+
+  
+  getline (in, line);
+  AssertThrow (line=="$ELM",
+              ExcInvalidGMSHInput(line));
+
+  in >> n_cells;
+                                  // set up array of cells
+  std::vector<CellData<dim> > cells;
+  SubCellData                 subcelldata;
+
+  for (unsigned int cell=0; cell<n_cells; ++cell) 
+    {
+                                      // note that since in the input
+                                      // file we found the number of
+                                      // cells at the top, there
+                                      // should still be input here,
+                                      // so check this:
+      AssertThrow (in, ExcIO());
+      
+/*  
+    $ENDNOD
+    $ELM
+    NUMBER-OF-ELEMENTS
+    ELM-NUMBER ELM-TYPE REG-PHYS REG-ELEM NUMBER-OF-NODES NODE-NUMBER-LIST
+    ...
+    $ENDELM
+*/
+      
+      unsigned int cell_type;
+      unsigned int material_id;
+      unsigned int nod_num;
+      
+      in >> dummy          // ELM-NUMBER
+        >> cell_type      // ELM-TYPE
+        >> material_id    // REG-PHYS
+        >> dummy          // reg_elm
+        >> nod_num;
+      
+/*       `ELM-TYPE'
+        defines the geometrical type of the N-th element:
+        `1'
+        Line (2 nodes, 1 edge).
+                                                                                                
+        `3'
+        Quadrangle (4 nodes, 4 edges).
+                                                                                                
+        `5'
+        Hexahedron (8 nodes, 12 edges, 6 faces).
+                                                                                    
+        `15'
+        Point (1 node).
+*/
+      
+      if (((cell_type == 1) && (dim == 1)) ||
+         ((cell_type == 3) && (dim == 2)) ||
+         ((cell_type == 5) && (dim == 3)))
+                                        // found a cell
+       {
+                                          // allocate and read indices
+         cells.push_back (CellData<dim>());
+         for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+           in >> cells.back().vertices[i];
+         cells.back().material_id = material_id;
+
+                                          // transform from ucd to
+                                          // consecutive numbering
+         for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+           if (vertex_indices.find (cells.back().vertices[i]) != vertex_indices.end())
+                                              // vertex with this index exists
+             cells.back().vertices[i] = vertex_indices[cells.back().vertices[i]];
+           else 
+             {
+                                                // no such vertex index
+               AssertThrow (false, ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
+               cells.back().vertices[i] = -1;
+             };
+       }
+      else
+       if ((cell_type == 1) && ((dim == 2) || (dim == 3)))
+                                          // boundary info
+         {
+           subcelldata.boundary_lines.push_back (CellData<1>());
+           in >> subcelldata.boundary_lines.back().vertices[0]
+              >> subcelldata.boundary_lines.back().vertices[1];
+           subcelldata.boundary_lines.back().material_id = material_id;
+
+                                            // transform from ucd to
+                                            // consecutive numbering
+           for (unsigned int i=0; i<2; ++i)
+             if (vertex_indices.find (subcelldata.boundary_lines.back().vertices[i]) !=
+                 vertex_indices.end())
+                                                // vertex with this index exists
+               subcelldata.boundary_lines.back().vertices[i]
+                 = vertex_indices[subcelldata.boundary_lines.back().vertices[i]];
+             else 
+               {
+                                                  // no such vertex index
+                 AssertThrow (false,
+                              ExcInvalidVertexIndex(cell,
+                                                    subcelldata.boundary_lines.back().vertices[i]));
+                 subcelldata.boundary_lines.back().vertices[i] = -1;
+               };
+         }
+       else
+         if ((cell_type == 3) && (dim == 3))
+                                            // boundary info
+           {
+             subcelldata.boundary_quads.push_back (CellData<2>());
+             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];
+             
+             subcelldata.boundary_quads.back().material_id = material_id;
+             
+                                              // transform from gmsh to
+                                              // consecutive numbering
+             for (unsigned int i=0; i<4; ++i)
+               if (vertex_indices.find (subcelldata.boundary_quads.back().vertices[i]) !=
+                   vertex_indices.end())
+                                                  // vertex with this index exists
+                 subcelldata.boundary_quads.back().vertices[i]
+                   = vertex_indices[subcelldata.boundary_quads.back().vertices[i]];
+               else
+                 {
+                                                    // no such vertex index
+                   Assert (false,
+                           ExcInvalidVertexIndex(cell,
+                                                 subcelldata.boundary_quads.back().vertices[i]));
+                   subcelldata.boundary_quads.back().vertices[i] = -1;
+                 };
+             
+           }
+         else
+                                            // cannot read this
+           AssertThrow (false, ExcGmshUnsupportedGeometry(cell_type));
+    };
+
+  
+                                  // check that no forbidden arrays are used
+  Assert (subcelldata.check_consistency(dim), ExcInternalError());
+
+  AssertThrow (in, ExcIO());
+
+                                  // do some clean-up on vertices...
+  delete_unused_vertices (vertices, cells, subcelldata);
+                                  // ... and cells
+  GridReordering<dim>::reorder_cells (cells);
+  tria->create_triangulation (vertices, cells, subcelldata);
+}
+
+
+
 
 template <int dim>
 void GridIn<dim>::skip_empty_lines (std::istream &in)
index 1726bb96c3f441b6b0c946c92b864a12404ec251..bbc4bfd050e362a3a132ffbe8ae0baa5f6793577 100644 (file)
@@ -58,7 +58,8 @@ tests_x = anna_? \
          dof_constraints_* \
          apply_boundary_values_* \
          error_estimator_* \
-         cylinder_shell_*
+         cylinder_shell_* \
+         grid_in_*
 
 ifeq ($(USE_CONTRIB_PETSC),yes)
   tests_x += petsc_*
diff --git a/tests/bits/grid_in_msh_01.2d.msh b/tests/bits/grid_in_msh_01.2d.msh
new file mode 100644 (file)
index 0000000..ef628fb
--- /dev/null
@@ -0,0 +1,15 @@
+$NOD
+4
+1 0 0 0
+2 1 0 0
+3 1 1 0
+4 0 1 0
+$ENDNOD
+$ELM
+5
+1 1 1 1 2 1 2
+2 1 1 2 2 2 3
+3 1 1 3 2 3 4
+4 1 1 4 2 4 1
+5 3 100 6 4 1 2 3 4
+$ENDELM
diff --git a/tests/bits/grid_in_msh_01.3d.msh b/tests/bits/grid_in_msh_01.3d.msh
new file mode 100644 (file)
index 0000000..489f312
--- /dev/null
@@ -0,0 +1,21 @@
+$NOD
+8
+1 0 0 0
+2 1 0 0
+3 1 1 0
+4 0 1 0
+5 1 0 1
+6 1 1 1
+10 0 1 1
+14 0 0 1
+$ENDNOD
+$ELM
+7
+1 3 1 6 4 1 2 3 4
+2 3 1 23 4 4 1 14 10
+3 3 1 19 4 3 4 10 6
+4 3 1 27 4 1 2 5 14
+5 3 1 15 4 2 3 6 5
+6 3 1 28 4 14 5 6 10
+7 5 100 1000 8 1 2 3 4 14 5 6 10
+$ENDELM
diff --git a/tests/bits/grid_in_msh_01.cc b/tests/bits/grid_in_msh_01.cc
new file mode 100644 (file)
index 0000000..775641f
--- /dev/null
@@ -0,0 +1,87 @@
+//----------------------------  grid_in_msh_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004 by the deal.II authors and Luca Heltai
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  grid_in_msh_01.cc  ---------------------------
+
+
+// check whether we can read in with the gmsh format
+
+#include "../tests.h"
+#include <base/logstream.h>
+
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+  
+#include <fstream>
+#include <cmath>
+  
+
+template <int dim>
+void gmsh_grid (const char *name)
+{
+  Triangulation<dim> tria;
+  GridIn<dim> grid_in;
+  grid_in.attach_triangulation (tria);
+  std::ifstream input_file(name);
+  grid_in.read_msh(input_file);
+  
+  deallog << "  " << tria.n_active_cells() << " active cells" << std::endl;
+
+  int hash = 0;
+  int index = 0;
+  for (typename Triangulation<dim>::active_cell_iterator c=tria.begin_active();
+       c!=tria.end(); ++c, ++index)
+    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+      hash += (index * i * c->vertex_index(i)) % (tria.n_active_cells()+1);
+  deallog << "  hash=" << hash << std::endl;
+}
+
+
+int main () 
+{
+  std::ofstream logfile("grid_in_msh_01.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+                               
+  try
+    {
+      gmsh_grid<2> ("grid_in_msh_01.2d.msh");
+      gmsh_grid<3> ("grid_in_msh_01.3d.msh");
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+  
+  return 0;
+}

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.