OGS
GeoLib Namespace Reference

Detailed Description

Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) Distributed under a Modified BSD License. See accompanying file LICENSE.txt or http://www.opengeosys.org/project/license

Contains all functionality concerned with defining and obtaining information about geometric objects such as points, (poly-)lines, polygones, surfaces, etc. Also includes additional algorithms for modifying such data.

Namespaces

 IO
 

Classes

class  QuadTree
 
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
 
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
 
struct  RasterHeader
 Contains the relevant information when storing a geoscientific raster data. More...
 
class  Raster
 Class Raster is used for managing 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 More...
 
using SurfaceVec = TemplateVec< GeoLib::Surface >
 

Enumerations

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

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)
 
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. More...
 
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)
 
bool operator== (Polygon const &lhs, Polygon const &rhs)
 
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)
 
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)
 

Typedef Documentation

◆ SurfaceVec

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

Definition at line 27 of file SurfaceVec.h.

Enumeration Type Documentation

◆ Location

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

Definition at line 30 of file Polyline.h.

◆ 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 705 of file AnalyticalGeometry.cpp.

706 {
707  Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
708  const double cos_theta = v[0];
709  const double sin_theta = v[1];
710  rot_mat(0, 0) = rot_mat(1, 1) = cos_theta;
711  rot_mat(0, 1) = sin_theta;
712  rot_mat(1, 0) = -sin_theta;
713  rot_mat(2, 2) = 1.0;
714  return rot_mat;
715 }
static const double v

References MathLib::v.

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 717 of file AnalyticalGeometry.cpp.

718 {
719  // a vector on the plane
720  Eigen::Vector3d yy = Eigen::Vector3d::Zero();
721  auto const eps = std::numeric_limits<double>::epsilon();
722  if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
723  {
724  yy[2] = 1.0;
725  }
726  else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
727  {
728  yy[0] = 1.0;
729  }
730  else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
731  {
732  yy[1] = 1.0;
733  }
734  else
735  {
736  for (unsigned i = 0; i < 3; i++)
737  {
738  if (std::abs(v[i]) > 0.0)
739  {
740  yy[i] = -v[i];
741  break;
742  }
743  }
744  }
745  // z"_vec
746  Eigen::Vector3d const zz = v.cross(yy).normalized();
747  // y"_vec
748  yy = zz.cross(v).normalized();
749 
750  Eigen::Matrix3d rot_mat;
751  rot_mat.row(0) = v;
752  rot_mat.row(1) = yy;
753  rot_mat.row(2) = zz;
754  return rot_mat;
755 }

References MathLib::v.

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 412 of file AnalyticalGeometry.cpp.

414 {
415  auto computeSegmentIntersections =
416  [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
417  {
418  for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
419  {
420  for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
421  {
422  GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
423  if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
424  {
425  std::size_t const id(
426  pnt_vec.push_back(new GeoLib::Point(s)));
427  poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
428  poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
429  }
430  }
431  }
432  };
433 
434  for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
435  {
436  auto it1(it0);
437  ++it1;
438  for (; it1 != plys.end(); ++it1)
439  {
440  computeSegmentIntersections(*(*it0), *(*it1));
441  }
442  }
443 }
std::size_t push_back(Point *pnt)
Definition: PointVec.cpp:129
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:51
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition: Polyline.cpp:54
SegmentIterator begin() const
Definition: Polyline.h:189
SegmentIterator end() const
Definition: Polyline.h:194
std::size_t size() const
Definition: TemplateVec.h:106
bool lineSegmentIntersect(GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)
static const double s

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), GeoLib::Polyline::insertPoint(), lineSegmentIntersect(), GeoLib::PointVec::push_back(), MathLib::s, 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 412 of file AnalyticalGeometry.cpp.

414 {
415  auto computeSegmentIntersections =
416  [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
417  {
418  for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
419  {
420  for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
421  {
422  GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
423  if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
424  {
425  std::size_t const id(
426  pnt_vec.push_back(new GeoLib::Point(s)));
427  poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
428  poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
429  }
430  }
431  }
432  };
433 
434  for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
435  {
436  auto it1(it0);
437  ++it1;
438  for (; it1 != plys.end(); ++it1)
439  {
440  computeSegmentIntersections(*(*it0), *(*it1));
441  }
442  }
443 }

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), GeoLib::Polyline::insertPoint(), lineSegmentIntersect(), GeoLib::PointVec::push_back(), MathLib::s, 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 303 of file AnalyticalGeometry.cpp.

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

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 519 of file Polyline.cpp.

