From f885e1ab5bcb9fb012cabed89cedd0b69280c01b Mon Sep 17 00:00:00 2001 From: Daniel Paukner Date: Tue, 14 Jul 2020 14:41:46 +0200 Subject: [PATCH] GridIn::read_msh() now supports triangular and tetrahedral elements --- doc/news/changes/minor/20200720DanielPaukner | 3 + include/deal.II/grid/grid_in.h | 2 + source/grid/grid_in.cc | 141 +++++++++++++----- source/grid/tria.cc | 2 +- tests/simplex/grid_in_msh.cc | 73 +++++++++ ...grid_in_msh.with_simplex_support=on.output | 135 +++++++++++++++++ tests/simplex/grid_in_msh/tet.msh | 107 +++++++++++++ tests/simplex/grid_in_msh/tri.msh | 41 +++++ 8 files changed, 467 insertions(+), 37 deletions(-) create mode 100644 doc/news/changes/minor/20200720DanielPaukner create mode 100644 tests/simplex/grid_in_msh.cc create mode 100644 tests/simplex/grid_in_msh.with_simplex_support=on.output create mode 100644 tests/simplex/grid_in_msh/tet.msh create mode 100644 tests/simplex/grid_in_msh/tri.msh diff --git a/doc/news/changes/minor/20200720DanielPaukner b/doc/news/changes/minor/20200720DanielPaukner new file mode 100644 index 0000000000..8a6e86ce3e --- /dev/null +++ b/doc/news/changes/minor/20200720DanielPaukner @@ -0,0 +1,3 @@ +New: GridIn::read_msh() now supports triangular and tetrahedral meshes. +
+(Daniel Paukner, 2020/07/20) diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h index b76d8a461e..74a6188e13 100644 --- 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 121f18783f..0fe7e29098 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -1731,6 +1731,8 @@ GridIn::read_msh(std::istream &in) std::vector> cells; SubCellData subcelldata; std::map boundary_ids_1d; + bool is_hex_mesh = false; + bool is_tet_mesh = false; { unsigned int global_cell = 0; @@ -1818,14 +1820,33 @@ GridIn::read_msh(std::istream &in) for (unsigned int i = 1; i < n_tags; ++i) in >> dummy; - nod_num = GeometryInfo::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::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::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::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::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::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::read_msh(std::istream &in) // transform from ucd to // consecutive numbering - for (const unsigned int i : GeometryInfo::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::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::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::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::invert_all_cells_of_negative_grid(vertices, - cells); - GridReordering::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::invert_all_cells_of_negative_grid( + vertices, cells); + GridReordering::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 eb0b44e86c..1134124599 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -9743,7 +9743,7 @@ Triangulation::create_triangulation( cell != this_round.end(); ++cell) { - for (const unsigned int i : GeometryInfo::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 index 0000000000..397ddab50f --- /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 +#include +#include + +#include + +#include "../tests.h" + +template +void +check_file(const std::string &file_name, const unsigned int vertices_per_cell) +{ + Triangulation tria; + + GridIn 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 index 0000000000..b1fa345b50 --- /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 index 0000000000..16e73d3584 --- /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 index 0000000000..822926ef4c --- /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 -- 2.39.5