17 #include <tclap/CmdLine.h>
19 #include <QCoreApplication>
50 std::size_t
const n_intervals)
52 auto points = std::make_unique<std::vector<GeoLib::Point*>>();
55 double const interval = length / n_intervals;
57 GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
59 for (std::size_t i = 1; i < n_intervals; ++i)
62 start[0] + ((i * interval) / length * vec[0]),
63 start[1] + ((i * interval) / length * vec[1]), 0, i));
73 std::size_t
const length = points.size();
74 for (std::size_t i = 0; i < length; ++i)
84 std::vector<std::string>
const& layer_names,
87 double const resolution)
89 std::vector<std::string> geo_name_list;
90 std::size_t
const n_layers = layer_names.size();
91 for (std::size_t i = 0; i < n_layers; ++i)
93 std::unique_ptr<MeshLib::Mesh>
const layer(
97 ERR(
"Could not read file {:s}. Skipping layer...", layer_names[i]);
100 if (layer->getDimension() != 2)
102 ERR(
"Layer {:d} is not a 2D mesh. Skipping layer...", i);
106 std::string geo_name(std::to_string(i));
107 std::unique_ptr<std::vector<GeoLib::Point*>> points(
111 auto lines = std::make_unique<std::vector<GeoLib::Polyline*>>();
117 geo_name_list.push_back(geo_name);
119 return geo_name_list;
126 std::vector<std::string>
const& geo_names,
127 std::string& merged_geo_name)
129 auto points = std::make_unique<std::vector<GeoLib::Point*>>();
130 auto lines = std::make_unique<std::vector<GeoLib::Polyline*>>();
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);
137 for (std::size_t i = 0; i < pnts_per_line; ++i)
139 std::size_t
const idx = pnts_per_line - i - 1;
141 last_line_idx[i] = idx;
143 for (std::size_t j = 1; j < n_layers; ++j)
146 for (std::size_t i = 0; i < pnts_per_line; ++i)
151 for (std::size_t i = 0; i < pnts_per_line; ++i)
155 std::size_t idx = last_line_idx[pnts_per_line - i - 1];
156 if ((*(*points)[idx])[2] > (*layer_pnts[i])[2])
158 idx = points->size();
160 last_line_idx[pnts_per_line - i - 1] = idx;
166 lines->push_back(line);
169 geo.
addPointVec(std::move(points), merged_geo_name);
175 std::vector<GeoLib::Point*>& points)
178 auto const [plane_normal, d] =
181 Eigen::Matrix3d
const rotation_matrix =
186 double const z_shift =
188 std::for_each(points.begin(), points.end(),
190 return {rotation_matrix, z_shift};
198 std::string
const& output_name,
199 std::string& merged_geo_name,
200 bool const keep_gml_file)
202 std::string
const filename(output_name +
".gml");
221 std::string
const& geo_name,
222 std::string
const& output_name,
double res)
224 std::string
const gmsh_geo_name(output_name +
".geo");
225 std::vector<std::string> gmsh_geo;
226 gmsh_geo.push_back(geo_name);
229 0, gmsh_geo,
false,
false);
233 ERR(
"Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
236 std::string
const gmsh_mesh_name = output_name +
".msh";
237 std::string gmsh_command =
"gmsh -2 -algo meshadapt " + gmsh_geo_name;
238 gmsh_command +=
" -o " + gmsh_mesh_name +
" -format msh22";
239 int const return_value = std::system(gmsh_command.c_str());
240 if (return_value != 0)
242 ERR(
"Execution of gmsh command returned non-zero status, %d",
250 double const z_shift)
252 std::vector<MeshLib::Node*>
const& nodes = mesh.
getNodes();
253 std::for_each(nodes.begin(), nodes.end(),
261 std::vector<std::size_t> line_idx;
262 std::vector<MeshLib::Element*>
const& elems = mesh.
getElements();
263 std::for_each(elems.begin(), elems.end(),
266 if (e->getGeomType() == MeshLib::MeshElemType::LINE)
268 line_idx.push_back(e->getID());
271 if (line_idx.size() == mesh.getNumberOfElements())
279 std::vector<std::size_t>
const& idx_array,
280 std::string
const& file_name)
282 std::unique_ptr<MeshLib::Mesh>
const boundary(
284 if (boundary ==
nullptr)
286 ERR(
"Error extracting boundary '{:s}'", file_name);
294 std::string
const& output_name,
298 std::unique_ptr<MeshLib::Mesh> boundary_mesh(
300 mesh,
"bulk_node_ids",
"bulk_element_ids",
"bulk_face_ids"));
302 auto const& elems = boundary_mesh->getElements();
303 std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
305 Eigen::Vector2d
const anchor(pnt_start[0], pnt_start[1]);
308 Eigen::Vector2d
const n1((*e->getNode(0))[0], (*e->getNode(0))[1]);
309 Eigen::Vector2d
const n2((*e->getNode(1))[0], (*e->getNode(1))[1]);
310 Eigen::Vector2d
const dist1(n1 - anchor);
311 Eigen::Vector2d
const dist2(n2 - anchor);
312 std::size_t
const id = e->getID();
314 if ((dist1 - dist2).squaredNorm() < eps)
316 top_bound_idx.push_back(
id);
317 bottom_bound_idx.push_back(
id);
319 if (dist1.squaredNorm() < eps)
321 right_bound_idx.push_back(
id);
325 left_bound_idx.push_back(
id);
330 if (dist2.squaredNorm() < dist1.squaredNorm())
332 top_bound_idx.push_back(
id);
333 left_bound_idx.push_back(
id);
334 right_bound_idx.push_back(
id);
338 bottom_bound_idx.push_back(
id);
339 left_bound_idx.push_back(
id);
340 right_bound_idx.push_back(
id);
344 writeBoundary(*boundary_mesh, left_bound_idx, output_name +
"_left");
345 writeBoundary(*boundary_mesh, right_bound_idx, output_name +
"_right");
346 writeBoundary(*boundary_mesh, top_bound_idx, output_name +
"_top");
347 writeBoundary(*boundary_mesh, bottom_bound_idx, output_name +
"_bottom");
350 int main(
int argc,
char* argv[])
352 QCoreApplication a(argc, argv);
354 "Creates a triangle-mesh of a vertical slice out of a list of input "
355 "layer meshes. The slice is defined by a start- and end-point. In "
356 "addition, the resolution for meshing the extracted slice (i.e. the "
357 "maximum edge length of the domain discretisation) needs to be "
358 "specified. The utility requires access to the meshing utility GMSH to "
359 "work correctly.\n\n"
360 "OpenGeoSys-6 software, version " +
363 "Copyright (c) 2012-2021, OpenGeoSys Community "
364 "(http://www.opengeosys.org)",
366 TCLAP::SwitchArg test_arg(
"t",
"testdata",
"keep test data",
false);
368 TCLAP::SwitchArg bound_arg(
"b",
"bounds",
"save mesh boundaries",
false);
370 TCLAP::ValueArg<double> res_arg(
372 "desired edge length of triangles in the resulting slice.",
true, 0,
373 "floating point number");
375 TCLAP::ValueArg<double> end_y_arg(
376 "",
"end-y",
"y-coordinates of the end point defining the slice",
true,
377 0,
"floating point number");
379 TCLAP::ValueArg<double> end_x_arg(
380 "",
"end-x",
"x-coordinates of the end point defining the slice",
true,
381 0,
"floating point number");
383 TCLAP::ValueArg<double> start_y_arg(
384 "",
"start-y",
"y-coordinates of the start point defining the slice",
385 true, 0,
"floating point number");
386 cmd.add(start_y_arg);
387 TCLAP::ValueArg<double> start_x_arg(
388 "",
"start-x",
"x-coordinates of the start point defining the slice",
389 true, 0,
"floating point number");
390 cmd.add(start_x_arg);
391 TCLAP::ValueArg<std::string> output_arg(
392 "o",
"output",
"name of output mesh (*.vtu)",
true,
"",
"string");
394 TCLAP::ValueArg<std::string> input_arg(
396 "name of the input file list containing the paths the all input layers "
397 "in correct order from top to bottom",
400 cmd.parse(argc, argv);
402 std::string
const input_name = input_arg.getValue();
403 std::string
const output_name = output_arg.getValue();
406 {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
408 {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
411 std::size_t
const res = std::ceil(length / res_arg.getValue());
412 double const interval_length = length / res;
414 std::vector<std::string>
const layer_names =
416 if (layer_names.size() < 2)
418 ERR(
"At least two layers are required to extract a slice.");
423 std::vector<std::string>
const geo_name_list =
426 if (geo_name_list.size() < 2)
428 ERR(
"Less than two geometries could be created from layers. Aborting "
433 std::string merged_geo_name =
"merged_geometries";
435 std::vector<GeoLib::Point*> points = *geo.
getPointVec(merged_geo_name);
439 std::unique_ptr<MeshLib::Mesh> mesh(
440 generateMesh(geo, merged_geo_name, output_name, interval_length));
441 if (!test_arg.getValue())
448 ERR(
"Error generating mesh... (GMSH was unable to output mesh)");
453 if (new_mesh ==
nullptr)
455 ERR(
"Error generating mesh... (GMSH created line mesh)");
462 std::unique_ptr<MeshLib::Mesh> revised_mesh(
463 rev.
simplifyMesh(
"RevisedMesh", new_mesh->getMinEdgeLength() / 100.0));
465 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(char const *fmt, Args const &... 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.
MeshLib::Mesh * removeLineElements(MeshLib::Mesh const &mesh)
removes line elements from mesh such that only triangles remain
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::unique_ptr< 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
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.
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 & getMinPoint() const
Eigen::Vector3d const & getMaxPoint() const
Container class for geometric objects.
const std::vector< Point * > * getPointVec(const std::string &name) const
bool removePointVec(const std::string &name)
void addPointVec(std::unique_ptr< std::vector< Point * >> points, std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> pnt_id_name_map=nullptr, double eps=std::sqrt(std::numeric_limits< double >::epsilon()))
void addPolylineVec(std::unique_ptr< std::vector< Polyline * >> lines, const std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> ply_names=nullptr)
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 i) const
virtual bool addPoint(std::size_t pnt_id)
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
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.
double getMinEdgeLength() const
Get the minimum edge length over all elements of the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
int writeStringToFile(std::string content, std::filesystem::path const &file_path)
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)
@ FixedMeshDensity
set the parameter with a fixed value
MeshLib::Mesh * readGMSHMesh(std::string const &fname)
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)
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.