From: Peter Munch Date: Thu, 23 Jul 2020 19:28:26 +0000 (+0200) Subject: GridIn::read_vtk() for simplex meshes X-Git-Tag: v9.3.0-rc1~1244^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10748%2Fhead;p=dealii.git GridIn::read_vtk() for simplex meshes --- diff --git a/doc/news/changes/minor/20200723Munch b/doc/news/changes/minor/20200723Munch new file mode 100644 index 0000000000..801696c752 --- /dev/null +++ b/doc/news/changes/minor/20200723Munch @@ -0,0 +1,3 @@ +New: GridIn::read_vtk() now supports triangular and tetrahedral meshes. +
+(Peter Munch, 2020/07/23) diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h index b76d8a461e..7bb542fc09 100644 --- a/include/deal.II/grid/grid_in.h +++ b/include/deal.II/grid/grid_in.h @@ -364,20 +364,22 @@ public: /** * Read grid data from a unstructured vtk file. The vtk file may contain - * the following VTK cell types: VTK_HEXAHEDRON, VTK_QUAD, and VTK_LINE. + * the following VTK cell types: VTK_HEXAHEDRON (12), VTK_TETRA (10), + * VTK_QUAD (9), VTK_TRIANGLE (5), and VTK_LINE (3). * * Depending on the template dimension, only some of the above are accepted. * * In particular, in three dimensions, this function expects the file to * contain * - * - VTK_HEXAHEDRON cell types - * - VTK_QUAD cell types, to specify optional boundary or interior quad faces + * - VTK_HEXAHEDRON/VTK_TETRA cell types + * - VTK_QUAD/VTK_TRIANGLE cell types, to specify optional boundary or + * interior quad faces * - VTK_LINE cell types, to specify optional boundary or interior edges * * In two dimensions: * - * - VTK_QUAD cell types + * - VTK_QUAD/VTK_TRIANGLE cell types * - VTK_LINE cell types, to specify optional boundary or interior edges * * In one dimension @@ -396,6 +398,8 @@ public: * * The companion GridOut::write_vtk function can be used to write VTK files * compatible with this method. + * + * @ingroup simplex */ void read_vtk(std::istream &in); diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 121f18783f..9f3db2828d 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -176,8 +176,33 @@ GridIn::read_vtk(std::istream &in) unsigned int n_geometric_objects = 0; unsigned int n_ints; + bool is_quad_or_hex_mesh = false; + bool is_tria_or_tet_mesh = false; + if (keyword == "CELLS") { + // jump to the `CELL_TYPES` section and read in cell types + std::vector cell_types; + { + std::streampos oldpos = in.tellg(); + + + while (in >> keyword) + if (keyword == "CELL_TYPES") + { + in >> n_ints; + + cell_types.resize(n_ints); + + for (unsigned int i = 0; i < n_ints; ++i) + in >> cell_types[i]; + + break; + } + + in.seekg(oldpos); + } + in >> n_geometric_objects; in >> n_ints; // Ignore this, since we don't need it. @@ -188,15 +213,20 @@ GridIn::read_vtk(std::istream &in) unsigned int type; in >> type; - if (type == 8) + if (cell_types[count] == 10 || cell_types[count] == 12) { + if (cell_types[count] == 10) + is_tria_or_tet_mesh = true; + if (cell_types[count] == 12) + is_quad_or_hex_mesh = true; + // we assume that the file contains first all cells, // and only then any faces or lines AssertThrow(subcelldata.boundary_quads.size() == 0 && subcelldata.boundary_lines.size() == 0, ExcNotImplemented()); - cells.emplace_back(); + cells.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed data in >> cells.back().vertices[j]; @@ -204,14 +234,19 @@ GridIn::read_vtk(std::istream &in) cells.back().material_id = 0; } - else if (type == 4) + else if (cell_types[count] == 5 || cell_types[count] == 9) { + if (cell_types[count] == 5) + is_tria_or_tet_mesh = true; + if (cell_types[count] == 9) + is_quad_or_hex_mesh = true; + // we assume that the file contains first all cells, // then all faces, and finally all lines AssertThrow(subcelldata.boundary_lines.size() == 0, ExcNotImplemented()); - subcelldata.boundary_quads.emplace_back(); + subcelldata.boundary_quads.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed the data to the boundary @@ -219,9 +254,9 @@ GridIn::read_vtk(std::istream &in) subcelldata.boundary_quads.back().material_id = 0; } - else if (type == 2) + else if (cell_types[count] == 3) { - subcelldata.boundary_lines.emplace_back(); + subcelldata.boundary_lines.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed the data to the boundary @@ -237,7 +272,6 @@ GridIn::read_vtk(std::istream &in) "While reading VTK file, unknown file type encountered")); } } - else if (dim == 2) { for (unsigned int count = 0; count < n_geometric_objects; count++) @@ -245,14 +279,19 @@ GridIn::read_vtk(std::istream &in) unsigned int type; in >> type; - if (type == 4) + if (cell_types[count] == 5 || cell_types[count] == 9) { // we assume that the file contains first all cells, // and only then any faces AssertThrow(subcelldata.boundary_lines.size() == 0, ExcNotImplemented()); - cells.emplace_back(); + if (cell_types[count] == 5) + is_tria_or_tet_mesh = true; + if (cell_types[count] == 9) + is_quad_or_hex_mesh = true; + + cells.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed data in >> cells.back().vertices[j]; @@ -260,11 +299,11 @@ GridIn::read_vtk(std::istream &in) cells.back().material_id = 0; } - else if (type == 2) + else if (cell_types[count] == 3) { // If this is encountered, the pointer comes out of the loop // and starts processing boundaries. - subcelldata.boundary_lines.emplace_back(); + subcelldata.boundary_lines.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed the data to the boundary @@ -290,10 +329,10 @@ GridIn::read_vtk(std::istream &in) in >> type; AssertThrow( - type == 2, + cell_types[count] == 3 && type == 2, ExcMessage( "While reading VTK file, unknown cell type encountered")); - cells.emplace_back(); + cells.emplace_back(type); for (unsigned int j = 0; j < type; j++) // loop to feed data in >> cells.back().vertices[j]; @@ -476,16 +515,36 @@ GridIn::read_vtk(std::istream &in) Assert(subcelldata.check_consistency(dim), ExcInternalError()); - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - if (dim == spacedim) - GridReordering::invert_all_cells_of_negative_grid( - vertices, cells); + // 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_tria_or_tet_mesh ^ is_quad_or_hex_mesh), + ExcNotImplemented()); + + if (dim == 1 || is_quad_or_hex_mesh) + { + GridTools::delete_unused_vertices(vertices, cells, subcelldata); + + if (dim == spacedim) + GridReordering::invert_all_cells_of_negative_grid( + vertices, cells); - GridReordering::reorder_cells(cells); - tria->create_triangulation_compatibility(vertices, cells, subcelldata); + GridReordering::reorder_cells(cells); + tria->create_triangulation_compatibility(vertices, + cells, + subcelldata); - return; + return; + } + else + { + tria->create_triangulation(vertices, cells, subcelldata); + } } else AssertThrow(false, diff --git a/source/grid/tria.cc b/source/grid/tria.cc index c4e207e428..3805f69cb1 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_vtk.cc b/tests/simplex/grid_in_vtk.cc new file mode 100644 index 0000000000..4b32c6d2fd --- /dev/null +++ b/tests/simplex/grid_in_vtk.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) +{ + Triangulation tria; + + GridIn grid_in; + grid_in.attach_triangulation(tria); + std::ifstream input_file(file_name); + grid_in.read_vtk(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_vtk/tri.vtk")); + deallog.pop(); + + // dim = 2, spacedim = 3 + deallog.push("triangluar_elements_dim2_spacedim3: "); + check_file<2, 3>(std::string(SOURCE_DIR "/grid_in_vtk/tri.vtk")); + deallog.pop(); + + + + // TETRAHEDRAL ELEMENTS + // dim = spacedim = 3 + deallog.push("tetrahedral_elements_dim3_spacedim3: "); + check_file<3, 3>(std::string(SOURCE_DIR "/grid_in_vtk/tet.vtk")); + deallog.pop(); +} diff --git a/tests/simplex/grid_in_vtk.with_simplex_support=on.output b/tests/simplex/grid_in_vtk.with_simplex_support=on.output new file mode 100644 index 0000000000..446eb9c86a --- /dev/null +++ b/tests/simplex/grid_in_vtk.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 +0 0 0 0 + + + +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 +0 0 0 0 + + + +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 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + + + +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_vtk/tet.vtk b/tests/simplex/grid_in_vtk/tet.vtk new file mode 100644 index 0000000000..661934ca7f --- /dev/null +++ b/tests/simplex/grid_in_vtk/tet.vtk @@ -0,0 +1,99 @@ +# vtk DataFile Version 3.0 +tet, Created by Gmsh +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 14 double +0 0 1 +0 0 0 +0 1 1 +0 1 0 +1 0 1 +1 0 0 +1 1 1 +1 1 0 +0 0.5 0.5 +1 0.5 0.5 +0.5 0 0.5 +0.5 1 0.5 +0.5 0.5 0 +0.5 0.5 1 + +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 CellEntityIds 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 diff --git a/tests/simplex/grid_in_vtk/tri.vtk b/tests/simplex/grid_in_vtk/tri.vtk new file mode 100644 index 0000000000..672448d625 --- /dev/null +++ b/tests/simplex/grid_in_vtk/tri.vtk @@ -0,0 +1,30 @@ +# vtk DataFile Version 3.0 +tri, Created by Gmsh +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 5 double +0 0 0 +1 0 0 +1 1 0 +0 1 0 +0.5 0.5 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 CellEntityIds int 1 +LOOKUP_TABLE default +1 +1 +1 +1