From f76948b12be52a8baa74dac09067adfb451df0ae Mon Sep 17 00:00:00 2001 From: David Schneller Date: Mon, 24 Jul 2023 11:12:57 +0200 Subject: [PATCH 1/2] Add PUMLcube --- pumlcube/CMakeLists.txt | 55 ++++++ pumlcube/src/main.cpp | 406 ++++++++++++++++++++++++++++++++++++++++ pumlcube/submodules | 1 + 3 files changed, 462 insertions(+) create mode 100644 pumlcube/CMakeLists.txt create mode 100644 pumlcube/src/main.cpp create mode 120000 pumlcube/submodules diff --git a/pumlcube/CMakeLists.txt b/pumlcube/CMakeLists.txt new file mode 100644 index 0000000..994edf8 --- /dev/null +++ b/pumlcube/CMakeLists.txt @@ -0,0 +1,55 @@ +cmake_minimum_required(VERSION 3.10) + +# set the project name +project(pumlcube) + +# specify the C++ standard +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED True) + +cmake_policy(SET CMP0074 NEW) + + +add_executable(pumlcube src/main.cpp) + +#logging +set(LOG_LEVEL "info" CACHE STRING "Log level for the code") +set(LOG_LEVEL_OPTIONS "debug" "info" "warning" "error") +set_property(CACHE LOG_LEVEL PROPERTY STRINGS ${LOG_LEVEL_OPTIONS}) +if("${LOG_LEVEL}" STREQUAL "debug") + target_compile_definitions(pumlcube PUBLIC LOG_LEVEL=3) +elseif("${LOG_LEVEL}" STREQUAL "info") + target_compile_definitions(pumlcube PUBLIC LOG_LEVEL=2) +elseif("${LOG_LEVEL}" STREQUAL "warning") + target_compile_definitions(pumlcube PUBLIC LOG_LEVEL=1) +elseif("${LOG_LEVEL}" STREQUAL "error") + target_compile_definitions(pumlcube PUBLIC LOG_LEVEL=0) +endif() + +#build and link libraries and executable +set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) + +find_package(HDF5 REQUIRED + COMPONENTS C HL) +target_include_directories(pumlcube PUBLIC ${HDF5_INCLUDE_DIRS}) +target_link_libraries(pumlcube PUBLIC ${HDF5_C_HL_LIBRARIES} ${HDF5_C_LIBRARIES}) + +find_package(MPI REQUIRED) +target_link_libraries(pumlcube PUBLIC MPI::MPI_CXX) + +find_package(OpenMP REQUIRED) +target_link_libraries(pumlcube PUBLIC OpenMP::OpenMP_CXX) + +#add some compiler specific flags +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + # using GCC + target_compile_options(pumlcube PUBLIC -fopenmp -pedantic $<$,$>:-Wall -Wextra -Wno-unused-parameter -Wno-unknown-pragmas>) + target_link_libraries(pumlcube PUBLIC "-fopenmp") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + target_compile_options(pumlcube PUBLIC -qopenmp -pedantic $<$,$>:-Wall -w3 -diag-disable:remark>) + target_link_libraries(pumlcube PUBLIC "-qopenmp") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp=libomp -Wall -Wextra -pedantic") +endif() + +target_include_directories(pumlcube PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/submodules) diff --git a/pumlcube/src/main.cpp b/pumlcube/src/main.cpp new file mode 100644 index 0000000..5c14355 --- /dev/null +++ b/pumlcube/src/main.cpp @@ -0,0 +1,406 @@ +/* (c) 2023 SeisSol, David Schneller. BSD-3 License */ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "utils/args.h" +#include "utils/logger.h" + +// cf. https://en.wikipedia.org/wiki/Schl%C3%A4fli_orthoscheme +static std::array, 4>, 6> cube2tet = { + std::array, 4>{ + std::array{0,0,0}, + std::array{1,0,0}, + std::array{0,1,0}, + std::array{0,1,1}, + }, + std::array, 4>{ + std::array{0,0,0}, + std::array{1,0,0}, + std::array{0,1,1}, + std::array{0,0,1}, + }, + std::array, 4>{ + std::array{1,0,0}, + std::array{0,1,1}, + std::array{0,0,1}, + std::array{1,0,1}, + }, + std::array, 4>{ + std::array{1,1,1}, + std::array{0,1,1}, + std::array{1,0,1}, + std::array{1,0,0}, + }, + std::array, 4>{ + std::array{1,1,1}, + std::array{0,1,1}, + std::array{1,0,0}, + std::array{1,1,0}, + }, + std::array, 4>{ + std::array{0,1,1}, + std::array{1,0,0}, + std::array{1,1,0}, + std::array{0,1,0}, + }, +}; + +// I don't know why this order works. But it seems like to does. +static std::array faceorder = { 3,2,0,1 }; + +static std::array, 4>, 6> tet2face(const std::array, 4>, 6>& data) { + std::array, 4>, 6> output; + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 4; ++j) { + for (int c = 0; c < 3; ++c) { + std::array values; + int d = 0; + for (int k = 0; k < 4; ++k) { + if (faceorder[j] != k) { + values[d] = data[i][k][c]; + ++d; + } + } + if (values[0] == values[1] && values[1] == values[2]) { + output[i][j][c] = values[0]; + } + else { + output[i][j][c] = -1; + } + } + } + } + return output; +} + +static std::array, 4>, 6> cube2tetface = tet2face(cube2tet); + +static std::array dim2str{"x", "y", "z"}; + +struct CubeDimension { + bool periodic; + int minBoundary, maxBoundary; + std::size_t numCells; + double scaling; + double translation; + + // (NOTE: utils::Args does not support const& right now for what we want to do) + CubeDimension(utils::Args& args, int dimension) { + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + logInfo(rank) << "Reading dimension" << dim2str[dimension]; + + periodic = args.getArgument("bper" + dim2str[dimension], true); + if (periodic) { + // periodic boundary values + minBoundary = 6; + maxBoundary = 6; + } + else { + minBoundary = args.getArgument("bmin" + dim2str[dimension], 1); + maxBoundary = args.getArgument("bmax" + dim2str[dimension], 1); + } + + numCells = args.getArgument("n" + dim2str[dimension]); + + scaling = args.getArgument("s" + dim2str[dimension], 1.0); + translation = args.getArgument("t" + dim2str[dimension], 0.0); + } + + std::size_t cellCount() const { + return numCells; + } + + std::size_t faceCount() const { + return cellCount() + 1; + } + + int getFaceType(std::size_t id) const { + if (id == 0) { + // left boundary + return minBoundary; + } + else if (id == cellCount()) { + // right boundary + return maxBoundary; + } + else { + // interior faces + return 0; + } + } + + double facePosition(std::size_t id) const { + return (scaling / cellCount()) * id + translation; + } + + std::size_t getFaceIdentification(std::size_t id) const { + // for periodicity, identify with the opposite-side vertex + if (periodic && id == cellCount()) { + return 0; + } + + // otherwise, identify with itself + return id; + } +}; + +struct Cubes { + CubeDimension dimx; + CubeDimension dimy; + CubeDimension dimz; + + Cubes(utils::Args& args) + : dimx(args, 0), dimy(args, 1), dimz(args, 2) + { + + } + + std::size_t vertexCount() const { + return dimx.faceCount() * dimy.faceCount() * dimz.faceCount(); + } + + std::size_t cubeCount() const { + return dimx.cellCount() * dimy.cellCount() * dimz.cellCount(); + } + + std::array unpackCubeId(std::size_t id) const { + auto z = id / (dimy.cellCount() * dimx.cellCount()); + auto xy = id % (dimy.cellCount() * dimx.cellCount()); + auto y = xy / dimx.cellCount(); + auto x = xy % dimx.cellCount(); + return {x,y,z}; + } + + std::array unpackVertexId(std::size_t id) const { + auto z = id / (dimy.faceCount() * dimx.faceCount()); + auto xy = id % (dimy.faceCount() * dimx.faceCount()); + auto y = xy / dimx.faceCount(); + auto x = xy % dimx.faceCount(); + return {x,y,z}; + } + + std::size_t packCubeId(const std::array& id) const { + return dimx.cellCount() * dimy.cellCount() * id[2] + dimx.cellCount() * id[1] + id[0]; + } + + std::size_t packVertexId(const std::array& id) const { + return dimx.faceCount() * dimy.faceCount() * id[2] + dimx.faceCount() * id[1] + id[0]; + } + + std::array, 6> cubeVertices(std::size_t id) const { + auto cid = unpackCubeId(id); + std::array, 6> indices; + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 4; ++j) { + auto offset = cube2tet[i][j]; + std::array vid = {cid[0] + offset[0], cid[1] + offset[1], cid[2] + offset[2]}; + indices[i][j] = packVertexId(vid); + } + } + return indices; + } + + std::array cubeBoundary(std::size_t id) const { + auto cid = unpackCubeId(id); + std::array indices; + for (int i = 0; i < 6; ++i) { + indices[i] = 0; + for (int j = 0; j < 4; ++j) { + auto offset = cube2tetface[i][j]; + int type = 0; + if (offset[0] >= 0) { + type = dimx.getFaceType(cid[0] + offset[0]); + } + if (offset[1] >= 0) { + type = dimy.getFaceType(cid[1] + offset[1]); + } + if (offset[2] >= 0) { + type = dimz.getFaceType(cid[2] + offset[2]); + } + indices[i] |= ((type & 0xff) << (8*j)); + } + } + return indices; + } + + std::array vertex(std::size_t id) const { + auto vid = unpackVertexId(id); + return {dimx.facePosition(vid[0]), dimy.facePosition(vid[1]), dimz.facePosition(vid[2])}; + } + + std::size_t vertexIdentification(std::size_t id) const { + auto vid = unpackVertexId(id); + return packVertexId({dimx.getFaceIdentification(vid[0]), dimy.getFaceIdentification(vid[1]), dimz.getFaceIdentification(vid[2])}); + } +}; + +template +static void writeData(hid_t fileHandle, hid_t xfer, const std::string& name, std::size_t count, std::size_t chunksize, const std::vector& elemsize, hid_t datatype, F&& accessor) { + int commSize; + int commRank; + MPI_Comm_size(MPI_COMM_WORLD, &commSize); + MPI_Comm_rank(MPI_COMM_WORLD, &commRank); + + logInfo(commRank) << "Writing" << name; + + std::size_t localCount = count / commSize; + std::size_t localOffset = localCount * commRank; + auto rest = count % commSize; + if ((std::size_t)commRank < rest) { + ++localCount; + localOffset += commRank; + } + else { + localOffset += rest; + } + + std::vector globalDims(elemsize.size() + 1); + globalDims[0] = chunksize * count; + for (std::size_t i = 0; i < elemsize.size(); ++i) { + globalDims[i + 1] = elemsize[i]; + } + + std::size_t localFlattened = 0; + std::size_t localChunk = 0; + std::vector localDims(elemsize.size() + 1); + localDims[0] = chunksize * localCount; + localFlattened = localDims[0]; + localChunk = chunksize; + for (std::size_t i = 0; i < elemsize.size(); ++i) { + localDims[i + 1] = elemsize[i]; + localFlattened *= elemsize[i]; + localChunk *= elemsize[i]; + } + + std::vector localStart(elemsize.size() + 1); + localStart[0] = chunksize * localOffset; + for (std::size_t i = 0; i < elemsize.size(); ++i) { + localStart[i + 1] = 0; + } + + hid_t h5p = H5Pcreate(H5P_DATASET_CREATE); + H5Pset_layout(h5p, H5D_CONTIGUOUS); + H5Pset_alloc_time(h5p, H5D_ALLOC_TIME_EARLY); + + hid_t globalSpace = H5Screate_simple(globalDims.size(), globalDims.data(), nullptr); + hid_t localSpace = H5Screate_simple(localDims.size(), localDims.data(), nullptr); + hid_t datasetHandle = H5Dcreate(fileHandle, name.c_str(), datatype, globalSpace, H5P_DEFAULT, h5p, H5P_DEFAULT); + + H5Pclose(h5p); + + H5Sselect_hyperslab(globalSpace, H5S_SELECT_SET, localStart.data(), nullptr, localDims.data(), nullptr); + H5Sselect_all(localSpace); + + std::vector buffer(localFlattened); + + #pragma omp parallel for + for (std::size_t i = 0; i < localCount; ++i) { + std::invoke(accessor, i + localOffset, i * localChunk, buffer); + } + + H5Dwrite(datasetHandle, datatype, localSpace, globalSpace, xfer, buffer.data()); + + H5Sclose(localSpace); + H5Sclose(globalSpace); + + H5Dclose(datasetHandle); +} + +int main(int argc, char** argv) { + int provided; + MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided); + + int commSize; + int commRank; + MPI_Comm_size(MPI_COMM_WORLD, &commSize); + MPI_Comm_rank(MPI_COMM_WORLD, &commRank); + + utils::Args args; + args.addOption("bperx", 0, "periodicity in x dimension", utils::Args::Required, false); + args.addOption("bpery", 0, "periodicity in y dimension", utils::Args::Required, false); + args.addOption("bperz", 0, "periodicity in z dimension", utils::Args::Required, false); + args.addOption("bminx", 0, "boundary condition at the start of x dimension (ignored, if bperx is true)", utils::Args::Required, false); + args.addOption("bmaxx", 0, "boundary condition at the end of x dimension (ignored, if bperx is true)", utils::Args::Required, false); + args.addOption("bminy", 0, "boundary condition at the start of y dimension (ignored, if bpery is true)", utils::Args::Required, false); + args.addOption("bmaxy", 0, "boundary condition at the end of y dimension (ignored, if bpery is true)", utils::Args::Required, false); + args.addOption("bminz", 0, "boundary condition at the start of z dimension (ignored, if bperz is true)", utils::Args::Required, false); + args.addOption("bmaxz", 0, "boundary condition at the end of z dimension (ignored, if bperz is true)", utils::Args::Required, false); + args.addOption("nx", 'x', "number of cubes in x dimension", utils::Args::Required, true); + args.addOption("ny", 'y', "number of cubes in y dimension", utils::Args::Required, true); + args.addOption("nz", 'z', "number of cubes in z dimension", utils::Args::Required, true); + args.addOption("output", 'o', "output file for resulting PUML mesh", utils::Args::Required, true); + args.addOption("sx", 0, "scale in x direction", utils::Args::Required, false); + args.addOption("sy", 0, "scale in y direction", utils::Args::Required, false); + args.addOption("sz", 0, "scale in z direction", utils::Args::Required, false); + args.addOption("tx", 0, "mesh translation in x direction (after scaling)", utils::Args::Required, false); + args.addOption("ty", 0, "mesh translation in y direction (after scaling)", utils::Args::Required, false); + args.addOption("tz", 0, "mesh translation in z direction (after scaling)", utils::Args::Required, false); + + logInfo(commRank) << "Using" << omp_get_max_threads() << "threads, and " << commSize << "ranks"; + + if (args.parse(argc, argv)) { + logError() << "Error parsing arguments"; + } + + Cubes cubes(args); + + auto filename = args.getArgument("output"); + + logInfo(commRank) << "Output to file" << filename; + + hid_t h5plist = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_libver_bounds(h5plist, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST); + H5Pset_fapl_mpio(h5plist, MPI_COMM_WORLD, MPI_INFO_NULL); + + hid_t fileHandle = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, h5plist); + H5Pclose(h5plist); + + hid_t xfer = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(xfer, H5FD_MPIO_COLLECTIVE); + + writeData(fileHandle, xfer, "geometry", cubes.vertexCount(), 1, {3}, H5T_NATIVE_DOUBLE, [&](std::size_t id, std::size_t datapos, std::vector& target) { + auto data = cubes.vertex(id); + target[datapos+0] = data[0]; + target[datapos+1] = data[1]; + target[datapos+2] = data[2]; + }); + writeData(fileHandle, xfer, "connect", cubes.cubeCount(), 6, {4}, H5T_NATIVE_ULONG, [&](std::size_t id, std::size_t datapos, std::vector& target){ + auto data = cubes.cubeVertices(id); + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 4; ++j) { + target[datapos + 4*i + j] = data[i][j]; + } + } + }); + writeData(fileHandle, xfer, "group", cubes.cubeCount(), 6, {}, H5T_NATIVE_INT, [&](std::size_t id, std::size_t datapos, std::vector& target){ + for (int i = 0; i < 6; ++i) { + target[datapos + i] = 0; + } + }); + writeData(fileHandle, xfer, "boundary", cubes.cubeCount(), 6, {}, H5T_NATIVE_UINT, [&](std::size_t id, std::size_t datapos, std::vector& target){ + auto data = cubes.cubeBoundary(id); + for (int i = 0; i < 6; ++i) { + target[datapos + i] = data[i]; + } + }); + writeData(fileHandle, xfer, "identify", cubes.vertexCount(), 1, {}, H5T_NATIVE_ULONG, [&](std::size_t id, std::size_t datapos, std::vector& target){ + target[datapos] = cubes.vertexIdentification(id); + }); + + H5Pclose(xfer); + H5Fclose(fileHandle); + + MPI_Finalize(); + return 0; +} diff --git a/pumlcube/submodules b/pumlcube/submodules new file mode 120000 index 0000000..1bc410d --- /dev/null +++ b/pumlcube/submodules @@ -0,0 +1 @@ +../submodules \ No newline at end of file From 542e5f6fb3669ab406a036d5f417c6225b733614 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 6 Nov 2024 20:29:08 +0100 Subject: [PATCH 2/2] Address review * add a slim Xdmf and XML writer * add HDF5 error checking functions * update boundary format info --- pumlcube/src/main.cpp | 409 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 336 insertions(+), 73 deletions(-) diff --git a/pumlcube/src/main.cpp b/pumlcube/src/main.cpp index 5c14355..1910e37 100644 --- a/pumlcube/src/main.cpp +++ b/pumlcube/src/main.cpp @@ -1,20 +1,41 @@ -/* (c) 2023 SeisSol, David Schneller. BSD-3 License */ +// SPDX-FileCopyrightText: 2023-2024 SeisSol Group +// +// SPDX-License-Identifier: BSD-3-Clause -#include -#include +#include +#include #include #include #include #include +#include #include #include #include #include +#include +#include +#include #include "utils/args.h" #include "utils/logger.h" +template static TT _hw(TT&& status, const char* file, int line) { + if (status < 0) { + logError() << utils::nospace << "An HDF5 error occurred (" << file << ": " << line << ")"; + } + return std::forward(status); +} + +#define hw(...) _hw(__VA_ARGS__, __FILE__, __LINE__) + +enum class BoundaryFormat { + Int32, + Int64, + Int32x4 +}; + // cf. https://en.wikipedia.org/wiki/Schl%C3%A4fli_orthoscheme static std::array, 4>, 6> cube2tet = { std::array, 4>{ @@ -100,15 +121,12 @@ struct CubeDimension { MPI_Comm_rank(MPI_COMM_WORLD, &rank); logInfo(rank) << "Reading dimension" << dim2str[dimension]; - periodic = args.getArgument("bper" + dim2str[dimension], true); - if (periodic) { - // periodic boundary values - minBoundary = 6; - maxBoundary = 6; - } - else { - minBoundary = args.getArgument("bmin" + dim2str[dimension], 1); - maxBoundary = args.getArgument("bmax" + dim2str[dimension], 1); + minBoundary = args.getArgument("bmin" + dim2str[dimension], 1); + maxBoundary = args.getArgument("bmax" + dim2str[dimension], 1); + + periodic = minBoundary == 6 && maxBoundary == 6; + if ((minBoundary == 6 || maxBoundary == 6) && !periodic) { + logError() << "One-sided periodic boundary conditions (i.e. value 6 for bmin or bmax) are not supported."; } numCells = args.getArgument("n" + dim2str[dimension]); @@ -174,6 +192,10 @@ struct Cubes { return dimx.cellCount() * dimy.cellCount() * dimz.cellCount(); } + std::size_t tetrahedronCount() const { + return cubeCount() * 6; + } + std::array unpackCubeId(std::size_t id) const { auto z = id / (dimy.cellCount() * dimx.cellCount()); auto xy = id % (dimy.cellCount() * dimx.cellCount()); @@ -211,24 +233,42 @@ struct Cubes { return indices; } - std::array cubeBoundary(std::size_t id) const { - auto cid = unpackCubeId(id); - std::array indices; + int faceBoundary(std::size_t id, int i, int j) const { + const auto cid = unpackCubeId(id); + const auto offset = cube2tetface[i][j]; + int type = 0; + if (offset[0] >= 0) { + type = dimx.getFaceType(cid[0] + offset[0]); + } + if (offset[1] >= 0) { + type = dimy.getFaceType(cid[1] + offset[1]); + } + if (offset[2] >= 0) { + type = dimz.getFaceType(cid[2] + offset[2]); + } + return type; + } + + template + std::array cubeBoundary(std::size_t id, int typesize) const { + std::array indices; + const T mask = (T(1) << typesize) - 1; for (int i = 0; i < 6; ++i) { indices[i] = 0; for (int j = 0; j < 4; ++j) { - auto offset = cube2tetface[i][j]; - int type = 0; - if (offset[0] >= 0) { - type = dimx.getFaceType(cid[0] + offset[0]); - } - if (offset[1] >= 0) { - type = dimy.getFaceType(cid[1] + offset[1]); - } - if (offset[2] >= 0) { - type = dimz.getFaceType(cid[2] + offset[2]); - } - indices[i] |= ((type & 0xff) << (8*j)); + const auto type = faceBoundary(id, i, j); + indices[i] |= ((type & mask) << (typesize*j)); + } + } + return indices; + } + + std::array cubeBoundary4x(std::size_t id) const { + std::array indices; + for (int i = 0; i < 6; ++i) { + for (int j = 0; j < 4; ++j) { + const auto type = faceBoundary(id, i, j); + indices[i * 4 + j] = type; } } return indices; @@ -245,6 +285,158 @@ struct Cubes { } }; +class XmlWriter { +public: + XmlWriter(const std::string& file) : stream(file) { + + } + + void writeIndent() { + for (std::size_t i = 0; i < nodes.size() * indent; ++i) { + stream << ' '; + } + } + + void writeHeader() { + stream << "\n"; + } + + void writeDoctype() { + stream << "\n"; + } + + void writeNode(const std::string& name) { + writeNodeHeaderEnd(); + writeIndent(); + stream << "<" << name; + nodes.push(name); + inHeader = true; + } + + void writeNodeEnd() { + const auto node = nodes.top(); + nodes.pop(); + if (inHeader) { + stream << " />\n"; + } + else { + writeIndent(); + stream << "\n"; + } + } + + void writeNodeHeaderEnd() { + if (inHeader) { + stream << ">\n"; + inHeader = false; + } + } + + template + void writeNodeAttribute(const std::string& name, const T& value) { + assert(inHeader); + stream << " " << name << "=\"" << value << "\""; + } + + template + void writeData(const T& data) { + writeNodeHeaderEnd(); + writeIndent(); + stream << data << "\n"; + } + + void closeDocument() { + while (!nodes.empty()) { + writeNodeEnd(); + } + stream.flush(); + } + + void writeComment(const std::string& text) { + writeNodeHeaderEnd(); + writeIndent(); + stream << "\n"; + } +private: + std::ofstream stream; + std::stack nodes; + bool inHeader{false}; + int indent{2}; +}; + +class XdmfXml { +public: + XdmfXml(const std::string& fileName, const std::string& dataFile) : xml(fileName), dataFile(dataFile) { + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + localWrite = rank == 0; + + logInfo(rank) << "Writing Xdmf file to" << fileName; + } + + void begin() { + if (localWrite) { + xml.writeHeader(); + xml.writeDoctype(); + xml.writeNode("Xdmf"); + xml.writeNodeAttribute("Version", "2.0"); + xml.writeNode("Domain"); + xml.writeNode("Grid"); + xml.writeNodeAttribute("Name", "cubemesh"); + xml.writeNodeAttribute("GridType", "Uniform"); + } + } + + void addData(const std::string& dataset, const std::vector& dimensions, const std::string& type, const std::string& name, const std::string& datatype, std::size_t size) { + if (localWrite) { + std::ostringstream dimensionStream; + dimensionStream << dimensions[0]; + for (std::size_t i = 1; i < dimensions.size(); ++i) { + dimensionStream << " " << dimensions[i]; + } + if (type == "Topology") { + xml.writeNode("Topology"); + xml.writeNodeAttribute("TopologyType", "Tetrahedron"); + xml.writeNodeAttribute("NumberOfElements", dimensions[0]); + } + else if (type == "Geometry") { + xml.writeNode("Geometry"); + xml.writeNodeAttribute("GeometryType", "XYZ"); + xml.writeNodeAttribute("NumberOfElements", dimensions[0]); + } + else if (type == "AttributeCell") { + xml.writeNode("Attribute"); + xml.writeNodeAttribute("Name", name); + xml.writeNodeAttribute("Center", "Cell"); + } + else if (type == "AttributeNode") { + xml.writeNode("Attribute"); + xml.writeNodeAttribute("Name", name); + xml.writeNodeAttribute("Center", "Node"); + } + + xml.writeNode("DataItem"); + xml.writeNodeAttribute("NumberType", datatype); + xml.writeNodeAttribute("Precision", size); + xml.writeNodeAttribute("Format", "HDF"); + xml.writeNodeAttribute("Dimensions", dimensionStream.str()); + xml.writeData(dataFile + ":/" + dataset); + xml.writeNodeEnd(); + xml.writeNodeEnd(); + } + } + + void end() { + if (localWrite) { + xml.closeDocument(); + } + } +private: + XmlWriter xml; + bool localWrite; + std::string dataFile; +}; + template static void writeData(hid_t fileHandle, hid_t xfer, const std::string& name, std::size_t count, std::size_t chunksize, const std::vector& elemsize, hid_t datatype, F&& accessor) { int commSize; @@ -289,18 +481,18 @@ static void writeData(hid_t fileHandle, hid_t xfer, const std::string& name, std localStart[i + 1] = 0; } - hid_t h5p = H5Pcreate(H5P_DATASET_CREATE); - H5Pset_layout(h5p, H5D_CONTIGUOUS); - H5Pset_alloc_time(h5p, H5D_ALLOC_TIME_EARLY); + hid_t h5p = hw(H5Pcreate(H5P_DATASET_CREATE)); + hw(H5Pset_layout(h5p, H5D_CONTIGUOUS)); + hw(H5Pset_alloc_time(h5p, H5D_ALLOC_TIME_EARLY)); - hid_t globalSpace = H5Screate_simple(globalDims.size(), globalDims.data(), nullptr); - hid_t localSpace = H5Screate_simple(localDims.size(), localDims.data(), nullptr); - hid_t datasetHandle = H5Dcreate(fileHandle, name.c_str(), datatype, globalSpace, H5P_DEFAULT, h5p, H5P_DEFAULT); + hid_t globalSpace = hw(H5Screate_simple(globalDims.size(), globalDims.data(), nullptr)); + hid_t localSpace = hw(H5Screate_simple(localDims.size(), localDims.data(), nullptr)); + hid_t datasetHandle = hw(H5Dcreate(fileHandle, name.c_str(), datatype, globalSpace, H5P_DEFAULT, h5p, H5P_DEFAULT)); - H5Pclose(h5p); + hw(H5Pclose(h5p)); - H5Sselect_hyperslab(globalSpace, H5S_SELECT_SET, localStart.data(), nullptr, localDims.data(), nullptr); - H5Sselect_all(localSpace); + hw(H5Sselect_hyperslab(globalSpace, H5S_SELECT_SET, localStart.data(), nullptr, localDims.data(), nullptr)); + hw(H5Sselect_all(localSpace)); std::vector buffer(localFlattened); @@ -309,12 +501,12 @@ static void writeData(hid_t fileHandle, hid_t xfer, const std::string& name, std std::invoke(accessor, i + localOffset, i * localChunk, buffer); } - H5Dwrite(datasetHandle, datatype, localSpace, globalSpace, xfer, buffer.data()); + hw(H5Dwrite(datasetHandle, datatype, localSpace, globalSpace, xfer, buffer.data())); - H5Sclose(localSpace); - H5Sclose(globalSpace); + hw(H5Sclose(localSpace)); + hw(H5Sclose(globalSpace)); - H5Dclose(datasetHandle); + hw(H5Dclose(datasetHandle)); } int main(int argc, char** argv) { @@ -327,27 +519,26 @@ int main(int argc, char** argv) { MPI_Comm_rank(MPI_COMM_WORLD, &commRank); utils::Args args; - args.addOption("bperx", 0, "periodicity in x dimension", utils::Args::Required, false); - args.addOption("bpery", 0, "periodicity in y dimension", utils::Args::Required, false); - args.addOption("bperz", 0, "periodicity in z dimension", utils::Args::Required, false); - args.addOption("bminx", 0, "boundary condition at the start of x dimension (ignored, if bperx is true)", utils::Args::Required, false); - args.addOption("bmaxx", 0, "boundary condition at the end of x dimension (ignored, if bperx is true)", utils::Args::Required, false); - args.addOption("bminy", 0, "boundary condition at the start of y dimension (ignored, if bpery is true)", utils::Args::Required, false); - args.addOption("bmaxy", 0, "boundary condition at the end of y dimension (ignored, if bpery is true)", utils::Args::Required, false); - args.addOption("bminz", 0, "boundary condition at the start of z dimension (ignored, if bperz is true)", utils::Args::Required, false); - args.addOption("bmaxz", 0, "boundary condition at the end of z dimension (ignored, if bperz is true)", utils::Args::Required, false); - args.addOption("nx", 'x', "number of cubes in x dimension", utils::Args::Required, true); - args.addOption("ny", 'y', "number of cubes in y dimension", utils::Args::Required, true); - args.addOption("nz", 'z', "number of cubes in z dimension", utils::Args::Required, true); - args.addOption("output", 'o', "output file for resulting PUML mesh", utils::Args::Required, true); + args.addOption("bminx", 0, "boundary condition at the start of x dimension", utils::Args::Required, false); + args.addOption("bmaxx", 0, "boundary condition at the end of x dimension", utils::Args::Required, false); + args.addOption("bminy", 0, "boundary condition at the start of y dimension", utils::Args::Required, false); + args.addOption("bmaxy", 0, "boundary condition at the end of y dimension", utils::Args::Required, false); + args.addOption("bminz", 0, "boundary condition at the start of z dimension", utils::Args::Required, false); + args.addOption("bmaxz", 0, "boundary condition at the end of z dimension", utils::Args::Required, false); + args.addOption("nx", 'x', "number of cubes in x dimension", utils::Args::Required, true); + args.addOption("ny", 'y', "number of cubes in y dimension", utils::Args::Required, true); + args.addOption("nz", 'z', "number of cubes in z dimension", utils::Args::Required, true); + args.addOption("output", 'o', "output file for resulting PUML mesh", utils::Args::Required, true); args.addOption("sx", 0, "scale in x direction", utils::Args::Required, false); args.addOption("sy", 0, "scale in y direction", utils::Args::Required, false); args.addOption("sz", 0, "scale in z direction", utils::Args::Required, false); args.addOption("tx", 0, "mesh translation in x direction (after scaling)", utils::Args::Required, false); args.addOption("ty", 0, "mesh translation in y direction (after scaling)", utils::Args::Required, false); args.addOption("tz", 0, "mesh translation in z direction (after scaling)", utils::Args::Required, false); + args.addOption("xdmf", 0, "create an XDMF XML file", utils::Args::Required, false); + args.addOption("boundary", 0, "boundary format to use", utils::Args::Required, false); - logInfo(commRank) << "Using" << omp_get_max_threads() << "threads, and " << commSize << "ranks"; + logInfo(commRank) << "Using" << omp_get_max_threads() << "threads, and " << commSize << "ranks"; if (args.parse(argc, argv)) { logError() << "Error parsing arguments"; @@ -355,19 +546,37 @@ int main(int argc, char** argv) { Cubes cubes(args); - auto filename = args.getArgument("output"); + const auto filename = args.getArgument("output"); + + const auto boundaryRaw = args.getArgument("boundary", "i32"); + + BoundaryFormat boundaryFormat = BoundaryFormat::Int32; + if (boundaryRaw == "i32") { + boundaryFormat = BoundaryFormat::Int32; + } + else if (boundaryRaw == "i64") { + boundaryFormat = BoundaryFormat::Int64; + } + else if (boundaryRaw == "i32x4") { + boundaryFormat = BoundaryFormat::Int32x4; + } + else { + logError() << "Unknown boundary format:" << boundaryRaw; + } + + const bool periodic = cubes.dimx.periodic || cubes.dimy.periodic || cubes.dimz.periodic; logInfo(commRank) << "Output to file" << filename; - hid_t h5plist = H5Pcreate(H5P_FILE_ACCESS); - H5Pset_libver_bounds(h5plist, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST); - H5Pset_fapl_mpio(h5plist, MPI_COMM_WORLD, MPI_INFO_NULL); + hid_t h5plist = hw(H5Pcreate(H5P_FILE_ACCESS)); + hw(H5Pset_libver_bounds(h5plist, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST)); + hw(H5Pset_fapl_mpio(h5plist, MPI_COMM_WORLD, MPI_INFO_NULL)); - hid_t fileHandle = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, h5plist); - H5Pclose(h5plist); + hid_t fileHandle = hw(H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, h5plist)); + hw(H5Pclose(h5plist)); - hid_t xfer = H5Pcreate(H5P_DATASET_XFER); - H5Pset_dxpl_mpio(xfer, H5FD_MPIO_COLLECTIVE); + hid_t xfer = hw(H5Pcreate(H5P_DATASET_XFER)); + hw(H5Pset_dxpl_mpio(xfer, H5FD_MPIO_COLLECTIVE)); writeData(fileHandle, xfer, "geometry", cubes.vertexCount(), 1, {3}, H5T_NATIVE_DOUBLE, [&](std::size_t id, std::size_t datapos, std::vector& target) { auto data = cubes.vertex(id); @@ -375,7 +584,7 @@ int main(int argc, char** argv) { target[datapos+1] = data[1]; target[datapos+2] = data[2]; }); - writeData(fileHandle, xfer, "connect", cubes.cubeCount(), 6, {4}, H5T_NATIVE_ULONG, [&](std::size_t id, std::size_t datapos, std::vector& target){ + writeData(fileHandle, xfer, "connect", cubes.cubeCount(), 6, {4}, H5T_NATIVE_UINT64, [&](std::size_t id, std::size_t datapos, std::vector& target){ auto data = cubes.cubeVertices(id); for (int i = 0; i < 6; ++i) { for (int j = 0; j < 4; ++j) { @@ -388,18 +597,72 @@ int main(int argc, char** argv) { target[datapos + i] = 0; } }); - writeData(fileHandle, xfer, "boundary", cubes.cubeCount(), 6, {}, H5T_NATIVE_UINT, [&](std::size_t id, std::size_t datapos, std::vector& target){ - auto data = cubes.cubeBoundary(id); - for (int i = 0; i < 6; ++i) { - target[datapos + i] = data[i]; + if (boundaryFormat == BoundaryFormat::Int32) { + writeData(fileHandle, xfer, "boundary", cubes.cubeCount(), 6, {}, H5T_NATIVE_INT32, [&](std::size_t id, std::size_t datapos, std::vector& target){ + auto data = cubes.cubeBoundary(id, 8); + for (int i = 0; i < 6; ++i) { + target[datapos + i] = data[i]; + } + }); + } + if (boundaryFormat == BoundaryFormat::Int64) { + writeData(fileHandle, xfer, "boundary", cubes.cubeCount(), 6, {}, H5T_NATIVE_INT64, [&](std::size_t id, std::size_t datapos, std::vector& target){ + auto data = cubes.cubeBoundary(id, 16); + for (int i = 0; i < 6; ++i) { + target[datapos + i] = data[i]; + } + }); + } + if (boundaryFormat == BoundaryFormat::Int32x4) { + writeData(fileHandle, xfer, "boundary", cubes.cubeCount(), 6, {4}, H5T_NATIVE_INT, [&](std::size_t id, std::size_t datapos, std::vector& target){ + auto data = cubes.cubeBoundary4x(id); + for (int i = 0; i < 24; ++i) { + target[datapos + i] = data[i]; + } + }); + } + if (periodic) { + writeData(fileHandle, xfer, "identify", cubes.vertexCount(), 1, {}, H5T_NATIVE_UINT64, [&](std::size_t id, std::size_t datapos, std::vector& target){ + target[datapos] = cubes.vertexIdentification(id); + }); + } + + hid_t attrSpace = hw(H5Screate(H5S_SCALAR)); + hid_t attrType = hw(H5Tcopy(H5T_C_S1)); + hw(H5Tset_size(attrType, H5T_VARIABLE)); + hid_t attrBoundary = hw( + H5Acreate(fileHandle, "boundary-format", attrType, attrSpace, H5P_DEFAULT, H5P_DEFAULT)); + const void* stringData = boundaryRaw.data(); + hw(H5Awrite(attrBoundary, attrType, &stringData)); + hw(H5Aclose(attrBoundary)); + hw(H5Sclose(attrSpace)); + hw(H5Tclose(attrType)); + + if (args.getArgument("xdmf", false)) { + auto xdmf = XdmfXml(filename + ".xdmf", filename); + xdmf.begin(); + xdmf.addData("geometry", {cubes.vertexCount(), 3}, "Geometry", "", "Float", 8); + xdmf.addData("connect", {cubes.tetrahedronCount(), 4}, "Topology", "", "Int", 8); + xdmf.addData("group", {cubes.tetrahedronCount()}, "AttributeCell", "group", "Int", 4); + if (boundaryFormat == BoundaryFormat::Int32) { + xdmf.addData("boundary", {cubes.tetrahedronCount()}, "AttributeCell", "boundary", "Int", 4); } - }); - writeData(fileHandle, xfer, "identify", cubes.vertexCount(), 1, {}, H5T_NATIVE_ULONG, [&](std::size_t id, std::size_t datapos, std::vector& target){ - target[datapos] = cubes.vertexIdentification(id); - }); + if (boundaryFormat == BoundaryFormat::Int64) { + xdmf.addData("boundary", {cubes.tetrahedronCount()}, "AttributeCell", "boundary", "Int", 8); + } + if (boundaryFormat == BoundaryFormat::Int32x4) { + xdmf.addData("boundary", {cubes.tetrahedronCount(), 4}, "AttributeCell", "boundary", "Int", 4); + } + if (periodic) { + xdmf.addData("identify", {cubes.vertexCount()}, "AttributeNode", "identify", "Int", 8); + } + xdmf.end(); + } + + hw(H5Pclose(xfer)); + hw(H5Fclose(fileHandle)); - H5Pclose(xfer); - H5Fclose(fileHandle); + logInfo(commRank) << "All done."; MPI_Finalize(); return 0;