OGS
VerticalSliceFromLayers.cpp File Reference

Detailed Description

Definition in file VerticalSliceFromLayers.cpp.

#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/MPI.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 206 of file VerticalSliceFromLayers.cpp.

210{
211 std::string const filename(output_name + ".gml");
213 xml.export_name = merged_geo_name;
214 BaseLib::IO::writeStringToFile(xml.writeToString(), filename);
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}
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:45
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 80 of file VerticalSliceFromLayers.cpp.

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}
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
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:41
A set of tools for mapping the elevation of geometric objects.
Definition GeoMapper.h:40
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 46 of file VerticalSliceFromLayers.cpp.

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}
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26

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 68 of file VerticalSliceFromLayers.cpp.

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}
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:40
virtual bool addPoint(std::size_t pnt_id)
Definition Polyline.cpp:35

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 303 of file VerticalSliceFromLayers.cpp.

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}
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)
Definition Properties.h:191
std::pair< double, double > minMaxEdgeLength(std::vector< Element * > const &elements)
Returns the minimum and maximum edge length for given elements.
Definition Mesh.cpp:190
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 229 of file VerticalSliceFromLayers.cpp.

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}
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()

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

Definition at line 121 of file VerticalSliceFromLayers.cpp.

122{
123 return ((*a)[2] < (*b)[2]) ? a : b;
124}

Referenced by mergeGeometries().

◆ main()

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

Definition at line 364 of file VerticalSliceFromLayers.cpp.

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}
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:57
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.
GITINFOLIB_EXPORT const std::string ogs_version

References consolidateGeometry(), createGeometries(), ERR(), extractBoundaries(), generateMesh(), GeoLib::GEOObjects::getPointVec(), 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 129 of file VerticalSliceFromLayers.cpp.

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}
static GeoLib::Point * getMinElevationPoint(GeoLib::Point *a, GeoLib::Point *b)
std::size_t getPointID(std::size_t const i) const
Definition Polyline.cpp:160

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 269 of file VerticalSliceFromLayers.cpp.

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}
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 183 of file VerticalSliceFromLayers.cpp.

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}
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
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 259 of file VerticalSliceFromLayers.cpp.

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}
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108

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 288 of file VerticalSliceFromLayers.cpp.

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}

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

Referenced by extractBoundaries().