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