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