diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index dfec676..a336d4d 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -87,6 +87,7 @@ jobs: - name: clang-format run: | sudo apt-get update - sudo apt-get install -qq python3 clang-format + sudo apt-get install -qq python3 python3-pip + pip3 install clang-format==18.1.5 git submodule update --init ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format --exclude src/third_party -r src diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index cb0066d..0000000 --- a/.travis.yml +++ /dev/null @@ -1,37 +0,0 @@ -sudo: required - -notifications: - email: - on_success: always - on_failure: always - -language: cpp - -compiler: gcc - -dist: bionic - -before_install: - - sudo add-apt-repository -y ppa:ubuntu-toolchain-r/test - - sudo apt-get update -qq - -install: - - sudo apt-get install -qq g++-5 openmpi-bin openmpi-common libopenmpi-dev hdf5-tools libhdf5-openmpi-100 libhdf5-openmpi-dev libnetcdf-dev clang-format-10 - - cd .. -# Download simmetrix libraries -# check license -#build pumi - - git clone https://github.com/SCOREC/core.git - - mkdir pumi - - mkdir core/build - - cd core/build - - cmake .. -DCMAKE_C_COMPILER="`which mpicc`" -DCMAKE_CXX_COMPILER="`which mpiCC`" -DCMAKE_C_FLAGS="-O2 -g -Wall" -DCMAKE_CXX_FLAGS="-O2 -g -Wall" -DCMAKE_INSTALL_PREFIX=../../pumi - - make -j2 && make install - - make -j2 && make install -#Go back to PUMGen - - cd ../../PUMGen - -script: - - ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format-10 --exclude src/third_party -r src - - mkdir build && cd build - - CC=mpicc CXX=mpiCC cmake .. -DCMAKE_PREFIX_PATH=../pumi && make -j2 diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e407e9..0efd978 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ set(SIMMETRIX_ROOT "" CACHE STRING "Root directory of simmetrix dev files") set(SIM_MPI "mpich3" CACHE STRING "MPI version used by simmetrix") option(NETCDF "Use netcdf" OFF) option(PARASOLID "enable support for parasolid files" OFF) +option(SCOREC "Use PUMI/APF" OFF) set(LOG_LEVEL "info" CACHE STRING "Log level for the code") set_property(CACHE LOG_LEVEL PROPERTY STRINGS "debug" "info" "warning" "error") @@ -24,20 +25,21 @@ set(CMAKE_CXX_EXTENSIONS OFF) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "RelWithDebInfo") # MinSizeRel is useless for us if(NOT CMAKE_BUILD_TYPE) - set(CMAKE_BUILD_TYPE Release) + set(CMAKE_BUILD_TYPE Release CACHE STRING "" FORCE) message(STATUS "Set build type to Release as none was supplied.") endif() set(CMAKE_CXX_FLAGS_RELEASE "-O2") add_executable(pumgen ${CMAKE_CURRENT_SOURCE_DIR}/src/pumgen.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/FidapReader.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GambitReader.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GMSH4Parser.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/input/ParallelVertexFilter.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSHLexer.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSHParser.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSH2Parser.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/aux/InsphereCalculator.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/aux/MPIConvenience.cpp ) target_include_directories(pumgen PUBLIC src @@ -74,9 +76,12 @@ target_link_libraries(pumgen PUBLIC MPI::MPI_CXX) target_compile_definitions(pumgen PUBLIC PARALLEL) -find_package(APF REQUIRED) -target_include_directories(pumgen PUBLIC ${APF_INCLUDE_DIR}) -target_link_libraries(pumgen PUBLIC ${APF_LIBRARIES}) +if (SCOREC) + find_package(APF REQUIRED) + target_include_directories(pumgen PUBLIC ${APF_INCLUDE_DIR}) + target_link_libraries(pumgen PUBLIC ${APF_LIBRARIES}) + target_compile_definitions(pumgen PUBLIC USE_SCOREC) +endif() add_library(tinyxml2 submodules/tinyxml2/tinyxml2.cpp) target_link_libraries(pumgen PUBLIC tinyxml2) diff --git a/azure-pipelines-build-template.yml b/azure-pipelines-build-template.yml deleted file mode 100644 index bf89ddb..0000000 --- a/azure-pipelines-build-template.yml +++ /dev/null @@ -1,26 +0,0 @@ -jobs: - - job: - displayName: pumgen-build - pool: - vmImage: 'ubuntu-22.04' - steps: - - bash: | - set -euo pipefail - export IFS=$'\n\t' - sudo apt-get update - sudo apt-get install -qq g++ openmpi-bin openmpi-common libopenmpi-dev hdf5-tools libhdf5-openmpi-103 libhdf5-openmpi-dev libnetcdf-dev - cd .. - # build pumi - git clone --recursive --branch v2.2.7 https://github.com/SCOREC/core.git - cd core - git --no-pager log --pretty=oneline -1 - mkdir build - cd build - cmake .. -DCMAKE_C_COMPILER="`which mpicc`" -DCMAKE_CXX_COMPILER="`which mpiCC`" -DCMAKE_C_FLAGS="-O2 -g -Wall" -DCMAKE_CXX_FLAGS="-O2 -g -Wall" -DSCOREC_CXX_WARNINGS=OFF -DCMAKE_INSTALL_PREFIX=../../pumi - make -j $(nproc) && make install - #Go back to PUMGen and build it - cd ../../s - git submodule update --init - mkdir build && cd build - CC=mpicc CXX=mpiCC cmake .. -DCMAKE_PREFIX_PATH=../pumi && make -j $(nproc) - diff --git a/azure-pipelines-clang-format-template.yml b/azure-pipelines-clang-format-template.yml deleted file mode 100644 index 9cb6996..0000000 --- a/azure-pipelines-clang-format-template.yml +++ /dev/null @@ -1,13 +0,0 @@ -jobs: - - job: - displayName: pumgen-clang-format - pool: - vmImage: 'ubuntu-22.04' - steps: - - bash: | - set -euo pipefail - export IFS=$'\n\t' - sudo apt-get update - sudo apt-get install -qq python3 clang-format-14 - git submodule update --init - ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format-14 --exclude src/third_party -r src diff --git a/azure-pipelines.yml b/azure-pipelines.yml deleted file mode 100644 index 997c5bf..0000000 --- a/azure-pipelines.yml +++ /dev/null @@ -1,3 +0,0 @@ -jobs: - - template: azure-pipelines-clang-format-template.yml - - template: azure-pipelines-build-template.yml diff --git a/cmake/FindSIMMETRIX.cmake b/cmake/FindSIMMETRIX.cmake index 55601a0..e8fdbef 100644 --- a/cmake/FindSIMMETRIX.cmake +++ b/cmake/FindSIMMETRIX.cmake @@ -1,18 +1,24 @@ include(FindPackageHandleStandardArgs) -find_path(GMI_SIM_INCLUDE_DIR gmi_sim.h) -find_path(APF_SIM_INCLUDE_DIR apfSIM.h) +if (SCOREC) + find_path(GMI_SIM_INCLUDE_DIR gmi_sim.h) + find_path(APF_SIM_INCLUDE_DIR apfSIM.h) + + list(APPEND SIMMETRIX_INCLUDE_DIR + ${GMI_SIM_INCLUDE_DIR} + ${APF_SIM_INCLUDE_DIR} + ) + + find_library(GMI_SIM_LIB gmi_sim) + find_library(APF_SIM_LIB apf_sim) +endif() + find_path(MESH_SIM_INCLUDE_DIR MeshSim.h) list(APPEND SIMMETRIX_INCLUDE_DIR - ${GMI_SIM_INCLUDE_DIR} - ${APF_SIM_INCLUDE_DIR} ${MESH_SIM_INCLUDE_DIR} ) -find_library(GMI_SIM_LIB gmi_sim) -find_library(APF_SIM_LIB apf_sim) - set(SIM_LIB_HINT ${SIMMETRIX_ROOT}/lib/x64_rhel7_gcc48) find_library(SIM_DISCRETE_LIB SimDiscrete ${SIM_LIB_HINT}) @@ -30,8 +36,6 @@ find_library(SIM_PS_KRNL_LIB pskernel ${SIM_LIB_HINT}/psKrnl) get_filename_component(SIM_PS_KRNL_LIB_DIR ${SIM_PS_KRNL_LIB} DIRECTORY) list(APPEND SIMMETRIX_LIBRARIES - "${GMI_SIM_LIB}" - "${APF_SIM_LIB}" "${SIM_DISCRETE_LIB}" "${SIM_EXPORT_LIB}" "${SIM_MESHING_LIB}" @@ -42,6 +46,13 @@ list(APPEND SIMMETRIX_LIBRARIES "${SIM_PS_KRNL_LIB}" ) +if (SCOREC) +list(APPEND SIMMETRIX_LIBRARIES + "${GMI_SIM_LIB}" + "${APF_SIM_LIB}" +) +endif() + if (PARASOLID) list(APPEND SIMMETRIX_LIBRARIES "${SIM_PARASOLID_LIB}" diff --git a/src/aux/Distributor.h b/src/aux/Distributor.h new file mode 100644 index 0000000..b5918a9 --- /dev/null +++ b/src/aux/Distributor.h @@ -0,0 +1,12 @@ +#ifndef PUMGEN_AUX_DISTRIBUTOR_H_ +#define PUMGEN_AUX_DISTRIBUTOR_H_ + +#include + +// no namespace for now + +constexpr std::size_t getChunksize(std::size_t total, int rank, int size) { + return (total / size) + (total % size > rank); +} + +#endif // PUMGEN_AUX_DISTRIBUTOR_H_ diff --git a/src/aux/InsphereCalculator.cpp b/src/aux/InsphereCalculator.cpp new file mode 100644 index 0000000..9f0735b --- /dev/null +++ b/src/aux/InsphereCalculator.cpp @@ -0,0 +1,165 @@ +#include "InsphereCalculator.h" +#include "third_party/MPITraits.h" +#include +#include +#include +#include +#include +#include + +#include "MPIConvenience.h" + +std::vector calculateInsphere(const std::vector& connectivity, + const std::vector& geometry, MPI_Comm comm) { + int commsize; + int commrank; + + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + + MPI_Datatype vertexType; + + MPI_Type_contiguous(3, MPI_DOUBLE, &vertexType); + MPI_Type_commit(&vertexType); + + std::vector outrequests(commsize); + std::vector inrequests(commsize); + + std::size_t localVertices = geometry.size() / 3; + + std::vector vertexDist(commsize + 1); + + MPI_Allgather(&localVertices, 1, tndm::mpi_type_t(), vertexDist.data() + 1, 1, + tndm::mpi_type_t(), comm); + + for (std::size_t i = 0; i < commsize; ++i) { + vertexDist[i + 1] += vertexDist[i]; + } + + std::vector> outidxmap(commsize); + + for (const auto& vertex : connectivity) { + auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); + auto position = std::distance(vertexDist.begin(), itPosition) - 1; + auto localVertex = vertex - vertexDist[position]; + + // only transfer each vertex coordinate once, and only if it's not already on the same rank + if (position != commrank && + outidxmap[position].find(localVertex) == outidxmap[position].end()) { + outidxmap[position][localVertex] = outrequests[position]; + ++outrequests[position]; + } + } + + MPI_Alltoall(outrequests.data(), 1, tndm::mpi_type_t(), inrequests.data(), 1, + tndm::mpi_type_t(), comm); + + std::vector outdisp(commsize + 1); + std::vector indisp(commsize + 1); + for (std::size_t i = 1; i < commsize + 1; ++i) { + outdisp[i] = outdisp[i - 1] + outrequests[i - 1]; + indisp[i] = indisp[i - 1] + inrequests[i - 1]; + } + + std::vector inidx(indisp[commsize]); + + { + std::vector outidx(outdisp[commsize]); + + std::vector counter(commsize); + + for (std::size_t i = 0; i < commsize; ++i) { + for (const auto& [localVertex, j] : outidxmap[i]) { + outidx[j + outdisp[i]] = localVertex; + } + } + + // transfer indices + + sparseAlltoallv(outidx.data(), outrequests.data(), outdisp.data(), + tndm::mpi_type_t(), inidx.data(), inrequests.data(), indisp.data(), + tndm::mpi_type_t(), comm); + } + + std::vector outvertices(3 * connectivity.size()); + + { + std::vector invertices(3 * indisp[commsize]); + + for (std::size_t i = 0; i < inidx.size(); ++i) { + for (int j = 0; j < 3; ++j) { + invertices[3 * i + j] = geometry[3 * inidx[i] + j]; + } + } + + // transfer vertices + + sparseAlltoallv(invertices.data(), inrequests.data(), indisp.data(), vertexType, + outvertices.data(), outrequests.data(), outdisp.data(), vertexType, comm); + } + + std::vector counter(commsize); + std::vector inspheres(connectivity.size() / 4); + + for (std::size_t i = 0; i < connectivity.size() / 4; ++i) { + std::array, 4> vertices; + for (int j = 0; j < 4; ++j) { + auto vertex = connectivity[i * 4 + j]; + auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); + auto position = std::distance(vertexDist.begin(), itPosition) - 1; + auto localVertex = vertex - vertexDist[position]; + if (position == commrank) { + vertices[j][0] = geometry[localVertex * 3 + 0]; + vertices[j][1] = geometry[localVertex * 3 + 1]; + vertices[j][2] = geometry[localVertex * 3 + 2]; + } else { + auto transferidx = outidxmap[position][localVertex] + outdisp[position]; + vertices[j][0] = outvertices[transferidx * 3 + 0]; + vertices[j][1] = outvertices[transferidx * 3 + 1]; + vertices[j][2] = outvertices[transferidx * 3 + 2]; + } + } + double a11 = vertices[1][0] - vertices[0][0]; + double a12 = vertices[1][1] - vertices[0][1]; + double a13 = vertices[1][2] - vertices[0][2]; + double a21 = vertices[2][0] - vertices[0][0]; + double a22 = vertices[2][1] - vertices[0][1]; + double a23 = vertices[2][2] - vertices[0][2]; + double a31 = vertices[3][0] - vertices[0][0]; + double a32 = vertices[3][1] - vertices[0][1]; + double a33 = vertices[3][2] - vertices[0][2]; + double b11 = vertices[2][0] - vertices[1][0]; + double b12 = vertices[2][1] - vertices[1][1]; + double b13 = vertices[2][2] - vertices[1][2]; + double b21 = vertices[3][0] - vertices[1][0]; + double b22 = vertices[3][1] - vertices[1][1]; + double b23 = vertices[3][2] - vertices[1][2]; + double det = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a13 * a22 * a31 - + a12 * a21 * a33 - a11 * a23 * a32; + double gram = std::abs(det); + + double a1a2x1 = a12 * a23 - a13 * a22; + double a1a2x2 = a13 * a21 - a11 * a23; + double a1a2x3 = a11 * a22 - a12 * a21; + double a1a3x1 = a12 * a33 - a13 * a32; + double a1a3x2 = a13 * a31 - a11 * a33; + double a1a3x3 = a11 * a32 - a12 * a31; + double a2a3x1 = a22 * a33 - a23 * a32; + double a2a3x2 = a23 * a31 - a21 * a33; + double a2a3x3 = a21 * a32 - a22 * a31; + double b1b2x1 = b12 * b23 - b13 * b22; + double b1b2x2 = b13 * b21 - b11 * b23; + double b1b2x3 = b11 * b22 - b12 * b21; + + double faces = std::sqrt(a1a2x1 * a1a2x1 + a1a2x2 * a1a2x2 + a1a2x3 * a1a2x3) + + std::sqrt(a1a3x1 * a1a3x1 + a1a3x2 * a1a3x2 + a1a3x3 * a1a3x3) + + std::sqrt(a2a3x1 * a2a3x1 + a2a3x2 * a2a3x2 + a2a3x3 * a2a3x3) + + std::sqrt(b1b2x1 * b1b2x1 + b1b2x2 * b1b2x2 + b1b2x3 * b1b2x3); + + inspheres[i] = gram / faces; + } + + MPI_Type_free(&vertexType); + + return inspheres; +} diff --git a/src/aux/InsphereCalculator.h b/src/aux/InsphereCalculator.h new file mode 100644 index 0000000..acd4b04 --- /dev/null +++ b/src/aux/InsphereCalculator.h @@ -0,0 +1,13 @@ +#ifndef PUMGEN_AUX_INSPHERE_CALCULATOR_H_ +#define PUMGEN_AUX_INSPHERE_CALCULATOR_H_ + +#include +#include +#include + +// assumes a contiguous distribution of all vertices over all processes + +std::vector calculateInsphere(const std::vector& connectivity, + const std::vector& geometry, MPI_Comm comm); + +#endif // PUMGEN_AUX_INSPHERE_CALCULATOR_H_ diff --git a/src/aux/MPIConvenience.cpp b/src/aux/MPIConvenience.cpp new file mode 100644 index 0000000..168a650 --- /dev/null +++ b/src/aux/MPIConvenience.cpp @@ -0,0 +1,108 @@ +#include "MPIConvenience.h" +#include +#include +#include + +namespace { +constexpr std::size_t DataIncrement = std::numeric_limits::max(); +constexpr int Tag = 1000; + +std::vector largeIsend(const void* sendbuf, std::size_t sendsize, std::size_t senddisp, + MPI_Datatype sendtype, int dest, int tag, MPI_Comm comm) { + int sendtypesizePre; + std::vector requests; + MPI_Type_size(sendtype, &sendtypesizePre); + std::size_t sendtypesize = sendtypesizePre; + for (std::size_t position = 0; position < sendsize; position += DataIncrement) { + int localSize = static_cast(std::min(sendsize - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Isend(reinterpret_cast(sendbuf) + (senddisp + position) * sendtypesize, + localSize, sendtype, dest, tag, comm, &requests.back()); + } + return requests; +} + +std::vector largeIrecv(void* recvbuf, std::size_t recvsize, std::size_t recvdisp, + MPI_Datatype recvtype, int dest, int tag, MPI_Comm comm) { + int recvtypesizePre; + std::vector requests; + MPI_Type_size(recvtype, &recvtypesizePre); + std::size_t recvtypesize = recvtypesizePre; + for (std::size_t position = 0; position < recvsize; position += DataIncrement) { + int localSize = static_cast(std::min(recvsize - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Irecv(reinterpret_cast(recvbuf) + (recvdisp + position) * recvtypesize, localSize, + recvtype, dest, tag, comm, &requests.back()); + } + return requests; +} +} // namespace + +void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const std::size_t* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + int commrank; + int sendtypesizePre; + int recvtypesizePre; + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + + std::vector requests; + requests.reserve(commsize * 2); + + // (no special handling of self-to-self comm at the moment) + + for (int i = 0; i < commsize; ++i) { + auto localRequests = largeIsend(sendbuf, sendsize[i], senddisp[i], sendtype, i, Tag, comm); + requests.insert(requests.end(), localRequests.begin(), localRequests.end()); + } + for (int i = 0; i < commsize; ++i) { + auto localRequests = largeIrecv(recvbuf, recvsize[i], recvdisp[i], recvtype, i, Tag, comm); + requests.insert(requests.end(), localRequests.begin(), localRequests.end()); + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); +} + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + MPI_Comm_size(comm, &commsize); + std::vector sendsizeSize(sendsize, sendsize + commsize); + std::vector recvsizeSize(recvsize, recvsize + commsize); + sparseAlltoallv(sendbuf, sendsizeSize.data(), senddisp, sendtype, recvbuf, recvsizeSize.data(), + recvdisp, recvtype, comm); +} + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, + MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + MPI_Comm_size(comm, &commsize); + std::vector senddispSize(senddisp, senddisp + commsize); + std::vector recvdispSize(recvdisp, recvdisp + commsize); + sparseAlltoallv(sendbuf, sendsize, senddispSize.data(), sendtype, recvbuf, recvsize, + recvdispSize.data(), recvtype, comm); +} + +void largeScatterv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, std::size_t recvsize, + MPI_Datatype recvtype, int root, MPI_Comm comm) { + std::vector requests; + int commrank; + int commsize; + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + if (commrank == root) { + requests.reserve(commsize + 1); + for (int i = 0; i < commsize; ++i) { + auto sendRequests = largeIsend(sendbuf, sendsize[i], senddisp[i], sendtype, i, Tag, comm); + requests.insert(requests.end(), sendRequests.begin(), sendRequests.end()); + } + } + auto recvRequests = largeIrecv(recvbuf, recvsize, 0, recvtype, root, Tag, comm); + requests.insert(requests.end(), recvRequests.begin(), recvRequests.end()); + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); +} diff --git a/src/aux/MPIConvenience.h b/src/aux/MPIConvenience.h new file mode 100644 index 0000000..b757e24 --- /dev/null +++ b/src/aux/MPIConvenience.h @@ -0,0 +1,23 @@ +#ifndef PUMGEN_AUX_MPI_CONVENIENCE_H_ +#define PUMGEN_AUX_MPI_CONVENIENCE_H_ + +#include +#include + +void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const std::size_t* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm); + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm); + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, + MPI_Datatype recvtype, MPI_Comm comm); + +void largeScatterv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, std::size_t recvsize, + MPI_Datatype recvtype, int root, MPI_Comm comm); + +#endif // PUMGEN_AUX_MPI_CONVENIENCE_H_ diff --git a/src/input/AnalysisAttributes.h b/src/input/AnalysisAttributes.h index e7d52b0..ed91eee 100644 --- a/src/input/AnalysisAttributes.h +++ b/src/input/AnalysisAttributes.h @@ -12,7 +12,7 @@ using namespace tinyxml2; struct faceBoundary { - faceBoundary(int faceID, int bcType) : faceID(faceID), bcType(bcType){}; + faceBoundary(int faceID, int bcType) : faceID(faceID), bcType(bcType) {}; int faceID; int bcType; }; diff --git a/src/input/ApfConvertWrapper.h b/src/input/ApfConvertWrapper.h deleted file mode 100644 index dc503c7..0000000 --- a/src/input/ApfConvertWrapper.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef APF_CONVERT_WRAPPER_H -#define APF_CONVERT_WRAPPER_H - -#include -#include - -// This is a small fix, because some time after PUMI (i.e. SCOREC/core) v2.2.7, -// the element ID was switched from a mere int to a named type apf::GiD, a typedef to long. -// Thus, PUMgen did not compile anymore, but to my knowledge, there was no flag to check the -// version. (of course, I could have missed some) To make PUMgen work with newer versions of PUMI, -// we therefore read out the type of the only function that actually was affected by the switch to -// apf::GiD and use that one in PUMgen. - -// inspired by https://stackoverflow.com/a/70954691 - -namespace internal { -template struct GetSecondFunctionArgument {}; - -template -struct GetSecondFunctionArgument { - using Type = Arg2; -}; - -using ElementIDConstPtr = GetSecondFunctionArgument::Type; -using ElementID = std::remove_const_t>; - -} // namespace internal - -using ElementID = internal::ElementID; - -#endif diff --git a/src/input/ApfNative.h b/src/input/ApfNative.h index 46ab3c2..a9794bf 100644 --- a/src/input/ApfNative.h +++ b/src/input/ApfNative.h @@ -21,9 +21,10 @@ #include "MeshInput.h" -class ApfNative : public MeshInput { +class ApfNative : public ApfMeshInput { public: - ApfNative(const char* mesh, const char* model = 0L) { + ApfNative(const char* mesh, int boundarySize, const char* model = 0L) + : ApfMeshInput(boundarySize) { if (model) gmi_register_mesh(); else { diff --git a/src/input/MeshAttributes.h b/src/input/MeshAttributes.h index 591e99c..86fe757 100644 --- a/src/input/MeshAttributes.h +++ b/src/input/MeshAttributes.h @@ -33,7 +33,7 @@ struct VelocityRefinementCube { VelocityRefinementCube(SimpleCuboid cuboid, double targetedFrequency, int bypassFindRegionAndUseGroup) : cuboid(cuboid), targetedFrequency(targetedFrequency), - bypassFindRegionAndUseGroup(bypassFindRegionAndUseGroup){}; + bypassFindRegionAndUseGroup(bypassFindRegionAndUseGroup) {}; SimpleCuboid cuboid; double targetedFrequency; diff --git a/src/input/MeshData.h b/src/input/MeshData.h new file mode 100644 index 0000000..230c209 --- /dev/null +++ b/src/input/MeshData.h @@ -0,0 +1,88 @@ + +#ifndef MESH_DATA_H +#define MESH_DATA_H + +#include +#include +#include + +#include "utils/logger.h" + +// TODO: to save space, the connectivity, geometry, group, and boundary arrays can also be +// constructed on the fly (for that, change the MeshData data structure to work on an iterator +// instead) + +/** + * Interface for mesh input + */ +class MeshData { + public: + virtual ~MeshData() = default; + + virtual std::size_t cellCount() = 0; + virtual std::size_t vertexCount() = 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; +}; + +class FullStorageMeshData : public MeshData { + public: + virtual ~FullStorageMeshData() = default; + + virtual std::size_t cellCount() { return cellCountValue; } + virtual std::size_t vertexCount() { return vertexCountValue; } + + virtual const std::vector& connectivity() { return connectivityData; } + virtual const std::vector& geometry() { return geometryData; } + virtual const std::vector& group() { return groupData; } + virtual const std::vector& boundary() { return boundaryData; } + + protected: + FullStorageMeshData(int boundarySize) : bndShift(boundarySize) {} + + std::size_t cellCountValue = 0; + std::size_t vertexCountValue = 0; + + std::vector connectivityData; + std::vector geometryData; + std::vector groupData; + std::vector boundaryData; + + int bndShift = 8; + + void setBoundary(std::size_t cell, int face, int value) { + if (bndShift < 0) { + boundaryData[cell * 4 + face] = value; + } else { + constexpr auto i32limit = std::numeric_limits::max(); + constexpr auto i64limit = std::numeric_limits::max(); + if (value < 0 || (value > i32limit && bndShift == 8) || + (value > i64limit && bndShift == 16)) { + logError() << "Cannot handle boundary condition" << value; + } + + boundaryData[cell] |= static_cast(value) << (face * bndShift); + } + } + + void setup(std::size_t cellCount, std::size_t vertexCount) { + cellCountValue = cellCount; + vertexCountValue = vertexCount; + + connectivityData.resize(cellCount * 4); + geometryData.resize(vertexCount * 3); + groupData.resize(cellCount); + if (bndShift < 0) { + boundaryData.resize(cellCount * 4); + } else { + boundaryData.resize(cellCount); + } + + // TODO: reconsider the times 4 or 3 + } +}; + +#endif // MESH_INPUT_H diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index b30a3e2..fae39a2 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -10,7 +10,11 @@ * @author Sebastian Rettenberger */ +#include "MeshData.h" +#include "utils/logger.h" #include +#include +#include #ifndef MESH_INTPUT_H #define MESH_INTPUT_H @@ -18,14 +22,98 @@ /** * Interface for mesh input */ -class MeshInput { +class ApfMeshInput : public FullStorageMeshData { + private: + template static void iterate(apf::Mesh* mesh, int dim, F&& function) { + apf::MeshIterator* it = mesh->begin(dim); + std::size_t index = 0; + while (apf::MeshEntity* element = mesh->iterate(it)) { + std::invoke(function, index, element); + ++index; + } + mesh->end(it); + } + protected: - apf::Mesh2* m_mesh; + apf::Mesh2* m_mesh = nullptr; - public: - virtual ~MeshInput() {} + ApfMeshInput(int boundarySize) : FullStorageMeshData(boundarySize) {} + public: apf::Mesh2* getMesh() { return m_mesh; } + + virtual ~ApfMeshInput() { + delete m_mesh; + m_mesh = nullptr; + } + + void generate() { + if (alignMdsMatches(m_mesh)) { + logWarning() << "Fixed misaligned matches"; + } + m_mesh->verify(); + apf::GlobalNumbering* vertexNum = apf::makeGlobal(apf::numberOwnedNodes(m_mesh, "vertices")); + apf::synchronize(vertexNum); + + apf::Sharing* sharing = apf::getSharing(m_mesh); + + setup(apf::countOwned(m_mesh, 3), apf::countOwned(m_mesh, 0)); + + // connectivity + iterate(m_mesh, 3, [&](auto index, auto* element) { + apf::NewArray vn; + apf::getElementNumbers(vertexNum, element, vn); + + for (int i = 0; i < 4; i++) { + connectivityData[index * 4 + i] = vn[i]; + } + }); + + std::size_t vertexIndex = 0; + + // geometry + iterate(m_mesh, 0, [&](auto index, auto* element) { + if (!sharing->isOwned(element)) { + return; + } + + apf::Vector3 point; + m_mesh->getPoint(element, 0, point); + double geometry[3]; + point.toArray(geometry); + + for (int i = 0; i < 3; ++i) { + geometryData[vertexIndex * 3 + i] = geometry[i]; + } + ++vertexIndex; + }); + + // groups + apf::MeshTag* groupTag = m_mesh->findTag("group"); + iterate(m_mesh, 3, [&](auto index, auto* element) { + int group; + m_mesh->getIntTag(element, groupTag, &group); + groupData[index] = group; + }); + + // boundary + apf::MeshTag* boundaryTag = m_mesh->findTag("boundary condition"); + iterate(m_mesh, 3, [&](auto index, auto* element) { + apf::Downward faces; + m_mesh->getDownward(element, 2, faces); + + for (int i = 0; i < 4; i++) { + if (m_mesh->hasTag(faces[i], boundaryTag)) { + int b; + m_mesh->getIntTag(faces[i], boundaryTag, &b); + + setBoundary(index, i, b); + } + } + }); + + delete sharing; + } }; #endif // MESH_INPUT_H diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index 934b6fb..fac05f7 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -18,33 +18,27 @@ #include #include -#include "ApfConvertWrapper.h" -#include -#include -#include -#include - +#include "MeshData.h" #include "utils/logger.h" -#include "MeshInput.h" +#include "MeshData.h" #include "NetCDFPartition.h" #include "ParallelVertexFilter.h" /** * Read PUMGen generated mesh files */ -class NetCDFMesh : public MeshInput { +class NetCDFMesh : public FullStorageMeshData { public: - NetCDFMesh(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) { + virtual ~NetCDFMesh() = default; + + NetCDFMesh(const char* meshFile, int boundarySize, MPI_Comm comm = MPI_COMM_WORLD) + : FullStorageMeshData(boundarySize) { int rank = 0; int nProcs = 1; MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &nProcs); - gmi_register_null(); - gmi_model* model = gmi_load(".null"); - m_mesh = apf::makeEmptyMdsMesh(model, 3, false); - int ncFile; checkNcError(nc_open_par(meshFile, NC_MPIIO, comm, MPI_INFO_NULL, &ncFile)); @@ -71,14 +65,8 @@ class NetCDFMesh : public MeshInput { checkNcError(nc_open_par(meshFile, NC_MPIIO, commIO, MPI_INFO_NULL, &ncFile)); } - PCU_Switch_Comm(commIO); - std::size_t nElements = 0; std::size_t nVertices = 0; - std::vector elements; - std::vector vertices; - std::vector boundaries; - std::vector groups; if (nLocalPart > 0) { std::vector partitions(nLocalPart); @@ -100,7 +88,7 @@ class NetCDFMesh : public MeshInput { bool useGroups = true; if (nc_inq_varid(ncFile, "element_group", &ncVarElemGroup) != NC_NOERR) { useGroups = false; - logWarning() << "No group found, using group 0 for all elements"; + logWarning(rank) << "No group found, using group 0 for all elements"; } else { collectiveAccess(ncFile, ncVarElemGroup); } @@ -161,94 +149,61 @@ class NetCDFMesh : public MeshInput { // Copy to the buffer std::vector elementsLocal(nElements * 4); - elements.resize(nElements * 4); - vertices.resize(nVertices * 3); + std::vector verticesLocal(nVertices * 3); - boundaries.resize(nElements * 4); - groups.resize(nElements); + std::size_t vertexOffset = 0; + for (std::size_t i = 0; i < nLocalPart; i++) { + std::copy_n(partitions[i].vertices(), partitions[i].nVertices() * 3, + verticesLocal.begin() + vertexOffset * 3); + vertexOffset += partitions[i].nVertices(); + } + + logInfo(rank) << "Running vertex filter"; + ParallelVertexFilter filter(commIO); + filter.filter(nVertices, verticesLocal); + + nVertices = filter.numLocalVertices(); + + setup(nElements, nVertices); + + std::copy_n(filter.localVertices().begin(), nVertices * 3, geometryData.begin()); std::size_t elementOffset = 0; - std::size_t vertexOffset = 0; + vertexOffset = 0; for (std::size_t i = 0; i < nLocalPart; i++) { #ifdef _OPENMP -#pragma omp parallel schedule(static) +#pragma omp parallel for schedule(static) #endif for (std::size_t j = 0; j < partitions[i].nElements() * 4; j++) { elementsLocal[elementOffset * 4 + j] = partitions[i].elements()[j] + vertexOffset; } - std::copy_n(partitions[i].vertices(), partitions[i].nVertices() * 3, - vertices.begin() + vertexOffset * 3); - partitions[i].convertBoundary(); - std::copy_n(partitions[i].boundaries(), partitions[i].nElements() * 4, - boundaries.begin() + elementOffset * 4); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t j = 0; j < partitions[i].nElements(); ++j) { + for (int k = 0; k < 4; ++k) { + setBoundary(elementOffset + j, k, partitions[i].boundaries()[j * 4 + k]); + } + } + std::copy_n(partitions[i].groups(), partitions[i].nElements(), - groups.begin() + elementOffset); + groupData.begin() + elementOffset); elementOffset += partitions[i].nElements(); vertexOffset += partitions[i].nVertices(); } - logInfo(rank) << "Running vertex filter"; - ParallelVertexFilter filter(commIO); - filter.filter(nVertices, vertices); - - // Create filtered vertex list - - nVertices = filter.numLocalVertices(); - vertices.resize(nVertices * 3); - std::copy_n(filter.localVertices(), nVertices * 3, vertices.begin()); - logInfo(rank) << "Converting local to global vertex identifier"; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel for #endif for (std::size_t i = 0; i < nElements * 4; i++) { - elements[i] = filter.globalIds()[elementsLocal[i]]; - } - } - - logInfo(rank) << "Constructing the mesh"; - apf::GlobalToVert vertMap; - apf::construct(m_mesh, elements.data(), nElements, apf::Mesh::TET, vertMap); - - apf::alignMdsRemotes(m_mesh); - apf::deriveMdsModel(m_mesh); - - logInfo(rank) << "Set coordinates in APF"; - apf::setCoords(m_mesh, vertices.data(), nVertices, vertMap); - - // Set boundaries - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::Adjacent adjacent; - m_mesh->getAdjacent(element, 2, adjacent); - - for (int j = 0; j < 4; j++) { - if (!boundaries[i * 4 + j]) - continue; - - m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); + connectivityData[i] = filter.globalIds()[elementsLocal[i]]; } - - i++; } - m_mesh->end(it); - - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - it = m_mesh->begin(3); - i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - m_mesh->setIntTag(element, groupTag, &groups[i]); - i++; - } - m_mesh->end(it); - - PCU_Switch_Comm(MPI_COMM_WORLD); } private: diff --git a/src/input/ParallelVertexFilter.h b/src/input/ParallelVertexFilter.h index 96b1136..144058e 100644 --- a/src/input/ParallelVertexFilter.h +++ b/src/input/ParallelVertexFilter.h @@ -30,6 +30,8 @@ #include "third_party/MPITraits.h" +#include "aux/MPIConvenience.h" + /** * Filters duplicate vertices in parallel */ @@ -162,15 +164,16 @@ class ParallelVertexFilter { } // Determine the (local and total) bucket size - std::vector bucketSize(m_numProcs); + std::vector bucketSize(m_numProcs); for (std::size_t i = 0; i < numVertices; i++) { ++bucketSize[bucket[i]]; } // Tell all processes what we are going to send them - std::vector recvSize(m_numProcs); + std::vector recvSize(m_numProcs); - MPI_Alltoall(bucketSize.data(), 1, MPI_INT, recvSize.data(), 1, MPI_INT, m_comm); + MPI_Alltoall(bucketSize.data(), 1, tndm::mpi_type_t(), recvSize.data(), 1, + tndm::mpi_type_t(), m_comm); std::size_t numSortVertices = 0; #ifdef _OPENMP @@ -192,16 +195,16 @@ class ParallelVertexFilter { // Allocate buffer for the vertices and exchange them std::vector sortVertices(3 * numSortVertices); - std::vector sDispls(m_numProcs); - std::vector rDispls(m_numProcs); + std::vector sDispls(m_numProcs); + std::vector rDispls(m_numProcs); sDispls[0] = 0; rDispls[0] = 0; for (int i = 1; i < m_numProcs; i++) { sDispls[i] = sDispls[i - 1] + bucketSize[i - 1]; rDispls[i] = rDispls[i - 1] + recvSize[i - 1]; } - MPI_Alltoallv(sendVertices.data(), bucketSize.data(), sDispls.data(), vertexType, - sortVertices.data(), recvSize.data(), rDispls.data(), vertexType, m_comm); + sparseAlltoallv(sendVertices.data(), bucketSize.data(), sDispls.data(), vertexType, + sortVertices.data(), recvSize.data(), rDispls.data(), vertexType, m_comm); // Chop the last 4 bits to avoid numerical errors roundVertices.resize(numSortVertices * 3); @@ -253,9 +256,9 @@ class ParallelVertexFilter { // Send result back std::vector globalIds(numVertices); - MPI_Alltoallv(gids.data(), recvSize.data(), rDispls.data(), tndm::mpi_type_t(), - globalIds.data(), bucketSize.data(), sDispls.data(), - tndm::mpi_type_t(), m_comm); + sparseAlltoallv(gids.data(), recvSize.data(), rDispls.data(), tndm::mpi_type_t(), + globalIds.data(), bucketSize.data(), sDispls.data(), + tndm::mpi_type_t(), m_comm); // Assign the global ids to the correct vertices m_globalIds.resize(numVertices); diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 5377314..362b26e 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -17,19 +17,20 @@ #include #endif // PARALLEL -#include "ApfConvertWrapper.h" -#include -#include -#include -#include +#include "aux/Distributor.h" +#include #include -#include "MeshInput.h" +#include "MeshData.h" +#include "utils/logger.h" /** * Read a mesh from a serial file */ -template class SerialMeshFile : public MeshInput { +template class SerialMeshFile : public FullStorageMeshData { + public: + virtual ~SerialMeshFile() = default; + private: #ifdef PARALLEL MPI_Comm m_comm; @@ -42,13 +43,13 @@ template class SerialMeshFile : public MeshInput { public: #ifdef PARALLEL - SerialMeshFile(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) - : m_comm(comm), m_meshReader(comm) { + SerialMeshFile(const char* meshFile, int boundarySize, MPI_Comm comm = MPI_COMM_WORLD) + : FullStorageMeshData(boundarySize), m_comm(comm), m_meshReader(comm) { init(); open(meshFile); } #else // PARALLEL - SerialMeshFile(const char* meshFile) { + SerialMeshFile(const char* meshFile, int boundarySize) : FullStorageMeshData(boundarySize) { init(); open(meshFile); } @@ -73,77 +74,27 @@ template class SerialMeshFile : public MeshInput { const std::size_t nVertices = m_meshReader.nVertices(); const std::size_t nElements = m_meshReader.nElements(); - std::size_t nLocalVertices = (nVertices + m_nProcs - 1) / m_nProcs; - std::size_t nLocalElements = (nElements + m_nProcs - 1) / m_nProcs; - if (m_rank == m_nProcs - 1) { - nLocalVertices = nVertices - (m_nProcs - 1) * nLocalVertices; - nLocalElements = nElements - (m_nProcs - 1) * nLocalElements; - } + const std::size_t nLocalVertices = getChunksize(nVertices, m_rank, m_nProcs); + const std::size_t nLocalElements = getChunksize(nElements, m_rank, m_nProcs); - gmi_register_null(); - gmi_model* model = gmi_load(".null"); - m_mesh = apf::makeEmptyMdsMesh(model, 3, false); - - // Create elements - apf::GlobalToVert vertMap; - { - std::vector elements; - { - std::vector elementsInt(nLocalElements * 4); - m_meshReader.readElements(elementsInt.data()); - elements = std::vector(elementsInt.begin(), elementsInt.end()); - } - logInfo(m_rank) << "Create APF connectivity"; - apf::construct(m_mesh, elements.data(), nLocalElements, apf::Mesh::TET, vertMap); - } + setup(nLocalElements, nLocalVertices); - apf::alignMdsRemotes(m_mesh); - apf::deriveMdsModel(m_mesh); + logInfo(m_rank) << "Read vertex coordinates"; + m_meshReader.readVertices(geometryData.data()); - // Set vertices - { - std::vector vertices(nLocalVertices * 3); - m_meshReader.readVertices(vertices.data()); - logInfo(m_rank) << "Set coordinates in APF"; - apf::setCoords(m_mesh, vertices.data(), nLocalVertices, vertMap); - } + logInfo(m_rank) << "Read cell vertices"; + m_meshReader.readElements(connectivityData.data()); - // Set boundaries - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - { - std::vector boundaries(nLocalElements * 4); - m_meshReader.readBoundaries(boundaries.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::Adjacent adjacent; - m_mesh->getAdjacent(element, 2, adjacent); - - for (int j = 0; j < 4; j++) { - if (!boundaries[i * 4 + j]) { - continue; - } - - m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); - } - - i++; - } - m_mesh->end(it); - } + logInfo(m_rank) << "Read cell groups"; + m_meshReader.readGroups(groupData.data()); - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - { - std::vector groups(nLocalElements); - m_meshReader.readGroups(groups.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - m_mesh->setIntTag(element, groupTag, &groups[i]); - i++; + logInfo(m_rank) << "Read boundary conditions"; + std::vector preBoundaryData(nLocalElements * 4); + 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]); } - m_mesh->end(it); } } }; diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index b28350a..416448a 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -6,8 +6,9 @@ * notice in the file 'COPYING' at the root directory of this package * and the copyright notice at https://github.com/SeisSol/PUMGen * - * @copyright 2017 Technical University of Munich + * @copyright 2017,2023 Technical University of Munich * @author Sebastian Rettenberger + * @author David Schneller */ #ifndef SIM_MOD_SUITE_H @@ -20,12 +21,6 @@ #include #include -#include -#include -#include -#include -#include - #include #include #ifdef BEFORE_SIM_18 @@ -48,14 +43,17 @@ #include "utils/path.h" #include "utils/progress.h" -#include "MeshInput.h" +#include "MeshData.h" #include "AnalysisAttributes.h" #include "EasiMeshSize.h" #include "MeshAttributes.h" +#include "ParallelVertexFilter.h" + #include #include +#include #include #include #include @@ -72,7 +70,7 @@ pAManager SModel_attManager(pModel model); * of this class * @todo Maybe add MS_setMaxEntities to limit the number of elements */ -class SimModSuite : public MeshInput { +class SimModSuite : public FullStorageMeshData { private: EasiMeshSize easiMeshSize; pGModel m_model; @@ -83,10 +81,11 @@ class SimModSuite : public MeshInput { bool m_log; public: - SimModSuite(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, - const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", - int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, - const char* logFile = 0L) { + SimModSuite(const char* modFile, int boundarySize, const char* cadFile = 0L, + const char* licenseFile = 0L, const char* meshCaseName = "mesh", + const char* analysisCaseName = "analysis", int enforceSize = 0, + const char* xmlFile = 0L, const bool analyseAR = false, const char* logFile = 0L) + : FullStorageMeshData(boundarySize) { // Init SimModSuite SimModel_start(); SimPartitionedMesh_start(0L, 0L); @@ -174,51 +173,171 @@ class SimModSuite : public MeshInput { analyse_mesh(); } - // Convert to APF mesh - apf::Mesh* tmpMesh = apf::createMesh(m_simMesh); - gmi_register_sim(); - gmi_model* model = gmi_import_sim(m_model); + logInfo(PMU_rank()) << "Iterating over mesh to get data..."; + int parts = PM_numParts(m_simMesh); + std::size_t vertexCount = 0; + std::size_t vertexIDCount = 0; + std::size_t cellCount = 0; + + std::vector vertexStart(parts + 1); + std::vector vertexIDStart(parts + 1); + std::vector cellStart(parts + 1); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Counting part" << i << "/" << parts; + vertexStart[i] = vertexCount; + vertexIDStart[i] = vertexIDCount; + cellStart[i] = cellCount; + { + size_t totalCount = 0; + size_t ownedCount = 0; + VIter vert_it = M_vertexIter(PM_mesh(m_simMesh, i)); + while (pVertex vertex = VIter_next(vert_it)) { + ++totalCount; + if (EN_isOwnerProc((pEntity)vertex)) { + ++ownedCount; + } + } + VIter_delete(vert_it); + + vertexCount += ownedCount; + vertexIDCount += totalCount; + } + { + size_t totalCount = 0; + size_t ownedCount = 0; + RIter reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + while (pRegion reg = RIter_next(reg_it)) { + ++totalCount; + if (EN_isOwnerProc((pEntity)reg)) { + ++ownedCount; + } + } + RIter_delete(reg_it); + + cellCount += ownedCount; + assert(totalCount == ownedCount); + } + } + vertexStart.back() = vertexCount; + cellStart.back() = cellCount; + vertexIDStart.back() = vertexIDCount; - logInfo(PMU_rank()) << "Converting mesh to APF"; - m_mesh = apf::createMdsMesh(model, tmpMesh); - apf::destroyMesh(tmpMesh); + logInfo(PMU_rank()) << "Local cells:" << cellCount; + logInfo(PMU_rank()) << "Local vertices:" << vertexCount; + logInfo(PMU_rank()) << "Local vertices (with duplicates):" << vertexIDCount; - // Set the boundary conditions from the geometric model AttCase_associate(analysisCase, 0L); - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - apf::MeshIterator* it = m_mesh->begin(2); - while (apf::MeshEntity* face = m_mesh->iterate(it)) { - apf::ModelEntity* modelFace = m_mesh->toModel(face); - if (m_mesh->getModelType(modelFace) != 2) - continue; - - pGEntity simFace = reinterpret_cast(modelFace); - - pAttribute attr = GEN_attrib(simFace, "boundaryCondition"); - if (attr) { - char* image = Attribute_imageClass(attr); - int boundary = parseBoundary(image); - Sim_deleteString(image); - - m_mesh->setIntTag(face, boundaryTag, &boundary); + + ParallelVertexFilter vertexFilter; + + { + std::vector vertexData(vertexIDCount * 3); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Processing part" << i << "/" << parts; + logInfo(PMU_rank()) << "Vertices:" << vertexStart[i] << "to" << vertexStart[i + 1]; + { + std::size_t indexID = vertexIDStart[i]; + VIter reg_it = M_vertexIter(PM_mesh(m_simMesh, i)); + while (pVertex vertex = VIter_next(reg_it)) { + EN_setID((pEntity)vertex, indexID); + double xyz[3]; + V_coord(vertex, xyz); + vertexData[indexID * 3 + 0] = xyz[0]; + vertexData[indexID * 3 + 1] = xyz[1]; + vertexData[indexID * 3 + 2] = xyz[2]; + ++indexID; + } + VIter_delete(reg_it); + } } + vertexFilter.filter(vertexIDCount, vertexData); } - m_mesh->end(it); + vertexCount = vertexFilter.numLocalVertices(); + + logInfo(PMU_rank()) << "Local vertices (after filtering):" << vertexCount; + + setup(cellCount, vertexCount); + std::copy(vertexFilter.localVertices().begin(), vertexFilter.localVertices().end(), + geometryData.begin()); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Processing part" << i << "/" << parts; + logInfo(PMU_rank()) << "Connectivity:" << cellStart[i] << "to" << cellStart[i + 1]; + { + std::size_t index = cellStart[i]; + RIter reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + pRegion reg; + while ((reg = RIter_next(reg_it))) { + if (EN_isOwnerProc((pEntity)reg)) { + EN_setID((pEntity)reg, index); + + pPList vertices = R_vertices(reg, 1); + void* iter = nullptr; + int j = 0; + while (pVertex vertex = (pVertex)PList_next(vertices, &iter)) { + connectivityData[index * 4 + j] = vertexFilter.globalIds()[EN_id((pEntity)vertex)]; + ++j; + } + PList_delete(vertices); + ++index; + } + } + RIter_delete(reg_it); + } - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - it = m_mesh->begin(3); - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::ModelEntity* modelRegion = m_mesh->toModel(element); + logInfo(PMU_rank()) << "Groups"; + { + int groupIdx = 0; + GRIter modelRegionIter = GM_regionIter(m_model); + while (pGRegion mreg = GRIter_next(modelRegionIter)) { + RIter reg_it = M_classifiedRegionIter(PM_mesh(m_simMesh, i), mreg); + while (pRegion reg = RIter_next(reg_it)) { + if (EN_isOwnerProc((pEntity)reg)) { + groupData[EN_id((pEntity)reg)] = groupIdx; + } + } + RIter_delete(reg_it); - pGRegion simRegion = reinterpret_cast(modelRegion); - std::unordered_map::const_iterator i = groupMap.find(simRegion); - if (i == groupMap.end()) - logError() << "Mesh element with unknown region found."; + ++groupIdx; + } + GRIter_delete(modelRegionIter); + } - m_mesh->setIntTag(element, groupTag, &i->second); + logInfo(PMU_rank()) << "Boundaries"; + { + GFIter modelFaceIter = GM_faceIter(m_model); + while (pGFace mface = GFIter_next(modelFaceIter)) { + pAttribute attr = GEN_attrib(mface, "boundaryCondition"); + if (attr != nullptr) { + char* image = Attribute_imageClass(attr); + int boundary = parseBoundary(image); + FIter face_it = M_classifiedFaceIter(PM_mesh(m_simMesh, i), mface, 1); + while (pFace face = FIter_next(face_it)) { + pRegion r1 = F_region(face, 0); + pRegion r2 = F_region(face, 1); + + // suboptimal: we go over all faces here... + if (r1 != nullptr && EN_isOwnerProc((pEntity)r1)) { + for (int k = 0; k < 4; ++k) { + if (R_face(r1, k) == face) { + setBoundary(EN_id((pEntity)r1), k, boundary); + } + } + } + if (r2 != nullptr && EN_isOwnerProc((pEntity)r2)) { + for (int k = 0; k < 4; ++k) { + if (R_face(r2, k) == face) { + setBoundary(EN_id((pEntity)r2), k, boundary); + } + } + } + } + FIter_delete(face_it); + Sim_deleteString(image); + } + } + GFIter_delete(modelFaceIter); + } } - m_mesh->end(it); AttCase_unassociate(analysisCase); diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h new file mode 100644 index 0000000..ad07f30 --- /dev/null +++ b/src/input/SimModSuiteApf.h @@ -0,0 +1,689 @@ +/** + * @file + * This file is part of PUMGen + * + * For conditions of distribution and use, please see the copyright + * notice in the file 'COPYING' at the root directory of this package + * and the copyright notice at https://github.com/SeisSol/PUMGen + * + * @copyright 2017 Technical University of Munich + * @author Sebastian Rettenberger + */ + +#ifndef SIM_MOD_SUITE_APF_H +#define SIM_MOD_SUITE_APF_H + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#ifdef BEFORE_SIM_18 +#include +#include +#include +#else +#include +#include +#include +#endif +#include +#include +#ifdef PARASOLID +#include +#endif +#include + +#include "utils/logger.h" +#include "utils/path.h" +#include "utils/progress.h" + +#include "MeshInput.h" + +#include "AnalysisAttributes.h" +#include "EasiMeshSize.h" +#include "MeshAttributes.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// forward declare +pAManager SModel_attManager(pModel model); + +/** + * @todo Currently it is not supported to create more than one instance + * of this class + * @todo Maybe add MS_setMaxEntities to limit the number of elements + */ +class SimModSuiteApf : public ApfMeshInput { + private: + EasiMeshSize easiMeshSize; + pGModel m_model; + + pParMesh m_simMesh; + + /** Enable Simmetrix logging file */ + bool m_log; + + public: + SimModSuiteApf(const char* modFile, int boundarySize, const char* cadFile = 0L, + const char* licenseFile = 0L, const char* meshCaseName = "mesh", + const char* analysisCaseName = "analysis", int enforceSize = 0, + const char* xmlFile = 0L, const bool analyseAR = false, const char* logFile = 0L) + : ApfMeshInput(boundarySize) { + // Init SimModSuite + SimModel_start(); + SimPartitionedMesh_start(0L, 0L); + if (logFile) { + m_log = true; + Sim_logOn(logFile); + } else + m_log = false; + Sim_readLicenseFile(licenseFile); + MS_init(); + SimDiscrete_start(0); +#ifdef PARASOLID + SimParasolid_start(1); +#endif + Sim_setMessageHandler(messageHandler); + + // Load CAD + logInfo(PMU_rank()) << "Loading model"; + + std::string smodFile = modFile; + if (cadFile != 0L) { + loadCAD(modFile, cadFile); + } else if (smodFile.substr(smodFile.find_last_of(".") + 1) == "smd") { + loadCAD(modFile, cadFile); + } else { + logError() << "unsupported model file"; + } + + // Create group lookup table + std::unordered_map groupMap; + GRIter regionIt = GM_regionIter(m_model); + while (pGRegion region = GRIter_next(regionIt)) { + groupMap[region] = static_cast(groupMap.size() + 1); + } + + // Extract cases + logInfo(PMU_rank()) << "Extracting cases"; + pACase meshCase, analysisCase; + MeshAttributes MeshAtt; + + if (xmlFile != nullptr) { + + // Read mesh Attributes from xml file + int numFaces = GM_numFaces(m_model); + MeshAtt.init(xmlFile); + AnalysisAttributes AnalysisAtt(xmlFile, numFaces); + setCases(m_model, meshCase, analysisCase, MeshAtt, AnalysisAtt, groupMap); + } else { + extractCases(m_model, meshCase, meshCaseName, analysisCase, analysisCaseName); + } + + m_simMesh = PM_new(0, m_model, PMU_size()); + + pProgress prog = Progress_new(); + Progress_setCallback(prog, progressHandler); + + // create the mesh + logInfo(PMU_rank()) << "Starting the surface mesher"; + pSurfaceMesher surfaceMesher = SurfaceMesher_new(meshCase, m_simMesh); + if (xmlFile != nullptr) { + SurfaceMesher_setSmoothing(surfaceMesher, MeshAtt.surfaceSmoothingLevel); + SurfaceMesher_setSmoothType(surfaceMesher, static_cast(MeshAtt.surfaceSmoothingType)); + SurfaceMesher_setFaceRotationLimit(surfaceMesher, MeshAtt.surfaceFaceRotationLimit); + SurfaceMesher_setSnapForDiscrete(surfaceMesher, MeshAtt.surfaceSnap); + } + SurfaceMesher_execute(surfaceMesher, prog); + SurfaceMesher_delete(surfaceMesher); + + PM_setTotalNumParts(m_simMesh, PMU_size()); + + logInfo(PMU_rank()) << "Starting the volume mesher"; + pVolumeMesher volumeMesher = VolumeMesher_new(meshCase, m_simMesh); + if (xmlFile != nullptr) { + VolumeMesher_setSmoothing(volumeMesher, MeshAtt.volumeSmoothingLevel); + VolumeMesher_setSmoothType(volumeMesher, static_cast(MeshAtt.volumeSmoothingType)); + VolumeMesher_setOptimization(volumeMesher, MeshAtt.VolumeMesherOptimization); + } + VolumeMesher_setEnforceSize(volumeMesher, enforceSize); + VolumeMesher_execute(volumeMesher, prog); + VolumeMesher_delete(volumeMesher); + + Progress_delete(prog); + + if (analyseAR) { + analyse_mesh(); + } + + // Convert to APF mesh + apf::Mesh* tmpMesh = apf::createMesh(m_simMesh); + gmi_register_sim(); + gmi_model* model = gmi_import_sim(m_model); + + logInfo(PMU_rank()) << "Converting mesh to APF"; + m_mesh = apf::createMdsMesh(model, tmpMesh); + apf::destroyMesh(tmpMesh); + + // Set the boundary conditions from the geometric model + AttCase_associate(analysisCase, 0L); + apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); + apf::MeshIterator* it = m_mesh->begin(2); + while (apf::MeshEntity* face = m_mesh->iterate(it)) { + apf::ModelEntity* modelFace = m_mesh->toModel(face); + if (m_mesh->getModelType(modelFace) != 2) + continue; + + pGEntity simFace = reinterpret_cast(modelFace); + + pAttribute attr = GEN_attrib(simFace, "boundaryCondition"); + if (attr) { + char* image = Attribute_imageClass(attr); + int boundary = parseBoundary(image); + Sim_deleteString(image); + + m_mesh->setIntTag(face, boundaryTag, &boundary); + } + } + m_mesh->end(it); + + // Set groups + apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); + it = m_mesh->begin(3); + while (apf::MeshEntity* element = m_mesh->iterate(it)) { + apf::ModelEntity* modelRegion = m_mesh->toModel(element); + + pGRegion simRegion = reinterpret_cast(modelRegion); + std::unordered_map::const_iterator i = groupMap.find(simRegion); + if (i == groupMap.end()) + logError() << "Mesh element with unknown region found."; + + m_mesh->setIntTag(element, groupTag, &i->second); + } + m_mesh->end(it); + + AttCase_unassociate(analysisCase); + + // Delete cases + MS_deleteMeshCase(meshCase); + MS_deleteMeshCase(analysisCase); + } + + virtual ~SimModSuiteApf() { + M_release(m_simMesh); + // We cannot delete the model here because it is still + // connected to the mesh + // GM_release(m_model); + + // Finalize SimModSuite +#ifdef PARASOLID + SimParasolid_stop(1); +#endif + SimDiscrete_stop(0); + MS_exit(); + Sim_unregisterAllKeys(); + if (m_log) + Sim_logOff(); + SimPartitionedMesh_stop(); + SimModel_stop(); + } + + private: + pACase extractCase(pAManager attMngr, const char* name) { + pACase acase = AMAN_findCase(attMngr, name); + if (!acase) + logError() << "Case" << std::string(name) << "not found."; + + AttCase_setModel(acase, m_model); + + return acase; + } + + private: + /** + * @todo Make this private as soon as APF does no longer need it + */ + static unsigned int parseBoundary(const char* boundaryCondition) { + if (strcmp(boundaryCondition, "freeSurface") == 0) + return 1; + if (strcmp(boundaryCondition, "dynamicRupture") == 0) + return 3; + if (strcmp(boundaryCondition, "absorbing") == 0) + return 5; + + // check if boundaryCondition starts with pattern + std::string pattern = "BC"; + std::string sboundaryCondition(boundaryCondition, 11); + if (sboundaryCondition.find(pattern) == 0) { + std::string sNumber = sboundaryCondition.substr(pattern.length()); + return stoi(sNumber); + } + logError() << "Unknown boundary condition" << boundaryCondition; + return -1; + } + + static void messageHandler(int type, const char* msg) { + switch (type) { + case Sim_InfoMsg: + // Show sim info messages as debug messages + logDebug(PMU_rank()) << "SimModeler:" << msg; + break; + case Sim_DebugMsg: + // Ignore sim debug messages + break; + case Sim_WarningMsg: + logWarning(PMU_rank()) << "SimModeler:" << msg; + break; + case Sim_ErrorMsg: + // Use warning because error will abort the program + logWarning() << "SimModeler:" << msg; + break; + } + } + + static void progressHandler(const char* what, int level, int startVal, int endVal, int currentVal, + void* ignore) { + if (PMU_rank() != 0) + return; + if (endVal != 0) { + switch (currentVal) { + case -2: + // task is started, do nothing + logInfo(PMU_rank()) << "Progress:" << what << ", 0" + << "/" << endVal; + break; + case -1: + // end of the task + logInfo(PMU_rank()) << "Progress:" << what << ", done"; + break; + default: + logInfo(PMU_rank()) << "Progress:" << what << "," << currentVal << "/" << endVal; + break; + } + } else { + if (currentVal == -2) + logInfo(PMU_rank()) << "Progress:" << what; + } + logDebug() << what << level << startVal << endVal << currentVal; + } + + private: + void extractCases(pGModel m_model, pACase& meshCase, const char* meshCaseName, + pACase& analysisCase, const char* analysisCaseName) { + logInfo(PMU_rank()) << "Extracting cases"; + pAManager attMngr = SModel_attManager(m_model); + + MeshingOptions meshingOptions; + meshCase = MS_newMeshCase(m_model); + + pACase meshCaseFile = extractCase(attMngr, meshCaseName); + AttCase_associate(meshCaseFile, nullptr); + MS_processSimModelerMeshingAtts(meshCaseFile, meshCase, &meshingOptions); + AttCase_setModel(meshCase, m_model); + + analysisCase = extractCase(attMngr, analysisCaseName); + pPList children = AttNode_children(analysisCase); + void* iter = 0L; + while (pANode child = static_cast(PList_next(children, &iter))) + AttCase_setModel(static_cast(child), m_model); + PList_delete(children); + } + + private: + void setMeshSize(pGModel model, pACase& meshCase, MeshAttributes& MeshAtt) { + // Set mesh size on vertices, edges, surfaces and regions + const char* element_name[4] = {"vertex", "edge", "face", "region"}; + const ElementType element_type[4] = {ElementType::vertex, ElementType::edge, ElementType::face, + ElementType::region}; + + for (int element_type_id = 0; element_type_id < 4; element_type_id++) { + std::list, double>> lpair_lId_MSize = + MeshAtt.getMSizeList(element_type[element_type_id]); + for (const auto& pair_lId_MSize : lpair_lId_MSize) { + double MSize = pair_lId_MSize.second; + std::list lElements = pair_lId_MSize.first; + for (const auto& element_id : lElements) { + pGEntity entity = GM_entityByTag(model, element_type_id, element_id); + if (entity == nullptr) { + logError() << element_name[element_type_id] << "id:" << element_id + << "not found in model."; + } else { + MS_setMeshSize(meshCase, entity, 1, MSize, nullptr); + logInfo(PMU_rank()) << element_name[element_type_id] << "id:" << element_id + << ", MSize =" << MSize; + } + } + } + } + } + + void setCases(pGModel model, pACase& meshCase, pACase& analysisCase, MeshAttributes& MeshAtt, + AnalysisAttributes& AnalysisAtt, std::unordered_map groupMap) { + + logInfo(PMU_rank()) << "Setting cases"; + // ------------------------------ Set boundary conditions + // ------------------------------ + + // Create a new manager. The manager is responsible for creating + // new attributes, saving/retrieving to/from file, and to look up + // cases, AttModels etc. + pAManager attMngr = AMAN_new(); + + // Create a case. A case serves as a grouping mechanism. It will (should) + // contain all the attributes that make up a complete analysis of + // some kind, although the system has no way of verifying this. + analysisCase = AMAN_newCase(attMngr, "analysis", "", (pModel)model); + + // Now we create an attribute information node. This can be seen + // as an "attribute generator object", as we will later on create an + // actual attribute for each geometric face from this attribute + // information node. We name the attribute T1, and give it the + // information type "boundaryCondition". + + // With this code we can tag any BC, not only 1,3,5 + int numFaces = GM_numFaces(model); + std::set UniquefaceBound; + for (auto const& f : AnalysisAtt.faceBound) { + UniquefaceBound.insert(f.bcType); + } + int numBC = UniquefaceBound.size(); + std::map aBC; + + for (auto const& ufb : UniquefaceBound) { + std::string strboundaryCondition = "BC"; + std::ostringstream stringstream; + stringstream << std::setw(9) << std::setfill('0') << ufb; + std::string sNumber = stringstream.str(); + + pAttInfoVoid iBC = AMAN_newAttInfoVoid(attMngr, "BC", "boundaryCondition"); + + strboundaryCondition.append(sNumber); + AttNode_setImageClass((pANode)iBC, strboundaryCondition.c_str()); + AttCase_addNode(analysisCase, (pANode)iBC); + aBC[ufb] = AttCase_newModelAssoc(analysisCase, (pANode)iBC); + } + + for (auto const& fb : AnalysisAtt.faceBound) { + // Get the face + pGEntity face = GM_entityByTag(model, 2, fb.faceID + 1); + if (face == nullptr) { + logError() << "faceid:" << fb.faceID + 1 << "not found in model."; + } else { + // Add the face to the model association. Note that we passed + // the Attribute Information Node into the Model Association + // at the time when the Model Association was created. That prepares + // the creation of the AttributeVoid on the face as soon as the + // association process is started + if (fb.bcType != 0) { + logInfo(PMU_rank()) << "faceBound[" << fb.faceID + 1 << "] =" << fb.bcType; + AMA_addGEntity(aBC[fb.bcType], face); + } + } + } + + // ------------------------------ Set meshing parameters + // ------------------------------ + + meshCase = MS_newMeshCase(model); + + // Set global mesh size + pModelItem modelDomain = GM_domain(model); + if (MeshAtt.globalMSize > 0) { + logInfo(PMU_rank()) << "globalMSize =" << MeshAtt.globalMSize; + // ( , , <1=absolute, 2=relative>, , ) + MS_setMeshSize(meshCase, modelDomain, 1, MeshAtt.globalMSize, nullptr); + } + + // Set mesh size on vertices, edges, surfaces and regions + setMeshSize(model, meshCase, MeshAtt); + + if (MeshAtt.velocityAwareRefinementSettings.isVelocityAwareRefinementOn()) { + logInfo(PMU_rank()) << "Enabling velocity aware meshing"; + easiMeshSize = EasiMeshSize(MeshAtt.velocityAwareRefinementSettings, model, groupMap); + auto easiMeshSizeFunc = [](pSizeAttData sadata, void* userdata) { + auto* easiMeshSize = static_cast(userdata); + std::array pt{}; + int haspt = SizeAttData_point(sadata, pt.data()); + if (!haspt) { + logError() << "!haspt"; + } + return easiMeshSize->getMeshSize(pt); + }; + // set the user-defined function for isotropic size + MS_setSizeAttFunc(meshCase, "setCustomMeshSize", easiMeshSizeFunc, &easiMeshSize); + // Relative anisotropic size is set for the entire model through the + // function anisoSize + MS_setMeshSize(meshCase, GM_domain(model), MS_userDefinedType | 1, 0, "setCustomMeshSize"); + } + + if (MeshAtt.gradation > 0) { + // Set gradation relative + logInfo(PMU_rank()) << "Gradation rate =" << MeshAtt.gradation; + MS_setGlobalSizeGradationRate(meshCase, MeshAtt.gradation); + } + if (MeshAtt.vol_AspectRatio > 0) { + // Set target equivolume AspectRatio + logInfo(PMU_rank()) << "Target equivolume AspectRatio =" << MeshAtt.vol_AspectRatio; + MS_setVolumeShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + MeshAtt.vol_AspectRatio); + } + if (MeshAtt.area_AspectRatio > 0) { + // Set target equiarea AspectRatio + logInfo(PMU_rank()) << "Target equiarea AspectRatio =" << MeshAtt.area_AspectRatio; +#ifdef BEFORE_SIM_11 + MS_setSurfaceShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + MeshAtt.area_AspectRatio); +#else + MS_setSurfaceShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + ElementType_Triangle, MeshAtt.area_AspectRatio); +#endif + } + for (auto& mycube : MeshAtt.lCube) { + logInfo(PMU_rank()) << "Cube mesh refinement: " << mycube.CubeMSize; + logInfo(PMU_rank()) << "Center" << mycube.CubeCenter[0] << " " << mycube.CubeCenter[1] << " " + << mycube.CubeCenter[2]; + logInfo(PMU_rank()) << "Width" << mycube.CubeWidth[0] << " " << mycube.CubeWidth[1] << " " + << mycube.CubeWidth[2]; + logInfo(PMU_rank()) << "Height" << mycube.CubeHeight[0] << " " << mycube.CubeHeight[1] << " " + << mycube.CubeHeight[2]; + logInfo(PMU_rank()) << "Depth" << mycube.CubeDepth[0] << " " << mycube.CubeDepth[1] << " " + << mycube.CubeDepth[2]; + MS_addCubeRefinement(meshCase, mycube.CubeMSize, &mycube.CubeCenter[0], &mycube.CubeWidth[0], + &mycube.CubeHeight[0], &mycube.CubeDepth[0]); + } + + if (MeshAtt.MeshSizePropagationDistance > 0.0) { + for (auto& iElem : MeshAtt.lFaceIdMeshSizePropagation) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "MeshSizeProp faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "MeshSizeProp faceid:" << iElem + << ", distance =" << MeshAtt.MeshSizePropagationDistance + << ", scaling factor =" << MeshAtt.MeshSizePropagationScalingFactor; +#ifdef BEFORE_SIM_15 + MS_setMeshSizePropagation(meshCase, face, 1, MeshAtt.MeshSizePropagationDistance, + MeshAtt.MeshSizePropagationScalingFactor); +#else + MS_setMeshSizePropagation(meshCase, face, 2, 1, MeshAtt.MeshSizePropagationDistance, + MeshAtt.MeshSizePropagationScalingFactor); +#endif + } + } + } + for (auto& iElem : MeshAtt.lFaceIdUseDiscreteMesh) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "UseDiscreteMesh; faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "UseDiscreteMesh; faceid, noModification:" << iElem + << MeshAtt.UseDiscreteMesh_noModification; + MS_useDiscreteGeometryMesh(meshCase, face, MeshAtt.UseDiscreteMesh_noModification); + // MS_limitSurfaceMeshModification(meshCase,face,UseDiscreteMesh_noModification); + } + } + for (auto& iElem : MeshAtt.lFaceIdNoMesh) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "No Mesh; faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "No Mesh; faceid:" << iElem; + MS_setNoMesh(meshCase, face, 1); + } + } + for (auto& iElem : MeshAtt.lRegionIdNoMesh) { + pGEntity region = GM_entityByTag(model, 3, iElem); + if (region == nullptr) { + logError() << "No Mesh; regionid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "No Mesh; regionid:" << iElem; + MS_setNoMesh(meshCase, region, 1); + } + } + } + + private: + void loadCAD(const char* modFile, const char* cadFile) { + // Load CAD + std::string sCadFile; + if (cadFile) + sCadFile = cadFile; + else { + sCadFile = modFile; + utils::StringUtils::replaceLast(sCadFile, ".smd", "_nat.x_t"); + } + pNativeModel nativeModel = 0L; +#ifdef PARASOLID + if (utils::Path(sCadFile).exists()) + nativeModel = ParasolidNM_createFromFile(sCadFile.c_str(), 0); +#endif + m_model = GM_load(modFile, nativeModel, 0L); + nativeModel = GM_nativeModel(m_model); + + if (nativeModel) + NM_release(nativeModel); + +#ifdef BEFORE_SIM_18 + // check for model errors + pPList modelErrors = PList_new(); + if (!GM_isValid(m_model, 1, modelErrors)) { + void* iter = 0L; + while (pSimError err = static_cast(PList_next(modelErrors, &iter))) { + logInfo(PMU_rank()) << " Error code: " << SimError_code(err) << std::endl; + logInfo(PMU_rank()) << " Error string: " << SimError_toString(err) << std::endl; + } + logError() << "Input model is not valid"; + } + PList_delete(modelErrors); +#else + // check for model infos + pPList modelInfos = PList_new(); + if (!GM_isValid(m_model, 1, modelInfos)) { + void* iter = 0L; + while (pSimInfo info = static_cast(PList_next(modelInfos, &iter))) { + logInfo(PMU_rank()) << " Info code: " << SimInfo_code(info) << std::endl; + logInfo(PMU_rank()) << " Info string: " << SimInfo_toString(info) << std::endl; + } + logError() << "Input model is not valid"; + } + PList_delete(modelInfos); +#endif + + // check for self-intersecting CAD mesh + pProgress prog = Progress_new(); + Progress_setCallback(prog, progressHandler); + pMesh mesh = DM_getMesh(static_cast(m_model)); + pPList entities = PList_new(); + if (MS_checkMeshIntersections(mesh, entities, prog) > 0) { + + const char* element_name[4] = {"vertex", "edge", "face", "region"}; + + for (int i = 0; i < PList_size(entities); i += 2) { + pEntity firstEnt = (pEntity)PList_item(entities, i); + pEntity secondEnt = (pEntity)PList_item(entities, i + 1); + logInfo(PMU_rank()) << "Self-intersection between" << element_name[EN_whatInType(firstEnt)] + << GEN_tag(EN_whatIn(firstEnt)) << "and" + << element_name[EN_whatInType(secondEnt)] + << GEN_tag(EN_whatIn(secondEnt)); + double centroid[3], centroid2[3]; + EN_centroid(firstEnt, ¢roid[0]); + EN_centroid(secondEnt, ¢roid2[0]); + logInfo(PMU_rank()) << "centroids" << centroid[0] << centroid[1] << centroid[2] << "and " + << centroid2[0] << centroid2[1] << centroid2[2] << std::flush; + } + logError() << PList_size(entities) / 2 << "Self-intersection(s) detected in CAD mesh"; + } + PList_delete(entities); + Progress_delete(prog); + } + + private: + void analyse_mesh() { + int num_bins = 8; + double AR[8] = {0, 2, 4, 6, 10, 20, 40, 100}; + long int AR_vol_bins[num_bins]; + long int skew_area_bins[num_bins]; + memset(AR_vol_bins, 0, num_bins * sizeof(long int)); + memset(skew_area_bins, 0, num_bins * sizeof(long int)); + + // Accumulate equivolume AspectRatio data + RIter reg_it; + int num_partMeshes = PM_numParts(m_simMesh); + pRegion reg; + double maxAR = 0, AR_global = 0, elAR; + int kAR; + for (int i = 0; i < num_partMeshes; i++) { + reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + while ((reg = RIter_next(reg_it))) { + kAR = num_bins; + for (int k = 1; k < num_bins; k++) { + elAR = R_aspectRatio(reg); + if (elAR < AR[k]) { + kAR = k; + break; + } + } + AR_vol_bins[kAR - 1]++; + maxAR = std::max(maxAR, elAR); + } + } + RIter_delete(reg_it); + + // Print the statistics + logInfo(PMU_rank()) << "AR statistics:"; + MPI_Allreduce(&maxAR, &AR_global, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + logInfo(PMU_rank()) << "AR max:" << AR_global; + logInfo(PMU_rank()) << "AR (target: < ~10):"; + long int bin_global; + for (int i = 0; i < num_bins - 1; i++) { + MPI_Allreduce(&AR_vol_bins[i], &bin_global, 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD); + logInfo(PMU_rank()) << std::fixed << std::setprecision(2) << "[" << AR[i] << "," << AR[i + 1] + << "):" << bin_global; + } + MPI_Allreduce(&AR_vol_bins[num_bins - 1], &bin_global, 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD); + logInfo(PMU_rank()) << std::fixed << std::setprecision(2) << "[" << AR[num_bins - 1] + << ",inf):" << bin_global; + } +}; + +#endif // SIM_MOD_SUITE_APF_H diff --git a/src/meshreader/FidapReader.cpp b/src/meshreader/FidapReader.cpp deleted file mode 100644 index fe39208..0000000 --- a/src/meshreader/FidapReader.cpp +++ /dev/null @@ -1,24 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#include "FidapReader.h" - -const char* FidapReader::FIDAP_FILE_ID = "** FIDAP NEUTRAL FILE"; -const char* FidapReader::NODAL_COORDINATES = "NODAL COORDINATES"; -const char* FidapReader::BOUNDARY_CONDITIONS = "BOUNDARY CONDITIONS"; -const char* FidapReader::ELEMENT_GROUPS = "ELEMENT GROUPS"; - -const char* FidapReader::ZONE_GROUP = "GROUP:"; -const char* FidapReader::ZONE_ELEMENTS = "ELEMENTS:"; -const char* FidapReader::ZONE_NODES = "NODES:"; -const char* FidapReader::ZONE_GEOMETRY = "GEOMETRY:"; -const char* FidapReader::ZONE_TYPE = "TYPE:"; diff --git a/src/meshreader/FidapReader.h b/src/meshreader/FidapReader.h deleted file mode 100644 index cbc39cc..0000000 --- a/src/meshreader/FidapReader.h +++ /dev/null @@ -1,414 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#ifndef FIDAP_READER_H -#define FIDAP_READER_H - -#include -#include -#include -#include -#include -#include - -#include "utils/stringutils.h" - -#include "MeshReader.h" - -struct FidapBoundaryFace { - /** The vertices of the face */ - std::array vertices; - /** The type of the boundary */ - int type; -}; - -/** Defines a face by three vertices */ -struct FaceVertex { - FaceVertex() {} - - FaceVertex(std::size_t v1, std::size_t v2, std::size_t v3) { - vertices[0] = v1; - vertices[1] = v2; - vertices[2] = v3; - std::sort(vertices.begin(), vertices.end()); - } - - FaceVertex(const std::array& v) { - std::copy(v.begin(), v.end(), vertices.begin()); - std::sort(vertices.begin(), vertices.end()); - } - - bool operator<(const FaceVertex& other) const { - return vertices[0] < other.vertices[0] || - (vertices[0] == other.vertices[0] && vertices[1] < other.vertices[1]) || - (vertices[0] == other.vertices[0] && vertices[1] == other.vertices[1] && - vertices[2] < other.vertices[2]); - } - - std::array vertices; -}; - -/** Defines a face by element and face side */ -struct FaceElement { - FaceElement() { - element = -1; - side = -1; - } - - FaceElement(std::size_t e, int s) { - element = e; - side = s; - } - - std::size_t element; - int side; -}; - -class FidapReader : public MeshReader { - private: - struct ElementSection : FileSection { - /** Group id for elements in this section */ - int group; - }; - - struct BoundarySection : FileSection { - /** Type of the boundary */ - int type; - }; - - // Section descriptions - FileSection m_vertices; - std::vector m_elements; - std::vector m_boundaries; - - public: - void open(const char* meshFile) { - MeshReader::open(meshFile); - - std::string line; - - // Read header information - getline(m_mesh, line); - utils::StringUtils::trim(line); - if (line != FIDAP_FILE_ID) - logError() << "Not an Fidap mesh file:" << meshFile; - - std::string name; - getline(m_mesh, name); // Internal name - - getline(m_mesh, line); // Skip line: - // VERSION: x.y - getline(m_mesh, line); // Date - getline(m_mesh, line); // Empty line - - std::size_t nElements, nZones, t, dimensions; - m_mesh >> m_vertices.nLines; - m_mesh >> nElements; - m_mesh >> nZones; - m_mesh >> t; // I don't know what this variable contains - m_mesh >> dimensions; - if (dimensions != 3) - logError() << "Fidap file does not contain a 3 dimensional mesh"; - getline(m_mesh, line); // Read newline - - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - - // Find seek positions - // Vertices - getline(m_mesh, line); - if (line.find(NODAL_COORDINATES) == std::string::npos) - logError() << "Invalid Fidap format: Coordinates expected, found" << line; - m_vertices.seekPosition = m_mesh.tellg(); - getline(m_mesh, line); - m_vertices.lineSize = line.length() + 1; - m_mesh.seekg(m_vertices.seekPosition + m_vertices.nLines * m_vertices.lineSize); - - // Boundary conditions (currently ignored) - getline(m_mesh, line); - if (line.find(BOUNDARY_CONDITIONS) == std::string::npos) - logError() << "Invalid Fidap format: Boundary conditions expected, found" << line; - getline(m_mesh, line); - - // Element groups - getline(m_mesh, line); - if (line.find(ELEMENT_GROUPS) == std::string::npos) - logError() << "Invalid Fidap format: Element groups expected, found" << line; - - for (std::size_t i = 0; i < nZones; i++) { - FileSection sec; - - // Group (unused) - m_mesh >> name; - if (name.find(ZONE_GROUP) == std::string::npos) - logError() << "Invalid Fidap format: Expected group, found" << name; - unsigned int id; - m_mesh >> id; - // Elements - m_mesh >> name; - if (name.find(ZONE_ELEMENTS) == std::string::npos) - logError() << "Invalid Fidap format: Expected elements, found" << name; - m_mesh >> sec.nLines; - // Nodes - m_mesh >> name; - if (name.find(ZONE_NODES) == std::string::npos) - logError() << "Invalid Fidap format: Expected nodes, found" << name; - unsigned int nodeType; - m_mesh >> nodeType; - // Geometry (unused) - m_mesh >> name; - if (name.find(ZONE_GEOMETRY) == std::string::npos) - logError() << "Invalid Fidap format: Expected geometry, found" << name; - unsigned int geoType; - m_mesh >> geoType; - // Type (unused) - m_mesh >> name; - if (name.find(ZONE_TYPE) == std::string::npos) - logError() << "Invalid Fidap format: Expected type, found" << name; - unsigned int type; - m_mesh >> type; - getline(m_mesh, line); // Read newline - - // Entity name - m_mesh >> name; - m_mesh >> name; - std::string entityName; - getline(m_mesh, entityName); - utils::StringUtils::trim(entityName); - - sec.seekPosition = m_mesh.tellg(); - getline(m_mesh, line); - sec.lineSize = line.size() + 1; - - switch (nodeType) { - case 3: { - BoundarySection boundary; - boundary.seekPosition = sec.seekPosition; - boundary.nLines = sec.nLines; - boundary.lineSize = sec.lineSize; - - boundary.type = convertType(entityName); - - m_boundaries.push_back(boundary); - break; - } - case 4: { - ElementSection element; - element.seekPosition = sec.seekPosition; - element.nLines = sec.nLines; - element.lineSize = sec.lineSize; - - element.group = m_elements.size() + 1; - - m_elements.push_back(element); - break; - } - default: - logError() << "Unsupported node type" << nodeType << "found"; - } - - m_mesh.seekg(sec.seekPosition + sec.nLines * sec.lineSize); - } - } - - std::size_t nVertices() const { return m_vertices.nLines; } - - std::size_t nElements() const { - std::size_t count = 0; - for (const auto& element : m_elements) { - count += element.nLines; - } - - return count; - } - - std::size_t nBoundaries() const { - std::size_t count = 0; - for (const auto& boundary : m_boundaries) { - count += boundary.nLines; - } - - return count; - } - - void readVertices(std::size_t start, std::size_t count, double* vertices) { - m_mesh.clear(); - - m_mesh.seekg(m_vertices.seekPosition + start * m_vertices.lineSize + m_vertices.lineSize - - 3 * COORDINATE_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - for (int j = 0; j < 3; j++) - m_mesh >> vertices[i * 3 + j]; - - m_mesh.seekg(m_vertices.lineSize - 3 * COORDINATE_SIZE, std::fstream::cur); - } - } - - void readElements(std::size_t start, std::size_t count, std::size_t* elements) { - m_mesh.clear(); - - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_elements.begin(); section != m_elements.end() && section->nLines < start; - section++) - start -= section->nLines; - - m_mesh.seekg(section->seekPosition + start * section->lineSize + section->lineSize - - 4 * VERTEX_ID_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - - m_mesh.seekg(section->seekPosition + section->lineSize - 4 * VERTEX_ID_SIZE - 1); - } - - for (int j = 0; j < 4; j++) { - m_mesh >> elements[i * 4 + j]; - elements[i * 4 + j]--; - } - - m_mesh.seekg(section->lineSize - 4 * VERTEX_ID_SIZE, std::fstream::cur); - - start++; // Line in the current section - } - } - - /** - * Read group ids from start to start+count. The group ids are sorted - * in the same order as the elements. - */ - void readGroups(std::size_t start, std::size_t count, int* groups) { - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_elements.begin(); section != m_elements.end() && section->nLines < start; - section++) - start -= section->nLines; - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - } - - groups[i] = section->group; - - start++; // "Line" in the current section - } - } - - /** - * Read all groups - */ - void readGroups(int* groups) { readGroups(0, nElements(), groups); } - - void readBoundaries(std::size_t start, std::size_t count, FidapBoundaryFace* boundaries) { - m_mesh.clear(); - - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_boundaries.begin(); section != m_boundaries.end() && section->nLines < start; - section++) - start -= section->nLines; - - m_mesh.seekg(section->seekPosition + start * section->lineSize + section->lineSize - - 3 * VERTEX_ID_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - - m_mesh.seekg(section->seekPosition + section->lineSize - 3 * VERTEX_ID_SIZE - 1); - } - - for (int j = 0; j < 3; j++) { - m_mesh >> boundaries[i].vertices[j]; - boundaries[i].vertices[j]--; - } - boundaries[i].type = section->type; - - // Seek to next position - m_mesh.seekg(section->lineSize - 3 * VERTEX_ID_SIZE, std::fstream::cur); - - start++; // Line in the current section - } - } - - // TODO void readBoundaries(std::size_t* boundaries) {} - - public: - /** - * Creates a map from vertex definition of a face to a element/side definition - * - * @param elements Elements (4 vertices per element) - * @param n Number of elements - * @param[out] map The map - */ - static void createFaceMap(const std::size_t* elements, std::size_t n, - std::map& map) { - for (std::size_t i = 0; i < n; i++) { - FaceVertex v1(elements[i * 4], elements[i * 4 + 1], elements[i * 4 + 2]); - FaceElement e1(i, 0); - map[v1] = e1; - - FaceVertex v2(elements[i * 4], elements[i * 4 + 1], elements[i * 4 + 3]); - FaceElement e2(i, 1); - map[v2] = e2; - - FaceVertex v3(elements[i * 4 + 1], elements[i * 4 + 2], elements[i * 4 + 3]); - FaceElement e3(i, 2); - map[v3] = e3; - - FaceVertex v4(elements[i * 4], elements[i * 4 + 2], elements[i * 4 + 3]); - FaceElement e4(i, 3); - map[v4] = e4; - } - } - - private: - static int convertType(const std::string& name) { - std::istringstream ss(name.substr(2)); - - int type; - ss >> type; - - return type % 100; - } - - /** Number of characters required to store a coordinate */ - static const size_t COORDINATE_SIZE = 20ul; - /** Number of characters required to store a vertex id */ - static const size_t VERTEX_ID_SIZE = 8ul; - - static const char* FIDAP_FILE_ID; - static const char* NODAL_COORDINATES; - static const char* BOUNDARY_CONDITIONS; - static const char* ELEMENT_GROUPS; - - static const char* ZONE_GROUP; - static const char* ZONE_ELEMENTS; - static const char* ZONE_NODES; - static const char* ZONE_GEOMETRY; - static const char* ZONE_TYPE; -}; - -#endif // FIDAP_READER_H diff --git a/src/meshreader/GMSHBuilder.h b/src/meshreader/GMSHBuilder.h index 3c547ac..f4c35c9 100644 --- a/src/meshreader/GMSHBuilder.h +++ b/src/meshreader/GMSHBuilder.h @@ -40,7 +40,7 @@ template class GMSHBuilder : public tndm::GMSHMeshBuilder { std::vector bcs; void setNumVertices(std::size_t numVertices) { vertices.resize(numVertices); } - void setVertex(long id, std::array const& x) { + void setVertex(long id, std::array const& x) { for (std::size_t i = 0; i < D; ++i) { vertices[id][i] = x[i]; } diff --git a/src/meshreader/ParallelFidapReader.h b/src/meshreader/ParallelFidapReader.h deleted file mode 100644 index dea856a..0000000 --- a/src/meshreader/ParallelFidapReader.h +++ /dev/null @@ -1,237 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#ifndef PARALLEL_FIDAP_READER -#define PARALLEL_FIDAP_READER - -#include - -#include - -#include "FidapReader.h" -#include "ParallelMeshReader.h" -#include "third_party/MPITraits.h" - -class ParallelFidapReader : public ParallelMeshReader { - private: - /** Face map required for boundaries */ - std::map m_faceMap; - - public: - ParallelFidapReader(MPI_Comm comm = MPI_COMM_WORLD) : ParallelMeshReader(comm) {} - - ParallelFidapReader(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) - : ParallelMeshReader(meshFile, comm) {} - - void readElements(std::size_t* elements) { - ParallelMeshReader::readElements(elements); - - // Create the face map - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - const std::size_t maxChunkSize = chunkSize; - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - - FidapReader::createFaceMap(elements, chunkSize, m_faceMap); - for (auto& face : m_faceMap) { - face.second.element += m_rank * maxChunkSize; - } - } - - /** - * @copydoc ParallelGambitReader::readGroups - */ - void readGroups(int* groups) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - - if (m_rank == 0) { - // Allocate second buffer so we can read and send in parallel - std::vector tempGroups(chunkSize); - int* groups2 = tempGroups.data(); - if (m_nProcs % 2 == 0) - // swap once so we have the correct buffer at the end - swap(groups, groups2); - - MPI_Request request = MPI_REQUEST_NULL; - - for (int i = 1; i < m_nProcs - 1; i++) { - logInfo() << "Reading group information part" << i << "of" << m_nProcs; - m_serialReader.readGroups(i * chunkSize, chunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - MPI_Isend(groups, chunkSize, MPI_INT, i, 0, m_comm, &request); - swap(groups, groups2); - } - - if (m_nProcs > 1) { - // Read last one - const std::size_t lastChunkSize = nElements() - (m_nProcs - 1) * chunkSize; - logInfo() << "Reading group information part" << (m_nProcs - 1) << "of" << m_nProcs; - m_serialReader.readGroups((m_nProcs - 1) * chunkSize, lastChunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - MPI_Isend(groups, lastChunkSize, MPI_INT, m_nProcs - 1, 0, m_comm, &request); - swap(groups, groups2); - } - - // Finally read the first part - logInfo() << "Reading group information part" << m_nProcs << "of" << m_nProcs; - m_serialReader.readGroups(0, chunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - } else { - if (m_rank == m_nProcs - 1) - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - - MPI_Recv(groups, chunkSize, MPI_INT, 0, 0, m_comm, MPI_STATUS_IGNORE); - } - } - - /** - * @copydoc ParallelGambitReader::readBoundaries - */ - void readBoundaries(int* boundaries) { - // Create the distributed face map - int vertexChunk = (nVertices() + m_nProcs - 1) / m_nProcs; - PCU_Comm_Begin(); - for (std::map::iterator i = m_faceMap.begin(); i != m_faceMap.end();) { - int proc = i->first.vertices[1] / vertexChunk; - if (proc == m_rank) - i++; - else { - PCU_COMM_PACK(proc, i->first); - PCU_COMM_PACK(proc, i->second); - m_faceMap.erase(i++); - } - } - PCU_Comm_Send(); - while (PCU_Comm_Receive()) { - FaceVertex v; - PCU_COMM_UNPACK(v); - FaceElement e; - PCU_COMM_UNPACK(e); - m_faceMap[v] = e; - } - - // Distribute the faces - std::vector facePos; - std::vector faceType; - - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - const std::size_t maxChunkSize = chunkSize; - - if (m_rank == 0) { - std::vector faces(maxChunkSize); - - std::vector> aggregator(m_nProcs - 1); - std::vector sizes(m_nProcs - 1); - std::vector requests((m_nProcs - 1) * 2, MPI_REQUEST_NULL); - - const std::size_t nChunks = (nBoundaries() + chunkSize - 1) / chunkSize; - for (std::size_t i = 0; i < nChunks; i++) { - logInfo() << "Reading boundary conditions part" << (i + 1) << "of" << nChunks; - - if (i == nChunks - 1) { - chunkSize = nBoundaries() - (nChunks - 1) * chunkSize; - } - - m_serialReader.readBoundaries(i * maxChunkSize, chunkSize, faces.data()); - - // Wait for all sending from last iteration - MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); - for (int j = 0; j < m_nProcs - 1; j++) { - aggregator[j].clear(); - } - - // Sort boundary conditions into the corresponding aggregator - for (std::size_t j = 0; j < chunkSize; j++) { - FaceVertex v(faces[j].vertices); - - if (v.vertices[1] < vertexChunk) { - facePos.push_back(m_faceMap.at(v)); - faceType.push_back(faces[j].type); - } else { - // Face for another processor - // Serials the struct to make sending easier - int proc = v.vertices[1] / vertexChunk; - aggregator[proc - 1].insert(aggregator[proc - 1].end(), v.vertices.begin(), - v.vertices.end()); - aggregator[proc - 1].push_back(faces[j].type); - } - } - - // Send send aggregated values - for (int j = 0; j < m_nProcs - 1; j++) { - if (aggregator[j].empty()) - continue; - - sizes[j] = aggregator[j].size() / 4; // 3 vertices + face type - MPI_Isend(&sizes[j], 1, tndm::mpi_type_t(), j + 1, 0, m_comm, - &requests[j * 2]); - MPI_Isend(&aggregator[j][0], aggregator[j].size(), MPI_INT, j + 1, 0, m_comm, - &requests[j * 2 + 1]); - } - } - - MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); - - // Send finish signal to all other processes - std::fill(sizes.begin(), sizes.end(), 0); - for (int i = 0; i < m_nProcs - 1; i++) { - MPI_Isend(&sizes[i], 1, tndm::mpi_type_t(), i + 1, 0, m_comm, &requests[i]); - } - MPI_Waitall(m_nProcs - 1, requests.data(), MPI_STATUSES_IGNORE); - } else { - if (m_rank == m_nProcs - 1) - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - - // Allocate enough memory - std::vector buf(chunkSize * 4); - - while (true) { - std::size_t size; - MPI_Recv(&size, 1, tndm::mpi_type_t(), 0, 0, m_comm, MPI_STATUS_IGNORE); - if (size == 0) - // Finished - break; - - MPI_Recv(buf.data(), size * 4, tndm::mpi_type_t(), 0, 0, m_comm, - MPI_STATUS_IGNORE); - - for (std::size_t i = 0; i < size * 4; i += 4) { - FaceVertex v(buf[i], buf[i + 1], buf[i + 2]); - facePos.push_back(m_faceMap.at(v)); - faceType.push_back(buf[i + 3]); - } - } - } - - // Distribute the faces to the final ranks - PCU_Comm_Begin(); - for (std::size_t i = 0; i < facePos.size(); i++) { - if (facePos[i].element / maxChunkSize == static_cast(m_rank)) - boundaries[(facePos[i].element % maxChunkSize) * 4 + facePos[i].side] = faceType[i]; - else { - PCU_COMM_PACK(facePos[i].element / maxChunkSize, facePos[i]); - PCU_COMM_PACK(facePos[i].element / maxChunkSize, faceType[i]); - } - } - PCU_Comm_Send(); - while (PCU_Comm_Receive()) { - FaceElement e; - PCU_COMM_UNPACK(e); - int type; - PCU_COMM_UNPACK(type); - boundaries[(e.element % maxChunkSize) * 4 + e.side] = type; - } - } -}; - -#endif // PARALLEL_FIDAP_READER diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index 66b2902..99388b8 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -11,6 +11,9 @@ #include #include +#include "aux/Distributor.h" +#include "aux/MPIConvenience.h" + namespace puml { template class ParallelGMSHReader { @@ -40,7 +43,7 @@ template class ParallelGMSHReader { } MPI_Bcast(&nVertices_, 1, tndm::mpi_type_t(), 0, comm_); - MPI_Bcast(&nElements_, 1, tndm::mpi_type_t(), 0, comm_); + MPI_Bcast(&nElements_, 1, tndm::mpi_type_t(), 0, comm_); } [[nodiscard]] std::size_t nVertices() const { return nVertices_; } @@ -140,19 +143,17 @@ template class ParallelGMSHReader { MPI_Comm_rank(comm_, &rank); MPI_Comm_size(comm_, &procs); - const std::size_t chunkSize = (numElements + procs - 1) / procs; - auto sendcounts = std::vector(procs); - auto displs = std::vector(procs + 1); - std::fill(sendcounts.begin(), sendcounts.end(), numPerElement * chunkSize); - sendcounts.back() = numPerElement * (numElements - (procs - 1u) * chunkSize); + auto sendcounts = std::vector(procs); + auto displs = std::vector(procs + 1); displs[0] = 0; for (int i = 0; i < procs; ++i) { + sendcounts[i] = numPerElement * getChunksize(numElements, i, procs); displs[i + 1] = displs[i] + sendcounts[i]; } auto recvcount = sendcounts[rank]; - MPI_Scatterv(sendbuf, sendcounts.data(), displs.data(), tndm::mpi_type_t(), recvbuf, - recvcount, tndm::mpi_type_t(), 0, comm_); + largeScatterv(sendbuf, sendcounts.data(), displs.data(), tndm::mpi_type_t(), recvbuf, + recvcount, tndm::mpi_type_t(), 0, comm_); } MPI_Comm comm_; diff --git a/src/meshreader/ParallelGambitReader.h b/src/meshreader/ParallelGambitReader.h index 2bfc970..b5b3d8c 100644 --- a/src/meshreader/ParallelGambitReader.h +++ b/src/meshreader/ParallelGambitReader.h @@ -17,6 +17,7 @@ #include "GambitReader.h" #include "ParallelMeshReader.h" +#include "aux/Distributor.h" #include "third_party/MPITraits.h" @@ -37,7 +38,7 @@ class ParallelGambitReader : public ParallelMeshReader { * This is a collective operation. */ void readGroups(int* groups) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; + std::size_t chunkSize = getChunksize(nElements(), m_rank, m_nProcs); if (m_rank == 0) { std::size_t maxChunkSize = chunkSize; @@ -50,9 +51,7 @@ class ParallelGambitReader : public ParallelMeshReader { for (int i = 0; i < m_nProcs; i++) { logInfo() << "Reading group information part" << (i + 1) << "of" << m_nProcs; - if (i == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } + chunkSize = getChunksize(nElements(), i, m_nProcs); m_serialReader.readGroups(i * maxChunkSize, chunkSize, map.data()); @@ -92,10 +91,6 @@ class ParallelGambitReader : public ParallelMeshReader { MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); } else { - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - // Allocate enough memory std::vector buf(chunkSize * 2); @@ -126,7 +121,7 @@ class ParallelGambitReader : public ParallelMeshReader { * @todo Only tetrahedral meshes are supported */ void readBoundaries(int* boundaries) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; + std::size_t chunkSize = getChunksize(nBoundaries(), m_rank, m_nProcs); if (m_rank == 0) { std::size_t maxChunkSize = chunkSize; @@ -140,9 +135,7 @@ class ParallelGambitReader : public ParallelMeshReader { for (std::size_t i = 0; i < nChunks; i++) { logInfo() << "Reading boundary conditions part" << (i + 1) << "of" << nChunks; - if (i == nChunks - 1) { - chunkSize = nBoundaries() - (nChunks - 1) * chunkSize; - } + chunkSize = getChunksize(nBoundaries(), i, m_nProcs); m_serialReader.readBoundaries(i * maxChunkSize, chunkSize, faces.data()); @@ -188,10 +181,6 @@ class ParallelGambitReader : public ParallelMeshReader { } MPI_Waitall(m_nProcs - 1, requests.data(), MPI_STATUSES_IGNORE); } else { - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - // Allocate enough memory std::vector buf(chunkSize * 2); diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 885c747..98a5a97 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -10,12 +10,12 @@ * @author Sebastian Rettenberger */ -#include -#include +#include #include #include #include +#include #include #include #include @@ -25,10 +25,17 @@ #include +#ifdef USE_SCOREC +#include #include #include // #include +#include "input/ApfNative.h" #include +#ifdef USE_SIMMOD +#include "input/SimModSuiteApf.h" +#endif +#endif #include #include "third_party/MPITraits.h" @@ -39,16 +46,16 @@ #ifdef USE_NETCDF #include "input/NetCDFMesh.h" #endif // USE_NETCDF -#include "input/ApfNative.h" #ifdef USE_SIMMOD #include "input/SimModSuite.h" #endif // USE_SIMMOD #include "meshreader/GMSH4Parser.h" -#include "meshreader/ParallelFidapReader.h" #include "meshreader/ParallelGMSHReader.h" #include "meshreader/ParallelGambitReader.h" #include "third_party/GMSH2Parser.h" +#include "aux/InsphereCalculator.h" + template static TT _checkH5Err(TT&& status, const char* file, int line) { if (status < 0) { logError() << utils::nospace << "An HDF5 error occurred (" << file << ": " << line << ")"; @@ -67,37 +74,29 @@ static int ilog(std::size_t value, int expbase = 1) { return count; } -template static void iterate(apf::Mesh* mesh, int dim, F&& function) { - apf::MeshIterator* it = mesh->begin(dim); - while (apf::MeshEntity* element = mesh->iterate(it)) { - std::invoke(function, element); - } - mesh->end(it); -} - constexpr std::size_t NoSecondDim = 0; -template -static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf::Mesh* mesh, +template +static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, void* mesh, int meshdim, hid_t h5memtype, hid_t h5outtype, std::size_t chunk, std::size_t localSize, std::size_t globalSize, bool reduceInts, - int filterEnable, std::size_t filterChunksize) { - constexpr std::size_t SecondSize = std::max(SecondDim, static_cast(1)); - constexpr std::size_t Dimensions = SecondDim == 0 ? 1 : 2; - const std::size_t chunkSize = chunk / SecondSize / sizeof(T); + int filterEnable, std::size_t filterChunksize, std::size_t secondDim) { + const std::size_t secondSize = std::max(secondDim, static_cast(1)); + const std::size_t dimensions = secondDim == 0 ? 1 : 2; + const std::size_t chunkSize = chunk / secondSize / sizeof(T); const std::size_t bufferSize = std::min(localSize, chunkSize); - std::vector data(SecondSize * bufferSize); + std::vector data(secondSize * bufferSize); std::size_t rounds = (localSize + chunkSize - 1) / chunkSize; MPI_Allreduce(MPI_IN_PLACE, &rounds, 1, tndm::mpi_type_t(), MPI_MAX, MPI_COMM_WORLD); - hsize_t globalSizes[2] = {globalSize, SecondDim}; - hid_t h5space = H5Screate_simple(Dimensions, globalSizes, nullptr); + hsize_t globalSizes[2] = {globalSize, secondDim}; + hid_t h5space = H5Screate_simple(dimensions, globalSizes, nullptr); - hsize_t sizes[2] = {bufferSize, SecondDim}; - hid_t h5memspace = H5Screate_simple(Dimensions, sizes, nullptr); + hsize_t sizes[2] = {bufferSize, secondDim}; + hid_t h5memspace = H5Screate_simple(dimensions, sizes, nullptr); std::size_t offset = localSize; @@ -106,7 +105,7 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: offset -= localSize; hsize_t start[2] = {offset, 0}; - hsize_t count[2] = {bufferSize, SecondDim}; + hsize_t count[2] = {bufferSize, secondDim}; hid_t h5dxlist = H5Pcreate(H5P_DATASET_XFER); checkH5Err(h5dxlist); @@ -126,7 +125,7 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: hid_t h5filter = H5P_DEFAULT; if (filterEnable > 0) { h5filter = checkH5Err(H5Pcreate(H5P_DATASET_CREATE)); - hsize_t chunk[2] = {std::min(filterChunksize, bufferSize), SecondDim}; + hsize_t chunk[2] = {std::min(filterChunksize, bufferSize), secondDim}; checkH5Err(H5Pset_chunk(h5filter, 2, chunk)); if (filterEnable == 1 && std::is_integral_v) { checkH5Err(H5Pset_scaleoffset(h5filter, H5Z_SO_INT, H5Z_SO_INT_MINBITS_DEFAULT)); @@ -142,8 +141,6 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: std::size_t written = 0; std::size_t index = 0; - apf::MeshIterator* it = mesh->begin(meshdim); - hsize_t nullstart[2] = {0, 0}; for (std::size_t i = 0; i < rounds; ++i) { @@ -151,15 +148,7 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: start[0] = offset + written; count[0] = std::min(localSize - written, bufferSize); - while (apf::MeshEntity* element = mesh->iterate(it)) { - if (handler(element, data.begin() + index * SecondSize)) { - ++index; - } - - if (index >= count[0]) { - break; - } - } + std::copy_n(handler.begin() + written * secondSize, count[0] * secondSize, data.begin()); checkH5Err(H5Sselect_hyperslab(h5memspace, H5S_SELECT_SET, nullstart, nullptr, count, nullptr)); @@ -170,8 +159,6 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: written += count[0]; } - mesh->end(it); - if (filterEnable > 0) { checkH5Err(H5Pclose(h5filter)); } @@ -187,21 +174,22 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: int main(int argc, char* argv[]) { int rank = 0; int processes = 1; + int mpithreadstate = MPI_THREAD_SINGLE; - MPI_Init(&argc, &argv); + MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &mpithreadstate); MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &processes); +#ifdef USE_SCOREC PCU_Comm_Init(); +#endif // Parse command line arguments utils::Args args; - const char* source[] = {"gambit", "fidap", "msh2", "msh4", "netcdf", "apf", "simmodsuite"}; + const char* source[] = {"gambit", "msh2", "msh4", "netcdf", + "apf", "simmodsuite", "simmodsuite-apf"}; args.addEnumOption("source", source, 's', "Mesh source (default: gambit)", false); - args.addOption("dump", 'd', "Dump APF mesh before partitioning it", utils::Args::Required, false); - args.addOption("model", 0, "Dump/Load a specific model file", utils::Args::Required, false); - args.addOption("vtk", 0, "Dump mesh to VTK files", utils::Args::Required, false); const char* filters[] = {"none", "scaleoffset", "deflate1", "deflate2", "deflate3", "deflate4", "deflate5", "deflate6", @@ -215,9 +203,9 @@ int main(int argc, char* argv[]) { utils::Args::Required, false); args.addOption("chunksize", 0, "Chunksize for writing (default=1 GiB).", utils::Args::Required, false); - const char* boundarytypes[] = {"int32", "int64"}; + const char* boundarytypes[] = {"int32", "int64", "i32", "i64", "i32x4"}; args.addEnumOption("boundarytype", boundarytypes, 0, - "Type for writing out boundary data (default: int32).", false); + "Type for writing out boundary data (default: i32).", false); args.addOption("license", 'l', "License file (only used by SimModSuite)", utils::Args::Required, false); @@ -285,14 +273,32 @@ int main(int argc, char* argv[]) { int boundaryType = args.getArgument("boundarytype", 0); hid_t boundaryDatatype; int faceOffset; - if (boundaryType == 0) { + int secondShape = NoSecondDim; + std::string boundaryFormatAttr = ""; + int boundaryPrecision = 0; + std::string secondDimBoundary = ""; + if (boundaryType == 0 || boundaryType == 2) { boundaryDatatype = H5T_STD_I32LE; faceOffset = 8; - logInfo(rank) << "Using 32-bit integer boundary type conditions."; - } else if (boundaryType == 1) { + secondShape = NoSecondDim; + boundaryFormatAttr = "i32"; + boundaryPrecision = 4; + logInfo(rank) << "Using 32-bit integer boundary type conditions, or 8 bit per face (i32)."; + } else if (boundaryType == 1 || boundaryType == 3) { boundaryDatatype = H5T_STD_I64LE; faceOffset = 16; - logInfo(rank) << "Using 64-bit integer boundary type conditions."; + secondShape = NoSecondDim; + boundaryFormatAttr = "i64"; + boundaryPrecision = 8; + logInfo(rank) << "Using 64-bit integer boundary type conditions, or 16 bit per face (i64)."; + } else if (boundaryType == 4) { + boundaryDatatype = H5T_STD_I32LE; + secondShape = 4; + faceOffset = -1; + boundaryFormatAttr = "i32x4"; + boundaryPrecision = 4; + secondDimBoundary = " 4"; + logInfo(rank) << "Using 32-bit integer per boundary face (i32x4)."; } std::string xdmfFile = outputFile; @@ -306,116 +312,94 @@ int main(int argc, char* argv[]) { } // Create/read the mesh - MeshInput* meshInput = 0L; - apf::Mesh2* mesh = 0L; + MeshData* meshInput = nullptr; switch (args.getArgument("source", 0)) { case 0: logInfo(rank) << "Using Gambit mesh"; - meshInput = new SerialMeshFile(inputFile); + meshInput = new SerialMeshFile(inputFile, faceOffset); break; case 1: - logInfo(rank) << "Using Fidap mesh"; - meshInput = new SerialMeshFile(inputFile); - break; - case 2: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; - meshInput = new SerialMeshFile>(inputFile); + meshInput = + new SerialMeshFile>(inputFile, faceOffset); break; - case 3: + case 2: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; - meshInput = new SerialMeshFile>(inputFile); + meshInput = + new SerialMeshFile>(inputFile, faceOffset); break; - case 4: + case 3: #ifdef USE_NETCDF logInfo(rank) << "Using netCDF mesh"; - meshInput = new NetCDFMesh(inputFile); + meshInput = new NetCDFMesh(inputFile, faceOffset); #else // USE_NETCDF logError() << "netCDF is not supported in this version"; #endif // USE_NETCDF break; - case 5: + case 4: logInfo(rank) << "Using APF native format"; - meshInput = new ApfNative(inputFile, args.getArgument("model", 0L)); +#ifdef USE_SCOREC + meshInput = new ApfNative(inputFile, faceOffset, args.getArgument("input", 0L)); + (dynamic_cast(meshInput))->generate(); +#else + logError() << "This version of PUMgen has been compiled without SCOREC. Hence, the APF format " + "is not available."; +#endif break; - case 6: + case 5: #ifdef USE_SIMMOD logInfo(rank) << "Using SimModSuite"; meshInput = new SimModSuite( - inputFile, args.getArgument("cad", 0L), + inputFile, faceOffset, args.getArgument("cad", 0L), args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), args.getArgument("analysis", "analysis"), args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), args.isSet("analyseAR"), args.getArgument("sim_log", 0L)); #else // USE_SIMMOD - logError() << "SimModSuite is not supported in this version"; + logError() << "SimModSuite is not supported in this version."; #endif // USE_SIMMOD break; + case 6: +#ifdef USE_SIMMOD +#ifdef USE_SCOREC + logInfo(rank) << "Using SimModSuite with APF (deprecated)"; + + meshInput = new SimModSuiteApf( + inputFile, faceOffset, args.getArgument("cad", 0L), + args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), + args.getArgument("analysis", "analysis"), + args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), + args.isSet("analyseAR"), args.getArgument("sim_log", 0L)); + (dynamic_cast(meshInput))->generate(); +#else + logError() << "This version of PUMgen has been compiled without SCOREC. Hence, this reader for " + "the SimModSuite is not available here."; +#endif +#else + logError() << "SimModSuite is not supported in this version."; +#endif + break; default: - logError() << "Unknown source"; + logError() << "Unknown source."; } - mesh = meshInput->getMesh(); - logInfo(rank) << "Parsed mesh successfully, writing output..."; - // Check mesh - if (alignMdsMatches(mesh)) - logWarning() << "Fixed misaligned matches"; - mesh->verify(); - - // Dump mesh for later usage? - const char* dumpFile = args.getArgument("dump", 0L); - if (dumpFile) { - logInfo(PCU_Comm_Self()) << "Writing native APF mesh"; - mesh->writeNative(dumpFile); - - const char* modelFile = args.getArgument("model", 0L); - if (modelFile) - gmi_write_dmg(mesh->getModel(), modelFile); - } - - // Dump VTK mesh - const char* vtkPrefix = args.getArgument("vtk", 0L); - if (vtkPrefix) { - logInfo(PCU_Comm_Self()) << "Writing VTK mesh"; - apf::writeVtkFiles(vtkPrefix, mesh); - } - - apf::Sharing* sharing = apf::getSharing(mesh); - - // oriented at the apf::countOwned method ... But with size_t - auto countOwnedLong = [sharing, mesh](int dim) { - std::size_t counter = 0; - - apf::MeshIterator* it = mesh->begin(dim); - while (apf::MeshEntity* element = mesh->iterate(it)) { - if (sharing->isOwned(element)) { - ++counter; - } - } - mesh->end(it); - - return counter; - }; - - // TODO(David): replace by apf::countOwned again once extended + void* mesh = nullptr; // Get local/global size - std::size_t localSize[2] = {countOwnedLong(3), countOwnedLong(0)}; + std::size_t localSize[2] = {meshInput->cellCount(), meshInput->vertexCount()}; std::size_t globalSize[2] = {localSize[0], localSize[1]}; MPI_Allreduce(MPI_IN_PLACE, globalSize, 2, tndm::mpi_type_t(), MPI_SUM, MPI_COMM_WORLD); - logInfo(rank) << "Mesh size:" << globalSize[0]; + logInfo(rank) << "Total cell count:" << globalSize[0]; + logInfo(rank) << "Total vertex count:" << globalSize[1]; - // Compute min insphere radius - double min = std::numeric_limits::max(); - apf::MeshIterator* it = mesh->begin(3); - while (apf::MeshEntity* element = mesh->iterate(it)) { - min = std::min(min, ma::getInsphere(mesh, element)); - } - mesh->end(it); + auto inspheres = + calculateInsphere(meshInput->connectivity(), meshInput->geometry(), MPI_COMM_WORLD); + double min = *std::min_element(inspheres.begin(), inspheres.end()); MPI_Reduce((rank == 0 ? MPI_IN_PLACE : &min), &min, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD); logInfo(rank) << "Minimum insphere found:" << min; @@ -425,10 +409,6 @@ int main(int argc, char* argv[]) { offsets[0] -= localSize[0]; offsets[1] -= localSize[1]; - // Create numbering for the vertices/elements - apf::GlobalNumbering* vertexNum = apf::makeGlobal(apf::numberOwnedNodes(mesh, "vertices")); - apf::synchronize(vertexNum); - // Create the H5 file hid_t h5falist = H5Pcreate(H5P_FILE_ACCESS); checkH5Err(h5falist); @@ -445,99 +425,41 @@ int main(int argc, char* argv[]) { // Write cells std::size_t connectBytesPerData = 8; logInfo(rank) << "Writing cells"; - writeH5Data( - [&](auto& element, auto&& data) { - apf::NewArray vn; - apf::getElementNumbers(vertexNum, element, vn); - - for (int i = 0; i < 4; i++) { - *(data + i) = vn[i]; - } - return true; - }, - h5file, "connect", mesh, 3, H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->connectivity(), h5file, "connect", mesh, 3, + H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], + globalSize[0], reduceInts, filterEnable, filterChunksize, 4); // Vertices logInfo(rank) << "Writing vertices"; - writeH5Data( - [&](auto& element, auto&& data) { - if (!sharing->isOwned(element)) { - return false; - } - - /*long gid = apf::getNumber(vertexNum, apf::Node(element, 0)); - - if (gid != static_cast(offsets[1] + index)) { - logError() << "Global vertex numbering is incorrect"; - }*/ - - apf::Vector3 point; - mesh->getPoint(element, 0, point); - double geometry[3]; - point.toArray(geometry); - - *(data + 0) = geometry[0]; - *(data + 1) = geometry[1]; - *(data + 2) = geometry[2]; - - return true; - }, - h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, H5T_IEEE_F64LE, chunksize, localSize[1], - globalSize[1], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->geometry(), h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, + H5T_IEEE_F64LE, chunksize, localSize[1], globalSize[1], reduceInts, + filterEnable, filterChunksize, 3); // Group information - apf::MeshTag* groupTag = mesh->findTag("group"); std::size_t groupBytesPerData = 4; - if (groupTag) { - logInfo(rank) << "Writing group information"; - writeH5Data( - [&](auto& element, auto&& data) { - int group; - mesh->getIntTag(element, groupTag, &group); - *data = group; - return true; - }, - h5file, "group", mesh, 3, H5T_NATIVE_INT, H5T_STD_I32LE, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); - } else { - logInfo() << "No group information found in mesh"; - } + logInfo(rank) << "Writing group information"; + writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, H5T_STD_I32LE, + chunksize, localSize[0], globalSize[0], reduceInts, filterEnable, + filterChunksize, NoSecondDim); // Write boundary condition logInfo(rank) << "Writing boundary condition"; - apf::MeshTag* boundaryTag = mesh->findTag("boundary condition"); - assert(boundaryTag); auto i32limit = std::numeric_limits::max(); auto i64limit = std::numeric_limits::max(); - writeH5Data( - [&](auto& element, auto&& data) { - long boundary = 0; - apf::Downward faces; - mesh->getDownward(element, 2, faces); - - for (int i = 0; i < 4; i++) { - if (mesh->hasTag(faces[i], boundaryTag)) { - int b; - mesh->getIntTag(faces[i], boundaryTag, &b); - - if (b <= 0 || (b > i32limit && boundaryType == 0) || - (b > i64limit && boundaryType == 1)) { - logError() << "Cannot handle boundary condition" << b; - } - - // NOTE: this fails for boundary values per face larger than 2**31 - 1 due to sign - // extension - boundary |= static_cast(b) << (i * faceOffset); - } - } - - *data = boundary; - return true; - }, - h5file, "boundary", mesh, 3, H5T_NATIVE_LONG, boundaryDatatype, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->boundary(), h5file, "boundary", mesh, 3, H5T_NATIVE_LONG, + boundaryDatatype, chunksize, localSize[0], globalSize[0], reduceInts, + filterEnable, filterChunksize, secondShape); + + hid_t attrSpace = checkH5Err(H5Screate(H5S_SCALAR)); + hid_t attrType = checkH5Err(H5Tcopy(H5T_C_S1)); + checkH5Err(H5Tset_size(attrType, boundaryFormatAttr.size())); + hid_t attrBoundary = checkH5Err( + H5Acreate(h5file, "boundary-format", attrType, attrSpace, H5P_DEFAULT, H5P_DEFAULT)); + checkH5Err(H5Awrite(attrBoundary, attrType, boundaryFormatAttr.data())); + checkH5Err(H5Aclose(attrBoundary)); + checkH5Err(H5Sclose(attrSpace)); + checkH5Err(H5Tclose(attrType)); // Writing XDMF file if (rank == 0) { @@ -579,9 +501,9 @@ int main(int argc, char* argv[]) { << globalSize[0] << "\">" << basename << ":/group" << std::endl << " " << std::endl << " " << std::endl - << " " << basename - << ":/boundary" << std::endl + << " " + << basename << ":/boundary" << std::endl << " " << std::endl << " " << std::endl << " " << std::endl @@ -590,13 +512,13 @@ int main(int argc, char* argv[]) { checkH5Err(H5Fclose(h5file)); - delete sharing; - delete mesh; delete meshInput; logInfo(rank) << "Finished successfully"; +#ifdef USE_SCOREC PCU_Comm_Free(); +#endif MPI_Finalize(); return 0;