OGS
VerticalSliceFromLayers.cpp
Go to the documentation of this file.
1
9
10#include <tclap/CmdLine.h>
11
12#include <QCoreApplication>
13#include <algorithm>
14#include <cmath>
15#include <memory>
16#include <string>
17#include <vector>
18
21#include "BaseLib/FileTools.h"
23#include "BaseLib/MPI.h"
24#include "GeoLib/AABB.h"
26#include "GeoLib/GEOObjects.h"
28#include "GeoLib/Point.h"
29#include "GeoLib/Polygon.h"
30#include "GeoLib/Polyline.h"
31#include "InfoLib/GitInfo.h"
32#include "MathLib/MathTools.h"
33#include "MathLib/Point3d.h"
38#include "MeshLib/Mesh.h"
40#include "MeshLib/Node.h"
44
46std::vector<GeoLib::Point*> createPoints(MathLib::Point3d const& start,
47 MathLib::Point3d const& end,
48 std::size_t const n_intervals)
49{
50 std::vector<GeoLib::Point*> points;
51 points.push_back(new GeoLib::Point(start, 0));
52 double const length = std::sqrt(MathLib::sqrDist(start, end));
53 double const interval = length / n_intervals;
54
55 GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
56 (end[2] - start[2]));
57 for (std::size_t i = 1; i < n_intervals; ++i)
58 {
59 points.push_back(new GeoLib::Point(
60 start[0] + ((i * interval) / length * vec[0]),
61 start[1] + ((i * interval) / length * vec[1]), 0, i));
62 }
63 points.push_back(new GeoLib::Point(end, n_intervals));
64 return points;
65}
66
68GeoLib::Polyline* createPolyline(std::vector<GeoLib::Point*> const& points)
69{
70 GeoLib::Polyline* line = new GeoLib::Polyline(points);
71 std::size_t const length = points.size();
72 for (std::size_t i = 0; i < length; ++i)
73 {
74 line->addPoint(i);
75 }
76 return line;
77}
78
80std::vector<std::string> createGeometries(
82 std::vector<std::string> const& layer_names,
83 MathLib::Point3d const& pnt_start,
84 MathLib::Point3d const& pnt_end,
85 double const resolution)
86{
87 std::vector<std::string> geo_name_list;
88 std::size_t const n_layers = layer_names.size();
89 for (std::size_t i = 0; i < n_layers; ++i)
90 {
91 std::unique_ptr<MeshLib::Mesh> const layer(
92 MeshLib::IO::readMeshFromFile(layer_names[i]));
93 if (layer == nullptr)
94 {
95 ERR("Could not read file {:s}. Skipping layer...", layer_names[i]);
96 continue;
97 }
98 if (layer->getDimension() != 2)
99 {
100 ERR("Layer {:d} is not a 2D mesh. Skipping layer...", i);
101 continue;
102 }
103
104 std::string geo_name(std::to_string(i));
105 auto points(createPoints(pnt_start, pnt_end, resolution));
106 geo.addPointVec(std::move(points), geo_name,
108
109 std::vector<GeoLib::Polyline*> lines;
110 lines.push_back(createPolyline(*geo.getPointVec(geo_name)));
111 geo.addPolylineVec(std::move(lines), geo_name,
113
114 MeshGeoToolsLib::GeoMapper mapper(geo, geo_name);
115 mapper.mapOnMesh(layer.get());
116 geo_name_list.push_back(geo_name);
117 }
118 return geo_name_list;
119}
120
122{
123 return ((*a)[2] < (*b)[2]) ? a : b;
124}
125
130 std::vector<std::string> const& geo_names,
131 std::string& merged_geo_name)
132{
133 std::vector<GeoLib::Point*> points;
134 std::vector<GeoLib::Polyline*> lines;
135
136 auto DEM_pnts = *geo.getPointVec(geo_names[0]);
137 std::size_t const pnts_per_line = DEM_pnts.size();
138 std::size_t const n_layers = geo_names.size();
139 std::vector<std::size_t> last_line_idx(pnts_per_line, 0);
140
141 auto layer_pnts = *geo.getPointVec(geo_names.back());
142 for (std::size_t i = 0; i < pnts_per_line; ++i)
143 {
144 points.push_back(new GeoLib::Point(
145 *getMinElevationPoint(layer_pnts[i], DEM_pnts[i]), i));
146 last_line_idx[i] = i;
147 }
148 for (int j = n_layers - 2; j >= 0; --j)
149 {
150 GeoLib::Polyline* line = new GeoLib::Polyline(points);
151 for (std::size_t i = 0; i < pnts_per_line; ++i)
152 {
153 line->addPoint(last_line_idx[pnts_per_line - i - 1]);
154 }
155 layer_pnts = *geo.getPointVec(geo_names[j]);
156 for (std::size_t i = 0; i < pnts_per_line; ++i)
157 {
158 // check if for current point the upper layer boundary is actually
159 // located above the upper boundary
160 std::size_t idx = last_line_idx[i];
161 if ((*points[idx])[2] < (*layer_pnts[i])[2])
162 {
163 idx = points.size();
164 // check current point against DEM
165 points.push_back(new GeoLib::Point(
166 *getMinElevationPoint(layer_pnts[i], DEM_pnts[i]), idx));
167 last_line_idx[i] = idx;
168 }
169 line->addPoint(idx);
170 }
171 // close polygon
172 line->addPoint(line->getPointID(0));
173 lines.push_back(line);
174 }
175
176 geo.addPointVec(std::move(points), merged_geo_name,
178 geo.addPolylineVec(std::move(lines), merged_geo_name,
180}
181
183std::pair<Eigen::Matrix3d, double> rotateGeometryToXY(
184 std::vector<GeoLib::Point*>& points)
185{
186 // compute the plane normal
187 auto const [plane_normal, d] =
188 GeoLib::getNewellPlane(points.begin(), points.end());
189 // rotate points into x-y-plane
190 Eigen::Matrix3d const rotation_matrix =
192 GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
193
194 GeoLib::AABB aabb(points.begin(), points.end());
195 double const z_shift =
196 (aabb.getMinPoint()[2] + aabb.getMaxPoint()[2]) / 2.0;
197 std::for_each(points.begin(), points.end(),
198 [z_shift](GeoLib::Point* p) { (*p)[2] -= z_shift; });
199 return {rotation_matrix, z_shift};
200}
201
207 std::string const& output_name,
208 std::string& merged_geo_name,
209 bool const keep_gml_file)
210{
211 std::string const filename(output_name + ".gml");
213 xml.export_name = merged_geo_name;
215
216 geo.removePolylineVec(merged_geo_name);
217 geo.removePointVec(merged_geo_name);
218
219 xml.readFile(filename);
220
221 if (!keep_gml_file)
222 {
223 BaseLib::removeFile(filename);
224 BaseLib::removeFile(filename + ".md5");
225 }
226}
227
230 std::string const& geo_name,
231 std::string const& output_name, double res)
232{
233 std::string const gmsh_geo_name(output_name + ".geo");
234 std::vector<std::string> gmsh_geo;
235 gmsh_geo.push_back(geo_name);
238 0, gmsh_geo, false, false);
239 gmsh_io.writePhysicalGroups(true);
240 if (!BaseLib::IO::writeStringToFile(gmsh_io.writeToString(), gmsh_geo_name))
241 {
242 ERR("Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
243 }
244
245 std::string const gmsh_mesh_name = output_name + ".msh";
246 std::string gmsh_command = "gmsh -2 -algo meshadapt " + gmsh_geo_name;
247 gmsh_command += " -o " + gmsh_mesh_name + " -format msh22";
248 int const return_value = std::system(gmsh_command.c_str());
249 if (return_value != 0)
250 {
251 ERR("Execution of gmsh command returned non-zero status, %d",
252 return_value);
253 }
254 return FileIO::GMSH::readGMSHMesh(gmsh_mesh_name,
255 false /*is_created_with_gmsh2*/);
256}
257
259void rotateMesh(MeshLib::Mesh& mesh, Eigen::Matrix3d const& rot_mat,
260 double const z_shift)
261{
262 std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
263 std::for_each(nodes.begin(), nodes.end(),
264 [z_shift](MeshLib::Node* n) { (*n)[2] += z_shift; });
265 GeoLib::rotatePoints(rot_mat.transpose(), nodes.begin(), nodes.end());
266}
267
270{
271 std::vector<std::size_t> line_idx;
272 std::vector<MeshLib::Element*> const& elems = mesh.getElements();
273 std::for_each(elems.begin(), elems.end(),
274 [&](auto e)
275 {
276 if (e->getGeomType() == MeshLib::MeshElemType::LINE)
277 {
278 line_idx.push_back(e->getID());
279 }
280 });
281 if (line_idx.size() == mesh.getNumberOfElements())
282 {
283 return nullptr;
284 }
285 return MeshToolsLib::removeElements(mesh, line_idx, "mesh");
286}
287
289 std::vector<std::size_t> const& idx_array,
290 std::string const& file_name)
291{
292 std::unique_ptr<MeshLib::Mesh> const boundary(
293 MeshToolsLib::removeElements(mesh, idx_array, "mesh"));
294 if (boundary == nullptr)
295 {
296 ERR("Error extracting boundary '{:s}'", file_name);
297 return;
298 }
299 MeshLib::IO::VtuInterface vtu(boundary.get());
300 vtu.writeToFile(file_name + ".vtu");
301}
302
304 std::string const& output_name,
305 MathLib::Point3d const& pnt_start)
306{
307 auto const edge_length = minMaxEdgeLength(mesh.getElements());
308 double const eps = edge_length.first / 100.0;
309 std::unique_ptr<MeshLib::Mesh> boundary_mesh(
311 mesh,
315
316 auto const& elems = boundary_mesh->getElements();
317 std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
318 bottom_bound_idx;
319 Eigen::Vector2d const anchor(pnt_start[0], pnt_start[1]);
320 for (auto e : elems)
321 {
322 Eigen::Vector2d const n1((*e->getNode(0))[0], (*e->getNode(0))[1]);
323 Eigen::Vector2d const n2((*e->getNode(1))[0], (*e->getNode(1))[1]);
324 Eigen::Vector2d const dist1(n1 - anchor);
325 Eigen::Vector2d const dist2(n2 - anchor);
326 std::size_t const id = e->getID();
327 // elements located at left or right side
328 if ((dist1 - dist2).squaredNorm() < eps)
329 {
330 top_bound_idx.push_back(id);
331 bottom_bound_idx.push_back(id);
332
333 if (dist1.squaredNorm() < eps)
334 {
335 right_bound_idx.push_back(id);
336 }
337 else
338 {
339 left_bound_idx.push_back(id);
340 }
341 continue;
342 }
343 // elements located at top or bottom
344 if (dist2.squaredNorm() < dist1.squaredNorm())
345 {
346 top_bound_idx.push_back(id);
347 left_bound_idx.push_back(id);
348 right_bound_idx.push_back(id);
349 }
350 else
351 {
352 bottom_bound_idx.push_back(id);
353 left_bound_idx.push_back(id);
354 right_bound_idx.push_back(id);
355 }
356 }
357
358 writeBoundary(*boundary_mesh, left_bound_idx, output_name + "_left");
359 writeBoundary(*boundary_mesh, right_bound_idx, output_name + "_right");
360 writeBoundary(*boundary_mesh, top_bound_idx, output_name + "_top");
361 writeBoundary(*boundary_mesh, bottom_bound_idx, output_name + "_bottom");
362}
363
364int main(int argc, char* argv[])
365{
366 QCoreApplication a(argc, argv);
367 TCLAP::CmdLine cmd(
368 "Creates a triangle-mesh of a vertical slice out of a list of input "
369 "layer meshes. The slice is defined by a start- and end-point. In "
370 "addition, the resolution for meshing the extracted slice (i.e. the "
371 "maximum edge length of the domain discretisation) needs to be "
372 "specified. The utility requires access to the meshing utility GMSH to "
373 "work correctly.\n\n"
374 "OpenGeoSys-6 software, version " +
376 ".\n"
377 "Copyright (c) 2012-2025, OpenGeoSys Community "
378 "(http://www.opengeosys.org)",
380 TCLAP::SwitchArg test_arg("t", "testdata", "keep test data", false);
381 cmd.add(test_arg);
382 TCLAP::SwitchArg bound_arg("b", "bounds", "save mesh boundaries", false);
383 cmd.add(bound_arg);
384 TCLAP::ValueArg<double> res_arg(
385 "r", "resolution",
386 "desired edge length of triangles in the resulting slice, (min = 0)",
387 true, 0, "RESOLUTION");
388 cmd.add(res_arg);
389 TCLAP::ValueArg<double> end_y_arg(
390 "", "end-y", "y-coordinates of the end point defining the slice", true,
391 0, "END_Y");
392 cmd.add(end_y_arg);
393 TCLAP::ValueArg<double> end_x_arg(
394 "", "end-x", "x-coordinates of the end point defining the slice", true,
395 0, "END_X");
396 cmd.add(end_x_arg);
397 TCLAP::ValueArg<double> start_y_arg(
398 "", "start-y", "y-coordinates of the start point defining the slice",
399 true, 0, "START_Y");
400 cmd.add(start_y_arg);
401 TCLAP::ValueArg<double> start_x_arg(
402 "", "start-x", "x-coordinates of the start point defining the slice",
403 true, 0, "START_X");
404 cmd.add(start_x_arg);
405 TCLAP::ValueArg<std::string> output_arg(
406 "o", "output", "Output (.vtu). Name of output mesh file", true, "",
407 "OUTPUT_FILE");
408 cmd.add(output_arg);
409 TCLAP::ValueArg<std::string> input_arg(
410 "i", "input",
411 "Input (.vtu | .msh). Name of the input file list containing the paths "
412 "the all input layers "
413 "in correct order from top to bottom",
414 true, "", "INPUT_FILE_LIST");
415 cmd.add(input_arg);
416 cmd.parse(argc, argv);
417
418 BaseLib::MPI::Setup mpi_setup(argc, argv);
419
420 std::string const input_name = input_arg.getValue();
421 std::string const output_name = output_arg.getValue();
422
423 MathLib::Point3d const pnt_start{
424 {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
425 MathLib::Point3d const pnt_end{
426 {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
427 double const length = std::sqrt(MathLib::sqrDist(pnt_start, pnt_end));
428
429 std::size_t const res = std::ceil(length / res_arg.getValue());
430 double const interval_length = length / res;
431
432 std::vector<std::string> const layer_names =
434 if (layer_names.size() < 2)
435 {
436 ERR("At least two layers are required to extract a slice.");
437 return EXIT_FAILURE;
438 }
439
441 std::vector<std::string> const geo_name_list =
442 createGeometries(geo, layer_names, pnt_start, pnt_end, res);
443
444 if (geo_name_list.size() < 2)
445 {
446 ERR("Less than two geometries could be created from layers. Aborting "
447 "extraction...");
448 return EXIT_FAILURE;
449 }
450
451 std::string merged_geo_name = "merged_geometries";
452 mergeGeometries(geo, geo_name_list, merged_geo_name);
453 std::vector<GeoLib::Point*> points = *geo.getPointVec(merged_geo_name);
454 auto const [rot_mat, z_shift] = rotateGeometryToXY(points);
455 consolidateGeometry(geo, output_name, merged_geo_name, test_arg.getValue());
456
457 std::unique_ptr<MeshLib::Mesh> mesh(
458 generateMesh(geo, merged_geo_name, output_name, interval_length));
459 if (!test_arg.getValue())
460 {
461 BaseLib::removeFile(output_name + ".geo");
462 BaseLib::removeFile(output_name + ".msh");
463 }
464 if (mesh == nullptr)
465 {
466 ERR("Error generating mesh... (GMSH was unable to output mesh)");
467 return EXIT_FAILURE;
468 }
469 rotateMesh(*mesh, rot_mat, z_shift);
470 std::unique_ptr<MeshLib::Mesh> new_mesh(removeLineElements(*mesh));
471 if (new_mesh == nullptr)
472 {
473 ERR("Error generating mesh... (GMSH created line mesh)");
474 return EXIT_FAILURE;
475 }
476
477 // collapse all nodes that might have been created due to gmsh physically
478 // separating layers
479 MeshToolsLib::MeshRevision rev(*new_mesh);
480 auto const edge_length = minMaxEdgeLength(new_mesh->getElements());
481 double const minimum_node_distance = edge_length.first / 100.0;
482
483 std::unique_ptr<MeshLib::Mesh> revised_mesh(
484 rev.simplifyMesh("RevisedMesh", minimum_node_distance));
485
486 if (bound_arg.getValue())
487 {
488 extractBoundaries(*revised_mesh, output_name, pnt_start);
489 }
490 MeshLib::IO::VtuInterface vtu(revised_mesh.get());
491 vtu.writeToFile(output_name + ".vtu");
492 return EXIT_SUCCESS;
493}
Definition of the AABB class.
Definition of analytical geometry functions.
Definition of the Element class.
Filename manipulation routines.
Definition of the GEOObjects class.
Definition of the Point class.
Definition of the GeoMapper class.
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
Definition of the MeshRevision class.
Definition of the MeshSurfaceExtraction class.
Definition of the Mesh class.
Definition of the Node class.
Definition of the Point3d class.
Definition of the Polygon class.
Definition of the PolyLine class.
int main(int argc, char *argv[])
void rotateMesh(MeshLib::Mesh &mesh, Eigen::Matrix3d const &rot_mat, double const z_shift)
inverse rotation of the mesh, back into original position
void extractBoundaries(MeshLib::Mesh const &mesh, std::string const &output_name, MathLib::Point3d const &pnt_start)
std::vector< std::string > createGeometries(GeoLib::GEOObjects &geo, std::vector< std::string > const &layer_names, MathLib::Point3d const &pnt_start, MathLib::Point3d const &pnt_end, double const resolution)
creates a mapped line for each of the mesh layers
void writeBoundary(MeshLib::Mesh const &mesh, std::vector< std::size_t > const &idx_array, std::string const &file_name)
void mergeGeometries(GeoLib::GEOObjects &geo, std::vector< std::string > const &geo_names, std::string &merged_geo_name)
MeshLib::Mesh * removeLineElements(MeshLib::Mesh const &mesh)
removes line elements from mesh such that only triangles remain
GeoLib::Polyline * createPolyline(std::vector< GeoLib::Point * > const &points)
creates a polyline to be mapped on a mesh layer
std::pair< Eigen::Matrix3d, double > rotateGeometryToXY(std::vector< GeoLib::Point * > &points)
rotates the merged geometry into the XY-plane
std::vector< GeoLib::Point * > createPoints(MathLib::Point3d const &start, MathLib::Point3d const &end, std::size_t const n_intervals)
creates a vector of sampling points based on the specified resolution
static GeoLib::Point * getMinElevationPoint(GeoLib::Point *a, GeoLib::Point *b)
void consolidateGeometry(GeoLib::GEOObjects &geo, std::string const &output_name, std::string &merged_geo_name, bool const keep_gml_file)
MeshLib::Mesh * generateMesh(GeoLib::GEOObjects &geo, std::string const &geo_name, std::string const &output_name, double res)
converts geometry into GMSH format and creates mesh
Implementation of the VtuInterface class.
Definition of the XmlGmlInterface class.
std::string writeToString()
Writes the object to a string.
Definition Writer.cpp:31
Reads and writes GMSH-files to and from OGS data structures.
void writePhysicalGroups(bool flag)
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
Container class for geometric objects.
Definition GEOObjects.h:57
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
const std::vector< Point * > * getPointVec(const std::string &name) const
void addPointVec(std::vector< Point * > &&points, std::string &name, PointVec::NameIdMap &&pnt_id_name_map, double const eps=std::sqrt(std::numeric_limits< double >::epsilon()))
bool removePointVec(const std::string &name)
bool removePolylineVec(const std::string &name)
Reads and writes GeoObjects to and from XML files.
int readFile(const QString &fileName) override
Reads an xml-file containing geometric object definitions into the GEOObjects used in the constructor...
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:40
std::size_t getPointID(std::size_t const i) const
Definition Polyline.cpp:160
virtual bool addPoint(std::size_t pnt_id)
Definition Polyline.cpp:35
std::map< std::string, std::size_t > NameIdMap
Definition TemplateVec.h:41
A set of tools for mapping the elevation of geometric objects.
Definition GeoMapper.h:40
void mapOnMesh(MeshLib::Mesh const *const mesh)
Definition GeoMapper.cpp:73
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
bool writeToFile(std::filesystem::path const &file_path)
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
MeshLib::Mesh * simplifyMesh(const std::string &new_mesh_name, double eps, unsigned min_elem_dim=1) const
std::vector< std::string > readStringListFromFile(std::string const &filename)
Reads non-empty lines from a list of strings from a file into a vector.
int writeStringToFile(std::string_view content, std::filesystem::path const &file_path)
Definition Writer.cpp:45
void removeFile(std::string const &filename)
@ FixedMeshDensity
set the parameter with a fixed value
MeshLib::Mesh * readGMSHMesh(std::string const &fname, bool const is_created_with_gmsh2)
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const &n)
std::pair< Eigen::Vector3d, double > getNewellPlane(InputIterator pnts_begin, InputIterator pnts_end)
GITINFOLIB_EXPORT const std::string ogs_version
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition Properties.h:191
std::unique_ptr< MeshLib::Mesh > getBoundaryElementsAsMesh(MeshLib::Mesh const &bulk_mesh, std::string_view subsfc_node_id_prop_name, std::string_view subsfc_element_id_prop_name, std::string_view face_id_prop_name)
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)
Definition of readMeshFromFile function.