520 {
521  if (id0 == id1)
522  {
523  ERR("no valid edge id0 == id1 == {:d}.", id0);
524  return false;
525  }
526  if (id0 > id1)
527  {
528  std::swap(id0, id1);
529  }
530  const std::size_t n(ply.getNumberOfPoints() - 1);
531  for (std::size_t k(0); k < n; k++)
532  {
533  std::size_t ply_pnt0(ply.getPointID(k));
534  std::size_t ply_pnt1(ply.getPointID(k + 1));
535  if (ply_pnt0 > ply_pnt1)
536  {
537  std::swap(ply_pnt0, ply_pnt1);
538  }
539  if (ply_pnt0 == id0 && ply_pnt1 == id1)
540  {
541  return true;
542  }
543  }
544  return false;
545 }
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42

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 97 of file SimplePolygonTree.h.

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

◆ 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 681 of file GEOObjects.cpp.

683 {
684  GeoLib::PointVec const* const pnt_obj(geo_objects.getPointVecObj(geo_name));
685  if (pnt_obj == nullptr)
686  {
687  ERR("Point vector {:s} not found.", geo_name);
688  return -1;
689  }
690  std::vector<GeoLib::Point*> const& pnts = *pnt_obj->getVector();
691  if (pnts.empty())
692  {
693  ERR("Point vector {:s} is empty.", geo_name);
694  return -1;
695  }
696  std::size_t const n_pnts(pnts.size());
697  std::vector<bool> transfer_pnts(n_pnts, true);
698  if (only_unused_pnts)
699  {
700  markUnusedPoints(geo_objects, geo_name, transfer_pnts);
701  }
702 
703  auto stations = std::make_unique<std::vector<GeoLib::Point*>>();
704  for (std::size_t i = 0; i < n_pnts; ++i)
705  {
706  if (!transfer_pnts[i])
707  {
708  continue;
709  }
710  std::string name = pnt_obj->getItemNameByID(i);
711  if (name.empty())
712  {
713  name = "Station " + std::to_string(i);
714  }
715  stations->push_back(new GeoLib::Station(pnts[i], name));
716  }
717  if (!stations->empty())
718  {
719  geo_objects.addStationVec(std::move(stations), stn_name);
720  }
721  else
722  {
723  WARN("No points found to convert.");
724  return 1;
725  }
726  return 0;
727 }
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition: PointVec.h:38
A Station (observation site) is basically a Point with some additional information.
Definition: Station.h:37
void markUnusedPoints(GEOObjects const &geo_objects, std::string const &geo_name, std::vector< bool > &transfer_pnts)
Definition: GEOObjects.cpp:643

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

Referenced by MainWindow::convertPointsToStations().

◆ 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 [12] .

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 44 of file AnalyticalGeometry-impl.h.

46 {
47  return getNewellPlane(pnts.begin(), pnts.end());
48 }
std::pair< Eigen::Vector3d, double > getNewellPlane(const std::vector< T_POINT > &pnts)

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 51 of file AnalyticalGeometry-impl.h.

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

◆ 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 [12] .

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; i = j, ++j) {
20  auto &pt_i = *(*i);
21  auto &pt_j = *(*j);
22  plane_normal[0] += (pt_i[1] - pt_j[1])
23  * (pt_i[2] + pt_j[2]); // projection on yz
24  plane_normal[1] += (pt_i[2] - pt_j[2])
25  * (pt_i[0] + pt_j[0]); // projection on xz
26  plane_normal[2] += (pt_i[0] - pt_j[0])
27  * (pt_i[1] + pt_j[1]); // projection on xy
28 
29  centroid += Eigen::Map<Eigen::Vector3d const>(pt_j.getCoords());
30  }
31 
32  plane_normal.normalize();
33  auto n_pnts(std::distance(pnts_begin, pnts_end));
34  assert(n_pnts > 2);
35  double d = 0.0;
36  if (n_pnts > 0)
37  {
38  d = centroid.dot(plane_normal) / n_pnts;
39  }
40  return std::make_pair(plane_normal, d);
41 }

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 52 of file AnalyticalGeometry.cpp.

55 {
56  double const orientation = ExactPredicates::getOrientation2d(p0, p1, p2);
57  if (orientation > 0)
58  {
59  return CCW;
60  }
61  if (orientation < 0)
62  {
63  return CW;
64  }
65  return COLLINEAR;
66 }
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::clipEars(), GeoLib::Polygon::ensureCCWOrientation(), GeoLib::EarClippingTriangulation::ensureCWOrientation(), 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 68 of file AnalyticalGeometry.cpp.

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

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

Referenced by MeshLib::ProjectPointOnMesh::getProjectedElement().

◆ 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(), and VtkStationSource::RequestData().

◆ isStation()

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

Definition at line 76 of file Station.cpp.

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

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 141 of file AnalyticalGeometry.cpp.

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

