24#include <boost/algorithm/string.hpp>
25#include <boost/archive/binary_iarchive.hpp>
26#include <boost/io/ios_state.hpp>
27#include <boost/property_tree/ptree.hpp>
28#include <boost/property_tree/xml_parser.hpp>
29#include <boost/serialization/serialization.hpp>
31#ifdef DEAL_II_GMSH_WITH_API
42#ifdef DEAL_II_WITH_ASSIMP
43# include <assimp/Importer.hpp>
44# include <assimp/postprocess.h>
45# include <assimp/scene.h>
48#ifdef DEAL_II_TRILINOS_WITH_SEACAS
75 template <
int spacedim>
77 assign_1d_boundary_ids(
84 if (cell->face(face_no)->at_boundary())
85 for (const auto &pair : boundary_ids)
86 if (cell->face(face_no)->vertex(0) == pair.
first)
88 cell->face(face_no)->set_boundary_id(pair.second);
94 template <
int dim,
int spacedim>
96 assign_1d_boundary_ids(
108 template <
int dim,
int spacedim>
116 const auto n_hypercube_vertices =
117 ReferenceCells::get_hypercube<dim>().n_vertices();
118 bool is_only_hypercube =
true;
120 if (cell.
vertices.size() != n_hypercube_vertices)
122 is_only_hypercube =
false;
130 if (is_only_hypercube)
135template <
int dim,
int spacedim>
137 : tria(nullptr, typeid(*this).name())
138 , default_format(ucd)
143template <
int dim,
int spacedim>
145 : tria(&t, typeid(*this).name())
146 , default_format(ucd)
151template <
int dim,
int spacedim>
160template <
int dim,
int spacedim>
172 text[0] =
"# vtk DataFile Version 3.0";
175 text[3] =
"DATASET UNSTRUCTURED_GRID";
177 for (
unsigned int i = 0; i < 4; ++i)
182 line.compare(text[i]) == 0,
185 "While reading VTK file, failed to find a header line with text <") +
192 std::vector<Point<spacedim>>
vertices;
193 std::vector<CellData<dim>> cells;
202 if (keyword ==
"POINTS")
204 unsigned int n_vertices;
209 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
213 in >> x[0] >> x[1] >> x[2];
216 for (
unsigned int d = 0; d < spacedim; ++d)
224 "While reading VTK file, failed to find POINTS section"));
228 unsigned int n_geometric_objects = 0;
231 if (keyword ==
"CELLS")
234 std::vector<unsigned int> cell_types;
236 std::streampos oldpos = in.tellg();
239 while (in >> keyword)
240 if (keyword ==
"CELL_TYPES")
244 cell_types.resize(n_ints);
246 for (
unsigned int i = 0; i < n_ints; ++i)
255 in >> n_geometric_objects;
260 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
262 unsigned int n_vertices;
266 if (cell_types[count] == 10 || cell_types[count] == 12)
274 cells.emplace_back(n_vertices);
276 for (
unsigned int j = 0; j < n_vertices;
278 in >> cells.back().vertices[j];
282 if (cell_types[count] == 12)
284 std::swap(cells.back().vertices[2],
285 cells.back().vertices[3]);
286 std::swap(cells.back().vertices[6],
287 cells.back().vertices[7]);
290 cells.back().material_id = 0;
293 else if (cell_types[count] == 5 || cell_types[count] == 9)
302 for (
unsigned int j = 0; j < n_vertices;
309 else if (cell_types[count] == 3)
313 for (
unsigned int j = 0; j < n_vertices;
324 "While reading VTK file, unknown cell type encountered"));
329 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
331 unsigned int n_vertices;
335 if (cell_types[count] == 5 || cell_types[count] == 9)
342 cells.emplace_back(n_vertices);
344 for (
unsigned int j = 0; j < n_vertices;
346 in >> cells.back().vertices[j];
350 if (cell_types[count] == 9)
354 std::swap(cells.back().vertices[2],
355 cells.back().vertices[3]);
358 cells.back().material_id = 0;
361 else if (cell_types[count] == 3)
367 for (
unsigned int j = 0; j < n_vertices;
380 "While reading VTK file, unknown cell type encountered"));
385 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
391 cell_types[count] == 3 && type == 2,
393 "While reading VTK file, unknown cell type encountered"));
394 cells.emplace_back(type);
396 for (
unsigned int j = 0; j < type; ++j)
397 in >> cells.back().vertices[j];
399 cells.back().material_id = 0;
405 "While reading VTK file, failed to find CELLS section"));
412 keyword ==
"CELL_TYPES",
414 "While reading VTK file, missing CELL_TYPES section. Found <" +
415 keyword +
"> instead.")));
419 n_ints == n_geometric_objects,
420 ExcMessage(
"The VTK reader found a CELL_DATA statement "
421 "that lists a total of " +
423 " cell data objects, but this needs to "
424 "equal the number of cells (which is " +
426 ") plus the number of quads (" +
428 " in 3d or the number of lines (" +
433 for (
unsigned int i = 0; i < n_ints; ++i)
437 while (in >> keyword)
438 if (keyword ==
"CELL_DATA")
444 ExcMessage(
"The VTK reader found a CELL_DATA statement "
445 "that lists a total of " +
447 " cell data objects, but this needs to "
448 "equal the number of cells (which is " +
450 ") plus the number of quads (" +
453 " in 3d or the number of lines (" +
458 const std::vector<std::string> data_sets{
"MaterialID",
461 for (
unsigned int i = 0; i < data_sets.size(); ++i)
464 while (in >> keyword)
465 if (keyword ==
"SCALARS")
470 std::string field_name;
472 if (std::find(data_sets.begin(),
474 field_name) == data_sets.end())
486 std::getline(in, line);
489 std::min(
static_cast<std::size_t
>(3),
490 line.size() - 1)) ==
"int",
492 "While reading VTK file, material- and manifold IDs can only have type 'int'."));
496 keyword ==
"LOOKUP_TABLE",
498 "While reading VTK file, missing keyword 'LOOKUP_TABLE'."));
502 keyword ==
"default",
504 "While reading VTK file, missing keyword 'default'."));
511 for (
unsigned int i = 0; i < cells.size(); ++i)
515 if (field_name ==
"MaterialID")
516 cells[i].material_id =
518 else if (field_name ==
"ManifoldID")
519 cells[i].manifold_id =
531 if (field_name ==
"MaterialID")
532 boundary_quad.material_id =
534 else if (field_name ==
"ManifoldID")
535 boundary_quad.manifold_id =
544 if (field_name ==
"MaterialID")
545 boundary_line.material_id =
547 else if (field_name ==
"ManifoldID")
548 boundary_line.manifold_id =
560 if (field_name ==
"MaterialID")
561 boundary_line.material_id =
563 else if (field_name ==
"ManifoldID")
564 boundary_line.manifold_id =
574 apply_grid_fixup_functions(
vertices, cells, subcelldata);
575 tria->create_triangulation(
vertices, cells, subcelldata);
580 "While reading VTK file, failed to find CELLS section"));
585template <
int dim,
int spacedim>
589 namespace pt = boost::property_tree;
591 pt::read_xml(in, tree);
592 auto section = tree.get_optional<std::string>(
"VTKFile.dealiiData");
596 "While reading a VTU file, failed to find dealiiData section. "
597 "Notice that we can only read grid files in .vtu format that "
598 "were created by the deal.II library, using a call to "
599 "GridOut::write_vtu(), where the flag "
600 "GridOutFlags::Vtu::serialize_triangulation is set to true."));
604 const auto string_archive =
606 std::istringstream in_stream(string_archive);
607 boost::archive::binary_iarchive ia(in_stream);
612template <
int dim,
int spacedim>
616 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
620 skip_comment_lines(in,
'#');
631 "Expected '-1' before and after a section."));
643 std::getline(in, line);
645 boost::algorithm::trim(line);
646 if (line.compare(
"-1") == 0)
657 std::vector<Point<spacedim>>
vertices;
676 in >> dummy >> dummy >> dummy;
679 in >> x[0] >> x[1] >> x[2];
683 for (
unsigned int d = 0; d < spacedim; ++d)
698 AssertThrow(tmp == 2412, ExcUnknownSectionType(tmp));
700 std::vector<CellData<dim>> cells;
727 in >> type >> dummy >> dummy >> dummy >> dummy;
729 AssertThrow((type == 11) || (type == 44) || (type == 94) || (type == 115),
730 ExcUnknownElementType(type));
732 if ((((type == 44) || (type == 94)) && (dim == 2)) ||
733 ((type == 115) && (dim == 3)))
735 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
736 cells.emplace_back();
741 .vertices[reference_cell.unv_vertex_to_deal_vertex(v)];
743 cells.back().material_id = 0;
746 cells.back().vertices[v] =
vertex_indices[cells.back().vertices[v]];
748 cell_indices[object_index] = n_cells;
752 else if (((type == 11) && (dim == 2)) ||
753 ((type == 11) && (dim == 3)))
756 in >> dummy >> dummy >> dummy;
761 for (
unsigned int &vertex :
767 for (
unsigned int &vertex :
771 line_indices[object_index] = n_lines;
775 else if (((type == 44) || (type == 94)) && (dim == 3))
786 .vertices[reference_cell.unv_vertex_to_deal_vertex(v)];
790 for (
unsigned int &vertex :
794 quad_indices[object_index] = n_quads;
802 "> when running in dim=" +
821 AssertThrow((tmp == 2467) || (tmp == 2477), ExcUnknownSectionType(tmp));
837 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
848 std::getline(in, line);
851 "The line before the line containing an ID has too "
852 "many entries. This is not a valid UNV file."));
854 std::getline(in, line);
856 std::istringstream id_stream(line);
859 !id_stream.fail() && id_stream.eof(),
861 "The given UNV file contains a boundary or material id set to '" +
863 "', which cannot be parsed as a fixed-width integer, whereas "
864 "deal.II only supports integer boundary and material ids. To fix "
865 "this, ensure that all such ids are given integer values."));
868 id <=
long(std::numeric_limits<types::material_id>::max()),
869 ExcMessage(
"The provided integer id '" + std::to_string(
id) +
870 "' is not convertible to either types::material_id nor "
871 "types::boundary_id."));
873 const unsigned int n_lines =
874 (n_entities % 2 == 0) ? (n_entities / 2) : ((n_entities + 1) / 2);
876 for (
unsigned int line = 0; line < n_lines; ++line)
878 unsigned int n_fragments;
880 if (line == n_lines - 1)
881 n_fragments = (n_entities % 2 == 0) ? (2) : (1);
885 for (
unsigned int no_fragment = 0; no_fragment < n_fragments;
889 in >> dummy >> no >> dummy >> dummy;
891 if (cell_indices.count(no) > 0)
892 cells[cell_indices[no]].material_id = id;
894 if (line_indices.count(no) > 0)
898 if (quad_indices.count(no) > 0)
906 apply_grid_fixup_functions(
vertices, cells, subcelldata);
907 tria->create_triangulation(
vertices, cells, subcelldata);
912template <
int dim,
int spacedim>
915 const bool apply_all_indicators_to_manifolds)
917 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
921 skip_comment_lines(in,
'#');
924 unsigned int n_vertices;
925 unsigned int n_cells;
928 in >> n_vertices >> n_cells >> dummy
934 std::vector<Point<spacedim>>
vertices(n_vertices);
940 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
947 in >> vertex_number >> x[0] >> x[1] >> x[2];
950 for (
unsigned int d = 0; d < spacedim; ++d)
959 std::vector<CellData<dim>> cells;
962 for (
unsigned int cell = 0; cell < n_cells; ++cell)
971 std::string cell_type;
975 unsigned int material_id;
981 if (((cell_type ==
"line") && (dim == 1)) ||
982 ((cell_type ==
"quad") && (dim == 2)) ||
983 ((cell_type ==
"hex") && (dim == 3)))
987 cells.emplace_back();
992 Assert(material_id <= std::numeric_limits<types::material_id>::max(),
995 std::numeric_limits<types::material_id>::max()));
1000 if (apply_all_indicators_to_manifolds)
1001 cells.back().manifold_id =
1003 cells.back().material_id = material_id;
1011 cells.back().vertices[i] =
1017 ExcInvalidVertexIndex(cell,
1018 cells.back().vertices[i]));
1023 else if ((cell_type ==
"line") && ((dim == 2) || (dim == 3)))
1031 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
1034 std::numeric_limits<types::boundary_id>::max()));
1043 if (apply_all_indicators_to_manifolds)
1060 for (
unsigned int &vertex :
1068 AssertThrow(
false, ExcInvalidVertexIndex(cell, vertex));
1072 else if ((cell_type ==
"quad") && (dim == 3))
1081 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
1084 std::numeric_limits<types::boundary_id>::max()));
1093 if (apply_all_indicators_to_manifolds)
1110 for (
unsigned int &vertex :
1118 Assert(
false, ExcInvalidVertexIndex(cell, vertex));
1124 AssertThrow(
false, ExcUnknownIdentifier(cell_type));
1129 apply_grid_fixup_functions(
vertices, cells, subcelldata);
1130 tria->create_triangulation(
vertices, cells, subcelldata);
1135 template <
int dim,
int spacedim>
1142 read_in_abaqus(std::istream &in);
1144 write_out_avs_ucd(std::ostream &out)
const;
1147 const double tolerance;
1150 get_global_node_numbers(
const int face_cell_no,
1151 const int face_cell_face_no)
const;
1154 std::vector<std::vector<double>> node_list;
1157 std::vector<std::vector<double>> cell_list;
1159 std::vector<std::vector<double>> face_list;
1162 std::map<std::string, std::vector<int>> elsets_list;
1166template <
int dim,
int spacedim>
1169 const bool apply_all_indicators_to_manifolds)
1171 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1176 Assert((spacedim == 2 && dim == spacedim) ||
1177 (spacedim == 3 && (dim == spacedim || dim == spacedim - 1)),
1183 Abaqus_to_UCD<dim, spacedim> abaqus_to_ucd;
1184 abaqus_to_ucd.read_in_abaqus(in);
1186 std::stringstream in_ucd;
1187 abaqus_to_ucd.write_out_avs_ucd(in_ucd);
1195 read_ucd(in_ucd, apply_all_indicators_to_manifolds);
1197 catch (std::exception &exc)
1199 std::cerr <<
"Exception on processing internal UCD data: " << std::endl
1200 << exc.what() << std::endl;
1205 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1206 "More information is provided in an error message printed above. "
1207 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1208 "listed in the documentation?"));
1215 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1216 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1217 "listed in the documentation?"));
1222template <
int dim,
int spacedim>
1226 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1232 skip_comment_lines(in,
'#');
1238 AssertThrow(line ==
"MeshVersionFormatted 0", ExcInvalidDBMESHInput(line));
1240 skip_empty_lines(in);
1244 AssertThrow(line ==
"Dimension", ExcInvalidDBMESHInput(line));
1245 unsigned int dimension;
1247 AssertThrow(dimension == dim, ExcDBMESHWrongDimension(dimension));
1248 skip_empty_lines(in);
1261 while (getline(in, line), line.find(
"# END") == std::string::npos)
1263 skip_empty_lines(in);
1268 AssertThrow(line ==
"Vertices", ExcInvalidDBMESHInput(line));
1270 unsigned int n_vertices;
1274 std::vector<Point<spacedim>>
vertices(n_vertices);
1275 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1278 for (
unsigned int d = 0; d < dim; ++d)
1285 skip_empty_lines(in);
1291 AssertThrow(line ==
"Edges", ExcInvalidDBMESHInput(line));
1293 unsigned int n_edges;
1295 for (
unsigned int edge = 0; edge < n_edges; ++edge)
1298 in >> dummy >> dummy;
1304 skip_empty_lines(in);
1313 AssertThrow(line ==
"CrackedEdges", ExcInvalidDBMESHInput(line));
1316 for (
unsigned int edge = 0; edge < n_edges; ++edge)
1319 in >> dummy >> dummy;
1325 skip_empty_lines(in);
1331 AssertThrow(line ==
"Quadrilaterals", ExcInvalidDBMESHInput(line));
1333 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
1334 {0, 1, 5, 4, 2, 3, 7, 6}};
1335 std::vector<CellData<dim>> cells;
1337 unsigned int n_cells;
1339 for (
unsigned int cell = 0; cell < n_cells; ++cell)
1343 cells.emplace_back();
1347 cells.back().vertices[dim == 3 ? local_vertex_numbering[i] :
1351 (
static_cast<unsigned int>(cells.back().vertices[i]) <=
1353 ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
1355 --cells.back().vertices[i];
1363 skip_empty_lines(in);
1371 while (getline(in, line), ((line.find(
"End") == std::string::npos) && (in)))
1377 apply_grid_fixup_functions(
vertices, cells, subcelldata);
1378 tria->create_triangulation(
vertices, cells, subcelldata);
1383template <
int dim,
int spacedim>
1387 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1390 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1394 std::getline(in, line);
1396 unsigned int n_vertices;
1397 unsigned int n_cells;
1401 std::getline(in, line);
1404 std::getline(in, line);
1407 for (
unsigned int i = 0; i < 8; ++i)
1408 std::getline(in, line);
1411 std::vector<CellData<dim>> cells(n_cells);
1422 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1423 in >> cell.vertices[reference_cell.exodusii_vertex_to_deal_vertex(i)];
1427 std::vector<Point<spacedim>>
vertices(n_vertices);
1430 for (
unsigned int d = 0; d < spacedim; ++d)
1432 for (
unsigned int d = spacedim; d < 3; ++d)
1441 apply_grid_fixup_functions(
vertices, cells, subcelldata);
1442 tria->create_triangulation(
vertices, cells, subcelldata);
1447template <
int dim,
int spacedim>
1451 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1464 std::stringstream whole_file;
1469 std::getline(in, line);
1477 if ((line.size() > 0) && (line.back() ==
'\r'))
1478 line.erase(line.size() - 1);
1483 if (line.find(
'#') != std::string::npos)
1484 line.erase(line.find(
'#'), std::string::npos);
1485 while ((line.size() > 0) && (line.back() ==
' '))
1486 line.erase(line.size() - 1);
1488 if (line.size() > 0)
1489 whole_file <<
'\n' << line;
1508 unsigned int version_major, version_minor;
1509 whole_file >> version_major >> version_minor;
1510 AssertThrow((version_major == 0) && (version_minor == 1),
1511 ExcMessage(
"deal.II can currently only read version 0.1 "
1512 "of the mphtxt file format."));
1517 unsigned int n_tags;
1518 whole_file >> n_tags;
1519 for (
unsigned int i = 0; i < n_tags; ++i)
1522 while (whole_file.peek() ==
'\n')
1524 std::getline(whole_file, dummy);
1530 unsigned int n_types;
1531 whole_file >> n_types;
1532 for (
unsigned int i = 0; i < n_types; ++i)
1535 while (whole_file.peek() ==
'\n')
1537 std::getline(whole_file, dummy);
1557 whole_file >> dummy >> dummy >> dummy;
1561 while (whole_file.peek() ==
'\n')
1563 std::getline(whole_file, s);
1565 ExcMessage(
"Expected '4 Mesh', but got '" + s +
"'."));
1568 unsigned int version;
1569 whole_file >> version;
1573 unsigned int file_space_dim;
1574 whole_file >> file_space_dim;
1578 "The mesh file uses a different number of space dimensions "
1579 "than the triangulation you want to read it into."));
1581 unsigned int n_vertices;
1582 whole_file >> n_vertices;
1584 unsigned int starting_vertex_index;
1585 whole_file >> starting_vertex_index;
1587 std::vector<Point<spacedim>>
vertices(n_vertices);
1588 for (
unsigned int v = 0; v < n_vertices; ++v)
1619 std::vector<CellData<dim>> cells;
1622 unsigned int n_types;
1623 whole_file >> n_types;
1624 for (
unsigned int type = 0; type < n_types; ++type)
1631 whole_file >> dummy;
1635 std::string object_name;
1636 whole_file >> object_name;
1638 static const std::map<std::string, ReferenceCell> name_to_type = {
1648 AssertThrow(name_to_type.find(object_name) != name_to_type.end(),
1649 ExcMessage(
"The input file contains a cell type <" +
1651 "> that the reader does not "
1652 "current support"));
1653 const ReferenceCell object_type = name_to_type.at(object_name);
1655 unsigned int n_vertices_per_element;
1656 whole_file >> n_vertices_per_element;
1658 unsigned int n_elements;
1659 whole_file >> n_elements;
1673 ExcMessage(
"Triangles should not appear in input files "
1681 "Quadrilaterals should not appear in input files "
1688 ExcMessage(
"Tetrahedra should not appear in input files "
1689 "for 1d or 2d meshes."));
1696 "Prisms (wedges) should not appear in input files "
1697 "for 1d or 2d meshes."));
1716 std::vector<unsigned int> vertices_for_this_element(
1717 n_vertices_per_element);
1718 for (
unsigned int e = 0; e < n_elements; ++e)
1721 for (
unsigned int v = 0; v < n_vertices_per_element; ++v)
1723 whole_file >> vertices_for_this_element[v];
1724 vertices_for_this_element[v] -= starting_vertex_index;
1733 cells.emplace_back();
1734 cells.back().vertices = vertices_for_this_element;
1740 vertices_for_this_element;
1748 cells.emplace_back();
1749 cells.back().vertices = vertices_for_this_element;
1755 vertices_for_this_element;
1763 cells.emplace_back();
1764 cells.back().vertices = vertices_for_this_element;
1777 unsigned int n_geom_entity_indices;
1778 whole_file >> n_geom_entity_indices;
1780 (n_geom_entity_indices == n_elements),
1787 if (n_geom_entity_indices != 0)
1789 for (
unsigned int e = 0; e < n_geom_entity_indices; ++e)
1792 unsigned int geometric_entity_index;
1793 whole_file >> geometric_entity_index;
1799 cells[cells.size() - n_elements + e].material_id =
1800 geometric_entity_index;
1805 .boundary_id = geometric_entity_index;
1811 cells[cells.size() - n_elements + e].material_id =
1812 geometric_entity_index;
1817 .boundary_id = geometric_entity_index;
1823 cells[cells.size() - n_elements + e].material_id =
1824 geometric_entity_index;
1840 tria->create_triangulation(
vertices, cells, {});
1850 if (line.vertices[1] < line.vertices[0])
1851 std::swap(line.vertices[0], line.vertices[1]);
1856 return std::lexicographical_compare(a.vertices.begin(),
1878 Assert((face.vertices.size() == 3) || (face.vertices.size() == 4),
1880 std::sort(face.vertices.begin(), face.vertices.end());
1885 return std::lexicographical_compare(a.vertices.begin(),
1895 for (
const auto &cell : tria->active_cell_iterators())
1896 for (
const auto &face : cell->face_iterators())
1897 if (face->at_boundary())
1903 std::array<unsigned int, 2> face_vertex_indices = {
1904 {face->vertex_index(0), face->vertex_index(1)}};
1905 if (face_vertex_indices[0] > face_vertex_indices[1])
1906 std::swap(face_vertex_indices[0], face_vertex_indices[1]);
1912 face_vertex_indices,
1914 const std::array<unsigned int, 2>
1915 &face_vertex_indices) ->
bool {
1916 return std::lexicographical_compare(
1919 face_vertex_indices.begin(),
1920 face_vertex_indices.end());
1924 (p->vertices[0] == face_vertex_indices[0]) &&
1925 (p->vertices[1] == face_vertex_indices[1]))
1927 face->set_boundary_id(p->boundary_id);
1935 std::vector<unsigned int> face_vertex_indices(
1936 face->n_vertices());
1937 for (
unsigned int v = 0; v < face->n_vertices(); ++v)
1938 face_vertex_indices[v] = face->vertex_index(v);
1939 std::sort(face_vertex_indices.begin(),
1940 face_vertex_indices.end());
1946 face_vertex_indices,
1948 const std::vector<unsigned int>
1949 &face_vertex_indices) ->
bool {
1950 return std::lexicographical_compare(
1953 face_vertex_indices.begin(),
1954 face_vertex_indices.end());
1958 (p->vertices == face_vertex_indices))
1960 face->set_boundary_id(p->boundary_id);
1965 for (
unsigned int e = 0; e < face->n_lines(); ++e)
1967 const auto edge = face->line(e);
1969 std::array<unsigned int, 2> edge_vertex_indices = {
1970 {edge->vertex_index(0), edge->vertex_index(1)}};
1971 if (edge_vertex_indices[0] > edge_vertex_indices[1])
1972 std::swap(edge_vertex_indices[0],
1973 edge_vertex_indices[1]);
1979 edge_vertex_indices,
1981 const std::array<unsigned int, 2>
1982 &edge_vertex_indices) ->
bool {
1983 return std::lexicographical_compare(
1986 edge_vertex_indices.begin(),
1987 edge_vertex_indices.end());
1991 (p->vertices[0] == edge_vertex_indices[0]) &&
1992 (p->vertices[1] == edge_vertex_indices[1]))
1994 edge->set_boundary_id(p->boundary_id);
2003template <
int dim,
int spacedim>
2007 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
2010 unsigned int n_vertices;
2011 unsigned int n_cells;
2017 std::array<std::map<int, int>, 4> tag_maps;
2020 std::string stripped_file;
2024 while (std::getline(input_stream, line))
2026 if (line ==
"@f$Comments")
2028 while (std::getline(input_stream, line))
2030 if (line ==
"@f$EndComments")
2037 stripped_file += line +
'\n';
2041 std::istringstream in(stripped_file);
2046 unsigned int gmsh_file_format = 0;
2047 if (line ==
"@f$NOD")
2048 gmsh_file_format = 10;
2049 else if (line ==
"@f$MeshFormat")
2050 gmsh_file_format = 20;
2056 if (gmsh_file_format == 20)
2059 unsigned int file_type, data_size;
2061 in >> version >> file_type >> data_size;
2064 gmsh_file_format =
static_cast<unsigned int>(version * 10);
2072 AssertThrow(line ==
"@f$EndMeshFormat", ExcInvalidGMSHInput(line));
2076 if (line ==
"@f$PhysicalNames")
2082 while (line !=
"@f$EndPhysicalNames");
2087 if (line ==
"@f$Entities")
2089 unsigned long n_points, n_curves, n_surfaces, n_volumes;
2091 in >> n_points >> n_curves >> n_surfaces >> n_volumes;
2092 for (
unsigned int i = 0; i < n_points; ++i)
2096 unsigned int n_physicals;
2097 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2101 if (gmsh_file_format > 40)
2103 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2105 box_max_x = box_min_x;
2106 box_max_y = box_min_y;
2107 box_max_z = box_min_z;
2111 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2112 box_max_x >> box_max_y >> box_max_z >> n_physicals;
2117 ExcMessage(
"More than one tag is not supported!"));
2119 int physical_tag = 0;
2120 for (
unsigned int j = 0; j < n_physicals; ++j)
2122 tag_maps[0][tag] = physical_tag;
2124 for (
unsigned int i = 0; i < n_curves; ++i)
2128 unsigned int n_physicals;
2129 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2133 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2134 box_max_y >> box_max_z >> n_physicals;
2138 ExcMessage(
"More than one tag is not supported!"));
2140 int physical_tag = 0;
2141 for (
unsigned int j = 0; j < n_physicals; ++j)
2143 tag_maps[1][tag] = physical_tag;
2148 for (
unsigned int j = 0; j < n_points; ++j)
2152 for (
unsigned int i = 0; i < n_surfaces; ++i)
2156 unsigned int n_physicals;
2157 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2161 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2162 box_max_y >> box_max_z >> n_physicals;
2166 ExcMessage(
"More than one tag is not supported!"));
2168 int physical_tag = 0;
2169 for (
unsigned int j = 0; j < n_physicals; ++j)
2171 tag_maps[2][tag] = physical_tag;
2176 for (
unsigned int j = 0; j < n_curves; ++j)
2179 for (
unsigned int i = 0; i < n_volumes; ++i)
2183 unsigned int n_physicals;
2184 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2188 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2189 box_max_y >> box_max_z >> n_physicals;
2193 ExcMessage(
"More than one tag is not supported!"));
2195 int physical_tag = 0;
2196 for (
unsigned int j = 0; j < n_physicals; ++j)
2198 tag_maps[3][tag] = physical_tag;
2203 for (
unsigned int j = 0; j < n_surfaces; ++j)
2207 AssertThrow(line ==
"@f$EndEntities", ExcInvalidGMSHInput(line));
2212 if (line ==
"@f$PartitionedEntities")
2218 while (line !=
"@f$EndPartitionedEntities");
2225 AssertThrow(line ==
"@f$Nodes", ExcInvalidGMSHInput(line));
2229 int n_entity_blocks = 1;
2230 if (gmsh_file_format > 40)
2234 in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag;
2236 else if (gmsh_file_format == 40)
2238 in >> n_entity_blocks >> n_vertices;
2242 std::vector<Point<spacedim>>
vertices(n_vertices);
2249 unsigned int global_vertex = 0;
2250 for (
int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2253 unsigned long numNodes;
2255 if (gmsh_file_format < 40)
2257 numNodes = n_vertices;
2264 int tagEntity, dimEntity;
2265 in >> tagEntity >> dimEntity >> parametric >> numNodes;
2268 std::vector<int> vertex_numbers;
2270 if (gmsh_file_format > 40)
2271 for (
unsigned long vertex_per_entity = 0;
2272 vertex_per_entity < numNodes;
2273 ++vertex_per_entity)
2275 in >> vertex_number;
2276 vertex_numbers.push_back(vertex_number);
2279 for (
unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes;
2280 ++vertex_per_entity, ++global_vertex)
2286 if (gmsh_file_format > 40)
2288 vertex_number = vertex_numbers[vertex_per_entity];
2289 in >> x[0] >> x[1] >> x[2];
2292 in >> vertex_number >> x[0] >> x[1] >> x[2];
2294 for (
unsigned int d = 0; d < spacedim; ++d)
2300 if (parametric != 0)
2315 static const std::string end_nodes_marker[] = {
"@f$ENDNOD",
"@f$EndNodes"};
2316 AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1],
2317 ExcInvalidGMSHInput(line));
2321 static const std::string begin_elements_marker[] = {
"@f$ELM",
"@f$Elements"};
2322 AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2323 ExcInvalidGMSHInput(line));
2326 if (gmsh_file_format > 40)
2330 in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag;
2332 else if (gmsh_file_format == 40)
2334 in >> n_entity_blocks >> n_cells;
2338 n_entity_blocks = 1;
2345 std::vector<CellData<dim>> cells;
2347 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2353 std::map<unsigned int, unsigned int> vertex_counts;
2356 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
2357 {0, 1, 5, 4, 2, 3, 7, 6}};
2358 unsigned int global_cell = 0;
2359 for (
int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2361 unsigned int material_id;
2362 unsigned long numElements;
2365 if (gmsh_file_format < 40)
2369 numElements = n_cells;
2371 else if (gmsh_file_format == 40)
2373 int tagEntity, dimEntity;
2374 in >> tagEntity >> dimEntity >> cell_type >> numElements;
2375 material_id = tag_maps[dimEntity][tagEntity];
2380 int tagEntity, dimEntity;
2381 in >> dimEntity >> tagEntity >> cell_type >> numElements;
2382 material_id = tag_maps[dimEntity][tagEntity];
2385 for (
unsigned int cell_per_entity = 0; cell_per_entity < numElements;
2386 ++cell_per_entity, ++global_cell)
2395 unsigned int nod_num;
2416 unsigned int elm_number = 0;
2417 if (gmsh_file_format < 40)
2423 if (gmsh_file_format < 20)
2429 else if (gmsh_file_format < 40)
2434 unsigned int n_tags;
2441 for (
unsigned int i = 1; i < n_tags; ++i)
2446 else if (cell_type == 2)
2448 else if (cell_type == 3)
2450 else if (cell_type == 4)
2452 else if (cell_type == 5)
2463 else if (cell_type == 2)
2465 else if (cell_type == 3)
2467 else if (cell_type == 4)
2469 else if (cell_type == 5)
2495 if (((cell_type == 1) && (dim == 1)) ||
2496 ((cell_type == 2) && (dim == 2)) ||
2497 ((cell_type == 3) && (dim == 2)) ||
2498 ((cell_type == 4) && (dim == 3)) ||
2499 ((cell_type == 5) && (dim == 3)))
2502 unsigned int vertices_per_cell = 0;
2504 vertices_per_cell = 2;
2505 else if (cell_type == 2)
2506 vertices_per_cell = 3;
2507 else if (cell_type == 3)
2508 vertices_per_cell = 4;
2509 else if (cell_type == 4)
2510 vertices_per_cell = 4;
2511 else if (cell_type == 5)
2512 vertices_per_cell = 8;
2516 "Number of nodes does not coincide with the "
2517 "number required for this object"));
2520 cells.emplace_back();
2522 cell.
vertices.resize(vertices_per_cell);
2523 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2526 if (vertices_per_cell ==
2529 local_vertex_numbering[i] :
2537 std::numeric_limits<types::material_id>::max(),
2541 std::numeric_limits<types::material_id>::max()));
2549 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2553 ExcInvalidVertexIndexGmsh(cell_per_entity,
2559 vertex_counts[vertex] += 1u;
2563 else if ((cell_type == 1) &&
2564 ((dim == 2) || (dim == 3)))
2573 std::numeric_limits<types::boundary_id>::max(),
2577 std::numeric_limits<types::boundary_id>::max()));
2588 for (
unsigned int &vertex :
2597 ExcInvalidVertexIndex(cell_per_entity,
2602 else if ((cell_type == 2 || cell_type == 3) &&
2606 unsigned int vertices_per_cell = 0;
2609 vertices_per_cell = 3;
2610 else if (cell_type == 3)
2611 vertices_per_cell = 4;
2619 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2624 std::numeric_limits<types::boundary_id>::max(),
2628 std::numeric_limits<types::boundary_id>::max()));
2639 for (
unsigned int &vertex :
2648 ExcInvalidVertexIndex(cell_per_entity, vertex));
2652 else if (cell_type == 15)
2655 unsigned int node_index = 0;
2656 if (gmsh_file_format < 20)
2676 AssertThrow(
false, ExcGmshUnsupportedGeometry(cell_type));
2684 static const std::string end_elements_marker[] = {
"@f$ENDELM",
"@f$EndElements"};
2685 AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2686 ExcInvalidGMSHInput(line));
2699 for (
const auto &pair : vertex_counts)
2700 if (pair.second == 1u)
2701 boundary_id_pairs.emplace_back(
vertices[pair.first],
2702 boundary_ids_1d[pair.first]);
2704 apply_grid_fixup_functions(
vertices, cells, subcelldata);
2705 tria->create_triangulation(
vertices, cells, subcelldata);
2710 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2715#ifdef DEAL_II_GMSH_WITH_API
2716template <
int dim,
int spacedim>
2720 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
2722 const std::map<int, std::uint8_t> gmsh_to_dealii_type = {
2723 {{15, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {7, 5}, {6, 6}, {5, 7}}};
2726 const std::array<std::vector<unsigned int>, 8> gmsh_to_dealii = {
2733 {{0, 1, 2, 3, 4, 5}},
2734 {{0, 1, 3, 2, 4, 5, 7, 6}}}};
2736 std::vector<Point<spacedim>>
vertices;
2737 std::vector<CellData<dim>> cells;
2739 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2745 std::map<unsigned int, unsigned int> vertex_counts;
2748 gmsh::option::setNumber(
"General.Verbosity", 0);
2752 ExcMessage(
"You are trying to read a gmsh file with dimension " +
2753 std::to_string(gmsh::model::getDimension()) +
2754 " into a grid of dimension " + std::to_string(dim)));
2759 gmsh::model::mesh::removeDuplicateNodes();
2760 gmsh::model::mesh::renumberNodes();
2761 std::vector<std::size_t> node_tags;
2762 std::vector<double> coord;
2763 std::vector<double> parametricCoord;
2764 gmsh::model::mesh::getNodes(
2765 node_tags, coord, parametricCoord, -1, -1,
false,
false);
2767 for (
unsigned int i = 0; i < node_tags.size(); ++i)
2771 for (
unsigned int d = 0; d < spacedim; ++d)
2775 for (
unsigned int d = spacedim; d < 3; ++d)
2777 ExcMessage(
"The grid you are reading contains nodes that are "
2778 "nonzero in the coordinate with index " +
2780 ", but you are trying to save "
2781 "it on a grid embedded in a " +
2782 std::to_string(spacedim) +
" dimensional space."));
2789 std::vector<std::pair<int, int>> entities;
2790 gmsh::model::getEntities(entities);
2792 for (
const auto &e : entities)
2795 const int &entity_dim = e.first;
2796 const int &entity_tag = e.second;
2802 std::vector<int> physical_tags;
2803 gmsh::model::getPhysicalGroupsForEntity(entity_dim,
2808 if (physical_tags.size())
2809 for (
auto physical_tag : physical_tags)
2812 gmsh::model::getPhysicalName(entity_dim, physical_tag, name);
2819 std::map<std::string, int> id_names;
2821 bool found_unrecognized_tag =
false;
2822 bool found_boundary_id =
false;
2825 for (
const auto &it : id_names)
2827 const auto &name = it.first;
2828 const auto &
id = it.second;
2829 if (entity_dim == dim && name ==
"MaterialID")
2832 found_boundary_id =
true;
2834 else if (entity_dim < dim && name ==
"BoundaryID")
2837 found_boundary_id =
true;
2839 else if (name ==
"ManifoldID")
2845 found_unrecognized_tag =
true;
2850 if (found_unrecognized_tag && !found_boundary_id)
2851 boundary_id = physical_tag;
2859 boundary_id = physical_tag;
2865 std::vector<int> element_types;
2866 std::vector<std::vector<std::size_t>> element_ids, element_nodes;
2867 gmsh::model::mesh::getElements(
2868 element_types, element_ids, element_nodes, entity_dim, entity_tag);
2870 for (
unsigned int i = 0; i < element_types.size(); ++i)
2872 const auto &type = gmsh_to_dealii_type.at(element_types[i]);
2873 const auto n_vertices = gmsh_to_dealii[type].size();
2874 const auto &elements = element_ids[i];
2875 const auto &nodes = element_nodes[i];
2876 for (
unsigned int j = 0; j < elements.size(); ++j)
2878 if (entity_dim == dim)
2880 cells.emplace_back(n_vertices);
2881 auto &cell = cells.back();
2882 for (
unsigned int v = 0; v < n_vertices; ++v)
2885 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2887 vertex_counts[cell.vertices[v]] += 1u;
2889 cell.manifold_id = manifold_id;
2890 cell.material_id = boundary_id;
2892 else if (entity_dim == 2)
2896 for (
unsigned int v = 0; v < n_vertices; ++v)
2898 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2900 face.manifold_id = manifold_id;
2901 face.boundary_id = boundary_id;
2903 else if (entity_dim == 1)
2907 for (
unsigned int v = 0; v < n_vertices; ++v)
2909 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2911 line.manifold_id = manifold_id;
2912 line.boundary_id = boundary_id;
2914 else if (entity_dim == 0)
2918 for (
unsigned int j = 0; j < elements.size(); ++j)
2919 boundary_ids_1d[nodes[j] - 1] = boundary_id;
2930 for (
const auto &pair : vertex_counts)
2931 if (pair.second == 1u)
2932 boundary_id_pairs.emplace_back(
vertices[pair.first],
2933 boundary_ids_1d[pair.first]);
2935 apply_grid_fixup_functions(
vertices, cells, subcelldata);
2936 tria->create_triangulation(
vertices, cells, subcelldata);
2941 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2950template <
int dim,
int spacedim>
2953 std::string &header,
2954 std::vector<unsigned int> &tecplot2deal,
2955 unsigned int &n_vars,
2956 unsigned int &n_vertices,
2957 unsigned int &n_cells,
2958 std::vector<unsigned int> &IJK,
2983 std::transform(header.begin(), header.end(), header.begin(), ::toupper);
2987 std::replace(header.begin(), header.end(),
'\t',
' ');
2988 std::replace(header.begin(), header.end(),
',',
' ');
2989 std::replace(header.begin(), header.end(),
'\n',
' ');
2993 std::string::size_type pos = header.find(
'=');
2995 while (pos !=
static_cast<std::string::size_type
>(std::string::npos))
2996 if (header[pos + 1] ==
' ')
2997 header.erase(pos + 1, 1);
2998 else if (header[pos - 1] ==
' ')
3000 header.erase(pos - 1, 1);
3004 pos = header.find(
'=', ++pos);
3007 std::vector<std::string> entries =
3011 for (
unsigned int i = 0; i < entries.size(); ++i)
3020 tecplot2deal[0] = 0;
3023 while (entries[i][0] ==
'"')
3025 if (entries[i] ==
"\"X\"")
3026 tecplot2deal[0] = n_vars;
3027 else if (entries[i] ==
"\"Y\"")
3033 tecplot2deal[1] = n_vars;
3035 else if (entries[i] ==
"\"Z\"")
3041 tecplot2deal[2] = n_vars;
3053 "Tecplot file must contain at least one variable for each dimension"));
3054 for (
unsigned int d = 1; d < dim; ++d)
3056 tecplot2deal[d] > 0,
3058 "Tecplot file must contain at least one variable for each dimension."));
3063 "ZONETYPE=FELINESEG") &&
3067 "ZONETYPE=FEQUADRILATERAL") &&
3071 "ZONETYPE=FEBRICK") &&
3079 "The tecplot file contains an unsupported ZONETYPE."));
3082 "DATAPACKING=POINT"))
3085 "DATAPACKING=BLOCK"))
3108 "ET=QUADRILATERAL") &&
3120 "The tecplot file contains an unsupported ElementType."));
3128 dim > 1 || IJK[1] == 1,
3130 "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
3136 dim > 2 || IJK[2] == 1,
3138 "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
3153 for (
unsigned int d = 0; d < dim; ++d)
3158 "Tecplot file does not contain a complete and consistent set of parameters"));
3159 n_vertices *= IJK[d];
3160 n_cells *= (IJK[d] - 1);
3168 "Tecplot file does not contain a complete and consistent set of parameters"));
3174 n_cells = *std::max_element(IJK.begin(), IJK.end());
3178 "Tecplot file does not contain a complete and consistent set of parameters"));
3188 const unsigned int dim = 2;
3189 const unsigned int spacedim = 2;
3190 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
3194 skip_comment_lines(in,
'#');
3197 std::string line, header;
3204 std::string letters =
"abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
3207 while (line.find_first_of(letters) != std::string::npos)
3209 header +=
" " + line;
3216 std::vector<unsigned int> tecplot2deal(dim);
3217 std::vector<unsigned int> IJK(dim);
3218 unsigned int n_vars, n_vertices, n_cells;
3219 bool structured, blocked;
3221 parse_tecplot_header(header,
3236 std::vector<Point<spacedim>>
vertices(n_vertices + 1);
3239 std::vector<CellData<dim>> cells(n_cells);
3255 unsigned int next_index = 0;
3259 if (tecplot2deal[0] == 0)
3263 std::vector<std::string> first_var =
3266 for (
unsigned int i = 1; i < first_var.size() + 1; ++i)
3267 vertices[i][0] = std::strtod(first_var[i - 1].c_str(), &endptr);
3272 for (
unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
3280 for (
unsigned int i = 1; i < n_vars; ++i)
3289 if (next_index == dim && structured)
3292 if ((next_index < dim) && (i == tecplot2deal[next_index]))
3295 for (
unsigned int j = 1; j < n_vertices + 1; ++j)
3303 for (
unsigned int j = 1; j < n_vertices + 1; ++j)
3315 std::vector<double> vars(n_vars);
3320 std::vector<std::string> first_vertex =
3323 for (
unsigned int d = 0; d < dim; ++d)
3325 std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
3329 for (
unsigned int v = 2; v < n_vertices + 1; ++v)
3331 for (
unsigned int i = 0; i < n_vars; ++i)
3337 for (
unsigned int i = 0; i < dim; ++i)
3338 vertices[v][i] = vars[tecplot2deal[i]];
3346 unsigned int I = IJK[0], J = IJK[1];
3348 unsigned int cell = 0;
3350 for (
unsigned int j = 0; j < J - 1; ++j)
3351 for (
unsigned int i = 1; i < I; ++i)
3353 cells[cell].vertices[0] = i + j * I;
3354 cells[cell].vertices[1] = i + 1 + j * I;
3355 cells[cell].vertices[2] = i + (j + 1) * I;
3356 cells[cell].vertices[3] = i + 1 + (j + 1) * I;
3360 std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
3362 for (
unsigned int i = 1; i < I + 1; ++i)
3364 boundary_vertices[k] = i;
3366 boundary_vertices[k] = i + (J - 1) * I;
3369 for (
unsigned int j = 1; j < J - 1; ++j)
3371 boundary_vertices[k] = 1 + j * I;
3373 boundary_vertices[k] = I + j * I;
3392 for (
unsigned int i = 0; i < n_cells; ++i)
3409 apply_grid_fixup_functions(
vertices, cells, subcelldata);
3410 tria->create_triangulation(
vertices, cells, subcelldata);
3415template <
int dim,
int spacedim>
3424template <
int dim,
int spacedim>
3427 const unsigned int mesh_index,
3428 const bool remove_duplicates,
3430 const bool ignore_unsupported_types)
3432#ifdef DEAL_II_WITH_ASSIMP
3437 Assimp::Importer importer;
3440 const aiScene *scene =
3441 importer.ReadFile(filename.c_str(),
3442 aiProcess_RemoveComponent |
3443 aiProcess_JoinIdenticalVertices |
3444 aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
3445 aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
3451 ExcMessage(
"Input file contains no meshes."));
3454 (mesh_index < scene->mNumMeshes),
3457 unsigned int start_mesh =
3459 unsigned int end_mesh =
3464 std::vector<Point<spacedim>>
vertices;
3465 std::vector<CellData<dim>> cells;
3469 unsigned int v_offset = 0;
3470 unsigned int c_offset = 0;
3472 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
3473 {0, 1, 5, 4, 2, 3, 7, 6}};
3475 for (
unsigned int m = start_mesh; m < end_mesh; ++m)
3477 const aiMesh *mesh = scene->mMeshes[m];
3481 if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
3484 ExcMessage(
"Incompatible mesh " + std::to_string(m) +
3485 "/" + std::to_string(scene->mNumMeshes)));
3488 else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
3491 ExcMessage(
"Incompatible mesh " + std::to_string(m) +
3492 "/" + std::to_string(scene->mNumMeshes)));
3496 const unsigned int n_vertices = mesh->mNumVertices;
3497 const aiVector3D *mVertices = mesh->mVertices;
3500 const unsigned int n_faces = mesh->mNumFaces;
3501 const aiFace *mFaces = mesh->mFaces;
3503 vertices.resize(v_offset + n_vertices);
3504 cells.resize(c_offset + n_faces);
3506 for (
unsigned int i = 0; i < n_vertices; ++i)
3507 for (
unsigned int d = 0; d < spacedim; ++d)
3508 vertices[i + v_offset][d] = mVertices[i][d];
3510 unsigned int valid_cell = c_offset;
3511 for (
unsigned int i = 0; i < n_faces; ++i)
3518 .vertices[dim == 3 ? local_vertex_numbering[f] :
3520 mFaces[i].mIndices[f] + v_offset;
3522 cells[valid_cell].material_id = m;
3528 ExcMessage(
"Face " + std::to_string(i) +
" of mesh " +
3529 std::to_string(m) +
" has " +
3530 std::to_string(mFaces[i].mNumIndices) +
3531 " vertices. We expected only " +
3536 cells.resize(valid_cell);
3541 v_offset += n_vertices;
3542 c_offset = valid_cell;
3549 if (remove_duplicates)
3554 unsigned int n_verts = 0;
3558 std::vector<unsigned int> considered_vertices;
3560 vertices, cells, subcelldata, considered_vertices, tol);
3564 apply_grid_fixup_functions(
vertices, cells, subcelldata);
3565 tria->create_triangulation(
vertices, cells, subcelldata);
3570 (void)remove_duplicates;
3572 (void)ignore_unsupported_types;
3577#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3584 exodusii_name_to_type(
const std::string &type_name,
3585 const int n_nodes_per_element)
3591 std::string type_name_2 = type_name;
3592 std::transform(type_name_2.begin(),
3594 type_name_2.begin(),
3595 [](
unsigned char c) { return std::toupper(c); });
3596 const std::string
numbers =
"0123456789";
3597 type_name_2.erase(std::find_first_of(type_name_2.begin(),
3604 if (type_name_2 ==
"BAR" || type_name_2 ==
"BEAM" ||
3605 type_name_2 ==
"EDGE" || type_name_2 ==
"TRUSS")
3607 else if (type_name_2 ==
"TRI" || type_name_2 ==
"TRIANGLE")
3609 else if (type_name_2 ==
"QUAD" || type_name_2 ==
"QUADRILATERAL")
3611 else if (type_name_2 ==
"SHELL")
3613 if (n_nodes_per_element == 3)
3618 else if (type_name_2 ==
"TET" || type_name_2 ==
"TETRA" ||
3619 type_name_2 ==
"TETRAHEDRON")
3621 else if (type_name_2 ==
"PYRA" || type_name_2 ==
"PYRAMID")
3623 else if (type_name_2 ==
"WEDGE")
3625 else if (type_name_2 ==
"HEX" || type_name_2 ==
"HEXAHEDRON")
3635 template <
int dim,
int spacedim = dim>
3636 std::pair<SubCellData, std::vector<std::vector<int>>>
3637 read_exodusii_sidesets(
const int ex_id,
3638 const int n_side_sets,
3640 const bool apply_all_indicators_to_manifolds)
3643 std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
3645 b_or_m_id_to_sideset_ids.emplace_back();
3651 if (dim == spacedim && n_side_sets > 0)
3653 std::vector<int> side_set_ids(n_side_sets);
3654 int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
3655 AssertThrowExodusII(ierr);
3663 std::map<std::size_t, std::vector<int>> face_side_sets;
3664 for (
const int side_set_id : side_set_ids)
3667 int n_distribution_factors = -1;
3669 ierr = ex_get_set_param(ex_id,
3673 &n_distribution_factors);
3674 AssertThrowExodusII(ierr);
3677 std::vector<int> elements(n_sides);
3678 std::vector<int> faces(n_sides);
3679 ierr = ex_get_set(ex_id,
3684 AssertThrowExodusII(ierr);
3692 for (
int side_n = 0; side_n < n_sides; ++side_n)
3694 const long element_n = elements[side_n] - 1;
3695 const long face_n = faces[side_n] - 1;
3696 const std::size_t face_id =
3697 element_n * max_faces_per_cell + face_n;
3698 face_side_sets[face_id].push_back(side_set_id);
3704 std::vector<std::pair<std::size_t, std::vector<int>>>
3705 face_id_to_side_sets;
3706 face_id_to_side_sets.reserve(face_side_sets.size());
3707 for (
auto &pair : face_side_sets)
3710 face_id_to_side_sets.emplace_back(std::move(pair));
3714 std::sort(face_id_to_side_sets.begin(),
3715 face_id_to_side_sets.end(),
3716 [](
const auto &a,
const auto &b) {
3717 return std::lexicographical_compare(a.second.begin(),
3728 for (
const auto &pair : face_id_to_side_sets)
3730 const std::size_t face_id = pair.first;
3731 const std::vector<int> &face_sideset_ids = pair.second;
3732 if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
3737 ++current_b_or_m_id;
3738 b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
3739 Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
3743 const unsigned int local_face_n = face_id % max_faces_per_cell;
3744 const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
3747 const unsigned int deal_face_n =
3758 if (apply_all_indicators_to_manifolds)
3759 boundary_line.manifold_id = current_b_or_m_id;
3761 boundary_line.boundary_id = current_b_or_m_id;
3762 for (
unsigned int j = 0; j < face_reference_cell.
n_vertices();
3764 boundary_line.vertices[j] =
3766 deal_face_n, j, 0)];
3773 if (apply_all_indicators_to_manifolds)
3774 boundary_quad.manifold_id = current_b_or_m_id;
3776 boundary_quad.boundary_id = current_b_or_m_id;
3777 for (
unsigned int j = 0; j < face_reference_cell.
n_vertices();
3779 boundary_quad.vertices[j] =
3781 deal_face_n, j, 0)];
3788 return std::make_pair(std::move(subcelldata),
3789 std::move(b_or_m_id_to_sideset_ids));
3794template <
int dim,
int spacedim>
3797 const std::string &filename,
3798 const bool apply_all_indicators_to_manifolds)
3800#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3802 int component_word_size =
sizeof(double);
3804 int floating_point_word_size = 0;
3805 float ex_version = 0.0;
3807 const int ex_id = ex_open(filename.c_str(),
3809 &component_word_size,
3810 &floating_point_word_size,
3813 ExcMessage(
"ExodusII failed to open the specified input file."));
3816 std::vector<char> string_temp(MAX_LINE_LENGTH + 1,
'\0');
3817 int mesh_dimension = 0;
3820 int n_element_blocks = 0;
3821 int n_node_sets = 0;
3822 int n_side_sets = 0;
3824 int ierr = ex_get_init(ex_id,
3832 AssertThrowExodusII(ierr);
3843 std::vector<Point<spacedim>>
vertices;
3846 std::vector<double> xs(n_nodes);
3847 std::vector<double> ys(n_nodes);
3848 std::vector<double> zs(n_nodes);
3850 ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
3851 AssertThrowExodusII(ierr);
3853 for (
int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
3858 vertices.emplace_back(xs[vertex_n]);
3861 vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
3864 vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
3872 std::vector<int> element_block_ids(n_element_blocks);
3873 ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
3874 AssertThrowExodusII(ierr);
3876 std::vector<CellData<dim>> cells;
3877 cells.reserve(n_elements);
3881 for (
const int element_block_id : element_block_ids)
3883 std::fill(string_temp.begin(), string_temp.end(),
'\0');
3884 int n_block_elements = 0;
3885 int n_nodes_per_element = 0;
3886 int n_edges_per_element = 0;
3887 int n_faces_per_element = 0;
3888 int n_attributes_per_element = 0;
3891 ierr = ex_get_block(ex_id,
3896 &n_nodes_per_element,
3897 &n_edges_per_element,
3898 &n_faces_per_element,
3899 &n_attributes_per_element);
3900 AssertThrowExodusII(ierr);
3902 exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
3905 "The ExodusII block " + std::to_string(element_block_id) +
3906 " with element type " + std::string(string_temp.data()) +
3908 ", which does not match the topological mesh dimension " +
3909 std::to_string(dim) +
"."));
3916 std::vector<int> connection(n_nodes_per_element * n_block_elements);
3917 ierr = ex_get_conn(ex_id,
3923 AssertThrowExodusII(ierr);
3925 for (
unsigned int elem_n = 0; elem_n < connection.size();
3926 elem_n += n_nodes_per_element)
3932 connection[elem_n + i] - 1;
3935 cells.push_back(std::move(cell));
3940 auto pair = read_exodusii_sidesets<dim, spacedim>(
3941 ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
3942 ierr = ex_close(ex_id);
3943 AssertThrowExodusII(ierr);
3945 apply_grid_fixup_functions(
vertices, cells, pair.first);
3946 tria->create_triangulation(
vertices, cells, pair.first);
3952 (void)apply_all_indicators_to_manifolds;
3959template <
int dim,
int spacedim>
3973 if (std::find_if(line.begin(), line.end(), [](
const char c) {
3978 for (
int i = line.size() - 1; i >= 0; --i)
3979 in.putback(line[i]);
3989template <
int dim,
int spacedim>
3992 const char comment_start)
3997 while (in.get(c) && c == comment_start)
4000 while (in.get() !=
'\n')
4010 skip_empty_lines(in);
4015template <
int dim,
int spacedim>
4033 double min_x =
vertices[cells[0].vertices[0]][0],
4034 max_x =
vertices[cells[0].vertices[0]][0],
4035 min_y =
vertices[cells[0].vertices[0]][1],
4036 max_y =
vertices[cells[0].vertices[0]][1];
4038 for (
unsigned int i = 0; i < cells.size(); ++i)
4040 for (
const auto vertex : cells[i].vertices)
4054 out <<
"# cell " << i << std::endl;
4056 for (
const auto vertex : cells[i].vertices)
4060 out <<
"set label \"" << i <<
"\" at " <<
center[0] <<
',' <<
center[1]
4061 <<
" center" << std::endl;
4064 for (
unsigned int f = 0; f < 2; ++f)
4068 <<
vertices[cells[i].vertices[(f + 1) % 4]][1] << std::endl;
4070 for (
unsigned int f = 2; f < 4; ++f)
4072 <<
',' <<
vertices[cells[i].vertices[(f + 1) % 4]][1] <<
" to "
4080 <<
"set nokey" << std::endl
4081 <<
"pl [" << min_x <<
':' << max_x <<
"][" << min_y <<
':' << max_y
4082 <<
"] " << min_y << std::endl
4083 <<
"pause -1" << std::endl;
4094 for (
const auto &cell : cells)
4161template <
int dim,
int spacedim>
4168 if (format == Default)
4170 const std::string::size_type slashpos = filename.find_last_of(
'/');
4171 const std::string::size_type dotpos = filename.find_last_of(
'.');
4172 if (dotpos < filename.size() &&
4173 (dotpos > slashpos || slashpos == std::string::npos))
4175 std::string ext = filename.substr(dotpos + 1);
4176 format = parse_format(ext);
4180 if (format == assimp)
4182 read_assimp(filename);
4184 else if (format == exodusii)
4186 read_exodusii(filename);
4190 std::ifstream in(filename);
4196template <
int dim,
int spacedim>
4200 if (format == Default)
4201 format = default_format;
4243 ExcMessage(
"There is no read_assimp(istream &) function. "
4244 "Use the read_assimp(string &filename, ...) "
4245 "functions, instead."));
4250 ExcMessage(
"There is no read_exodusii(istream &) function. "
4251 "Use the read_exodusii(string &filename, ...) "
4252 "function, instead."));
4263template <
int dim,
int spacedim>
4292 return ".unknown_format";
4298template <
int dim,
int spacedim>
4302 if (format_name ==
"dbmesh")
4305 if (format_name ==
"exodusii")
4308 if (format_name ==
"msh")
4311 if (format_name ==
"unv")
4314 if (format_name ==
"vtk")
4317 if (format_name ==
"vtu")
4321 if (format_name ==
"inp")
4324 if (format_name ==
"ucd")
4327 if (format_name ==
"xda")
4330 if (format_name ==
"tecplot")
4333 if (format_name ==
"dat")
4336 if (format_name ==
"plt")
4355template <
int dim,
int spacedim>
4359 return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
4366 template <
int dim,
int spacedim>
4367 Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
4380 from_string(T &t,
const std::string &s, std::ios_base &(*f)(std::ios_base &))
4382 std::istringstream iss(s);
4383 return !(iss >> f >> t).fail();
4390 extract_int(
const std::string &s)
4393 for (
const char c : s)
4395 if (isdigit(c) != 0)
4402 from_string(number, tmp, std::dec);
4408 template <
int dim,
int spacedim>
4410 Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
4419 while (std::getline(input_stream, line))
4422 std::transform(line.begin(), line.end(), line.begin(), ::toupper);
4424 if (line.compare(
"*HEADING") == 0 || line.compare(0, 2,
"**") == 0 ||
4425 line.compare(0, 5,
"*PART") == 0)
4428 while (std::getline(input_stream, line))
4434 else if (line.compare(0, 5,
"*NODE") == 0)
4443 while (std::getline(input_stream, line))
4448 std::vector<double> node(spacedim + 1);
4450 std::istringstream iss(line);
4452 for (
unsigned int i = 0; i < spacedim + 1; ++i)
4453 iss >> node[i] >> comma;
4455 node_list.push_back(node);
4458 else if (line.compare(0, 8,
"*ELEMENT") == 0)
4473 const std::string before_material =
"ELSET=EB";
4474 const std::size_t idx = line.find(before_material);
4475 if (idx != std::string::npos)
4477 from_string(material,
4478 line.substr(idx + before_material.size()),
4484 while (std::getline(input_stream, line))
4489 std::istringstream iss(line);
4495 const unsigned int n_data_per_cell =
4497 std::vector<double> cell(n_data_per_cell);
4498 for (
unsigned int i = 0; i < n_data_per_cell; ++i)
4499 iss >> cell[i] >> comma;
4502 cell[0] =
static_cast<double>(material);
4503 cell_list.push_back(cell);
4506 else if (line.compare(0, 8,
"*SURFACE") == 0)
4517 const std::string name_key =
"NAME=";
4518 const std::size_t name_idx_start =
4519 line.find(name_key) + name_key.size();
4520 std::size_t name_idx_end = line.find(
',', name_idx_start);
4521 if (name_idx_end == std::string::npos)
4523 name_idx_end = line.size();
4525 const int b_indicator = extract_int(
4526 line.substr(name_idx_start, name_idx_end - name_idx_start));
4533 while (std::getline(input_stream, line))
4539 std::transform(line.begin(),
4547 std::istringstream iss(line);
4554 std::vector<double> quad_node_list;
4555 const std::string elset_name = line.substr(0, line.find(
','));
4556 if (elsets_list.count(elset_name) != 0)
4560 iss >> stmp >> temp >> face_number;
4562 const std::vector<int> cells = elsets_list[elset_name];
4563 for (
const int cell : cells)
4567 get_global_node_numbers(el_idx, face_number);
4568 quad_node_list.insert(quad_node_list.begin(),
4571 face_list.push_back(quad_node_list);
4578 iss >> el_idx >> comma >> temp >> face_number;
4580 get_global_node_numbers(el_idx, face_number);
4581 quad_node_list.insert(quad_node_list.begin(), b_indicator);
4583 face_list.push_back(quad_node_list);
4587 else if (line.compare(0, 6,
"*ELSET") == 0)
4591 std::string elset_name;
4593 const std::string elset_key =
"*ELSET, ELSET=";
4594 const std::size_t idx = line.find(elset_key);
4595 if (idx != std::string::npos)
4597 const std::string comma =
",";
4598 const std::size_t first_comma = line.find(comma);
4599 const std::size_t second_comma =
4600 line.find(comma, first_comma + 1);
4601 const std::size_t elset_name_start =
4602 line.find(elset_key) + elset_key.size();
4603 elset_name = line.substr(elset_name_start,
4604 second_comma - elset_name_start);
4614 std::vector<int> elements;
4615 const std::size_t generate_idx = line.find(
"GENERATE");
4616 if (generate_idx != std::string::npos)
4619 std::getline(input_stream, line);
4620 std::istringstream iss(line);
4630 iss >> elid_start >> comma >> elid_end;
4634 "While reading an ABAQUS file, the reader "
4635 "expected a comma but found a <") +
4636 comma +
"> in the line <" + line +
">."));
4638 elid_start <= elid_end,
4641 "While reading an ABAQUS file, the reader encountered "
4642 "a GENERATE statement in which the upper bound <") +
4644 "> for the element numbers is not larger or equal "
4645 "than the lower bound <" +
4649 if (iss.rdbuf()->in_avail() != 0)
4650 iss >> comma >> elis_step;
4654 "While reading an ABAQUS file, the reader "
4655 "expected a comma but found a <") +
4656 comma +
"> in the line <" + line +
">."));
4658 for (
int i = elid_start; i <= elid_end; i += elis_step)
4659 elements.push_back(i);
4660 elsets_list[elset_name] = elements;
4662 std::getline(input_stream, line);
4667 while (std::getline(input_stream, line))
4672 std::istringstream iss(line);
4677 iss >> elid >> comma;
4682 "While reading an ABAQUS file, the reader "
4683 "expected a comma but found a <") +
4684 comma +
"> in the line <" + line +
">."));
4686 elements.push_back(elid);
4690 elsets_list[elset_name] = elements;
4695 else if (line.compare(0, 5,
"*NSET") == 0)
4698 while (std::getline(input_stream, line))
4704 else if (line.compare(0, 14,
"*SOLID SECTION") == 0)
4708 const std::string elset_key =
"ELSET=";
4709 const std::size_t elset_start =
4710 line.find(
"ELSET=") + elset_key.size();
4711 const std::size_t elset_end = line.find(
',', elset_start + 1);
4712 const std::string elset_name =
4713 line.substr(elset_start, elset_end - elset_start);
4718 const std::string material_key =
"MATERIAL=";
4719 const std::size_t last_equal =
4720 line.find(
"MATERIAL=") + material_key.size();
4721 const std::size_t material_id_start = line.find(
'-', last_equal);
4723 from_string(material_id,
4724 line.substr(material_id_start + 1),
4728 const std::vector<int> &elset_cells = elsets_list[elset_name];
4729 for (
const int elset_cell : elset_cells)
4731 const int cell_id = elset_cell - 1;
4739 template <
int dim,
int spacedim>
4741 Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
4742 const int face_cell_no,
4743 const int face_cell_face_no)
const
4753 if (face_cell_face_no == 1)
4755 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4756 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4758 else if (face_cell_face_no == 2)
4760 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4761 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4763 else if (face_cell_face_no == 3)
4765 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4766 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4768 else if (face_cell_face_no == 4)
4770 quad_node_list[0] = cell_list[face_cell_no - 1][4];
4771 quad_node_list[1] = cell_list[face_cell_no - 1][1];
4781 if (face_cell_face_no == 1)
4783 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4784 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4785 quad_node_list[2] = cell_list[face_cell_no - 1][3];
4786 quad_node_list[3] = cell_list[face_cell_no - 1][2];
4788 else if (face_cell_face_no == 2)
4790 quad_node_list[0] = cell_list[face_cell_no - 1][5];
4791 quad_node_list[1] = cell_list[face_cell_no - 1][8];
4792 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4793 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4795 else if (face_cell_face_no == 3)
4797 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4798 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4799 quad_node_list[2] = cell_list[face_cell_no - 1][6];
4800 quad_node_list[3] = cell_list[face_cell_no - 1][5];
4802 else if (face_cell_face_no == 4)
4804 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4805 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4806 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4807 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4809 else if (face_cell_face_no == 5)
4811 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4812 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4813 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4814 quad_node_list[3] = cell_list[face_cell_no - 1][7];
4816 else if (face_cell_face_no == 6)
4818 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4819 quad_node_list[1] = cell_list[face_cell_no - 1][5];
4820 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4821 quad_node_list[3] = cell_list[face_cell_no - 1][4];
4834 return quad_node_list;
4837 template <
int dim,
int spacedim>
4839 Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output)
const
4848 const boost::io::ios_base_all_saver formatting_saver(output);
4852 output <<
"# Abaqus to UCD mesh conversion" << std::endl;
4853 output <<
"# Mesh type: AVS UCD" << std::endl;
4882 output << node_list.size() <<
"\t" << (cell_list.size() + face_list.size())
4883 <<
"\t0\t0\t0" << std::endl;
4886 output.precision(8);
4890 for (
const auto &node : node_list)
4893 output << node[0] <<
"\t";
4896 output.setf(std::ios::scientific, std::ios::floatfield);
4897 for (
unsigned int jj = 1; jj < spacedim + 1; ++jj)
4900 if (
std::abs(node[jj]) > tolerance)
4901 output << static_cast<double>(node[jj]) <<
"\t";
4903 output << 0.0 <<
"\t";
4906 output << 0.0 <<
"\t";
4908 output << std::endl;
4909 output.unsetf(std::ios::floatfield);
4913 for (
unsigned int ii = 0; ii < cell_list.size(); ++ii)
4915 output << ii + 1 <<
"\t" << cell_list[ii][0] <<
"\t"
4916 << (dim == 2 ?
"quad" :
"hex") <<
"\t";
4917 for (
unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
4919 output << cell_list[ii][jj] <<
"\t";
4921 output << std::endl;
4925 for (
unsigned int ii = 0; ii < face_list.size(); ++ii)
4927 output << ii + 1 <<
"\t" << face_list[ii][0] <<
"\t"
4928 << (dim == 2 ?
"line" :
"quad") <<
"\t";
4929 for (
unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
4931 output << face_list[ii][jj] <<
"\t";
4933 output << std::endl;
4940#include "grid_in.inst"
void read_vtk(std::istream &in)
static void skip_empty_lines(std::istream &in)
void read_assimp(const std::string &filename, const unsigned int mesh_index=numbers::invalid_unsigned_int, const bool remove_duplicates=true, const double tol=1e-12, const bool ignore_unsupported_element_types=true)
void read_abaqus(std::istream &in, const bool apply_all_indicators_to_manifolds=false)
static std::string default_suffix(const Format format)
void read_xda(std::istream &in)
void read_comsol_mphtxt(std::istream &in)
static void skip_comment_lines(std::istream &in, const char comment_start)
void attach_triangulation(Triangulation< dim, spacedim > &tria)
void read_msh(std::istream &in)
static Format parse_format(const std::string &format_name)
void read_vtu(std::istream &in)
void read_tecplot(std::istream &in)
ExodusIIData read_exodusii(const std::string &filename, const bool apply_all_indicators_to_manifolds=false)
void read_dbmesh(std::istream &in)
void read_ucd(std::istream &in, const bool apply_all_indicators_to_manifolds=false)
void read(std::istream &in, Format format=Default)
static void debug_output_grid(const std::vector< CellData< dim > > &cells, const std::vector< Point< spacedim > > &vertices, std::ostream &out)
static std::string get_format_names()
void read_unv(std::istream &in)
static void parse_tecplot_header(std::string &header, std::vector< unsigned int > &tecplot2deal, unsigned int &n_vars, unsigned int &n_vertices, unsigned int &n_cells, std::vector< unsigned int > &IJK, bool &structured, bool &blocked)
std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices() const
unsigned int n_vertices() const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const unsigned char face_orientation) const
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int get_dimension() const
ReferenceCell face_reference_cell(const unsigned int face_no) const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_FALLTHROUGH
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcFileNotOpen(std::string arg1)
static ::ExceptionBase & ExcNeedsAssimp()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcNeedsExodusII()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void consistently_order_cells(std::vector< CellData< dim > > &cells)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
std::pair< int, unsigned int > get_integer_at_position(const std::string &name, const unsigned int position)
std::vector< unsigned char > decode_base64(const std::string &base64_input)
std::vector< std::string > break_text_into_lines(const std::string &original_text, const unsigned int width, const char delimiter=' ')
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
bool match_at_string_start(const std::string &name, const std::string &pattern)
std::string decompress(const std::string &compressed_input)
const types::material_id invalid_material_id
const types::boundary_id internal_face_boundary_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< unsigned int > vertices
types::material_id material_id
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
std::vector< std::vector< int > > id_to_sideset_ids
std::vector< CellData< 2 > > boundary_quads
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines