OGS
GeoLib Namespace Reference

Detailed Description

This module consists of classes governing geometric objects and related algorithms.

Namespaces

namespace  IO
 

Classes

class  AABB
 Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type PNT_TYPE. More...
 
class  DuplicateGeometry
 
class  EarClippingTriangulation
 
struct  GeoObject
 
class  GEOObjects
 Container class for geometric objects. More...
 
class  Grid
 
class  LineSegment
 
class  MinimalBoundingSphere
 
struct  MinMaxPoints
 
struct  NamedRaster
 
class  OctTree
 
class  Point
 
class  PointVec
 This class manages pointers to Points in a std::vector along with a name. It also handles the deletion of points. Furthermore, typically, PointVec objects are managed by GEOObjects using the instance name for identification. For this reason PointVec must have a unique name. More...
 
class  Polygon
 
class  PolygonWithSegmentMarker
 
class  Polyline
 Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices in the point vector. A polyline consists of at least one line segment. The polyline is specified by the points of the line segments. The class Polyline stores ids of pointers to the points in the _ply_pnt_ids vector. More...
 
class  PolylineWithSegmentMarker
 
class  QuadTree
 
class  Raster
 Class Raster is used for managing raster data. More...
 
struct  RasterHeader
 Contains the relevant information when storing a geoscientific raster data. More...
 
class  SimplePolygonTree
 This class computes and stores the topological relations between polygons. Every node of the SimplePolygonTree represents a polygon. A child node c of a parent node p mean that the polygon represented by c is contained in the polygon represented by p. More...
 
class  Station
 A Station (observation site) is basically a Point with some additional information. More...
 
class  StationBorehole
 A borehole as a geometric object. More...
 
class  Surface
 A Surface is represented by Triangles. It consists of a reference to a vector of (pointers to) points (_sfc_pnts) and a vector that stores the Triangles consisting of points from _sfc_pnts. More...
 
class  SurfaceGrid
 
class  TemplateVec
 The class TemplateVec takes a unique name and manages a std::vector of pointers to data elements of type T. More...
 
class  Triangle
 Class Triangle consists of a reference to a point vector and a vector that stores the indices in the point vector. A surface is composed by triangles. The class Surface stores the position of pointers to the points of triangles in the _pnt_ids vector. More...
 

Typedefs

using PolylineVec = TemplateVec<GeoLib::Polyline>
 class PolylineVec encapsulate a std::vector of Polylines additional one can give the vector of polylines a name
 
using SurfaceVec = TemplateVec<GeoLib::Surface>
 

Enumerations

enum  Orientation { CW = -1 , COLLINEAR = 0 , CCW = 1 }
 
enum class  GEOTYPE { POINT , POLYLINE , SURFACE }
 
enum class  Location {
  LEFT , RIGHT , BEYOND , BEHIND ,
  BETWEEN , SOURCE , DESTINATION
}
 
enum class  EdgeType { TOUCHING , CROSSING , INESSENTIAL }
 

Functions

template<typename InputIterator >
std::pair< Eigen::Vector3d, double > getNewellPlane (InputIterator pnts_begin, InputIterator pnts_end)
 
template<class T_POINT >
std::pair< Eigen::Vector3d, double > getNewellPlane (const std::vector< T_POINT * > &pnts)
 
template<class T_POINT >
std::pair< Eigen::Vector3d, double > getNewellPlane (const std::vector< T_POINT > &pnts)
 