References MaterialPropertyLib::c, GeoLib::LineSegment::getBeginPoint(), MathLib::TemplatePoint< T, DIM >::getCoords(), GeoLib::LineSegment::getEndPoint(), MathLib::isCoplanar(), parallel(), MathLib::s, MathLib::sqrDist(), MathLib::u, and MathLib::v.

Referenced by FileIO::GMSH::GMSHPolygonTree::checkIntersectionsSegmentExistingPolylines(), computeAndInsertAllIntersectionPoints(), GeoLib::Polygon::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 469 of file AnalyticalGeometry.cpp.

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

References MaterialPropertyLib::c, GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), getOrientation(), OGS_FATAL, MathLib::q, MathLib::sqrDist2d(), and MathLib::t.

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 269 of file AnalyticalGeometry.cpp.

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

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 643 of file GEOObjects.cpp.

646 {
647  GeoLib::PolylineVec const* const ply_obj(
648  geo_objects.getPolylineVecObj(geo_name));
649  if (ply_obj)
650  {
651  std::vector<GeoLib::Polyline*> const& lines(*ply_obj->getVector());
652  for (auto* line : lines)
653  {
654  std::size_t const n_pnts(line->getNumberOfPoints());
655  for (std::size_t i = 0; i < n_pnts; ++i)
656  {
657  transfer_pnts[line->getPointID(i)] = false;
658  }
659  }
660  }
661 
662  GeoLib::SurfaceVec const* const sfc_obj(
663  geo_objects.getSurfaceVecObj(geo_name));
664  if (sfc_obj)
665  {
666  std::vector<GeoLib::Surface*> const& surfaces = *sfc_obj->getVector();
667  for (auto* sfc : surfaces)
668  {
669  std::size_t const n_tri(sfc->getNumberOfTriangles());
670  for (std::size_t i = 0; i < n_tri; ++i)
671  {
672  GeoLib::Triangle const& t = *(*sfc)[i];
673  transfer_pnts[t[0]] = false;
674  transfer_pnts[t[1]] = false;
675  transfer_pnts[t[2]] = false;
676  }
677  }
678  }
679 }
The class TemplateVec takes a unique name and manages a std::vector of pointers to data elements of t...
Definition: TemplateVec.h:40
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
Definition: Triangle.h:26

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

Referenced by geoPointsToStations().

◆ 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 MathLib::s.

◆ 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 521 of file Polygon.cpp.

522 {
523  if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
524  {
525  return false;
526  }
527 
528  const std::size_t n(lhs.getNumberOfPoints());
529  const std::size_t start_pnt(lhs.getPointID(0));
530 
531  // search start point of first polygon in second polygon
532  bool nfound(true);
533  std::size_t k(0);
534  for (; k < n - 1 && nfound; k++)
535  {
536  if (start_pnt == rhs.getPointID(k))
537  {
538  nfound = false;
539  break;
540  }
541  }
542 
543  // case: start point not found in second polygon
544  if (nfound)
545  {
546  return false;
547  }
548 
549  // *** determine direction
550  // opposite direction
551  if (k == n - 2)
552  {
553  for (k = 1; k < n - 1; k++)
554  {
555  if (lhs.getPointID(k) != rhs.getPointID(n - 1 - k))
556  {
557  return false;
558  }
559  }
560  return true;
561  }
562 
563  // same direction - start point of first polygon at arbitrary position in
564  // second polygon
565  if (lhs.getPointID(1) == rhs.getPointID(k + 1))
566  {
567  std::size_t j(k + 2);
568  for (; j < n - 1; j++)
569  {
570  if (lhs.getPointID(j - k) != rhs.getPointID(j))
571  {
572  return false;
573  }
574  }
575  j = 0; // new start point at second polygon
576  for (; j < k + 1; j++)
577  {
578  if (lhs.getPointID(n - (k + 2) + j + 1) != rhs.getPointID(j))
579  {
580  return false;
581  }
582  }
583  return true;
584  }
585  // opposite direction with start point of first polygon at arbitrary
586  // position
587  // *** ATTENTION
588  WARN(
589  "operator==(Polygon const& lhs, Polygon const& rhs) - not tested case "
590  "(implementation is probably buggy) - please contact "
591  "thomas.fischer@ufz.de mentioning the problem.");
592  // in second polygon
593  if (lhs.getPointID(1) == rhs.getPointID(k - 1))
594  {
595  std::size_t j(k - 2);
596  for (; j > 0; j--)
597  {
598  if (lhs.getPointID(k - 2 - j) != rhs.getPointID(j))
599  {
600  return false;
601  }
602  }
603  // new start point at second polygon - the point n-1 of a polygon is
604  // equal to the first point of the polygon (for this reason: n-2)
605  j = n - 2;
606  for (; j > k - 1; j--)
607  {
608  if (lhs.getPointID(n - 2 + j + k - 2) != rhs.getPointID(j))
609  {
610  return false;
611  }
612  }
613  return true;
614  }
615  return false;
616 }

