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>
173 text[0] =
"# vtk DataFile Version 3.0";
176 text[3] =
"DATASET UNSTRUCTURED_GRID";
178 for (
unsigned int i = 0; i < 4; ++i)
181 if (i == 2 || i == 3)
183 line.compare(text[i]) == 0,
186 "While reading VTK file, failed to find a header line with text <") +
193 std::vector<Point<spacedim>> vertices;
194 std::vector<CellData<dim>> cells;
203 if (keyword ==
"POINTS")
205 unsigned int n_vertices;
210 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
214 in >> x[0] >> x[1] >> x[2];
216 vertices.emplace_back();
217 for (
unsigned int d = 0; d < spacedim; ++d)
218 vertices.back()[d] = x[d];
225 "While reading VTK file, failed to find POINTS section"));
229 unsigned int n_geometric_objects = 0;
232 if (keyword ==
"CELLS")
235 std::vector<unsigned int> cell_types;
237 std::streampos oldpos = in.tellg();
240 while (in >> keyword)
241 if (keyword ==
"CELL_TYPES")
245 cell_types.resize(n_ints);
247 for (
unsigned int i = 0; i < n_ints; ++i)
256 in >> n_geometric_objects;
261 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
263 unsigned int n_vertices;
267 if (cell_types[count] == 10 || cell_types[count] == 12)
275 cells.emplace_back(n_vertices);
277 for (
unsigned int j = 0; j < n_vertices;
279 in >> cells.back().vertices[j];
283 if (cell_types[count] == 12)
285 std::swap(cells.back().vertices[2],
286 cells.back().vertices[3]);
287 std::swap(cells.back().vertices[6],
288 cells.back().vertices[7]);
291 cells.back().material_id = 0;
294 else if (cell_types[count] == 5 || cell_types[count] == 9)
303 for (
unsigned int j = 0; j < n_vertices;
310 else if (cell_types[count] == 3)
314 for (
unsigned int j = 0; j < n_vertices;
325 "While reading VTK file, unknown cell type encountered"));
330 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
332 unsigned int n_vertices;
336 if (cell_types[count] == 5 || cell_types[count] == 9)
343 cells.emplace_back(n_vertices);
345 for (
unsigned int j = 0; j < n_vertices;
347 in >> cells.back().vertices[j];
351 if (cell_types[count] == 9)
355 std::swap(cells.back().vertices[2],
356 cells.back().vertices[3]);
359 cells.back().material_id = 0;
362 else if (cell_types[count] == 3)
368 for (
unsigned int j = 0; j < n_vertices;
381 "While reading VTK file, unknown cell type encountered"));
386 for (
unsigned int count = 0; count < n_geometric_objects; ++count)
392 cell_types[count] == 3 && type == 2,
394 "While reading VTK file, unknown cell type encountered"));
395 cells.emplace_back(type);
397 for (
unsigned int j = 0; j < type; ++j)
398 in >> cells.back().vertices[j];
400 cells.back().material_id = 0;
406 "While reading VTK file, failed to find CELLS section"));
413 keyword ==
"CELL_TYPES",
415 "While reading VTK file, missing CELL_TYPES section. Found <" +
416 keyword +
"> instead.")));
420 n_ints == n_geometric_objects,
421 ExcMessage(
"The VTK reader found a CELL_DATA statement "
422 "that lists a total of " +
424 " cell data objects, but this needs to "
425 "equal the number of cells (which is " +
427 ") plus the number of quads (" +
429 " in 3d or the number of lines (" +
434 for (
unsigned int i = 0; i < n_ints; ++i)
441 while (in >> keyword)
443 if (keyword ==
"CELL_DATA")
449 n_ids == n_geometric_objects,
451 "The VTK reader found a CELL_DATA statement "
452 "that lists a total of " +
454 " cell data objects, but this needs to "
455 "equal the number of cells (which is " +
457 ") plus the number of quads (" +
459 " in 3d or the number of lines (" +
463 const std::vector<std::string> data_sets{
"MaterialID",
467 for (
unsigned int i = 0; i < data_sets.size(); ++i)
470 if (keyword ==
"SCALARS")
475 std::string field_name;
477 if (std::find(data_sets.begin(),
479 field_name) == data_sets.end())
491 std::getline(in, line);
494 std::min(
static_cast<std::size_t
>(3),
495 line.size() - 1)) ==
"int",
497 "While reading VTK file, material- and manifold IDs can only have type 'int'."));
501 keyword ==
"LOOKUP_TABLE",
503 "While reading VTK file, missing keyword 'LOOKUP_TABLE'."));
507 keyword ==
"default",
509 "While reading VTK file, missing keyword 'default'."));
516 for (
unsigned int i = 0; i < cells.size(); ++i)
520 if (field_name ==
"MaterialID")
521 cells[i].material_id =
523 else if (field_name ==
"ManifoldID")
524 cells[i].manifold_id =
536 if (field_name ==
"MaterialID")
537 boundary_quad.material_id =
539 else if (field_name ==
"ManifoldID")
540 boundary_quad.manifold_id =
549 if (field_name ==
"MaterialID")
550 boundary_line.material_id =
552 else if (field_name ==
"ManifoldID")
553 boundary_line.manifold_id =
565 if (field_name ==
"MaterialID")
566 boundary_line.material_id =
568 else if (field_name ==
"ManifoldID")
569 boundary_line.manifold_id =
579 std::streampos oldpos = in.tellg();
581 if (keyword ==
"SCALARS")
592 else if (keyword ==
"FIELD")
594 unsigned int n_fields;
597 keyword ==
"FieldData",
599 "While reading VTK file, missing keyword FieldData"));
603 for (
unsigned int i = 0; i < n_fields; ++i)
605 std::string section_name;
606 std::string data_type;
607 unsigned int temp, n_ids;
613 n_ids == n_geometric_objects,
615 "The VTK reader found a FIELD statement "
616 "that lists a total of " +
618 " cell data objects, but this needs to equal the number of cells (which is " +
620 ") plus the number of quads (" +
623 " in 3d or the number of lines (" +
630 for (
unsigned int j = 0; j < n_ids; ++j)
633 if (j < cells.size())
636 this->cell_data[section_name] = std::move(temp_data);
647 apply_grid_fixup_functions(vertices, cells, subcelldata);
648 tria->create_triangulation(vertices, cells, subcelldata);
653 "While reading VTK file, failed to find CELLS section"));
656template <
int dim,
int spacedim>
657const std::map<std::string, Vector<double>> &
660 return this->cell_data;
663template <
int dim,
int spacedim>
667 namespace pt = boost::property_tree;
669 pt::read_xml(in, tree);
670 auto section = tree.get_optional<std::string>(
"VTKFile.dealiiData");
674 "While reading a VTU file, failed to find dealiiData section. "
675 "Notice that we can only read grid files in .vtu format that "
676 "were created by the deal.II library, using a call to "
677 "GridOut::write_vtu(), where the flag "
678 "GridOutFlags::Vtu::serialize_triangulation is set to true."));
682 const auto string_archive =
684 std::istringstream in_stream(string_archive);
685 boost::archive::binary_iarchive ia(in_stream);
690template <
int dim,
int spacedim>
694 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
698 skip_comment_lines(in,
'#');
709 "Expected '-1' before and after a section."));
721 std::getline(in, line);
723 boost::algorithm::trim(line);
724 if (line.compare(
"-1") == 0)
735 std::vector<Point<spacedim>> vertices;
754 in >> dummy >> dummy >> dummy;
757 in >> x[0] >> x[1] >> x[2];
759 vertices.emplace_back();
761 for (
unsigned int d = 0; d < spacedim; ++d)
762 vertices.back()[d] = x[d];
776 AssertThrow(tmp == 2412, ExcUnknownSectionType(tmp));
778 std::vector<CellData<dim>> cells;
805 in >> type >> dummy >> dummy >> dummy >> dummy;
807 AssertThrow((type == 11) || (type == 44) || (type == 94) || (type == 115),
808 ExcUnknownElementType(type));
810 if ((((type == 44) || (type == 94)) && (dim == 2)) ||
811 ((type == 115) && (dim == 3)))
813 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
814 cells.emplace_back();
819 .vertices[reference_cell.unv_vertex_to_deal_vertex(v)];
821 cells.back().material_id = 0;
824 cells.back().vertices[v] =
vertex_indices[cells.back().vertices[v]];
826 cell_indices[object_index] = n_cells;
830 else if (((type == 11) && (dim == 2)) ||
831 ((type == 11) && (dim == 3)))
834 in >> dummy >> dummy >> dummy;
839 for (
unsigned int &vertex :
845 for (
unsigned int &vertex :
849 line_indices[object_index] = n_lines;
853 else if (((type == 44) || (type == 94)) && (dim == 3))
864 .vertices[reference_cell.unv_vertex_to_deal_vertex(v)];
868 for (
unsigned int &vertex :
872 quad_indices[object_index] = n_quads;
880 "> when running in dim=" +
899 AssertThrow((tmp == 2467) || (tmp == 2477), ExcUnknownSectionType(tmp));
915 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
926 std::getline(in, line);
929 "The line before the line containing an ID has too "
930 "many entries. This is not a valid UNV file."));
932 std::getline(in, line);
934 std::istringstream id_stream(line);
937 !id_stream.fail() && id_stream.eof(),
939 "The given UNV file contains a boundary or material id set to '" +
941 "', which cannot be parsed as a fixed-width integer, whereas "
942 "deal.II only supports integer boundary and material ids. To fix "
943 "this, ensure that all such ids are given integer values."));
946 id <=
long(std::numeric_limits<types::material_id>::max()),
947 ExcMessage(
"The provided integer id '" + std::to_string(
id) +
948 "' is not convertible to either types::material_id nor "
949 "types::boundary_id."));
951 const unsigned int n_lines =
952 (n_entities % 2 == 0) ? (n_entities / 2) : ((n_entities + 1) / 2);
954 for (
unsigned int line = 0; line < n_lines; ++line)
956 unsigned int n_fragments;
958 if (line == n_lines - 1)
959 n_fragments = (n_entities % 2 == 0) ? (2) : (1);
963 for (
unsigned int no_fragment = 0; no_fragment < n_fragments;
967 in >> dummy >> no >> dummy >> dummy;
969 if (cell_indices.count(no) > 0)
970 cells[cell_indices[no]].material_id = id;
972 if (line_indices.count(no) > 0)
976 if (quad_indices.count(no) > 0)
984 apply_grid_fixup_functions(vertices, cells, subcelldata);
985 tria->create_triangulation(vertices, cells, subcelldata);
990template <
int dim,
int spacedim>
993 const bool apply_all_indicators_to_manifolds)
995 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
999 skip_comment_lines(in,
'#');
1002 unsigned int n_vertices;
1003 unsigned int n_cells;
1006 in >> n_vertices >> n_cells >> dummy
1012 std::vector<Point<spacedim>> vertices(n_vertices);
1018 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1025 in >> vertex_number >> x[0] >> x[1] >> x[2];
1028 for (
unsigned int d = 0; d < spacedim; ++d)
1029 vertices[vertex][d] = x[d];
1037 std::vector<CellData<dim>> cells;
1040 for (
unsigned int cell = 0; cell < n_cells; ++cell)
1049 std::string cell_type;
1053 unsigned int material_id;
1059 if (((cell_type ==
"line") && (dim == 1)) ||
1060 ((cell_type ==
"quad") && (dim == 2)) ||
1061 ((cell_type ==
"hex") && (dim == 3)))
1065 cells.emplace_back();
1070 Assert(material_id <= std::numeric_limits<types::material_id>::max(),
1073 std::numeric_limits<types::material_id>::max()));
1078 if (apply_all_indicators_to_manifolds)
1079 cells.back().manifold_id =
1081 cells.back().material_id = material_id;
1089 cells.back().vertices[i] =
1095 ExcInvalidVertexIndex(cell,
1096 cells.back().vertices[i]));
1101 else if ((cell_type ==
"line") && ((dim == 2) || (dim == 3)))
1109 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
1112 std::numeric_limits<types::boundary_id>::max()));
1121 if (apply_all_indicators_to_manifolds)
1138 for (
unsigned int &vertex :
1146 AssertThrow(
false, ExcInvalidVertexIndex(cell, vertex));
1150 else if ((cell_type ==
"quad") && (dim == 3))
1159 Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
1162 std::numeric_limits<types::boundary_id>::max()));
1171 if (apply_all_indicators_to_manifolds)
1188 for (
unsigned int &vertex :
1196 Assert(
false, ExcInvalidVertexIndex(cell, vertex));
1202 AssertThrow(
false, ExcUnknownIdentifier(cell_type));
1207 apply_grid_fixup_functions(vertices, cells, subcelldata);
1208 tria->create_triangulation(vertices, cells, subcelldata);
1213 template <
int dim,
int spacedim>
1220 read_in_abaqus(std::istream &in);
1222 write_out_avs_ucd(std::ostream &out)
const;
1225 const double tolerance;
1228 get_global_node_numbers(
const int face_cell_no,
1229 const int face_cell_face_no)
const;
1232 std::vector<std::vector<double>> node_list;
1235 std::vector<std::vector<double>> cell_list;
1237 std::vector<std::vector<double>> face_list;
1240 std::map<std::string, std::vector<int>> elsets_list;
1244template <
int dim,
int spacedim>
1247 const bool apply_all_indicators_to_manifolds)
1249 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1254 Assert((spacedim == 2 && dim == spacedim) ||
1255 (spacedim == 3 && (dim == spacedim || dim == spacedim - 1)),
1261 Abaqus_to_UCD<dim, spacedim> abaqus_to_ucd;
1262 abaqus_to_ucd.read_in_abaqus(in);
1264 std::stringstream in_ucd;
1265 abaqus_to_ucd.write_out_avs_ucd(in_ucd);
1273 read_ucd(in_ucd, apply_all_indicators_to_manifolds);
1275 catch (std::exception &exc)
1277 std::cerr <<
"Exception on processing internal UCD data: " << std::endl
1278 << exc.what() << std::endl;
1283 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1284 "More information is provided in an error message printed above. "
1285 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1286 "listed in the documentation?"));
1293 "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
1294 "Are you sure that your ABAQUS mesh file conforms with the requirements "
1295 "listed in the documentation?"));
1300template <
int dim,
int spacedim>
1304 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1310 skip_comment_lines(in,
'#');
1316 AssertThrow(line ==
"MeshVersionFormatted 0", ExcInvalidDBMESHInput(line));
1318 skip_empty_lines(in);
1322 AssertThrow(line ==
"Dimension", ExcInvalidDBMESHInput(line));
1323 unsigned int dimension;
1325 AssertThrow(dimension == dim, ExcDBMESHWrongDimension(dimension));
1326 skip_empty_lines(in);
1339 while (getline(in, line), line.find(
"# END") == std::string::npos)
1341 skip_empty_lines(in);
1346 AssertThrow(line ==
"Vertices", ExcInvalidDBMESHInput(line));
1348 unsigned int n_vertices;
1352 std::vector<Point<spacedim>> vertices(n_vertices);
1353 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1356 for (
unsigned int d = 0; d < dim; ++d)
1357 in >> vertices[vertex][d];
1363 skip_empty_lines(in);
1369 AssertThrow(line ==
"Edges", ExcInvalidDBMESHInput(line));
1371 unsigned int n_edges;
1373 for (
unsigned int edge = 0; edge < n_edges; ++edge)
1376 in >> dummy >> dummy;
1382 skip_empty_lines(in);
1391 AssertThrow(line ==
"CrackedEdges", ExcInvalidDBMESHInput(line));
1394 for (
unsigned int edge = 0; edge < n_edges; ++edge)
1397 in >> dummy >> dummy;
1403 skip_empty_lines(in);
1409 AssertThrow(line ==
"Quadrilaterals", ExcInvalidDBMESHInput(line));
1411 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
1412 {0, 1, 5, 4, 2, 3, 7, 6}};
1413 std::vector<CellData<dim>> cells;
1415 unsigned int n_cells;
1417 for (
unsigned int cell = 0; cell < n_cells; ++cell)
1421 cells.emplace_back();
1425 cells.back().vertices[dim == 3 ? local_vertex_numbering[i] :
1429 (
static_cast<unsigned int>(cells.back().vertices[i]) <=
1431 ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
1433 --cells.back().vertices[i];
1441 skip_empty_lines(in);
1449 while (getline(in, line), ((line.find(
"End") == std::string::npos) && (in)))
1455 apply_grid_fixup_functions(vertices, cells, subcelldata);
1456 tria->create_triangulation(vertices, cells, subcelldata);
1461template <
int dim,
int spacedim>
1465 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1468 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
1472 std::getline(in, line);
1474 unsigned int n_vertices;
1475 unsigned int n_cells;
1479 std::getline(in, line);
1482 std::getline(in, line);
1485 for (
unsigned int i = 0; i < 8; ++i)
1486 std::getline(in, line);
1489 std::vector<CellData<dim>> cells(n_cells);
1500 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1501 in >> cell.vertices[reference_cell.exodusii_vertex_to_deal_vertex(i)];
1505 std::vector<Point<spacedim>> vertices(n_vertices);
1508 for (
unsigned int d = 0; d < spacedim; ++d)
1510 for (
unsigned int d = spacedim; d < 3; ++d)
1519 apply_grid_fixup_functions(vertices, cells, subcelldata);
1520 tria->create_triangulation(vertices, cells, subcelldata);
1525template <
int dim,
int spacedim>
1529 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
1542 std::stringstream whole_file;
1547 std::getline(in, line);
1555 if ((line.size() > 0) && (line.back() ==
'\r'))
1556 line.erase(line.size() - 1);
1561 if (line.find(
'#') != std::string::npos)
1562 line.erase(line.find(
'#'), std::string::npos);
1563 while ((line.size() > 0) && (line.back() ==
' '))
1564 line.erase(line.size() - 1);
1566 if (line.size() > 0)
1567 whole_file <<
'\n' << line;
1586 unsigned int version_major, version_minor;
1587 whole_file >> version_major >> version_minor;
1588 AssertThrow((version_major == 0) && (version_minor == 1),
1589 ExcMessage(
"deal.II can currently only read version 0.1 "
1590 "of the mphtxt file format."));
1595 unsigned int n_tags;
1596 whole_file >> n_tags;
1597 for (
unsigned int i = 0; i < n_tags; ++i)
1600 while (whole_file.peek() ==
'\n')
1602 std::getline(whole_file, dummy);
1608 unsigned int n_types;
1609 whole_file >> n_types;
1610 for (
unsigned int i = 0; i < n_types; ++i)
1613 while (whole_file.peek() ==
'\n')
1615 std::getline(whole_file, dummy);
1635 whole_file >> dummy >> dummy >> dummy;
1639 while (whole_file.peek() ==
'\n')
1641 std::getline(whole_file, s);
1643 ExcMessage(
"Expected '4 Mesh', but got '" + s +
"'."));
1646 unsigned int version;
1647 whole_file >> version;
1651 unsigned int file_space_dim;
1652 whole_file >> file_space_dim;
1656 "The mesh file uses a different number of space dimensions "
1657 "than the triangulation you want to read it into."));
1659 unsigned int n_vertices;
1660 whole_file >> n_vertices;
1662 unsigned int starting_vertex_index;
1663 whole_file >> starting_vertex_index;
1665 std::vector<Point<spacedim>> vertices(n_vertices);
1666 for (
unsigned int v = 0; v < n_vertices; ++v)
1667 whole_file >> vertices[v];
1697 std::vector<CellData<dim>> cells;
1700 unsigned int n_types;
1701 whole_file >> n_types;
1702 for (
unsigned int type = 0; type < n_types; ++type)
1709 whole_file >> dummy;
1713 std::string object_name;
1714 whole_file >> object_name;
1716 static const std::map<std::string, ReferenceCell> name_to_type = {
1726 AssertThrow(name_to_type.find(object_name) != name_to_type.end(),
1727 ExcMessage(
"The input file contains a cell type <" +
1729 "> that the reader does not "
1730 "current support"));
1731 const ReferenceCell object_type = name_to_type.at(object_name);
1733 unsigned int n_vertices_per_element;
1734 whole_file >> n_vertices_per_element;
1736 unsigned int n_elements;
1737 whole_file >> n_elements;
1751 ExcMessage(
"Triangles should not appear in input files "
1759 "Quadrilaterals should not appear in input files "
1766 ExcMessage(
"Tetrahedra should not appear in input files "
1767 "for 1d or 2d meshes."));
1774 "Prisms (wedges) should not appear in input files "
1775 "for 1d or 2d meshes."));
1794 std::vector<unsigned int> vertices_for_this_element(
1795 n_vertices_per_element);
1796 for (
unsigned int e = 0; e < n_elements; ++e)
1799 for (
unsigned int v = 0; v < n_vertices_per_element; ++v)
1801 whole_file >> vertices_for_this_element[v];
1802 vertices_for_this_element[v] -= starting_vertex_index;
1811 cells.emplace_back();
1812 cells.back().vertices = vertices_for_this_element;
1818 vertices_for_this_element;
1826 cells.emplace_back();
1827 cells.back().vertices = vertices_for_this_element;
1833 vertices_for_this_element;
1841 cells.emplace_back();
1842 cells.back().vertices = vertices_for_this_element;
1855 unsigned int n_geom_entity_indices;
1856 whole_file >> n_geom_entity_indices;
1858 (n_geom_entity_indices == n_elements),
1865 if (n_geom_entity_indices != 0)
1867 for (
unsigned int e = 0; e < n_geom_entity_indices; ++e)
1870 unsigned int geometric_entity_index;
1871 whole_file >> geometric_entity_index;
1877 cells[cells.size() - n_elements + e].material_id =
1878 geometric_entity_index;
1883 .boundary_id = geometric_entity_index;
1889 cells[cells.size() - n_elements + e].material_id =
1890 geometric_entity_index;
1895 .boundary_id = geometric_entity_index;
1901 cells[cells.size() - n_elements + e].material_id =
1902 geometric_entity_index;
1918 tria->create_triangulation(vertices, cells, {});
1928 if (line.vertices[1] < line.vertices[0])
1929 std::swap(line.vertices[0], line.vertices[1]);
1934 return std::lexicographical_compare(a.vertices.begin(),
1956 Assert((face.vertices.size() == 3) || (face.vertices.size() == 4),
1958 std::sort(face.vertices.begin(), face.vertices.end());
1963 return std::lexicographical_compare(a.vertices.begin(),
1973 for (
const auto &cell : tria->active_cell_iterators())
1974 for (
const auto &face : cell->face_iterators())
1975 if (face->at_boundary())
1981 std::array<unsigned int, 2> face_vertex_indices = {
1982 {face->vertex_index(0), face->vertex_index(1)}};
1983 if (face_vertex_indices[0] > face_vertex_indices[1])
1984 std::swap(face_vertex_indices[0], face_vertex_indices[1]);
1990 face_vertex_indices,
1992 const std::array<unsigned int, 2>
1993 &face_vertex_indices) ->
bool {
1994 return std::lexicographical_compare(
1997 face_vertex_indices.begin(),
1998 face_vertex_indices.end());
2002 (p->vertices[0] == face_vertex_indices[0]) &&
2003 (p->vertices[1] == face_vertex_indices[1]))
2005 face->set_boundary_id(p->boundary_id);
2013 std::vector<unsigned int> face_vertex_indices(
2014 face->n_vertices());
2015 for (
unsigned int v = 0; v < face->n_vertices(); ++v)
2016 face_vertex_indices[v] = face->vertex_index(v);
2017 std::sort(face_vertex_indices.begin(),
2018 face_vertex_indices.end());
2024 face_vertex_indices,
2026 const std::vector<unsigned int>
2027 &face_vertex_indices) ->
bool {
2028 return std::lexicographical_compare(
2031 face_vertex_indices.begin(),
2032 face_vertex_indices.end());
2036 (p->vertices == face_vertex_indices))
2038 face->set_boundary_id(p->boundary_id);
2043 for (
unsigned int e = 0; e < face->n_lines(); ++e)
2045 const auto edge = face->line(e);
2047 std::array<unsigned int, 2> edge_vertex_indices = {
2048 {edge->vertex_index(0), edge->vertex_index(1)}};
2049 if (edge_vertex_indices[0] > edge_vertex_indices[1])
2050 std::swap(edge_vertex_indices[0],
2051 edge_vertex_indices[1]);
2057 edge_vertex_indices,
2059 const std::array<unsigned int, 2>
2060 &edge_vertex_indices) ->
bool {
2061 return std::lexicographical_compare(
2064 edge_vertex_indices.begin(),
2065 edge_vertex_indices.end());
2069 (p->vertices[0] == edge_vertex_indices[0]) &&
2070 (p->vertices[1] == edge_vertex_indices[1]))
2072 edge->set_boundary_id(p->boundary_id);
2081template <
int dim,
int spacedim>
2085 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
2088 unsigned int n_vertices;
2089 unsigned int n_cells;
2095 std::array<std::map<int, int>, 4> tag_maps;
2098 std::string stripped_file;
2102 while (std::getline(input_stream, line))
2104 if (line ==
"@f$Comments")
2106 while (std::getline(input_stream, line))
2108 if (line ==
"@f$EndComments")
2115 stripped_file += line +
'\n';
2119 std::istringstream in(stripped_file);
2124 unsigned int gmsh_file_format = 0;
2125 if (line ==
"@f$NOD")
2126 gmsh_file_format = 10;
2127 else if (line ==
"@f$MeshFormat")
2128 gmsh_file_format = 20;
2134 if (gmsh_file_format == 20)
2137 unsigned int file_type, data_size;
2139 in >> version >> file_type >> data_size;
2142 gmsh_file_format =
static_cast<unsigned int>(version * 10);
2150 AssertThrow(line ==
"@f$EndMeshFormat", ExcInvalidGMSHInput(line));
2154 if (line ==
"@f$PhysicalNames")
2160 while (line !=
"@f$EndPhysicalNames");
2165 if (line ==
"@f$Entities")
2167 unsigned long n_points, n_curves, n_surfaces, n_volumes;
2169 in >> n_points >> n_curves >> n_surfaces >> n_volumes;
2170 for (
unsigned int i = 0; i < n_points; ++i)
2174 unsigned int n_physicals;
2175 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2179 if (gmsh_file_format > 40)
2181 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2183 box_max_x = box_min_x;
2184 box_max_y = box_min_y;
2185 box_max_z = box_min_z;
2189 in >> tag >> box_min_x >> box_min_y >> box_min_z >>
2190 box_max_x >> box_max_y >> box_max_z >> n_physicals;
2195 ExcMessage(
"More than one tag is not supported!"));
2197 int physical_tag = 0;
2198 for (
unsigned int j = 0; j < n_physicals; ++j)
2200 tag_maps[0][tag] = physical_tag;
2202 for (
unsigned int i = 0; i < n_curves; ++i)
2206 unsigned int n_physicals;
2207 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2211 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2212 box_max_y >> box_max_z >> n_physicals;
2216 ExcMessage(
"More than one tag is not supported!"));
2218 int physical_tag = 0;
2219 for (
unsigned int j = 0; j < n_physicals; ++j)
2221 tag_maps[1][tag] = physical_tag;
2226 for (
unsigned int j = 0; j < n_points; ++j)
2230 for (
unsigned int i = 0; i < n_surfaces; ++i)
2234 unsigned int n_physicals;
2235 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2239 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2240 box_max_y >> box_max_z >> n_physicals;
2244 ExcMessage(
"More than one tag is not supported!"));
2246 int physical_tag = 0;
2247 for (
unsigned int j = 0; j < n_physicals; ++j)
2249 tag_maps[2][tag] = physical_tag;
2254 for (
unsigned int j = 0; j < n_curves; ++j)
2257 for (
unsigned int i = 0; i < n_volumes; ++i)
2261 unsigned int n_physicals;
2262 double box_min_x, box_min_y, box_min_z, box_max_x, box_max_y,
2266 in >> tag >> box_min_x >> box_min_y >> box_min_z >> box_max_x >>
2267 box_max_y >> box_max_z >> n_physicals;
2271 ExcMessage(
"More than one tag is not supported!"));
2273 int physical_tag = 0;
2274 for (
unsigned int j = 0; j < n_physicals; ++j)
2276 tag_maps[3][tag] = physical_tag;
2281 for (
unsigned int j = 0; j < n_surfaces; ++j)
2285 AssertThrow(line ==
"@f$EndEntities", ExcInvalidGMSHInput(line));
2290 if (line ==
"@f$PartitionedEntities")
2296 while (line !=
"@f$EndPartitionedEntities");
2303 AssertThrow(line ==
"@f$Nodes", ExcInvalidGMSHInput(line));
2307 int n_entity_blocks = 1;
2308 if (gmsh_file_format > 40)
2312 in >> n_entity_blocks >> n_vertices >> min_node_tag >> max_node_tag;
2314 else if (gmsh_file_format == 40)
2316 in >> n_entity_blocks >> n_vertices;
2320 std::vector<Point<spacedim>> vertices(n_vertices);
2327 unsigned int global_vertex = 0;
2328 for (
int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2331 unsigned long numNodes;
2333 if (gmsh_file_format < 40)
2335 numNodes = n_vertices;
2342 int tagEntity, dimEntity;
2343 in >> tagEntity >> dimEntity >> parametric >> numNodes;
2346 std::vector<int> vertex_numbers;
2348 if (gmsh_file_format > 40)
2349 for (
unsigned long vertex_per_entity = 0;
2350 vertex_per_entity < numNodes;
2351 ++vertex_per_entity)
2353 in >> vertex_number;
2354 vertex_numbers.push_back(vertex_number);
2357 for (
unsigned long vertex_per_entity = 0; vertex_per_entity < numNodes;
2358 ++vertex_per_entity, ++global_vertex)
2364 if (gmsh_file_format > 40)
2366 vertex_number = vertex_numbers[vertex_per_entity];
2367 in >> x[0] >> x[1] >> x[2];
2370 in >> vertex_number >> x[0] >> x[1] >> x[2];
2372 for (
unsigned int d = 0; d < spacedim; ++d)
2373 vertices[global_vertex][d] = x[d];
2378 if (parametric != 0)
2393 static const std::string end_nodes_marker[] = {
"@f$ENDNOD",
"@f$EndNodes"};
2394 AssertThrow(line == end_nodes_marker[gmsh_file_format == 10 ? 0 : 1],
2395 ExcInvalidGMSHInput(line));
2399 static const std::string begin_elements_marker[] = {
"@f$ELM",
"@f$Elements"};
2400 AssertThrow(line == begin_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2401 ExcInvalidGMSHInput(line));
2404 if (gmsh_file_format > 40)
2408 in >> n_entity_blocks >> n_cells >> min_node_tag >> max_node_tag;
2410 else if (gmsh_file_format == 40)
2412 in >> n_entity_blocks >> n_cells;
2416 n_entity_blocks = 1;
2423 std::vector<CellData<dim>> cells;
2425 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2431 std::map<unsigned int, unsigned int> vertex_counts;
2434 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
2435 {0, 1, 5, 4, 2, 3, 7, 6}};
2436 unsigned int global_cell = 0;
2437 for (
int entity_block = 0; entity_block < n_entity_blocks; ++entity_block)
2439 unsigned int material_id;
2440 unsigned long numElements;
2443 if (gmsh_file_format < 40)
2447 numElements = n_cells;
2449 else if (gmsh_file_format == 40)
2451 int tagEntity, dimEntity;
2452 in >> tagEntity >> dimEntity >> cell_type >> numElements;
2453 material_id = tag_maps[dimEntity][tagEntity];
2458 int tagEntity, dimEntity;
2459 in >> dimEntity >> tagEntity >> cell_type >> numElements;
2460 material_id = tag_maps[dimEntity][tagEntity];
2463 for (
unsigned int cell_per_entity = 0; cell_per_entity < numElements;
2464 ++cell_per_entity, ++global_cell)
2473 unsigned int nod_num;
2494 unsigned int elm_number = 0;
2495 if (gmsh_file_format < 40)
2501 if (gmsh_file_format < 20)
2507 else if (gmsh_file_format < 40)
2512 unsigned int n_tags;
2519 for (
unsigned int i = 1; i < n_tags; ++i)
2524 else if (cell_type == 2)
2526 else if (cell_type == 3)
2528 else if (cell_type == 4)
2530 else if (cell_type == 5)
2541 else if (cell_type == 2)
2543 else if (cell_type == 3)
2545 else if (cell_type == 4)
2547 else if (cell_type == 5)
2573 if (((cell_type == 1) && (dim == 1)) ||
2574 ((cell_type == 2) && (dim == 2)) ||
2575 ((cell_type == 3) && (dim == 2)) ||
2576 ((cell_type == 4) && (dim == 3)) ||
2577 ((cell_type == 5) && (dim == 3)))
2580 unsigned int vertices_per_cell = 0;
2582 vertices_per_cell = 2;
2583 else if (cell_type == 2)
2584 vertices_per_cell = 3;
2585 else if (cell_type == 3)
2586 vertices_per_cell = 4;
2587 else if (cell_type == 4)
2588 vertices_per_cell = 4;
2589 else if (cell_type == 5)
2590 vertices_per_cell = 8;
2594 "Number of nodes does not coincide with the "
2595 "number required for this object"));
2598 cells.emplace_back();
2600 cell.
vertices.resize(vertices_per_cell);
2601 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2604 if (vertices_per_cell ==
2607 local_vertex_numbering[i] :
2615 std::numeric_limits<types::material_id>::max(),
2619 std::numeric_limits<types::material_id>::max()));
2627 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2631 ExcInvalidVertexIndexGmsh(cell_per_entity,
2637 vertex_counts[vertex] += 1u;
2641 else if ((cell_type == 1) &&
2642 ((dim == 2) || (dim == 3)))
2651 std::numeric_limits<types::boundary_id>::max(),
2655 std::numeric_limits<types::boundary_id>::max()));
2666 for (
unsigned int &vertex :
2675 ExcInvalidVertexIndex(cell_per_entity,
2680 else if ((cell_type == 2 || cell_type == 3) &&
2684 unsigned int vertices_per_cell = 0;
2687 vertices_per_cell = 3;
2688 else if (cell_type == 3)
2689 vertices_per_cell = 4;
2697 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
2702 std::numeric_limits<types::boundary_id>::max(),
2706 std::numeric_limits<types::boundary_id>::max()));
2717 for (
unsigned int &vertex :
2726 ExcInvalidVertexIndex(cell_per_entity, vertex));
2730 else if (cell_type == 15)
2733 unsigned int node_index = 0;
2734 if (gmsh_file_format < 20)
2754 AssertThrow(
false, ExcGmshUnsupportedGeometry(cell_type));
2762 static const std::string end_elements_marker[] = {
"@f$ENDELM",
"@f$EndElements"};
2763 AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1],
2764 ExcInvalidGMSHInput(line));
2777 for (
const auto &pair : vertex_counts)
2778 if (pair.second == 1u)
2779 boundary_id_pairs.emplace_back(vertices[pair.first],
2780 boundary_ids_1d[pair.first]);
2782 apply_grid_fixup_functions(vertices, cells, subcelldata);
2783 tria->create_triangulation(vertices, cells, subcelldata);
2788 assign_1d_boundary_ids(boundary_id_pairs, *tria);
2793#ifdef DEAL_II_GMSH_WITH_API
2794template <
int dim,
int spacedim>
2798 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
2800 const std::map<int, std::uint8_t> gmsh_to_dealii_type = {
2801 {{15, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {7, 5}, {6, 6}, {5, 7}}};
2804 const std::array<std::vector<unsigned int>, 8> gmsh_to_dealii = {
2811 {{0, 1, 2, 3, 4, 5}},
2812 {{0, 1, 3, 2, 4, 5, 7, 6}}}};
2814 std::vector<Point<spacedim>> vertices;
2815 std::vector<CellData<dim>> cells;
2817 std::map<unsigned int, types::boundary_id> boundary_ids_1d;
2823 std::map<unsigned int, unsigned int> vertex_counts;
2826 gmsh::option::setNumber(
"General.Verbosity", 0);
2830 ExcMessage(
"You are trying to read a gmsh file with dimension " +
2831 std::to_string(gmsh::model::getDimension()) +
2832 " into a grid of dimension " + std::to_string(dim)));
2837 gmsh::model::mesh::removeDuplicateNodes();
2838 gmsh::model::mesh::renumberNodes();
2839 std::vector<std::size_t> node_tags;
2840 std::vector<double> coord;
2841 std::vector<double> parametricCoord;
2842 gmsh::model::mesh::getNodes(
2843 node_tags, coord, parametricCoord, -1, -1,
false,
false);
2844 vertices.resize(node_tags.size());
2845 for (
unsigned int i = 0; i < node_tags.size(); ++i)
2849 for (
unsigned int d = 0; d < spacedim; ++d)
2850 vertices[i][d] = coord[i * 3 + d];
2854 for (
unsigned int d = spacedim; d < 3; ++d)
2857 "The grid you are reading contains nodes that are "
2858 "nonzero in the coordinate with index " +
2860 ", but you are trying to save "
2861 "it on a grid embedded in a " +
2862 std::to_string(spacedim) +
" dimensional space."));
2869 std::vector<std::pair<int, int>> entities;
2870 gmsh::model::getEntities(entities);
2872 for (
const auto &e : entities)
2875 const int &entity_dim = e.first;
2876 const int &entity_tag = e.second;
2882 std::vector<int> physical_tags;
2883 gmsh::model::getPhysicalGroupsForEntity(entity_dim,
2888 if (physical_tags.size())
2889 for (
auto physical_tag : physical_tags)
2892 gmsh::model::getPhysicalName(entity_dim, physical_tag, name);
2899 std::map<std::string, int> id_names;
2901 bool found_unrecognized_tag =
false;
2902 bool found_boundary_id =
false;
2905 for (
const auto &it : id_names)
2907 const auto &name = it.first;
2908 const auto &
id = it.second;
2909 if (entity_dim == dim && name ==
"MaterialID")
2912 found_boundary_id =
true;
2914 else if (entity_dim < dim && name ==
"BoundaryID")
2917 found_boundary_id =
true;
2919 else if (name ==
"ManifoldID")
2925 found_unrecognized_tag =
true;
2930 if (found_unrecognized_tag && !found_boundary_id)
2931 boundary_id = physical_tag;
2939 boundary_id = physical_tag;
2945 std::vector<int> element_types;
2946 std::vector<std::vector<std::size_t>> element_ids, element_nodes;
2947 gmsh::model::mesh::getElements(
2948 element_types, element_ids, element_nodes, entity_dim, entity_tag);
2950 for (
unsigned int i = 0; i < element_types.size(); ++i)
2952 const auto &type = gmsh_to_dealii_type.at(element_types[i]);
2953 const auto n_vertices = gmsh_to_dealii[type].size();
2954 const auto &elements = element_ids[i];
2955 const auto &nodes = element_nodes[i];
2956 for (
unsigned int j = 0; j < elements.size(); ++j)
2958 if (entity_dim == dim)
2960 cells.emplace_back(n_vertices);
2961 auto &cell = cells.back();
2962 for (
unsigned int v = 0; v < n_vertices; ++v)
2965 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2967 vertex_counts[cell.vertices[v]] += 1u;
2969 cell.manifold_id = manifold_id;
2970 cell.material_id = boundary_id;
2972 else if (entity_dim == 2)
2976 for (
unsigned int v = 0; v < n_vertices; ++v)
2978 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2980 face.manifold_id = manifold_id;
2981 face.boundary_id = boundary_id;
2983 else if (entity_dim == 1)
2987 for (
unsigned int v = 0; v < n_vertices; ++v)
2989 nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1;
2991 line.manifold_id = manifold_id;
2992 line.boundary_id = boundary_id;
2994 else if (entity_dim == 0)
2998 for (
unsigned int j = 0; j < elements.size(); ++j)
2999 boundary_ids_1d[nodes[j] - 1] = boundary_id;
3010 for (
const auto &pair : vertex_counts)
3011 if (pair.second == 1u)
3012 boundary_id_pairs.emplace_back(vertices[pair.first],
3013 boundary_ids_1d[pair.first]);
3015 apply_grid_fixup_functions(vertices, cells, subcelldata);
3016 tria->create_triangulation(vertices, cells, subcelldata);
3021 assign_1d_boundary_ids(boundary_id_pairs, *tria);
3030template <
int dim,
int spacedim>
3033 std::string &header,
3034 std::vector<unsigned int> &tecplot2deal,
3035 unsigned int &n_vars,
3036 unsigned int &n_vertices,
3037 unsigned int &n_cells,
3038 std::vector<unsigned int> &IJK,
3063 std::transform(header.begin(),
3066 static_cast<int (*)(
int)
>(std::toupper));
3070 std::replace(header.begin(), header.end(),
'\t',
' ');
3071 std::replace(header.begin(), header.end(),
',',
' ');
3072 std::replace(header.begin(), header.end(),
'\n',
' ');
3076 std::string::size_type pos = header.find(
'=');
3078 while (pos !=
static_cast<std::string::size_type
>(std::string::npos))
3079 if (header[pos + 1] ==
' ')
3080 header.erase(pos + 1, 1);
3081 else if (header[pos - 1] ==
' ')
3083 header.erase(pos - 1, 1);
3087 pos = header.find(
'=', ++pos);
3090 std::vector<std::string> entries =
3094 for (
unsigned int i = 0; i < entries.size(); ++i)
3103 tecplot2deal[0] = 0;
3106 while (entries[i][0] ==
'"')
3108 if (entries[i] ==
"\"X\"")
3109 tecplot2deal[0] = n_vars;
3110 else if (entries[i] ==
"\"Y\"")
3116 tecplot2deal[1] = n_vars;
3118 else if (entries[i] ==
"\"Z\"")
3124 tecplot2deal[2] = n_vars;
3136 "Tecplot file must contain at least one variable for each dimension"));
3137 for (
unsigned int d = 1; d < dim; ++d)
3139 tecplot2deal[d] > 0,
3141 "Tecplot file must contain at least one variable for each dimension."));
3146 "ZONETYPE=FELINESEG") &&
3150 "ZONETYPE=FEQUADRILATERAL") &&
3154 "ZONETYPE=FEBRICK") &&
3162 "The tecplot file contains an unsupported ZONETYPE."));
3165 "DATAPACKING=POINT"))
3168 "DATAPACKING=BLOCK"))
3191 "ET=QUADRILATERAL") &&
3203 "The tecplot file contains an unsupported ElementType."));
3211 dim > 1 || IJK[1] == 1,
3213 "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
3219 dim > 2 || IJK[2] == 1,
3221 "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
3236 for (
unsigned int d = 0; d < dim; ++d)
3241 "Tecplot file does not contain a complete and consistent set of parameters"));
3242 n_vertices *= IJK[d];
3243 n_cells *= (IJK[d] - 1);
3251 "Tecplot file does not contain a complete and consistent set of parameters"));
3257 n_cells = *std::max_element(IJK.begin(), IJK.end());
3261 "Tecplot file does not contain a complete and consistent set of parameters"));
3271 const unsigned int dim = 2;
3272 const unsigned int spacedim = 2;
3273 Assert(tria !=
nullptr, ExcNoTriangulationSelected());
3277 skip_comment_lines(in,
'#');
3280 std::string line, header;
3287 std::string letters =
"abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";
3290 while (line.find_first_of(letters) != std::string::npos)
3292 header +=
" " + line;
3299 std::vector<unsigned int> tecplot2deal(dim);
3300 std::vector<unsigned int> IJK(dim);
3301 unsigned int n_vars, n_vertices, n_cells;
3302 bool structured, blocked;
3304 parse_tecplot_header(header,
3319 std::vector<Point<spacedim>> vertices(n_vertices + 1);
3322 std::vector<CellData<dim>> cells(n_cells);
3338 unsigned int next_index = 0;
3342 if (tecplot2deal[0] == 0)
3346 std::vector<std::string> first_var =
3349 for (
unsigned int i = 1; i < first_var.size() + 1; ++i)
3350 vertices[i][0] = std::strtod(first_var[i - 1].c_str(), &endptr);
3355 for (
unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
3356 in >> vertices[j][next_index];
3363 for (
unsigned int i = 1; i < n_vars; ++i)
3372 if (next_index == dim && structured)
3375 if ((next_index < dim) && (i == tecplot2deal[next_index]))
3378 for (
unsigned int j = 1; j < n_vertices + 1; ++j)
3379 in >> vertices[j][next_index];
3386 for (
unsigned int j = 1; j < n_vertices + 1; ++j)
3398 std::vector<double> vars(n_vars);
3403 std::vector<std::string> first_vertex =
3406 for (
unsigned int d = 0; d < dim; ++d)
3408 std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);
3412 for (
unsigned int v = 2; v < n_vertices + 1; ++v)
3414 for (
unsigned int i = 0; i < n_vars; ++i)
3420 for (
unsigned int i = 0; i < dim; ++i)
3421 vertices[v][i] = vars[tecplot2deal[i]];
3429 unsigned int I = IJK[0], J = IJK[1];
3431 unsigned int cell = 0;
3433 for (
unsigned int j = 0; j < J - 1; ++j)
3434 for (
unsigned int i = 1; i < I; ++i)
3436 cells[cell].vertices[0] = i + j * I;
3437 cells[cell].vertices[1] = i + 1 + j * I;
3438 cells[cell].vertices[2] = i + (j + 1) * I;
3439 cells[cell].vertices[3] = i + 1 + (j + 1) * I;
3443 std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
3445 for (
unsigned int i = 1; i < I + 1; ++i)
3447 boundary_vertices[k] = i;
3449 boundary_vertices[k] = i + (J - 1) * I;
3452 for (
unsigned int j = 1; j < J - 1; ++j)
3454 boundary_vertices[k] = 1 + j * I;
3456 boundary_vertices[k] = I + j * I;
3475 for (
unsigned int i = 0; i < n_cells; ++i)
3492 apply_grid_fixup_functions(vertices, cells, subcelldata);
3493 tria->create_triangulation(vertices, cells, subcelldata);
3498template <
int dim,
int spacedim>
3507template <
int dim,
int spacedim>
3510 const unsigned int mesh_index,
3511 const bool remove_duplicates,
3513 const bool ignore_unsupported_types)
3515#ifdef DEAL_II_WITH_ASSIMP
3520 Assimp::Importer importer;
3523 const aiScene *scene =
3524 importer.ReadFile(filename.c_str(),
3525 aiProcess_RemoveComponent |
3526 aiProcess_JoinIdenticalVertices |
3527 aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
3528 aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);
3534 ExcMessage(
"Input file contains no meshes."));
3537 (mesh_index < scene->mNumMeshes),
3540 unsigned int start_mesh =
3542 unsigned int end_mesh =
3547 std::vector<Point<spacedim>> vertices;
3548 std::vector<CellData<dim>> cells;
3552 unsigned int v_offset = 0;
3553 unsigned int c_offset = 0;
3555 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
3556 {0, 1, 5, 4, 2, 3, 7, 6}};
3558 for (
unsigned int m = start_mesh; m < end_mesh; ++m)
3560 const aiMesh *mesh = scene->mMeshes[m];
3564 if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
3567 ExcMessage(
"Incompatible mesh " + std::to_string(m) +
3568 "/" + std::to_string(scene->mNumMeshes)));
3571 else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
3574 ExcMessage(
"Incompatible mesh " + std::to_string(m) +
3575 "/" + std::to_string(scene->mNumMeshes)));
3579 const unsigned int n_vertices = mesh->mNumVertices;
3580 const aiVector3D *mVertices = mesh->mVertices;
3583 const unsigned int n_faces = mesh->mNumFaces;
3584 const aiFace *mFaces = mesh->mFaces;
3586 vertices.resize(v_offset + n_vertices);
3587 cells.resize(c_offset + n_faces);
3589 for (
unsigned int i = 0; i < n_vertices; ++i)
3590 for (
unsigned int d = 0; d < spacedim; ++d)
3591 vertices[i + v_offset][d] = mVertices[i][d];
3593 unsigned int valid_cell = c_offset;
3594 for (
unsigned int i = 0; i < n_faces; ++i)
3601 .vertices[dim == 3 ? local_vertex_numbering[f] :
3603 mFaces[i].mIndices[f] + v_offset;
3605 cells[valid_cell].material_id = m;
3611 ExcMessage(
"Face " + std::to_string(i) +
" of mesh " +
3612 std::to_string(m) +
" has " +
3613 std::to_string(mFaces[i].mNumIndices) +
3614 " vertices. We expected only " +
3619 cells.resize(valid_cell);
3624 v_offset += n_vertices;
3625 c_offset = valid_cell;
3632 if (remove_duplicates)
3637 unsigned int n_verts = 0;
3638 while (n_verts != vertices.size())
3640 n_verts = vertices.size();
3641 std::vector<unsigned int> considered_vertices;
3643 vertices, cells, subcelldata, considered_vertices, tol);
3647 apply_grid_fixup_functions(vertices, cells, subcelldata);
3648 tria->create_triangulation(vertices, cells, subcelldata);
3653 (void)remove_duplicates;
3655 (void)ignore_unsupported_types;
3660#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3667 exodusii_name_to_type(
const std::string &type_name,
3668 const int n_nodes_per_element)
3674 std::string type_name_2 = type_name;
3675 std::transform(type_name_2.begin(),
3677 type_name_2.begin(),
3678 static_cast<int (*)(
int)
>(std::toupper));
3679 const std::string
numbers =
"0123456789";
3680 type_name_2.erase(std::find_first_of(type_name_2.begin(),
3687 if (type_name_2 ==
"BAR" || type_name_2 ==
"BEAM" ||
3688 type_name_2 ==
"EDGE" || type_name_2 ==
"TRUSS")
3690 else if (type_name_2 ==
"TRI" || type_name_2 ==
"TRIANGLE")
3692 else if (type_name_2 ==
"QUAD" || type_name_2 ==
"QUADRILATERAL")
3694 else if (type_name_2 ==
"SHELL")
3696 if (n_nodes_per_element == 3)
3701 else if (type_name_2 ==
"TET" || type_name_2 ==
"TETRA" ||
3702 type_name_2 ==
"TETRAHEDRON")
3704 else if (type_name_2 ==
"PYRA" || type_name_2 ==
"PYRAMID")
3706 else if (type_name_2 ==
"WEDGE")
3708 else if (type_name_2 ==
"HEX" || type_name_2 ==
"HEXAHEDRON")
3718 template <
int dim,
int spacedim = dim>
3719 std::pair<SubCellData, std::vector<std::vector<int>>>
3720 read_exodusii_sidesets(
const int ex_id,
3721 const int n_side_sets,
3723 const bool apply_all_indicators_to_manifolds)
3726 std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
3728 b_or_m_id_to_sideset_ids.emplace_back();
3734 if (dim == spacedim && n_side_sets > 0)
3736 std::vector<int> side_set_ids(n_side_sets);
3737 int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
3745 std::map<std::size_t, std::vector<int>> face_side_sets;
3746 for (
const int side_set_id : side_set_ids)
3749 int n_distribution_factors = -1;
3751 ierr = ex_get_set_param(ex_id,
3755 &n_distribution_factors);
3759 std::vector<int> elements(n_sides);
3760 std::vector<int> faces(n_sides);
3761 ierr = ex_get_set(ex_id,
3774 for (
int side_n = 0; side_n < n_sides; ++side_n)
3776 const long element_n = elements[side_n] - 1;
3777 const long face_n = faces[side_n] - 1;
3778 const std::size_t face_id =
3779 element_n * ReferenceCells::max_n_faces<dim>() + face_n;
3780 face_side_sets[face_id].push_back(side_set_id);
3786 std::vector<std::pair<std::size_t, std::vector<int>>>
3787 face_id_to_side_sets;
3788 face_id_to_side_sets.reserve(face_side_sets.size());
3789 for (
auto &pair : face_side_sets)
3792 face_id_to_side_sets.emplace_back(std::move(pair));
3796 std::sort(face_id_to_side_sets.begin(),
3797 face_id_to_side_sets.end(),
3798 [](
const auto &a,
const auto &b) {
3799 return std::lexicographical_compare(a.second.begin(),
3810 for (
const auto &pair : face_id_to_side_sets)
3812 const std::size_t face_id = pair.first;
3813 const std::vector<int> &face_sideset_ids = pair.second;
3814 if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
3819 ++current_b_or_m_id;
3820 b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
3821 Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
3825 const unsigned int local_face_n =
3826 face_id % ReferenceCells::max_n_faces<dim>();
3828 cells[face_id / ReferenceCells::max_n_faces<dim>()];
3831 const unsigned int deal_face_n =
3842 if (apply_all_indicators_to_manifolds)
3843 boundary_line.manifold_id = current_b_or_m_id;
3845 boundary_line.boundary_id = current_b_or_m_id;
3846 for (
unsigned int j = 0; j < face_reference_cell.
n_vertices();
3848 boundary_line.vertices[j] =
3857 if (apply_all_indicators_to_manifolds)
3858 boundary_quad.manifold_id = current_b_or_m_id;
3860 boundary_quad.boundary_id = current_b_or_m_id;
3861 for (
unsigned int j = 0; j < face_reference_cell.
n_vertices();
3863 boundary_quad.vertices[j] =
3872 return std::make_pair(std::move(subcelldata),
3873 std::move(b_or_m_id_to_sideset_ids));
3878template <
int dim,
int spacedim>
3881 const std::string &filename,
3882 const bool apply_all_indicators_to_manifolds)
3884#ifdef DEAL_II_TRILINOS_WITH_SEACAS
3886 int component_word_size =
sizeof(double);
3888 int floating_point_word_size = 0;
3889 float ex_version = 0.0;
3891 const int ex_id = ex_open(filename.c_str(),
3893 &component_word_size,
3894 &floating_point_word_size,
3897 ExcMessage(
"ExodusII failed to open the specified input file."));
3900 std::vector<char> string_temp(MAX_LINE_LENGTH + 1,
'\0');
3901 int mesh_dimension = 0;
3904 int n_element_blocks = 0;
3905 int n_node_sets = 0;
3906 int n_side_sets = 0;
3908 int ierr = ex_get_init(ex_id,
3927 std::vector<Point<spacedim>> vertices;
3928 vertices.reserve(n_nodes);
3930 std::vector<double> xs(n_nodes);
3931 std::vector<double> ys(n_nodes);
3932 std::vector<double> zs(n_nodes);
3934 ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
3937 for (
int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
3942 vertices.emplace_back(xs[vertex_n]);
3945 vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
3948 vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
3956 std::vector<int> element_block_ids(n_element_blocks);
3957 ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
3960 std::vector<CellData<dim>> cells;
3961 cells.reserve(n_elements);
3965 for (
const int element_block_id : element_block_ids)
3967 std::fill(string_temp.begin(), string_temp.end(),
'\0');
3968 int n_block_elements = 0;
3969 int n_nodes_per_element = 0;
3970 int n_edges_per_element = 0;
3971 int n_faces_per_element = 0;
3972 int n_attributes_per_element = 0;
3975 ierr = ex_get_block(ex_id,
3980 &n_nodes_per_element,
3981 &n_edges_per_element,
3982 &n_faces_per_element,
3983 &n_attributes_per_element);
3986 exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
3989 "The ExodusII block " + std::to_string(element_block_id) +
3990 " with element type " + std::string(string_temp.data()) +
3992 ", which does not match the topological mesh dimension " +
3993 std::to_string(dim) +
"."));
4000 std::vector<int> connection(n_nodes_per_element * n_block_elements);
4001 ierr = ex_get_conn(ex_id,
4009 for (
unsigned int elem_n = 0; elem_n < connection.size();
4010 elem_n += n_nodes_per_element)
4016 connection[elem_n + i] - 1;
4019 cells.push_back(std::move(cell));
4024 auto pair = read_exodusii_sidesets<dim, spacedim>(
4025 ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
4026 ierr = ex_close(ex_id);
4029 apply_grid_fixup_functions(vertices, cells, pair.first);
4030 tria->create_triangulation(vertices, cells, pair.first);
4036 (void)apply_all_indicators_to_manifolds;
4043template <
int dim,
int spacedim>
4057 if (std::find_if(line.begin(), line.end(), [](
const char c) {
4062 for (
int i = line.size() - 1; i >= 0; --i)
4063 in.putback(line[i]);
4073template <
int dim,
int spacedim>
4076 const char comment_start)
4081 while (in.get(c) && c == comment_start)
4084 while (in.get() !=
'\n')
4094 skip_empty_lines(in);
4099template <
int dim,
int spacedim>
4114 const std::vector<
Point<2>> &vertices,
4117 double min_x = vertices[cells[0].vertices[0]][0],
4118 max_x = vertices[cells[0].vertices[0]][0],
4119 min_y = vertices[cells[0].vertices[0]][1],
4120 max_y = vertices[cells[0].vertices[0]][1];
4122 for (
unsigned int i = 0; i < cells.size(); ++i)
4124 for (
const auto vertex : cells[i].vertices)
4126 const Point<2> &p = vertices[vertex];
4138 out <<
"# cell " << i << std::endl;
4140 for (
const auto vertex : cells[i].vertices)
4141 center += vertices[vertex];
4144 out <<
"set label \"" << i <<
"\" at " << center[0] <<
',' << center[1]
4145 <<
" center" << std::endl;
4148 for (
unsigned int f = 0; f < 2; ++f)
4149 out <<
"set arrow from " << vertices[cells[i].vertices[f]][0] <<
','
4150 << vertices[cells[i].vertices[f]][1] <<
" to "
4151 << vertices[cells[i].vertices[(f + 1) % 4]][0] <<
','
4152 << vertices[cells[i].vertices[(f + 1) % 4]][1] << std::endl;
4154 for (
unsigned int f = 2; f < 4; ++f)
4155 out <<
"set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]][0]
4156 <<
',' << vertices[cells[i].vertices[(f + 1) % 4]][1] <<
" to "
4157 << vertices[cells[i].vertices[f]][0] <<
','
4158 << vertices[cells[i].vertices[f]][1] << std::endl;
4164 <<
"set nokey" << std::endl
4165 <<
"pl [" << min_x <<
':' << max_x <<
"][" << min_y <<
':' << max_y
4166 <<
"] " << min_y << std::endl
4167 <<
"pause -1" << std::endl;
4175 const std::vector<
Point<3>> &vertices,
4178 for (
const auto &cell : cells)
4181 out << vertices[cell.
vertices[0]] << std::endl
4182 << vertices[cell.
vertices[1]] << std::endl
4186 out << vertices[cell.
vertices[1]] << std::endl
4187 << vertices[cell.
vertices[2]] << std::endl
4191 out << vertices[cell.
vertices[3]] << std::endl
4192 << vertices[cell.
vertices[2]] << std::endl
4196 out << vertices[cell.
vertices[0]] << std::endl
4197 << vertices[cell.
vertices[3]] << std::endl
4201 out << vertices[cell.
vertices[4]] << std::endl
4202 << vertices[cell.
vertices[5]] << std::endl
4206 out << vertices[cell.
vertices[5]] << std::endl
4207 << vertices[cell.
vertices[6]] << std::endl
4211 out << vertices[cell.
vertices[7]] << std::endl
4212 << vertices[cell.
vertices[6]] << std::endl
4216 out << vertices[cell.
vertices[4]] << std::endl
4217 << vertices[cell.
vertices[7]] << std::endl
4221 out << vertices[cell.
vertices[0]] << std::endl
4222 << vertices[cell.
vertices[4]] << std::endl
4226 out << vertices[cell.
vertices[1]] << std::endl
4227 << vertices[cell.
vertices[5]] << std::endl
4231 out << vertices[cell.
vertices[2]] << std::endl
4232 << vertices[cell.
vertices[6]] << std::endl
4236 out << vertices[cell.
vertices[3]] << std::endl
4237 << vertices[cell.
vertices[7]] << std::endl
4245template <
int dim,
int spacedim>
4252 if (format == Default)
4254 const std::string::size_type slashpos = filename.find_last_of(
'/');
4255 const std::string::size_type dotpos = filename.find_last_of(
'.');
4256 if (dotpos < filename.size() &&
4257 (dotpos > slashpos || slashpos == std::string::npos))
4259 std::string ext = filename.substr(dotpos + 1);
4260 format = parse_format(ext);
4264 if (format == assimp)
4266 read_assimp(filename);
4268 else if (format == exodusii)
4270 read_exodusii(filename);
4274 std::ifstream in(filename);
4280template <
int dim,
int spacedim>
4284 if (format == Default)
4285 format = default_format;
4327 ExcMessage(
"There is no read_assimp(istream &) function. "
4328 "Use the read_assimp(string &filename, ...) "
4329 "functions, instead."));
4334 ExcMessage(
"There is no read_exodusii(istream &) function. "
4335 "Use the read_exodusii(string &filename, ...) "
4336 "function, instead."));
4347template <
int dim,
int spacedim>
4376 return ".unknown_format";
4382template <
int dim,
int spacedim>
4386 if (format_name ==
"dbmesh")
4389 if (format_name ==
"exodusii")
4392 if (format_name ==
"msh")
4395 if (format_name ==
"unv")
4398 if (format_name ==
"vtk")
4401 if (format_name ==
"vtu")
4405 if (format_name ==
"inp")
4408 if (format_name ==
"ucd")
4411 if (format_name ==
"xda")
4414 if (format_name ==
"tecplot")
4417 if (format_name ==
"dat")
4420 if (format_name ==
"plt")
4439template <
int dim,
int spacedim>
4443 return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
4450 template <
int dim,
int spacedim>
4451 Abaqus_to_UCD<dim, spacedim>::Abaqus_to_UCD()
4464 from_string(T &t,
const std::string &s, std::ios_base &(*f)(std::ios_base &))
4466 std::istringstream iss(s);
4467 return !(iss >> f >> t).fail();
4474 extract_int(
const std::string &s)
4477 for (
const char c : s)
4479 if (std::isdigit(c) != 0)
4486 from_string(number, tmp, std::dec);
4492 template <
int dim,
int spacedim>
4494 Abaqus_to_UCD<dim, spacedim>::read_in_abaqus(std::istream &input_stream)
4503 while (std::getline(input_stream, line))
4506 std::transform(line.begin(),
4509 static_cast<int (*)(
int)
>(std::toupper));
4511 if (line.compare(
"*HEADING") == 0 || line.compare(0, 2,
"**") == 0 ||
4512 line.compare(0, 5,
"*PART") == 0)
4515 while (std::getline(input_stream, line))
4521 else if (line.compare(0, 5,
"*NODE") == 0)
4530 while (std::getline(input_stream, line))
4535 std::vector<double> node(spacedim + 1);
4537 std::istringstream iss(line);
4539 for (
unsigned int i = 0; i < spacedim + 1; ++i)
4540 iss >> node[i] >> comma;
4542 node_list.push_back(node);
4545 else if (line.compare(0, 8,
"*ELEMENT") == 0)
4560 const std::string before_material =
"ELSET=EB";
4561 const std::size_t idx = line.find(before_material);
4562 if (idx != std::string::npos)
4564 from_string(material,
4565 line.substr(idx + before_material.size()),
4571 while (std::getline(input_stream, line))
4576 std::istringstream iss(line);
4582 const unsigned int n_data_per_cell =
4584 std::vector<double> cell(n_data_per_cell);
4585 for (
unsigned int i = 0; i < n_data_per_cell; ++i)
4586 iss >> cell[i] >> comma;
4589 cell[0] =
static_cast<double>(material);
4590 cell_list.push_back(cell);
4593 else if (line.compare(0, 8,
"*SURFACE") == 0)
4604 const std::string name_key =
"NAME=";
4605 const std::size_t name_idx_start =
4606 line.find(name_key) + name_key.size();
4607 std::size_t name_idx_end = line.find(
',', name_idx_start);
4608 if (name_idx_end == std::string::npos)
4610 name_idx_end = line.size();
4612 const int b_indicator = extract_int(
4613 line.substr(name_idx_start, name_idx_end - name_idx_start));
4620 while (std::getline(input_stream, line))
4626 std::transform(line.begin(),
4629 static_cast<int (*)(
int)
>(std::toupper));
4634 std::istringstream iss(line);
4641 std::vector<double> quad_node_list;
4642 const std::string elset_name = line.substr(0, line.find(
','));
4643 if (elsets_list.count(elset_name) != 0)
4647 iss >> stmp >> temp >> face_number;
4649 const std::vector<int> cells = elsets_list[elset_name];
4650 for (
const int cell : cells)
4654 get_global_node_numbers(el_idx, face_number);
4655 quad_node_list.insert(quad_node_list.begin(),
4658 face_list.push_back(quad_node_list);
4665 iss >> el_idx >> comma >> temp >> face_number;
4667 get_global_node_numbers(el_idx, face_number);
4668 quad_node_list.insert(quad_node_list.begin(), b_indicator);
4670 face_list.push_back(quad_node_list);
4674 else if (line.compare(0, 6,
"*ELSET") == 0)
4678 std::string elset_name;
4680 const std::string elset_key =
"*ELSET, ELSET=";
4681 const std::size_t idx = line.find(elset_key);
4682 if (idx != std::string::npos)
4684 const std::string comma =
",";
4685 const std::size_t first_comma = line.find(comma);
4686 const std::size_t second_comma =
4687 line.find(comma, first_comma + 1);
4688 const std::size_t elset_name_start =
4689 line.find(elset_key) + elset_key.size();
4690 elset_name = line.substr(elset_name_start,
4691 second_comma - elset_name_start);
4701 std::vector<int> elements;
4702 const std::size_t generate_idx = line.find(
"GENERATE");
4703 if (generate_idx != std::string::npos)
4706 std::getline(input_stream, line);
4707 std::istringstream iss(line);
4717 iss >> elid_start >> comma >> elid_end;
4721 "While reading an ABAQUS file, the reader "
4722 "expected a comma but found a <") +
4723 comma +
"> in the line <" + line +
">."));
4725 elid_start <= elid_end,
4728 "While reading an ABAQUS file, the reader encountered "
4729 "a GENERATE statement in which the upper bound <") +
4731 "> for the element numbers is not larger or equal "
4732 "than the lower bound <" +
4736 if (iss.rdbuf()->in_avail() != 0)
4737 iss >> comma >> elis_step;
4741 "While reading an ABAQUS file, the reader "
4742 "expected a comma but found a <") +
4743 comma +
"> in the line <" + line +
">."));
4745 for (
int i = elid_start; i <= elid_end; i += elis_step)
4746 elements.push_back(i);
4747 elsets_list[elset_name] = elements;
4749 std::getline(input_stream, line);
4754 while (std::getline(input_stream, line))
4759 std::istringstream iss(line);
4764 iss >> elid >> comma;
4769 "While reading an ABAQUS file, the reader "
4770 "expected a comma but found a <") +
4771 comma +
"> in the line <" + line +
">."));
4773 elements.push_back(elid);
4777 elsets_list[elset_name] = elements;
4782 else if (line.compare(0, 5,
"*NSET") == 0)
4785 while (std::getline(input_stream, line))
4791 else if (line.compare(0, 14,
"*SOLID SECTION") == 0)
4795 const std::string elset_key =
"ELSET=";
4796 const std::size_t elset_start =
4797 line.find(
"ELSET=") + elset_key.size();
4798 const std::size_t elset_end = line.find(
',', elset_start + 1);
4799 const std::string elset_name =
4800 line.substr(elset_start, elset_end - elset_start);
4805 const std::string material_key =
"MATERIAL=";
4806 const std::size_t last_equal =
4807 line.find(
"MATERIAL=") + material_key.size();
4808 const std::size_t material_id_start = line.find(
'-', last_equal);
4810 from_string(material_id,
4811 line.substr(material_id_start + 1),
4815 const std::vector<int> &elset_cells = elsets_list[elset_name];
4816 for (
const int elset_cell : elset_cells)
4818 const int cell_id = elset_cell - 1;
4826 template <
int dim,
int spacedim>
4828 Abaqus_to_UCD<dim, spacedim>::get_global_node_numbers(
4829 const int face_cell_no,
4830 const int face_cell_face_no)
const
4840 if (face_cell_face_no == 1)
4842 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4843 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4845 else if (face_cell_face_no == 2)
4847 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4848 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4850 else if (face_cell_face_no == 3)
4852 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4853 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4855 else if (face_cell_face_no == 4)
4857 quad_node_list[0] = cell_list[face_cell_no - 1][4];
4858 quad_node_list[1] = cell_list[face_cell_no - 1][1];
4868 if (face_cell_face_no == 1)
4870 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4871 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4872 quad_node_list[2] = cell_list[face_cell_no - 1][3];
4873 quad_node_list[3] = cell_list[face_cell_no - 1][2];
4875 else if (face_cell_face_no == 2)
4877 quad_node_list[0] = cell_list[face_cell_no - 1][5];
4878 quad_node_list[1] = cell_list[face_cell_no - 1][8];
4879 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4880 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4882 else if (face_cell_face_no == 3)
4884 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4885 quad_node_list[1] = cell_list[face_cell_no - 1][2];
4886 quad_node_list[2] = cell_list[face_cell_no - 1][6];
4887 quad_node_list[3] = cell_list[face_cell_no - 1][5];
4889 else if (face_cell_face_no == 4)
4891 quad_node_list[0] = cell_list[face_cell_no - 1][2];
4892 quad_node_list[1] = cell_list[face_cell_no - 1][3];
4893 quad_node_list[2] = cell_list[face_cell_no - 1][7];
4894 quad_node_list[3] = cell_list[face_cell_no - 1][6];
4896 else if (face_cell_face_no == 5)
4898 quad_node_list[0] = cell_list[face_cell_no - 1][3];
4899 quad_node_list[1] = cell_list[face_cell_no - 1][4];
4900 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4901 quad_node_list[3] = cell_list[face_cell_no - 1][7];
4903 else if (face_cell_face_no == 6)
4905 quad_node_list[0] = cell_list[face_cell_no - 1][1];
4906 quad_node_list[1] = cell_list[face_cell_no - 1][5];
4907 quad_node_list[2] = cell_list[face_cell_no - 1][8];
4908 quad_node_list[3] = cell_list[face_cell_no - 1][4];
4921 return quad_node_list;
4924 template <
int dim,
int spacedim>
4926 Abaqus_to_UCD<dim, spacedim>::write_out_avs_ucd(std::ostream &output)
const
4935 const boost::io::ios_base_all_saver formatting_saver(output);
4939 output <<
"# Abaqus to UCD mesh conversion" << std::endl;
4940 output <<
"# Mesh type: AVS UCD" << std::endl;
4969 output << node_list.size() <<
"\t" << (cell_list.size() + face_list.size())
4970 <<
"\t0\t0\t0" << std::endl;
4973 output.precision(8);
4977 for (
const auto &node : node_list)
4980 output << node[0] <<
"\t";
4983 output.setf(std::ios::scientific, std::ios::floatfield);
4984 for (
unsigned int jj = 1; jj < spacedim + 1; ++jj)
4987 if (
std::abs(node[jj]) > tolerance)
4988 output << static_cast<double>(node[jj]) <<
"\t";
4990 output << 0.0 <<
"\t";
4993 output << 0.0 <<
"\t";
4995 output << std::endl;
4996 output.unsetf(std::ios::floatfield);
5000 for (
unsigned int ii = 0; ii < cell_list.size(); ++ii)
5002 output << ii + 1 <<
"\t" << cell_list[ii][0] <<
"\t"
5003 << (dim == 2 ?
"quad" :
"hex") <<
"\t";
5004 for (
unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
5006 output << cell_list[ii][jj] <<
"\t";
5008 output << std::endl;
5012 for (
unsigned int ii = 0; ii < face_list.size(); ++ii)
5014 output << ii + 1 <<
"\t" << face_list[ii][0] <<
"\t"
5015 << (dim == 2 ?
"line" :
"quad") <<
"\t";
5016 for (
unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
5018 output << face_list[ii][jj] <<
"\t";
5020 output << std::endl;
5027#include "grid/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)
const std::map< std::string, Vector< double > > & get_cell_data() const
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 exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const types::geometric_orientation face_orientation) 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)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_FALLTHROUGH
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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 AssertThrowExodusII(error_code)
#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)
std::vector< index_type > data
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Pyramid
constexpr ReferenceCell Wedge
constexpr ReferenceCell Vertex
constexpr ReferenceCell Invalid
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
constexpr 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)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
constexpr types::manifold_id flat_manifold_id
constexpr types::material_id invalid_material_id
constexpr types::geometric_orientation default_geometric_orientation
::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