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::Point >	triangleLineIntersection (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::Point3d >	lineSegmentIntersect2d (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::Point >	getAllIntersectionPoints (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< Polyline >	createPolyline (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

typedef TemplateVec< GeoLib::Polyline > GeoLib::PolylineVec = TemplateVec<GeoLib::Polyline>

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

Definition at line 16 of file PolylineVec.h.

SurfaceVec

using GeoLib::SurfaceVec = TemplateVec<GeoLib::Surface>

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

Definition at line 17 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 27 of file Polygon.cpp.

{
    TOUCHING,    
    CROSSING,    
    INESSENTIAL  
};

GEOTYPE

enum class GeoLib::GEOTYPE

strong

Enumerator
POINT
POLYLINE
SURFACE

Definition at line 11 of file GeoType.h.

{
    POINT,
    POLYLINE,
    SURFACE
};

Location

enum class GeoLib::Location

strong

Enumerator
LEFT
RIGHT
BEYOND
BEHIND
BETWEEN
SOURCE
DESTINATION

Definition at line 13 of file Polygon.cpp.

{
    LEFT,
    RIGHT,
    BEYOND,
    BEHIND,
    BETWEEN,
    SOURCE,
    DESTINATION
};

Orientation

enum GeoLib::Orientation

Enumerator
CW
COLLINEAR
CCW

Definition at line 13 of file AnalyticalGeometry.h.

{
    CW = -1,
    COLLINEAR = 0,
    CCW = 1
};

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

v	a 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 681 of file AnalyticalGeometry.cpp.

{
    Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
    const double cos_theta = v[0];
    const double sin_theta = v[1];
    rot_mat(0, 0) = rot_mat(1, 1) = cos_theta;
    rot_mat(0, 1) = sin_theta;
    rot_mat(1, 0) = -sin_theta;
    rot_mat(2, 2) = 1.0;
    return rot_mat;
}

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

v	a 3D normal vector to be rotated

Returns

a 3x3 rotation matrix

Definition at line 693 of file AnalyticalGeometry.cpp.

{
    // a vector on the plane
    Eigen::Vector3d yy = Eigen::Vector3d::Zero();
    auto const eps = std::numeric_limits<double>::epsilon();
    if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
    {
        yy[2] = 1.0;
    }
    else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
    {
        yy[0] = 1.0;
    }
    else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
    {
        yy[1] = 1.0;
    }
    else
    {
        for (unsigned i = 0; i < 3; i++)
        {
            if (std::abs(v[i]) > 0.0)
            {
                yy[i] = -v[i];
                break;
            }
        }
    }
    // z"_vec
    Eigen::Vector3d const zz = v.cross(yy).normalized();
    // y"_vec
    yy = zz.cross(v).normalized();
 
    Eigen::Matrix3d rot_mat;
    rot_mat.row(0) = v;
    rot_mat.row(1) = yy;
    rot_mat.row(2) = zz;
    return rot_mat;
}

Referenced by detail::getRotationMatrixToGlobal().

computeAndInsertAllIntersectionPoints() [1/2]

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

Definition at line 388 of file AnalyticalGeometry.cpp.

{
    auto computeSegmentIntersections =
        [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
    {
        for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
        {
            for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
            {
                GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
                if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
                {
                    std::size_t const id(
                        pnt_vec.push_back(new GeoLib::Point(s)));
                    poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
                    poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
                }
            }
        }
    };
 
    for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
    {
        auto it1(it0);
        ++it1;
        for (; it1 != plys.end(); ++it1)
        {
            computeSegmentIntersections(*(*it0), *(*it1));
        }
    }
}

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.

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

n	the (3d) vector that is rotated parallel to the \(z\) axis

Returns

rot_mat 3x3 rotation matrix

Definition at line 285 of file AnalyticalGeometry.cpp.

{
    Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
    // check if normal points already in the right direction
    if (n[0] == 0 && n[1] == 0)
    {
        rot_mat(1, 1) = 1.0;
 
        if (n[2] > 0)
        {
            // identity matrix
            rot_mat(0, 0) = 1.0;
            rot_mat(2, 2) = 1.0;
        }
        else
        {
            // rotate by pi about the y-axis
            rot_mat(0, 0) = -1.0;
            rot_mat(2, 2) = -1.0;
        }
 
        return rot_mat;
    }
 
    // sqrt (n_1^2 + n_3^2)
    double const h0(std::sqrt(n[0] * n[0] + n[2] * n[2]));
 
    // In case the x and z components of the normal are both zero the rotation
    // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
    // required. The angle is either pi/2 or 3/2*pi. Thus the components of
    // rot_mat are as follows.
    if (h0 < std::numeric_limits<double>::epsilon())
    {
        rot_mat(0, 0) = 1.0;
        if (n[1] > 0)
        {
            rot_mat(1, 2) = -1.0;
            rot_mat(2, 1) = 1.0;
        }
        else
        {
            rot_mat(1, 2) = 1.0;
            rot_mat(2, 1) = -1.0;
        }
        return rot_mat;
    }
 
    double const h1(1 / n.norm());
 
    // general case: calculate entries of rotation matrix
    rot_mat(0, 0) = n[2] / h0;
    rot_mat(0, 1) = 0;
    rot_mat(0, 2) = -n[0] / h0;
    rot_mat(1, 0) = -n[1] * n[0] / h0 * h1;
    rot_mat(1, 1) = h0 * h1;
    rot_mat(1, 2) = -n[1] * n[2] / h0 * h1;
    rot_mat(2, 0) = n[0] * h1;
    rot_mat(2, 1) = n[1] * h1;
    rot_mat(2, 2) = n[2] * h1;
 
    return rot_mat;
}

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

{
    if (id0 == id1)
    {
        ERR("no valid edge id0 == id1 == {:d}.", id0);
        return false;
    }
    if (id0 > id1)
    {
        std::swap(id0, id1);
    }
    const std::size_t n(ply.getNumberOfPoints() - 1);
    for (std::size_t k(0); k < n; k++)
    {
        std::size_t ply_pnt0(ply.getPointID(k));
        std::size_t ply_pnt1(ply.getPointID(k + 1));
        if (ply_pnt0 > ply_pnt1)
        {
            std::swap(ply_pnt0, ply_pnt1);
        }
        if (ply_pnt0 == id0 && ply_pnt1 == id1)
        {
            return true;
        }
    }
    return false;
}

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 12 of file GeoType.cpp.

{
    switch (geo_type)
    {
        case GEOTYPE::POINT:
            return "POINT";
        case GEOTYPE::POLYLINE:
            return "POLYLINE";
        case GEOTYPE::SURFACE:
            return "SURFACE";
    }
 
    // Cannot happen, because switch covers all cases.
    // Used to silence compiler warning.
    OGS_FATAL("convertGeoTypeToString(): Given geo type is not supported");
}

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

{
    for (auto it0 = list_of_simple_polygon_hierarchies.begin();
         it0 != list_of_simple_polygon_hierarchies.end();
         ++it0)
    {
        auto it1 = list_of_simple_polygon_hierarchies.begin();
        while (it1 != list_of_simple_polygon_hierarchies.end())
        {
            if (it0 == it1)
            {  // don't check same polygons
                ++it1;
                // skip test if it1 points to the end after increment
                if (it1 == list_of_simple_polygon_hierarchies.end())
                {
                    break;
                }
            }
            if ((*it0)->isPolygonInside(*it1))
            {
                (*it0)->insertSimplePolygonTree(*it1);
                it1 = list_of_simple_polygon_hierarchies.erase(it1);
            }
            else
            {
                ++it1;
            }
        }
    }
}

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

{
    auto const& point_id_map = points_vec.getIDMap();
    auto polyline = std::make_unique<Polyline>(points_vec.getVector());
    for (auto point_id : point_ids)
    {
        if (point_id >= point_id_map.size())
        {
            WARN("The point id {} doesn't exist in the underlying PointVec.",
                 point_id);
            continue;
        }
        polyline->addPoint(point_id_map[point_id]);
    }
    return polyline;
}

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

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

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 10 of file GeoLib/Utils.cpp.

{
    if (number_of_subdivisions < 0)
    {
        OGS_FATAL(
            "generateEquidistantPoints: number of subdivisions is required to "
            "be non-negative.");
    }
 
    auto const& start = begin.asEigenVector3d();
    auto const& stop = end.asEigenVector3d();
    auto const delta = (stop - start) / (number_of_subdivisions + 1);
 
    std::vector<GeoLib::Point*> points;
 
    for (int i = 0; i <= number_of_subdivisions; ++i)
    {
        auto const p = start + i * delta;
        points.push_back(new GeoLib::Point{p[0], p[1], p[2]});
    }
    points.push_back(new GeoLib::Point{stop[0], stop[1], stop[2]});
 
    return points;
}

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

{
    GeoLib::PointVec const* const pnt_obj(geo_objects.getPointVecObj(geo_name));
    if (pnt_obj == nullptr)
    {
        ERR("Point vector {:s} not found.", geo_name);
        return -1;
    }
    auto const& pnts = pnt_obj->getVector();
    if (pnts.empty())
    {
        ERR("Point vector {:s} is empty.", geo_name);
        return -1;
    }
    std::size_t const n_pnts(pnts.size());
    std::vector<bool> transfer_pnts(n_pnts, true);
    if (only_unused_pnts)
    {
        markUnusedPoints(geo_objects, geo_name, transfer_pnts);
    }
 
    std::vector<GeoLib::Point*> stations;
    for (std::size_t i = 0; i < n_pnts; ++i)
    {
        if (!transfer_pnts[i])
        {
            continue;
        }
        std::string name = pnt_obj->getItemNameByID(i);
        if (name.empty())
        {
            name = "Station " + std::to_string(i);
        }
        stations.push_back(new GeoLib::Station(pnts[i], name));
    }
    if (!stations.empty())
    {
        geo_objects.addStationVec(std::move(stations), stn_name);
    }
    else
    {
        WARN("No points found to convert.");
        return 1;
    }
    return 0;
}

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

polygon	the polygon the segment line segment that will be processed
segment	the line segment that will be processed

Returns

a possible empty vector containing the intersection points

Definition at line 178 of file Polygon.cpp.

{
    std::vector<GeoLib::Point> intersections;
    GeoLib::Point s;
    for (auto&& seg_it : polygon)
    {
        if (GeoLib::lineSegmentIntersect(seg_it, segment, s))
        {
            intersections.push_back(s);
        }
    }
 
    return intersections;
}

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

a	first point of line segment
b	last point of line segment
pnt	point that is edge type computed for

Returns

a value of enum EdgeType

Definition at line 91 of file Polygon.cpp.

{
    switch (getLocationOfPoint(a, b, pnt))
    {
        case Location::LEFT:
        {
            if (a[1] < pnt[1] && pnt[1] <= b[1])
            {
                return EdgeType::CROSSING;
            }
 
            return EdgeType::INESSENTIAL;
        }
        case Location::RIGHT:
        {
            if (b[1] < pnt[1] && pnt[1] <= a[1])
            {
                return EdgeType::CROSSING;
            }
 
            return EdgeType::INESSENTIAL;
        }
        case Location::BETWEEN:
        case Location::SOURCE:
        case Location::DESTINATION:
            return EdgeType::TOUCHING;
        default:
            return EdgeType::INESSENTIAL;
    }
}

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

source	the first point of the line segment
destination	the end point of the line segment
pnt	the test point

Returns

a value of enum Location

Definition at line 44 of file Polygon.cpp.

{
    long double const a[2] = {destination[0] - source[0],
                              destination[1] - source[1]};  // vector
    long double const b[2] = {pnt[0] - source[0],
                              pnt[1] - source[1]};  // vector
 
    long double const det_2x2(a[0] * b[1] - a[1] * b[0]);
    constexpr double eps = std::numeric_limits<double>::epsilon();
 
    if (det_2x2 > eps)
    {
        return Location::LEFT;
    }
    if (eps < std::abs(det_2x2))
    {
        return Location::RIGHT;
    }
    if (a[0] * b[0] < 0.0 || a[1] * b[1] < 0.0)
    {
        return Location::BEHIND;
    }
    if (a[0] * a[0] + a[1] * a[1] < b[0] * b[0] + b[1] * b[1])
    {
        return Location::BEYOND;
    }
    if (MathLib::sqrDist(pnt, source) < pow(eps, 2))
    {
        return Location::SOURCE;
    }
    if (MathLib::sqrDist(pnt, destination) < std::sqrt(eps))
    {
        return Location::DESTINATION;
    }
    return Location::BETWEEN;
}

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

pnts	points 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 39 of file AnalyticalGeometry-impl.h.

{
    return getNewellPlane(pnts.begin(), 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 46 of file AnalyticalGeometry-impl.h.

{
    Eigen::Vector3d plane_normal({0, 0, 0});
    Eigen::Vector3d centroid({0, 0, 0});
    std::size_t n_pnts(pnts.size());
    for (std::size_t i = n_pnts - 1, j = 0; j < n_pnts; i = j, j++)
    {
        plane_normal[0] += (pnts[i][1] - pnts[j][1]) *
                           (pnts[i][2] + pnts[j][2]);  // projection on yz
        plane_normal[1] += (pnts[i][2] - pnts[j][2]) *
                           (pnts[i][0] + pnts[j][0]);  // projection on xz
        plane_normal[2] += (pnts[i][0] - pnts[j][0]) *
                           (pnts[i][1] + pnts[j][1]);  // projection on xy
 
        centroid += pnts[j].asEigenVector3d();
    }
 
    plane_normal.normalize();
    double const d = centroid.dot(plane_normal) / n_pnts;
    return std::make_pair(plane_normal, d);
}

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_begin	Iterator pointing to the initial point of a closed polyline describing a polygon
pnts_end	Iterator 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 7 of file AnalyticalGeometry-impl.h.

{
    Eigen::Vector3d plane_normal({0, 0, 0});
    Eigen::Vector3d centroid({0, 0, 0});
    for (auto i = std::prev(pnts_end), j = pnts_begin; j != pnts_end;
         i = j, ++j)
    {
        auto const& pt_i = *(*i);
        auto const& pt_j = *(*j);
        plane_normal[0] +=
            (pt_i[1] - pt_j[1]) * (pt_i[2] + pt_j[2]);  // projection on yz
        plane_normal[1] +=
            (pt_i[2] - pt_j[2]) * (pt_i[0] + pt_j[0]);  // projection on xz
        plane_normal[2] +=
            (pt_i[0] - pt_j[0]) * (pt_i[1] + pt_j[1]);  // projection on xy
 
        centroid += pt_j.asEigenVector3d();
    }
 
    plane_normal.normalize();
    auto n_pnts(std::distance(pnts_begin, pnts_end));
    assert(n_pnts > 2);
    double d = 0.0;
    if (n_pnts > 0)
    {
        d = centroid.dot(plane_normal) / n_pnts;
    }
    return std::make_pair(plane_normal, d);
}

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

{
    double const orientation = getOrientation2d(p0, p1, p2);
    if (orientation > 0)
    {
        return CCW;
    }
    if (orientation < 0)
    {
        return CW;
    }
    return COLLINEAR;
}

References CCW, COLLINEAR, and CW.

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

{
    double const orientation = getOrientation2dFast(p0, p1, p2);
    if (orientation > 0)
    {
        return CCW;
    }
    if (orientation < 0)
    {
        return CW;
    }
    return COLLINEAR;
}

References CCW, COLLINEAR, and CW.

isBorehole()

bool GeoLib::isBorehole

(

GeoLib::Point const *

pnt

)

Definition at line 95 of file StationBorehole.cpp.

{
    auto const* bh(dynamic_cast<GeoLib::StationBorehole const*>(pnt));
    return bh != nullptr;
}

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

isStation()

bool GeoLib::isStation

(

GeoLib::Point const *

pnt

)

Definition at line 66 of file Station.cpp.

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

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

s0	the first line segment.
s1	the second line segment.
s	the intersection point if the segments do intersect

Returns

true, if the line segments intersect, else false

Definition at line 127 of file AnalyticalGeometry.cpp.

{
    GeoLib::Point const& pa{s0.getBeginPoint()};
    GeoLib::Point const& pb{s0.getEndPoint()};
    GeoLib::Point const& pc{s1.getBeginPoint()};
    GeoLib::Point const& pd{s1.getEndPoint()};
 
    if (!isCoplanar(pa, pb, pc, pd))
    {
        return false;
    }
 
    auto const& a = s0.getBeginPoint().asEigenVector3d();
    auto const& b = s0.getEndPoint().asEigenVector3d();
    auto const& c = s1.getBeginPoint().asEigenVector3d();
    auto const& d = s1.getEndPoint().asEigenVector3d();
 
    Eigen::Vector3d const v = b - a;
    Eigen::Vector3d const w = d - c;
    Eigen::Vector3d const qp = c - a;
    Eigen::Vector3d const pq = a - c;
 
    double const eps = std::numeric_limits<double>::epsilon();
    double const squared_eps = eps * eps;
    // handle special cases here to avoid computing intersection numerical
    if (qp.squaredNorm() < squared_eps || (d - a).squaredNorm() < squared_eps)
    {
        s = pa;
        return true;
    }
    if ((c - b).squaredNorm() < squared_eps ||
        (d - b).squaredNorm() < squared_eps)
    {
        s = pb;
        return true;
    }
 
    auto isLineSegmentIntersectingAB =
        [&v](Eigen::Vector3d const& ap, std::size_t i)
    {
        // check if p is located at v=(a,b): (ap = t*v, t in [0,1])
        return 0.0 <= ap[i] / v[i] && ap[i] / v[i] <= 1.0;
    };
 
    if (parallel(v, w))
    {  // original line segments (a,b) and (c,d) are parallel
        if (parallel(pq, v))
        {  // line segment (a,b) and (a,c) are also parallel
            // Here it is already checked that the line segments (a,b) and (c,d)
            // are parallel. At this point it is also known that the line
            // segment (a,c) is also parallel to (a,b). In that case it is
            // possible to express c as c(t) = a + t * (b-a) (analog for the
            // point d). Since the evaluation of all three coordinate equations
            // (x,y,z) have to lead to the same solution for the parameter t it
            // is sufficient to evaluate t only once.
 
            // Search id of coordinate with largest absolute value which is will
            // be used in the subsequent computations. This prevents division by
            // zero in case the line segments are parallel to one of the
            // coordinate axis.
            std::size_t i_max(std::abs(v[0]) <= std::abs(v[1]) ? 1 : 0);
            i_max = std::abs(v[i_max]) <= std::abs(v[2]) ? 2 : i_max;
            if (isLineSegmentIntersectingAB(qp, i_max))
            {
                s = pc;
                return true;
            }
            Eigen::Vector3d const ad = d - a;
            if (isLineSegmentIntersectingAB(ad, i_max))
            {
                s = pd;
                return true;
            }
            return false;
        }
        return false;
    }
 
    // general case
    const double sqr_len_v(v.squaredNorm());
    const double sqr_len_w(w.squaredNorm());
 
    Eigen::Matrix2d mat;
    mat(0, 0) = sqr_len_v;
    mat(0, 1) = -v.dot(w);
    mat(1, 1) = sqr_len_w;
    mat(1, 0) = mat(0, 1);
 
    Eigen::Vector2d rhs{v.dot(qp), w.dot(pq)};
 
    rhs = mat.partialPivLu().solve(rhs);
 
    // no theory for the following tolerances, determined by testing
    // lower tolerance: little bit smaller than zero
    const double l(-1.0 * std::numeric_limits<float>::epsilon());
    // upper tolerance a little bit greater than one
    const double u(1.0 + std::numeric_limits<float>::epsilon());
    if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1])
    {
        return false;
    }
 
    // compute points along line segments with minimal distance
    GeoLib::Point const p0(a[0] + rhs[0] * v[0], a[1] + rhs[0] * v[1],
                           a[2] + rhs[0] * v[2]);
    GeoLib::Point const p1(c[0] + rhs[1] * w[0], c[1] + rhs[1] * w[1],
                           c[2] + rhs[1] * w[2]);
 
    double const min_dist(std::sqrt(MathLib::sqrDist(p0, p1)));
    double const min_seg_len(
        std::min(std::sqrt(sqr_len_v), std::sqrt(sqr_len_w)));
    if (min_dist < min_seg_len * 1e-6)
    {
        s[0] = 0.5 * (p0[0] + p1[0]);
        s[1] = 0.5 * (p0[1] + p1[1]);
        s[2] = 0.5 * (p0[2] + p1[2]);
        return true;
    }
 
    return false;
}

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

ab	first line segment
cd	second 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 445 of file AnalyticalGeometry.cpp.

{
    GeoLib::Point const& a{ab.getBeginPoint()};
    GeoLib::Point const& b{ab.getEndPoint()};
    GeoLib::Point const& c{cd.getBeginPoint()};
    GeoLib::Point const& d{cd.getEndPoint()};
 
    double const orient_abc(getOrientation(a, b, c));
    double const orient_abd(getOrientation(a, b, d));
 
    // check if the segment (cd) lies on the left or on the right of (ab)
    if ((orient_abc > 0 && orient_abd > 0) ||
        (orient_abc < 0 && orient_abd < 0))
    {
        return std::vector<MathLib::Point3d>();
    }
 
    // check: (cd) and (ab) are on the same line
    if (orient_abc == 0.0 && orient_abd == 0.0)
    {
        double const eps(std::numeric_limits<double>::epsilon());
        if (MathLib::sqrDist2d(a, c) < eps && MathLib::sqrDist2d(b, d) < eps)
        {
            return {{a, b}};
        }
        if (MathLib::sqrDist2d(a, d) < eps && MathLib::sqrDist2d(b, c) < eps)
        {
            return {{a, b}};
        }
 
        // Since orient_ab and orient_abd vanish, a, b, c, d are on the same
        // line and for this reason it is enough to check the x-component.
        auto isPointOnSegment = [](double q, double p0, double p1)
        {
            double const t((q - p0) / (p1 - p0));
            return 0 <= t && t <= 1;
        };
 
        // check if c in (ab)
        if (isPointOnSegment(c[0], a[0], b[0]))
        {
            // check if a in (cd)
            if (isPointOnSegment(a[0], c[0], d[0]))
            {
                return {{a, c}};
            }
            // check b == c
            if (MathLib::sqrDist2d(b, c) < eps)
            {
                return {{b}};
            }
            // check if b in (cd)
            if (isPointOnSegment(b[0], c[0], d[0]))
            {
                return {{b, c}};
            }
            // check d in (ab)
            if (isPointOnSegment(d[0], a[0], b[0]))
            {
                return {{c, d}};
            }
            std::stringstream err;
            err.precision(std::numeric_limits<double>::max_digits10);
            err << ab << " x " << cd;
            OGS_FATAL(
                "The case of parallel line segments ({:s}) is not handled yet. "
                "Aborting.",
                err.str());
        }
 
        // check if d in (ab)
        if (isPointOnSegment(d[0], a[0], b[0]))
        {
            // check if a in (cd)
            if (isPointOnSegment(a[0], c[0], d[0]))
            {
                return {{a, d}};
            }
            // check if b==d
            if (MathLib::sqrDist2d(b, d) < eps)
            {
                return {{b}};
            }
            // check if b in (cd)
            if (isPointOnSegment(b[0], c[0], d[0]))
            {
                return {{b, d}};
            }
            // d in (ab), b not in (cd): check c in (ab)
            if (isPointOnSegment(c[0], a[0], b[0]))
            {
                return {{c, d}};
            }
 
            std::stringstream err;
            err.precision(std::numeric_limits<double>::max_digits10);
            err << ab << " x " << cd;
            OGS_FATAL(
                "The case of parallel line segments ({:s}) is not handled yet. "
                "Aborting.",
                err.str());
        }
        return std::vector<MathLib::Point3d>();
    }
 
    // precondition: points a, b, c are collinear
    // the function checks if the point c is onto the line segment (a,b)
    auto isCollinearPointOntoLineSegment = [](MathLib::Point3d const& a,
                                              MathLib::Point3d const& b,
                                              MathLib::Point3d const& c)
    {
        if (b[0] - a[0] != 0)
        {
            double const t = (c[0] - a[0]) / (b[0] - a[0]);
            return 0.0 <= t && t <= 1.0;
        }
        if (b[1] - a[1] != 0)
        {
            double const t = (c[1] - a[1]) / (b[1] - a[1]);
            return 0.0 <= t && t <= 1.0;
        }
        if (b[2] - a[2] != 0)
        {
            double const t = (c[2] - a[2]) / (b[2] - a[2]);
            return 0.0 <= t && t <= 1.0;
        }
        return false;
    };
 
    if (orient_abc == 0.0)
    {
        if (isCollinearPointOntoLineSegment(a, b, c))
        {
            return {{c}};
        }
        return std::vector<MathLib::Point3d>();
    }
 
    if (orient_abd == 0.0)
    {
        if (isCollinearPointOntoLineSegment(a, b, d))
        {
            return {{d}};
        }
        return std::vector<MathLib::Point3d>();
    }
 
    // check if the segment (ab) lies on the left or on the right of (cd)
    double const orient_cda(getOrientation(c, d, a));
    double const orient_cdb(getOrientation(c, d, b));
    if ((orient_cda > 0 && orient_cdb > 0) ||
        (orient_cda < 0 && orient_cdb < 0))
    {
        return std::vector<MathLib::Point3d>();
    }
 
    // at this point it is sure that there is an intersection and the system of
    // linear equations will be invertible
    // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously
    Eigen::Matrix2d mat;
    mat(0, 0) = b[0] - a[0];
    mat(0, 1) = c[0] - d[0];
    mat(1, 0) = b[1] - a[1];
    mat(1, 1) = c[1] - d[1];
    Eigen::Vector2d rhs{c[0] - a[0], c[1] - a[1]};
 
    rhs = mat.partialPivLu().solve(rhs);
    if (0 <= rhs[1] && rhs[1] <= 1.0)
    {
        return {MathLib::Point3d{std::array<double, 3>{
            {c[0] + rhs[1] * (d[0] - c[0]), c[1] + rhs[1] * (d[1] - c[1]),
             c[2] + rhs[1] * (d[2] - c[2])}}}};
    }
    return std::vector<MathLib::Point3d>();  // parameter s not in the valid
                                             // range
}

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

ply	the polyline
seg_it0	iterator pointing to the first segment that has an intersection
seg_it1	iterator pointing to the second segment that has an intersection
intersection_pnt	the intersection point if the segments intersect

Returns

true, if the polyline contains intersections

Definition at line 251 of file AnalyticalGeometry.cpp.

{
    std::size_t const n_segs(ply->getNumberOfSegments());
    // Neighbouring segments always intersects at a common vertex. The algorithm
    // checks for intersections of non-neighbouring segments.
    for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0)
    {
        seg_it1 = seg_it0 + 2;
        std::size_t const seg_num_0 = seg_it0.getSegmentNumber();
        for (; seg_it1 != ply->end(); ++seg_it1)
        {
            // Do not check first and last segment, because they are
            // neighboured.
            if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1))
            {
                if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt))
                {
                    return true;
                }
            }
        }
    }
    return false;
}

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