template<typename InputIterator >
void rotatePoints (Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
 
template<typename InputIterator1 , typename InputIterator2 >
Eigen::Matrix3d rotatePointsToXY (InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
 
template<typename P >
void rotatePoints (Eigen::Matrix3d const &rot_mat, std::vector< P * > const &pnts)
 
Orientation getOrientation (MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
 
Orientation getOrientationFast (MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
 
bool parallel (Eigen::Vector3d v, Eigen::Vector3d w)
 
bool lineSegmentIntersect (GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)
 
bool lineSegmentsIntersect (const GeoLib::Polyline *ply, GeoLib::Polyline::SegmentIterator &seg_it0, GeoLib::Polyline::SegmentIterator &seg_it1, GeoLib::Point &intersection_pnt)
 
void rotatePoints (Eigen::Matrix3d const &rot_mat, std::vector< GeoLib::Point * > &pnts)
 
Eigen::Matrix3d computeRotationMatrixToXY (Eigen::Vector3d const &n)
 
Eigen::Matrix3d rotatePointsToXY (std::vector< GeoLib::Point * > &pnts)
 
std::unique_ptr< GeoLib::PointtriangleLineIntersection (MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, MathLib::Point3d const &p, MathLib::Point3d const &q)
 
void computeAndInsertAllIntersectionPoints (GeoLib::PointVec &pnt_vec, std::vector< GeoLib::Polyline * > &plys)
 
std::tuple< std::vector< GeoLib::Point * >, Eigen::Vector3d > rotatePolygonPointsToXY (GeoLib::Polygon const &polygon_in)
 
std::vector< MathLib::Point3dlineSegmentIntersect2d (GeoLib::LineSegment const &ab, GeoLib::LineSegment const &cd)
 
void sortSegments (MathLib::Point3d const &seg_beg_pnt, std::vector< GeoLib::LineSegment > &sub_segments)
 
Eigen::Matrix3d compute2DRotationMatrixToX (Eigen::Vector3d const &v)
 
Eigen::Matrix3d compute3DRotationMatrixToX (Eigen::Vector3d const &v)
 
Eigen::Matrix3d rotatePointsToXY (std::vector< Point * > &pnts)
 
void computeAndInsertAllIntersectionPoints (PointVec &pnt_vec, std::vector< Polyline * > &plys)
 
void sortSegments (MathLib::Point3d const &seg_beg_pnt, std::vector< LineSegment > &sub_segments)
 
template<typename Container >
auto findVectorByName (Container const &container, std::string const &name)
 
void markUnusedPoints (GEOObjects const &geo_objects, std::string const &geo_name, std::vector< bool > &transfer_pnts)
 
int geoPointsToStations (GEOObjects &geo_objects, std::string const &geo_name, std::string &stn_name, bool const only_unused_pnts=true)
 Constructs a station-vector based on the points of a given geometry.
 
std::string convertGeoTypeToString (GEOTYPE geo_type)
 
std::ostream & operator<< (std::ostream &os, LineSegment const &s)
 
std::ostream & operator<< (std::ostream &os, std::pair< GeoLib::LineSegment const &, GeoLib::LineSegment const & > const &seg_pair)
 
bool operator== (LineSegment const &s0, LineSegment const &s1)
 
Location getLocationOfPoint (MathLib::Point3d const &source, MathLib::Point3d const &destination, MathLib::Point3d const &pnt)
 
EdgeType getEdgeType (MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &pnt)
 
std::vector< GeoLib::PointgetAllIntersectionPoints (Polygon const &polygon, GeoLib::LineSegment const &segment)
 
bool operator== (Polygon const &lhs, Polygon const &rhs)
 
void resetPointIDs (Polyline &polyline, std::vector< std::size_t > const &mapping)
 Resets the point IDs of the polyline corresponding to the mapping.
 
void markUsedPoints (Polyline const &polyline, std::vector< bool > &used_points)
 Resets the point IDs of the polyline corresponding to the mapping.
 
bool containsEdge (const Polyline &ply, std::size_t id0, std::size_t id1)
 
bool operator== (Polyline const &lhs, Polyline const &rhs)
 
bool pointsAreIdentical (const std::vector< Point * > &pnt_vec, std::size_t i, std::size_t j, double prox)
 
std::unique_ptr< PolylinecreatePolyline (GeoLib::PointVec const &points_vec, std::vector< std::size_t > &&point_ids)
 Create a polyline from given point ids.
 
template<typename POLYGONTREETYPE >
void createPolygonTrees (std::list< POLYGONTREETYPE * > &list_of_simple_polygon_hierarchies)
 
bool isStation (GeoLib::Point const *pnt)
 
bool isBorehole (GeoLib::Point const *pnt)
 
bool operator== (Surface const &lhs, Surface const &rhs)
 
void resetPointIDs (Surface &surface, std::vector< std::size_t > const &mapping)
 Resets the point IDs of the surface corresponding to the mapping.
 
void markUsedPoints (Surface const &surface, std::vector< bool > &used_points)
 
std::vector< GeoLib::Point * > generateEquidistantPoints (MathLib::Point3d const &begin, MathLib::Point3d const &end, int const number_of_subdivisions)
 

Typedef Documentation

◆ PolylineVec

class PolylineVec encapsulate a std::vector of Polylines additional one can give the vector of polylines a name

Definition at line 27 of file PolylineVec.h.

◆ SurfaceVec

Class SurfaceVec encapsulate a std::vector of Surfaces and a name.

Definition at line 28 of file SurfaceVec.h.

Enumeration Type Documentation

◆ EdgeType

enum class GeoLib::EdgeType
strong

edge classification

Enumerator
TOUCHING 

TOUCHING.

CROSSING 

CROSSING.

INESSENTIAL 

INESSENTIAL.

Definition at line 38 of file Polygon.cpp.

39{
40 TOUCHING,
41 CROSSING,
43};
@ INESSENTIAL
INESSENTIAL.
@ CROSSING
CROSSING.
@ TOUCHING
TOUCHING.

◆ GEOTYPE

enum class GeoLib::GEOTYPE
strong
Enumerator
POINT 
POLYLINE 
SURFACE 

Definition at line 22 of file GeoType.h.

◆ Location

enum class GeoLib::Location
strong
Enumerator
LEFT 
RIGHT 
BEYOND 
BEHIND 
BETWEEN 
SOURCE 
DESTINATION 

Definition at line 24 of file Polygon.cpp.

◆ Orientation

Enumerator
CW 
COLLINEAR 
CCW 

Definition at line 24 of file AnalyticalGeometry.h.

25{
26 CW = -1,
27 COLLINEAR = 0,
28 CCW = 1
29};

Function Documentation

◆ compute2DRotationMatrixToX()

Eigen::Matrix3d GeoLib::compute2DRotationMatrixToX ( Eigen::Vector3d const & v)

Computes a rotation matrix that rotates the given 2D normal vector parallel to X-axis

Parameters
va 2D normal vector to be rotated
Returns
a 3x3 rotation matrix where the upper, left, 2x2 block contains the entries necessary for the 2D rotation

Definition at line 693 of file AnalyticalGeometry.cpp.

694{
695 Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
696 const double cos_theta = v[0];
697 const double sin_theta = v[1];
698 rot_mat(0, 0) = rot_mat(1, 1) = cos_theta;
699 rot_mat(0, 1) = sin_theta;
700 rot_mat(1, 0) = -sin_theta;
701 rot_mat(2, 2) = 1.0;
702 return rot_mat;
703}

Referenced by detail::getRotationMatrixToGlobal().

◆ compute3DRotationMatrixToX()

Eigen::Matrix3d GeoLib::compute3DRotationMatrixToX ( Eigen::Vector3d const & v)

Computes a rotation matrix that rotates the given 3D normal vector parallel to X-axis.

Parameters
va 3D normal vector to be rotated
Returns
a 3x3 rotation matrix

Definition at line 705 of file AnalyticalGeometry.cpp.

706{
707 // a vector on the plane
708 Eigen::Vector3d yy = Eigen::Vector3d::Zero();
709 auto const eps = std::numeric_limits<double>::epsilon();
710 if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
711 {
712 yy[2] = 1.0;
713 }
714 else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
715 {
716 yy[0] = 1.0;
717 }
718 else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
719 {
720 yy[1] = 1.0;
721 }
722 else
723 {
724 for (unsigned i = 0; i < 3; i++)
725 {
726 if (std::abs(v[i]) > 0.0)
727 {
728 yy[i] = -v[i];
729 break;
730 }
731 }
732 }
733 // z"_vec
734 Eigen::Vector3d const zz = v.cross(yy).normalized();
735 // y"_vec
736 yy = zz.cross(v).normalized();
737
738 Eigen::Matrix3d rot_mat;
739 rot_mat.row(0) = v;
740 rot_mat.row(1) = yy;
741 rot_mat.row(2) = zz;
742 return rot_mat;
743}

Referenced by detail::getRotationMatrixToGlobal().

◆ computeAndInsertAllIntersectionPoints() [1/2]

void GeoLib::computeAndInsertAllIntersectionPoints ( PointVec & pnt_vec,
std::vector< Polyline * > & plys )

Method first computes the intersection points of line segments of Polyline objects (

See also
computeIntersectionPoints()) and pushes each intersection point in the PointVec pnt_vec. For each intersection an id is returned. This id is used to split the two intersecting straight line segments in four straight line segments.

Definition at line 400 of file AnalyticalGeometry.cpp.

402{
403 auto computeSegmentIntersections =
404 [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
405 {
406 for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
407 {
408 for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
409 {
410 GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
411 if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
412 {
413 std::size_t const id(
414 pnt_vec.push_back(new GeoLib::Point(s)));
415 poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
416 poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
417 }
418 }
419 }
420 };
421
422 for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
423 {
424 auto it1(it0);
425 ++it1;
426 for (; it1 != plys.end(); ++it1)
427 {
428 computeSegmentIntersections(*(*it0), *(*it1));
429 }
430 }
431}
std::size_t push_back(Point *pnt)
Definition PointVec.cpp:133
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:40
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition Polyline.cpp:55
SegmentIterator begin() const
Definition Polyline.h:188
SegmentIterator end() const
Definition Polyline.h:190
std::size_t size() const
Definition TemplateVec.h:99
bool lineSegmentIntersect(GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), GeoLib::Polyline::insertPoint(), lineSegmentIntersect(), GeoLib::PointVec::push_back(), and GeoLib::TemplateVec< T >::size().

Referenced by FileIO::GMSH::GMSHInterface::writeGMSHInputFile().

◆ computeAndInsertAllIntersectionPoints() [2/2]

void GeoLib::computeAndInsertAllIntersectionPoints ( PointVec & pnt_vec,
std::vector< Polyline * > & plys )

Method first computes the intersection points of line segments of Polyline objects (

See also
computeIntersectionPoints()) and pushes each intersection point in the PointVec pnt_vec. For each intersection an id is returned. This id is used to split the two intersecting straight line segments in four straight line segments.

Definition at line 400 of file AnalyticalGeometry.cpp.

402{
403 auto computeSegmentIntersections =
404 [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
405 {
406 for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
407 {
408 for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
409 {
410 GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
411 if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
412 {
413 std::size_t const id(
414 pnt_vec.push_back(new GeoLib::Point(s)));
415 poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
416 poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
417 }
418 }
419 }
420 };
421
422 for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
423 {
424 auto it1(it0);
425 ++it1;
426 for (; it1 != plys.end(); ++it1)
427 {
428 computeSegmentIntersections(*(*it0), *(*it1));
429 }
430 }
431}

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), GeoLib::Polyline::insertPoint(), lineSegmentIntersect(), GeoLib::PointVec::push_back(), and GeoLib::TemplateVec< T >::size().

Referenced by FileIO::GMSH::GMSHInterface::writeGMSHInputFile().

◆ computeRotationMatrixToXY()

Eigen::Matrix3d GeoLib::computeRotationMatrixToXY ( Eigen::Vector3d const & n)

Method computes the rotation matrix that rotates the given vector parallel to the \(z\) axis.

Parameters
nthe (3d) vector that is rotated parallel to the \(z\) axis
Returns
rot_mat 3x3 rotation matrix

Definition at line 297 of file AnalyticalGeometry.cpp.

298{
299 Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
300 // check if normal points already in the right direction
301 if (n[0] == 0 && n[1] == 0)
302 {
303 rot_mat(1, 1) = 1.0;
304
305 if (n[2] > 0)
306 {
307 // identity matrix
308 rot_mat(0, 0) = 1.0;
309 rot_mat(2, 2) = 1.0;
310 }
311 else
312 {
313 // rotate by pi about the y-axis
314 rot_mat(0, 0) = -1.0;
315 rot_mat(2, 2) = -1.0;
316 }
317
318 return rot_mat;
319 }
320
321 // sqrt (n_1^2 + n_3^2)
322 double const h0(std::sqrt(n[0] * n[0] + n[2] * n[2]));
323
324 // In case the x and z components of the normal are both zero the rotation
325 // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
326 // required. The angle is either pi/2 or 3/2*pi. Thus the components of
327 // rot_mat are as follows.
328 if (h0 < std::numeric_limits<double>::epsilon())
329 {
330 rot_mat(0, 0) = 1.0;
331 if (n[1] > 0)
332 {
333 rot_mat(1, 2) = -1.0;
334 rot_mat(2, 1) = 1.0;
335 }
336 else
337 {
338 rot_mat(1, 2) = 1.0;
339 rot_mat(2, 1) = -1.0;
340 }
341 return rot_mat;
342 }
343
344 double const h1(1 / n.norm());
345
346 // general case: calculate entries of rotation matrix
347 rot_mat(0, 0) = n[2] / h0;
348 rot_mat(0, 1) = 0;
349 rot_mat(0, 2) = -n[0] / h0;
350 rot_mat(1, 0) = -n[1] * n[0] / h0 * h1;
351 rot_mat(1, 1) = h0 * h1;
352 rot_mat(1, 2) = -n[1] * n[2] / h0 * h1;
353 rot_mat(2, 0) = n[0] * h1;
354 rot_mat(2, 1) = n[1] * h1;
355 rot_mat(2, 2) = n[2] * h1;
356
357 return rot_mat;
358}

Referenced by detail::getRotationMatrixToGlobal(), MeshGeoToolsLib::markNodesOutSideOfPolygon(), rotateGeometryToXY(), rotatePointsToXY(), and rotatePolygonPointsToXY().

◆ containsEdge()

bool GeoLib::containsEdge ( const Polyline & ply,
std::size_t id0,
std::size_t id1 )

Definition at line 505 of file Polyline.cpp.

506{
507 if (id0 == id1)
508 {
509 ERR("no valid edge id0 == id1 == {:d}.", id0);
510 return false;
511 }
512 if (id0 > id1)
513 {
514 std::swap(id0, id1);
515 }
516 const std::size_t n(ply.getNumberOfPoints() - 1);
517 for (std::size_t k(0); k < n; k++)
518 {
519 std::size_t ply_pnt0(ply.getPointID(k));
520 std::size_t ply_pnt1(ply.getPointID(k + 1));
521 if (ply_pnt0 > ply_pnt1)
522 {
523 std::swap(ply_pnt0, ply_pnt1);
524 }
525 if (ply_pnt0 == id0 && ply_pnt1 == id1)
526 {
527 return true;
528 }
529 }
530 return false;
531}
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
std::size_t getPointID(std::size_t const i) const
Definition Polyline.cpp:160
std::size_t getNumberOfPoints() const
Definition Polyline.cpp:109

References ERR(), GeoLib::Polyline::getNumberOfPoints(), and GeoLib::Polyline::getPointID().

Referenced by FileIO::GMSH::GMSHPolygonTree::markSharedSegments(), and FileIO::GMSH::GMSHPolygonTree::writeLineConstraints().

◆ convertGeoTypeToString()

std::string GeoLib::convertGeoTypeToString ( GEOTYPE geo_type)

Definition at line 23 of file GeoType.cpp.

24{
25 switch (geo_type)
26 {
27 case GEOTYPE::POINT:
28 return "POINT";
29 case GEOTYPE::POLYLINE:
30 return "POLYLINE";
31 case GEOTYPE::SURFACE:
32 return "SURFACE";
33 }
34
35 // Cannot happen, because switch covers all cases.
36 // Used to silence compiler warning.
37 OGS_FATAL("convertGeoTypeToString(): Given geo type is not supported");
38}
#define OGS_FATAL(...)
Definition Error.h:26

References OGS_FATAL, POINT, POLYLINE, and SURFACE.

Referenced by GEOModels::addNameForElement(), GEOModels::addNameForObjectPoints(), GeoLib::GEOObjects::getGeoObject(), and MainWindow::showGeoNameDialog().

◆ createPolygonTrees()

template<typename POLYGONTREETYPE >
void GeoLib::createPolygonTrees ( std::list< POLYGONTREETYPE * > & list_of_simple_polygon_hierarchies)

creates from a list of simple polygons a list of trees (SimplePolygonTrees)

Definition at line 90 of file SimplePolygonTree.h.

92{
93 for (auto it0 = list_of_simple_polygon_hierarchies.begin();
94 it0 != list_of_simple_polygon_hierarchies.end();
95 ++it0)
96 {
97 auto it1 = list_of_simple_polygon_hierarchies.begin();
98 while (it1 != list_of_simple_polygon_hierarchies.end())
99 {
100 if (it0 == it1)
101 { // don't check same polygons
102 ++it1;
103 // skip test if it1 points to the end after increment
104 if (it1 == list_of_simple_polygon_hierarchies.end())
105 {
106 break;
107 }
108 }
109 if ((*it0)->isPolygonInside(*it1))
110 {
111 (*it0)->insertSimplePolygonTree(*it1);
112 it1 = list_of_simple_polygon_hierarchies.erase(it1);
113 }
114 else
115 {
116 ++it1;
117 }
118 }
119 }
120}

Referenced by FileIO::GMSH::GMSHInterface::writeGMSHInputFile().

◆ createPolyline()

std::unique_ptr< Polyline > GeoLib::createPolyline ( GeoLib::PointVec const & points_vec,
std::vector< std::size_t > && point_ids )

Create a polyline from given point ids.

Definition at line 564 of file Polyline.cpp.

566{
567 auto const& point_id_map = points_vec.getIDMap();
568 auto polyline = std::make_unique<Polyline>(points_vec.getVector());
569 for (auto point_id : point_ids)
570 {
571 if (point_id >= point_id_map.size())
572 {
573 WARN("The point id {} doesn't exist in the underlying PointVec.",
574 point_id);
575 continue;
576 }
577 polyline->addPoint(point_id_map[point_id]);
578 }
579 return polyline;
580}
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40

References GeoLib::PointVec::getIDMap(), GeoLib::TemplateVec< T >::getVector(), and WARN().

Referenced by generatePolylineGeometry(), and generateQuadGeometry().

◆ findVectorByName()

template<typename Container >
auto GeoLib::findVectorByName ( Container const & container,
std::string const & name )

Function to find PointVec, PolylineVec, or SurfaceVec pointers that are stored in the container

Definition at line 28 of file GEOObjects.cpp.

29{
30 return std::find_if(container.begin(), container.end(),
31 [&name](auto const* vector)
32 { return vector->getName() == name; });
33}

Referenced by GeoLib::GEOObjects::appendPolylineVec(), GeoLib::GEOObjects::appendSurfaceVec(), GeoLib::GEOObjects::exists(), GeoLib::GEOObjects::getSurfaceVec(), GeoLib::GEOObjects::getSurfaceVecObj(), GeoLib::GEOObjects::isPntVecUsed(), GeoLib::GEOObjects::removeSurfaceVec(), and GeoLib::GEOObjects::renameGeometry().

◆ generateEquidistantPoints()

std::vector< GeoLib::Point * > GeoLib::generateEquidistantPoints ( MathLib::Point3d const & begin,
MathLib::Point3d const & end,
int const number_of_subdivisions )

Function generates equidistant points in the interval [begin, end] according to the given number of subdivisions. In case of zero subdivisions a vector containing the begin and the end point is returned.

Definition at line 18 of file Utils.cpp.

21{
22 if (number_of_subdivisions < 0)
23 {
25 "generateEquidistantPoints: number of subdivisions is required to "
26 "be non-negative.");
27 }
28
29 auto const& start = begin.asEigenVector3d();
30 auto const& stop = end.asEigenVector3d();
31 auto const delta = (stop - start) / (number_of_subdivisions + 1);
32
33 std::vector<GeoLib::Point*> points;
34
35 for (int i = 0; i <= number_of_subdivisions; ++i)
36 {
37 auto const p = start + i * delta;
38 points.push_back(new GeoLib::Point{p[0], p[1], p[2]});
39 }
40 points.push_back(new GeoLib::Point{stop[0], stop[1], stop[2]});
41
42 return points;
43}

References MathLib::Point3d::asEigenVector3d(), and OGS_FATAL.

Referenced by generatePolylineGeometry(), and generateQuadPoints().

◆ geoPointsToStations()

int GeoLib::geoPointsToStations ( GEOObjects & geo_objects,
std::string const & geo_name,
std::string & stn_name,
bool const only_unused_pnts )

Constructs a station-vector based on the points of a given geometry.

Definition at line 611 of file GEOObjects.cpp.

613{
614 GeoLib::PointVec const* const pnt_obj(geo_objects.getPointVecObj(geo_name));
615 if (pnt_obj == nullptr)
616 {
617 ERR("Point vector {:s} not found.", geo_name);
618 return -1;
619 }
620 auto const& pnts = pnt_obj->getVector();
621 if (pnts.empty())
622 {
623 ERR("Point vector {:s} is empty.", geo_name);
624 return -1;
625 }
626 std::size_t const n_pnts(pnts.size());
627 std::vector<bool> transfer_pnts(n_pnts, true);
628 if (only_unused_pnts)
629 {
630 markUnusedPoints(geo_objects, geo_name, transfer_pnts);
631 }
632
633 std::vector<GeoLib::Point*> stations;
634 for (std::size_t i = 0; i < n_pnts; ++i)
635 {
636 if (!transfer_pnts[i])
637 {
638 continue;
639 }
640 std::string name = pnt_obj->getItemNameByID(i);
641 if (name.empty())
642 {
643 name = "Station " + std::to_string(i);
644 }
645 stations.push_back(new GeoLib::Station(pnts[i], name));
646 }
647 if (!stations.empty())
648 {
649 geo_objects.addStationVec(std::move(stations), stn_name);
650 }
651 else
652 {
653 WARN("No points found to convert.");
654 return 1;
655 }
656 return 0;
657}
const PointVec * getPointVecObj(const std::string &name) const
void addStationVec(std::vector< Point * > &&stations, std::string &name)
Adds a vector of stations with the given name and colour to GEOObjects.
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 markUnusedPoints(GEOObjects const &geo_objects, std::string const &geo_name, std::vector< bool > &transfer_pnts)

References GeoLib::GEOObjects::addStationVec(), ERR(), GeoLib::PointVec::getItemNameByID(), GeoLib::GEOObjects::getPointVecObj(), GeoLib::TemplateVec< T >::getVector(), markUnusedPoints(), and WARN().

Referenced by MainWindow::convertPointsToStations().

◆ getAllIntersectionPoints()

std::vector< GeoLib::Point > GeoLib::getAllIntersectionPoints ( Polygon const & polygon,
GeoLib::LineSegment const & segment )

Computes all intersections of the straight line segment and the polyline boundary

Parameters
polygonthe polygon the segment line segment that will be processed
segmentthe line segment that will be processed
Returns
a possible empty vector containing the intersection points

Definition at line 189 of file Polygon.cpp.

191{
192 std::vector<GeoLib::Point> intersections;
194 for (auto&& seg_it : polygon)
195 {
196 if (GeoLib::lineSegmentIntersect(seg_it, segment, s))
197 {
198 intersections.push_back(s);
199 }
200 }
201
202 return intersections;
203}

References lineSegmentIntersect().

Referenced by GeoLib::Polygon::containsSegment().

◆ getEdgeType()

EdgeType GeoLib::getEdgeType ( MathLib::Point3d const & a,
MathLib::Point3d const & b,
MathLib::Point3d const & pnt )

from book: Computational Geometry and Computer Graphics in C++, page 119 get the type of edge with respect to the given point (2d method!)

Parameters
afirst point of line segment
blast point of line segment
pntpoint that is edge type computed for
Returns
a value of enum EdgeType

Definition at line 102 of file Polygon.cpp.

105{
106 switch (getLocationOfPoint(a, b, pnt))
107 {
108 case Location::LEFT:
109 {
110 if (a[1] < pnt[1] && pnt[1] <= b[1])
111 {
112 return EdgeType::CROSSING;
113 }
114
115 return EdgeType::INESSENTIAL;
116 }
117 case Location::RIGHT:
118 {
119 if (b[1] < pnt[1] && pnt[1] <= a[1])
120 {
121 return EdgeType::CROSSING;
122 }
123
124 return EdgeType::INESSENTIAL;
125 }
126 case Location::BETWEEN:
127 case Location::SOURCE:
128 case Location::DESTINATION:
129 return EdgeType::TOUCHING;
130 default:
131 return EdgeType::INESSENTIAL;
132 }
133}
Location getLocationOfPoint(MathLib::Point3d const &source, MathLib::Point3d const &destination, MathLib::Point3d const &pnt)
Definition Polygon.cpp:55

References BETWEEN, CROSSING, DESTINATION, getLocationOfPoint(), INESSENTIAL, LEFT, RIGHT, SOURCE, and TOUCHING.

Referenced by GeoLib::Polygon::isPntInPolygon().

◆ getLocationOfPoint()

Location GeoLib::getLocationOfPoint ( MathLib::Point3d const & source,
MathLib::Point3d const & destination,
MathLib::Point3d const & pnt )

2D method - ignores z coordinate. It calculates the location of the point relative to the line segment specified by the points source and destination. (literature reference: Computational Geometry and Computer Graphics in C++; Michael J. Laszlo)

Parameters
sourcethe first point of the line segment
destinationthe end point of the line segment
pntthe test point
Returns
a value of enum Location

Definition at line 55 of file Polygon.cpp.

58{
59 long double const a[2] = {destination[0] - source[0],
60 destination[1] - source[1]}; // vector
61 long double const b[2] = {pnt[0] - source[0],
62 pnt[1] - source[1]}; // vector
63
64 long double const det_2x2(a[0] * b[1] - a[1] * b[0]);
65 constexpr double eps = std::numeric_limits<double>::epsilon();
66
67 if (det_2x2 > eps)
68 {
69 return Location::LEFT;
70 }
71 if (eps < std::abs(det_2x2))
72 {
73 return Location::RIGHT;
74 }
75 if (a[0] * b[0] < 0.0 || a[1] * b[1] < 0.0)
76 {
77 return Location::BEHIND;
78 }
79 if (a[0] * a[0] + a[1] * a[1] < b[0] * b[0] + b[1] * b[1])
80 {
81 return Location::BEYOND;
82 }
83 if (MathLib::sqrDist(pnt, source) < pow(eps, 2))
84 {
85 return Location::SOURCE;
86 }
87 if (MathLib::sqrDist(pnt, destination) < std::sqrt(eps))
88 {
89 return Location::DESTINATION;
90 }
91 return Location::BETWEEN;
92}
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26

References BEHIND, BETWEEN, BEYOND, DESTINATION, LEFT, RIGHT, SOURCE, and MathLib::sqrDist().

Referenced by getEdgeType().

◆ getNewellPlane() [1/3]

template<class T_POINT >
std::pair< Eigen::Vector3d, double > GeoLib::getNewellPlane ( const std::vector< T_POINT * > & pnts)

Compute a supporting plane (represented by plane_normal and the value d) for the polygon. Let \(n\) be the plane normal and \(d\) a parameter. Then for all points \(p \in R^3\) of the plane it holds \( n \cdot p + d = 0\). The Newell algorithm is described in [14] .

Parameters
pntspoints of a closed polyline describing a polygon
Returns
pair of plane_normal and the parameter d: the normal of the plane the polygon is located in, parameter d from the plane equation

Definition at line 46 of file AnalyticalGeometry-impl.h.

48{
49 return getNewellPlane(pnts.begin(), pnts.end());
50}
std::pair< Eigen::Vector3d, double > getNewellPlane(InputIterator pnts_begin, InputIterator pnts_end)

References getNewellPlane().

◆ getNewellPlane() [2/3]

template<class T_POINT >
std::pair< Eigen::Vector3d, double > GeoLib::getNewellPlane ( const std::vector< T_POINT > & pnts)

Same as getNewellPlane(pnts).

Definition at line 53 of file AnalyticalGeometry-impl.h.

55{
56 Eigen::Vector3d plane_normal({0, 0, 0});
57 Eigen::Vector3d centroid({0, 0, 0});
58 std::size_t n_pnts(pnts.size());
59 for (std::size_t i = n_pnts - 1, j = 0; j < n_pnts; i = j, j++)
60 {
61 plane_normal[0] += (pnts[i][1] - pnts[j][1]) *
62 (pnts[i][2] + pnts[j][2]); // projection on yz
63 plane_normal[1] += (pnts[i][2] - pnts[j][2]) *
64 (pnts[i][0] + pnts[j][0]); // projection on xz
65 plane_normal[2] += (pnts[i][0] - pnts[j][0]) *
66 (pnts[i][1] + pnts[j][1]); // projection on xy
67
68 centroid += pnts[j].asEigenVector3d();
69 }
70
71 plane_normal.normalize();
72 double const d = centroid.dot(plane_normal) / n_pnts;
73 return std::make_pair(plane_normal, d);
74}

◆ getNewellPlane() [3/3]

template<typename InputIterator >
std::pair< Eigen::Vector3d, double > GeoLib::getNewellPlane ( InputIterator pnts_begin,
InputIterator pnts_end )

Compute a supporting plane (represented by plane_normal and the value d) for the polygon. Let \(n\) be the plane normal and \(d\) a parameter. Then for all points \(p \in R^3\) of the plane it holds \( n \cdot p + d = 0\). The Newell algorithm is described in [14] .

Parameters
pnts_beginIterator pointing to the initial point of a closed polyline describing a polygon
pnts_endIterator pointing to the element following the last point of a closed polyline describing a polygon
Returns
pair of plane_normal and the parameter d: the normal of the plane the polygon is located in, d parameter from the plane equation

Definition at line 14 of file AnalyticalGeometry-impl.h.

16{
17 Eigen::Vector3d plane_normal({0, 0, 0});
18 Eigen::Vector3d centroid({0, 0, 0});
19 for (auto i = std::prev(pnts_end), j = pnts_begin; j != pnts_end;
20 i = j, ++j)
21 {
22 auto const& pt_i = *(*i);
23 auto const& pt_j = *(*j);
24 plane_normal[0] +=
25 (pt_i[1] - pt_j[1]) * (pt_i[2] + pt_j[2]); // projection on yz
26 plane_normal[1] +=
27 (pt_i[2] - pt_j[2]) * (pt_i[0] + pt_j[0]); // projection on xz
28 plane_normal[2] +=
29 (pt_i[0] - pt_j[0]) * (pt_i[1] + pt_j[1]); // projection on xy
30
31 centroid += pt_j.asEigenVector3d();
32 }
33
34 plane_normal.normalize();
35 auto n_pnts(std::distance(pnts_begin, pnts_end));
36 assert(n_pnts > 2);
37 double d = 0.0;
38 if (n_pnts > 0)
39 {
40 d = centroid.dot(plane_normal) / n_pnts;
41 }
42 return std::make_pair(plane_normal, d);
43}

Referenced by getNewellPlane(), detail::getRotationMatrixToGlobal(), rotateGeometryToXY(), rotatePointsToXY(), and rotatePolygonPointsToXY().

◆ getOrientation()

Orientation GeoLib::getOrientation ( MathLib::Point3d const & p0,
MathLib::Point3d const & p1,
MathLib::Point3d const & p2 )

Computes the orientation of the three 2D-Points. This is a robust method.

Returns
CW (clockwise), CCW (counterclockwise) or COLLINEAR (points are on a line)

Definition at line 50 of file AnalyticalGeometry.cpp.

53{
54 double const orientation = ExactPredicates::getOrientation2d(p0, p1, p2);
55 if (orientation > 0)
56 {
57 return CCW;
58 }
59 if (orientation < 0)
60 {
61 return CW;
62 }
63 return COLLINEAR;
64}
double getOrientation2d(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)

References CCW, COLLINEAR, CW, and ExactPredicates::getOrientation2d().

Referenced by GeoLib::EarClippingTriangulation::addLastTriangle(), GeoLib::EarClippingTriangulation::clipEars(), GeoLib::Polygon::ensureCCWOrientation(), GeoLib::EarClippingTriangulation::ensureCWOrientation(), MeshToolsLib::ProjectPointOnMesh::getProjectedElement(), GeoLib::EarClippingTriangulation::initLists(), and lineSegmentIntersect2d().

◆ getOrientationFast()

Orientation GeoLib::getOrientationFast ( MathLib::Point3d const & p0,
MathLib::Point3d const & p1,
MathLib::Point3d const & p2 )

Computes the orientation of the three 2D-Points. This is a non-robust method.

Returns
CW (clockwise), CCW (counterclockwise) or COLLINEAR (points are on a line)

Definition at line 66 of file AnalyticalGeometry.cpp.

69{
70 double const orientation =
72 if (orientation > 0)
73 {
74 return CCW;
75 }
76 if (orientation < 0)
77 {
78 return CW;
79 }
80 return COLLINEAR;
81}
double getOrientation2dFast(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)

References CCW, COLLINEAR, CW, and ExactPredicates::getOrientation2dFast().

◆ isBorehole()

bool GeoLib::isBorehole ( GeoLib::Point const * pnt)

Definition at line 106 of file StationBorehole.cpp.

107{
108 auto const* bh(dynamic_cast<GeoLib::StationBorehole const*>(pnt));
109 return bh != nullptr;
110}
A borehole as a geometric object.

Referenced by MeshGeoToolsLib::GeoMapper::mapStationData().

◆ isStation()

bool GeoLib::isStation ( GeoLib::Point const * pnt)

Definition at line 77 of file Station.cpp.

78{
79 auto const* bh(dynamic_cast<GeoLib::Station const*>(pnt));
80 return bh != nullptr;
81}

Referenced by MeshGeoToolsLib::GeoMapper::mapOnDEM(), and MeshGeoToolsLib::GeoMapper::mapOnMesh().

◆ lineSegmentIntersect()

bool GeoLib::lineSegmentIntersect ( LineSegment const & s0,
LineSegment const & s1,
Point & s )

A line segment is given by its two end-points. The function checks, if the two line segments (ab) and (cd) intersects.

Parameters
s0the first line segment.
s1the second line segment.
sthe intersection point if the segments do intersect
Returns
true, if the line segments intersect, else false

Definition at line 139 of file AnalyticalGeometry.cpp.

142{
143 GeoLib::Point const& pa{s0.getBeginPoint()};
144 GeoLib::Point const& pb{s0.getEndPoint()};
145 GeoLib::Point const& pc{s1.getBeginPoint()};
146 GeoLib::Point const& pd{s1.getEndPoint()};
147
148 if (!isCoplanar(pa, pb, pc, pd))
149 {
150 return false;
151 }
152
153 auto const& a = s0.getBeginPoint().asEigenVector3d();
154 auto const& b = s0.getEndPoint().asEigenVector3d();
155 auto const& c = s1.getBeginPoint().asEigenVector3d();
156 auto const& d = s1.getEndPoint().asEigenVector3d();
157
158 Eigen::Vector3d const v = b - a;
159 Eigen::Vector3d const w = d - c;
160 Eigen::Vector3d const qp = c - a;
161 Eigen::Vector3d const pq = a - c;
162
163 double const eps = std::numeric_limits<double>::epsilon();
164 double const squared_eps = eps * eps;
165 // handle special cases here to avoid computing intersection numerical
166 if (qp.squaredNorm() < squared_eps || (d - a).squaredNorm() < squared_eps)
167 {
168 s = pa;
169 return true;
170 }
171 if ((c - b).squaredNorm() < squared_eps ||
172 (d - b).squaredNorm() < squared_eps)
173 {
174 s = pb;
175 return true;
176 }
177
178 auto isLineSegmentIntersectingAB =
179 [&v](Eigen::Vector3d const& ap, std::size_t i)
180 {
181 // check if p is located at v=(a,b): (ap = t*v, t in [0,1])
182 return 0.0 <= ap[i] / v[i] && ap[i] / v[i] <= 1.0;
183 };
184
185 if (parallel(v, w))
186 { // original line segments (a,b) and (c,d) are parallel
187 if (parallel(pq, v))
188 { // line segment (a,b) and (a,c) are also parallel
189 // Here it is already checked that the line segments (a,b) and (c,d)
190 // are parallel. At this point it is also known that the line
191 // segment (a,c) is also parallel to (a,b). In that case it is
192 // possible to express c as c(t) = a + t * (b-a) (analog for the
193 // point d). Since the evaluation of all three coordinate equations
194 // (x,y,z) have to lead to the same solution for the parameter t it
195 // is sufficient to evaluate t only once.
196
197 // Search id of coordinate with largest absolute value which is will
198 // be used in the subsequent computations. This prevents division by
199 // zero in case the line segments are parallel to one of the
200 // coordinate axis.
201 std::size_t i_max(std::abs(v[0]) <= std::abs(v[1]) ? 1 : 0);
202 i_max = std::abs(v[i_max]) <= std::abs(v[2]) ? 2 : i_max;
203 if (isLineSegmentIntersectingAB(qp, i_max))
204 {
205 s = pc;
206 return true;
207 }
208 Eigen::Vector3d const ad = d - a;
209 if (isLineSegmentIntersectingAB(ad, i_max))
210 {
211 s = pd;
212 return true;
213 }
214 return false;
215 }
216 return false;
217 }
218
219 // general case
220 const double sqr_len_v(v.squaredNorm());
221 const double sqr_len_w(w.squaredNorm());
222
223 Eigen::Matrix2d mat;
224 mat(0, 0) = sqr_len_v;
225 mat(0, 1) = -v.dot(w);
226 mat(1, 1) = sqr_len_w;
227 mat(1, 0) = mat(0, 1);
228
229 Eigen::Vector2d rhs{v.dot(qp), w.dot(pq)};
230
231 rhs = mat.partialPivLu().solve(rhs);
232
233 // no theory for the following tolerances, determined by testing
234 // lower tolerance: little bit smaller than zero
235 const double l(-1.0 * std::numeric_limits<float>::epsilon());
236 // upper tolerance a little bit greater than one
237 const double u(1.0 + std::numeric_limits<float>::epsilon());
238 if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1])
239 {
240 return false;
241 }
242
243 // compute points along line segments with minimal distance
244 GeoLib::Point const p0(a[0] + rhs[0] * v[0], a[1] + rhs[0] * v[1],
245 a[2] + rhs[0] * v[2]);
246 GeoLib::Point const p1(c[0] + rhs[1] * w[0], c[1] + rhs[1] * w[1],
247 c[2] + rhs[1] * w[2]);
248
249 double const min_dist(std::sqrt(MathLib::sqrDist(p0, p1)));
250 double const min_seg_len(
251 std::min(std::sqrt(sqr_len_v), std::sqrt(sqr_len_w)));
252 if (min_dist < min_seg_len * 1e-6)
253 {
254 s[0] = 0.5 * (p0[0] + p1[0]);
255 s[1] = 0.5 * (p0[1] + p1[1]);
256 s[2] = 0.5 * (p0[2] + p1[2]);
257 return true;
258 }
259
260 return false;
261}
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
bool parallel(Eigen::Vector3d v, Eigen::Vector3d w)
bool isCoplanar(const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane.

References MathLib::Point3d::asEigenVector3d(), GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), parallel(), and MathLib::sqrDist().

Referenced by FileIO::GMSH::GMSHPolygonTree::checkIntersectionsSegmentExistingPolylines(), computeAndInsertAllIntersectionPoints(), getAllIntersectionPoints(), GeoLib::Polygon::getNextIntersectionPointPolygonLine(), and lineSegmentsIntersect().

◆ lineSegmentIntersect2d()

std::vector< MathLib::Point3d > GeoLib::lineSegmentIntersect2d ( LineSegment const & ab,
LineSegment const & cd )

A line segment is given by its two end-points. The function checks, if the two line segments (ab) and (cd) intersects. This method checks the intersection only in 2d.

Parameters
abfirst line segment
cdsecond line segment
Returns
empty vector in case there isn't an intersection point, a vector containing one point if the line segments intersect in a single point, a vector containing two points describing the line segment the original line segments are interfering.

Definition at line 457 of file AnalyticalGeometry.cpp.

459{
460 GeoLib::Point const& a{ab.getBeginPoint()};
461 GeoLib::Point const& b{ab.getEndPoint()};
462 GeoLib::Point const& c{cd.getBeginPoint()};
463 GeoLib::Point const& d{cd.getEndPoint()};
464
465 double const orient_abc(getOrientation(a, b, c));
466 double const orient_abd(getOrientation(a, b, d));
467
468 // check if the segment (cd) lies on the left or on the right of (ab)
469 if ((orient_abc > 0 && orient_abd > 0) ||
470 (orient_abc < 0 && orient_abd < 0))
471 {
472 return std::vector<MathLib::Point3d>();
473 }
474
475 // check: (cd) and (ab) are on the same line
476 if (orient_abc == 0.0 && orient_abd == 0.0)
477 {
478 double const eps(std::numeric_limits<double>::epsilon());
479 if (MathLib::sqrDist2d(a, c) < eps && MathLib::sqrDist2d(b, d) < eps)
480 {
481 return {{a, b}};
482 }
483 if (MathLib::sqrDist2d(a, d) < eps && MathLib::sqrDist2d(b, c) < eps)
484 {
485 return {{a, b}};
486 }
487
488 // Since orient_ab and orient_abd vanish, a, b, c, d are on the same
489 // line and for this reason it is enough to check the x-component.
490 auto isPointOnSegment = [](double q, double p0, double p1)
491 {
492 double const t((q - p0) / (p1 - p0));
493 return 0 <= t && t <= 1;
494 };
495
496 // check if c in (ab)
497 if (isPointOnSegment(c[0], a[0], b[0]))
498 {
499 // check if a in (cd)
500 if (isPointOnSegment(a[0], c[0], d[0]))
501 {
502 return {{a, c}};
503 }
504 // check b == c
505 if (MathLib::sqrDist2d(b, c) < eps)
506 {
507 return {{b}};
508 }
509 // check if b in (cd)
510 if (isPointOnSegment(b[0], c[0], d[0]))
511 {
512 return {{b, c}};
513 }
514 // check d in (ab)
515 if (isPointOnSegment(d[0], a[0], b[0]))
516 {
517 return {{c, d}};
518 }
519 std::stringstream err;
520 err.precision(std::numeric_limits<double>::max_digits10);
521 err << ab << " x " << cd;
522 OGS_FATAL(
523 "The case of parallel line segments ({:s}) is not handled yet. "
524 "Aborting.",
525 err.str());
526 }
527
528 // check if d in (ab)
529 if (isPointOnSegment(d[0], a[0], b[0]))
530 {
531 // check if a in (cd)
532 if (isPointOnSegment(a[0], c[0], d[0]))
533 {
534 return {{a, d}};
535 }
536 // check if b==d
537 if (MathLib::sqrDist2d(b, d) < eps)
538 {
539 return {{b}};
540 }
541 // check if b in (cd)
542 if (isPointOnSegment(b[0], c[0], d[0]))
543 {
544 return {{b, d}};
545 }
546 // d in (ab), b not in (cd): check c in (ab)
547 if (isPointOnSegment(c[0], a[0], b[0]))
548 {
549 return {{c, d}};
550 }
551
552 std::stringstream err;
553 err.precision(std::numeric_limits<double>::max_digits10);
554 err << ab << " x " << cd;
555 OGS_FATAL(
556 "The case of parallel line segments ({:s}) is not handled yet. "
557 "Aborting.",
558 err.str());
559 }
560 return std::vector<MathLib::Point3d>();
561 }
562
563 // precondition: points a, b, c are collinear
564 // the function checks if the point c is onto the line segment (a,b)
565 auto isCollinearPointOntoLineSegment = [](MathLib::Point3d const& a,
566 MathLib::Point3d const& b,
567 MathLib::Point3d const& c)
568 {
569 if (b[0] - a[0] != 0)
570 {
571 double const t = (c[0] - a[0]) / (b[0] - a[0]);
572 return 0.0 <= t && t <= 1.0;
573 }
574 if (b[1] - a[1] != 0)
575 {
576 double const t = (c[1] - a[1]) / (b[1] - a[1]);
577 return 0.0 <= t && t <= 1.0;
578 }
579 if (b[2] - a[2] != 0)
580 {
581 double const t = (c[2] - a[2]) / (b[2] - a[2]);
582 return 0.0 <= t && t <= 1.0;
583 }
584 return false;
585 };
586
587 if (orient_abc == 0.0)
588 {
589 if (isCollinearPointOntoLineSegment(a, b, c))
590 {
591 return {{c}};
592 }
593 return std::vector<MathLib::Point3d>();
594 }
595
596 if (orient_abd == 0.0)
597 {
598 if (isCollinearPointOntoLineSegment(a, b, d))
599 {
600 return {{d}};
601 }
602 return std::vector<MathLib::Point3d>();
603 }
604
605 // check if the segment (ab) lies on the left or on the right of (cd)
606 double const orient_cda(getOrientation(c, d, a));
607 double const orient_cdb(getOrientation(c, d, b));
608 if ((orient_cda > 0 && orient_cdb > 0) ||
609 (orient_cda < 0 && orient_cdb < 0))
610 {
611 return std::vector<MathLib::Point3d>();
612 }
613
614 // at this point it is sure that there is an intersection and the system of
615 // linear equations will be invertible
616 // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously
617 Eigen::Matrix2d mat;
618 mat(0, 0) = b[0] - a[0];
619 mat(0, 1) = c[0] - d[0];
620 mat(1, 0) = b[1] - a[1];
621 mat(1, 1) = c[1] - d[1];
622 Eigen::Vector2d rhs{c[0] - a[0], c[1] - a[1]};
623
624 rhs = mat.partialPivLu().solve(rhs);
625 if (0 <= rhs[1] && rhs[1] <= 1.0)
626 {
627 return {MathLib::Point3d{std::array<double, 3>{
628 {c[0] + rhs[1] * (d[0] - c[0]), c[1] + rhs[1] * (d[1] - c[1]),
629 c[2] + rhs[1] * (d[2] - c[2])}}}};
630 }
631 return std::vector<MathLib::Point3d>(); // parameter s not in the valid
632 // range
633}
Orientation getOrientation(MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
double sqrDist2d(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.h:123

References GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), getOrientation(), OGS_FATAL, and MathLib::sqrDist2d().

Referenced by MeshGeoToolsLib::computeElementSegmentIntersections().

◆ lineSegmentsIntersect()

bool GeoLib::lineSegmentsIntersect ( const Polyline * ply,
Polyline::SegmentIterator & seg_it0,
Polyline::SegmentIterator & seg_it1,
Point & intersection_pnt )

test for intersections of the line segments of the Polyline

Parameters
plythe polyline
seg_it0iterator pointing to the first segment that has an intersection
seg_it1iterator pointing to the second segment that has an intersection
intersection_pntthe intersection point if the segments intersect
Returns
true, if the polyline contains intersections

Definition at line 263 of file AnalyticalGeometry.cpp.

267{
268 std::size_t const n_segs(ply->getNumberOfSegments());
269 // Neighbouring segments always intersects at a common vertex. The algorithm
270 // checks for intersections of non-neighbouring segments.
271 for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0)
272 {
273 seg_it1 = seg_it0 + 2;
274 std::size_t const seg_num_0 = seg_it0.getSegmentNumber();
275 for (; seg_it1 != ply->end(); ++seg_it1)
276 {
277 // Do not check first and last segment, because they are
278 // neighboured.
279 if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1))
280 {
281 if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt))
282 {
283 return true;
284 }
285 }
286 }
287 }
288 return false;
289}
std::size_t getSegmentNumber() const
Definition Polyline.cpp:387
std::size_t getNumberOfSegments() const
Definition Polyline.cpp:114

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), GeoLib::Polyline::getNumberOfSegments(), GeoLib::Polyline::SegmentIterator::getSegmentNumber(), and lineSegmentIntersect().

