OGS
VerticalSliceFromLayers.cpp File Reference
#include <tclap/CmdLine.h>
#include <QCoreApplication>
#include <algorithm>
#include <cmath>
#include <memory>
#include <string>
#include <vector>
#include "Applications/FileIO/Gmsh/GMSHInterface.h"
#include "Applications/FileIO/Gmsh/GmshReader.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/IO/readStringListFromFile.h"
#include "BaseLib/Logging.h"
#include "BaseLib/MPI.h"
#include "BaseLib/TCLAPArguments.h"
#include "GeoLib/AABB.h"
#include "GeoLib/AnalyticalGeometry.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
#include "GeoLib/Point.h"
#include "GeoLib/Polygon.h"
#include "GeoLib/Polyline.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/MathTools.h"
#include "MathLib/Point3d.h"
#include "MeshGeoToolsLib/GeoMapper.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/ElementSearch.h"
#include "MeshLib/Node.h"
#include "MeshToolsLib/MeshEditing/MeshRevision.h"
#include "MeshToolsLib/MeshEditing/RemoveMeshComponents.h"
#include "MeshToolsLib/MeshSurfaceExtraction.h"
Include dependency graph for VerticalSliceFromLayers.cpp:

Go to the source code of this file.

Functions

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
GeoLib::PolylinecreatePolyline (std::vector< GeoLib::Point * > const &points)
 creates a polyline to be mapped on a mesh layer
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
static GeoLib::PointgetMinElevationPoint (GeoLib::Point *a, GeoLib::Point *b)
void mergeGeometries (GeoLib::GEOObjects &geo, std::vector< std::string > const &geo_names, std::string &merged_geo_name)
std::pair< Eigen::Matrix3d, double > rotateGeometryToXY (std::vector< GeoLib::Point * > &points)
 rotates the merged geometry into the XY-plane
void consolidateGeometry (GeoLib::GEOObjects &geo, std::string const &output_name, std::string &merged_geo_name, bool const keep_gml_file)
MeshLib::MeshgenerateMesh (GeoLib::GEOObjects &geo, std::string const &geo_name, std::string const &output_name, double res)
 converts geometry into GMSH format and creates mesh
void rotateMesh (MeshLib::Mesh &mesh, Eigen::Matrix3d const &rot_mat, double const z_shift)
 inverse rotation of the mesh, back into original position
MeshLib::MeshremoveLineElements (MeshLib::Mesh const &mesh)
 removes line elements from mesh such that only triangles remain
void writeBoundary (MeshLib::Mesh const &mesh, std::vector< std::size_t > const &idx_array, std::string const &file_name)
void extractBoundaries (MeshLib::Mesh const &mesh, std::string const &output_name, MathLib::Point3d const &pnt_start)
int main (int argc, char *argv[])

Function Documentation

◆ consolidateGeometry()

void consolidateGeometry ( GeoLib::GEOObjects & geo,
std::string const & output_name,
std::string & merged_geo_name,
bool const keep_gml_file )

This encapsulates a workaround: For unknown reasons, the GML->GEO converter will not work correctly when inputting the merged geometry directly. However, if the geometry is saved to a file and immediately loaded again, everything works fine.

Definition at line 202 of file VerticalSliceFromLayers.cpp.

206{
207 std::string const filename(output_name + ".gml");
209 xml.export_name = merged_geo_name;
210 BaseLib::IO::writeStringToFile(xml.writeToString(), filename);
211
212 geo.removePolylineVec(merged_geo_name);
213 geo.removePointVec(merged_geo_name);
214
215 xml.readFile(filename);
216
217 if (!keep_gml_file)
218 {
219 BaseLib::removeFile(filename);
220 BaseLib::removeFile(filename + ".md5");
221 }
222}
bool removePointVec(const std::string &name)
bool removePolylineVec(const std::string &name)
Reads and writes GeoObjects to and from XML files.
int writeStringToFile(std::string_view content, std::filesystem::path const &file_path)
Definition Writer.cpp:34
void removeFile(std::string const &filename)

References BaseLib::IO::XMLInterface::export_name, GeoLib::IO::XmlGmlInterface::readFile(), BaseLib::removeFile(), GeoLib::GEOObjects::removePointVec(), GeoLib::GEOObjects::removePolylineVec(), BaseLib::IO::writeStringToFile(), and BaseLib::IO::Writer::writeToString().