{
    GeoLib::PolylineVec const* const ply_obj(
        geo_objects.getPolylineVecObj(geo_name));
    if (ply_obj)
    {
        std::vector<GeoLib::Polyline*> const& lines(ply_obj->getVector());
        for (auto* line : lines)
        {
            std::size_t const n_pnts(line->getNumberOfPoints());
            for (std::size_t i = 0; i < n_pnts; ++i)
            {
                transfer_pnts[line->getPointID(i)] = false;
            }
        }
    }
 
    GeoLib::SurfaceVec const* const sfc_obj(
        geo_objects.getSurfaceVecObj(geo_name));
    if (sfc_obj)
    {
        std::vector<GeoLib::Surface*> const& surfaces = sfc_obj->getVector();
        for (auto* sfc : surfaces)
        {
            std::size_t const n_tri(sfc->getNumberOfTriangles());
            for (std::size_t i = 0; i < n_tri; ++i)
            {
                GeoLib::Triangle const& t = *(*sfc)[i];
                transfer_pnts[t[0]] = false;
                transfer_pnts[t[1]] = false;
                transfer_pnts[t[2]] = false;
            }
        }
    }
}

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

{
    if (polyline.getPointsVec().size() != used_points.size())
    {
        OGS_FATAL(
            "internal error in markUsedPoints(): polyline based on point "
            "vector of size {}, given used_points has size {}",
            polyline.getPointsVec().size(), used_points.size());
    }
    for (std::size_t i = 0; i < polyline.getNumberOfPoints(); ++i)
    {
        used_points[polyline.getPointID(i)] = true;
    }
}

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 126 of file Surface.cpp.