Referenced by GeoLib::Polygon::splitPolygonAtIntersection().

◆ markUnusedPoints()

void GeoLib::markUnusedPoints ( GEOObjects const & geo_objects,
std::string const & geo_name,
std::vector< bool > & transfer_pnts )

Definition at line 573 of file GEOObjects.cpp.

576{
577 GeoLib::PolylineVec const* const ply_obj(
578 geo_objects.getPolylineVecObj(geo_name));
579 if (ply_obj)
580 {
581 std::vector<GeoLib::Polyline*> const& lines(ply_obj->getVector());
582 for (auto* line : lines)
583 {
584 std::size_t const n_pnts(line->getNumberOfPoints());
585 for (std::size_t i = 0; i < n_pnts; ++i)
586 {
587 transfer_pnts[line->getPointID(i)] = false;
588 }
589 }
590 }
591
592 GeoLib::SurfaceVec const* const sfc_obj(
593 geo_objects.getSurfaceVecObj(geo_name));
594 if (sfc_obj)
595 {
596 std::vector<GeoLib::Surface*> const& surfaces = sfc_obj->getVector();
597 for (auto* sfc : surfaces)
598 {
599 std::size_t const n_tri(sfc->getNumberOfTriangles());
600 for (std::size_t i = 0; i < n_tri; ++i)
601 {
602 GeoLib::Triangle const& t = *(*sfc)[i];
603 transfer_pnts[t[0]] = false;
604 transfer_pnts[t[1]] = false;
605 transfer_pnts[t[2]] = false;
606 }
607 }
608 }
609}
The class TemplateVec takes a unique name and manages a std::vector of pointers to data elements of t...
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
Definition Triangle.h:27

