diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index 19dcb99..4579d3f 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -13,7 +13,9 @@ #include #include +#include #include +#include #include #include "helper/Distributor.h" @@ -43,6 +45,7 @@ template class ParallelGMSHReader { logError() << meshFile << std::endl << parser.getErrorMessage(); } builder_.postprocess(); + validateIdentifyTopology(); convertBoundaryConditions(); nVertices_ = builder_.vertices.size(); @@ -77,6 +80,60 @@ template class ParallelGMSHReader { } private: + /** + * Checks that periodic vertex identification does not create non-manifold + * topology (i.e. identified faces shared by more than two cells). + */ + void validateIdentifyTopology() const { + if (builder_.identify.empty()) { + return; + } + + using FaceKey = std::array; + struct FaceKeyHash { + std::size_t operator()(const FaceKey& key) const { + return (key[0] * 73856093u) ^ (key[1] * 19349663u) ^ (key[2] * 83492791u); + } + }; + + std::unordered_map faceMultiplicity; + faceMultiplicity.reserve(builder_.elements.size() * (Dim + 1u)); + + bool hasNonManifoldFace = false; + FaceKey exampleFace{0, 0, 0}; + std::size_t exampleCount = 0; + + for (const auto& element : builder_.elements) { + for (int localFctNo = 0; localFctNo < Dim + 1u; ++localFctNo) { + FaceKey face{}; + for (int localNo = 0; localNo < Dim; ++localNo) { + const auto localNoElement = Facet2Nodes[localFctNo][localNo]; + const auto vertex = element[localNoElement]; + face[localNo] = builder_.identify[vertex]; + } + std::sort(face.begin(), face.end()); + + const auto count = ++faceMultiplicity[face]; + if (count > 2) { + hasNonManifoldFace = true; + exampleFace = face; + exampleCount = count; + } + } + } + + if (hasNonManifoldFace) { + logError() + << "Invalid periodic identify topology: multiple distinct faces collapse to the same" + << "canonical face under indentify, yielding a face shared by more than two cells." + << "Example canonical face ids: (" << exampleFace[0] << ", " << exampleFace[1] << ", " + << exampleFace[2] << ") with " << exampleCount + << ". Verify periodic entity pairing and the periodic section in the mesh file, " + << "and regenerate the mesh with consistent periodic definitions"; + MPI_Abort(comm_, EXIT_FAILURE); + } + } + /** * GMSH stores boundary conditions on a surface mesh whereas SeisSol expects * boundary conditions to be stored per element. In the following we convert