{
    if (surface.getPointVec()->size() != used_points.size())
    {
        OGS_FATAL(
            "internal error in markUsedPoints(): surface based on point vector "
            "of size {}, given used_points has size {}",
            surface.getPointVec()->size(), used_points.size());
    }
    for (std::size_t i = 0; i < surface.getNumberOfTriangles(); ++i)
    {
        auto const& triangle = *surface[i];
        for (std::size_t k = 0; k < 3; ++k)
        {
            used_points[triangle[k]] = true;
        }
    }
}

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 80 of file LineSegment.cpp.

{
    os << "{(" << s.getBeginPoint() << "), (" << s.getEndPoint() << ")}";
    return os;
}

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 86 of file LineSegment.cpp.

{
    os << seg_pair.first << " x " << seg_pair.second;
    return os;
}

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 94 of file LineSegment.cpp.

{
    double const tol(std::numeric_limits<double>::epsilon());
    return (MathLib::sqrDist(s0.getBeginPoint(), s1.getBeginPoint()) < tol &&
            MathLib::sqrDist(s0.getEndPoint(), s1.getEndPoint()) < tol) ||
           (MathLib::sqrDist(s0.getBeginPoint(), s1.getEndPoint()) < tol &&
            MathLib::sqrDist(s0.getEndPoint(), s1.getBeginPoint()) < tol);
}

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