References GeoLib::GEOObjects::getPolylineVecObj(), GeoLib::GEOObjects::getSurfaceVecObj(), and GeoLib::TemplateVec< T >::getVector().

Referenced by geoPointsToStations().

◆ markUsedPoints() [1/2]

void GeoLib::markUsedPoints ( Polyline const & polyline,
std::vector< bool > & used_points )

Resets the point IDs of the polyline corresponding to the mapping.

Definition at line 490 of file Polyline.cpp.

491{
492 if (polyline.getPointsVec().size() != used_points.size())
493 {
494 OGS_FATAL(
495 "internal error in markUsedPoints(): polyline based on point "
496 "vector of size {}, given used_points has size {}",
497 polyline.getPointsVec().size(), used_points.size());
498 }
499 for (std::size_t i = 0; i < polyline.getNumberOfPoints(); ++i)
500 {
501 used_points[polyline.getPointID(i)] = true;
502 }
503}

References GeoLib::Polyline::getNumberOfPoints(), GeoLib::Polyline::getPointID(), GeoLib::Polyline::getPointsVec(), and OGS_FATAL.

Referenced by main().

◆ markUsedPoints() [2/2]

void GeoLib::markUsedPoints ( Surface const & surface,
std::vector< bool > & used_points )

Definition at line 132 of file Surface.cpp.

