template <int dim, int spacedim>
void
-GridIn<dim, spacedim>::read_msh(std::istream &in)
+GridIn<dim, spacedim>::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<Point<spacedim>> vertices;
-
- // set up mapping between numbering in msh-file (nod) and in the vertices
- // vector
- std::map<int, int> 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<CellData<dim>> cells;
- SubCellData subcelldata;
- std::map<unsigned int, types::boundary_id> 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<unsigned int, unsigned int> 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<std::map<int, int>, 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)
// 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")
{
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<Point<spacedim>> vertices(n_vertices);
+ // set up mapping between numbering
+ // in msh-file (nod) and in the
+ // vertices vector
+ std::map<int, int> 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<int> 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<int> 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<CellData<dim>> cells; SubCellData
- // subcelldata; std::map<unsigned int, types::boundary_id>
- // 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<unsigned int, unsigned int> 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<unsigned int, 8>
- 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<CellData<dim>> cells;
+ SubCellData subcelldata;
+ std::map<unsigned int, types::boundary_id> 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<unsigned int, unsigned int> 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<unsigned int, 8> 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<dim> &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<dim>::vertices_per_cell)
- in >> cell.vertices
- [dim == 3 ?
- local_vertex_numbering[i] :
- GeometryInfo<dim>::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<types::material_id>::max(),
- ExcIndexRange(
- material_id,
- 0,
- std::numeric_limits<types::material_id>::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<types::boundary_id>::max(),
- ExcIndexRange(
- material_id,
- 0,
- std::numeric_limits<types::boundary_id>::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<types::boundary_id>(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<types::boundary_id>::max(),
- ExcIndexRange(
- material_id,
- 0,
- std::numeric_limits<types::boundary_id>::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<types::boundary_id>(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<dim> &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<dim>::vertices_per_cell)
+ in >> cell.vertices[dim == 3 ?
+ local_vertex_numbering[i] :
+ GeometryInfo<dim>::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<types::material_id>::max(),
+ ExcIndexRange(
+ material_id,
+ 0,
+ std::numeric_limits<types::material_id>::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<types::boundary_id>::max(),
+ ExcIndexRange(
+ material_id,
+ 0,
+ std::numeric_limits<types::boundary_id>::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<types::boundary_id>(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<types::boundary_id>::max(),
+ ExcIndexRange(
+ material_id,
+ 0,
+ std::numeric_limits<types::boundary_id>::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<types::boundary_id>(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,
projects / dealii.git / commitdiff