OGS
generateGeometry.cpp
Go to the documentation of this file.
1
11#include <tclap/CmdLine.h>
12
13#include <numeric>
14
15#include "BaseLib/Logging.h"
16#include "BaseLib/MPI.h"
18#include "GeoLib/GEOObjects.h"
20#include "GeoLib/Point.h"
21#include "GeoLib/Polyline.h"
22#include "GeoLib/Utils.h"
23#include "InfoLib/GitInfo.h"
24
25std::tuple<std::vector<GeoLib::Polyline*>, GeoLib::PolylineVec::NameIdMap>
26appendNamedPolyline(std::unique_ptr<GeoLib::Polyline> polyline,
27 std::string&& polyline_name)
28{
29 std::vector<GeoLib::Polyline*> lines;
31
32 lines.push_back(polyline.release());
33 name_map[std::move(polyline_name)] = lines.size() - 1;
34
35 return {lines, name_map};
36}
37
39 std::string&& point_name,
40 std::string&& geometry_name,
41 GeoLib::GEOObjects& geometry)
42{
43 std::vector<GeoLib::Point*> points;
44 points.push_back(new GeoLib::Point{point});
45
46 GeoLib::PointVec::NameIdMap name_map{{std::move(point_name), 0}};
47
48 geometry.addPointVec(std::move(points), geometry_name, std::move(name_map));
49}
50
52 GeoLib::Point const& point1,
53 int const number_of_subdivisions,
54 std::string&& polyline_name,
55 std::string&& geometry_name,
56 GeoLib::GEOObjects& geometry)
57{
58 auto intermediate_points = GeoLib::generateEquidistantPoints(
59 point0, point1, number_of_subdivisions);
60 std::vector<GeoLib::Point*> points(intermediate_points.begin(),
61 intermediate_points.end());
62 geometry.addPointVec(std::move(points), geometry_name,
64 auto const& point_vec = *geometry.getPointVecObj(geometry_name);
65
66 std::vector<std::size_t> polyline_point_ids(point_vec.getVector().size());
67 std::iota(begin(polyline_point_ids), end(polyline_point_ids), 0);
68
69 auto [lines, name_map] = appendNamedPolyline(
70 GeoLib::createPolyline(point_vec, std::move(polyline_point_ids)),
71 std::move(polyline_name));
72
73 geometry.addPolylineVec(std::move(lines), geometry_name,
74 std::move(name_map));
75}
76
77std::vector<GeoLib::Point*> generateQuadPoints(
78 std::array<GeoLib::Point, 4> const& points,
79 std::array<int, 4> const& number_of_subdivisions_per_edge)
80{
81 std::vector<GeoLib::Point*> quad_points;
82
83 auto addPointsOnLine = [&quad_points](auto const& begin, auto const& end,
84 auto const number_of_subdivisions)
85 {
86 auto intermediate_points = GeoLib::generateEquidistantPoints(
87 begin, end, number_of_subdivisions);
88 quad_points.insert(quad_points.end(), intermediate_points.begin(),
89 --intermediate_points.end());
90 delete intermediate_points.back(); // Release last point, other points
91 // are managed by GEOObjects.
92 };
93
94 addPointsOnLine(points[0], points[1], number_of_subdivisions_per_edge[0]);
95 addPointsOnLine(points[1], points[2], number_of_subdivisions_per_edge[1]);
96 addPointsOnLine(points[2], points[3], number_of_subdivisions_per_edge[2]);
97 addPointsOnLine(points[3], points[0], number_of_subdivisions_per_edge[3]);
98
99 return quad_points;
100}
101
103 GeoLib::Point const& point1,
104 int const number_of_subdivisions_first_x,
105 int const number_of_subdivisions_second_x,
106 int const number_of_subdivisions_first_y,
107 int const number_of_subdivisions_second_y,
108 int const number_of_subdivisions_first_z,
109 int const number_of_subdivisions_second_z,
110 std::string&& quad_name, std::string&& geometry_name,
111 GeoLib::GEOObjects& geometry)
112{
113 std::array<GeoLib::Point, 4> edge_points;
114 edge_points[0] = point0;
115 edge_points[2] = point1;
116 std::array<int, 4> number_of_subdivisions;
117 if (point0[0] != point1[0] && point0[1] != point1[1] &&
118 point0[2] == point1[2])
119 {
120 // quad in xy plane
121 edge_points[1] =
122 GeoLib::Point{point1[0], point0[1], point0[2]}; // right front
123 edge_points[3] =
124 GeoLib::Point{point0[0], point1[1], point0[2]}; // left back
125 number_of_subdivisions = {
126 number_of_subdivisions_first_x, number_of_subdivisions_first_y,
127 number_of_subdivisions_second_x, number_of_subdivisions_second_y};
128 }
129 else if (point0[0] != point1[0] && point0[1] == point1[1] &&
130 point0[2] != point1[2])
131 {
132 // quad in xz plane
133 edge_points[1] =
134 GeoLib::Point{point1[0], point1[1], point0[2]}; // lower right
135 edge_points[3] =
136 GeoLib::Point{point0[0], point0[1], point1[2]}; // upper left
137 number_of_subdivisions = {
138 number_of_subdivisions_first_x, number_of_subdivisions_first_z,
139 number_of_subdivisions_second_x, number_of_subdivisions_second_z};
140 }
141 else if (point0[0] == point1[0] && point0[1] != point1[1] &&
142 point0[2] != point1[2])
143 {
144 // quad in yz plane
145 edge_points[1] =
146 GeoLib::Point{point1[0], point1[1], point0[2]}; // lower back
147 edge_points[3] =
148 GeoLib::Point{point0[0], point0[1], point1[2]}; // upper front
149 number_of_subdivisions = {
150 number_of_subdivisions_first_y, number_of_subdivisions_first_z,
151 number_of_subdivisions_second_y, number_of_subdivisions_second_z};
152 }
153 else
154 {
155 ERR("Input coordinates don't describe an axis aligned polyline or "
156 "quad.");
157 return EXIT_FAILURE;
158 }
159
160 geometry.addPointVec(
161 generateQuadPoints(edge_points, number_of_subdivisions), geometry_name,
163 auto const& point_vec = *geometry.getPointVecObj(geometry_name);
164
165 std::vector<std::size_t> polyline_point_ids(point_vec.getVector().size() +
166 1);
167 std::iota(begin(polyline_point_ids), end(polyline_point_ids), 0);
168 polyline_point_ids.back() = polyline_point_ids.front();
169
170 auto [lines, name_map] = appendNamedPolyline(
171 GeoLib::createPolyline(point_vec, std::move(polyline_point_ids)),
172 std::move(quad_name));
173
174 geometry.addPolylineVec(std::move(lines), geometry_name,
175 std::move(name_map));
176 return EXIT_SUCCESS;
177}
178
179int main(int argc, char* argv[])
180{
181 TCLAP::CmdLine cmd(
182 "Generate point, axis parallel polyline or axis parallel quad "
183 "geometry. \n\n"
184 "OpenGeoSys-6 software, version " +
186 ".\n"
187 "Copyright (c) 2012-2025, OpenGeoSys Community "
188 "(http://www.opengeosys.org)",
190 TCLAP::ValueArg<double> z0("", "z0", "z coordinate of the first point",
191 true, 0.0, "Z0");
192 cmd.add(z0);
193 TCLAP::ValueArg<double> y0("", "y0", "y coordinate of the first point",
194 true, 0.0, "Y0");
195 cmd.add(y0);
196 TCLAP::ValueArg<double> x0("", "x0", "x coordinate of the first point",
197 true, 0.0, "X0");
198 cmd.add(x0);
199 TCLAP::ValueArg<double> z1("", "z1", "z coordinate of the first point",
200 true, 1.0, "Z1");
201 cmd.add(z1);
202 TCLAP::ValueArg<double> y1("", "y1", "y coordinate of the first point",
203 true, 1.0, "Y1");
204 cmd.add(y1);
205 TCLAP::ValueArg<double> x1("", "x1", "x coordinate of the first point",
206 true, 1.0, "X1");
207 cmd.add(x1);
208 TCLAP::ValueArg<int> nx("", "nx", "number of subdivisions in x direction",
209 false, 0, "NX");
210 cmd.add(nx);
211 TCLAP::ValueArg<int> nx1("", "nx1", "number of subdivisions in x direction",
212 false, -1, "NX1");
213 cmd.add(nx1);
214 TCLAP::ValueArg<int> ny("", "ny", "number of subdivisions in y direction",
215 false, 0, "NY");
216 cmd.add(ny);
217 TCLAP::ValueArg<int> ny1("", "ny1", "number of subdivisions in y direction",
218 false, -1, "NY1");
219 cmd.add(ny1);
220 TCLAP::ValueArg<int> nz("", "nz", "number of subdivisions in z direction",
221 false, 0, "number of subdivisions in z direction");
222 cmd.add(nz);
223 TCLAP::ValueArg<int> nz1("", "nz1", "number of subdivisions in z direction",
224 false, -1, "NZ1");
225 cmd.add(nz1);
226 TCLAP::ValueArg<std::string> geometry_name(
227 "", "geometry_name", "name of the generated geometry", false,
228 "conceptual model", "GEOMETRY_NAME");
229 cmd.add(geometry_name);
230 TCLAP::ValueArg<std::string> polyline_name(
231 "", "polyline_name", "name of the generated polyline", false,
232 "polyline", "POLYLINE_NAME");
233 cmd.add(polyline_name);
234 TCLAP::ValueArg<std::string> geo_output_arg(
235 "o", "output", "Output (.gml) geometry file", true, "", "OUTPUT_FILE");
236 cmd.add(geo_output_arg);
237 auto log_level_arg = BaseLib::makeLogLevelArg();
238 cmd.add(log_level_arg);
239 cmd.parse(argc, argv);
240
241 BaseLib::MPI::Setup mpi_setup(argc, argv);
242 BaseLib::initOGSLogger(log_level_arg.getValue());
243
244 auto const p0 = GeoLib::Point{x0.getValue(), y0.getValue(), z0.getValue()};
245 auto const p1 = GeoLib::Point{x1.getValue(), y1.getValue(), z1.getValue()};
246
247 GeoLib::GEOObjects geometry;
248 auto constexpr eps = std::numeric_limits<double>::epsilon();
249 if (p1[0] - p0[0] < eps && p1[1] - p0[1] < eps && p1[2] - p0[2] < eps)
250 {
251 generateSinglePointGeometry(p0, std::string(polyline_name.getValue()),
252 std::string(geometry_name.getValue()),
253 geometry);
254 }
255 else if ((p1[0] - p0[0] >= eps && p1[1] - p0[1] < eps &&
256 p1[2] - p0[2] < eps) ||
257 (p1[0] - p0[0] < eps && p1[1] - p0[1] >= eps &&
258 p1[2] - p0[2] < eps) ||
259 (p1[0] - p0[0] < eps && p1[1] - p0[1] < eps &&
260 p1[2] - p0[2] >= eps))
261 {
263 p0, p1, nx.getValue(), std::string(polyline_name.getValue()),
264 std::string(geometry_name.getValue()), geometry);
265 }
266 else
267 {
268 auto eval = [](int v, int v1)
269 {
270 if (v1 == -1)
271 {
272 return v;
273 }
274 else
275 {
276 return v1;
277 }
278 };
280 p0, p1, nx.getValue(), eval(nx.getValue(), nx1.getValue()),
281 ny.getValue(), eval(ny.getValue(), ny1.getValue()),
282 nz.getValue(), eval(nz.getValue(), nz1.getValue()),
283 std::string(polyline_name.getValue()),
284 std::string(geometry_name.getValue()),
285 geometry) == EXIT_FAILURE)
286 {
287 return EXIT_FAILURE;
288 }
289 }
290
292 xml.export_name = geometry.getGeometryNames()[0];
293 BaseLib::IO::writeStringToFile(xml.writeToString(),
294 geo_output_arg.getValue());
295
296 return EXIT_SUCCESS;
297}
Definition of the BoostXmlGmlInterface class.
Definition of the GEOObjects class.
Definition of the Point class.
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
Definition of the PolyLine class.
Container class for geometric objects.
Definition GEOObjects.h:57
std::vector< std::string > getGeometryNames() const
Returns the names of all geometry vectors.
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
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()))
const PointVec * getPointVecObj(const std::string &name) const
std::map< std::string, std::size_t > NameIdMap
Definition TemplateVec.h:41
int main(int argc, char *argv[])
void generateSinglePointGeometry(GeoLib::Point const &point, std::string &&point_name, std::string &&geometry_name, GeoLib::GEOObjects &geometry)
void generatePolylineGeometry(GeoLib::Point const &point0, GeoLib::Point const &point1, int const number_of_subdivisions, std::string &&polyline_name, std::string &&geometry_name, GeoLib::GEOObjects &geometry)
std::tuple< std::vector< GeoLib::Polyline * >, GeoLib::PolylineVec::NameIdMap > appendNamedPolyline(std::unique_ptr< GeoLib::Polyline > polyline, std::string &&polyline_name)
std::vector< GeoLib::Point * > generateQuadPoints(std::array< GeoLib::Point, 4 > const &points, std::array< int, 4 > const &number_of_subdivisions_per_edge)
int generateQuadGeometry(GeoLib::Point const &point0, GeoLib::Point const &point1, int const number_of_subdivisions_first_x, int const number_of_subdivisions_second_x, int const number_of_subdivisions_first_y, int const number_of_subdivisions_second_y, int const number_of_subdivisions_first_z, int const number_of_subdivisions_second_z, std::string &&quad_name, std::string &&geometry_name, GeoLib::GEOObjects &geometry)
int writeStringToFile(std::string_view content, std::filesystem::path const &file_path)
Definition Writer.cpp:45
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
std::unique_ptr< Polyline > createPolyline(GeoLib::PointVec const &points_vec, std::vector< std::size_t > &&point_ids)
Create a polyline from given point ids.
Definition Polyline.cpp:564
std::vector< GeoLib::Point * > generateEquidistantPoints(MathLib::Point3d const &begin, MathLib::Point3d const &end, int const number_of_subdivisions)
Definition Utils.cpp:18
GITINFOLIB_EXPORT const std::string ogs_version