Referenced by main().

◆ createGeometries()

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

Definition at line 76 of file VerticalSliceFromLayers.cpp.

82{
83 std::vector<std::string> geo_name_list;
84 std::size_t const n_layers = layer_names.size();
85 for (std::size_t i = 0; i < n_layers; ++i)
86 {
87 std::unique_ptr<MeshLib::Mesh> const layer(
88 MeshLib::IO::readMeshFromFile(layer_names[i]));
89 if (layer == nullptr)
90 {
91 ERR("Could not read file {:s}. Skipping layer...", layer_names[i]);
92 continue;
93 }
94 if (layer->getDimension() != 2)
95 {
96 ERR("Layer {:d} is not a 2D mesh. Skipping layer...", i);
97 continue;
98 }
99
100 std::string geo_name(std::to_string(i));
101 auto points(createPoints(pnt_start, pnt_end, resolution));
102 geo.addPointVec(std::move(points), geo_name,
104
105 std::vector<GeoLib::Polyline*> lines;
106 lines.push_back(createPolyline(*geo.getPointVec(geo_name)));
107 geo.addPolylineVec(std::move(lines), geo_name,
109
110 MeshGeoToolsLib::GeoMapper mapper(geo, geo_name);
111 mapper.mapOnMesh(layer.get());
112 geo_name_list.push_back(geo_name);
113 }
114 return geo_name_list;
115}
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
GeoLib::Polyline * createPolyline(std::vector< GeoLib::Point * > const &points)
creates a polyline to be mapped on a mesh layer
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
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()))
std::map< std::string, std::size_t > NameIdMap
Definition TemplateVec.h:30
A set of tools for mapping the elevation of geometric objects.
Definition GeoMapper.h:29
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)

References GeoLib::GEOObjects::addPointVec(), GeoLib::GEOObjects::addPolylineVec(), createPoints(), createPolyline(), ERR(), GeoLib::GEOObjects::getPointVec(), MeshGeoToolsLib::GeoMapper::mapOnMesh(), and MeshLib::IO::readMeshFromFile().

Referenced by main().

◆ createPoints()

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

Definition at line 42 of file VerticalSliceFromLayers.cpp.

45{
46 std::vector<GeoLib::Point*> points;
47 points.push_back(new GeoLib::Point(start, 0));
48 double const length = std::sqrt(MathLib::sqrDist(start, end));
49 double const interval = length / n_intervals;
50
51 GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
52 (end[2] - start[2]));
53 for (std::size_t i = 1; i < n_intervals; ++i)
54 {
55 points.push_back(new GeoLib::Point(
56 start[0] + ((i * interval) / length * vec[0]),
57 start[1] + ((i * interval) / length * vec[1]), 0, i));
58 }
59 points.push_back(new GeoLib::Point(end, n_intervals));
60 return points;
61}
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:19

References MathLib::sqrDist().

Referenced by createGeometries().

◆ createPolyline()

GeoLib::Polyline * createPolyline ( std::vector< GeoLib::Point * > const & points)

creates a polyline to be mapped on a mesh layer

Definition at line 64 of file VerticalSliceFromLayers.cpp.

65{
66 GeoLib::Polyline* line = new GeoLib::Polyline(points);
67 std::size_t const length = points.size();
68 for (std::size_t i = 0; i < length; ++i)
69 {
70 line->addPoint(i);
71 }
72 return line;
73}
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:29
virtual bool addPoint(std::size_t pnt_id)
Definition Polyline.cpp:24

References GeoLib::Polyline::addPoint().

Referenced by createGeometries().

◆ extractBoundaries()

void extractBoundaries ( MeshLib::Mesh const & mesh,
std::string const & output_name,
MathLib::Point3d const & pnt_start )

Definition at line 299 of file VerticalSliceFromLayers.cpp.

