From: Wolfgang Bangerth Date: Tue, 16 Nov 2021 00:42:06 +0000 (-0700) Subject: Add the ability to read COMSOL .mphtxt files. X-Git-Tag: v9.4.0-rc1~820^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7cdde2c706bda0b656a4965f8ee63355292e6923;p=dealii.git Add the ability to read COMSOL .mphtxt files. --- diff --git a/doc/doxygen/images/comsol-mesh-boundary-lines.png b/doc/doxygen/images/comsol-mesh-boundary-lines.png new file mode 100644 index 0000000000..c1665cf59c Binary files /dev/null and b/doc/doxygen/images/comsol-mesh-boundary-lines.png differ diff --git a/doc/doxygen/images/comsol-mesh-boundary-triangles.png b/doc/doxygen/images/comsol-mesh-boundary-triangles.png new file mode 100644 index 0000000000..69901509f0 Binary files /dev/null and b/doc/doxygen/images/comsol-mesh-boundary-triangles.png differ diff --git a/doc/doxygen/images/comsol-mesh-boundary-volume-mesh.png b/doc/doxygen/images/comsol-mesh-boundary-volume-mesh.png new file mode 100644 index 0000000000..4c8fc908fe Binary files /dev/null and b/doc/doxygen/images/comsol-mesh-boundary-volume-mesh.png differ diff --git a/doc/doxygen/images/comsol-mesh-marked-lines.png b/doc/doxygen/images/comsol-mesh-marked-lines.png new file mode 100644 index 0000000000..784d616467 Binary files /dev/null and b/doc/doxygen/images/comsol-mesh-marked-lines.png differ diff --git a/doc/doxygen/images/comsol-mesh-marked-triangles.png b/doc/doxygen/images/comsol-mesh-marked-triangles.png new file mode 100644 index 0000000000..cd00cc367f Binary files /dev/null and b/doc/doxygen/images/comsol-mesh-marked-triangles.png differ diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h index a2d0b42366..3015ecd240 100644 --- a/include/deal.II/grid/grid_in.h +++ b/include/deal.II/grid/grid_in.h @@ -570,7 +570,7 @@ public: * value, i.e., zero) * - some elements have manifold id 1 and material id equal to 2 * - * If the physical groups are not named, then the behaviour is the same as + * If the physical groups are not named, then the behavior is the same as * the other read_msh() function, i.e., the physical tag itself is interpreted * as a boundary or material id. * @@ -580,6 +580,44 @@ public: read_msh(const std::string &filename); #endif + /** + * Read grid data from a `.mphtxt` file. `.mphtxt` is one of the file formats + * typically generated by COMSOL. The file format is described at + * http://victorsndvg.github.io/FEconv/formats/mphtxt.xhtml . + * + * The reader interprets the "geometric entity indicators" that COMSOL + * writes into these files as either boundary indicators (for edges and faces + * of cells) or as material ids (for cells). See the glossary for a + * description of + * @ref GlossBoundaryIndicator "boundary indicators". + * and + * @ref GlossMaterialId "material indicators" + * + * COMSOL has a habit of assigning "geometric entity indicators" not only + * to edges and faces on the actual boundary, but also to interior faces + * and edges. For example, for the following volume mesh generated by + * COMSOL, + * @image html "comsol-mesh-boundary-volume-mesh.png" + * the marked edges and faces are as follows: + * @image html "comsol-mesh-marked-lines.png" + * @image html "comsol-mesh-marked-triangles.png" + * Here, some of the marked lines and faces with explicitly given + * geometric entity indicators are in the *interior* of the domain -- an + * artifact of the geometry description that was used to describe + * the mesh. However, we can of course not assign boundary indicators to + * interior edges and faces. As a consequence, this reader function simply + * ignores the geometric entity indicator for edges and faces that + * are not in fact on the boundary of the domain. The result is then a mesh + * in which only the following edges and faces are explicitly assigned + * boundary indicators: + * @image html "comsol-mesh-boundary-lines.png" + * @image html "comsol-mesh-boundary-triangles.png" + * + * @ingroup simplex + */ + void + read_comsol_mphtxt(std::istream &in); + /** * Read grid data from a file containing tecplot ASCII data. This also works * in the absence of any tecplot installation. diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index f6ec8f94ce..d64e89cd1f 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -1446,6 +1446,485 @@ GridIn::read_xda(std::istream &in) +template +void +GridIn::read_comsol_mphtxt(std::istream &in) +{ + Assert(tria != nullptr, ExcNoTriangulationSelected()); + AssertThrow(in, ExcIO()); + + // Start by making our life a bit easier: The file format + // allows for comments in a whole bunch of places, including + // on separate lines, at line ends, and that's just a hassle to + // parse because we will have to check in every line whether there + // is a comment. To make things easier, just read it all in up + // front, strip comments, eat trailing whitespace, and + // concatenate it all into one big string from which we will + // then read. We lose the ability to output error messages tied + // to individual lines of the input, but none of the other + // readers does that either. + std::stringstream whole_file; + while (in) + { + // read one line + std::string line; + std::getline(in, line); + + // Strip trailing comments, then strip whatever spaces are at the end + // of the line, and if anything is left, concatenate that to the previous + // content of the file: + if (line.find('#') != std::string::npos) + line.erase(line.find('#'), std::string::npos); + while ((line.size() > 0) && (line.back() == ' ')) + line.erase(line.size() - 1); + + if (line.size() > 0) + whole_file << '\n' << line; + } + + // Now start to read the contents of this so-simplified file. A typical + // header of these files will look like this: + // # Created by COMSOL Multiphysics. + // + // # Major & minor version + // 0 1 + // 1 # number of tags + // # Tags + // 5 mesh1 + // 1 # number of types + // # Types + // 3 obj + + AssertThrow(whole_file, ExcIO()); + + { + unsigned int version_major, version_minor; + whole_file >> version_major >> version_minor; + AssertThrow((version_major == 0) && (version_minor == 1), + ExcMessage("deal.II can currently only read version 0.1 " + "of the mphtxt file format.")); + } + + // It's not clear what the 'tags' are, but read them and discard them + { + unsigned int n_tags; + whole_file >> n_tags; + for (unsigned int i = 0; i < n_tags; ++i) + { + std::string dummy; + while (whole_file.peek() == '\n') + whole_file.get(); + std::getline(whole_file, dummy); + } + } + + // Do the same with the 'types' + { + unsigned int n_types; + whole_file >> n_types; + for (unsigned int i = 0; i < n_types; ++i) + { + std::string dummy; + while (whole_file.peek() == '\n') + whole_file.get(); + std::getline(whole_file, dummy); + } + } + + // Then move on to the actual mesh. A typical header of this part will + // look like this: + // # --------- Object 0 ---------- + // + // 0 0 1 + // 4 Mesh # class + // 4 # version + // 3 # sdim + // 1204 # number of mesh vertices + // 0 # lowest mesh vertex index + // + // # Mesh vertex coordinates + // ... + AssertThrow(whole_file, ExcIO()); + { + unsigned int dummy; + whole_file >> dummy >> dummy >> dummy; + } + { + std::string s; + while (whole_file.peek() == '\n') + whole_file.get(); + std::getline(whole_file, s); + AssertThrow(s == "4 Mesh", ExcNotImplemented()); + } + { + unsigned int version; + whole_file >> version; + AssertThrow(version == 4, ExcNotImplemented()); + } + { + unsigned int file_space_dim; + whole_file >> file_space_dim; + + AssertThrow(file_space_dim == spacedim, + ExcMessage( + "The mesh file uses a different number of space dimensions " + "than the triangulation you want to read it into.")); + } + unsigned int n_vertices; + whole_file >> n_vertices; + + unsigned int starting_vertex_index; + whole_file >> starting_vertex_index; + + std::vector> vertices(n_vertices); + for (unsigned int v = 0; v < n_vertices; ++v) + whole_file >> vertices[v]; + + // Then comes a block that looks like this: + // 4 # number of element types + // + // # Type #0 + // 3 vtx # type name + // + // + // 1 # number of vertices per element + // 18 # number of elements + // # Elements + // 4 + // 12 + // 19 + // 80 + // 143 + // [...] + // 1203 + // + // 18 # number of geometric entity indices + // # Geometric entity indices + // 2 + // 0 + // 11 + // 6 + // 3 + // [...] + AssertThrow(whole_file, ExcIO()); + + std::vector> cells; + SubCellData subcelldata; + + unsigned int n_types; + whole_file >> n_types; + for (unsigned int type = 0; type < n_types; ++type) + { + std::string object_type; + { + unsigned int dummy; + whole_file >> dummy; + } + whole_file >> object_type; + + unsigned int n_vertices_per_element; + whole_file >> n_vertices_per_element; + + unsigned int n_elements; + whole_file >> n_elements; + + if (object_type == "vtx") + { + AssertThrow(n_vertices_per_element == 1, ExcInternalError()); + } + else if (object_type == "edg") + { + AssertThrow(n_vertices_per_element == 2, ExcInternalError()); + if (dim == 1) + cells.resize(n_elements); + else + subcelldata.boundary_lines.resize(n_elements); + } + else if (object_type == "tri") + { + AssertThrow(dim >= 2, + ExcMessage("Triangles should not appear in input files " + "for 1d meshes.")); + AssertThrow(n_vertices_per_element == 3, ExcInternalError()); + if (dim == 2) + cells.resize(n_elements); + else + subcelldata.boundary_quads.resize(n_elements); + } + else if (object_type == "tet") + { + AssertThrow(dim >= 3, + ExcMessage("Tetrahedra should not appear in input files " + "for 1d or 2d meshes.")); + AssertThrow(n_vertices_per_element == 4, ExcInternalError()); + if (dim == 3) + cells.resize(n_elements); + else + Assert(false, ExcInternalError()); + } + else + AssertThrow(false, ExcNotImplemented()); + + + // Next, for each element read the vertex numbers. Then we have to decide + // what to do with it. If it is a vertex, we ignore the information. + // If it is a cell, we have to put it into the appropriate object, and the + // same if it is an edge or face. + // + // In any case, we adjust vertex indices right after reading them based on + // the starting index read above + std::vector vertices_for_this_element( + n_vertices_per_element); + for (unsigned int e = 0; e < n_elements; ++e) + { + AssertThrow(whole_file, ExcIO()); + for (unsigned int v = 0; v < n_vertices_per_element; ++v) + { + whole_file >> vertices_for_this_element[v]; + vertices_for_this_element[v] -= starting_vertex_index; + } + + if (object_type == "vtx") + ; // do nothing + else if (object_type == "edg") + { + if (spacedim == 1) + cells[e].vertices = vertices_for_this_element; + else + subcelldata.boundary_lines[e].vertices = + vertices_for_this_element; + } + else if (object_type == "tri") + { + if (spacedim == 2) + cells[e].vertices = vertices_for_this_element; + else + subcelldata.boundary_quads[e].vertices = + vertices_for_this_element; + } + else if (object_type == "tet") + { + if (spacedim == 3) + cells[e].vertices = vertices_for_this_element; + else + Assert(false, ExcInternalError()); + } + else + Assert(false, ExcNotImplemented()); + } + + // Then also read the "geometric entity indices". There need to be as + // many as there were elements to begin with + { + unsigned int dummy; + whole_file >> dummy; + AssertThrow(dummy == n_elements, ExcInternalError()); + } + + for (unsigned int e = 0; e < n_elements; ++e) + { + AssertThrow(whole_file, ExcIO()); + unsigned int geometric_entity_index; + whole_file >> geometric_entity_index; + if (object_type == "vtx") + ; // do nothing + else if (object_type == "edg") + { + if (spacedim == 1) + cells[e].boundary_id = geometric_entity_index; + else + subcelldata.boundary_lines[e].boundary_id = + geometric_entity_index; + } + else if (object_type == "tri") + { + if (spacedim == 2) + cells[e].boundary_id = geometric_entity_index; + else + subcelldata.boundary_quads[e].boundary_id = + geometric_entity_index; + } + else if (object_type == "tet") + { + if (spacedim == 3) + cells[e].boundary_id = geometric_entity_index; + else + Assert(false, ExcInternalError()); + } + else + Assert(false, ExcNotImplemented()); + } + } + AssertThrow(whole_file, ExcIO()); + + // Now finally create the mesh. Because of the quirk with boundary + // edges and faces described in the documentation of this function, + // we can't pass 'subcelldata' as third argument to this function. + // Rather, we then have to fix up the generated triangulation + // after the fact :-( + tria->create_triangulation(vertices, cells, {}); + + // Now for the "fixing up" step mentioned above. To make things a bit + // simpler, let us sort first normalize the order of vertices in edges + // and triangles/quads, and then sort lexicographically: + if (dim >= 2) + { + for (auto &line : subcelldata.boundary_lines) + { + Assert(line.vertices.size() == 2, ExcInternalError()); + if (line.vertices[1] < line.vertices[0]) + std::swap(line.vertices[0], line.vertices[1]); + } + std::sort(subcelldata.boundary_lines.begin(), + subcelldata.boundary_lines.end(), + [](const CellData<1> &a, const CellData<1> &b) { + return std::lexicographical_compare(a.vertices.begin(), + a.vertices.end(), + b.vertices.begin(), + b.vertices.end()); + }); + } + + // Now for boundary faces. For triangles, we can sort the vertices in + // ascending vertex index order because every order corresponds to a circular + // order either seen from one side or the other. For quads, the situation is + // more difficult. But fortunately, we do not actually need to keep the + // vertices in any specific order because there can be no two quads with the + // same vertices but listed in different orders that actually correspond to + // different things. If we had given this information to + // Triangulation::create_triangulation(), we would probably have wanted to + // keep things in a specific order so that the vertices define a proper + // coordinate system on the quad, but that's not our goal here so we just + // sort. + if (dim >= 3) + { + for (auto &face : subcelldata.boundary_quads) + { + Assert((face.vertices.size() == 3) || (face.vertices.size() == 4), + ExcInternalError()); + std::sort(face.vertices.begin(), face.vertices.end()); + } + std::sort(subcelldata.boundary_quads.begin(), + subcelldata.boundary_quads.end(), + [](const CellData<2> &a, const CellData<2> &b) { + return std::lexicographical_compare(a.vertices.begin(), + a.vertices.end(), + b.vertices.begin(), + b.vertices.end()); + }); + } + + // OK, now we can finally go about fixing up edges and faces. + if (dim >= 2) + { + for (const auto &cell : tria->active_cell_iterators()) + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + // We found a face at the boundary. Let us look up whether it + // was listed in subcelldata + if (dim == 2) + { + std::array face_vertex_indices = { + {face->vertex_index(0), face->vertex_index(1)}}; + if (face_vertex_indices[0] > face_vertex_indices[1]) + std::swap(face_vertex_indices[0], face_vertex_indices[1]); + + // See if we can find an edge with these indices: + const auto p = + std::lower_bound(subcelldata.boundary_lines.begin(), + subcelldata.boundary_lines.end(), + face_vertex_indices, + [](const CellData<1> &a, + const std::array + &face_vertex_indices) -> bool { + return std::lexicographical_compare( + a.vertices.begin(), + a.vertices.end(), + face_vertex_indices.begin(), + face_vertex_indices.end()); + }); + + if ((p != subcelldata.boundary_lines.end()) && + (p->vertices[0] == face_vertex_indices[0]) && + (p->vertices[1] == face_vertex_indices[1])) + { + face->set_boundary_id(p->boundary_id); + } + } + else if (dim == 3) + { + // In 3d, we need to look things up in the boundary_quads + // structure (which also stores boundary triangles) as well as + // for the edges + std::vector face_vertex_indices( + face->n_vertices()); + for (unsigned int v = 0; v < face->n_vertices(); ++v) + face_vertex_indices[v] = face->vertex_index(v); + std::sort(face_vertex_indices.begin(), + face_vertex_indices.end()); + + // See if we can find a face with these indices: + const auto p = + std::lower_bound(subcelldata.boundary_quads.begin(), + subcelldata.boundary_quads.end(), + face_vertex_indices, + [](const CellData<2> &a, + const std::vector + &face_vertex_indices) -> bool { + return std::lexicographical_compare( + a.vertices.begin(), + a.vertices.end(), + face_vertex_indices.begin(), + face_vertex_indices.end()); + }); + + if ((p != subcelldata.boundary_quads.end()) && + (p->vertices == face_vertex_indices)) + { + face->set_boundary_id(p->boundary_id); + } + + + // Now do the same for the edges + for (unsigned int e = 0; e < face->n_lines(); ++e) + { + const auto edge = face->line(e); + + std::array edge_vertex_indices = { + {edge->vertex_index(0), edge->vertex_index(1)}}; + if (edge_vertex_indices[0] > edge_vertex_indices[1]) + std::swap(edge_vertex_indices[0], + edge_vertex_indices[1]); + + // See if we can find an edge with these indices: + const auto p = + std::lower_bound(subcelldata.boundary_lines.begin(), + subcelldata.boundary_lines.end(), + edge_vertex_indices, + [](const CellData<1> &a, + const std::array + &edge_vertex_indices) -> bool { + return std::lexicographical_compare( + a.vertices.begin(), + a.vertices.end(), + edge_vertex_indices.begin(), + edge_vertex_indices.end()); + }); + + if ((p != subcelldata.boundary_lines.end()) && + (p->vertices[0] == edge_vertex_indices[0]) && + (p->vertices[1] == edge_vertex_indices[1])) + { + edge->set_boundary_id(p->boundary_id); + } + } + } + } + } +} + + template void GridIn::read_msh(std::istream &in) @@ -1473,8 +1952,8 @@ GridIn::read_msh(std::istream &in) 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 file format is 2.0 or greater then we also have to read the rest of + // the header if (gmsh_file_format == 20) { double version; @@ -1488,8 +1967,8 @@ GridIn::read_msh(std::istream &in) Assert(file_type == 0, ExcNotImplemented()); Assert(data_size == sizeof(double), ExcNotImplemented()); - // read the end of the header and the first line of the nodes description - // to synch ourselves with the format 1 handling above + // read the end of the header and the first line of the nodes + // description to synch ourselves with the format 1 handling above in >> line; AssertThrow(line == "$EndMeshFormat", ExcInvalidGMSHInput(line)); @@ -1533,7 +2012,8 @@ GridIn::read_msh(std::istream &in) in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >> box_max_y >> box_max_z >> n_physicals; } - // if there is a physical tag, we will use it as boundary id below + // if there is a physical tag, we will use it as boundary id + // below AssertThrow(n_physicals < 2, ExcMessage("More than one tag is not supported!")); // if there is no physical tag, use 0 as default @@ -1553,7 +2033,8 @@ GridIn::read_msh(std::istream &in) // we only care for 'tag' as key for tag_maps[1] in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >> box_max_y >> box_max_z >> n_physicals; - // if there is a physical tag, we will use it as boundary id below + // if there is a physical tag, we will use it as boundary id + // below AssertThrow(n_physicals < 2, ExcMessage("More than one tag is not supported!")); // if there is no physical tag, use 0 as default @@ -1561,8 +2042,9 @@ GridIn::read_msh(std::istream &in) for (unsigned int j = 0; j < n_physicals; ++j) in >> physical_tag; tag_maps[1][tag] = physical_tag; - // we don't care about the points associated to a curve, but have - // to parse them anyway because their format is unstructured + // we don't care about the points associated to a curve, but + // have to parse them anyway because their format is + // unstructured in >> n_points; for (unsigned int j = 0; j < n_points; ++j) in >> tag; @@ -1579,7 +2061,8 @@ GridIn::read_msh(std::istream &in) // we only care for 'tag' as key for tag_maps[2] in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >> box_max_y >> box_max_z >> n_physicals; - // if there is a physical tag, we will use it as boundary id below + // if there is a physical tag, we will use it as boundary id + // below AssertThrow(n_physicals < 2, ExcMessage("More than one tag is not supported!")); // if there is no physical tag, use 0 as default @@ -1588,7 +2071,8 @@ GridIn::read_msh(std::istream &in) in >> physical_tag; tag_maps[2][tag] = physical_tag; // we don't care about the curves associated to a surface, but - // have to parse them anyway because their format is unstructured + // have to parse them anyway because their format is + // unstructured in >> n_curves; for (unsigned int j = 0; j < n_curves; ++j) in >> tag; @@ -1604,7 +2088,8 @@ GridIn::read_msh(std::istream &in) // we only care for 'tag' as key for tag_maps[3] in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >> box_max_y >> box_max_z >> n_physicals; - // if there is a physical tag, we will use it as boundary id below + // if there is a physical tag, we will use it as boundary id + // below AssertThrow(n_physicals < 2, ExcMessage("More than one tag is not supported!")); // if there is no physical tag, use 0 as default @@ -1613,7 +2098,8 @@ GridIn::read_msh(std::istream &in) in >> physical_tag; tag_maps[3][tag] = physical_tag; // we don't care about the surfaces associated to a volume, but - // have to parse them anyway because their format is unstructured + // have to parse them anyway because their format is + // unstructured in >> n_surfaces; for (unsigned int j = 0; j < n_surfaces; ++j) in >> tag; @@ -1755,8 +2241,8 @@ GridIn::read_msh(std::istream &in) } // 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 + // 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; @@ -1804,21 +2290,22 @@ GridIn::read_msh(std::istream &in) 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 + 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 + 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. + 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; @@ -2093,8 +2580,9 @@ GridIn::read_msh(std::istream &in) in >> node_index; } - // we only care about boundary indicators assigned to individual - // vertices in 1d (because otherwise the vertices are not faces) + // 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; } @@ -2271,8 +2759,8 @@ GridIn::read_msh(const std::string &fname) } // If we didn't find a BoundaryID:XX or MaterialID:XX, and // something was found but not recognized, then we set the - // material id or boundary id in the catch block below, using - // directly the physical tag + // material id or boundary id in the catch block below, + // using directly the physical tag if (throw_anyway && !found_boundary_id) throw; } @@ -2999,16 +3487,16 @@ GridIn::read_assimp(const std::string &filename, // Namespace containing some extra functions for reading ExodusII files namespace { - // Convert ExodusII strings to cell types. Use the number of nodes per element - // to disambiguate some cases. + // Convert ExodusII strings to cell types. Use the number of nodes per + // element to disambiguate some cases. ReferenceCell exodusii_name_to_type(const std::string &type_name, const int n_nodes_per_element) { Assert(type_name.size() > 0, ExcInternalError()); - // Try to canonify the name by switching to upper case and removing trailing - // numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and PYRAMID13 all - // equal. + // Try to canonify the name by switching to upper case and removing + // trailing numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and + // PYRAMID13 all equal. std::string type_name_2 = type_name; std::transform(type_name_2.begin(), type_name_2.end(), @@ -3060,10 +3548,11 @@ namespace std::vector> b_or_m_id_to_sideset_ids; // boundary id 0 is the default b_or_m_id_to_sideset_ids.emplace_back(); - // deal.II does not support assigning boundary ids with nonzero codimension - // meshes so completely skip this information in that case. + // deal.II does not support assigning boundary ids with nonzero + // codimension meshes so completely skip this information in that case. // - // Exodus prints warnings if we try to get empty sets so always check first + // Exodus prints warnings if we try to get empty sets so always check + // first if (dim == spacedim && n_side_sets > 0) { std::vector side_set_ids(n_side_sets); @@ -3071,9 +3560,10 @@ namespace AssertThrowExodusII(ierr); // First collect all side sets on all boundary faces (indexed here as - // max_faces_per_cell * cell_n + face_n). We then sort and uniquify the - // side sets so that we can convert a set of side set indices into a - // single deal.II boundary or manifold id (and save the correspondence). + // max_faces_per_cell * cell_n + face_n). We then sort and uniquify + // the side sets so that we can convert a set of side set indices into + // a single deal.II boundary or manifold id (and save the + // correspondence). constexpr auto max_faces_per_cell = GeometryInfo::faces_per_cell; std::map> face_side_sets; for (const int side_set_id : side_set_ids) @@ -3100,10 +3590,10 @@ namespace // According to the manual (subsection 4.8): "The internal // number of an element numbering is defined implicitly by the - // order in which it appears in the file. Elements are numbered - // internally (beginning with 1) consecutively across all - // element blocks." Hence element i in Exodus numbering is entry - // i - 1 in the cells array. + // order in which it appears in the file. Elements are + // numbered internally (beginning with 1) consecutively across + // all element blocks." Hence element i in Exodus numbering is + // entry i - 1 in the cells array. for (int side_n = 0; side_n < n_sides; ++side_n) { const long element_n = elements[side_n] - 1; @@ -3141,8 +3631,9 @@ namespace const std::vector &face_sideset_ids = pair.second; if (face_sideset_ids != b_or_m_id_to_sideset_ids.back()) { - // Since we sorted by sideset ids we are guaranteed that if this - // doesn't match the last set then it has not yet been seen + // Since we sorted by sideset ids we are guaranteed that if + // this doesn't match the last set then it has not yet been + // seen ++current_b_or_m_id; b_or_m_id_to_sideset_ids.push_back(face_sideset_ids); Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1, @@ -3158,9 +3649,9 @@ namespace const ReferenceCell face_reference_cell = cell_type.face_reference_cell(deal_face_n); - // The orientation we pick doesn't matter here since when we create - // the Triangulation we will sort the vertices for each CellData - // object created here. + // The orientation we pick doesn't matter here since when we + // create the Triangulation we will sort the vertices for each + // CellData object created here. if (dim == 2) { CellData<1> boundary_line(face_reference_cell.n_vertices()); @@ -3251,9 +3742,9 @@ GridIn::read_exodusii( // Even if there is a node numbering array the values stored inside the // ExodusII file must use the contiguous, internal ordering (see Section 4.5 - // of the manual - "Internal (contiguously numbered) node and element IDs must - // be used for all data structures that contain node or element numbers (IDs), - // including node set node lists, side set element lists, and element + // of the manual - "Internal (contiguously numbered) node and element IDs + // must be used for all data structures that contain node or element numbers + // (IDs), including node set node lists, side set element lists, and element // connectivity.") std::vector> vertices; vertices.reserve(n_nodes); @@ -3961,8 +4452,8 @@ namespace ::toupper); // Surface can be created from ELSET, or directly from cells - // If elsets_list contains a key with specific name - refers to - // that ELSET, otherwise refers to cell + // If elsets_list contains a key with specific name - refers + // to that ELSET, otherwise refers to cell std::istringstream iss(line); int el_idx; int face_number; @@ -4043,7 +4534,8 @@ namespace int elis_step = 1; // Default if case stride not provided // Some files don't have the stride size - // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp to + // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp + // to // ./grids/abaqus/2d/2d_test_abaqus.inp iss >> elid_start >> comma >> elid_end; AssertThrow(comma == ',', @@ -4121,7 +4613,8 @@ namespace } else if (line.compare(0, 14, "*SOLID SECTION") == 0) { - // The ELSET name, which describes a section for particular material + // The ELSET name, which describes a section for particular + // material const std::string elset_key = "ELSET="; const std::size_t elset_start = line.find("ELSET=") + elset_key.size(); @@ -4264,8 +4757,8 @@ namespace // save old formatting options const boost::io::ios_base_all_saver formatting_saver(output); - // Write out title - Note: No other commented text can be inserted below the - // title in a UCD file + // Write out title - Note: No other commented text can be inserted below + // the title in a UCD file output << "# Abaqus to UCD mesh conversion" << std::endl; output << "# Mesh type: AVS UCD" << std::endl; @@ -4275,21 +4768,25 @@ namespace // Comments, if present, must precede all data in the file. // Comments within the data will cause read errors. // The general order of the data is as follows: - // 1. Numbers defining the overall structure, including the number of nodes, + // 1. Numbers defining the overall structure, including the number of + // nodes, // the number of cells, and the length of the vector of data associated // with the nodes, cells, and the model. // e.g. 1: // // e.g. 2: - // n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges - // outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n") - // 2. For each node, its node id and the coordinates of that node in space. + // n_elements = n_hex_cells + n_bc_quads + n_quad_cells + + // n_bc_edges outfile.write(str(n_nodes) + " " + str(n_elements) + + // " 0 0 0\n") + // 2. For each node, its node id and the coordinates of that node in + // space. // Node-ids must be integers, but any number including non sequential // numbers can be used. Mid-edge nodes are treated like any other node. - // 3. For each cell: its cell-id, material, cell type (hexahedral, pyramid, - // etc.), and the list of node-ids that correspond to each of the cell's - // vertices. The below table specifies the different cell types and the - // keyword used to represent them in the file. + // 3. For each cell: its cell-id, material, cell type (hexahedral, + // pyramid, + // etc.), and the list of node-ids that correspond to each of the + // cell's vertices. The below table specifies the different cell types + // and the keyword used to represent them in the file. // Write out header output << node_list.size() << "\t" << (cell_list.size() + face_list.size())