diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp index 60c38bfb84fb975c0e2cdc3b07d1d97e39eb0021..d83e56ebe2b2e9e761fab8d5221da13b50473d32 100644 --- a/src/mesh/ConnectivityBuilderBase.cpp +++ b/src/mesh/ConnectivityBuilderBase.cpp @@ -116,6 +116,30 @@ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityD face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed())); break; } + case CellType::Prism: { + // face 0 + Face f0({cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector); + face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed())); + + // face 1 + Face f1({cell_nodes[3], cell_nodes[4], cell_nodes[5]}, node_number_vector); + face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed())); + + // face 2 + Face f2({cell_nodes[1], cell_nodes[2], cell_nodes[5], cell_nodes[4]}, node_number_vector); + face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed())); + + // face 3 + Face f3({cell_nodes[0], cell_nodes[1], cell_nodes[4], cell_nodes[3]}, node_number_vector); + face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed())); + + // face 4 + Face f4({cell_nodes[2], cell_nodes[0], cell_nodes[3], cell_nodes[5]}, node_number_vector); + face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed())); + + cell_nb_faces[j] = 5; + break; + } case CellType::Pyramid: { cell_nb_faces[j] = cell_nodes.size(); std::vector<unsigned int> base_nodes; @@ -332,6 +356,22 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector); break; } + case CellType::Prism: { + constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}}; + std::vector<unsigned int> cell_edge_vector; + cell_edge_vector.reserve(9); + for (int i_edge = 0; i_edge < 9; ++i_edge) { + const auto e = local_edge[i_edge]; + Edge edge{{cell_nodes[e[0]], cell_nodes[e[1]]}, node_number_vector}; + auto i = edge_id_map.find(edge); + if (i == edge_id_map.end()) { + throw NormalError("could not find this edge"); + } + cell_edge_vector.push_back(i->second); + } + descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector); + break; + } case CellType::Pyramid: { const size_t number_of_edges = 2 * cell_nodes.size(); std::vector<unsigned int> base_nodes; diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp index 5ecd752dccc8cc5aade741c3fd058bf1c8fec9c3..73853e01ad48a6d75a12d62771c14f25549978e8 100644 --- a/src/mesh/GmshReader.cpp +++ b/src/mesh/GmshReader.cpp @@ -248,9 +248,20 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData& descriptor.cell_number_vector[jh] = gmsh_data.__hexahedra_number[j]; } + const size_t nb_prisms = gmsh_data.__prisms.size(); + for (size_t j = 0; j < nb_prisms; ++j) { + const size_t jp = nb_tetrahedra + nb_hexahedra + j; + descriptor.cell_to_node_vector[jp].resize(6); + for (int r = 0; r < 6; ++r) { + descriptor.cell_to_node_vector[jp][r] = gmsh_data.__prisms[j][r]; + } + descriptor.cell_type_vector[jp] = CellType::Prism; + descriptor.cell_number_vector[jp] = gmsh_data.__prisms_number[j]; + } + const size_t nb_pyramids = gmsh_data.__pyramids.size(); for (size_t j = 0; j < nb_pyramids; ++j) { - const size_t jh = nb_tetrahedra + nb_hexahedra + j; + const size_t jh = nb_tetrahedra + nb_hexahedra + nb_prisms + j; descriptor.cell_to_node_vector[jh].resize(5); for (int r = 0; r < 5; ++r) { descriptor.cell_to_node_vector[jh][r] = gmsh_data.__pyramids[j][r]; @@ -503,6 +514,8 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename) __keywordList["$EndElements"] = ENDELEMENTS; __keywordList["$PhysicalNames"] = PHYSICALNAMES; __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES; + __keywordList["$Periodic"] = PERIODIC; + __keywordList["$EndPeriodic"] = ENDPERIODIC; __numberOfPrimitiveNodes.resize(16); __numberOfPrimitiveNodes[0] = 2; // edge @@ -532,7 +545,7 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename) __primitivesNames[4] = "hexahedra"; __supportedPrimitives[4] = true; __primitivesNames[5] = "prisms"; - __supportedPrimitives[5] = false; + __supportedPrimitives[5] = true; __primitivesNames[6] = "pyramids"; __supportedPrimitives[6] = true; __primitivesNames[7] = "second order edges"; @@ -718,6 +731,24 @@ GmshReader::__readPhysicalNames2_2() } } +void +GmshReader::__readPeriodic2_2() +{ + // This is just a compatibility reading, periodic information is not + // used + const int number_of_periodic = this->_getInteger(); + for (int i = 0; i < number_of_periodic; ++i) { + [[maybe_unused]] const int dimension = this->_getInteger(); + [[maybe_unused]] const int id = this->_getInteger(); + [[maybe_unused]] const int master_id = this->_getInteger(); + [[maybe_unused]] const int nb_corresponding_nodes = this->_getInteger(); + for (int i_node = 0; i_node < nb_corresponding_nodes; ++i_node) { + [[maybe_unused]] const int node_id = this->_getInteger(); + [[maybe_unused]] const int master_id = this->_getInteger(); + } + } +} + // std::shared_ptr<IConnectivity> // GmshReader::_buildConnectivity3D(const size_t nb_cells) // { @@ -1026,6 +1057,12 @@ GmshReader::__proceedData() m_mesh_data.__hexahedra_number.resize(elementNumber[i]); break; } + case 5: { // prism + m_mesh_data.__prisms = Array<GmshData::Prism>(elementNumber[i]); + m_mesh_data.__prisms_ref.resize(elementNumber[i]); + m_mesh_data.__prisms_number.resize(elementNumber[i]); + break; + } case 6: { // pyramid m_mesh_data.__pyramids = Array<GmshData::Pyramid>(elementNumber[i]); m_mesh_data.__pyramids_ref.resize(elementNumber[i]); @@ -1039,7 +1076,6 @@ GmshReader::__proceedData() break; } // Unsupported entities - case 5: // prism case 7: // second order edge case 8: // second order triangle case 9: // second order quadrangle @@ -1169,6 +1205,24 @@ GmshReader::__proceedData() hexahedronNumber++; // one more hexahedron break; } + case 5: { // prism + int& prism_number = elementNumber[5]; + const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]]; + const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]]; + const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]]; + const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]]; + const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]]; + const int f = m_mesh_data.__verticesCorrepondance[elementVertices[5]]; + if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0)) { + throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) + + " [bad vertices definition]"); + } + m_mesh_data.__prisms[prism_number] = GmshData::Prism(a, b, c, d, e, f); + m_mesh_data.__prisms_ref[prism_number] = m_mesh_data.__references[i]; + m_mesh_data.__prisms_number[prism_number] = m_mesh_data.__elementNumber[i]; + prism_number++; // one more prism + break; + } case 6: { // pyramid int& pyramid_number = elementNumber[6]; @@ -1197,7 +1251,6 @@ GmshReader::__proceedData() point_number++; break; } - case 5: // prism case 7: // second order edge case 8: // second order triangle case 9: // second order quadrangle @@ -1352,7 +1405,15 @@ GmshReader::__readGmshFormat2_2() case PHYSICALNAMES: { this->__readPhysicalNames2_2(); if (this->__nextKeyword().second != ENDPHYSICALNAMES) { - throw NormalError("reading file '" + m_filename + "': expecting $EndNodes, '" + kw.first + "' was found"); + throw NormalError("reading file '" + m_filename + "': expecting $EndPhysicalNames, '" + kw.first + + "' was found"); + } + break; + } + case PERIODIC: { + this->__readPeriodic2_2(); + if (this->__nextKeyword().second != ENDPERIODIC) { + throw NormalError("reading file '" + m_filename + "': expecting $EndPeriodic, '" + kw.first + "' was found"); } break; } diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp index 36e21508f10eda492ec1128816215d5a4d8c624d..860695ef57143b1b67575f54da95c12f8c7c2a01 100644 --- a/src/mesh/GmshReader.hpp +++ b/src/mesh/GmshReader.hpp @@ -109,6 +109,11 @@ class GmshReader : public MeshBuilderBase std::vector<int> __pyramids_ref; std::vector<int> __pyramids_number; + using Prism = TinyVector<6, unsigned int>; + Array<Prism> __prisms; + std::vector<int> __prisms_ref; + std::vector<int> __prisms_number; + using Hexahedron = TinyVector<8, unsigned int>; Array<Hexahedron> __hexahedra; std::vector<int> __hexahedra_ref; @@ -212,6 +217,8 @@ class GmshReader : public MeshBuilderBase ENDELEMENTS, PHYSICALNAMES, ENDPHYSICALNAMES, + PERIODIC, + ENDPERIODIC, Unknown, EndOfFile @@ -268,6 +275,12 @@ class GmshReader : public MeshBuilderBase */ void __readPhysicalNames2_2(); + /** + * Reads periodic + * + */ + void __readPeriodic2_2(); + public: GmshReader(const std::string& filename); ~GmshReader() = default;