302{
303 auto const edge_length = minMaxEdgeLength(mesh.getElements());
304 double const eps = edge_length.first / 100.0;
305 std::unique_ptr<MeshLib::Mesh> boundary_mesh(
307 mesh,
311
312 auto const& elems = boundary_mesh->getElements();
313 std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
314 bottom_bound_idx;
315 Eigen::Vector2d const anchor(pnt_start[0], pnt_start[1]);
316 for (auto e : elems)
317 {
318 Eigen::Vector2d const n1((*e->getNode(0))[0], (*e->getNode(0))[1]);
319 Eigen::Vector2d const n2((*e->getNode(1))[0], (*e->getNode(1))[1]);
320 Eigen::Vector2d const dist1(n1 - anchor);
321 Eigen::Vector2d const dist2(n2 - anchor);
322 std::size_t const id = e->getID();
323 // elements located at left or right side
324 if ((dist1 - dist2).squaredNorm() < eps)
325 {
326 top_bound_idx.push_back(id);
327 bottom_bound_idx.push_back(id);
328
329 if (dist1.squaredNorm() < eps)
330 {
331 right_bound_idx.push_back(id);
332 }
333 else
334 {
335 left_bound_idx.push_back(id);
336 }
337 continue;
338 }
339 // elements located at top or bottom
340 if (dist2.squaredNorm() < dist1.squaredNorm())
341 {
342 top_bound_idx.push_back(id);
343 left_bound_idx.push_back(id);
344 right_bound_idx.push_back(id);
345 }
346 else
347 {
348 bottom_bound_idx.push_back(id);
349 left_bound_idx.push_back(id);
350 right_bound_idx.push_back(id);
351 }
352 }
353
354 writeBoundary(*boundary_mesh, left_bound_idx, output_name + "_left");
355 writeBoundary(*boundary_mesh, right_bound_idx, output_name + "_right");
356 writeBoundary(*boundary_mesh, top_bound_idx, output_name + "_top");
357 writeBoundary(*boundary_mesh, bottom_bound_idx, output_name + "_bottom");
358}
void writeBoundary(MeshLib::Mesh const &mesh, std::vector< std::size_t > const &idx_array, std::string const &file_name)
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.
Definition Mesh.cpp:179
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)

References MeshLib::Cell, MeshLib::Face, MeshToolsLib::BoundaryExtraction::getBoundaryElementsAsMesh(), MeshLib::getBulkIDString(), MeshLib::Mesh::getElements(), MeshLib::Node, and writeBoundary().

Referenced by main().

◆ generateMesh()

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

Definition at line 225 of file VerticalSliceFromLayers.cpp.

228{
229 std::string const gmsh_geo_name(output_name + ".geo");
230 std::vector<std::string> gmsh_geo;
231 gmsh_geo.push_back(geo_name);
234 0, gmsh_geo, false, false);
235 gmsh_io.writePhysicalGroups(true);
236 if (!BaseLib::IO::writeStringToFile(gmsh_io.writeToString(), gmsh_geo_name))
237 {
238 ERR("Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
239 }
240
241 std::string const gmsh_mesh_name = output_name + ".msh";
242 std::string gmsh_command = "gmsh -2 -algo meshadapt " + gmsh_geo_name;
243 gmsh_command += " -o " + gmsh_mesh_name + " -format msh22";
244 int const return_value = std::system(gmsh_command.c_str());
245 if (return_value != 0)
246 {
247 ERR("Execution of gmsh command returned non-zero status, %d",
248 return_value);
249 }
250 return FileIO::GMSH::readGMSHMesh(gmsh_mesh_name,
251 false /*is_created_with_gmsh2*/);
252}
Reads and writes GMSH-files to and from OGS data structures.
@ FixedMeshDensity
set the parameter with a fixed value
MeshLib::Mesh * readGMSHMesh(std::string const &fname, bool const is_created_with_gmsh2)

References ERR(), FileIO::GMSH::FixedMeshDensity, FileIO::GMSH::readGMSHMesh(), FileIO::GMSH::GMSHInterface::writePhysicalGroups(), BaseLib::IO::writeStringToFile(), and BaseLib::IO::Writer::writeToString().

Referenced by main().

◆ getMinElevationPoint()

GeoLib::Point * getMinElevationPoint ( GeoLib::Point * a,
GeoLib::Point * b )
static

Definition at line 117 of file VerticalSliceFromLayers.cpp.

118{
119 return ((*a)[2] < (*b)[2]) ? a : b;
120}

Referenced by mergeGeometries().

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 360 of file VerticalSliceFromLayers.cpp.

