OGS
GMSHPolygonTree.cpp
Go to the documentation of this file.
1 
11 #include "GMSHPolygonTree.h"
12 
13 #include <sstream>
14 
16 #include "GMSHFixedMeshDensity.h"
18 #include "GeoLib/GEOObjects.h"
19 #include "GeoLib/Point.h"
20 #include "GeoLib/Polygon.h"
23 #include "GeoLib/Station.h"
24 
25 namespace FileIO
26 {
27 namespace GMSH
28 {
30  GMSHPolygonTree* parent,
31  GeoLib::GEOObjects& geo_objs,
32  std::string const& geo_name,
33  GMSHMeshDensityStrategy& mesh_density_strategy)
34  : GeoLib::SimplePolygonTree(polygon, parent),
35  _geo_objs(geo_objs),
36  _geo_name(geo_name),
37  _mesh_density_strategy(mesh_density_strategy)
38 {
39 }
40 
42 {
43  // the polylines are processed also by the children, but the root is
44  // responsible to cleanup up
45  if (isRoot())
46  { // root
47  for (auto* polyline : _plys)
48  {
49  delete polyline;
50  }
51  }
52 }
53 
55 {
56  if (isRoot())
57  {
58  return;
59  }
60 
61  for (auto& child : *this)
62  {
63  std::size_t const n_pnts(child->polygon().getNumberOfPoints());
64  for (std::size_t k(1); k < n_pnts; k++)
65  {
67  polygon().getPointID(k - 1),
68  polygon().getPointID(k)))
69  {
71  .markSegment(k, true);
72  }
73  }
74  }
75 }
76 
78 {
79  if (polygon().isPntInPolygon(*station))
80  {
81  // try to insert station into the child nodes
82  for (auto* child : *this)
83  {
84  if (child->polygon().isPntInPolygon(*station))
85  {
86  bool rval(dynamic_cast<GMSHPolygonTree*>(child)->insertStation(
87  station));
88  // stop recursion if sub SimplePolygonTree is a leaf
89  if (rval && child->getNumberOfChildren() == 0)
90  {
91  _stations.push_back(station);
92  }
93  return rval;
94  }
95  }
96  // station did not fit into child nodes -> insert the station into this
97  // node
98  _stations.push_back(station);
99  return true;
100  }
101  return false;
102 }
103 
105 {
106  if (!polygon().isPartOfPolylineInPolygon(*ply))
107  {
108  return;
109  }
110 
111  // check if polyline segments are inside of the polygon, intersect the
112  // polygon or are part of the boundary of the polygon
113  for (auto* child : *this)
114  {
115  dynamic_cast<GMSHPolygonTree*>(child)->insertPolyline(ply);
116  }
117 
118  // calculate possible intersection points between the node polygon and the
119  // given polyline ply
120  // pay attention: loop bound is not fix!
121  GeoLib::Point tmp_pnt;
123  for (auto segment_it(ply->begin()); segment_it != ply->end(); ++segment_it)
124  {
125  if (ply->isSegmentMarked(segment_it.getSegmentNumber()))
126  {
127  continue;
128  }
129 
130  if (polygon().containsSegment(*segment_it))
131  {
132  ply->markSegment(segment_it.getSegmentNumber(), true);
133  continue;
134  }
135 
136  std::size_t seg_num(0);
137  GeoLib::Point intersection_pnt;
138  while (polygon().getNextIntersectionPointPolygonLine(
139  *segment_it, intersection_pnt, seg_num))
140  {
141  // insert the intersection point to point vector of GEOObjects
142  // instance
143  const std::size_t pnt_vec_size(pnt_vec.size());
144  std::size_t pnt_id(
145  pnt_vec.push_back(new GeoLib::Point(intersection_pnt)));
146  if (pnt_vec_size < pnt_vec.size())
147  { // case: new point
148  // modify the polygon
149  polygon().insertPoint(seg_num + 1, pnt_id);
150  // modify the polyline
151  ply->insertPoint(segment_it.getSegmentNumber(), pnt_id);
152  }
153  else
154  { // case: existing point
155  // check if point id is within the polygon
156  if (!polygon().isPointIDInPolyline(pnt_id))
157  {
158  polygon().insertPoint(seg_num + 1, pnt_id);
159  }
160 
161  // check if point id is in polyline
162  if (!ply->isPointIDInPolyline(pnt_id))
163  {
164  ply->insertPoint(segment_it.getSegmentNumber() + 1, pnt_id);
165  }
166  }
167 
168  std::size_t tmp_seg_num(seg_num + 1);
169  if (!polygon().getNextIntersectionPointPolygonLine(
170  *segment_it, tmp_pnt, tmp_seg_num))
171  {
172  // check a point of the segment except the end points
173  for (std::size_t i(0); i < 3; i++)
174  {
175  tmp_pnt[i] = ((*segment_it).getBeginPoint()[i] +
176  (*segment_it).getEndPoint()[i]) /
177  2;
178  }
179  if (polygon().isPntInPolygon(tmp_pnt))
180  {
181  ply->markSegment(segment_it.getSegmentNumber(), true);
182  // insert line segment as constraint
183  _gmsh_lines_for_constraints.push_back(
184  new GMSHLine((*segment_it).getBeginPoint().getID(),
185  (*segment_it).getEndPoint().getID()));
186  }
187  }
188  seg_num++;
189 
190  // check a point of the segment except the end points
191  for (std::size_t i(0); i < 3; i++)
192  {
193  tmp_pnt[i] = ((*segment_it).getBeginPoint()[i] +
194  (*segment_it).getEndPoint()[i]) /
195  2;
196  }
197 
199 
200  if (polygon().isPntInPolygon(tmp_pnt))
201  {
202  ply->markSegment(segment_it.getSegmentNumber(), true);
203  // insert line segment as constraint
204  _gmsh_lines_for_constraints.push_back(
205  new GMSHLine((*segment_it).getBeginPoint().getID(),
206  (*segment_it).getEndPoint().getID()));
207  }
208  }
209  }
210 
211  _plys.push_back(ply);
212 }
213 
216  GeoLib::Polyline::SegmentIterator const& seg_it)
217 {
218  std::size_t const ply_segment_number(seg_it.getSegmentNumber());
220  {
222  for (auto seg_it_p(p->begin()); seg_it_p != p->end(); ++seg_it_p)
223  {
224  GeoLib::Point s; // intersection point
225  if (GeoLib::lineSegmentIntersect(*seg_it, *seg_it_p, s))
226  {
227  const std::size_t pnt_vec_size(pnt_vec.size());
228  // point id of new point in GEOObjects instance
229  const std::size_t pnt_id(
230  pnt_vec.push_back(new GeoLib::Point(s)));
231  if (pnt_vec_size < pnt_vec.size())
232  { // case: new point
233  // modify polyline already in this node
234  p->insertPoint(seg_it_p.getSegmentNumber() + 1, pnt_id);
235  // modify polyline
236  ply->insertPoint(ply_segment_number + 1, pnt_id);
237  }
238  else
239  { // case: point exists already in geometry
240  // check if point is not already in polyline p
241  std::size_t const k(seg_it_p.getSegmentNumber());
242  if (p->getPointID(k) != pnt_id &&
243  p->getPointID(k + 1) != pnt_id)
244  {
245  p->insertPoint(k + 1, pnt_id);
246  }
247  // check if point is not already in polyline ply
248  if (ply->getPointID(ply_segment_number) != pnt_id &&
249  ply->getPointID(ply_segment_number + 1) != pnt_id)
250  {
251  ply->insertPoint(ply_segment_number + 1, pnt_id);
252  }
253  }
254  }
255  }
256  }
257 }
258 
260 {
261  if (auto* adaptive_mesh_density =
263  {
264  // collect points
265  std::vector<GeoLib::Point const*> pnts;
266  const std::size_t n_pnts_polygon(polygon().getNumberOfPoints());
267  for (std::size_t k(0); k < n_pnts_polygon; k++)
268  {
269  pnts.push_back(polygon().getPoint(k));
270  }
272 
273  const std::size_t n_plys(_plys.size());
274  for (std::size_t k(0); k < n_plys; k++)
275  {
276  const std::size_t n_pnts_in_kth_ply(_plys[k]->getNumberOfPoints());
277  for (std::size_t j(0); j < n_pnts_in_kth_ply; j++)
278  {
279  pnts.push_back(_plys[k]->getPoint(j));
280  }
281  }
282 
283  // give collected points to the mesh density strategy
284  adaptive_mesh_density->initialize(pnts);
285  // insert constraints
286  adaptive_mesh_density->addPoints(_stations);
287  std::vector<GeoLib::Point const*> stations;
289  adaptive_mesh_density->addPoints(stations);
290  }
291 }
292 
293 void GMSHPolygonTree::createGMSHPoints(std::vector<GMSHPoint*>& gmsh_pnts) const
294 {
295  const std::size_t n_pnts_polygon(polygon().getNumberOfPoints());
296  for (std::size_t k(0); k < n_pnts_polygon - 1; k++)
297  {
298  const std::size_t id(polygon().getPointID(k));
299  GeoLib::Point const* const pnt(polygon().getPoint(k));
300  // if this point was already part of another polyline
301  if (gmsh_pnts[id] != nullptr)
302  {
303  continue;
304  }
305  gmsh_pnts[id] = new GMSHPoint(
307  }
308 
309  const std::size_t n_plys(_plys.size());
310  std::stringstream error_messages;
311  error_messages.precision(std::numeric_limits<double>::digits10);
312  for (std::size_t k(0); k < n_plys; k++)
313  {
314  const std::size_t n_pnts_in_ply(_plys[k]->getNumberOfPoints());
315  for (std::size_t j(0); j < n_pnts_in_ply; j++)
316  {
317  if (polygon().isPntInPolygon(*(_plys[k]->getPoint(j))))
318  {
319  const std::size_t id(_plys[k]->getPointID(j));
320  // if this point was already part of another polyline
321  if (gmsh_pnts[id] != nullptr)
322  {
323  continue;
324  }
325  GeoLib::Point const* const pnt(_plys[k]->getPoint(j));
326  gmsh_pnts[id] = new GMSHPoint(
327  *pnt, id,
329  }
330  else
331  {
332  auto const& p = *(_plys[k]->getPoint(j));
333  error_messages << "\n\tpoint with id " << p.getID()
334  << " and coordinates (" << p[0] << ", " << p[1]
335  << ", " << p[2]
336  << ") is outside of the polygon.";
337  }
338  }
339  }
340  if (!parent())
341  {
342  auto const error_message = error_messages.str();
343  if (!error_message.empty())
344  {
345  OGS_FATAL(error_message);
346  }
347  }
348 
349  // walk through children
350  for (auto* child : *this)
351  {
352  dynamic_cast<GMSHPolygonTree*>(child)->createGMSHPoints(gmsh_pnts);
353  }
354 }
355 
356 void GMSHPolygonTree::writeLineLoop(std::size_t& line_offset,
357  std::size_t& sfc_offset, std::ostream& out,
358  bool const write_physical) const
359 {
360  const std::size_t n_pnts(polygon().getNumberOfPoints());
361  for (std::size_t k(1), first_pnt_id(polygon().getPointID(0)); k < n_pnts;
362  k++)
363  {
364  std::size_t const second_pnt_id = polygon().getPointID(k);
365  out << "Line(" << line_offset + k - 1 << ") = {" << first_pnt_id << ","
366  << second_pnt_id << "};\n";
367  first_pnt_id = second_pnt_id;
368  }
369  out << "Line Loop(" << line_offset + n_pnts - 1 << ") = {";
370  for (std::size_t k(0); k < n_pnts - 2; k++)
371  {
372  out << line_offset + k << ",";
373  }
374  out << line_offset + n_pnts - 2 << "};\n";
375  out << "Plane Surface(" << sfc_offset << ") = {" << line_offset + n_pnts - 1
376  << "};\n";
377  if (write_physical)
378  {
379  out << "Physical Curve(" << sfc_offset << ") = {";
380  for (std::size_t k(0); k < n_pnts - 2; k++)
381  {
382  out << line_offset + k << ",";
383  }
384  out << line_offset + n_pnts - 2 << "};\n";
385  out << "Physical Surface(" << sfc_offset << ") = {" << sfc_offset
386  << "};\n";
387  }
388  line_offset += n_pnts;
389  sfc_offset++;
390 }
391 
392 void GMSHPolygonTree::writeLineConstraints(std::size_t& line_offset,
393  std::size_t sfc_number,
394  std::ostream& out) const
395 {
396  for (auto polyline : _plys)
397  {
398  const std::size_t n_pnts(polyline->getNumberOfPoints());
399  std::size_t first_pnt_id(polyline->getPointID(0));
400  for (std::size_t k(1); k < n_pnts; k++)
401  {
402  auto const second_pnt_id = polyline->getPointID(k);
403  if (polyline->isSegmentMarked(k - 1) &&
404  polygon().isPntInPolygon(*(polyline->getPoint(k))) &&
405  !GeoLib::containsEdge(polygon(), first_pnt_id, second_pnt_id))
406  {
407  out << "Line(" << line_offset + k - 1 << ") = {" << first_pnt_id
408  << "," << second_pnt_id << "};\n";
409  out << "Line { " << line_offset + k - 1 << " } In Surface { "
410  << sfc_number << " };\n";
411  }
412  first_pnt_id = second_pnt_id;
413  }
414  line_offset += n_pnts;
415  }
416 }
417 
419  std::size_t& line_offset, std::size_t sfc_number, std::ostream& out) const
420 {
421  for (auto* child : *this)
422  {
423  dynamic_cast<GMSHPolygonTree*>(child)
424  ->writeSubPolygonsAsLineConstraints(line_offset, sfc_number, out);
425  }
426 
427  if (!isRoot())
428  {
429  const std::size_t n_pnts(polygon().getNumberOfPoints());
430  std::size_t first_pnt_id(polygon().getPointID(0));
431  for (std::size_t k(1); k < n_pnts; k++)
432  {
433  auto const second_pnt_id = polygon().getPointID(k);
434  out << "Line(" << line_offset + k - 1 << ") = {" << first_pnt_id
435  << "," << second_pnt_id << "};\n";
436  first_pnt_id = second_pnt_id;
437  out << "Line { " << line_offset + k - 1 << " } In Surface { "
438  << sfc_number << " };\n";
439  }
440  line_offset += n_pnts;
441  }
442 }
443 
444 void GMSHPolygonTree::writeStations(std::size_t& pnt_id_offset,
445  std::size_t sfc_number,
446  std::ostream& out) const
447 {
448  for (auto const* station : _stations)
449  {
450  out << "Point(" << pnt_id_offset << ") = {" << (*station)[0] << ", "
451  << (*station)[1] << ", 0.0, "
453  << "}; // Station "
454  << static_cast<GeoLib::Station const*>(station)->getName() << " \n";
455  out << "Point { " << pnt_id_offset << " } In Surface { " << sfc_number
456  << " };\n";
457  ++pnt_id_offset;
458  }
459 }
460 
461 void GMSHPolygonTree::writeAdditionalPointData(std::size_t& pnt_id_offset,
462  std::size_t sfc_number,
463  std::ostream& out) const
464 {
465  if (auto* adaptive_mesh_density =
467  {
468  std::vector<GeoLib::Point*> steiner_pnts;
469  adaptive_mesh_density->getSteinerPoints(steiner_pnts, 0);
470  const std::size_t n(steiner_pnts.size());
471  for (std::size_t k(0); k < n; k++)
472  {
473  if (polygon().isPntInPolygon(*(steiner_pnts[k])))
474  {
475  out << "Point(" << pnt_id_offset + k << ") = {"
476  << (*(steiner_pnts[k]))[0] << "," << (*(steiner_pnts[k]))[1]
477  << ", 0.0, ";
479  steiner_pnts[k])
480  << "};\n";
481  out << "Point { " << pnt_id_offset + k << " } In Surface { "
482  << sfc_number << " };\n";
483  }
484  delete steiner_pnts[k];
485  }
486  pnt_id_offset += n;
487  }
488 
489 #ifndef NDEBUG
490  if (auto* adaptive_mesh_density =
492  {
493  std::vector<GeoLib::Point*> pnts;
494  std::vector<GeoLib::Polyline*> plys;
495  adaptive_mesh_density->getQuadTreeGeometry(pnts, plys);
496  std::string quad_tree_geo("QuadTree");
497  _geo_objs.addPointVec(std::move(pnts), quad_tree_geo,
499  std::vector<std::size_t> const& id_map(
500  (_geo_objs.getPointVecObj(quad_tree_geo))->getIDMap());
501  for (std::size_t k(0); k < plys.size(); k++)
502  {
503  for (std::size_t j(0); j < plys[k]->getNumberOfPoints(); j++)
504  {
505  (plys[k])->setPointID(j, id_map[(plys[k])->getPointID(j)]);
506  }
507  }
508  _geo_objs.addPolylineVec(std::move(plys), quad_tree_geo,
510  }
511 #endif
512 }
513 
515  std::vector<GeoLib::Point const*>& pnts) const
516 {
517  for (auto const* child : *this)
518  {
519  dynamic_cast<GMSHPolygonTree const*>(child)->getPointsFromSubPolygons(
520  pnts);
521  }
522 }
523 
525  std::vector<GeoLib::Point const*>& stations) const
526 {
527  const std::size_t n_stations(_stations.size());
528  for (std::size_t k(0); k < n_stations; k++)
529  {
530  stations.push_back(_stations[k]);
531  }
532 
533  for (auto const* child : *this)
534  {
535  dynamic_cast<GMSHPolygonTree const*>(child)
536  ->getStationsInsideSubPolygons(stations);
537  }
538 }
539 
540 } // end namespace GMSH
541 } // end namespace FileIO
Definition of analytical geometry functions.
#define OGS_FATAL(...)
Definition: Error.h:26
Definition of the GEOObjects class.
Definition of the Point class.
Definition of the Polygon class.
Definition of the PolylineWithSegmentMarker class.
Definition of the Station class.
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
virtual double getMeshDensityAtStation(GeoLib::Point const *const) const =0
virtual double getMeshDensityAtPoint(GeoLib::Point const *const) const =0
void writeSubPolygonsAsLineConstraints(std::size_t &line_offset, std::size_t sfc_number, std::ostream &out) const
std::vector< GeoLib::PolylineWithSegmentMarker * > _plys
GMSHMeshDensityStrategy & _mesh_density_strategy
void insertPolyline(GeoLib::PolylineWithSegmentMarker *ply)
virtual void writeLineConstraints(std::size_t &line_offset, std::size_t sfc_number, std::ostream &out) const
std::string const & _geo_name
void writeAdditionalPointData(std::size_t &pnt_id_offset, std::size_t sfc_number, std::ostream &out) const
std::vector< GeoLib::Point const * > _stations
GMSHPolygonTree(GeoLib::PolygonWithSegmentMarker *polygon, GMSHPolygonTree *parent, GeoLib::GEOObjects &geo_objs, std::string const &geo_name, GMSHMeshDensityStrategy &mesh_density_strategy)
void getStationsInsideSubPolygons(std::vector< GeoLib::Point const * > &stations) const
void writeStations(std::size_t &pnt_id_offset, std::size_t sfc_number, std::ostream &out) const
void createGMSHPoints(std::vector< GMSHPoint * > &gmsh_pnts) const
void getPointsFromSubPolygons(std::vector< GeoLib::Point const * > &pnts) const
void checkIntersectionsSegmentExistingPolylines(GeoLib::PolylineWithSegmentMarker *ply, GeoLib::Polyline::SegmentIterator const &seg_it)
GeoLib::GEOObjects & _geo_objs
bool insertStation(GeoLib::Point const *station)
std::vector< GMSHLine * > _gmsh_lines_for_constraints
virtual void writeLineLoop(std::size_t &line_offset, std::size_t &sfc_offset, std::ostream &out, bool const write_physical) const
Container class for geometric objects.
Definition: GEOObjects.h:61
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
Definition: GEOObjects.cpp:152
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()))
Definition: GEOObjects.cpp:47
const PointVec * getPointVecObj(const std::string &name) const
Definition: GEOObjects.cpp:85
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition: PointVec.h:38
std::size_t push_back(Point *pnt)
Definition: PointVec.cpp:133
bool isPntInPolygon(const MathLib::Point3d &pnt) const
Definition: Polygon.cpp:71
bool isSegmentMarked(std::size_t seg_num) const
void markSegment(std::size_t seg_num, bool mark_val=true)
bool insertPoint(std::size_t pos, std::size_t pnt_id) override
std::size_t getSegmentNumber() const
Definition: Polyline.cpp:432
std::size_t getPointID(std::size_t i) const
Definition: Polyline.cpp:160
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition: Polyline.cpp:55
SegmentIterator begin() const
Definition: Polyline.h:191
bool isPointIDInPolyline(std::size_t pnt_id) const
Definition: Polyline.cpp:154
SegmentIterator end() const
Definition: Polyline.h:196
const SimplePolygonTree * parent() const
Polygon const & polygon() const
A Station (observation site) is basically a Point with some additional information.
Definition: Station.h:37
std::map< std::string, std::size_t > NameIdMap
Definition: TemplateVec.h:44
std::size_t size() const
Definition: TemplateVec.h:103
bool containsEdge(const Polyline &ply, std::size_t id0, std::size_t id1)
Definition: Polyline.cpp:520
bool lineSegmentIntersect(GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)
static const double s