diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 51e46fefd1777995f3d435a90c5a3cf76ce8c094..334434e499d844520a8dd93e88c5f010a9a3bb3a 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -234,93 +234,6 @@ bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
return false;
}
-bool isPointInTriangle(MathLib::Point3d const& p,
- MathLib::Point3d const& a,
- MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane,
- double eps_pnt_out_of_tri,
- GeoLib::TriangleTest algorithm)
-{
- switch (algorithm)
- {
- case GeoLib::GAUSS:
- return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
- case GeoLib::BARYCENTRIC:
- return barycentricPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
- default:
- ERR ("Selected algorithm for point in triangle testing not found, falling back on default.");
- }
- return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane, eps_pnt_out_of_tri);
-}
-
-bool gaussPointInTriangle(MathLib::Point3d const& q,
- MathLib::Point3d const& a,
- MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane,
- double eps_pnt_out_of_tri)
-{
- MathLib::Vector3 const v(a, b);
- MathLib::Vector3 const w(a, c);
-
- MathLib::DenseMatrix<double> mat (2,2);
- mat(0,0) = v.getSqrLength();
- mat(0,1) = v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
- mat(1,0) = mat(0,1);
- mat(1,1) = w.getSqrLength();
- double y[2] = {
- v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
- w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2])
- };
-
- MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
- gauss.solve(mat, y);
-
- const double lower (eps_pnt_out_of_tri);
- const double upper (1 + lower);
-
- if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper && y[0] + y[1] <=
- upper) {
- MathLib::Point3d const q_projected(std::array<double,3>{{
- a[0] + y[0] * v[0] + y[1] * w[0],
- a[1] + y[0] * v[1] + y[1] * w[1],
- a[2] + y[0] * v[2] + y[1] * w[2]
- }});
- if (MathLib::sqrDist(q, q_projected) <= eps_pnt_out_of_plane)
- return true;
- }
-
- return false;
-}
-
-bool barycentricPointInTriangle(MathLib::Point3d const& p,
- MathLib::Point3d const& a,
- MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane,
- double eps_pnt_out_of_tri)
-{
- if (std::abs(MathLib::orientation3d(p, a, b, c)) > eps_pnt_out_of_plane)
- return false;
-
- MathLib::Vector3 const pa (p,a);
- MathLib::Vector3 const pb (p,b);
- MathLib::Vector3 const pc (p,c);
- double const area_x_2 (calcTriangleArea(a,b,c) * 2);
-
- double const alpha ((MathLib::crossProduct(pb,pc).getLength()) / area_x_2);
- if (alpha < -eps_pnt_out_of_tri || alpha > 1+eps_pnt_out_of_tri)
- return false;
- double const beta ((MathLib::crossProduct(pc,pa).getLength()) / area_x_2);
- if (beta < -eps_pnt_out_of_tri || beta > 1+eps_pnt_out_of_tri)
- return false;
- double const gamma (1 - alpha - beta);
- if (gamma < -eps_pnt_out_of_tri || gamma > 1+eps_pnt_out_of_tri)
- return false;
- return true;
-}
-
void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat)
{
// *** some frequently used terms ***
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index a02c7cdd8e9f9f319a6a0c15d4075ce37931ec9e..83bac893dd2d963925ed8846d255d18fd6bcfed9 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -32,11 +32,6 @@ namespace GeoLib
class Polyline;
class LineSegment;
-enum TriangleTest
-{
- GAUSS, BARYCENTRIC
-};
-
enum Orientation
{
CW = 1, CCW = 2, COLLINEAR = 3
@@ -176,65 +171,6 @@ MathLib::DenseMatrix<double> rotatePointsToXY(InputIterator1 p_pnts_begin,
*/
void rotatePointsToXZ(std::vector<GeoLib::Point*> &pnts);
-/**
- * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
- * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
- * as well as an 'epsilon' outside the triangles plane.
- * @param p test point
- * @param a edge node of triangle
- * @param b edge node of triangle
- * @param c edge node of triangle
- * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
- * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
- * @param algorithm defines the method to use
- * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
- */
-bool isPointInTriangle(MathLib::Point3d const& p,
- MathLib::Point3d const& a,
- MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
- double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon(),
- GeoLib::TriangleTest algorithm = GeoLib::GAUSS);
-
-/**
- * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
- * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
- * as well as an 'epsilon' outside the triangles plane.
- * @param p test point
- * @param a edge node of triangle
- * @param b edge node of triangle
- * @param c edge node of triangle
- * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
- * ((orthogonal distance to the plane spaned by triangle)
- * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
- * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
- */
-bool gaussPointInTriangle(MathLib::Point3d const& p,
- MathLib::Point3d const& a, MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
- double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
-
-/**
- * Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.
- * Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle
- * as well as an 'epsilon' outside the triangles plane.
- * Algorithm based on "Fundamentals of Computer Graphics" by Peter Shirley.
- * @param p test point
- * @param a edge node of triangle
- * @param b edge node of triangle
- * @param c edge node of triangle
- * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc
- * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
- * @return true if the test point p is within the 'epsilon'-neighbourhood of the triangle
- */
-bool barycentricPointInTriangle(MathLib::Point3d const& p,
- MathLib::Point3d const& a, MathLib::Point3d const& b,
- MathLib::Point3d const& c,
- double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
- double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
-
/**
* test for intersections of the line segments of the Polyline
* @param ply the polyline
diff --git a/GeoLib/EarClippingTriangulation.cpp b/GeoLib/EarClippingTriangulation.cpp
index 523a684fa378f9d378fc9ccd9e901b5124907b43..6cbd549cafa166e2a8c41299d5461ef661f785c7 100644
--- a/GeoLib/EarClippingTriangulation.cpp
+++ b/GeoLib/EarClippingTriangulation.cpp
@@ -16,6 +16,8 @@
#include "BaseLib/uniqueInsert.h"
+#include "MathLib/GeometricBasics.h"
+
#include "Polygon.h"
#include "Triangle.h"
#include "Point.h"
@@ -115,7 +117,7 @@ bool EarClippingTriangulation::isEar(std::size_t v0, std::size_t v1, std::size_t
for (std::list<std::size_t>::const_iterator it (_vertex_list.begin ());
it != _vertex_list.end(); ++it) {
if (*it != v0 && *it != v1 && *it != v2) {
- if (isPointInTriangle (*_pnts[*it], *_pnts[v0], *_pnts[v1], *_pnts[v2])) {
+ if (MathLib::isPointInTriangle (*_pnts[*it], *_pnts[v0], *_pnts[v1], *_pnts[v2])) {
return false;
}
}
diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp
index 5480a9c28499593c3bfcf3ebc23b75826d885c3f..c28af316aaa0b09b8da8c2a1348a9e523de35271 100644
--- a/GeoLib/Triangle.cpp
+++ b/GeoLib/Triangle.cpp
@@ -14,6 +14,7 @@
#include "AnalyticalGeometry.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+#include "MathLib/GeometricBasics.h"
namespace GeoLib {
@@ -38,7 +39,7 @@ bool Triangle::containsPoint(MathLib::Point3d const& q, double eps) const
GeoLib::Point const& a(*(_pnts[_pnt_ids[0]]));
GeoLib::Point const& b(*(_pnts[_pnt_ids[1]]));
GeoLib::Point const& c(*(_pnts[_pnt_ids[2]]));
- return GeoLib::isPointInTriangle(q, a, b, c, eps);
+ return MathLib::isPointInTriangle(q, a, b, c, eps);
}
bool Triangle::containsPoint2D (Point const& pnt) const
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index 3d26f8e877f08863d8d200d1562be7bf475514d6..11e2fcf4e25cebb751a64cd231ab055b6336b30c 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -8,10 +8,15 @@
*
*/
-#include "GeometricBasics.h"
+#include <logog/include/logog.hpp>
+
+#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
#include "Point3d.h"
#include "Vector3.h"
+#include "GeometricBasics.h"
+
namespace MathLib
{
double orientation3d(MathLib::Point3d const& p,
@@ -75,4 +80,93 @@ bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a,
return false;
}
+bool isPointInTriangle(MathLib::Point3d const& p,
+ MathLib::Point3d const& a,
+ MathLib::Point3d const& b,
+ MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri,
+ MathLib::TriangleTest algorithm)
+{
+ switch (algorithm)
+ {
+ case MathLib::GAUSS:
+ return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
+ eps_pnt_out_of_tri);
+ case MathLib::BARYCENTRIC:
+ return barycentricPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
+ eps_pnt_out_of_tri);
+ default:
+ ERR("Selected algorithm for point in triangle testing not found, "
+ "falling back on default.");
+ }
+ return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
+ eps_pnt_out_of_tri);
+}
+
+bool gaussPointInTriangle(MathLib::Point3d const& q,
+ MathLib::Point3d const& a,
+ MathLib::Point3d const& b,
+ MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri)
+{
+ MathLib::Vector3 const v(a, b);
+ MathLib::Vector3 const w(a, c);
+
+ MathLib::DenseMatrix<double> mat(2, 2);
+ mat(0, 0) = v.getSqrLength();
+ mat(0, 1) = v[0] * w[0] + v[1] * w[1] + v[2] * w[2];
+ mat(1, 0) = mat(0, 1);
+ mat(1, 1) = w.getSqrLength();
+ double y[2] = {
+ v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
+ w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2])};
+
+ MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
+ gauss.solve(mat, y);
+
+ const double lower(eps_pnt_out_of_tri);
+ const double upper(1 + lower);
+
+ if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper &&
+ y[0] + y[1] <= upper)
+ {
+ MathLib::Point3d const q_projected(std::array<double, 3>{
+ {a[0] + y[0] * v[0] + y[1] * w[0], a[1] + y[0] * v[1] + y[1] * w[1],
+ a[2] + y[0] * v[2] + y[1] * w[2]}});
+ if (MathLib::sqrDist(q, q_projected) <= eps_pnt_out_of_plane)
+ return true;
+ }
+
+ return false;
+}
+
+bool barycentricPointInTriangle(MathLib::Point3d const& p,
+ MathLib::Point3d const& a,
+ MathLib::Point3d const& b,
+ MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane,
+ double eps_pnt_out_of_tri)
+{
+ if (std::abs(MathLib::orientation3d(p, a, b, c)) > eps_pnt_out_of_plane)
+ return false;
+
+ MathLib::Vector3 const pa(p, a);
+ MathLib::Vector3 const pb(p, b);
+ MathLib::Vector3 const pc(p, c);
+ double const area_x_2(calcTriangleArea(a, b, c) * 2);
+
+ double const alpha((MathLib::crossProduct(pb, pc).getLength()) / area_x_2);
+ if (alpha < -eps_pnt_out_of_tri || alpha > 1 + eps_pnt_out_of_tri)
+ return false;
+ double const beta((MathLib::crossProduct(pc, pa).getLength()) / area_x_2);
+ if (beta < -eps_pnt_out_of_tri || beta > 1 + eps_pnt_out_of_tri)
+ return false;
+ double const gamma(1 - alpha - beta);
+ if (gamma < -eps_pnt_out_of_tri || gamma > 1 + eps_pnt_out_of_tri)
+ return false;
+ return true;
+}
+
} // end namespace MathLib
diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h
index ba569dddbf31abce433629e5c7399c61d2db1828..d1c45460ad380e541fc2eda10816b137d77de7ea 100644
--- a/MathLib/GeometricBasics.h
+++ b/MathLib/GeometricBasics.h
@@ -20,6 +20,11 @@ namespace MathLib
template <typename T, std::size_t DIM> class TemplatePoint;
typedef MathLib::TemplatePoint<double,3> Point3d;
+enum TriangleTest
+{
+ GAUSS, BARYCENTRIC
+};
+
/**
* Checks if a point p is on the left or right side of a plane spanned by three
* points a, b, c.
@@ -76,6 +81,79 @@ bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a,
MathLib::Point3d const& d,
double eps = std::numeric_limits<double>::epsilon());
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes
+ * a, b and c.Using two eps-values it is possible to test an 'epsilon'
+ * neighbourhood around the triangle as well as an 'epsilon' outside the
+ * triangles plane.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane
+ * spanned by abc
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of
+ * abc
+ * @param algorithm defines the method to use
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the
+ * triangle
+ */
+bool isPointInTriangle(
+ MathLib::Point3d const& p,
+ MathLib::Point3d const& a,
+ MathLib::Point3d const& b,
+ MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon(),
+ MathLib::TriangleTest algorithm = MathLib::GAUSS);
+
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes
+ * a, b and c.
+ * Using two eps-values it is possible to test an 'epsilon' neighbourhood around
+ * the triangle
+ * as well as an 'epsilon' outside the triangles plane.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane
+ * spanned by abc ((orthogonal distance to the plane spaned by triangle)
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of
+ * abc
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the
+ * triangle
+ */
+bool gaussPointInTriangle(
+ MathLib::Point3d const& p, MathLib::Point3d const& a,
+ MathLib::Point3d const& b, MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
+
+/**
+ * Tests if the given point p is within the triangle, defined by its edge nodes
+ * a, b and c.
+ * Using two eps-values it is possible to test an 'epsilon' neighbourhood around
+ * the triangle
+ * as well as an 'epsilon' outside the triangles plane.
+ * Algorithm based on "Fundamentals of Computer Graphics" by Peter Shirley.
+ * @param p test point
+ * @param a edge node of triangle
+ * @param b edge node of triangle
+ * @param c edge node of triangle
+ * @param eps_pnt_out_of_plane eps allowing for p to be slightly off the plane
+ * spanned by abc
+ * @param eps_pnt_out_of_tri eps allowing for p to be slightly off outside of
+ * abc
+ * @return true if the test point p is within the 'epsilon'-neighbourhood of the
+ * triangle
+ */
+bool barycentricPointInTriangle(
+ MathLib::Point3d const& p, MathLib::Point3d const& a,
+ MathLib::Point3d const& b, MathLib::Point3d const& c,
+ double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
+ double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
+
} // end namespace MathLib
#endif
diff --git a/MeshLib/Elements/QuadRule4.cpp b/MeshLib/Elements/QuadRule4.cpp
index fec6900320583cb0297e6e37b226278b0d8e1450..d370adc593c608e819105708c6644428225400e0 100644
--- a/MeshLib/Elements/QuadRule4.cpp
+++ b/MeshLib/Elements/QuadRule4.cpp
@@ -12,6 +12,7 @@
#include "logog/include/logog.hpp"
#include "GeoLib/AnalyticalGeometry.h"
+
#include "MathLib/GeometricBasics.h"
#include "MeshLib/Node.h"
@@ -34,8 +35,8 @@ double QuadRule4::computeVolume(Node const* const* _nodes)
bool QuadRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
{
- return (GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps) ||
- GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
+ return (MathLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps) ||
+ MathLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
}
unsigned QuadRule4::identifyFace(Node const* const* _nodes, Node* nodes[3])
diff --git a/MeshLib/Elements/TriRule3.cpp b/MeshLib/Elements/TriRule3.cpp
index 0927588134cac1a6860d07766bff4918e184bfab..07d4fed4383793037271ea167b72ef7e8f7c4cf6 100644
--- a/MeshLib/Elements/TriRule3.cpp
+++ b/MeshLib/Elements/TriRule3.cpp
@@ -11,7 +11,6 @@
#include "logog/include/logog.hpp"
-#include "GeoLib/AnalyticalGeometry.h"
#include "MathLib/GeometricBasics.h"
#include "MeshLib/Node.h"
@@ -32,7 +31,7 @@ double TriRule3::computeVolume(Node const* const* _nodes)
bool TriRule3::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps)
{
- return GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps);
+ return MathLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps);
}
unsigned TriRule3::identifyFace(Node const* const* _nodes, Node* nodes[3])