lhs	the first polygon
rhs	the second polygon

Returns

true, if the polygons describe the same geometrical object

Definition at line 585 of file Polygon.cpp.

{
    if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
    {
        return false;
    }
 
    const std::size_t n(lhs.getNumberOfPoints());
    const std::size_t start_pnt(lhs.getPointID(0));
 
    // search start point of first polygon in second polygon
    bool nfound(true);
    std::size_t k(0);
    for (; k < n - 1 && nfound; k++)
    {
        if (start_pnt == rhs.getPointID(k))
        {
            nfound = false;
            break;
        }
    }
 
    // case: start point not found in second polygon
    if (nfound)
    {
        return false;
    }
 
    // *** determine direction
    // opposite direction
    if (k == n - 2)
    {
        for (k = 1; k < n - 1; k++)
        {
            if (lhs.getPointID(k) != rhs.getPointID(n - 1 - k))
            {
                return false;
            }
        }
        return true;
    }
 
    // same direction - start point of first polygon at arbitrary position in
    // second polygon
    if (lhs.getPointID(1) == rhs.getPointID(k + 1))
    {
        std::size_t j(k + 2);
        for (; j < n - 1; j++)
        {
            if (lhs.getPointID(j - k) != rhs.getPointID(j))
            {
                return false;
            }
        }
        j = 0;  // new start point at second polygon
        for (; j < k + 1; j++)
        {
            if (lhs.getPointID(n - (k + 2) + j + 1) != rhs.getPointID(j))
            {
                return false;
            }
        }
        return true;
    }
    // opposite direction with start point of first polygon at arbitrary
    // position
    // *** ATTENTION
    WARN(
        "operator==(Polygon const& lhs, Polygon const& rhs) - not tested case "
        "(implementation is probably buggy) - please contact "
        "thomas.fischer@ufz.de mentioning the problem.");
    // in second polygon
    if (lhs.getPointID(1) == rhs.getPointID(k - 1))
    {
        std::size_t j(k - 2);
        for (; j > 0; j--)
        {
            if (lhs.getPointID(k - 2 - j) != rhs.getPointID(j))
            {
                return false;
            }
        }
        // new start point at second polygon - the point n-1 of a polygon is
        // equal to the first point of the polygon (for this reason: n-2)
        j = n - 2;
        for (; j > k - 1; j--)
        {
            if (lhs.getPointID(n - 2 + j + k - 2) != rhs.getPointID(j))
            {
                return false;
            }
        }
        return true;
    }
    return false;
}