133{
134 if (surface.getPointVec()->size() != used_points.size())
135 {
136 OGS_FATAL(
137 "internal error in markUsedPoints(): surface based on point vector "
138 "of size {}, given used_points has size {}",
139 surface.getPointVec()->size(), used_points.size());
140 }
141 for (std::size_t i = 0; i < surface.getNumberOfTriangles(); ++i)
142 {
143 auto const& triangle = *surface[i];
144 for (std::size_t k = 0; k < 3; ++k)
145 {
146 used_points[triangle[k]] = true;
147 }
148 }
149}

References GeoLib::Surface::getNumberOfTriangles(), GeoLib::Surface::getPointVec(), and OGS_FATAL.

◆ operator<<() [1/2]

std::ostream & GeoLib::operator<< ( std::ostream & os,
LineSegment const & s )

Definition at line 86 of file LineSegment.cpp.

87{
88 os << "{(" << s.getBeginPoint() << "), (" << s.getEndPoint() << ")}";
89 return os;
90}

References GeoLib::LineSegment::getBeginPoint(), and GeoLib::LineSegment::getEndPoint().

◆ operator<<() [2/2]

std::ostream & GeoLib::operator<< ( std::ostream & os,
std::pair< GeoLib::LineSegment const &, GeoLib::LineSegment const & > const & seg_pair )

