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