operator==() [3/4]

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

comparison operator

Parameters

lhs	first polyline
rhs	second polyline

Returns

true, if the polylines consists of the same sequence of line segments

Definition at line 522 of file Polyline.cpp.

{
    if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
    {
        return false;
    }
 
    const std::size_t n(lhs.getNumberOfPoints());
    for (std::size_t k(0); k < n; k++)
    {
        if (lhs.getPointID(k) != rhs.getPointID(k))
        {
            return false;
        }
    }
 
    return true;
}

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

operator==() [4/4]

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

Definition at line 103 of file Surface.cpp.

{
    return &lhs == &rhs;
}

parallel()

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

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

Parameters

v	first vector
w	second vector

Returns

true if the vectors are in parallel, else false

Definition at line 71 of file AnalyticalGeometry.cpp.

{
    const double eps(std::numeric_limits<double>::epsilon());
    double const eps_squared = eps * eps;
 
    // check degenerated cases
    if (v.squaredNorm() < eps_squared)
    {
        return false;
    }
 
    if (w.squaredNorm() < eps_squared)
    {
        return false;
    }
 
    v.normalize();
    w.normalize();
 
    bool parallel(true);
    if (std::abs(v[0] - w[0]) > eps)
    {
        parallel = false;
    }
    if (std::abs(v[1] - w[1]) > eps)
    {
        parallel = false;
    }
    if (std::abs(v[2] - w[2]) > eps)
    {
        parallel = false;
    }
 
    if (!parallel)
    {
        parallel = true;
        // change sense of direction of v_normalised
        v *= -1.0;
        // check again
        if (std::abs(v[0] - w[0]) > eps)
        {
            parallel = false;
        }
        if (std::abs(v[1] - w[1]) > eps)
        {
            parallel = false;
        }
        if (std::abs(v[2] - w[2]) > eps)
        {
            parallel = false;
        }
    }
 
    return parallel;
}

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