361{
362 QCoreApplication a(argc, argv);
363 TCLAP::CmdLine cmd(
364 "Creates a triangle-mesh of a vertical slice out of a list of input "
365 "layer meshes. The slice is defined by a start- and end-point. In "
366 "addition, the resolution for meshing the extracted slice (i.e. the "
367 "maximum edge length of the domain discretisation) needs to be "
368 "specified. The utility requires access to the meshing utility GMSH to "
369 "work correctly.\n\n"
370 "OpenGeoSys-6 software, version " +
372 ".\n"
373 "Copyright (c) 2012-2026, OpenGeoSys Community "
374 "(http://www.opengeosys.org)",
376 TCLAP::SwitchArg test_arg("t", "testdata", "keep test data", false);
377 cmd.add(test_arg);
378 TCLAP::SwitchArg bound_arg("b", "bounds", "save mesh boundaries", false);
379 cmd.add(bound_arg);
380 TCLAP::ValueArg<double> res_arg(
381 "r", "resolution",
382 "desired edge length of triangles in the resulting slice, (min = 0)",
383 true, 0, "RESOLUTION");
384 cmd.add(res_arg);
385 TCLAP::ValueArg<double> end_y_arg(
386 "", "end-y", "y-coordinates of the end point defining the slice", true,
387 0, "END_Y");
388 cmd.add(end_y_arg);
389 TCLAP::ValueArg<double> end_x_arg(
390 "", "end-x", "x-coordinates of the end point defining the slice", true,
391 0, "END_X");
392 cmd.add(end_x_arg);
393 TCLAP::ValueArg<double> start_y_arg(
394 "", "start-y", "y-coordinates of the start point defining the slice",
395 true, 0, "START_Y");
396 cmd.add(start_y_arg);
397 TCLAP::ValueArg<double> start_x_arg(
398 "", "start-x", "x-coordinates of the start point defining the slice",
399 true, 0, "START_X");
400 cmd.add(start_x_arg);
401 TCLAP::ValueArg<std::string> output_arg(
402 "o", "output", "Output (.vtu). Name of output mesh file", true, "",
403 "OUTPUT_FILE");
404 cmd.add(output_arg);
405 TCLAP::ValueArg<std::string> input_arg(
406 "i", "input",
407 "Input (.vtu | .msh). Name of the input file list containing the paths "
408 "the all input layers "
409 "in correct order from top to bottom",
410 true, "", "INPUT_FILE_LIST");
411 cmd.add(input_arg);
412 auto log_level_arg = BaseLib::makeLogLevelArg();
413 cmd.add(log_level_arg);
414 cmd.parse(argc, argv);
415
416 BaseLib::MPI::Setup mpi_setup(argc, argv);
417 BaseLib::initOGSLogger(log_level_arg.getValue());
418
419 std::string const input_name = input_arg.getValue();
420 std::string const output_name = output_arg.getValue();
421
422 MathLib::Point3d const pnt_start{
423 {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
424 MathLib::Point3d const pnt_end{
425 {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
426 double const length = std::sqrt(MathLib::sqrDist(pnt_start, pnt_end));
427
428 std::size_t const res = std::ceil(length / res_arg.getValue());
429 double const interval_length = length / res;
430
431 std::vector<std::string> const layer_names =
433 if (layer_names.size() < 2)
434 {
435 ERR("At least two layers are required to extract a slice.");
436 return EXIT_FAILURE;
437 }
438
440 std::vector<std::string> const geo_name_list =
441 createGeometries(geo, layer_names, pnt_start, pnt_end, res);
442
443 if (geo_name_list.size() < 2)
444 {
445 ERR("Less than two geometries could be created from layers. Aborting "
446 "extraction...");
447 return EXIT_FAILURE;
448 }
449
450 std::string merged_geo_name = "merged_geometries";
451 mergeGeometries(geo, geo_name_list, merged_geo_name);
452 std::vector<GeoLib::Point*> points = *geo.getPointVec(merged_geo_name);
453 auto const [rot_mat, z_shift] = rotateGeometryToXY(points);
454 consolidateGeometry(geo, output_name, merged_geo_name, test_arg.getValue());
455
456 std::unique_ptr<MeshLib::Mesh> mesh(
457 generateMesh(geo, merged_geo_name, output_name, interval_length));
458 if (!test_arg.getValue())
459 {
460 BaseLib::removeFile(output_name + ".geo");
461 BaseLib::removeFile(output_name + ".msh");
462 }
463 if (mesh == nullptr)
464 {
465 ERR("Error generating mesh... (GMSH was unable to output mesh)");
466 return EXIT_FAILURE;
467 }
468 rotateMesh(*mesh, rot_mat, z_shift);
469 std::unique_ptr<MeshLib::Mesh> new_mesh(removeLineElements(*mesh));
470 if (new_mesh == nullptr)
471 {
472 ERR("Error generating mesh... (GMSH created line mesh)");
473 return EXIT_FAILURE;
474 }
475
476 // collapse all nodes that might have been created due to gmsh physically
477 // separating layers
478 MeshToolsLib::MeshRevision rev(*new_mesh);
479 auto const edge_length = minMaxEdgeLength(new_mesh->getElements());
480 double const minimum_node_distance = edge_length.first / 100.0;
481
482 std::unique_ptr<MeshLib::Mesh> revised_mesh(
483 rev.simplifyMesh("RevisedMesh", minimum_node_distance));
484
485 if (bound_arg.getValue())
486 {
487 extractBoundaries(*revised_mesh, output_name, pnt_start);
488 }
489 MeshLib::IO::VtuInterface vtu(revised_mesh.get());
490 vtu.writeToFile(output_name + ".vtu");
491 return EXIT_SUCCESS;
492}
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 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
std::pair< Eigen::Matrix3d, double > rotateGeometryToXY(std::vector< GeoLib::Point * > &points)
rotates the merged geometry into the XY-plane
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
Container class for geometric objects.
Definition GEOObjects.h:46
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
std::vector< std::string > readStringListFromFile(std::string const &filename)
Reads non-empty lines from a list of strings from a file into a vector.
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:56
GITINFOLIB_EXPORT const std::string ogs_version

References consolidateGeometry(), createGeometries(), ERR(), extractBoundaries(), generateMesh(), GeoLib::GEOObjects::getPointVec(), BaseLib::initOGSLogger(), BaseLib::makeLogLevelArg(), mergeGeometries(), GitInfoLib::GitInfo::ogs_version, BaseLib::IO::readStringListFromFile(), BaseLib::removeFile(), removeLineElements(), rotateGeometryToXY(), rotateMesh(), MeshToolsLib::MeshRevision::simplifyMesh(), MathLib::sqrDist(), and MeshLib::IO::VtuInterface::writeToFile().

◆ mergeGeometries()

void mergeGeometries ( GeoLib::GEOObjects & geo,
std::vector< std::string > const & geo_names,
std::string & merged_geo_name )

Merges all layer geometries into one. Each layer is specified by one polygon. This will ensure that GMSH deals with the subsequent assignment of material groups automatically.

Definition at line 125 of file VerticalSliceFromLayers.cpp.

128{
129 std::vector<GeoLib::Point*> points;
130 std::vector<GeoLib::Polyline*> lines;
131
132 auto DEM_pnts = *geo.getPointVec(geo_names[0]);
133 std::size_t const pnts_per_line = DEM_pnts.size();
134 std::size_t const n_layers = geo_names.size();
135 std::vector<std::size_t> last_line_idx(pnts_per_line, 0);
136
137 auto layer_pnts = *geo.getPointVec(geo_names.back());
138 for (std::size_t i = 0; i < pnts_per_line; ++i)
139 {
140 points.push_back(new GeoLib::Point(
141 *getMinElevationPoint(layer_pnts[i], DEM_pnts[i]), i));
142 last_line_idx[i] = i;
143 }
144 for (int j = n_layers - 2; j >= 0; --j)
145 {
146 GeoLib::Polyline* line = new GeoLib::Polyline(points);
147 for (std::size_t i = 0; i < pnts_per_line; ++i)
148 {
149 line->addPoint(last_line_idx[pnts_per_line - i - 1]);
150 }
151 layer_pnts = *geo.getPointVec(geo_names[j]);
152 for (std::size_t i = 0; i < pnts_per_line; ++i)
153 {
154 // check if for current point the upper layer boundary is actually
155 // located above the upper boundary
156 std::size_t idx = last_line_idx[i];
157 if ((*points[idx])[2] < (*layer_pnts[i])[2])
158 {
159 idx = points.size();
160 // check current point against DEM
161 points.push_back(new GeoLib::Point(
162 *getMinElevationPoint(layer_pnts[i], DEM_pnts[i]), idx));
163 last_line_idx[i] = idx;
164 }
165 line->addPoint(idx);
166 }
167 // close polygon
168 line->addPoint(line->getPointID(0));
169 lines.push_back(line);
170 }
171
172 geo.addPointVec(std::move(points), merged_geo_name,
174 geo.addPolylineVec(std::move(lines), merged_geo_name,
176}
static GeoLib::Point * getMinElevationPoint(GeoLib::Point *a, GeoLib::Point *b)
std::size_t getPointID(std::size_t const i) const
Definition Polyline.cpp:149

References GeoLib::Polyline::addPoint(), GeoLib::GEOObjects::addPointVec(), GeoLib::GEOObjects::addPolylineVec(), getMinElevationPoint(), GeoLib::Polyline::getPointID(), and GeoLib::GEOObjects::getPointVec().

Referenced by main().

◆ removeLineElements()

MeshLib::Mesh * removeLineElements ( MeshLib::Mesh const & mesh)

removes line elements from mesh such that only triangles remain

Definition at line 265 of file VerticalSliceFromLayers.cpp.

266{
267 std::vector<std::size_t> line_idx;
268 std::vector<MeshLib::Element*> const& elems = mesh.getElements();
269 std::for_each(elems.begin(), elems.end(),
270 [&](auto e)
271 {
272 if (e->getGeomType() == MeshLib::MeshElemType::LINE)
273 {
274 line_idx.push_back(e->getID());
275 }
276 });
277 if (line_idx.size() == mesh.getNumberOfElements())
278 {
279 return nullptr;
280 }
281 return MeshToolsLib::removeElements(mesh, line_idx, "mesh");
282}
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)

References MeshLib::Mesh::getElements().

Referenced by main().

◆ rotateGeometryToXY()

std::pair< Eigen::Matrix3d, double > rotateGeometryToXY ( std::vector< GeoLib::Point * > & points)

rotates the merged geometry into the XY-plane

Definition at line 179 of file VerticalSliceFromLayers.cpp.

181{
182 // compute the plane normal
183 auto const [plane_normal, d] =
184 GeoLib::getNewellPlane(points.begin(), points.end());
185 // rotate points into x-y-plane
186 Eigen::Matrix3d const rotation_matrix =
188 GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
189
190 GeoLib::AABB aabb(points.begin(), points.end());
191 double const z_shift =
192 (aabb.getMinPoint()[2] + aabb.getMaxPoint()[2]) / 2.0;
193 std::for_each(points.begin(), points.end(),
194 [z_shift](GeoLib::Point* p) { (*p)[2] -= z_shift; });
195 return {rotation_matrix, z_shift};
196}
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:45
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)

References GeoLib::computeRotationMatrixToXY(), GeoLib::AABB::getMaxPoint(), GeoLib::AABB::getMinPoint(), GeoLib::getNewellPlane(), and GeoLib::rotatePoints().

Referenced by main().

◆ rotateMesh()

void rotateMesh ( MeshLib::Mesh & mesh,
Eigen::Matrix3d const & rot_mat,
double const z_shift )

inverse rotation of the mesh, back into original position

Definition at line 255 of file VerticalSliceFromLayers.cpp.

257{
258 std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
259 std::for_each(nodes.begin(), nodes.end(),
260 [z_shift](MeshLib::Node* n) { (*n)[2] += z_shift; });
261 GeoLib::rotatePoints(rot_mat.transpose(), nodes.begin(), nodes.end());
262}
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:97

References MeshLib::Mesh::getNodes(), and GeoLib::rotatePoints().

Referenced by main().

◆ writeBoundary()

void writeBoundary ( MeshLib::Mesh const & mesh,
std::vector< std::size_t > const & idx_array,
std::string const & file_name )

Definition at line 284 of file VerticalSliceFromLayers.cpp.

287{
288 std::unique_ptr<MeshLib::Mesh> const boundary(
289 MeshToolsLib::removeElements(mesh, idx_array, "mesh"));
290 if (boundary == nullptr)
291 {
292 ERR("Error extracting boundary '{:s}'", file_name);
293 return;
294 }
295 MeshLib::IO::VtuInterface vtu(boundary.get());
296 vtu.writeToFile(file_name + ".vtu");
297}

References ERR(), MeshToolsLib::removeElements(), and MeshLib::IO::VtuInterface::writeToFile().

Referenced by extractBoundaries().