Definition at line 92 of file LineSegment.cpp.

95{
96 os << seg_pair.first << " x " << seg_pair.second;
97 return os;
98}

◆ operator==() [1/4]

bool GeoLib::operator== ( LineSegment const & s0,
LineSegment const & s1 )

Compares two objects of type LineSegment.

Note
The comparison is not strict. Two line segments \(s_0 = (b_0, e_0)\) and \(s_1 = (b_1, e_1)\) are equal, if \(\|b_0-b_1\|^2 < \epsilon\) and \(\|e_0-e_1\|^2 < \epsilon\) or \(\|b_0-e_1\|^2 < \epsilon\) and \(\|e_0-b_1\|^2 < \epsilon\) where \(\epsilon\) is the machine precision of double.

Definition at line 100 of file LineSegment.cpp.

101{
102 double const tol(std::numeric_limits<double>::epsilon());
103 return (MathLib::sqrDist(s0.getBeginPoint(), s1.getBeginPoint()) < tol &&
104 MathLib::sqrDist(s0.getEndPoint(), s1.getEndPoint()) < tol) ||
105 (MathLib::sqrDist(s0.getBeginPoint(), s1.getEndPoint()) < tol &&
106 MathLib::sqrDist(s0.getEndPoint(), s1.getBeginPoint()) < tol);
107}

References GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), and MathLib::sqrDist().

◆ operator==() [2/4]

bool GeoLib::operator== ( Polygon const & lhs,
Polygon const & rhs )

comparison operator for polygons

Parameters
lhsthe first polygon
rhsthe second polygon
Returns
true, if the polygons describe the same geometrical object

Definition at line 596 of file Polygon.cpp.

