OGS
VerticalSliceFromLayers.cpp
Go to the documentation of this file.
1 
10 #include <algorithm>
11 #include <cmath>
12 #include <memory>
13 #include <string>
14 #include <vector>
15 
16 // ThirdParty
17 #include <tclap/CmdLine.h>
18 
19 #include <QCoreApplication>
20 
23 #include "BaseLib/FileTools.h"
25 #include "GeoLib/AABB.h"
27 #include "GeoLib/GEOObjects.h"
29 #include "GeoLib/Point.h"
30 #include "GeoLib/Polygon.h"
31 #include "GeoLib/Polyline.h"
32 #include "InfoLib/GitInfo.h"
33 #include "MathLib/MathTools.h"
34 #include "MathLib/Point3d.h"
39 #include "MeshLib/Mesh.h"
44 #include "MeshLib/Node.h"
45 
47 std::vector<GeoLib::Point*> createPoints(MathLib::Point3d const& start,
48  MathLib::Point3d const& end,
49  std::size_t const n_intervals)
50 {
51  std::vector<GeoLib::Point*> points;
52  points.push_back(new GeoLib::Point(start, 0));
53  double const length = std::sqrt(MathLib::sqrDist(start, end));
54  double const interval = length / n_intervals;
55 
56  GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
57  (end[2] - start[2]));
58  for (std::size_t i = 1; i < n_intervals; ++i)
59  {
60  points.push_back(new GeoLib::Point(
61  start[0] + ((i * interval) / length * vec[0]),
62  start[1] + ((i * interval) / length * vec[1]), 0, i));
63  }
64  points.push_back(new GeoLib::Point(end, n_intervals));
65  return points;
66 }
67 
69 GeoLib::Polyline* createPolyline(std::vector<GeoLib::Point*> const& points)
70 {
71  GeoLib::Polyline* line = new GeoLib::Polyline(points);
72  std::size_t const length = points.size();
73  for (std::size_t i = 0; i < length; ++i)
74  {
75  line->addPoint(i);
76  }
77  return line;
78 }
79 
81 std::vector<std::string> createGeometries(
82  GeoLib::GEOObjects& geo,
83  std::vector<std::string> const& layer_names,
84  MathLib::Point3d const& pnt_start,
85  MathLib::Point3d const& pnt_end,
86  double const resolution)
87 {
88  std::vector<std::string> geo_name_list;
89  std::size_t const n_layers = layer_names.size();
90  for (std::size_t i = 0; i < n_layers; ++i)
91  {
92  std::unique_ptr<MeshLib::Mesh> const layer(
93  MeshLib::IO::readMeshFromFile(layer_names[i]));
94  if (layer == nullptr)
95  {
96  ERR("Could not read file {:s}. Skipping layer...", layer_names[i]);
97  continue;
98  }
99  if (layer->getDimension() != 2)
100  {
101  ERR("Layer {:d} is not a 2D mesh. Skipping layer...", i);
102  continue;
103  }
104 
105  std::string geo_name(std::to_string(i));
106  auto points(createPoints(pnt_start, pnt_end, resolution));
107  geo.addPointVec(std::move(points), geo_name,
109 
110  std::vector<GeoLib::Polyline*> lines;
111  lines.push_back(createPolyline(*geo.getPointVec(geo_name)));
112  geo.addPolylineVec(std::move(lines), geo_name,
114 
115  MeshGeoToolsLib::GeoMapper mapper(geo, geo_name);
116  mapper.mapOnMesh(layer.get());
117  geo_name_list.push_back(geo_name);
118  }
119  return geo_name_list;
120 }
121 
126  std::vector<std::string> const& geo_names,
127  std::string& merged_geo_name)
128 {
129  std::vector<GeoLib::Point*> points;
130  std::vector<GeoLib::Polyline*> lines;
131 
132  auto layer_pnts = *geo.getPointVec(geo_names[0]);
133  std::size_t const pnts_per_line = layer_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  for (std::size_t i = 0; i < pnts_per_line; ++i)
138  {
139  std::size_t const idx = pnts_per_line - i - 1;
140  points.push_back(new GeoLib::Point(*layer_pnts[i], idx));
141  last_line_idx[i] = idx;
142  }
143  for (std::size_t j = 1; j < n_layers; ++j)
144  {
145  GeoLib::Polyline* line = new GeoLib::Polyline(points);
146  for (std::size_t i = 0; i < pnts_per_line; ++i)
147  {
148  line->addPoint(last_line_idx[i]);
149  }
150  layer_pnts = *geo.getPointVec(geo_names[j]);
151  for (std::size_t i = 0; i < pnts_per_line; ++i)
152  {
153  // check if for current point the lower layer boundary is actually
154  // located below the upper boundary
155  std::size_t idx = last_line_idx[pnts_per_line - i - 1];
156  if ((*points[idx])[2] > (*layer_pnts[i])[2])
157  {
158  idx = points.size();
159  points.push_back(new GeoLib::Point(*layer_pnts[i], idx));
160  last_line_idx[pnts_per_line - i - 1] = idx;
161  }
162  line->addPoint(idx);
163  }
164  // close polygon
165  line->addPoint(line->getPointID(0));
166  lines.push_back(line);
167  }
168 
169  geo.addPointVec(std::move(points), merged_geo_name,
171  geo.addPolylineVec(std::move(lines), merged_geo_name,
173 }
174 
176 std::pair<Eigen::Matrix3d, double> rotateGeometryToXY(
177  std::vector<GeoLib::Point*>& points)
178 {
179  // compute the plane normal
180  auto const [plane_normal, d] =
181  GeoLib::getNewellPlane(points.begin(), points.end());
182  // rotate points into x-y-plane
183  Eigen::Matrix3d const rotation_matrix =
184  GeoLib::computeRotationMatrixToXY(plane_normal);
185  GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
186 
187  GeoLib::AABB aabb(points.begin(), points.end());
188  double const z_shift =
189  (aabb.getMinPoint()[2] + aabb.getMaxPoint()[2]) / 2.0;
190  std::for_each(points.begin(), points.end(),
191  [z_shift](GeoLib::Point* p) { (*p)[2] -= z_shift; });
192  return {rotation_matrix, z_shift};
193 }
194 
200  std::string const& output_name,
201  std::string& merged_geo_name,
202  bool const keep_gml_file)
203 {
204  std::string const filename(output_name + ".gml");
206  xml.export_name = merged_geo_name;
208 
209  geo.removePolylineVec(merged_geo_name);
210  geo.removePointVec(merged_geo_name);
211 
212  xml.readFile(filename);
213 
214  if (!keep_gml_file)
215  {
216  BaseLib::removeFile(filename);
217  BaseLib::removeFile(filename + ".md5");
218  }
219 }
220 
223  std::string const& geo_name,
224  std::string const& output_name, double res)
225 {
226  std::string const gmsh_geo_name(output_name + ".geo");
227  std::vector<std::string> gmsh_geo;
228  gmsh_geo.push_back(geo_name);
231  0, gmsh_geo, false, false);
232  gmsh_io.writePhysicalGroups(true);
233  if (!BaseLib::IO::writeStringToFile(gmsh_io.writeToString(), gmsh_geo_name))
234  {
235  ERR("Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
236  }
237 
238  std::string const gmsh_mesh_name = output_name + ".msh";
239  std::string gmsh_command = "gmsh -2 -algo meshadapt " + gmsh_geo_name;
240  gmsh_command += " -o " + gmsh_mesh_name + " -format msh22";
241  int const return_value = std::system(gmsh_command.c_str());
242  if (return_value != 0)
243  {
244  ERR("Execution of gmsh command returned non-zero status, %d",
245  return_value);
246  }
247  return FileIO::GMSH::readGMSHMesh(gmsh_mesh_name);
248 }
249 
251 void rotateMesh(MeshLib::Mesh& mesh, Eigen::Matrix3d const& rot_mat,
252  double const z_shift)
253 {
254  std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
255  std::for_each(nodes.begin(), nodes.end(),
256  [z_shift](MeshLib::Node* n) { (*n)[2] += z_shift; });
257  GeoLib::rotatePoints(rot_mat.transpose(), nodes.begin(), nodes.end());
258 }
259 
262 {
263  std::vector<std::size_t> line_idx;
264  std::vector<MeshLib::Element*> const& elems = mesh.getElements();
265  std::for_each(elems.begin(), elems.end(),
266  [&](auto e)
267  {
268  if (e->getGeomType() == MeshLib::MeshElemType::LINE)
269  {
270  line_idx.push_back(e->getID());
271  }
272  });
273  if (line_idx.size() == mesh.getNumberOfElements())
274  {
275  return nullptr;
276  }
277  return MeshLib::removeElements(mesh, line_idx, "mesh");
278 }
279 
280 void writeBoundary(MeshLib::Mesh const& mesh,
281  std::vector<std::size_t> const& idx_array,
282  std::string const& file_name)
283 {
284  std::unique_ptr<MeshLib::Mesh> const boundary(
285  MeshLib::removeElements(mesh, idx_array, "mesh"));
286  if (boundary == nullptr)
287  {
288  ERR("Error extracting boundary '{:s}'", file_name);
289  return;
290  }
291  MeshLib::IO::VtuInterface vtu(boundary.get());
292  vtu.writeToFile(file_name + ".vtu");
293 }
294 
296  std::string const& output_name,
297  MathLib::Point3d const& pnt_start)
298 {
299  double const eps = mesh.getMinEdgeLength() / 100.0;
300  std::unique_ptr<MeshLib::Mesh> boundary_mesh(
302  mesh, "bulk_node_ids", "bulk_element_ids", "bulk_face_ids"));
303 
304  auto const& elems = boundary_mesh->getElements();
305  std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
306  bottom_bound_idx;
307  Eigen::Vector2d const anchor(pnt_start[0], pnt_start[1]);
308  for (auto e : elems)
309  {
310  Eigen::Vector2d const n1((*e->getNode(0))[0], (*e->getNode(0))[1]);
311  Eigen::Vector2d const n2((*e->getNode(1))[0], (*e->getNode(1))[1]);
312  Eigen::Vector2d const dist1(n1 - anchor);
313  Eigen::Vector2d const dist2(n2 - anchor);
314  std::size_t const id = e->getID();
315  // elements located at left or right side
316  if ((dist1 - dist2).squaredNorm() < eps)
317  {
318  top_bound_idx.push_back(id);
319  bottom_bound_idx.push_back(id);
320 
321  if (dist1.squaredNorm() < eps)
322  {
323  right_bound_idx.push_back(id);
324  }
325  else
326  {
327  left_bound_idx.push_back(id);
328  }
329  continue;
330  }
331  // elements located at top or bottom
332  if (dist2.squaredNorm() < dist1.squaredNorm())
333  {
334  top_bound_idx.push_back(id);
335  left_bound_idx.push_back(id);
336  right_bound_idx.push_back(id);
337  }
338  else
339  {
340  bottom_bound_idx.push_back(id);
341  left_bound_idx.push_back(id);
342  right_bound_idx.push_back(id);
343  }
344  }
345 
346  writeBoundary(*boundary_mesh, left_bound_idx, output_name + "_left");
347  writeBoundary(*boundary_mesh, right_bound_idx, output_name + "_right");
348  writeBoundary(*boundary_mesh, top_bound_idx, output_name + "_top");
349  writeBoundary(*boundary_mesh, bottom_bound_idx, output_name + "_bottom");
350 }
351 
352 int main(int argc, char* argv[])
353 {
354  QCoreApplication a(argc, argv);
355  TCLAP::CmdLine cmd(
356  "Creates a triangle-mesh of a vertical slice out of a list of input "
357  "layer meshes. The slice is defined by a start- and end-point. In "
358  "addition, the resolution for meshing the extracted slice (i.e. the "
359  "maximum edge length of the domain discretisation) needs to be "
360  "specified. The utility requires access to the meshing utility GMSH to "
361  "work correctly.\n\n"
362  "OpenGeoSys-6 software, version " +
364  ".\n"
365  "Copyright (c) 2012-2022, OpenGeoSys Community "
366  "(http://www.opengeosys.org)",
368  TCLAP::SwitchArg test_arg("t", "testdata", "keep test data", false);
369  cmd.add(test_arg);
370  TCLAP::SwitchArg bound_arg("b", "bounds", "save mesh boundaries", false);
371  cmd.add(bound_arg);
372  TCLAP::ValueArg<double> res_arg(
373  "r", "resolution",
374  "desired edge length of triangles in the resulting slice.", true, 0,
375  "floating point number");
376  cmd.add(res_arg);
377  TCLAP::ValueArg<double> end_y_arg(
378  "", "end-y", "y-coordinates of the end point defining the slice", true,
379  0, "floating point number");
380  cmd.add(end_y_arg);
381  TCLAP::ValueArg<double> end_x_arg(
382  "", "end-x", "x-coordinates of the end point defining the slice", true,
383  0, "floating point number");
384  cmd.add(end_x_arg);
385  TCLAP::ValueArg<double> start_y_arg(
386  "", "start-y", "y-coordinates of the start point defining the slice",
387  true, 0, "floating point number");
388  cmd.add(start_y_arg);
389  TCLAP::ValueArg<double> start_x_arg(
390  "", "start-x", "x-coordinates of the start point defining the slice",
391  true, 0, "floating point number");
392  cmd.add(start_x_arg);
393  TCLAP::ValueArg<std::string> output_arg(
394  "o", "output", "name of output mesh (*.vtu)", true, "", "string");
395  cmd.add(output_arg);
396  TCLAP::ValueArg<std::string> input_arg(
397  "i", "input",
398  "name of the input file list containing the paths the all input layers "
399  "in correct order from top to bottom",
400  true, "", "string");
401  cmd.add(input_arg);
402  cmd.parse(argc, argv);
403 
404  std::string const input_name = input_arg.getValue();
405  std::string const output_name = output_arg.getValue();
406 
407  MathLib::Point3d const pnt_start{
408  {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
409  MathLib::Point3d const pnt_end{
410  {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
411  double const length = std::sqrt(MathLib::sqrDist(pnt_start, pnt_end));
412 
413  std::size_t const res = std::ceil(length / res_arg.getValue());
414  double const interval_length = length / res;
415 
416  std::vector<std::string> const layer_names =
418  if (layer_names.size() < 2)
419  {
420  ERR("At least two layers are required to extract a slice.");
421  return EXIT_FAILURE;
422  }
423 
424  GeoLib::GEOObjects geo;
425  std::vector<std::string> const geo_name_list =
426  createGeometries(geo, layer_names, pnt_start, pnt_end, res);
427 
428  if (geo_name_list.size() < 2)
429  {
430  ERR("Less than two geometries could be created from layers. Aborting "
431  "extraction...");
432  return EXIT_FAILURE;
433  }
434 
435  std::string merged_geo_name = "merged_geometries";
436  mergeGeometries(geo, geo_name_list, merged_geo_name);
437  std::vector<GeoLib::Point*> points = *geo.getPointVec(merged_geo_name);
438  auto const [rot_mat, z_shift] = rotateGeometryToXY(points);
439  consolidateGeometry(geo, output_name, merged_geo_name, test_arg.getValue());
440 
441  std::unique_ptr<MeshLib::Mesh> mesh(
442  generateMesh(geo, merged_geo_name, output_name, interval_length));
443  if (!test_arg.getValue())
444  {
445  BaseLib::removeFile(output_name + ".geo");
446  BaseLib::removeFile(output_name + ".msh");
447  }
448  if (mesh == nullptr)
449  {
450  ERR("Error generating mesh... (GMSH was unable to output mesh)");
451  return EXIT_FAILURE;
452  }
453  rotateMesh(*mesh, rot_mat, z_shift);
454  std::unique_ptr<MeshLib::Mesh> new_mesh(removeLineElements(*mesh));
455  if (new_mesh == nullptr)
456  {
457  ERR("Error generating mesh... (GMSH created line mesh)");
458  return EXIT_FAILURE;
459  }
460 
461  // collapse all nodes that might have been created due to gmsh physically
462  // separating layers
463  MeshLib::MeshRevision rev(*new_mesh);
464  std::unique_ptr<MeshLib::Mesh> revised_mesh(
465  rev.simplifyMesh("RevisedMesh", new_mesh->getMinEdgeLength() / 100.0));
466 
467  if (bound_arg.getValue())
468  {
469  extractBoundaries(*revised_mesh, output_name, pnt_start);
470  }
471  MeshLib::IO::VtuInterface vtu(revised_mesh.get());
472  vtu.writeToFile(output_name + ".vtu");
473  return EXIT_SUCCESS;
474 }
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(char const *fmt, Args const &... args)
Definition: Logging.h:42
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.
MeshLib::Mesh * removeLineElements(MeshLib::Mesh const &mesh)
removes line elements from mesh such that only triangles remain
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
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)
void writeBoundary(MeshLib::Mesh const &mesh, std::vector< std::size_t > const &idx_array, std::string const &file_name)
std::pair< Eigen::Matrix3d, double > rotateGeometryToXY(std::vector< GeoLib::Point * > &points)
rotates the merged geometry into the XY-plane
void mergeGeometries(GeoLib::GEOObjects &geo, std::vector< std::string > const &geo_names, std::string &merged_geo_name)
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
GeoLib::Polyline * createPolyline(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
void consolidateGeometry(GeoLib::GEOObjects &geo, std::string const &output_name, std::string &merged_geo_name, bool const keep_gml_file)
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.
Definition: GMSHInterface.h:46
void writePhysicalGroups(bool flag)
Definition: GMSHInterface.h:83
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition: AABB.h:49
Eigen::Vector3d const & getMinPoint() const
Definition: AABB.h:172
Eigen::Vector3d const & getMaxPoint() const
Definition: AABB.h:179
Container class for geometric objects.
Definition: GEOObjects.h:61
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
Definition: GEOObjects.cpp:152
const std::vector< Point * > * getPointVec(const std::string &name) const
Definition: GEOObjects.cpp:72
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()))
Definition: GEOObjects.cpp:47
bool removePointVec(const std::string &name)
Definition: GEOObjects.cpp:98
bool removePolylineVec(const std::string &name)
Definition: GEOObjects.cpp:222
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:53
std::size_t getPointID(std::size_t 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:44
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....
Definition: VtuInterface.h:38
bool writeToFile(std::filesystem::path const &file_path)
MeshLib::Mesh * simplifyMesh(const std::string &new_mesh_name, double eps, unsigned min_elem_dim=1) const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:95
double getMinEdgeLength() const
Get the minimum edge length over all elements of the mesh.
Definition: Mesh.h:80
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
int writeStringToFile(std::string content, std::filesystem::path const &file_path)
Definition: Writer.cpp:45
std::vector< std::string > readStringListFromFile(std::string const &filename)
Reads non-empty lines from a list of strings from a file into a vector.
void removeFile(std::string const &filename)
Definition: FileTools.cpp:236
@ FixedMeshDensity
set the parameter with a fixed value
MeshLib::Mesh * readGMSHMesh(std::string const &fname)
Definition: GmshReader.cpp:270
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
static const double p
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition: Point3d.h:48
std::unique_ptr< MeshLib::Mesh > getBoundaryElementsAsMesh(MeshLib::Mesh const &bulk_mesh, std::string const &subsfc_node_id_prop_name, std::string const &subsfc_element_id_prop_name, std::string const &face_id_prop_name)
MeshLib::Mesh * readMeshFromFile(const std::string &file_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.