diff --git a/CMakeLists.txt b/CMakeLists.txt index f969e0c6..c37da694 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,6 +44,8 @@ endif() find_package(MPI REQUIRED) +option(ASTE_USE_ASTE_JOINER_CPP "Use C++ implementation of precice-aste-join" OFF) + add_executable(precice-aste-run src/precice-aste-run.cpp src/common.cpp src/mesh.cpp src/configreader.cpp src/modes.cpp src/utilities.cpp src/logger.cpp) target_include_directories(precice-aste-run PRIVATE src thirdparty) target_link_libraries(precice-aste-run @@ -62,6 +64,14 @@ if(METIS_FOUND) target_link_libraries(precice-aste-run metisAPI) endif() +if(ASTE_USE_ASTE_JOINER_CPP) +add_executable(precice-aste-join src/precice-aste-join.cpp) +target_include_directories(precice-aste-join PRIVATE src thirdparty) +target_link_libraries(precice-aste-join + ${Boost_LIBRARIES} + ${VTK_LIBRARIES} +) +endif() add_executable(testing tests/testing.cpp tests/read_test.cpp tests/write_test.cpp src/mesh.cpp src/logger.cpp) target_include_directories(testing PRIVATE src thirdparty) @@ -77,13 +87,15 @@ if (VTK_VERSION VERSION_LESS "8.90.0") else () # vtk_module_autoinit is needed vtk_module_autoinit( - TARGETS precice-aste-run testing + TARGETS precice-aste-run testing precice-aste-join MODULES ${VTK_LIBRARIES} ) endif() file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-partition DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +if(NOT ASTE_USE_ASTE_JOINER_CPP) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-join DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +endif() file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-evaluate DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) include(GNUInstallDirs) @@ -91,7 +103,12 @@ install(TARGETS precice-aste-run DESTINATION ${CMAKE_INSTALL_BINDIR}) if(METIS_FOUND) install(TARGETS metisAPI DESTINATION ${CMAKE_INSTALL_LIBDIR}) endif() +if(ASTE_USE_ASTE_JOINER_CPP) +install(TARGETS precice-aste-join DESTINATION ${CMAKE_INSTALL_BINDIR}) +install(PROGRAMS ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-partition ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-evaluate DESTINATION ${CMAKE_INSTALL_BINDIR}) +else() install(PROGRAMS ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-partition ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-join ${CMAKE_CURRENT_SOURCE_DIR}/src/precice-aste-evaluate DESTINATION ${CMAKE_INSTALL_BINDIR}) +endif() enable_testing() diff --git a/src/precice-aste-join.cpp b/src/precice-aste-join.cpp new file mode 100644 index 00000000..3fa0a4c8 --- /dev/null +++ b/src/precice-aste-join.cpp @@ -0,0 +1,323 @@ +#include "precice-aste-join.hpp" + +auto getOptions(int argc, char *argv[]) -> OptionMap +{ + namespace po = boost::program_options; + + // Declare the supported options. + po::options_description desc("Allowed options"); + desc.add_options()("help", "produce help message")("mesh,m", po::value(), "The partitioned mesh prefix used as input (only VTU format is accepted)(Looking for _<#filerank>.vtu).")("output,o", po::value(), "The output mesh. Can be VTK or VTU format. If it is not given _joined.vtk will be used.")("recovery,r", po::value(), "The path to the recovery file to fully recover it's state.")("numparts,n", po::value()->default_value(0), "The number of parts to read from the input mesh. By default the entire mesh is read.")("directory,dir", po::value()->default_value("."), "Directory for output files (optional)"); + + po::variables_map vm; + try { + po::store(parse_command_line(argc, argv, desc), vm); + + if (vm.count("help")) { + std::cout << desc << std::endl; + std::exit(EXIT_SUCCESS); + } + // Needs to be called + po::notify(vm); + + if (!vm.count("mesh")) { + std::cout << "You must specify a mesh file to read from." << std::endl; + std::exit(EXIT_SUCCESS); + } + + } catch (const std::exception &e) { + std::cerr << "ERROR: " << e.what() << "\n"; + std::cerr << desc << std::endl; + std::exit(EXIT_FAILURE); + } + return vm; +} + +void readRecoveryFile(const std::string &recoveryFile, int &size, std::vector &cellTypes, std::vector> &cells) +{ + // Parse the file + std::ifstream ifs(recoveryFile); + json recoveryData = json::parse(ifs); + // Get content + try { + size = recoveryData["size"].get(); + } catch (nlohmann::detail::parse_error &) { + std::cerr << "Error while parsing the recovery file: \"size\" attribute is missing"; + std::exit(EXIT_FAILURE); + } catch (nlohmann::detail::type_error &) { + std::cerr << "Error while parsing recovery file \"size\" is missing"; + std::exit(EXIT_FAILURE); + } + + try { + cellTypes = recoveryData["cell_types"].get>(); + } catch (nlohmann::detail::parse_error &) { + std::cerr << "Error while parsing the recovery file: \"cell_types\" attribute is missing"; + std::exit(EXIT_FAILURE); + } catch (nlohmann::detail::type_error &) { + std::cerr << "Error while parsing recovery file \"cell_types\" is missing"; + std::exit(EXIT_FAILURE); + } + + try { + cells = recoveryData["cells"].get>>(); + } catch (nlohmann::detail::parse_error &) { + std::cerr << "Error while parsing recovery file \"cells\" is missing"; + std::exit(EXIT_FAILURE); + } catch (nlohmann::detail::type_error &) { + std::cerr << "Error while parsing recovery file \"cells\" is missing"; + std::exit(EXIT_FAILURE); + } +} + +auto countPartitions(const std::string &prefix) -> size_t +{ + namespace fs = boost::filesystem; + + std::string filename; + size_t count = 0; + + while (true) { + filename = prefix + "_" + std::to_string(count) + ".vtu"; + if (!fs::exists(filename)) { + break; + } + ++count; + } + return count; +} + +void writeMesh(const std::string &filename, const std::string &directory, vtkSmartPointer mesh) +{ + namespace fs = boost::filesystem; + if (fs::exists(directory)) { + if (!fs::is_directory(directory)) { + std::cerr << "Error: " << directory << " is not a directory." << std::endl; + std::exit(EXIT_FAILURE); + } + } else { + if (!fs::create_directory(directory)) { + std::cerr << "Error: Could not create directory " << directory << std::endl; + std::exit(EXIT_FAILURE); + } + } + + auto output_path = fs::current_path() / fs::path(directory) / fs::path(filename); + auto output_filename = output_path.c_str(); + + if (fs::extension(output_filename) == ".vtu") { + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetFileName(output_filename); + writer->SetInputData(mesh); + writer->Write(); + } else if (fs::extension(filename) == ".vtk") { + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetFileName(output_filename); + writer->SetInputData(mesh); + writer->SetFileTypeToBinary(); + writer->Write(); + } else { + std::cerr << "Error: " << filename << " is not a valid output file." << std::endl; + std::exit(EXIT_FAILURE); + } +} + +auto partitionwiseMerge(const std::string &prefix, size_t numparts) -> vtkSmartPointer +{ + auto joinedMesh = vtkSmartPointer::New(); + auto joinedPoints = vtkSmartPointer::New(); + joinedPoints->SetDataTypeToDouble(); + auto joinedCells = vtkSmartPointer::New(); + std::vector joinedCellTypes; + + std::vector> joinedDataVec; + std::vector joinedDatanames; + + auto reader = vtkSmartPointer::New(); + for (size_t i = 0; i < numparts; ++i) { + // Read mesh + auto partname = prefix + "_" + std::to_string(i) + ".vtu"; + reader->SetFileName(partname.c_str()); + reader->Update(); + // Extract mesh + auto grid = reader->GetOutput(); + // Cells + const auto offset = joinedPoints->GetNumberOfPoints(); + auto numCells = grid->GetCells()->GetNumberOfCells(); + joinedCellTypes.reserve(joinedCellTypes.size() + numCells); + auto cellIds = vtkSmartPointer::New(); + for (int j = 0; j < numCells; ++j) { + cellIds->Reset(); + std::for_each(grid->GetCell(j)->GetPointIds()->begin(), grid->GetCell(j)->GetPointIds()->end(), [&cellIds, &offset](auto &pointId) { cellIds->InsertNextId(pointId + offset); }); + joinedCells->InsertNextCell(cellIds); + joinedCellTypes.push_back(grid->GetCellType(j)); + } + // Points + auto points = grid->GetPoints(); + joinedPoints->InsertPoints(joinedPoints->GetNumberOfPoints(), grid->GetNumberOfPoints(), 0, points); + // Point Data + auto partPointData = grid->GetPointData(); + auto numArrays = partPointData->GetNumberOfArrays(); + for (int j = 0; j < numArrays; ++j) { + auto partData = partPointData->GetArray(j); + auto name = partData->GetName(); + if (std::find(joinedDatanames.begin(), joinedDatanames.end(), name) == joinedDatanames.end()) { + joinedDatanames.emplace_back(name); + auto newJoinedData = vtkSmartPointer::New(); + newJoinedData->SetName(name); + newJoinedData->SetNumberOfComponents(partData->GetNumberOfComponents()); + joinedDataVec.push_back(newJoinedData); + } + auto joinedData = joinedDataVec[j]; + joinedData->InsertTuples(joinedData->GetNumberOfTuples(), partData->GetNumberOfTuples(), 0, partData); + } + } + + joinedMesh->SetPoints(joinedPoints); + for (const auto &data : joinedDataVec) { + joinedMesh->GetPointData()->AddArray(data); + } + joinedMesh->SetCells(joinedCellTypes.data(), joinedCells); + + return joinedMesh; +} + +auto recoveryMerge(const std::string &prefix, std::size_t numparts, int size, const std::vector &cellTypes, const std::vector> &cells) -> vtkSmartPointer +{ + auto joinedMesh = vtkSmartPointer::New(); + auto joinedPoints = vtkSmartPointer::New(); + joinedPoints->SetDataTypeToDouble(); + auto joinedCells = vtkSmartPointer::New(); + std::vector joinedCellTypes; + + joinedPoints->SetNumberOfPoints(size); + + std::vector> joinedDataVec; + std::vector joinedDataNames; + + auto reader = vtkSmartPointer::New(); + auto globalIds = vtkSmartPointer::New(); + auto localIds = vtkSmartPointer::New(); + + for (size_t i = 0; i < numparts; ++i) { + globalIds->Reset(); + // Read mesh + auto partname = prefix + "_" + std::to_string(i) + ".vtu"; + reader->SetFileName(partname.c_str()); + reader->Update(); + // Extract mesh + auto grid = reader->GetOutput(); + // Points + auto points = grid->GetPoints(); + // Set local Ids + localIds->SetNumberOfIds(points->GetNumberOfPoints()); + std::iota(localIds->begin(), localIds->end(), 0); + // Extract Global Ids + auto partPointData = grid->GetPointData(); + auto globalIdsArray = partPointData->GetArray("GlobalIDs"); + if (globalIdsArray == nullptr) { + std::cerr << "GlobalIDs not found in " << partname << std::endl; + std::cout << " Fall back to partitionwise merge" << std::endl; + return partitionwiseMerge(prefix, numparts); + } else { + globalIds->Allocate(globalIdsArray->GetNumberOfTuples()); + for (vtkIdType j = 0; j < globalIdsArray->GetNumberOfTuples(); ++j) { + globalIds->InsertNextId(static_cast(globalIdsArray->GetTuple1(j))); + } + joinedPoints->InsertPoints(globalIds, localIds, points); + } + + // Cells + auto numCells = grid->GetCells()->GetNumberOfCells(); + joinedCellTypes.reserve(joinedCellTypes.size() + numCells); + auto cellIds = vtkSmartPointer::New(); + for (vtkIdType j = 0; j < numCells; ++j) { + cellIds->Reset(); + std::for_each(grid->GetCell(j)->GetPointIds()->begin(), grid->GetCell(j)->GetPointIds()->end(), [&cellIds, &globalIds](auto &localPointId) { cellIds->InsertNextId(globalIds->GetId(localPointId)); }); + joinedCells->InsertNextCell(cellIds); + joinedCellTypes.push_back(grid->GetCellType(j)); + } + + // Point Data + auto numArrays = partPointData->GetNumberOfArrays(); + for (int j = 0; j < numArrays; ++j) { + auto partData = partPointData->GetArray(j); + auto name = partData->GetName(); + if (std::find(joinedDataNames.begin(), joinedDataNames.end(), name) == joinedDataNames.end()) { + joinedDataNames.emplace_back(name); + auto newJoinedData = vtkSmartPointer::New(); + newJoinedData->SetName(name); + newJoinedData->SetNumberOfComponents(partData->GetNumberOfComponents()); + newJoinedData->Allocate(size); + joinedDataVec.push_back(newJoinedData); + } + auto joinedData = joinedDataVec[j]; + joinedData->InsertTuples(globalIds, localIds, partData); + } + } + + // Add Recovery cells + auto numCells = cells.size(); + joinedCellTypes.reserve(joinedCellTypes.size() + numCells); + auto cellIds = vtkSmartPointer::New(); + for (std::size_t i = 0; i < numCells; ++i) { + cellIds->Reset(); + std::for_each(cells[i].begin(), cells[i].end(), [&cellIds](auto &pointId) { cellIds->InsertNextId(pointId); }); + joinedCells->InsertNextCell(cellIds); + joinedCellTypes.push_back(cellTypes[i]); + } + + // Assemble final mesh + for (const auto &data : joinedDataVec) { + joinedMesh->GetPointData()->AddArray(data); + } + joinedMesh->SetPoints(joinedPoints); + joinedMesh->SetCells(joinedCellTypes.data(), joinedCells); + + return joinedMesh; +} + +void join(int argc, char *argv[]) +{ + namespace fs = boost::filesystem; + auto options = getOptions(argc, argv); + std::string prefix = options["mesh"].as(); + std::string output{}; + std::string recovery{}; + if (options.find("output") == options.end()) { + output = prefix + "_joined.vtk"; + } else { + output = options["output"].as(); + } + if (options.find("recovery") != options.end()) { + recovery = options["recovery"].as(); + } else { + recovery = prefix + "_recovery.json"; + } + std::string directory = options["directory"].as(); + size_t numparts = options["numparts"].as(); + + if (numparts == 0) { + numparts = countPartitions(prefix); + } + + vtkSmartPointer joinedMesh = nullptr; + if (fs::exists(recovery)) { + std::cout << "Recovery file found. Recovering the connectivity information of the mesh." << std::endl; + int size; + std::vector cellTypes; + std::vector> cells; + readRecoveryFile(recovery, size, cellTypes, cells); + joinedMesh = recoveryMerge(prefix, numparts, size, cellTypes, cells); + } else { + std::cout << "Recovery file not found. Partition-wise merging will be done." << std::endl; + joinedMesh = partitionwiseMerge(prefix, numparts); + } + + writeMesh(output, options["directory"].as(), joinedMesh); +} + +auto main(int argc, char *argv[]) -> int +{ + join(argc, argv); + return EXIT_SUCCESS; +} diff --git a/src/precice-aste-join.hpp b/src/precice-aste-join.hpp new file mode 100644 index 00000000..0a21b0eb --- /dev/null +++ b/src/precice-aste-join.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "json.hpp" + +using OptionMap = boost::program_options::variables_map; +using json = nlohmann::json; + +/** + * @brief Get the Options object from command line arguments + * + * @param argc + * @param argv + * @return OptionMap + */ +auto getOptions(int argc, char *argv[]) -> OptionMap; + +void readRecoveryFile(const std::string &recoveryFile, int &size, std::vector &cellTypes, std::vector> &cells); + +/** + * @brief Count the number of partitioned mesh files for given prefix + * + * @param prefix + * @return size_t + */ +auto countPartitions(const std::string &prefix) -> size_t; + +/** + * @brief Write the joined mesh to a VTK file + * + * @param filename + * @param directory + * @param mesh + */ +void writeMesh(const std::string &filename, const std::string &directory, vtkSmartPointer mesh); + +/** + * @brief Merge the meshes from the partitioned files + * @details The meshes are merged in the order of the partitioned files + * @details The cells betwen the meshes are not connected + * @details The point numbering is not preserved between unpartitioned mesh and the merged mesh + * + * @param prefix + * @param numparts + * @return vtkSmartPointer + */ +auto partitionwiseMerge(const std::string &prefix, size_t numparts) -> vtkSmartPointer; + +/** + * @brief Merge the meshes from the partitioned files + * @details The meshes are merged to recover the original mesh + * @details All the cells are preserved and connected + * @details The point numbering is preserved between unpartitioned mesh and the merged mesh + * + * @param prefix + * @param numparts + * @param size + * @param cellTypes + * @param cells + * @return vtkSmartPointer + */ +auto recoveryMerge(const std::string &prefix, std::size_t numparts, int size, const std::vector &cellTypes, const std::vector> &cells) -> vtkSmartPointer; + +/** + * @brief Read the commandline arguments and merge the meshes to a single mesh + */ +void join(int argc, char *argv[]);