597{
598 if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
599 {
600 return false;
601 }
602
603 const std::size_t n(lhs.getNumberOfPoints());
604 const std::size_t start_pnt(lhs.getPointID(0));
605
606 // search start point of first polygon in second polygon
607 bool nfound(true);
608 std::size_t k(0);
609 for (; k < n - 1 && nfound; k++)
610 {
611 if (start_pnt == rhs.getPointID(k))
612 {
613 nfound = false;
614 break;
615 }
616 }
617
618 // case: start point not found in second polygon
619 if (nfound)
620 {
621 return false;
622 }
623
624 // *** determine direction
625 // opposite direction
626 if (k == n - 2)
627 {
628 for (k = 1; k < n - 1; k++)
629 {
630 if (lhs.getPointID(k) != rhs.getPointID(n - 1 - k))
631 {
632 return false;
633 }
634 }
635 return true;
636 }
637
638 // same direction - start point of first polygon at arbitrary position in
639 // second polygon
640 if (lhs.getPointID(1) == rhs.getPointID(k + 1))
641 {
642 std::size_t j(k + 2);
643 for (; j < n - 1; j++)
644 {
645 if (lhs.getPointID(j - k) != rhs.getPointID(j))
646 {
647 return false;
648 }
649 }
650 j = 0; // new start point at second polygon
651 for (; j < k + 1; j++)
652 {
653 if (lhs.getPointID(n - (k + 2) + j + 1) != rhs.getPointID(j))
654 {
655 return false;
656 }
657 }
658 return true;
659 }
660 // opposite direction with start point of first polygon at arbitrary
661 // position
662 // *** ATTENTION
663 WARN(
664 "operator==(Polygon const& lhs, Polygon const& rhs) - not tested case "
665 "(implementation is probably buggy) - please contact "
666 "thomas.fischer@ufz.de mentioning the problem.");
667 // in second polygon
668 if (lhs.getPointID(1) == rhs.getPointID(k - 1))
669 {
670 std::size_t j(k - 2);
671 for (; j > 0; j--)
672 {
673 if (lhs.getPointID(k - 2 - j) != rhs.getPointID(j))
674 {
675 return false;
676 }
677 }
678 // new start point at second polygon - the point n-1 of a polygon is
679 // equal to the first point of the polygon (for this reason: n-2)
680 j = n - 2;
681 for (; j > k - 1; j--)
682 {
683 if (lhs.getPointID(n - 2 + j + k - 2) != rhs.getPointID(j))
684 {
685 return false;
686 }
687 }
688 return true;
689 }
690 return false;
691}

◆ operator==() [3/4]

bool GeoLib::operator== ( Polyline const & lhs,
Polyline const & rhs )

comparison operator

Parameters
lhsfirst polyline
rhssecond polyline
Returns
true, if the polylines consists of the same sequence of line segments

Definition at line 533 of file Polyline.cpp.

534{
535 if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
536 {
537 return false;
538 }
539
540 const std::size_t n(lhs.getNumberOfPoints());
541 for (std::size_t k(0); k < n; k++)
542 {
543 if (lhs.getPointID(k) != rhs.getPointID(k))
544 {
545 return false;
546 }
547 }
548
549 return true;
550}

References GeoLib::Polyline::getNumberOfPoints(), and GeoLib::Polyline::getPointID().

◆ operator==() [4/4]

bool GeoLib::operator== ( Surface const & lhs,
Surface const & rhs )

Definition at line 109 of file Surface.cpp.

110{
111 return &lhs == &rhs;
112}

◆ parallel()

bool GeoLib::parallel ( Eigen::Vector3d v,
Eigen::Vector3d w )

Check if the two vectors \(v, w \in R^3\) are in parallel

Parameters
vfirst vector
wsecond vector
Returns
true if the vectors are in parallel, else false

Definition at line 83 of file AnalyticalGeometry.cpp.

84{
85 const double eps(std::numeric_limits<double>::epsilon());
86 double const eps_squared = eps * eps;
87
88 // check degenerated cases
89 if (v.squaredNorm() < eps_squared)
90 {
91 return false;
92 }
93
94 if (w.squaredNorm() < eps_squared)
95 {
96 return false;
97 }
98
99 v.normalize();
100 w.normalize();
101
102 bool parallel(true);
103 if (std::abs(v[0] - w[0]) > eps)
104 {
105 parallel = false;
106 }
107 if (std::abs(v[1] - w[1]) > eps)
108 {
109 parallel = false;
110 }
111 if (std::abs(v[2] - w[2]) > eps)
112 {
113 parallel = false;
114 }
115
116 if (!parallel)
117 {
118 parallel = true;
119 // change sense of direction of v_normalised
120 v *= -1.0;
121 // check again
122 if (std::abs(v[0] - w[0]) > eps)
123 {
124 parallel = false;
125 }
126 if (std::abs(v[1] - w[1]) > eps)
127 {
128 parallel = false;
129 }
130 if (std::abs(v[2] - w[2]) > eps)
131 {
132 parallel = false;
133 }
134 }
135
136 return parallel;
137}

References parallel().

Referenced by lineSegmentIntersect(), and parallel().

◆ pointsAreIdentical()

bool GeoLib::pointsAreIdentical ( const std::vector< Point * > & pnt_vec,
std::size_t i,
std::size_t j,
double prox )

Definition at line 552 of file Polyline.cpp.

556{
557 if (i == j)
558 {
559 return true;
560 }
561 return MathLib::sqrDist(*pnt_vec[i], *pnt_vec[j]) < prox;
562}

References MathLib::sqrDist().

Referenced by GeoLib::Polyline::constructPolylineFromSegments().

◆ resetPointIDs() [1/2]

void GeoLib::resetPointIDs ( Polyline & polyline,
std::vector< std::size_t > const & mapping )

Resets the point IDs of the polyline corresponding to the mapping.

Definition at line 475 of file Polyline.cpp.

476{
477 if (polyline.getPointsVec().size() != mapping.size())
478 {
479 OGS_FATAL(
480 "internal error in resetPointIDs(): polyline based on point vector "
481 "of size {}, given mapping has size {}",
482 polyline.getPointsVec().size(), mapping.size());
483 }
484 for (std::size_t i = 0; i < polyline.getNumberOfPoints(); ++i)
485 {
486 polyline.setPointID(i, mapping[polyline.getPointID(i)]);
487 }
488}
void setPointID(std::size_t idx, std::size_t id)
Definition Polyline.cpp:173
std::vector< Point * > const & getPointsVec() const
Definition Polyline.cpp:185

References GeoLib::Polyline::getNumberOfPoints(), GeoLib::Polyline::getPointID(), GeoLib::Polyline::getPointsVec(), OGS_FATAL, and GeoLib::Polyline::setPointID().

Referenced by main(), and FileIO::GMSH::GMSHPolygonTree::writeAdditionalPointData().

◆ resetPointIDs() [2/2]

void GeoLib::resetPointIDs ( Surface & surface,
std::vector< std::size_t > const & mapping )

Resets the point IDs of the surface corresponding to the mapping.

Definition at line 114 of file Surface.cpp.

115{
116 if (surface.getPointVec()->size() != mapping.size())
117 {
118 OGS_FATAL(
119 "internal error in resetPointIDs(): surface based on point vector "
120 "of size {}, given mapping vector has size {}",
121 surface.getPointVec()->size(), mapping.size());
122 }
123 for (std::size_t i = 0; i < surface.getNumberOfTriangles(); ++i)
124 {
125 auto& triangle = *surface[i];
126 const_cast<std::size_t&>(triangle[0]) = mapping[triangle[0]];
127 const_cast<std::size_t&>(triangle[1]) = mapping[triangle[1]];
128 const_cast<std::size_t&>(triangle[2]) = mapping[triangle[2]];
129 }
130}
std::size_t getNumberOfTriangles() const
Definition Surface.cpp:82
const std::vector< Point * > * getPointVec() const
Definition Surface.h:73

References GeoLib::Surface::getNumberOfTriangles(), GeoLib::Surface::getPointVec(), and OGS_FATAL.

◆ rotatePoints() [1/3]

template<typename InputIterator >
void GeoLib::rotatePoints ( Eigen::Matrix3d const & rot_mat,
InputIterator pnts_begin,
InputIterator pnts_end )

rotate points according to the rotation matrix

Parameters
rot_mat3x3 dimensional rotation matrix
pnts_beginIterator pointing to the initial element in a vector of points to be rotated
pnts_endIterator pointing to the element following the last element in a vector of points to be rotated

Definition at line 77 of file AnalyticalGeometry-impl.h.

79{
80 for (auto it = pnts_begin; it != pnts_end; ++it)
81 {
82 (*it)->asEigenVector3d() = rot_mat * (*it)->asEigenVector3d();
83 }
84}

Referenced by MeshGeoToolsLib::markNodesOutSideOfPolygon(), rotateGeometryToXY(), rotateMesh(), rotatePoints(), rotatePoints(), rotatePointsToXY(), rotatePolygonPointsToXY(), and FileIO::GMSH::GMSHInterface::writeGMSHInputFile().

◆ rotatePoints() [2/3]

void GeoLib::rotatePoints ( Eigen::Matrix3d const & rot_mat,
std::vector< GeoLib::Point * > & pnts )

Definition at line 291 of file AnalyticalGeometry.cpp.

293{
294 rotatePoints(rot_mat, pnts.begin(), pnts.end());
295}
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)

References rotatePoints().

◆ rotatePoints() [3/3]

template<typename P >
void GeoLib::rotatePoints ( Eigen::Matrix3d const & rot_mat,
std::vector< P * > const & pnts )

rotate points according to the rotation matrix

Parameters
rot_mat3x3 dimensional rotation matrix
pntsvector of points

Definition at line 111 of file AnalyticalGeometry-impl.h.

112{
113 rotatePoints(rot_mat, pnts.begin(), pnts.end());
114}

References rotatePoints().

◆ rotatePointsToXY() [1/3]

template<typename InputIterator1 , typename InputIterator2 >
Eigen::Matrix3d GeoLib::rotatePointsToXY ( InputIterator1 p_pnts_begin,
InputIterator1 p_pnts_end,
InputIterator2 r_pnts_begin,
InputIterator2 r_pnts_end )

rotate points to X-Y plane

Parameters
p_pnts_beginIterator pointing to the initial element in a vector of points used for computing a rotation matrix
p_pnts_endIterator pointing to the element following the last point in a vector of points used for computing a rotation matrix
r_pnts_beginIterator pointing to the initial element in a vector of points to be rotated
r_pnts_endIterator pointing to the element following the last point in a vector of points to be rotated Points are rotated using a rotation matrix computed from the first three points in the vector. Point coordinates are modified as a result of the rotation.

