From def16b3fa6b5a2f9a190e5fca7413523de0f2be2 Mon Sep 17 00:00:00 2001 From: Mikael Vaillant Date: Mon, 23 Sep 2024 13:06:32 -0400 Subject: [PATCH] Modified code logic to delete all the comments before parsing --- source/grid/grid_in.cc | 991 ++++++++++++++++++++--------------------- 1 file changed, 485 insertions(+), 506 deletions(-) diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 17f4d90f6f..e8bf9b1e0b 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -2002,70 +2002,55 @@ GridIn::read_comsol_mphtxt(std::istream &in) template void -GridIn::read_msh(std::istream &in) +GridIn::read_msh(std::istream &input_stream) { Assert(tria != nullptr, ExcNoTriangulationSelected()); - AssertThrow(in.fail() == false, ExcIO()); + AssertThrow(input_stream.fail() == false, ExcIO()); unsigned int n_vertices; unsigned int n_cells; - int n_entity_blocks = 1; unsigned int dummy; std::string line; - - // set up vector of vertices. This vector is filled using information from the - // $Elements or $ELM section - std::vector> vertices; - - // set up mapping between numbering in msh-file (nod) and in the vertices - // vector - std::map vertex_indices; - - // set up array of cells and subcells (faces). In 1d, there is currently no - // standard way in deal.II to pass boundary indicators attached to - // individual vertices, so do this by hand via the boundary_ids_1d array - std::vector> cells; - SubCellData subcelldata; - std::map boundary_ids_1d; - - // Track the number of times each vertex is used in 1D. This determines - // whether or not we can assign a boundary id to a vertex. This is necessary - // because sometimes gmsh saves internal vertices in the $ELEM list in codim - // 1 or codim 2. - std::map vertex_counts; - // This array stores maps from the 'entities' to the 'physical tags' for // points, curves, surfaces and volumes. We use this information later to // assign boundary ids. std::array, 4> tag_maps; - // Section markers according to mesh format - static const std::string begin_nodes_marker[] = {"$NOD", "$Nodes"}; - static const std::string begin_elements_marker[] = {"$ELM", "$Elements"}; - - in >> line; + // contain the the file stripped of the comments + std::string stripped_file; - // first determine file format and skip all preliminary information such as - // comments - unsigned int gmsh_file_format = 0; - const std::string first_line = line; - while (true) + // Comments can be included by mesh generating softwares and must be deleted, + // a string is filed with the content of the file stripped of the comments + while (std::getline(input_stream, line)) { - if (line == "$NOD") - { - gmsh_file_format = 10; - break; - } - else if (line == "$MeshFormat") + if (line == "$Comments") { - gmsh_file_format = 20; - break; + while (std::getline(input_stream, line)) + { + if (line == "$EndComments") + { + break; + } + } + continue; } - else if (in.eof()) - AssertThrow(false, ExcInvalidGMSHInput(first_line)); - in >> line; + stripped_file += line + '\n'; } + // Restart file reading + std::istringstream in(stripped_file); + + in >> line; + + // first determine file format + unsigned int gmsh_file_format = 0; + if (line == "$NOD") + gmsh_file_format = 10; + else if (line == "$MeshFormat") + gmsh_file_format = 20; + else + AssertThrow(false, ExcInvalidGMSHInput(line)); + // if file format is 2.0 or greater then we also have to read the rest of // the header if (gmsh_file_format == 20) @@ -2085,13 +2070,19 @@ GridIn::read_msh(std::istream &in) // description to synch ourselves with the format 1 handling above in >> line; AssertThrow(line == "$EndMeshFormat", ExcInvalidGMSHInput(line)); - } - // Loop through every line of the file. Processes information in $Entities, - // $Nodes/$NOD and $Elements/$ELLM section of the .msh file to build the - // triangulation - while (!in.eof()) - { + in >> line; + // if the next block is of kind $PhysicalNames, ignore it + if (line == "$PhysicalNames") + { + do + { + in >> line; + } + while (line != "$EndPhysicalNames"); + in >> line; + } + // if the next block is of kind $Entities, parse it if (line == "$Entities") { @@ -2217,495 +2208,483 @@ GridIn::read_msh(std::istream &in) in >> line; } - // if the next block is of kind $Nodes or $NOD, parse it - if (line == begin_nodes_marker[gmsh_file_format == 10 ? 0 : 1]) + // if the next block is of kind $PartitionedEntities, ignore it + if (line == "$PartitionedEntities") { - // now read the nodes list - int n_entity_blocks = 1; - if (gmsh_file_format > 40) + do { - int min_node_tag; - int max_node_tag; - in >> n_entity_blocks >> n_vertices >> min_node_tag >> - max_node_tag; + in >> line; } - else if (gmsh_file_format == 40) - { - in >> n_entity_blocks >> n_vertices; - } - else - in >> n_vertices; + while (line != "$EndPartitionedEntities"); + in >> line; + } - vertices.resize(n_vertices); + // but the next thing should, + // in any case, be the list of + // nodes: + AssertThrow(line == "$Nodes", ExcInvalidGMSHInput(line)); + } - { - unsigned int global_vertex = 0; - for (int entity_block = 0; entity_block < n_entity_blocks; - ++entity_block) - { - int parametric; - unsigned long numNodes; + // now read the nodes list + int n_entity_blocks = 1; + if (gmsh_file_format > 40) + { + int min_node_tag; + int max_node_tag; + in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag; + } + else if (gmsh_file_format == 40) + { + in >> n_entity_blocks >> n_vertices; + } + else + 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; - if (gmsh_file_format < 40) - { - numNodes = n_vertices; - parametric = 0; - } - else - { - // for gmsh_file_format 4.1 the order of tag and dim is - // reversed, but we are ignoring both anyway. - int tagEntity, dimEntity; - in >> tagEntity >> dimEntity >> parametric >> numNodes; - } + { + unsigned int global_vertex = 0; + for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block) + { + int parametric; + unsigned long numNodes; - std::vector vertex_numbers; - int vertex_number; - if (gmsh_file_format > 40) - for (unsigned long vertex_per_entity = 0; - vertex_per_entity < numNodes; - ++vertex_per_entity) - { - in >> vertex_number; - vertex_numbers.push_back(vertex_number); - } + if (gmsh_file_format < 40) + { + numNodes = n_vertices; + parametric = 0; + } + else + { + // for gmsh_file_format 4.1 the order of tag and dim is reversed, + // but we are ignoring both anyway. + int tagEntity, dimEntity; + in >> tagEntity >> dimEntity >> parametric >> numNodes; + } - for (unsigned long vertex_per_entity = 0; - vertex_per_entity < numNodes; - ++vertex_per_entity, ++global_vertex) - { - int vertex_number; - double x[3]; + std::vector vertex_numbers; + int vertex_number; + if (gmsh_file_format > 40) + for (unsigned long vertex_per_entity = 0; + vertex_per_entity < numNodes; + ++vertex_per_entity) + { + in >> vertex_number; + vertex_numbers.push_back(vertex_number); + } - // read vertex - if (gmsh_file_format > 40) - { - vertex_number = vertex_numbers[vertex_per_entity]; - in >> x[0] >> x[1] >> x[2]; - } - else - in >> vertex_number >> x[0] >> x[1] >> x[2]; + for (unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes; + ++vertex_per_entity, ++global_vertex) + { + int vertex_number; + double x[3]; - for (unsigned int d = 0; d < spacedim; ++d) - vertices[global_vertex][d] = x[d]; - // store mapping - vertex_indices[vertex_number] = global_vertex; + // read vertex + if (gmsh_file_format > 40) + { + vertex_number = vertex_numbers[vertex_per_entity]; + in >> x[0] >> x[1] >> x[2]; + } + else + in >> vertex_number >> x[0] >> x[1] >> x[2]; - // ignore parametric coordinates - if (parametric != 0) - { - double u = 0.; - double v = 0.; - in >> u >> v; - (void)u; - (void)v; - } - } + for (unsigned int d = 0; d < spacedim; ++d) + vertices[global_vertex][d] = x[d]; + // store mapping + vertex_indices[vertex_number] = global_vertex; + + // ignore parametric coordinates + if (parametric != 0) + { + double u = 0.; + double v = 0.; + in >> u >> v; + (void)u; + (void)v; } - AssertDimension(global_vertex, n_vertices); } + } + AssertDimension(global_vertex, n_vertices); + } - // Assert we reached the end of the block - in >> line; - static const std::string end_nodes_marker[] = {"$ENDNOD", - "$EndNodes"}; - AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1], - ExcInvalidGMSHInput(line)); - } + // Assert we reached the end of the block + in >> line; + static const std::string end_nodes_marker[] = {"$ENDNOD", "$EndNodes"}; + AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1], + ExcInvalidGMSHInput(line)); - // if the next block is of kind $Elements or $ELM, parse it - if (line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1]) - { - // now read the cell list - if (gmsh_file_format > 40) - { - int min_node_tag; - int max_node_tag; - in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag; - } - else if (gmsh_file_format == 40) - { - in >> n_entity_blocks >> n_cells; - } - else - { - n_entity_blocks = 1; - in >> n_cells; - } + // Now read in next bit + in >> line; + static const std::string begin_elements_marker[] = {"$ELM", "$Elements"}; + AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1], + ExcInvalidGMSHInput(line)); - // // set up array of cells and subcells (faces). In 1d, there is - // currently no - // // standard way in deal.II to pass boundary indicators attached to - // // individual vertices, so do this by hand via the boundary_ids_1d - // array std::vector> cells; SubCellData - // subcelldata; std::map - // boundary_ids_1d; - - // // Track the number of times each vertex is used in 1D. This - // determines - // // whether or not we can assign a boundary id to a vertex. This is - // necessary - // // because sometimes gmsh saves internal vertices in the $ELEM list - // in codim - // // 1 or codim 2. - // std::map vertex_counts; + // now read the cell list + if (gmsh_file_format > 40) + { + int min_node_tag; + int max_node_tag; + in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag; + } + else if (gmsh_file_format == 40) + { + in >> n_entity_blocks >> n_cells; + } + else + { + n_entity_blocks = 1; + in >> n_cells; + } - { - static constexpr std::array - local_vertex_numbering = {{0, 1, 5, 4, 2, 3, 7, 6}}; - unsigned int global_cell = 0; - for (int entity_block = 0; entity_block < n_entity_blocks; - ++entity_block) - { - unsigned int material_id; - unsigned long numElements; - int cell_type; + // set up array of cells and subcells (faces). In 1d, there is currently no + // standard way in deal.II to pass boundary indicators attached to + // individual vertices, so do this by hand via the boundary_ids_1d array + std::vector> cells; + SubCellData subcelldata; + std::map boundary_ids_1d; - if (gmsh_file_format < 40) - { - material_id = 0; - cell_type = 0; - numElements = n_cells; - } - else if (gmsh_file_format == 40) - { - int tagEntity, dimEntity; - in >> tagEntity >> dimEntity >> cell_type >> numElements; - material_id = tag_maps[dimEntity][tagEntity]; - } - else - { - // for gmsh_file_format 4.1 the order of tag and dim is - // reversed, - int tagEntity, dimEntity; - in >> dimEntity >> tagEntity >> cell_type >> numElements; - material_id = tag_maps[dimEntity][tagEntity]; - } + // Track the number of times each vertex is used in 1D. This determines + // whether or not we can assign a boundary id to a vertex. This is necessary + // because sometimes gmsh saves internal vertices in the $ELEM list in codim + // 1 or codim 2. + std::map vertex_counts; - for (unsigned int cell_per_entity = 0; - cell_per_entity < numElements; - ++cell_per_entity, ++global_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.fail() == false, ExcIO()); - - unsigned int nod_num; - - /* - For file format version 1, the format of each cell is as - follows: elm-number elm-type reg-phys reg-elem - number-of-nodes node-number-list - - However, for version 2, the format reads like this: - elm-number elm-type number-of-tags < tag > ... - node-number-list - - For version 4, we have: - tag(int) numVert(int) ... - - In the following, we will ignore the element number (we - simply enumerate them in the order in which we read them, - and we will take reg-phys (version 1) or the first tag - (version 2, if any tag is given at all) as material id. - For version 4, we already read the material and the cell - type in above. - */ - - unsigned int elm_number = 0; - if (gmsh_file_format < 40) - { - in >> elm_number // ELM-NUMBER - >> cell_type; // ELM-TYPE - } + { + static constexpr std::array local_vertex_numbering = { + {0, 1, 5, 4, 2, 3, 7, 6}}; + unsigned int global_cell = 0; + for (int entity_block = 0; entity_block < n_entity_blocks; ++entity_block) + { + unsigned int material_id; + unsigned long numElements; + int cell_type; - if (gmsh_file_format < 20) - { - in >> material_id // REG-PHYS - >> dummy // reg_elm - >> nod_num; - } - else if (gmsh_file_format < 40) - { - // read the tags; ignore all but the first one which we - // will interpret as the material_id (for cells) or - // boundary_id (for faces) - unsigned int n_tags; - in >> n_tags; - if (n_tags > 0) - in >> material_id; - else - material_id = 0; - - for (unsigned int i = 1; i < n_tags; ++i) - in >> dummy; - - 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 // file format version 4.0 and later - { - // ignore tag - int tag; - in >> tag; - - 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; - } + if (gmsh_file_format < 40) + { + material_id = 0; + cell_type = 0; + numElements = n_cells; + } + else if (gmsh_file_format == 40) + { + int tagEntity, dimEntity; + in >> tagEntity >> dimEntity >> cell_type >> numElements; + material_id = tag_maps[dimEntity][tagEntity]; + } + else + { + // for gmsh_file_format 4.1 the order of tag and dim is reversed, + int tagEntity, dimEntity; + in >> dimEntity >> tagEntity >> cell_type >> numElements; + material_id = tag_maps[dimEntity][tagEntity]; + } + for (unsigned int cell_per_entity = 0; cell_per_entity < numElements; + ++cell_per_entity, ++global_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.fail() == false, ExcIO()); + + unsigned int nod_num; + + /* + For file format version 1, the format of each cell is as + follows: elm-number elm-type reg-phys reg-elem number-of-nodes + node-number-list + + However, for version 2, the format reads like this: + elm-number elm-type number-of-tags < tag > ... + node-number-list + + For version 4, we have: + tag(int) numVert(int) ... + + In the following, we will ignore the element number (we simply + enumerate them in the order in which we read them, and we will + take reg-phys (version 1) or the first tag (version 2, if any + tag is given at all) as material id. For version 4, we already + read the material and the cell type in above. + */ + + unsigned int elm_number = 0; + if (gmsh_file_format < 40) + { + in >> elm_number // ELM-NUMBER + >> cell_type; // ELM-TYPE + } - /* `ELM-TYPE' - defines the geometrical type of the N-th element: - `1' - Line (2 nodes, 1 edge). + if (gmsh_file_format < 20) + { + in >> material_id // REG-PHYS + >> dummy // reg_elm + >> nod_num; + } + else if (gmsh_file_format < 40) + { + // read the tags; ignore all but the first one which we will + // interpret as the material_id (for cells) or boundary_id + // (for faces) + unsigned int n_tags; + in >> n_tags; + if (n_tags > 0) + in >> material_id; + else + material_id = 0; + + for (unsigned int i = 1; i < n_tags; ++i) + in >> dummy; + + 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 // file format version 4.0 and later + { + // ignore tag + int tag; + in >> tag; - `2' - Triangle (3 nodes, 3 edges). + 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; + } - `3' - Quadrangle (4 nodes, 4 edges). - `4' - Tetrahedron (4 nodes, 6 edges, 6 faces). + /* `ELM-TYPE' + defines the geometrical type of the N-th element: + `1' + Line (2 nodes, 1 edge). - `5' - Hexahedron (8 nodes, 12 edges, 6 faces). + `2' + Triangle (3 nodes, 3 edges). - `15' - Point (1 node). - */ + `3' + Quadrangle (4 nodes, 4 edges). - if (((cell_type == 1) && (dim == 1)) || // a line in 1d - ((cell_type == 2) && (dim == 2)) || // a triangle in 2d - ((cell_type == 3) && - (dim == 2)) || // a quadrilateral in 2d - ((cell_type == 4) && (dim == 3)) || // a tet in 3d - ((cell_type == 5) && (dim == 3))) // a hex in 3d - // found a 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; - else if (cell_type == 3) // quad - vertices_per_cell = 4; - else if (cell_type == 4) // tet - vertices_per_cell = 4; - else if (cell_type == 5) // hex - vertices_per_cell = 8; + `4' + Tetrahedron (4 nodes, 6 edges, 6 faces). - 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(); - CellData &cell = cells.back(); - cell.vertices.resize(vertices_per_cell); - for (unsigned int i = 0; i < vertices_per_cell; ++i) - { - // hypercube cells need to be reordered - if (vertices_per_cell == - GeometryInfo::vertices_per_cell) - in >> cell.vertices - [dim == 3 ? - local_vertex_numbering[i] : - GeometryInfo::ucd_to_deal[i]]; - else - in >> cell.vertices[i]; - } + `5' + Hexahedron (8 nodes, 12 edges, 6 faces). - // to make sure that the cast won't fail - Assert( - material_id <= - std::numeric_limits::max(), - ExcIndexRange( - material_id, - 0, - std::numeric_limits::max())); - // we use only material_ids in the range from 0 to - // numbers::invalid_material_id-1 - AssertIndexRange(material_id, - numbers::invalid_material_id); - - cell.material_id = material_id; - - // transform from gmsh to consecutive numbering - for (unsigned int i = 0; i < vertices_per_cell; ++i) - { - AssertThrow( - vertex_indices.find(cell.vertices[i]) != - vertex_indices.end(), - ExcInvalidVertexIndexGmsh(cell_per_entity, - elm_number, - cell.vertices[i])); - - const auto vertex = - vertex_indices[cell.vertices[i]]; - if (dim == 1) - vertex_counts[vertex] += 1u; - cell.vertices[i] = vertex; - } - } - else if ((cell_type == 1) && - ((dim == 2) || (dim == 3))) // a line in 2d or 3d - // boundary info - { - subcelldata.boundary_lines.emplace_back(); - in >> subcelldata.boundary_lines.back().vertices[0] >> - subcelldata.boundary_lines.back().vertices[1]; - - // to make sure that the cast won't fail - Assert( - material_id <= - std::numeric_limits::max(), - ExcIndexRange( - material_id, - 0, - std::numeric_limits::max())); - // we use only boundary_ids in the range from 0 to - // numbers::internal_face_boundary_id-1 - AssertIndexRange(material_id, - numbers::internal_face_boundary_id); - - subcelldata.boundary_lines.back().boundary_id = - static_cast(material_id); - - // transform from ucd to - // consecutive numbering - for (unsigned int &vertex : - subcelldata.boundary_lines.back().vertices) - if (vertex_indices.find(vertex) != - vertex_indices.end()) - // vertex with this index exists - vertex = vertex_indices[vertex]; - else - { - // no such vertex index - AssertThrow(false, - ExcInvalidVertexIndex(cell_per_entity, - vertex)); - vertex = numbers::invalid_unsigned_int; - } - } - else if ((cell_type == 2 || cell_type == 3) && - (dim == 3)) // a triangle or a quad in 3d - // boundary info - { - unsigned int vertices_per_cell = 0; - // check cell type - if (cell_type == 2) // tri - vertices_per_cell = 3; - else if (cell_type == 3) // quad - vertices_per_cell = 4; - - subcelldata.boundary_quads.emplace_back(); - - // 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 <= - std::numeric_limits::max(), - ExcIndexRange( - material_id, - 0, - std::numeric_limits::max())); - // we use only boundary_ids in the range from 0 to - // numbers::internal_face_boundary_id-1 - AssertIndexRange(material_id, - numbers::internal_face_boundary_id); - - subcelldata.boundary_quads.back().boundary_id = - static_cast(material_id); - - // transform from gmsh to - // consecutive numbering - for (unsigned int &vertex : - subcelldata.boundary_quads.back().vertices) - if (vertex_indices.find(vertex) != - vertex_indices.end()) - // vertex with this index exists - vertex = vertex_indices[vertex]; - else - { - // no such vertex index - Assert(false, - ExcInvalidVertexIndex(cell_per_entity, - vertex)); - vertex = numbers::invalid_unsigned_int; - } - } - else if (cell_type == 15) - { - // read the indices of nodes given - unsigned int node_index = 0; - if (gmsh_file_format < 20) - { - // For points (cell_type==15), we can only ever - // list one node index. - AssertThrow(nod_num == 1, ExcInternalError()); - in >> node_index; - } - else - { - in >> node_index; - } + `15' + Point (1 node). + */ - // we only care about boundary indicators assigned to - // individual vertices in 1d (because otherwise the - // vertices are not faces) - if (dim == 1) - boundary_ids_1d[vertex_indices[node_index]] = - material_id; - } + if (((cell_type == 1) && (dim == 1)) || // a line in 1d + ((cell_type == 2) && (dim == 2)) || // a triangle in 2d + ((cell_type == 3) && (dim == 2)) || // a quadrilateral in 2d + ((cell_type == 4) && (dim == 3)) || // a tet in 3d + ((cell_type == 5) && (dim == 3))) // a hex in 3d + // found a 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; + else if (cell_type == 3) // quad + vertices_per_cell = 4; + else if (cell_type == 4) // tet + vertices_per_cell = 4; + else if (cell_type == 5) // hex + vertices_per_cell = 8; + + 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(); + CellData &cell = cells.back(); + cell.vertices.resize(vertices_per_cell); + for (unsigned int i = 0; i < vertices_per_cell; ++i) + { + // hypercube cells need to be reordered + if (vertices_per_cell == + GeometryInfo::vertices_per_cell) + in >> cell.vertices[dim == 3 ? + local_vertex_numbering[i] : + GeometryInfo::ucd_to_deal[i]]; else - { - AssertThrow(false, - ExcGmshUnsupportedGeometry(cell_type)); - } + in >> cell.vertices[i]; + } + + // to make sure that the cast won't fail + Assert(material_id <= + std::numeric_limits::max(), + ExcIndexRange( + material_id, + 0, + std::numeric_limits::max())); + // we use only material_ids in the range from 0 to + // numbers::invalid_material_id-1 + AssertIndexRange(material_id, numbers::invalid_material_id); + + cell.material_id = material_id; + + // transform from gmsh to consecutive numbering + for (unsigned int i = 0; i < vertices_per_cell; ++i) + { + AssertThrow(vertex_indices.find(cell.vertices[i]) != + vertex_indices.end(), + ExcInvalidVertexIndexGmsh(cell_per_entity, + elm_number, + cell.vertices[i])); + + const auto vertex = vertex_indices[cell.vertices[i]]; + if (dim == 1) + vertex_counts[vertex] += 1u; + cell.vertices[i] = vertex; + } + } + else if ((cell_type == 1) && + ((dim == 2) || (dim == 3))) // a line in 2d or 3d + // boundary info + { + subcelldata.boundary_lines.emplace_back(); + in >> subcelldata.boundary_lines.back().vertices[0] >> + subcelldata.boundary_lines.back().vertices[1]; + + // to make sure that the cast won't fail + Assert(material_id <= + std::numeric_limits::max(), + ExcIndexRange( + material_id, + 0, + std::numeric_limits::max())); + // we use only boundary_ids in the range from 0 to + // numbers::internal_face_boundary_id-1 + AssertIndexRange(material_id, + numbers::internal_face_boundary_id); + + subcelldata.boundary_lines.back().boundary_id = + static_cast(material_id); + + // transform from ucd to + // consecutive numbering + for (unsigned int &vertex : + subcelldata.boundary_lines.back().vertices) + if (vertex_indices.find(vertex) != vertex_indices.end()) + // vertex with this index exists + vertex = vertex_indices[vertex]; + else + { + // no such vertex index + AssertThrow(false, + ExcInvalidVertexIndex(cell_per_entity, + vertex)); + vertex = numbers::invalid_unsigned_int; + } + } + else if ((cell_type == 2 || cell_type == 3) && + (dim == 3)) // a triangle or a quad in 3d + // boundary info + { + unsigned int vertices_per_cell = 0; + // check cell type + if (cell_type == 2) // tri + vertices_per_cell = 3; + else if (cell_type == 3) // quad + vertices_per_cell = 4; + + subcelldata.boundary_quads.emplace_back(); + + // 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 <= + std::numeric_limits::max(), + ExcIndexRange( + material_id, + 0, + std::numeric_limits::max())); + // we use only boundary_ids in the range from 0 to + // numbers::internal_face_boundary_id-1 + AssertIndexRange(material_id, + numbers::internal_face_boundary_id); + + subcelldata.boundary_quads.back().boundary_id = + static_cast(material_id); + + // transform from gmsh to + // consecutive numbering + for (unsigned int &vertex : + subcelldata.boundary_quads.back().vertices) + if (vertex_indices.find(vertex) != vertex_indices.end()) + // vertex with this index exists + vertex = vertex_indices[vertex]; + else + { + // no such vertex index + Assert(false, + ExcInvalidVertexIndex(cell_per_entity, vertex)); + vertex = numbers::invalid_unsigned_int; + } + } + else if (cell_type == 15) + { + // read the indices of nodes given + unsigned int node_index = 0; + if (gmsh_file_format < 20) + { + // For points (cell_type==15), we can only ever + // list one node index. + AssertThrow(nod_num == 1, ExcInternalError()); + in >> node_index; + } + else + { + in >> node_index; } + + // we only care about boundary indicators assigned to + // individual vertices in 1d (because otherwise the vertices + // are not faces) + if (dim == 1) + boundary_ids_1d[vertex_indices[node_index]] = material_id; + } + else + { + AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type)); } - AssertDimension(global_cell, n_cells); } - // Assert that we reached the end of the block - in >> line; - static const std::string end_elements_marker[] = {"$ENDELM", - "$EndElements"}; - AssertThrow(line == - end_elements_marker[gmsh_file_format == 10 ? 0 : 1], - ExcInvalidGMSHInput(line)); - } - // skip the line if not recognized - in >> line; - } + } + AssertDimension(global_cell, n_cells); + } + // Assert that we reached the end of the block + in >> line; + static const std::string end_elements_marker[] = {"$ENDELM", "$EndElements"}; + AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1], + ExcInvalidGMSHInput(line)); + AssertThrow(in.fail() == false, ExcIO()); // check that we actually read some cells. AssertThrow(cells.size() > 0, -- 2.39.5