10#include <tclap/CmdLine.h>
12#include <QCoreApplication>
48 std::size_t
const n_intervals)
50 std::vector<GeoLib::Point*> points;
53 double const interval = length / n_intervals;
55 GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
57 for (std::size_t i = 1; i < n_intervals; ++i)
60 start[0] + ((i * interval) / length * vec[0]),
61 start[1] + ((i * interval) / length * vec[1]), 0, i));
71 std::size_t
const length = points.size();
72 for (std::size_t i = 0; i < length; ++i)
82 std::vector<std::string>
const& layer_names,
85 double const resolution)
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)
91 std::unique_ptr<MeshLib::Mesh>
const layer(
95 ERR(
"Could not read file {:s}. Skipping layer...", layer_names[i]);
98 if (layer->getDimension() != 2)
100 ERR(
"Layer {:d} is not a 2D mesh. Skipping layer...", i);
104 std::string geo_name(std::to_string(i));
105 auto points(
createPoints(pnt_start, pnt_end, resolution));
109 std::vector<GeoLib::Polyline*> lines;
116 geo_name_list.push_back(geo_name);
118 return geo_name_list;
123 return ((*a)[2] < (*b)[2]) ? a : b;
130 std::vector<std::string>
const& geo_names,
131 std::string& merged_geo_name)
133 std::vector<GeoLib::Point*> points;
134 std::vector<GeoLib::Polyline*> lines;
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);
141 auto layer_pnts = *geo.
getPointVec(geo_names.back());
142 for (std::size_t i = 0; i < pnts_per_line; ++i)
146 last_line_idx[i] = i;
148 for (
int j = n_layers - 2; j >= 0; --j)
151 for (std::size_t i = 0; i < pnts_per_line; ++i)
153 line->
addPoint(last_line_idx[pnts_per_line - i - 1]);
156 for (std::size_t i = 0; i < pnts_per_line; ++i)
160 std::size_t idx = last_line_idx[i];
161 if ((*points[idx])[2] < (*layer_pnts[i])[2])
167 last_line_idx[i] = idx;
173 lines.push_back(line);
176 geo.
addPointVec(std::move(points), merged_geo_name,
184 std::vector<GeoLib::Point*>& points)
187 auto const [plane_normal, d] =
190 Eigen::Matrix3d
const rotation_matrix =
195 double const z_shift =
197 std::for_each(points.begin(), points.end(),
199 return {rotation_matrix, z_shift};
207 std::string
const& output_name,
208 std::string& merged_geo_name,
209 bool const keep_gml_file)
211 std::string
const filename(output_name +
".gml");
230 std::string
const& geo_name,
231 std::string
const& output_name,
double res)
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);
242 ERR(
"Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
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)
251 ERR(
"Execution of gmsh command returned non-zero status, %d",
260 double const z_shift)
262 std::vector<MeshLib::Node*>
const& nodes = mesh.
getNodes();
263 std::for_each(nodes.begin(), nodes.end(),
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(),
276 if (e->getGeomType() == MeshLib::MeshElemType::LINE)
278 line_idx.push_back(e->getID());
281 if (line_idx.size() == mesh.getNumberOfElements())
289 std::vector<std::size_t>
const& idx_array,
290 std::string
const& file_name)
292 std::unique_ptr<MeshLib::Mesh>
const boundary(
294 if (boundary ==
nullptr)
296 ERR(
"Error extracting boundary '{:s}'", file_name);
304 std::string
const& output_name,
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(
316 auto const& elems = boundary_mesh->getElements();
317 std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
319 Eigen::Vector2d
const anchor(pnt_start[0], pnt_start[1]);
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();
328 if ((dist1 - dist2).squaredNorm() < eps)
330 top_bound_idx.push_back(
id);
331 bottom_bound_idx.push_back(
id);
333 if (dist1.squaredNorm() < eps)
335 right_bound_idx.push_back(
id);
339 left_bound_idx.push_back(
id);
344 if (dist2.squaredNorm() < dist1.squaredNorm())
346 top_bound_idx.push_back(
id);
347 left_bound_idx.push_back(
id);
348 right_bound_idx.push_back(
id);
352 bottom_bound_idx.push_back(
id);
353 left_bound_idx.push_back(
id);
354 right_bound_idx.push_back(
id);
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");
364int main(
int argc,
char* argv[])
366 QCoreApplication a(argc, argv);
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 " +
377 "Copyright (c) 2012-2025, OpenGeoSys Community "
378 "(http://www.opengeosys.org)",
380 TCLAP::SwitchArg test_arg(
"t",
"testdata",
"keep test data",
false);
382 TCLAP::SwitchArg bound_arg(
"b",
"bounds",
"save mesh boundaries",
false);
384 TCLAP::ValueArg<double> res_arg(
386 "desired edge length of triangles in the resulting slice, (min = 0)",
387 true, 0,
"RESOLUTION");
389 TCLAP::ValueArg<double> end_y_arg(
390 "",
"end-y",
"y-coordinates of the end point defining the slice",
true,
393 TCLAP::ValueArg<double> end_x_arg(
394 "",
"end-x",
"x-coordinates of the end point defining the slice",
true,
397 TCLAP::ValueArg<double> start_y_arg(
398 "",
"start-y",
"y-coordinates of the start point defining the slice",
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",
404 cmd.add(start_x_arg);
405 TCLAP::ValueArg<std::string> output_arg(
406 "o",
"output",
"Output (.vtu). Name of output mesh file",
true,
"",
409 TCLAP::ValueArg<std::string> input_arg(
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");
416 cmd.parse(argc, argv);
420 std::string
const input_name = input_arg.getValue();
421 std::string
const output_name = output_arg.getValue();
424 {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
426 {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
429 std::size_t
const res = std::ceil(length / res_arg.getValue());
430 double const interval_length = length / res;
432 std::vector<std::string>
const layer_names =
434 if (layer_names.size() < 2)
436 ERR(
"At least two layers are required to extract a slice.");
441 std::vector<std::string>
const geo_name_list =
444 if (geo_name_list.size() < 2)
446 ERR(
"Less than two geometries could be created from layers. Aborting "
451 std::string merged_geo_name =
"merged_geometries";
453 std::vector<GeoLib::Point*> points = *geo.
getPointVec(merged_geo_name);
457 std::unique_ptr<MeshLib::Mesh> mesh(
458 generateMesh(geo, merged_geo_name, output_name, interval_length));
459 if (!test_arg.getValue())
466 ERR(
"Error generating mesh... (GMSH was unable to output mesh)");
471 if (new_mesh ==
nullptr)
473 ERR(
"Error generating mesh... (GMSH created line mesh)");
480 auto const edge_length = minMaxEdgeLength(new_mesh->getElements());
481 double const minimum_node_distance = edge_length.first / 100.0;
483 std::unique_ptr<MeshLib::Mesh> revised_mesh(
484 rev.
simplifyMesh(
"RevisedMesh", minimum_node_distance));
486 if (bound_arg.getValue())
Definition of the AABB class.
Definition of analytical geometry functions.
Definition of the Element class.
Definition of the GEOObjects class.
Definition of the Point class.
Definition of the GeoMapper class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the MeshRevision 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.
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 ...
Eigen::Vector3d const & getMaxPoint() const
Eigen::Vector3d const & getMinPoint() const
Container class for geometric objects.
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 ...
std::size_t getPointID(std::size_t const i) const
virtual bool addPoint(std::size_t pnt_id)
std::map< std::string, std::size_t > NameIdMap
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.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
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)
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)
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
constexpr std::string_view getBulkIDString(MeshItemType mesh_item_type)
Definition of readMeshFromFile function.