Definition at line 87 of file AnalyticalGeometry-impl.h.

91{
92 assert(std::distance(p_pnts_begin, p_pnts_end) > 2);
93
94 // compute the plane normal
95 auto const [plane_normal, d] =
96 GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
97
98 // rotate points into x-y-plane
99 Eigen::Matrix3d const rot_mat = computeRotationMatrixToXY(plane_normal);
100 rotatePoints(rot_mat, r_pnts_begin, r_pnts_end);
101
102 for (auto it = r_pnts_begin; it != r_pnts_end; ++it)
103 {
104 (*(*it))[2] = 0.0; // should be -= d but there are numerical errors
105 }
106
107 return rot_mat;
108}
Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const &n)

References computeRotationMatrixToXY(), getNewellPlane(), and rotatePoints().

Referenced by GeoLib::EarClippingTriangulation::EarClippingTriangulation(), GeoLib::Polygon::ensureCCWOrientation(), rotatePointsToXY(), and FileIO::GMSH::GMSHInterface::writeGMSHInputFile().

◆ rotatePointsToXY() [2/3]

Eigen::Matrix3d GeoLib::rotatePointsToXY ( std::vector< Point * > & pnts)

rotate points to X-Y plane

Parameters
pntsa vector of points with a minimum length of three. Points are rotated using a rotation matrix computed from the first three points in the vector. Point coordinates are modified as a result of the rotation.

Definition at line 360 of file AnalyticalGeometry.cpp.

361{
362 return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
363}
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)

References rotatePointsToXY().

◆ rotatePointsToXY() [3/3]

Eigen::Matrix3d GeoLib::rotatePointsToXY ( std::vector< Point * > & pnts)

rotate points to X-Y plane

Parameters
pntsa vector of points with a minimum length of three. Points are rotated using a rotation matrix computed from the first three points in the vector. Point coordinates are modified as a result of the rotation.

Definition at line 360 of file AnalyticalGeometry.cpp.

361{
362 return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
363}

References rotatePointsToXY().

◆ rotatePolygonPointsToXY()

std::tuple< std::vector< GeoLib::Point * >, Eigen::Vector3d > GeoLib::rotatePolygonPointsToXY ( GeoLib::Polygon const & polygon_in)

Function rotates a polygon to the xy plane. For this reason, (1) the points of the given polygon are copied, (2) a so called Newell plane is computed (getNewellPlane()) and the points are rotated, (3) for accuracy reasons the \(z\) coordinates of the rotated points are set to zero

See also
getNewellPlane()
Parameters
polygon_ina copy of the polygon_in polygon will be rotated
Returns
vector of rotated points and normal based on the original Newell plane

Definition at line 434 of file AnalyticalGeometry.cpp.

435{
436 // 1 copy all points
437 std::vector<GeoLib::Point*> polygon_points;
438 polygon_points.reserve(polygon_in.getNumberOfPoints());
439 for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
440 {
441 polygon_points.push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
442 }
443
444 // 2 rotate points
445 auto [plane_normal, d_polygon] = GeoLib::getNewellPlane(polygon_points);
446 Eigen::Matrix3d const rot_mat =
448 GeoLib::rotatePoints(rot_mat, polygon_points);
449
450 // 3 set z coord to zero
451 std::for_each(polygon_points.begin(), polygon_points.end(),
452 [](GeoLib::Point* p) { (*p)[2] = 0.0; });
453
454 return {polygon_points, plane_normal};
455}

References computeRotationMatrixToXY(), getNewellPlane(), GeoLib::Polyline::getNumberOfPoints(), GeoLib::Polyline::getPoint(), and rotatePoints().

Referenced by MeshGeoToolsLib::markNodesOutSideOfPolygon().

◆ sortSegments() [1/2]

void GeoLib::sortSegments ( MathLib::Point3d const & seg_beg_pnt,
std::vector< LineSegment > & sub_segments )

Sorts the vector of segments such that the \(i\)-th segment is connected with the \(i+1\)st segment, i.e. the end point of the \(i\)-th segment is the start point of the \(i+1\)st segment. The current implementation requires that all segments have to be connectable. In order to obtain a unique result the segments are sorted such that the begin point of the first segment is seg_beg_pnt.

Definition at line 635 of file AnalyticalGeometry.cpp.

637{
638 double const eps(std::numeric_limits<double>::epsilon());
639
640 auto findNextSegment =
641 [&eps](MathLib::Point3d const& seg_beg_pnt,
642 std::vector<GeoLib::LineSegment>& sub_segments,
643 std::vector<GeoLib::LineSegment>::iterator& sub_seg_it)
644 {
645 if (sub_seg_it == sub_segments.end())
646 {
647 return;
648 }
649 // find appropriate segment for the given segment begin point
650 auto act_beg_seg_it = std::find_if(
651 sub_seg_it, sub_segments.end(),
652 [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
653 {
654 return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) <
655 eps ||
656 MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
657 });
658 if (act_beg_seg_it == sub_segments.end())
659 {
660 return;
661 }
662 // if necessary correct orientation of segment, i.e. swap beg and
663 // end
664 if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
665 MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
666 {
667 std::swap(act_beg_seg_it->getBeginPoint(),
668 act_beg_seg_it->getEndPoint());
669 }
670 assert(sub_seg_it != sub_segments.end());
671 // exchange segments within the container
672 if (sub_seg_it != act_beg_seg_it)
673 {
674 std::swap(*sub_seg_it, *act_beg_seg_it);
675 }
676 };
677
678 // find start segment
679 auto seg_it = sub_segments.begin();
680 findNextSegment(seg_beg_pnt, sub_segments, seg_it);
681
682 while (seg_it != sub_segments.end())
683 {
684 MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint());
685 seg_it++;
686 if (seg_it != sub_segments.end())
687 {
688 findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
689 }
690 }
691}

References MathLib::sqrDist().

Referenced by MeshGeoToolsLib::mapLineSegment().

◆ sortSegments() [2/2]

void GeoLib::sortSegments ( MathLib::Point3d const & seg_beg_pnt,
std::vector< LineSegment > & sub_segments )

Sorts the vector of segments such that the \(i\)-th segment is connected with the \(i+1\)st segment, i.e. the end point of the \(i\)-th segment is the start point of the \(i+1\)st segment. The current implementation requires that all segments have to be connectable. In order to obtain a unique result the segments are sorted such that the begin point of the first segment is seg_beg_pnt.

Definition at line 635 of file AnalyticalGeometry.cpp.

637{
638 double const eps(std::numeric_limits<double>::epsilon());
639
640 auto findNextSegment =
641 [&eps](MathLib::Point3d const& seg_beg_pnt,
642 std::vector<GeoLib::LineSegment>& sub_segments,
643 std::vector<GeoLib::LineSegment>::iterator& sub_seg_it)
644 {
645 if (sub_seg_it == sub_segments.end())
646 {
647 return;
648 }
649 // find appropriate segment for the given segment begin point
650 auto act_beg_seg_it = std::find_if(
651 sub_seg_it, sub_segments.end(),
652 [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
653 {
654 return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) <
655 eps ||
656 MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
657 });
658 if (act_beg_seg_it == sub_segments.end())
659 {
660 return;
661 }
662 // if necessary correct orientation of segment, i.e. swap beg and
663 // end
664 if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
665 MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
666 {
667 std::swap(act_beg_seg_it->getBeginPoint(),
668 act_beg_seg_it->getEndPoint());
669 }
670 assert(sub_seg_it != sub_segments.end());
671 // exchange segments within the container
672 if (sub_seg_it != act_beg_seg_it)
673 {
674 std::swap(*sub_seg_it, *act_beg_seg_it);
675 }
676 };
677
678 // find start segment
679 auto seg_it = sub_segments.begin();
680 findNextSegment(seg_beg_pnt, sub_segments, seg_it);
681
682 while (seg_it != sub_segments.end())
683 {
684 MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint());
685 seg_it++;
686 if (seg_it != sub_segments.end())
687 {
688 findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
689 }
690 }
691}

References MathLib::sqrDist().

Referenced by MeshGeoToolsLib::mapLineSegment().

◆ triangleLineIntersection()

std::unique_ptr< Point > GeoLib::triangleLineIntersection ( MathLib::Point3d const & a,
MathLib::Point3d const & b,
MathLib::Point3d const & c,
MathLib::Point3d const & p,
MathLib::Point3d const & q )

Calculates the intersection points of a line PQ and a triangle ABC. This method requires ABC to be counterclockwise and PQ to point downward.

Returns
Intersection point or nullptr if there is no intersection.

Definition at line 365 of file AnalyticalGeometry.cpp.

369{
370 Eigen::Vector3d const pq = q.asEigenVector3d() - p.asEigenVector3d();
371 Eigen::Vector3d const pa = a.asEigenVector3d() - p.asEigenVector3d();
372 Eigen::Vector3d const pb = b.asEigenVector3d() - p.asEigenVector3d();
373 Eigen::Vector3d const pc = c.asEigenVector3d() - p.asEigenVector3d();
374
375 double u = pq.cross(pc).dot(pb);
376 if (u < 0)
377 {
378 return nullptr;
379 }
380 double v = pq.cross(pa).dot(pc);
381 if (v < 0)
382 {
383 return nullptr;
384 }
385 double w = pq.cross(pb).dot(pa);
386 if (w < 0)
387 {
388 return nullptr;
389 }
390
391 const double denom(1.0 / (u + v + w));
392 u *= denom;
393 v *= denom;
394 w *= denom;
395 return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
396 u * a[1] + v * b[1] + w * c[1],
397 u * a[2] + v * b[2] + w * c[2]);
398}

References MathLib::Point3d::asEigenVector3d().

Referenced by MeshGeoToolsLib::GeoMapper::getMeshElevation().