{
    if (i == j)
    {
        return true;
    }
    return MathLib::sqrDist(*pnt_vec[i], *pnt_vec[j]) < prox;
}

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

{
    if (polyline.getPointsVec().size() != mapping.size())
    {
        OGS_FATAL(
            "internal error in resetPointIDs(): polyline based on point vector "
            "of size {}, given mapping has size {}",
            polyline.getPointsVec().size(), mapping.size());
    }
    for (std::size_t i = 0; i < polyline.getNumberOfPoints(); ++i)
    {
        polyline.setPointID(i, mapping[polyline.getPointID(i)]);
    }
}

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 108 of file Surface.cpp.

{
    if (surface.getPointVec()->size() != mapping.size())
    {
        OGS_FATAL(
            "internal error in resetPointIDs(): surface based on point vector "
            "of size {}, given mapping vector has size {}",
            surface.getPointVec()->size(), mapping.size());
    }
    for (std::size_t i = 0; i < surface.getNumberOfTriangles(); ++i)
    {
        auto& triangle = *surface[i];
        const_cast<std::size_t&>(triangle[0]) = mapping[triangle[0]];
        const_cast<std::size_t&>(triangle[1]) = mapping[triangle[1]];
        const_cast<std::size_t&>(triangle[2]) = mapping[triangle[2]];
    }
}

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_mat	3x3 dimensional rotation matrix
pnts_begin	Iterator pointing to the initial element in a vector of points to be rotated
pnts_end	Iterator pointing to the element following the last element in a vector of points to be rotated

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

