diff --git a/src/input/MeshData.h b/src/input/MeshData.h index 230c209..286c764 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -19,21 +19,25 @@ class MeshData { public: virtual ~MeshData() = default; - virtual std::size_t cellCount() = 0; - virtual std::size_t vertexCount() = 0; + virtual std::size_t cellCount() const = 0; + virtual std::size_t vertexCount() const = 0; virtual const std::vector& connectivity() = 0; virtual const std::vector& geometry() = 0; virtual const std::vector& group() = 0; virtual const std::vector& boundary() = 0; + + // TODO: generalize? + virtual std::size_t vertexSize() const { return 3; } + virtual std::size_t cellSize() const { return 4; } }; class FullStorageMeshData : public MeshData { public: virtual ~FullStorageMeshData() = default; - virtual std::size_t cellCount() { return cellCountValue; } - virtual std::size_t vertexCount() { return vertexCountValue; } + virtual std::size_t cellCount() const { return cellCountValue; } + virtual std::size_t vertexCount() const { return vertexCountValue; } virtual const std::vector& connectivity() { return connectivityData; } virtual const std::vector& geometry() { return geometryData; } @@ -55,7 +59,7 @@ class FullStorageMeshData : public MeshData { void setBoundary(std::size_t cell, int face, int value) { if (bndShift < 0) { - boundaryData[cell * 4 + face] = value; + boundaryData[cell * (vertexSize() + 1) + face] = value; } else { constexpr auto i32limit = std::numeric_limits::max(); constexpr auto i64limit = std::numeric_limits::max(); @@ -72,11 +76,11 @@ class FullStorageMeshData : public MeshData { cellCountValue = cellCount; vertexCountValue = vertexCount; - connectivityData.resize(cellCount * 4); - geometryData.resize(vertexCount * 3); + connectivityData.resize(cellCount * cellSize()); + geometryData.resize(vertexCount * vertexSize()); groupData.resize(cellCount); if (bndShift < 0) { - boundaryData.resize(cellCount * 4); + boundaryData.resize(cellCount * (vertexSize() + 1)); } else { boundaryData.resize(cellCount); } diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 362b26e..bab36a9 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -31,6 +31,9 @@ template class SerialMeshFile : public FullStorageMeshData { public: virtual ~SerialMeshFile() = default; + std::size_t vertexSize() const override { return T::Dim; } + std::size_t cellSize() const override { return puml::nodeCount(T::Dim, T::Order); } + private: #ifdef PARALLEL MPI_Comm m_comm; @@ -89,11 +92,11 @@ template class SerialMeshFile : public FullStorageMeshData { m_meshReader.readGroups(groupData.data()); logInfo(m_rank) << "Read boundary conditions"; - std::vector preBoundaryData(nLocalElements * 4); + std::vector preBoundaryData(nLocalElements * (vertexSize() + 1)); m_meshReader.readBoundaries(preBoundaryData.data()); for (std::size_t i = 0; i < nLocalElements; ++i) { - for (int j = 0; j < 4; ++j) { - setBoundary(i, j, preBoundaryData[4 * i + j]); + for (int j = 0; j < (vertexSize() + 1); ++j) { + setBoundary(i, j, preBoundaryData[(vertexSize() + 1) * i + j]); } } } diff --git a/src/meshreader/GMSHBuilder.h b/src/meshreader/GMSHBuilder.h index f4c35c9..778edbe 100644 --- a/src/meshreader/GMSHBuilder.h +++ b/src/meshreader/GMSHBuilder.h @@ -1,37 +1,147 @@ #ifndef GMSHBUILDER_20201014_H #define GMSHBUILDER_20201014_H +#include #include #include #include +#include "utils/logger.h" + #include "third_party/GMSHMeshBuilder.h" namespace puml { -template struct GMSHSimplexType {}; +constexpr std::size_t nodeCount(std::size_t D, std::size_t Order) { + std::size_t count = 1; + for (std::size_t d = 1; d <= D; ++d) { + count *= (Order + D - d + 1); + count /= d; + } + return count; +} + +template typename ParamT, typename... Args> +BaseT* makePointerI(std::size_t param, Args&&... args) { + if constexpr (P < 11) { + if (param == P) { + return new ParamT

(std::forward(args)...); + } else { + return makePointerI

(param, std::forward(args)...); + } + } else { + logError() << "Order too high:" << param; + return nullptr; + } +} + +template typename ParamT, typename... Args> +BaseT* makePointer(std::size_t param, Args&&... args) { + return makePointerI<1, BaseT, ParamT>(param, std::forward(args)...); +} + +template struct GMSHSimplexType {}; // clang format seems to mess up the following lines from time to time, hence we disable it -// clang-format off -template <> struct GMSHSimplexType<0u> { +template struct GMSHSimplexType<0U, Order> { static constexpr long type = 15; -}; // 15 -> point -template <> struct GMSHSimplexType<1u> { +}; +template <> struct GMSHSimplexType<1U, 1U> { static constexpr long type = 1; -}; // 1 -> line -template <> struct GMSHSimplexType<2u> { +}; +template <> struct GMSHSimplexType<2U, 1U> { static constexpr long type = 2; -}; // 2 -> triangle -template <> struct GMSHSimplexType<3u> { +}; +template <> struct GMSHSimplexType<3U, 1U> { static constexpr long type = 4; -}; // 4 -> tetrahedron -// clang-format on +}; +template <> struct GMSHSimplexType<1U, 2U> { + static constexpr long type = 8; +}; +template <> struct GMSHSimplexType<2U, 2U> { + static constexpr long type = 9; +}; +template <> struct GMSHSimplexType<3U, 2U> { + static constexpr long type = 11; +}; +template <> struct GMSHSimplexType<1U, 3U> { + static constexpr long type = 26; +}; +template <> struct GMSHSimplexType<1U, 4U> { + static constexpr long type = 27; +}; +template <> struct GMSHSimplexType<1U, 5U> { + static constexpr long type = 28; +}; +template <> struct GMSHSimplexType<1U, 6U> { + static constexpr long type = 62; +}; +template <> struct GMSHSimplexType<1U, 7U> { + static constexpr long type = 63; +}; +template <> struct GMSHSimplexType<1U, 8U> { + static constexpr long type = 64; +}; +template <> struct GMSHSimplexType<1U, 9U> { + static constexpr long type = 65; +}; +template <> struct GMSHSimplexType<1U, 10U> { + static constexpr long type = 66; +}; +template <> struct GMSHSimplexType<2U, 3U> { + static constexpr long type = 21; +}; +template <> struct GMSHSimplexType<2U, 4U> { + static constexpr long type = 23; +}; +template <> struct GMSHSimplexType<2U, 5U> { + static constexpr long type = 25; +}; +template <> struct GMSHSimplexType<2U, 6U> { + static constexpr long type = 42; +}; +template <> struct GMSHSimplexType<2U, 7U> { + static constexpr long type = 43; +}; +template <> struct GMSHSimplexType<2U, 8U> { + static constexpr long type = 44; +}; +template <> struct GMSHSimplexType<2U, 9U> { + static constexpr long type = 45; +}; +template <> struct GMSHSimplexType<2U, 10U> { + static constexpr long type = 46; +}; +template <> struct GMSHSimplexType<3U, 3U> { + static constexpr long type = 29; +}; +template <> struct GMSHSimplexType<3U, 4U> { + static constexpr long type = 30; +}; +template <> struct GMSHSimplexType<3U, 5U> { + static constexpr long type = 31; +}; +template <> struct GMSHSimplexType<3U, 6U> { + static constexpr long type = 71; +}; +template <> struct GMSHSimplexType<3U, 7U> { + static constexpr long type = 72; +}; +template <> struct GMSHSimplexType<3U, 8U> { + static constexpr long type = 73; +}; +template <> struct GMSHSimplexType<3U, 9U> { + static constexpr long type = 74; +}; +template <> struct GMSHSimplexType<3U, 10U> { + static constexpr long type = 75; +}; -template class GMSHBuilder : public tndm::GMSHMeshBuilder { +template class GMSHBuilder : public tndm::GMSHMeshBuilder { public: using vertex_t = std::array; - using element_t = std::array; - using facet_t = std::array; + using element_t = std::array; + using facet_t = std::array; std::vector vertices; std::vector elements; @@ -56,16 +166,16 @@ template class GMSHBuilder : public tndm::GMSHMeshBuilder { bcs.reserve(numElements); } void addElement(long type, long tag, long* node, std::size_t numNodes) { - if (type == GMSHSimplexType::type) { - assert(numNodes == D + 1u); + if (type == GMSHSimplexType::type) { element_t elem; - std::copy(node, node + D + 1u, elem.begin()); + assert(numNodes == elem.size()); + std::copy_n(node, elem.size(), elem.begin()); elements.emplace_back(elem); groups.emplace_back(static_cast(tag)); - } else if (type == GMSHSimplexType::type) { - assert(numNodes == D); + } else if (type == GMSHSimplexType::type) { facet_t elem; - std::copy(node, node + D, elem.begin()); + assert(numNodes == elem.size()); + std::copy_n(node, elem.size(), elem.begin()); facets.emplace_back(elem); bcs.emplace_back(static_cast(tag)); } diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index 99388b8..53f8031 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -9,6 +9,7 @@ #include #include +#include #include #include "aux/Distributor.h" @@ -16,8 +17,9 @@ namespace puml { -template class ParallelGMSHReader { +template class ParallelGMSHReader { public: + constexpr static std::size_t Order = OrderP; constexpr static std::size_t Dim = 3; using bc_t = std::array; constexpr static std::size_t Facet2Nodes[Dim + 1][Dim] = { @@ -49,11 +51,12 @@ template class ParallelGMSHReader { [[nodiscard]] std::size_t nVertices() const { return nVertices_; } [[nodiscard]] std::size_t nElements() const { return nElements_; } void readElements(std::size_t* elements) const { - static_assert(sizeof(GMSHBuilder::element_t) == (Dim + 1) * sizeof(std::size_t)); - scatter(builder_.elements.data()->data(), elements, nElements(), Dim + 1); + static_assert(sizeof(typename GMSHBuilder::element_t) == + nodeCount(Dim, Order) * sizeof(std::size_t)); + scatter(builder_.elements.data()->data(), elements, nElements(), nodeCount(Dim, Order)); } void readVertices(double* vertices) const { - static_assert(sizeof(GMSHBuilder::vertex_t) == Dim * sizeof(double)); + static_assert(sizeof(typename GMSHBuilder::vertex_t) == Dim * sizeof(double)); scatter(builder_.vertices.data()->data(), vertices, nVertices(), Dim); } void readBoundaries(int* boundaries) const { @@ -157,7 +160,7 @@ template class ParallelGMSHReader { } MPI_Comm comm_; - GMSHBuilder builder_; + GMSHBuilder builder_; std::vector bcs_; std::size_t nVertices_ = 0; std::size_t nElements_ = 0; diff --git a/src/meshreader/ParallelMeshReader.h b/src/meshreader/ParallelMeshReader.h index 8d52036..f7b92c1 100644 --- a/src/meshreader/ParallelMeshReader.h +++ b/src/meshreader/ParallelMeshReader.h @@ -36,6 +36,8 @@ template class ParallelMeshReader { R m_serialReader; public: + constexpr static std::size_t Order = 1; + constexpr static std::size_t Dim = 3; ParallelMeshReader(MPI_Comm comm = MPI_COMM_WORLD) : m_nVertices(0), m_nElements(0), m_nBoundaries(0), m_comm(comm) { init(); diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 5a46c04..d396cc1 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -10,6 +10,7 @@ * @author Sebastian Rettenberger */ +#include "meshreader/GMSHBuilder.h" #include #include @@ -169,6 +170,11 @@ static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, checkH5Err(H5Pclose(h5dxlist)); } +template +using SMF2 = SerialMeshFile>; +template +using SMF4 = SerialMeshFile>; + int main(int argc, char* argv[]) { int rank = 0; int processes = 1; @@ -227,6 +233,7 @@ int main(int argc, char* argv[]) { args.addOption("sim_log", 0, "Create SimModSuite log", utils::Args::Required, false); args.addAdditionalOption("input", "Input file (mesh or model)"); args.addAdditionalOption("output", "Output parallel unstructured mesh file", false); + args.addOption("order", 'o', "Mesh order", utils::Args::Optional, false); if (args.parse(argc, argv, rank == 0) != utils::Args::Success) return 1; @@ -308,6 +315,8 @@ int main(int argc, char* argv[]) { xdmfFile.append(".xdmf"); } + const auto meshOrder = args.getArgument("order", 1); + // Create/read the mesh MeshData* meshInput = nullptr; switch (args.getArgument("source", 0)) { @@ -317,13 +326,11 @@ int main(int argc, char* argv[]) { break; case 1: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; - meshInput = - new SerialMeshFile>(inputFile, faceOffset); + meshInput = puml::makePointer(meshOrder, inputFile, faceOffset); break; case 2: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; - meshInput = - new SerialMeshFile>(inputFile, faceOffset); + meshInput = puml::makePointer(meshOrder, inputFile, faceOffset); break; case 3: #ifdef USE_NETCDF @@ -391,6 +398,8 @@ int main(int argc, char* argv[]) { MPI_Allreduce(MPI_IN_PLACE, globalSize, 2, tndm::mpi_type_t(), MPI_SUM, MPI_COMM_WORLD); + logInfo(rank) << "Coordinates in vertex:" << meshInput->vertexSize(); + logInfo(rank) << "Vertices in cell:" << meshInput->cellSize(); logInfo(rank) << "Total cell count:" << globalSize[0]; logInfo(rank) << "Total vertex count:" << globalSize[1]; @@ -424,13 +433,13 @@ int main(int argc, char* argv[]) { logInfo(rank) << "Writing cells"; writeH5Data(meshInput->connectivity(), h5file, "connect", mesh, 3, H5T_NATIVE_UINT64, H5T_STD_U64LE, chunksize, localSize[0], globalSize[0], reduceInts, - filterEnable, filterChunksize, 4); + filterEnable, filterChunksize, meshInput->cellSize()); // Vertices logInfo(rank) << "Writing vertices"; writeH5Data(meshInput->geometry(), h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, H5T_IEEE_F64LE, chunksize, localSize[1], globalSize[1], reduceInts, - filterEnable, filterChunksize, 3); + filterEnable, filterChunksize, meshInput->vertexSize()); // Group information @@ -442,6 +451,10 @@ int main(int argc, char* argv[]) { // Write boundary condition logInfo(rank) << "Writing boundary condition"; + if (secondShape != NoSecondDim) { + // TODO: a bit ugly, but it works + secondShape = meshInput->vertexSize() + 1; + } if (boundaryFormatAttr == "i32") { writeH5Data(meshInput->boundary(), h5file, "boundary", mesh, 3, H5T_NATIVE_INT32, boundaryDatatype, chunksize, localSize[0], globalSize[0], reduceInts, diff --git a/src/third_party/GMSHParser.h b/src/third_party/GMSHParser.h index 855159d..99f598e 100644 --- a/src/third_party/GMSHParser.h +++ b/src/third_party/GMSHParser.h @@ -7,100 +7,196 @@ #include #include -#include "GMSHMeshBuilder.h" #include "GMSHLexer.h" +#include "GMSHMeshBuilder.h" namespace tndm { class GMSHParser { -public: - // Look-up table from gmsh type to number of nodes - static constexpr std::size_t NumNodes[] = { - 2, // line - 3, // triangle, - 4, // quadrangle, - 4, // tetrahedron - 8, // hexahedron - 6, // prism - 5, // pyramid - 3, // P1 line - 6, // P1 triangle - 9, // P1 quadrangle - 10, // P1 tetrahedron - 27, // P2 hexahedron - 18, // P2 prism - 14, // P2 pyramid - 1, // point - }; - static constexpr const char* ElementTypes[] = { - "line", - "triangle", - "quadrangle", - "tetrahedron", - "hexahedron", - "prism", - "pyramid", - "P1 line", - "P1 triangle", - "P1 quadrangle", - "P1 tetrahedron", - "P2 hexahedron", - "P2 prism", - "P2 pyramid", - "point", - }; + public: + // Look-up table from gmsh type to number of nodes + static constexpr std::size_t NumNodes[] = { + 2, // MSH_LIN_2 + 3, // MSH_TRI_3 + 4, // MSH_QUA_4 + 4, // MSH_TET_4 + 8, // MSH_HEX_8 + 6, // MSH_PRI_6 + 5, // MSH_PYR_5 + 3, // MSH_LIN_3 + 6, // MSH_TRI_6 + 9, // MSH_QUA_9 + 10, // MSH_TET_10 + 27, // MSH_HEX_27 + 18, // MSH_PRI_18 + 14, // MSH_PYR_14 + 1, // MSH_PNT + 8, // MSH_QUA_8 + 20, // MSH_HEX_20 + 15, // MSH_PRI_15 + 13, // MSH_PYR_13 + 9, // MSH_TRI_9 + 10, // MSH_TRI_10 + 12, // MSH_TRI_12 + 15, // MSH_TRI_15 + 15, // MSH_TRI_15I + 21, // MSH_TRI_21 + 4, // MSH_LIN_4 + 5, // MSH_LIN_5 + 6, // MSH_LIN_6 + 20, // MSH_TET_20 + 35, // MSH_TET_35 + 56, // MSH_TET_56 + 22, // MSH_TET_22 + 28, // MSH_TET_28 + 0, // MSH_POLYG_ + 0, // MSH_POLYH_ + 16, // MSH_QUA_16 + 25, // MSH_QUA_25 + 36, // MSH_QUA_36 + 12, // MSH_QUA_12 + 16, // MSH_QUA_16I + 20, // MSH_QUA_20 + 28, // MSH_TRI_28 + 36, // MSH_TRI_36 + 45, // MSH_TRI_45 + 55, // MSH_TRI_55 + 66, // MSH_TRI_66 + 49, // MSH_QUA_49 + 64, // MSH_QUA_64 + 81, // MSH_QUA_81 + 100, // MSH_QUA_100 + 121, // MSH_QUA_121 + 18, // MSH_TRI_18 + 21, // MSH_TRI_21I + 24, // MSH_TRI_24 + 27, // MSH_TRI_27 + 30, // MSH_TRI_30 + 24, // MSH_QUA_24 + 28, // MSH_QUA_28 + 32, // MSH_QUA_32 + 36, // MSH_QUA_36I + 40, // MSH_QUA_40 + 7, // MSH_LIN_7 + 8, // MSH_LIN_8 + 9, // MSH_LIN_9 + 10, // MSH_LIN_10 + 11, // MSH_LIN_11 + 0, // MSH_LIN_B + 0, // MSH_TRI_B + 0, // MSH_POLYG_B + 0, // MSH_LIN_C + 84, // MSH_TET_84 + 120, // MSH_TET_120 + 165, // MSH_TET_165 + 220, // MSH_TET_220 + 286, // MSH_TET_286 + 34, // MSH_TET_34 + 40, // MSH_TET_40 + 46, // MSH_TET_46 + 52, // MSH_TET_52 + 58, // MSH_TET_58 + 1, // MSH_LIN_1 + 1, // MSH_TRI_1 + 1, // MSH_QUA_1 + 1, // MSH_TET_1 + 1, // MSH_HEX_1 + 1, // MSH_PRI_1 + 40, // MSH_PRI_40 + 75, // MSH_PRI_75 + 64, // MSH_HEX_64 + 125, // MSH_HEX_125 + 216, // MSH_HEX_216 + 343, // MSH_HEX_343 + 512, // MSH_HEX_512 + 729, // MSH_HEX_729 + 1000, // MSH_HEX_1000 + 32, // MSH_HEX_32 + 44, // MSH_HEX_44 + 56, // MSH_HEX_56 + 68, // MSH_HEX_68 + 80, // MSH_HEX_80 + 92, // MSH_HEX_92 + 104, // MSH_HEX_104 + 126, // MSH_PRI_126 + 196, // MSH_PRI_196 + 288, // MSH_PRI_288 + 405, // MSH_PRI_405 + 550, // MSH_PRI_550 + 24, // MSH_PRI_24 + 33, // MSH_PRI_33 + 42, // MSH_PRI_42 + 51, // MSH_PRI_51 + 60, // MSH_PRI_60 + 69, // MSH_PRI_69 + 78, // MSH_PRI_78 + 30, // MSH_PYR_30 + 55, // MSH_PYR_55 + 91, // MSH_PYR_91 + 140, // MSH_PYR_140 + 204, // MSH_PYR_204 + 285, // MSH_PYR_285 + 385, // MSH_PYR_385 + 21, // MSH_PYR_21 + 29, // MSH_PYR_29 + 37, // MSH_PYR_37 + 45, // MSH_PYR_45 + 53, // MSH_PYR_53 + 61, // MSH_PYR_61 + 69, // MSH_PYR_69 + }; - explicit GMSHParser(GMSHMeshBuilder* builder) : builder(builder){}; - [[nodiscard]] std::string_view getErrorMessage() const { return errorMsg; } + explicit GMSHParser(GMSHMeshBuilder* builder) : builder(builder) {}; + [[nodiscard]] std::string_view getErrorMessage() const { return errorMsg; } - bool parseFile(std::string const& fileName) { - std::ifstream in(fileName); - if (!in.is_open()) { - return logError("Unable to open MSH file"); - } - lexer.setIStream(&in); - return parse_(); + bool parseFile(std::string const& fileName) { + std::ifstream in(fileName); + if (!in.is_open()) { + return logError("Unable to open MSH file"); } + lexer.setIStream(&in); + return parse_(); + } -protected: - template T logErrorAnnotated(std::string_view msg) { - std::stringstream ss; - ss << "GMSH parser error in line " << curLoc.line << " in column " << curLoc.col << ":\n"; - ss << '\t' << msg << '\n'; - errorMsg += ss.str(); - return {}; - } + protected: + template T logErrorAnnotated(std::string_view msg) { + std::stringstream ss; + ss << "GMSH parser error in line " << curLoc.line << " in column " << curLoc.col << ":\n"; + ss << '\t' << msg << '\n'; + errorMsg += ss.str(); + return {}; + } - template T logError(std::string_view msg) { - errorMsg += "GMSH parser error:\n\t"; - errorMsg += msg; - errorMsg += '\n'; - return {}; - } + template T logError(std::string_view msg) { + errorMsg += "GMSH parser error:\n\t"; + errorMsg += msg; + errorMsg += '\n'; + return {}; + } - GMSHMeshBuilder* builder; - GMSHSourceLocation curLoc = {0, 0}; - GMSHToken curTok; - GMSHLexer lexer; - std::string errorMsg; + GMSHMeshBuilder* builder; + GMSHSourceLocation curLoc = {0, 0}; + GMSHToken curTok; + GMSHLexer lexer; + std::string errorMsg; - GMSHToken getNextToken() { - curLoc = lexer.getSourceLoc(); - return curTok = lexer.getToken(); - } + GMSHToken getNextToken() { + curLoc = lexer.getSourceLoc(); + return curTok = lexer.getToken(); + } - std::optional getNumber() { - if (curTok == GMSHToken::integer) { - return {lexer.getInteger()}; - } else if (curTok == GMSHToken::real) { - return {lexer.getReal()}; - } - return std::nullopt; + std::optional getNumber() { + if (curTok == GMSHToken::integer) { + return {lexer.getInteger()}; + } else if (curTok == GMSHToken::real) { + return {lexer.getReal()}; } + return std::nullopt; + } - virtual bool parse_() = 0; + virtual bool parse_() = 0; - double parseMeshFormat(); + double parseMeshFormat(); }; } // namespace tndm diff --git a/src/third_party/README.md b/src/third_party/README.md index 6b1e238..3db1044 100644 --- a/src/third_party/README.md +++ b/src/third_party/README.md @@ -1,3 +1,4 @@ # Third party sources -Taken from https://github.com/TEAR-ERC/tandem/commit/555bbfa06aad688d609cfa2fae6b5bc2e2315589 +Taken from ; +one further addition from .