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 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)
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