From: wolf Date: Wed, 21 Apr 2004 21:28:06 +0000 (+0000) Subject: Add Luca Heltai's reader for gmsh format. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fead4ade354b6bc28e523c751c5b904de415a3c7;p=dealii-svn.git Add Luca Heltai's reader for gmsh format. git-svn-id: https://svn.dealii.org/trunk@9080 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_in.h b/deal.II/deal.II/include/grid/grid_in.h index 35199ae262..b13588b715 100644 --- a/deal.II/deal.II/include/grid/grid_in.h +++ b/deal.II/deal.II/include/grid/grid_in.h @@ -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. + * + *
  • Gmsh mesh 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). * * * @@ -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 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 diff --git a/deal.II/deal.II/source/grid/grid_in.cc b/deal.II/deal.II/source/grid/grid_in.cc index f65b079cf5..4af96c5da7 100644 --- a/deal.II/deal.II/source/grid/grid_in.cc +++ b/deal.II/deal.II/source/grid/grid_in.cc @@ -527,6 +527,207 @@ void GridIn<3>::read_xda (std::istream &in) +template +void GridIn::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 > vertices (n_vertices); + // set up mapping between numbering + // in msh-file (nod) and in the + // vertices vector + std::map vertex_indices; + + for (unsigned int vertex=0; vertex> vertex_number + >> x[0] >> x[1] >> x[2]; + + for (unsigned int d=0; d> n_cells; + // set up array of cells + std::vector > cells; + SubCellData subcelldata; + + for (unsigned int cell=0; cell> 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()); + for (unsigned int i=0; i::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::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::reorder_cells (cells); + tria->create_triangulation (vertices, cells, subcelldata); +} + + + template void GridIn::skip_empty_lines (std::istream &in) diff --git a/tests/bits/Makefile b/tests/bits/Makefile index 1726bb96c3..bbc4bfd050 100644 --- a/tests/bits/Makefile +++ b/tests/bits/Makefile @@ -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 index 0000000000..ef628fb294 --- /dev/null +++ b/tests/bits/grid_in_msh_01.2d.msh @@ -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 index 0000000000..489f3127c8 --- /dev/null +++ b/tests/bits/grid_in_msh_01.3d.msh @@ -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 index 0000000000..775641fef1 --- /dev/null +++ b/tests/bits/grid_in_msh_01.cc @@ -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 + +#include +#include +#include +#include +#include + +#include +#include + + +template +void gmsh_grid (const char *name) +{ + Triangulation tria; + GridIn 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::active_cell_iterator c=tria.begin_active(); + c!=tria.end(); ++c, ++index) + for (unsigned int i=0; i::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; +}