diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index 9b02181d3934986562479e783bbb942bdb1d1919..b8006bd6c43f744aa8af341580431739bbe5b913 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -321,36 +321,6 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p, return true; } -bool isPointInTetrahedron(MathLib::Point3d const& p, - MathLib::Point3d const& a, MathLib::Point3d const& b, - MathLib::Point3d const& c, MathLib::Point3d const& d, double eps) -{ - double const d0 (MathLib::orientation3d(d,a,b,c)); - // if tetrahedron is not coplanar - if (std::abs(d0) > std::numeric_limits<double>::epsilon()) - { - bool const d0_sign (d0>0); - // if p is on the same side of bcd as a - double const d1 (MathLib::orientation3d(d, p, b, c)); - if (!(d0_sign == (d1>=0) || std::abs(d1) < eps)) - return false; - // if p is on the same side of acd as b - double const d2 (MathLib::orientation3d(d, a, p, c)); - if (!(d0_sign == (d2>=0) || std::abs(d2) < eps)) - return false; - // if p is on the same side of abd as c - double const d3 (MathLib::orientation3d(d, a, b, p)); - if (!(d0_sign == (d3>=0) || std::abs(d3) < eps)) - return false; - // if p is on the same side of abc as d - double const d4 (MathLib::orientation3d(p, a, b, c)); - if (!(d0_sign == (d4>=0) || std::abs(d4) < eps)) - return false; - return true; - } - return false; -} - double calcTriangleArea(MathLib::Point3d const& a, MathLib::Point3d const& b, MathLib::Point3d const& c) { diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h index 709cdadfd6c075aa173fd39caa0833242723d624..cf5f1d83c5a974578792c1b3302a16b85853275a 100644 --- a/GeoLib/AnalyticalGeometry.h +++ b/GeoLib/AnalyticalGeometry.h @@ -243,23 +243,6 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p, 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 located within a tetrahedron spanned by points a, b, c, d. - * If the tet specified by a, b, c, d is degenerated (i.e. all points are coplanar) the function - * will return false because there is no meaningful concept of "inside" for such elements. - * @param p test point - * @param a edge node of tetrahedron - * @param b edge node of tetrahedron - * @param c edge node of tetrahedron - * @param d edge node of tetrahedron - * @param eps Accounts for numerical inaccuracies by allowing a point to be slightly outside of the element and still be regarded as inside. - * @return true if the test point p is not located outside of abcd (i.e. inside or on a plane/edge). - */ -bool isPointInTetrahedron(MathLib::Point3d const& p, - MathLib::Point3d const& a, MathLib::Point3d const& b, - MathLib::Point3d const& c, MathLib::Point3d const& d, - double eps = std::numeric_limits<double>::epsilon()); - /** * test for intersections of the line segments of the Polyline * @param ply the polyline diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp index 5e43c3cc61f204e093e6066827fe3c9e7f1f710d..caab906d2d31b68da916f1d376d3763a2fd05d53 100644 --- a/MathLib/GeometricBasics.cpp +++ b/MathLib/GeometricBasics.cpp @@ -36,4 +36,34 @@ double calcTetrahedronVolume(MathLib::Point3d const& a, return std::abs(MathLib::scalarTriple(ac, ad, ab)) / 6.0; } +bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a, + MathLib::Point3d const& b, MathLib::Point3d const& c, + MathLib::Point3d const& d, double eps) +{ + double const d0 (MathLib::orientation3d(d,a,b,c)); + // if tetrahedron is not coplanar + if (std::abs(d0) > std::numeric_limits<double>::epsilon()) + { + bool const d0_sign (d0>0); + // if p is on the same side of bcd as a + double const d1 (MathLib::orientation3d(d, p, b, c)); + if (!(d0_sign == (d1>=0) || std::abs(d1) < eps)) + return false; + // if p is on the same side of acd as b + double const d2 (MathLib::orientation3d(d, a, p, c)); + if (!(d0_sign == (d2>=0) || std::abs(d2) < eps)) + return false; + // if p is on the same side of abd as c + double const d3 (MathLib::orientation3d(d, a, b, p)); + if (!(d0_sign == (d3>=0) || std::abs(d3) < eps)) + return false; + // if p is on the same side of abc as d + double const d4 (MathLib::orientation3d(p, a, b, c)); + if (!(d0_sign == (d4>=0) || std::abs(d4) < eps)) + return false; + return true; + } + return false; +} + } // end namespace MathLib diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h index d5164d5365c58b55637dde8963c6c5e454c1e15a..e97a9b5e407b2ba461c298cbb9c9084994676bb5 100644 --- a/MathLib/GeometricBasics.h +++ b/MathLib/GeometricBasics.h @@ -12,6 +12,7 @@ #define GEOMETRIC_BASICS_H_ #include <cstddef> +#include <limits> namespace MathLib { @@ -44,6 +45,27 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1, MathLib::Point3d const& x2, MathLib::Point3d const& x3, MathLib::Point3d const& x4); +/** + * Tests if the given point p is located within a tetrahedron spanned by points + * a, b, c, d. + * If the tet specified by a, b, c, d is degenerated (i.e. all points are + * coplanar) the function + * will return false because there is no meaningful concept of "inside" for such + * elements. + * @param p test point + * @param a edge node of tetrahedron + * @param b edge node of tetrahedron + * @param c edge node of tetrahedron + * @param d edge node of tetrahedron + * @param eps Accounts for numerical inaccuracies by allowing a point to be + * slightly outside of the element and still be regarded as inside. + * @return true if the test point p is not located outside of abcd (i.e. inside + * or on a plane/edge). + */ +bool isPointInTetrahedron(MathLib::Point3d const& p, MathLib::Point3d const& a, + MathLib::Point3d const& b, MathLib::Point3d const& c, + MathLib::Point3d const& d, + double eps = std::numeric_limits<double>::epsilon()); } // end namespace MathLib diff --git a/MeshLib/Elements/HexRule8.cpp b/MeshLib/Elements/HexRule8.cpp index db669c01b53c10ac0ea4482d9d063b97db6e9a7f..8a4a9e4ce0d0e2cab4a1c415310331f28c98e249 100644 --- a/MeshLib/Elements/HexRule8.cpp +++ b/MeshLib/Elements/HexRule8.cpp @@ -13,7 +13,6 @@ #include "logog/include/logog.hpp" -#include "GeoLib/AnalyticalGeometry.h" #include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" @@ -73,12 +72,12 @@ double HexRule8::computeVolume(Node const* const* _nodes) bool HexRule8::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) { - return (GeoLib::isPointInTetrahedron(pnt, *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2], eps)); + return (MathLib::isPointInTetrahedron(pnt, *_nodes[4], *_nodes[7], *_nodes[5], *_nodes[0], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[3], *_nodes[1], *_nodes[0], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[3], *_nodes[0], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[5], *_nodes[7], *_nodes[6], *_nodes[2], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[3], *_nodes[5], *_nodes[2], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[3], *_nodes[7], *_nodes[5], *_nodes[2], eps)); } unsigned HexRule8::identifyFace(Node const* const* _nodes, Node* nodes[3]) diff --git a/MeshLib/Elements/PrismRule6.cpp b/MeshLib/Elements/PrismRule6.cpp index 8a1f72546f5d29e870035c50e69bb5af8802c8d2..33dcd3688605eac7e5430f7665f0fc048f6370e4 100644 --- a/MeshLib/Elements/PrismRule6.cpp +++ b/MeshLib/Elements/PrismRule6.cpp @@ -11,7 +11,6 @@ #include "logog/include/logog.hpp" -#include "GeoLib/AnalyticalGeometry.h" #include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" @@ -79,9 +78,9 @@ double PrismRule6::computeVolume(Node const* const* _nodes) bool PrismRule6::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) { - return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3], eps)); + return (MathLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[1], *_nodes[4], *_nodes[2], *_nodes[3], eps) || + MathLib::isPointInTetrahedron(pnt, *_nodes[2], *_nodes[4], *_nodes[5], *_nodes[3], eps)); } unsigned PrismRule6::identifyFace(Node const* const* _nodes, Node* nodes[3]) diff --git a/MeshLib/Elements/PyramidRule5.cpp b/MeshLib/Elements/PyramidRule5.cpp index 98b3edcfdace8225dc088e41996b199976dcac29..476aad78092d513778b38857b04095c8435518ef 100644 --- a/MeshLib/Elements/PyramidRule5.cpp +++ b/MeshLib/Elements/PyramidRule5.cpp @@ -11,7 +11,6 @@ #include "logog/include/logog.hpp" -#include "GeoLib/AnalyticalGeometry.h" #include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" @@ -75,10 +74,14 @@ double PyramidRule5::computeVolume(Node const* const* _nodes) + MathLib::calcTetrahedronVolume(*_nodes[2], *_nodes[3], *_nodes[0], *_nodes[4]); } -bool PyramidRule5::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) +bool PyramidRule5::isPntInElement(Node const* const* _nodes, + MathLib::Point3d const& pnt, + double eps) { - return (GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4], eps) || - GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[2], *_nodes[3], *_nodes[4], eps)); + return (MathLib::isPointInTetrahedron( + pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[4], eps) || + MathLib::isPointInTetrahedron( + pnt, *_nodes[0], *_nodes[2], *_nodes[3], *_nodes[4], eps)); } unsigned PyramidRule5::identifyFace(Node const* const* _nodes, Node* nodes[3]) diff --git a/MeshLib/Elements/TetRule4.cpp b/MeshLib/Elements/TetRule4.cpp index daeee3683381dd24d33567be34a4b91cb56f7bba..7f545b5454f1823d62fa432e81613d98e0bbbf27 100644 --- a/MeshLib/Elements/TetRule4.cpp +++ b/MeshLib/Elements/TetRule4.cpp @@ -13,7 +13,6 @@ #include "logog/include/logog.hpp" -#include "GeoLib/AnalyticalGeometry.h" #include "MathLib/GeometricBasics.h" #include "MeshLib/Node.h" @@ -60,7 +59,7 @@ double TetRule4::computeVolume(Node const* const* _nodes) bool TetRule4::isPntInElement(Node const* const* _nodes, MathLib::Point3d const& pnt, double eps) { - return GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps); + return MathLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps); } unsigned TetRule4::identifyFace(Node const* const* _nodes, Node* nodes[3])