GeoLib Namespace Reference

Detailed Description

Copyright

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

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.

Copyright

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

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 deleting of points. Additionally, each vector of points is identified by a unique name from class GEOObject. For this reason PointVec should have a 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::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)
GeoLib::Polygon	rotatePolygonToXY (GeoLib::Polygon const &polygon_in, Eigen::Vector3d &plane_normal)
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)
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

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

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.

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

Orientation

enum GeoLib::Orientation

Enumerator
CW
COLLINEAR
CCW

Definition at line 24 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 712 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 724 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	(	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.

{
    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.

Definition at line 412 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().

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 303 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 rotatePolygonToXY().

containsEdge()

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

Definition at line 510 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 23 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 97 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;
            }
        }
    }
}

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 698 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;
    }
    std::vector<GeoLib::Point*> 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);
    }
 
    auto stations = std::make_unique<std::vector<GeoLib::Point*>>();
    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(), 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 [10] .

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 44 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 51 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 += Eigen::Map<Eigen::Vector3d const>(pnts[j].getCoords());
    }
 
    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 [10] .

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 14 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 &pt_i = *(*i);
        auto &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 += Eigen::Map<Eigen::Vector3d const>(pt_j.getCoords());
    }
 
    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 rotatePolygonToXY().

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.

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

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.

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

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.

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

Referenced by MeshGeoToolsLib::GeoMapper::mapStationData(), and VtkStationSource::RequestData().

isStation()

bool GeoLib::isStation

(

GeoLib::Point const *

pnt

)

Definition at line 76 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 141 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 =
        Eigen::Map<Eigen::Vector3d const>(s0.getBeginPoint().getCoords());
    auto const b =
        Eigen::Map<Eigen::Vector3d const>(s0.getEndPoint().getCoords());
    auto const c =
        Eigen::Map<Eigen::Vector3d const>(s1.getBeginPoint().getCoords());
    auto const d =
        Eigen::Map<Eigen::Vector3d const>(s1.getEndPoint().getCoords());
 
    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 MaterialPropertyLib::c, GeoLib::LineSegment::getBeginPoint(), MathLib::TemplatePoint< T, DIM >::getCoords(), GeoLib::LineSegment::getEndPoint(), MathLib::isCoplanar(), parallel(), and MathLib::sqrDist().

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

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 476 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>::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>::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 MaterialPropertyLib::c, GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), getOrientation(), OGS_FATAL, MathLib::q, 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 269 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 660 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().

operator<<() [1/2]

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

Definition at line 86 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 92 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 100 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 521 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 538 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 109 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 85 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;
}

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 557 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().

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

{
    for (auto it = pnts_begin; it != pnts_end; ++it)
    {
        Eigen::Map<Eigen::Vector3d>((*it)->getCoords()) =
            rot_mat * Eigen::Map<Eigen::Vector3d const>((*it)->getCoords());
    }
}

Referenced by MeshGeoToolsLib::markNodesOutSideOfPolygon(), rotateGeometryToXY(), rotateMesh(), rotatePoints(), rotatePointsToXY(), rotatePolygonToXY(), 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.

{
    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 110 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 86 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< 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.

Definition at line 366 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.

Definition at line 366 of file AnalyticalGeometry.cpp.

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

References rotatePointsToXY().

rotatePolygonToXY()

Polygon GeoLib::rotatePolygonToXY	(	Polygon const &	polygon_in,
		Eigen::Vector3d &	plane_normal
	)

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 security the $z$ coordinates of the rotated points are set to zero and finally, (4) a new polygon is constructed using the rotated points.

See also

getNewellPlane()

Parameters

polygon_in	a copy of the polygon_in polygon will be rotated
plane_normal	the normal of the original Newell plane

Returns

a rotated polygon

Definition at line 445 of file AnalyticalGeometry.cpp.

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

References GeoLib::Polyline::addPoint(), 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 654 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.

Definition at line 654 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().

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.

{
    auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
    auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
    auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
    auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
    auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
 
    Eigen::Vector3d const pq = vq - vp;
    Eigen::Vector3d const pa = va - vp;
    Eigen::Vector3d const pb = vb - vp;
    Eigen::Vector3d const pc = vc - vp;
 
    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 MaterialPropertyLib::c, MathLib::TemplatePoint< T, DIM >::getCoords(), and MathLib::q.

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