12#include <tclap/CmdLine.h>
15#include <boost/property_tree/ptree.hpp>
16#include <boost/property_tree/xml_parser.hpp>
21#include <range/v3/algorithm/copy.hpp>
22#include <range/v3/algorithm/transform.hpp>
24#include <unordered_set>
52 std::vector<std::unique_ptr<MeshLib::Mesh>>
const& partitioned_meshes,
55 std::vector<MeshEntityMapInfo>
const& merged_node_map,
56 std::vector<MeshEntityMapInfo>
const& merged_element_map)
75 if (pv_name ==
"OGS_VERSION" || pv_name ==
"IntegrationPointMetaData")
77 auto new_pv = MeshLib::getOrCreateMeshProperty<T>(
78 merged_mesh, pv_name, item_type, pv_num_components);
79 new_pv->resize(pv->
size());
81 std::copy(pv->begin(), pv->end(), new_pv->begin());
86 std::vector<MeshLib::PropertyVector<T>*> partition_property_vectors;
87 partition_property_vectors.reserve(partitioned_meshes.size());
88 for (
auto const& mesh : partitioned_meshes)
90 partition_property_vectors.emplace_back(
91 mesh->getProperties().getPropertyVector<T>(pv_name, item_type,
95 auto createNewCellOrNodePropertyVector =
96 [&](std::vector<MeshEntityMapInfo>
const& mesh_entity_map)
98 auto new_pv = MeshLib::getOrCreateMeshProperty<T>(
99 merged_mesh, pv_name, item_type, pv_num_components);
100 std::size_t counter = 0;
101 for (
auto const& entity_info : mesh_entity_map)
103 auto const& partition_pv =
104 partition_property_vectors[entity_info.partition_id];
105 for (
int i_com = 0; i_com < pv_num_components; i_com++)
107 (*new_pv)[counter * pv_num_components + i_com] =
108 (*partition_pv)[entity_info.original_id *
118 createNewCellOrNodePropertyVector(merged_node_map);
124 createNewCellOrNodePropertyVector(merged_element_map);
130 std::vector<std::vector<std::size_t>> partition_element_offsets;
131 partition_element_offsets.reserve(partitioned_meshes.size());
132 for (
auto const& mesh : partitioned_meshes)
134 partition_element_offsets.emplace_back(
136 mesh->getElements(), *pv, properties));
139 auto new_pv = MeshLib::getOrCreateMeshProperty<T>(
140 merged_mesh, pv_name, item_type, pv_num_components);
143 std::size_t counter = 0;
144 auto const ip_meta_data =
147 for (
auto const element : merged_mesh.
getElements())
149 int const number_of_integration_points =
152 counter += number_of_integration_points;
154 new_pv->resize(counter * pv_num_components);
156 auto const global_ip_offsets =
160 std::size_t element_counter = 0;
161 for (
auto const& element_info : merged_element_map)
164 *(partition_property_vectors[element_info.partition_id]);
166 auto const& offsets =
167 partition_element_offsets[element_info.partition_id];
169 int const begin_pos = offsets[element_info.original_id];
170 int const end_pos = offsets[element_info.original_id + 1];
172 std::copy(partition_pv.begin() + begin_pos,
173 partition_pv.begin() + end_pos,
174 new_pv->begin() + global_ip_offsets[element_counter]);
185 std::ifstream ins(pvtu_file_name);
189 OGS_FATAL(
"Could not open pvtu file {:s}.", pvtu_file_name);
192 using boost::property_tree::ptree;
196 auto root = pt.get_child(
"VTKFile");
198 std::vector<std::string> vtu_file_names;
202 for (ptree::value_type
const& v : root.get_child(
"PUnstructuredGrid"))
204 if (v.first ==
"Piece")
209 std::filesystem::path(v.second.get(
"<xmlattr>.Source",
""))
215 if (vtu_file_names.empty())
217 OGS_FATAL(
"PVTU file {:s} does not contain any vtu piece",
221 return vtu_file_names;
225std::tuple<std::vector<MeshLib::Node*>, std::vector<std::size_t>>
228 std::vector<std::size_t> number_of_nodes_per_partition;
229 ranges::transform(meshes, std::back_inserter(number_of_nodes_per_partition),
231 {
return mesh->getNumberOfNodes(); });
232 std::vector<std::size_t> offsets(meshes.size() + 1, 0);
233 std::partial_sum(number_of_nodes_per_partition.begin(),
234 number_of_nodes_per_partition.end(), offsets.begin() + 1);
236 std::vector<MeshLib::Node*> all_nodes;
237 all_nodes.reserve(offsets.back());
238 for (
auto const& mesh : meshes)
240 ranges::copy(mesh->getNodes(), std::back_inserter(all_nodes));
242 return {all_nodes, offsets};
245std::tuple<std::vector<MeshLib::Element*>, std::vector<MeshEntityMapInfo>>
248 std::vector<MeshLib::Element*> regular_elements;
250 std::size_t partition_counter = 0;
251 std::vector<MeshEntityMapInfo> merged_element_map;
252 for (
auto const& mesh : meshes)
256 std::vector<unsigned char>
const& ghost_id_vector =
259 auto const& mesh_elements = mesh->getElements();
261 auto const last_element_id_of_previous_partition =
262 regular_elements.size();
264 std::copy_if(mesh_elements.begin(), mesh_elements.end(),
265 std::back_inserter(regular_elements),
266 [&ghost_id_vector](
auto const element)
267 { return ghost_id_vector[element->getID()] == 0; });
269 for (
auto element_id = last_element_id_of_previous_partition;
270 element_id < regular_elements.size();
273 merged_element_map.push_back(
274 {partition_counter, regular_elements[element_id]->getID()});
280 return {regular_elements, merged_element_map};
287 std::vector<MeshLib::Node*> query_nodes;
289 std::numeric_limits<double>::epsilon() * extent.squaredNorm();
294 if (query_nodes.empty())
297 "query_nodes for node [{}], ({}, {}, {}) of element "
298 "[{}] are empty, eps is {}",
299 node.
getID(), node[0], node[1], node[2], element_id, eps);
302 std::find_if(query_nodes.begin(), query_nodes.end(),
303 [&node, eps](
auto const* p)
305 return (p->asEigenVector3d() - node.asEigenVector3d())
306 .squaredNorm() < eps;
308 if (it == query_nodes.end())
311 "did not find node: [{}] ({}, {}, {}) of element [{}] in "
313 node.
getID(), node[0], node[1], node[2], element_id);
319 std::vector<MeshLib::Element*>
const& regular_elements,
321 Eigen::Vector3d
const& extent)
323 for (
auto& e : regular_elements)
325 for (
unsigned i = 0; i < e->getNumberOfNodes(); i++)
327 auto*
const node = e->getNode(i);
329 if (!oct_tree.
addPoint(node, node_ptr))
332 oct_tree, extent, *node, e->getID());
333 e->setNode(i, node_ptr);
339std::pair<std::vector<MeshLib::Node*>, std::vector<MeshEntityMapInfo>>
341 std::vector<std::size_t>
const& partition_offsets,
344 std::vector<MeshLib::Node*> unique_merged_nodes;
345 unique_merged_nodes.reserve(all_merged_nodes_tmp.size());
347 std::vector<MeshEntityMapInfo> merged_node_map;
348 merged_node_map.reserve(all_merged_nodes_tmp.size());
350 for (std::size_t i = 0; i < partition_offsets.size() - 1; ++i)
352 for (std::size_t pos = partition_offsets[i];
353 pos < partition_offsets[i + 1];
356 auto* node = all_merged_nodes_tmp[pos];
358 if (oct_tree.
addPoint(node, node_ptr))
360 unique_merged_nodes.push_back(node);
361 merged_node_map.push_back({i, pos - partition_offsets[i]});
365 return {unique_merged_nodes, merged_node_map};
368int main(
int argc,
char* argv[])
371 "This tool merges VTU files of PVTU into one single VTU file. Apart "
372 "from the mesh data, all property data are merged as well"
373 "\n\nOpenGeoSys-6 software, version " +
376 "Copyright (c) 2012-2024, OpenGeoSys Community "
377 "(http://www.opengeosys.org)",
380 TCLAP::ValueArg<std::string> output_arg(
381 "o",
"output",
"the output mesh (*.vtu)",
true,
"",
"output.vtu");
384 TCLAP::ValueArg<std::string> input_arg(
385 "i",
"input",
"the partitioned input mesh (*.pvtu)",
true,
"",
388 cmd.parse(argc, argv);
392 OGS_FATAL(
"The extension of input file name {:s} is not \"pvtu\"",
393 input_arg.getValue());
397 OGS_FATAL(
"The extension of output file name {:s} is not \"vtu\"",
398 output_arg.getValue());
403 std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
404 meshes.reserve(vtu_file_names.size());
408 for (
auto const& file_name : vtu_file_names)
410 auto mesh = std::unique_ptr<MeshLib::Mesh>(
418 "Property vector vtkGhostType does not exist in mesh {:s}.",
422 meshes.emplace_back(std::move(mesh));
424 INFO(
"Reading meshes took {} s", io_timer.
elapsed());
427 merged_element_timer.
start();
432 "Collection of {} regular elements and computing element map took {} s",
433 regular_elements.size(), merged_element_timer.
elapsed());
437 collect_nodes_timer.
start();
438 auto [all_merged_nodes_tmp, partition_offsets] =
440 INFO(
"Collection of {} nodes and computing offsets took {} s",
441 all_merged_nodes_tmp.size(), collect_nodes_timer.
elapsed());
444 merged_nodes_timer.
start();
445 GeoLib::AABB aabb(all_merged_nodes_tmp.begin(), all_merged_nodes_tmp.end());
446 auto oct_tree = std::unique_ptr<GeoLib::OctTree<MeshLib::Node, 16>>(
450 auto [unique_merged_nodes, merged_node_map] =
452 INFO(
"Make nodes unique ({} unique nodes) / computing map took {} s",
453 unique_merged_nodes.size(), merged_nodes_timer.
elapsed());
456 reset_nodes_in_elements_timer.
start();
459 INFO(
"Reset nodes in regular elements took {} s",
460 reset_nodes_in_elements_timer.
elapsed());
463 mesh_creation_timer.
start();
467 MeshLib::Mesh(
"pvtu_merged_mesh", unique_merged_nodes, regular_elements,
469 INFO(
"creation of merged mesh took {} s", mesh_creation_timer.
elapsed());
471 auto const& properties = meshes[0]->getProperties();
474 property_timer.
start();
477 [&, &merged_node_map = merged_node_map,
478 &merged_element_map = merged_element_map](
auto type,
479 auto const& property)
481 return createPropertyVector<decltype(type)>(
485 properties, merged_node_map, merged_element_map);
487 INFO(
"merge properties into merged mesh took {} s",
493 writing_timer.
start();
494 auto const result = writer.
writeToFile(output_arg.getValue());
497 ERR(
"Could not write mesh to '{:s}'.", output_arg.getValue());
500 INFO(
"writing mesh took {} s", writing_timer.
elapsed());
Definition of the AABB class.
Definition of the Element class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the class Properties that implements a container of properties.
Definition of the Mesh class.
Definition of the Node class.
Implementation of the OctTree class.
int main(int argc, char *argv[])
std::pair< std::vector< MeshLib::Node * >, std::vector< MeshEntityMapInfo > > makeNodesUnique(std::vector< MeshLib::Node * > const &all_merged_nodes_tmp, std::vector< std::size_t > const &partition_offsets, GeoLib::OctTree< MeshLib::Node, 16 > &oct_tree)
MeshLib::Node * getExistingNodeFromOctTree(GeoLib::OctTree< MeshLib::Node, 16 > &oct_tree, Eigen::Vector3d const &extent, MeshLib::Node const &node, std::size_t const element_id)
void resetNodesInRegularElements(std::vector< MeshLib::Element * > const ®ular_elements, GeoLib::OctTree< MeshLib::Node, 16 > &oct_tree, Eigen::Vector3d const &extent)
std::tuple< std::vector< MeshLib::Element * >, std::vector< MeshEntityMapInfo > > getRegularElements(std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes)
std::vector< std::string > readVtuFileNames(std::string const &pvtu_file_name)
bool createPropertyVector(MeshLib::Mesh &merged_mesh, std::vector< std::unique_ptr< MeshLib::Mesh > > const &partitioned_meshes, MeshLib::PropertyVector< T > const *const pv, MeshLib::Properties const &properties, std::vector< MeshEntityMapInfo > const &merged_node_map, std::vector< MeshEntityMapInfo > const &merged_element_map)
std::tuple< std::vector< MeshLib::Node * >, std::vector< std::size_t > > getMergedNodesVector(std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes)
Definition of the Point3d class.
void applyToPropertyVectors(Properties const &properties, Function f)
Definition of the RunTime class.
Implementation of the VtuInterface class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Eigen::Vector3d const & getMaxPoint() const
Eigen::Vector3d const & getMinPoint() const
bool addPoint(POINT *pnt, POINT *&ret_pnt)
void getPointsInRange(T const &min, T const &max, std::vector< POINT * > &pnts) const
static OctTree< POINT, MAX_POINTS > * createOctTree(Eigen::Vector3d ll, Eigen::Vector3d ur, double eps=std::numeric_limits< double >::epsilon())
std::size_t getID() const
Eigen::Vector3d const & asEigenVector3d() const
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
static MeshLib::Mesh * readVTUFile(std::string const &file_name, bool const compute_element_neighbors=false)
bool writeToFile(std::filesystem::path const &file_path)
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
bool existsPropertyVector(std::string_view name) const
PropertyVector< T > const * getPropertyVector(std::string_view name) const
MeshItemType getMeshItemType() const
int getNumberOfGlobalComponents() const
std::string const & getPropertyName() const
std::string getFileExtension(const std::string &path)
std::string extractPath(std::string const &pathname)
std::string joinPaths(std::string const &pathA, std::string const &pathB)
GITINFOLIB_EXPORT const std::string ogs_version
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Properties const &properties, std::string const &name)