From: Luca Heltai Date: Wed, 7 Nov 2018 21:01:17 +0000 (+0100) Subject: New GridIn::read_vtk version. X-Git-Tag: v9.1.0-rc1~556^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7fa3a7550a71a112f1064f273c25591ea709341a;p=dealii.git New GridIn::read_vtk version. This version works also for dim == 1, and is more robust w.r.t. the VTK file format. --- diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h index 58c3dab83f..d646ead8f4 100644 --- a/include/deal.II/grid/grid_in.h +++ b/include/deal.II/grid/grid_in.h @@ -361,9 +361,42 @@ public: read(const std::string &in, Format format = Default); /** - * Read grid data from an vtk file. Numerical data is ignored. + * 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. + * + * 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_LINE cell types, to specify optional boundary or interior edges + * + * In two dimensions: + * + * - VTK_QUAD cell types + * - VTK_LINE cell types, to specify optional boundary or interior edges + * + * In one dimension + * + * - VTK_LINE cell types + * + * The input file may specify boundary ids, material ids, and manifold ids + * using the CELL_DATA section of the + * [VTK file format](http://www.vtk.org/VTK/img/file-formats.pdf). + * + * This function interprets two types of CELL_DATA contained in the input + * file: `SCALARS MaterialID`, used to specify the material_id of the cells, + * or the boundary_id of the faces and edges, and `SCALARS ManifoldID`, that + * can be used to specify the manifold id of any Triangulation object (cell, + * face, or edge). + * + * The companion GridOut::write_vtk function can be used to write VTK files + * compatible with this method. * * @author Mayank Sabharwal, Andreas Putz, 2013 + * @author Luca Heltai, 2018 */ void read_vtk(std::istream &in); diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 37ba68e4ac..e003ac7c08 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -113,7 +113,6 @@ template void GridIn::read_vtk(std::istream &in) { - Assert((dim == 2) || (dim == 3), ExcNotImplemented()); std::string line; // verify that the first, third and fourth lines match @@ -157,8 +156,7 @@ GridIn::read_vtk(std::istream &in) unsigned int n_vertices; in >> n_vertices; - in.ignore(256, - '\n'); // ignoring the number beside the total no. of points. + in >> keyword; // float, double, int, char, etc. for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) { @@ -177,30 +175,15 @@ GridIn::read_vtk(std::istream &in) ExcMessage( "While reading VTK file, failed to find POINTS section")); - - //////////////////ignoring space between points and cells - /// sections//////////////////// - std::string checkline; - int no; - in.ignore(256, - '\n'); // this move pointer to the next line ignoring unwanted no. - no = in.tellg(); - getline(in, checkline); - if (checkline.compare("") != 0) - { - in.seekg(no); - } - in >> keyword; - ///////////////////Processing the CELLS section that contains cells(cells) and - /// bound_quads(subcelldata)/////////////////////// + unsigned int n_geometric_objects = 0; + unsigned int n_ints; if (keyword == "CELLS") { - unsigned int n_geometric_objects; in >> n_geometric_objects; - in.ignore(256, '\n'); + in >> n_ints; // Ignore this, since we don't need it. if (dim == 3) { @@ -227,6 +210,11 @@ GridIn::read_vtk(std::istream &in) else if (type == 4) { + // 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(); for (unsigned int j = 0; j < type; @@ -235,6 +223,16 @@ GridIn::read_vtk(std::istream &in) subcelldata.boundary_quads.back().material_id = 0; } + else if (type == 2) + { + subcelldata.boundary_lines.emplace_back(); + + for (unsigned int j = 0; j < type; + j++) // loop to feed the data to the boundary + in >> subcelldata.boundary_lines.back().vertices[j]; + + subcelldata.boundary_lines.back().material_id = 0; + } else AssertThrow( @@ -285,7 +283,30 @@ GridIn::read_vtk(std::istream &in) AssertThrow( false, ExcMessage( - "While reading VTK file, unknown file type encountered")); + "While reading VTK file, unknown cell type encountered")); + } + } + else if (dim == 1) + { + for (unsigned int count = 0; count < n_geometric_objects; count++) + { + unsigned int type; + in >> type; + + if (type == 2) + { + cells.emplace_back(); + + for (unsigned int j = 0; j < type; j++) // loop to feed data + in >> cells.back().vertices[j]; + + cells.back().material_id = 0; + } + else + AssertThrow( + false, + ExcMessage( + "While reading VTK file, unknown cell type encountered")); } } else @@ -301,98 +322,144 @@ GridIn::read_vtk(std::istream &in) if (keyword == "CELL_TYPES") // Entering the cell_types section and ignoring data. { - in.ignore(256, '\n'); - - while (!in.fail() && !in.eof()) - { - in >> keyword; - if (keyword != "12" && keyword != "9") - { - break; - } - } + in >> n_ints; + int tmp_int; + for (unsigned int i = 0; i < n_ints; ++i) + in >> tmp_int; } + else + AssertThrow( + false, + ExcMessage(std::string( + "While reading VTK file, missing CELL_TYPES section. Found <" + + keyword + "> instead."))); + + // Ignore everything up to CELL_DATA + while (in >> keyword) + if (keyword == "CELL_DATA") + { + unsigned int n_ids; + in >> n_ids; + + AssertThrow(n_ids == n_geometric_objects, + ExcMessage("The VTK reader found a CELL_DATA statement " + "that lists a total of " + + Utilities::int_to_string(n_ids) + + " cell data objects, but this needs to " + "equal the number of cells (which is " + + Utilities::int_to_string(cells.size()) + + ") plus the number of quads (" + + Utilities::int_to_string( + subcelldata.boundary_quads.size()) + + " in 3d or the number of lines (" + + Utilities::int_to_string( + subcelldata.boundary_lines.size()) + + ") in 2d.")); + + const std::vector data_sets{"MaterialID", + "ManifoldID"}; + + for (unsigned int i = 0; i < data_sets.size(); ++i) + { + // Ignore everything until we get to a SCALARS data set + while (in >> keyword) + if (keyword == "SCALARS") + { + // Now see if we know about this type of data set, + // if not, just ignore everything till the next SCALARS + // keyword + std::string set = ""; + in >> keyword; + for (auto set_cmp : data_sets) + if (keyword == set_cmp) + set = keyword; + if (set == "") + // keep ignoring everything until the next SCALARS + // keyword + continue; + + // Now we got somewhere. Proceed from here. + // Ignore everything till the end of the line. + // SCALARS MaterialID 1 // this "1" is optional + in.ignore(256, '\n'); + + in >> keyword; + AssertThrow( + keyword == "LOOKUP_TABLE", + ExcMessage( + "While reading VTK file, missing keyword LOOKUP_TABLE")); - ////////////////////////Processing the CELL_DATA - /// section///////////////////////////// - - if (keyword == "CELL_DATA") - { - unsigned int n_ids; - in >> n_ids; - - AssertThrow( - n_ids == cells.size() + - (dim == 3 ? - subcelldata.boundary_quads.size() : - (dim == 2 ? subcelldata.boundary_lines.size() : 0)), - ExcMessage( - "The VTK reader found a CELL_DATA statement " - "that lists a total of " + - Utilities::int_to_string(n_ids) + - " cell data objects, but this needs to " - "equal the number of cells (which is " + - Utilities::int_to_string(cells.size()) + - ") plus the number of quads (" + - Utilities::int_to_string(subcelldata.boundary_quads.size()) + - " in 3d or the number of lines (" + - Utilities::int_to_string(subcelldata.boundary_lines.size()) + - ") in 2d.")); - - - std::string linenew; - std::string textnew[2]; - textnew[0] = "SCALARS MaterialID double"; - textnew[1] = "LOOKUP_TABLE default"; - - in.ignore(256, '\n'); - - for (unsigned int i = 0; i < 2; i++) - { - getline(in, linenew); - if (i == 0) - if (linenew.size() > textnew[0].size()) - linenew.resize(textnew[0].size()); - - AssertThrow(linenew.compare(textnew[i]) == 0, - ExcMessage( - std::string( - "While reading VTK file, failed to find <") + - textnew[i] + "> section")); - } - - // read material ids first for all cells, then for all - // faces. the assumption that cells come before all faces - // has been verified above via an assertion, so the order - // used in the following blocks makes sense - for (unsigned int i = 0; i < cells.size(); i++) - { - double id; - in >> id; - cells[i].material_id = id; - } - - if (dim == 3) - { - for (unsigned int i = 0; i < subcelldata.boundary_quads.size(); - i++) - { - double id; - in >> id; - subcelldata.boundary_quads[i].material_id = id; - } - } - else if (dim == 2) - { - for (unsigned int i = 0; i < subcelldata.boundary_lines.size(); - i++) - { - double id; - in >> id; - subcelldata.boundary_lines[i].material_id = id; - } - } - } + in >> keyword; + AssertThrow( + keyword == "default", + ExcMessage( + "While reading VTK file, missing keyword default")); + + // read material or manifold ids first for all cells, + // then for all faces, and finally for all lines. the + // assumption that cells come before all faces and + // lines has been verified above via an assertion, so + // the order used in the following blocks makes sense + for (unsigned int i = 0; i < cells.size(); i++) + { + double id; + in >> id; + if (set == "MaterialID") + cells[i].material_id = id; + else if (set == "ManifoldID") + cells[i].manifold_id = id; + else + Assert(false, ExcInternalError()); + } + + if (dim == 3) + { + for (unsigned int i = 0; + i < subcelldata.boundary_quads.size(); + i++) + { + double id; + in >> id; + if (set == "MaterialID") + subcelldata.boundary_quads[i].material_id = id; + else if (set == "ManifoldID") + subcelldata.boundary_quads[i].manifold_id = id; + else + Assert(false, ExcInternalError()); + } + for (unsigned int i = 0; + i < subcelldata.boundary_lines.size(); + i++) + { + double id; + in >> id; + if (set == "MaterialID") + subcelldata.boundary_lines[i].material_id = id; + else if (set == "ManifoldID") + subcelldata.boundary_lines[i].manifold_id = id; + else + Assert(false, ExcInternalError()); + } + } + else if (dim == 2) + { + for (unsigned int i = 0; + i < subcelldata.boundary_lines.size(); + i++) + { + double id; + in >> id; + if (set == "MaterialID") + subcelldata.boundary_lines[i].material_id = id; + else if (set == "ManifoldID") + subcelldata.boundary_lines[i].manifold_id = id; + else + Assert(false, ExcInternalError()); + } + } + } + } + } Assert(subcelldata.check_consistency(dim), ExcInternalError()); @@ -414,7 +481,6 @@ GridIn::read_vtk(std::istream &in) } - template void GridIn::read_unv(std::istream &in) diff --git a/tests/grid/grid_in_vtk_2d/mesh b/tests/grid/grid_in_vtk_2d/mesh index 6e1d90b7ed..71ed09807b 100644 --- a/tests/grid/grid_in_vtk_2d/mesh +++ b/tests/grid/grid_in_vtk_2d/mesh @@ -7,6 +7,8 @@ POINTS 4 float 0.5 0.5 0 CELLS 1 5 4 0 1 3 2 +CELL_TYPES 1 +9 CELL_DATA 1 SCALARS MaterialID double LOOKUP_TABLE default