{
    for (auto it = pnts_begin; it != pnts_end; ++it)
    {
        (*it)->asEigenVector3d() = rot_mat * (*it)->asEigenVector3d();
    }
}

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

{
    rotatePoints(rot_mat, pnts.begin(), 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_mat	3x3 dimensional rotation matrix
pnts	vector of points

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

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

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_begin	Iterator pointing to the initial element in a vector of points used for computing a rotation matrix
p_pnts_end	Iterator pointing to the element following the last point in a vector of points used for computing a rotation matrix
r_pnts_begin	Iterator pointing to the initial element in a vector of points to be rotated
r_pnts_end	Iterator 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 80 of file AnalyticalGeometry-impl.h.

{
    assert(std::distance(p_pnts_begin, p_pnts_end) > 2);
 
    // compute the plane normal
    auto const [plane_normal, d] =
        GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
 
    // rotate points into x-y-plane
    Eigen::Matrix3d const rot_mat = computeRotationMatrixToXY(plane_normal);
    rotatePoints(rot_mat, r_pnts_begin, r_pnts_end);
 
    for (auto it = r_pnts_begin; it != r_pnts_end; ++it)
    {
        (*(*it))[2] = 0.0;  // should be -= d but there are numerical errors
    }
 
    return rot_mat;
}

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< GeoLib::Point * > &

pnts

)

Definition at line 348 of file AnalyticalGeometry.cpp.

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

References rotatePointsToXY().

rotatePointsToXY() [3/3]

Eigen::Matrix3d GeoLib::rotatePointsToXY

(

std::vector< Point * > &

pnts

)

rotate points to X-Y plane

Parameters

pnts	a 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.

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_in

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

{
    // 1 copy all points
    std::vector<GeoLib::Point*> polygon_points;
    polygon_points.reserve(polygon_in.getNumberOfPoints());
    for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
    {
        polygon_points.push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
    }
 
    // 2 rotate points
    auto [plane_normal, d_polygon] = GeoLib::getNewellPlane(polygon_points);
    Eigen::Matrix3d const rot_mat =
        GeoLib::computeRotationMatrixToXY(plane_normal);
    GeoLib::rotatePoints(rot_mat, polygon_points);
 
    // 3 set z coord to zero
    std::for_each(polygon_points.begin(), polygon_points.end(),
                  [](GeoLib::Point* p) { (*p)[2] = 0.0; });
 
    return {polygon_points, plane_normal};
}

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< GeoLib::LineSegment > &	sub_segments )