◆ 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 547 of file Polyline.cpp.

548 {
549  if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
550  {
551  return false;
552  }
553 
554  const std::size_t n(lhs.getNumberOfPoints());
555  for (std::size_t k(0); k < n; k++)
556  {
557  if (lhs.getPointID(k) != rhs.getPointID(k))
558  {
559  return false;
560  }
561  }
562 
563  return true;
564 }

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 85 of file AnalyticalGeometry.cpp.

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

References MathLib::v.

Referenced by lineSegmentIntersect().

◆ pointsAreIdentical()

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

Definition at line 566 of file Polyline.cpp.

570 {
571  if (i == j)
572  {
573  return true;
574  }
575  return MathLib::sqrDist(*pnt_vec[i], *pnt_vec[j]) < prox;
576 }

References MathLib::sqrDist().

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

◆ 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 75 of file AnalyticalGeometry-impl.h.

77 {
78  for (auto it = pnts_begin; it != pnts_end; ++it)
79  {
80  Eigen::Map<Eigen::Vector3d>((*it)->getCoords()) =
81  rot_mat * Eigen::Map<Eigen::Vector3d const>((*it)->getCoords());
82  }
83 }

Referenced by MeshGeoToolsLib::markNodesOutSideOfPolygon(), rotateGeometryToXY(), rotateMesh(), 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 297 of file AnalyticalGeometry.cpp.

299 {
300  rotatePoints(rot_mat, pnts.begin(), pnts.end());
301 }
void rotatePoints(Eigen::Matrix3d const &rot_mat, std::vector< GeoLib::Point * > &pnts)

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 110 of file AnalyticalGeometry-impl.h.

111 {
112  rotatePoints(rot_mat, pnts.begin(), pnts.end());
113 }
void rotatePoints(Eigen::Matrix3d const &rot_mat, std::vector< P * > const &pnts)

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 86 of file AnalyticalGeometry-impl.h.

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

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 366 of file AnalyticalGeometry.cpp.

367 {
368  return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
369 }
Eigen::Matrix3d rotatePointsToXY(std::vector< GeoLib::Point * > &pnts)

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 366 of file AnalyticalGeometry.cpp.

367 {
368  return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
369 }

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 446 of file AnalyticalGeometry.cpp.

447 {
448  // 1 copy all points
449  std::vector<GeoLib::Point*> polygon_points;
450  polygon_points.reserve(polygon_in.getNumberOfPoints());
451  for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
452  {
453  polygon_points.push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
454  }
455 
456  // 2 rotate points
457  auto [plane_normal, d_polygon] = GeoLib::getNewellPlane(polygon_points);
458  Eigen::Matrix3d const rot_mat =
459  GeoLib::computeRotationMatrixToXY(plane_normal);
460  GeoLib::rotatePoints(rot_mat, polygon_points);
461 
462  // 3 set z coord to zero
463  std::for_each(polygon_points.begin(), polygon_points.end(),
464  [](GeoLib::Point* p) { (*p)[2] = 0.0; });
465 
466  return {polygon_points, plane_normal};
467 }
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
static const double p

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 647 of file AnalyticalGeometry.cpp.

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

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 647 of file AnalyticalGeometry.cpp.

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

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 371 of file AnalyticalGeometry.cpp.

375 {
376  auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
377  auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
378  auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
379  auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
380  auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
381 
382  Eigen::Vector3d const pq = vq - vp;
383  Eigen::Vector3d const pa = va - vp;
384  Eigen::Vector3d const pb = vb - vp;
385  Eigen::Vector3d const pc = vc - vp;
386 
387  double u = pq.cross(pc).dot(pb);
388  if (u < 0)
389  {
390  return nullptr;
391  }
392  double v = pq.cross(pa).dot(pc);
393  if (v < 0)
394  {
395  return nullptr;
396  }
397  double w = pq.cross(pb).dot(pa);
398  if (w < 0)
399  {
400  return nullptr;
401  }
402 
403  const double denom(1.0 / (u + v + w));
404  u *= denom;
405  v *= denom;
406  w *= denom;
407  return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
408  u * a[1] + v * b[1] + w * c[1],
409  u * a[2] + v * b[2] + w * c[2]);
410 }

References MaterialPropertyLib::c, MathLib::TemplatePoint< T, DIM >::getCoords(), MathLib::q, MathLib::u, and MathLib::v.

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