OGS
PVTU2VTU.cpp
Go to the documentation of this file.
1
12#include <tclap/CmdLine.h>
13
14#include <algorithm>
15#include <boost/property_tree/ptree.hpp>
16#include <boost/property_tree/xml_parser.hpp>
17#include <filesystem>
18#include <fstream>
19#include <memory>
20#include <numeric>
21#include <range/v3/algorithm/copy.hpp>
22#include <range/v3/algorithm/transform.hpp>
23#include <string>
24#include <unordered_set>
25#include <vector>
26
27#include "BaseLib/FileTools.h"
28#include "BaseLib/RunTime.h"
29#include "GeoLib/AABB.h"
30#include "GeoLib/OctTree.h"
31#include "InfoLib/GitInfo.h"
32#include "MathLib/Point3d.h"
35#include "MeshLib/Location.h"
36#include "MeshLib/Mesh.h"
37#include "MeshLib/Node.h"
38#include "MeshLib/Properties.h"
42
44{
45 std::size_t partition_id;
46 std::size_t original_id;
47};
48
49template <typename T>
51 MeshLib::Mesh& merged_mesh,
52 std::vector<std::unique_ptr<MeshLib::Mesh>> const& partitioned_meshes,
53 MeshLib::PropertyVector<T> const* const pv,
54 MeshLib::Properties const& properties,
55 std::vector<MeshEntityMapInfo> const& merged_node_map,
56 std::vector<MeshEntityMapInfo> const& merged_element_map)
57{
58 if (pv == nullptr)
59 {
60 return false;
61 }
62
63 if (pv->getPropertyName() == "vtkGhostType")
64 {
65 // Do nothing
66 return true;
67 }
68
69 auto const item_type = pv->getMeshItemType();
70
71 auto const pv_name = pv->getPropertyName();
72
73 auto const pv_num_components = pv->getNumberOfGlobalComponents();
74
75 if (pv_name == "OGS_VERSION" || pv_name == "IntegrationPointMetaData")
76 {
77 auto new_pv = MeshLib::getOrCreateMeshProperty<T>(
78 merged_mesh, pv_name, item_type, pv_num_components);
79 new_pv->resize(pv->size());
80
81 std::copy(pv->begin(), pv->end(), new_pv->begin());
82
83 return true;
84 }
85
86 std::vector<MeshLib::PropertyVector<T>*> partition_property_vectors;
87 partition_property_vectors.reserve(partitioned_meshes.size());
88 for (auto const& mesh : partitioned_meshes)
89 {
90 partition_property_vectors.emplace_back(
91 mesh->getProperties().getPropertyVector<T>(pv_name, item_type,
92 pv_num_components));
93 }
94
95 auto createNewCellOrNodePropertyVector =
96 [&](std::vector<MeshEntityMapInfo> const& mesh_entity_map)
97 {
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)
102 {
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++)
106 {
107 (*new_pv)[counter * pv_num_components + i_com] =
108 (*partition_pv)[entity_info.original_id *
109 pv_num_components +
110 i_com];
111 }
112 counter++;
113 }
114 };
115
116 if (item_type == MeshLib::MeshItemType::Node)
117 {
118 createNewCellOrNodePropertyVector(merged_node_map);
119 return true;
120 }
121
122 if (item_type == MeshLib::MeshItemType::Cell)
123 {
124 createNewCellOrNodePropertyVector(merged_element_map);
125 return true;
126 }
127
129 {
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)
133 {
134 partition_element_offsets.emplace_back(
136 mesh->getElements(), *pv, properties));
137 }
138
139 auto new_pv = MeshLib::getOrCreateMeshProperty<T>(
140 merged_mesh, pv_name, item_type, pv_num_components);
141
142 // Count the integration points
143 std::size_t counter = 0;
144 auto const ip_meta_data =
145 MeshLib::getIntegrationPointMetaData(properties, pv_name);
146
147 for (auto const element : merged_mesh.getElements())
148 {
149 int const number_of_integration_points =
151 *element);
152 counter += number_of_integration_points;
153 }
154 new_pv->resize(counter * pv_num_components);
155
156 auto const global_ip_offsets =
158 merged_mesh.getElements(), *pv, properties);
159
160 std::size_t element_counter = 0;
161 for (auto const& element_info : merged_element_map)
162 {
163 MeshLib::PropertyVector<T> const& partition_pv =
164 *(partition_property_vectors[element_info.partition_id]);
165
166 auto const& offsets =
167 partition_element_offsets[element_info.partition_id];
168
169 int const begin_pos = offsets[element_info.original_id];
170 int const end_pos = offsets[element_info.original_id + 1];
171
172 std::copy(partition_pv.begin() + begin_pos,
173 partition_pv.begin() + end_pos,
174 new_pv->begin() + global_ip_offsets[element_counter]);
175
176 element_counter++;
177 }
178 }
179
180 return true;
181}
182
183std::vector<std::string> readVtuFileNames(std::string const& pvtu_file_name)
184{
185 std::ifstream ins(pvtu_file_name);
186
187 if (!ins)
188 {
189 OGS_FATAL("Could not open pvtu file {:s}.", pvtu_file_name);
190 }
191
192 using boost::property_tree::ptree;
193 ptree pt;
194 read_xml(ins, pt);
195
196 auto root = pt.get_child("VTKFile");
197
198 std::vector<std::string> vtu_file_names;
199
200 std::string file_path = BaseLib::extractPath(pvtu_file_name);
201
202 for (ptree::value_type const& v : root.get_child("PUnstructuredGrid"))
203 {
204 if (v.first == "Piece")
205 {
206 vtu_file_names.push_back(BaseLib::joinPaths(
207 file_path,
208 // only gets the vtu file name:
209 std::filesystem::path(v.second.get("<xmlattr>.Source", ""))
210 .filename()
211 .string()));
212 }
213 }
214
215 if (vtu_file_names.empty())
216 {
217 OGS_FATAL("PVTU file {:s} does not contain any vtu piece",
218 pvtu_file_name);
219 }
220
221 return vtu_file_names;
222}
223
224// all nodes (also 'ghost' nodes) of all meshes sorted by partition
225std::tuple<std::vector<MeshLib::Node*>, std::vector<std::size_t>>
226getMergedNodesVector(std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
227{
228 std::vector<std::size_t> number_of_nodes_per_partition;
229 ranges::transform(meshes, std::back_inserter(number_of_nodes_per_partition),
230 [](auto const& mesh)
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);
235
236 std::vector<MeshLib::Node*> all_nodes;
237 all_nodes.reserve(offsets.back());
238 for (auto const& mesh : meshes)
239 {
240 ranges::copy(mesh->getNodes(), std::back_inserter(all_nodes));
241 }
242 return {all_nodes, offsets};
243}
244
245std::tuple<std::vector<MeshLib::Element*>, std::vector<MeshEntityMapInfo>>
246getRegularElements(std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
247{
248 std::vector<MeshLib::Element*> regular_elements;
249
250 std::size_t partition_counter = 0;
251 std::vector<MeshEntityMapInfo> merged_element_map;
252 for (auto const& mesh : meshes)
253 {
254 MeshLib::Properties const& properties = mesh->getProperties();
255
256 std::vector<unsigned char> const& ghost_id_vector =
257 *(properties.getPropertyVector<unsigned char>("vtkGhostType"));
258
259 auto const& mesh_elements = mesh->getElements();
260
261 auto const last_element_id_of_previous_partition =
262 regular_elements.size();
263
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; });
268
269 for (auto element_id = last_element_id_of_previous_partition;
270 element_id < regular_elements.size();
271 element_id++)
272 {
273 merged_element_map.push_back(
274 {partition_counter, regular_elements[element_id]->getID()});
275 }
276
277 partition_counter++;
278 }
279
280 return {regular_elements, merged_element_map};
281}
282
284 GeoLib::OctTree<MeshLib::Node, 16>& oct_tree, Eigen::Vector3d const& extent,
285 MeshLib::Node const& node, std::size_t const element_id)
286{
287 std::vector<MeshLib::Node*> query_nodes;
288 auto const eps =
289 std::numeric_limits<double>::epsilon() * extent.squaredNorm();
290 Eigen::Vector3d const min = node.asEigenVector3d().array() - eps;
291 Eigen::Vector3d const max = node.asEigenVector3d().array() + eps;
292 oct_tree.getPointsInRange(min, max, query_nodes);
293
294 if (query_nodes.empty())
295 {
296 OGS_FATAL(
297 "query_nodes for node [{}], ({}, {}, {}) of element "
298 "[{}] are empty, eps is {}",
299 node.getID(), node[0], node[1], node[2], element_id, eps);
300 }
301 auto const it =
302 std::find_if(query_nodes.begin(), query_nodes.end(),
303 [&node, eps](auto const* p)
304 {
305 return (p->asEigenVector3d() - node.asEigenVector3d())
306 .squaredNorm() < eps;
307 });
308 if (it == query_nodes.end())
309 {
310 OGS_FATAL(
311 "did not find node: [{}] ({}, {}, {}) of element [{}] in "
312 "query_nodes",
313 node.getID(), node[0], node[1], node[2], element_id);
314 }
315 return *it;
316}
317
319 std::vector<MeshLib::Element*> const& regular_elements,
321 Eigen::Vector3d const& extent)
322{
323 for (auto& e : regular_elements)
324 {
325 for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
326 {
327 auto* const node = e->getNode(i);
328 MeshLib::Node* node_ptr = nullptr;
329 if (!oct_tree.addPoint(node, node_ptr))
330 {
331 auto const node_ptr = getExistingNodeFromOctTree(
332 oct_tree, extent, *node, e->getID());
333 e->setNode(i, node_ptr);
334 }
335 }
336 }
337}
338
339std::pair<std::vector<MeshLib::Node*>, std::vector<MeshEntityMapInfo>>
340makeNodesUnique(std::vector<MeshLib::Node*> const& all_merged_nodes_tmp,
341 std::vector<std::size_t> const& partition_offsets,
343{
344 std::vector<MeshLib::Node*> unique_merged_nodes;
345 unique_merged_nodes.reserve(all_merged_nodes_tmp.size());
346
347 std::vector<MeshEntityMapInfo> merged_node_map;
348 merged_node_map.reserve(all_merged_nodes_tmp.size());
349
350 for (std::size_t i = 0; i < partition_offsets.size() - 1; ++i)
351 {
352 for (std::size_t pos = partition_offsets[i];
353 pos < partition_offsets[i + 1];
354 ++pos)
355 {
356 auto* node = all_merged_nodes_tmp[pos];
357 MeshLib::Node* node_ptr = nullptr;
358 if (oct_tree.addPoint(node, node_ptr))
359 {
360 unique_merged_nodes.push_back(node);
361 merged_node_map.push_back({i, pos - partition_offsets[i]});
362 }
363 }
364 }
365 return {unique_merged_nodes, merged_node_map};
366}
367
368int main(int argc, char* argv[])
369{
370 TCLAP::CmdLine cmd(
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 " +
375 ".\n"
376 "Copyright (c) 2012-2024, OpenGeoSys Community "
377 "(http://www.opengeosys.org)",
379
380 TCLAP::ValueArg<std::string> output_arg(
381 "o", "output", "the output mesh (*.vtu)", true, "", "output.vtu");
382 cmd.add(output_arg);
383
384 TCLAP::ValueArg<std::string> input_arg(
385 "i", "input", "the partitioned input mesh (*.pvtu)", true, "",
386 "input.pvtu");
387 cmd.add(input_arg);
388 cmd.parse(argc, argv);
389
390 if (BaseLib::getFileExtension(input_arg.getValue()) != ".pvtu")
391 {
392 OGS_FATAL("The extension of input file name {:s} is not \"pvtu\"",
393 input_arg.getValue());
394 }
395 if (BaseLib::getFileExtension(output_arg.getValue()) != ".vtu")
396 {
397 OGS_FATAL("The extension of output file name {:s} is not \"vtu\"",
398 output_arg.getValue());
399 }
400
401 auto const vtu_file_names = readVtuFileNames(input_arg.getValue());
402
403 std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
404 meshes.reserve(vtu_file_names.size());
405
406 BaseLib::RunTime io_timer;
407 io_timer.start();
408 for (auto const& file_name : vtu_file_names)
409 {
410 auto mesh = std::unique_ptr<MeshLib::Mesh>(
412
413 MeshLib::Properties const& properties = mesh->getProperties();
414
415 if (!properties.existsPropertyVector<unsigned char>("vtkGhostType"))
416 {
417 OGS_FATAL(
418 "Property vector vtkGhostType does not exist in mesh {:s}.",
419 file_name);
420 }
421
422 meshes.emplace_back(std::move(mesh));
423 }
424 INFO("Reading meshes took {} s", io_timer.elapsed());
425
426 BaseLib::RunTime merged_element_timer;
427 merged_element_timer.start();
428 // If structured binding is used for the returned tuple, Mac compiler gives
429 // an error in reference to local binding in calling applyToPropertyVectors.
430 auto [regular_elements, merged_element_map] = getRegularElements(meshes);
431 INFO(
432 "Collection of {} regular elements and computing element map took {} s",
433 regular_elements.size(), merged_element_timer.elapsed());
434
435 // alternative implementation of getNodesOfRegularElements
436 BaseLib::RunTime collect_nodes_timer;
437 collect_nodes_timer.start();
438 auto [all_merged_nodes_tmp, partition_offsets] =
439 getMergedNodesVector(meshes);
440 INFO("Collection of {} nodes and computing offsets took {} s",
441 all_merged_nodes_tmp.size(), collect_nodes_timer.elapsed());
442
443 BaseLib::RunTime merged_nodes_timer;
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>>(
448 aabb.getMinPoint(), aabb.getMaxPoint(), 1e-16));
449
450 auto [unique_merged_nodes, merged_node_map] =
451 makeNodesUnique(all_merged_nodes_tmp, partition_offsets, *oct_tree);
452 INFO("Make nodes unique ({} unique nodes) / computing map took {} s",
453 unique_merged_nodes.size(), merged_nodes_timer.elapsed());
454
455 BaseLib::RunTime reset_nodes_in_elements_timer;
456 reset_nodes_in_elements_timer.start();
457 auto const extent = aabb.getMaxPoint() - aabb.getMinPoint();
458 resetNodesInRegularElements(regular_elements, *oct_tree, extent);
459 INFO("Reset nodes in regular elements took {} s",
460 reset_nodes_in_elements_timer.elapsed());
461
462 BaseLib::RunTime mesh_creation_timer;
463 mesh_creation_timer.start();
464 // The Node pointers of 'merged_nodes' and Element pointers of
465 // 'regular_elements' are shared with 'meshes', the partitioned meshes.
466 MeshLib::Mesh merged_mesh =
467 MeshLib::Mesh("pvtu_merged_mesh", unique_merged_nodes, regular_elements,
468 false /* compute_element_neighbors */);
469 INFO("creation of merged mesh took {} s", mesh_creation_timer.elapsed());
470
471 auto const& properties = meshes[0]->getProperties();
472
473 BaseLib::RunTime property_timer;
474 property_timer.start();
476 properties,
477 [&, &merged_node_map = merged_node_map,
478 &merged_element_map = merged_element_map](auto type,
479 auto const& property)
480 {
481 return createPropertyVector<decltype(type)>(
482 merged_mesh, meshes,
483 dynamic_cast<MeshLib::PropertyVector<decltype(type)> const*>(
484 property),
485 properties, merged_node_map, merged_element_map);
486 });
487 INFO("merge properties into merged mesh took {} s",
488 property_timer.elapsed());
489
490 MeshLib::IO::VtuInterface writer(&merged_mesh);
491
492 BaseLib::RunTime writing_timer;
493 writing_timer.start();
494 auto const result = writer.writeToFile(output_arg.getValue());
495 if (!result)
496 {
497 ERR("Could not write mesh to '{:s}'.", output_arg.getValue());
498 return EXIT_FAILURE;
499 }
500 INFO("writing mesh took {} s", writing_timer.elapsed());
501
502 // Since the Node pointers of 'merged_nodes' and Element pointers of
503 // 'regular_elements' are held by 'meshes', the partitioned meshes, the
504 // memory by these pointers are released by 'meshes' automatically.
505 // Therefore, only node vector and element vector of merged_mesh should be
506 // cleaned.
507 merged_mesh.shallowClean();
508
509 return EXIT_SUCCESS;
510}
Definition of the AABB class.
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
Filename manipulation routines.
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
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[])
Definition PVTU2VTU.cpp:368
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)
Definition PVTU2VTU.cpp:340
MeshLib::Node * getExistingNodeFromOctTree(GeoLib::OctTree< MeshLib::Node, 16 > &oct_tree, Eigen::Vector3d const &extent, MeshLib::Node const &node, std::size_t const element_id)
Definition PVTU2VTU.cpp:283
void resetNodesInRegularElements(std::vector< MeshLib::Element * > const &regular_elements, GeoLib::OctTree< MeshLib::Node, 16 > &oct_tree, Eigen::Vector3d const &extent)
Definition PVTU2VTU.cpp:318
std::tuple< std::vector< MeshLib::Element * >, std::vector< MeshEntityMapInfo > > getRegularElements(std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes)
Definition PVTU2VTU.cpp:246
std::vector< std::string > readVtuFileNames(std::string const &pvtu_file_name)
Definition PVTU2VTU.cpp:183
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)
Definition PVTU2VTU.cpp:50
std::tuple< std::vector< MeshLib::Node * >, std::vector< std::size_t > > getMergedNodesVector(std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes)
Definition PVTU2VTU.cpp:226
Definition of the Point3d class.
void applyToPropertyVectors(Properties const &properties, Function f)
Definition of the RunTime class.
Implementation of the VtuInterface class.
Count the running time.
Definition RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition RunTime.h:42
void start()
Start the timer.
Definition RunTime.h:32
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
Eigen::Vector3d const & getMaxPoint() const
Definition AABB.h:187
Eigen::Vector3d const & getMinPoint() const
Definition AABB.h:180
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
Definition Point3d.h:64
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.
Definition Mesh.h:109
void shallowClean()
Definition Mesh.cpp:133
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition Properties.h:36
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::size_t size() 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)
std::vector< std::size_t > getIntegrationPointDataOffsetsOfMeshElements(std::vector< MeshLib::Element * > const &mesh_elements, MeshLib::PropertyVectorBase const &pv, MeshLib::Properties const &properties)
int getNumberOfElementIntegrationPoints(MeshLib::IntegrationPointMetaData const &ip_meta_data, MeshLib::Element const &e)
std::size_t partition_id
Definition PVTU2VTU.cpp:45
std::size_t original_id
Definition PVTU2VTU.cpp:46