OGS
GMSHInterface.cpp
Go to the documentation of this file.
1
12
13#include <fstream>
14#include <memory>
15#include <vector>
16
22#include "BaseLib/Algorithm.h"
23#include "BaseLib/FileTools.h"
24#include "BaseLib/Logging.h"
26#include "GeoLib/GEOObjects.h"
27#include "GeoLib/Polygon.h"
30#include "GeoLib/Station.h"
31#include "InfoLib/GitInfo.h"
32
33namespace FileIO
34{
35namespace GMSH
36{
37static std::ostream& operator<<(std::ostream& os,
38 std::vector<GMSHPoint*> const& points)
39{
40 for (auto const* const point : points)
41 {
42 if (point)
43 {
44 os << *point << "\n";
45 }
46 }
47 return os;
48}
49
51 GeoLib::GEOObjects& geo_objs, bool /*include_stations_as_constraints*/,
52 GMSH::MeshDensityAlgorithm const mesh_density_algorithm,
53 double const pnt_density, double const station_density,
54 std::size_t const max_pnts_per_leaf,
55 std::vector<std::string> const& selected_geometries, bool const rotate,
56 bool const keep_preprocessed_geometry)
57 : _n_lines(0),
58 _n_plane_sfc(0),
59 _geo_objs(geo_objs),
60 _selected_geometries(selected_geometries),
61 _rotate(rotate),
62 _keep_preprocessed_geometry(keep_preprocessed_geometry)
63{
64 switch (mesh_density_algorithm)
65 {
68 std::make_unique<GMSH::GMSHFixedMeshDensity>(pnt_density);
69 break;
72 std::make_unique<GMSH::GMSHAdaptiveMeshDensity>(
73 pnt_density, station_density, max_pnts_per_leaf);
74 break;
75 }
76}
77
79{
81 for (auto const* polygon_tree : _polygon_tree_list)
82 {
83 delete polygon_tree;
84 }
85}
86
88{
89 out << "// GMSH input file created by OpenGeoSys "
91 out << "\n\n";
92
93 return writeGMSHInputFile(out) <= 0;
94}
95
97{
98 DBUG("GMSHInterface::writeGMSHInputFile(): get data from GEOObjects.");
99
100 if (_selected_geometries.empty())
101 {
102 return 1;
103 }
104
105 // *** get and merge data from _geo_objs
106 if (_selected_geometries.size() > 1)
107 {
108 _gmsh_geo_name = "GMSHGeometry";
110 {
111 return 2;
112 }
113 }
114 else
115 {
118 }
119
120 auto* merged_pnts(const_cast<std::vector<GeoLib::Point*>*>(
122 if (!merged_pnts)
123 {
124 ERR("GMSHInterface::writeGMSHInputFile(): Did not found any points.");
125 return 2;
126 }
127
128 if (_rotate)
129 {
130 // Rotate points to the x-y-plane.
132 // Compute inverse rotation matrix to reverse the rotation later on.
133 _inverse_rot_mat.transposeInPlace();
134 }
135 else
136 {
137 // project data on the x-y-plane
138 _inverse_rot_mat = Eigen::Matrix3d::Identity();
139 for (auto pnt : *merged_pnts)
140 {
141 (*pnt)[2] = 0.0;
142 }
143 }
144
145 std::vector<GeoLib::Polyline*> const* merged_plys(
147 DBUG("GMSHInterface::writeGMSHInputFile(): Obtained data.");
148
149 if (!merged_plys)
150 {
151 ERR("GMSHInterface::writeGMSHInputFile(): Did not find any polylines.");
152 return 2;
153 }
154
155 // *** compute and insert all intersection points between polylines
156 GeoLib::PointVec& pnt_vec(*const_cast<GeoLib::PointVec*>(
159 pnt_vec, *(const_cast<std::vector<GeoLib::Polyline*>*>(merged_plys)));
160
161 std::vector<GeoLib::Polyline*> polygons;
162 // for each closed polyline add a PolygonWithSegmentMarker object into
163 // polygons
164 for (auto polyline : *merged_plys)
165 {
166 if (!polyline->isClosed())
167 {
168 continue;
169 }
170 polygons.push_back(new GeoLib::PolygonWithSegmentMarker(*polyline));
171 }
172 if (polygons.empty())
173 {
174 OGS_FATAL("GMSHInterface::writeGMSHInputFile(): no polygons found.");
175 }
176 // let the polygon memory management be done by GEOObjects
178 // create for each polygon a PolygonTree
179 std::transform(
180 polygons.begin(), polygons.end(),
181 std::back_inserter(_polygon_tree_list),
182 [this](auto const& polygon)
183 {
184 return new GMSH::GMSHPolygonTree(
185 dynamic_cast<GeoLib::PolygonWithSegmentMarker*>(polygon),
186 nullptr, _geo_objs, _gmsh_geo_name, *_mesh_density_strategy);
187 });
188 DBUG(
189 "GMSHInterface::writeGMSHInputFile(): Computed topological hierarchy - "
190 "detected {:d} polygons.",
191 _polygon_tree_list.size());
192 // compute topological hierarchy of polygons
193 GeoLib::createPolygonTrees<GMSH::GMSHPolygonTree>(_polygon_tree_list);
194 DBUG(
195 "GMSHInterface::writeGMSHInputFile(): Computed topological hierarchy - "
196 "calculated {:d} polygon trees.",
197 _polygon_tree_list.size());
198
199 // *** Mark in each polygon tree the segments shared by two polygons.
200 for (auto* polygon_tree : _polygon_tree_list)
201 {
202 polygon_tree->markSharedSegments();
203 }
204
205 // *** insert stations and polylines (except polygons) in the appropriate
206 // object of
207 // class GMSHPolygonTree
208 // *** insert stations
209 std::vector<GeoLib::Point*> gmsh_stations{};
210 for (auto const& geometry_name : _selected_geometries)
211 {
212 auto const* stations(_geo_objs.getStationVec(geometry_name));
213 if (stations)
214 {
215 for (auto* station : *stations)
216 {
217 bool found(false);
218 for (auto it(_polygon_tree_list.begin());
219 it != _polygon_tree_list.end() && !found;
220 ++it)
221 {
222 gmsh_stations.emplace_back(new GeoLib::Station(
223 *static_cast<GeoLib::Station*>(station)));
224 if ((*it)->insertStation(gmsh_stations.back()))
225 {
226 found = true;
227 }
228 }
229 }
230 }
231 }
232 std::string gmsh_stations_name(_gmsh_geo_name + "-Stations");
233 if (!gmsh_stations.empty())
234 {
235 _geo_objs.addStationVec(std::move(gmsh_stations), gmsh_stations_name);
236 }
237
238 // *** insert polylines
239 for (auto polyline : *merged_plys)
240 {
241 if (!polyline->isClosed())
242 {
243 for (auto* polygon_tree : _polygon_tree_list)
244 {
245 auto polyline_with_segment_marker =
247 polygon_tree->insertPolyline(polyline_with_segment_marker);
248 }
249 }
250 }
251
252 // *** init mesh density strategies
253 for (auto& polygon_tree : _polygon_tree_list)
254 {
255 polygon_tree->initMeshDensityStrategy();
256 }
257
258 // *** create GMSH data structures
259 const std::size_t n_merged_pnts(merged_pnts->size());
260 _gmsh_pnts.resize(n_merged_pnts);
261 for (std::size_t k(0); k < n_merged_pnts; k++)
262 {
263 _gmsh_pnts[k] = nullptr;
264 }
265 for (auto& polygon_tree : _polygon_tree_list)
266 {
267 polygon_tree->createGMSHPoints(_gmsh_pnts);
268 }
269
270 std::stringstream error_messages;
271 error_messages.precision(std::numeric_limits<double>::max_digits10);
272 for (std::size_t k = 0; k < _gmsh_pnts.size(); ++k)
273 {
274 if (_gmsh_pnts[k] == nullptr)
275 {
276 error_messages
277 << "The point at (" << *(*merged_pnts)[k]
278 << ") is not part of a polyline, and won't be used in the "
279 "meshing as a constraint. If you want to include it in the "
280 "mesh please create a observation/measurement station for "
281 "the point and include it additional in the meshing "
282 "process.\n";
283 }
284 }
285 auto const error_message = error_messages.str();
286 if (!error_message.empty())
287 {
288 OGS_FATAL("{}", error_message);
289 }
290
291 // *** finally write data :-)
293 out << _gmsh_pnts;
294
295 std::size_t pnt_id_offset(_gmsh_pnts.size());
296 for (auto* polygon_tree : _polygon_tree_list)
297 {
298 polygon_tree->writeLineLoop(_n_lines, _n_plane_sfc, out,
300 polygon_tree->writeSubPolygonsAsLineConstraints(_n_lines,
301 _n_plane_sfc - 1, out);
302 polygon_tree->writeLineConstraints(_n_lines, _n_plane_sfc - 1, out);
303 polygon_tree->writeStations(pnt_id_offset, _n_plane_sfc - 1, out);
304 polygon_tree->writeAdditionalPointData(pnt_id_offset, _n_plane_sfc - 1,
305 out);
306 }
307
309 {
313 _geo_objs.removeStationVec(gmsh_stations_name);
314 }
315
316 return 0;
317}
318
319} // end namespace GMSH
320} // end namespace FileIO
Definition of analytical geometry functions.
#define OGS_FATAL(...)
Definition Error.h:26
Filename manipulation routines.
Definition of the GEOObjects class.
Git information.
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the Polygon class.
Definition of the PolylineWithSegmentMarker class.
Definition of the Station class.
std::ostringstream out
The stream to write to.
Definition Writer.h:47
GMSHInterface(GeoLib::GEOObjects &geo_objs, bool include_stations_as_constraints, GMSH::MeshDensityAlgorithm mesh_density_algorithm, double pnt_density, double station_density, std::size_t max_pnts_per_leaf, std::vector< std::string > const &selected_geometries, bool rotate, bool keep_preprocessed_geometry)
std::unique_ptr< GMSH::GMSHMeshDensityStrategy > _mesh_density_strategy
Eigen::Matrix3d _inverse_rot_mat
std::vector< GMSH::GMSHPoint * > _gmsh_pnts
int writeGMSHInputFile(std::ostream &out)
std::vector< std::string > const & _selected_geometries
GeoLib::GEOObjects & _geo_objs
bool write() override
Writes the object to the internal stream. This method must be implemented by a subclass....
std::list< GMSH::GMSHPolygonTree * > _polygon_tree_list
Container class for geometric objects.
Definition GEOObjects.h:57
const std::vector< Point * > * getPointVec(const std::string &name) const
const PointVec * getPointVecObj(const std::string &name) const
bool removePointVec(const std::string &name)
bool removeSurfaceVec(const std::string &name)
void addStationVec(std::vector< Point * > &&stations, std::string &name)
Adds a vector of stations with the given name and colour to GEOObjects.
bool removeStationVec(const std::string &name)
Removes the station vector with the given name from GEOObjects.
Definition GEOObjects.h:135
int mergeGeometries(std::vector< std::string > const &geo_names, std::string &merged_geo_name)
const std::vector< Polyline * > * getPolylineVec(const std::string &name) const
bool appendPolylineVec(const std::vector< Polyline * > &polylines, const std::string &name)
const std::vector< GeoLib::Point * > * getStationVec(const std::string &name) const
Returns the station vector with the given name.
bool removePolylineVec(const std::string &name)
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition PointVec.h:36
A Station (observation site) is basically a Point with some additional information.
Definition Station.h:37
void cleanupVectorElements(std::vector< T * > &items)
Definition Algorithm.h:252
@ AdaptiveMeshDensity
computing the mesh density employing a QuadTree
@ FixedMeshDensity
set the parameter with a fixed value
static std::ostream & operator<<(std::ostream &os, std::vector< GMSHPoint * > const &points)
void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec, std::vector< GeoLib::Polyline * > &plys)
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
GITINFOLIB_EXPORT const std::string ogs_version