Definition at line 623 of file AnalyticalGeometry.cpp.

{
    double const eps(std::numeric_limits<double>::epsilon());
 
    auto findNextSegment =
        [&eps](MathLib::Point3d const& seg_beg_pnt,
               std::vector<GeoLib::LineSegment>& sub_segments,
               std::vector<GeoLib::LineSegment>::iterator& sub_seg_it)
    {
        if (sub_seg_it == sub_segments.end())
        {
            return;
        }
        // find appropriate segment for the given segment begin point
        auto act_beg_seg_it = std::find_if(
            sub_seg_it, sub_segments.end(),
            [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
            {
                return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) <
                           eps ||
                       MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
            });
        if (act_beg_seg_it == sub_segments.end())
        {
            return;
        }
        // if necessary correct orientation of segment, i.e. swap beg and
        // end
        if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
            MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
        {
            std::swap(act_beg_seg_it->getBeginPoint(),
                      act_beg_seg_it->getEndPoint());
        }
        assert(sub_seg_it != sub_segments.end());
        // exchange segments within the container
        if (sub_seg_it != act_beg_seg_it)
        {
            std::swap(*sub_seg_it, *act_beg_seg_it);
        }
    };
 
    // find start segment
    auto seg_it = sub_segments.begin();
    findNextSegment(seg_beg_pnt, sub_segments, seg_it);
 
    while (seg_it != sub_segments.end())
    {
        MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint());
        seg_it++;
        if (seg_it != sub_segments.end())
        {
            findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
        }
    }
}

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.

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

{
    Eigen::Vector3d const pq = q.asEigenVector3d() - p.asEigenVector3d();
    Eigen::Vector3d const pa = a.asEigenVector3d() - p.asEigenVector3d();
    Eigen::Vector3d const pb = b.asEigenVector3d() - p.asEigenVector3d();
    Eigen::Vector3d const pc = c.asEigenVector3d() - p.asEigenVector3d();
 
    double u = pq.cross(pc).dot(pb);
    if (u < 0)
    {
        return nullptr;
    }
    double v = pq.cross(pa).dot(pc);
    if (v < 0)
    {
        return nullptr;
    }
    double w = pq.cross(pb).dot(pa);
    if (w < 0)
    {
        return nullptr;
    }
 
    const double denom(1.0 / (u + v + w));
    u *= denom;
    v *= denom;
    w *= denom;
    return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
                                           u * a[1] + v * b[1] + w * c[1],
                                           u * a[2] + v * b[2] + w * c[2]);
}

References MathLib::Point3d::asEigenVector3d().

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