diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 7e5773a2767246415ac21a515c32e058eeaf0873..c9d9a01740cf7e77f5cd071da670f006e0c01803 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -242,6 +242,35 @@ bool isPointInTriangle(GeoLib::Point const& q,
return false;
}
+bool isPointInTetrahedron(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
+ GeoLib::Point const& c, GeoLib::Point const& d, double eps)
+{
+ double const d0 (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 (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 (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 (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 (orientation3d(p, a, b, c));
+ if (!(d0_sign == (d4>=0) || std::abs(d4) < eps))
+ return false;
+ return true;
+ }
+ return false;
+}
+
double calcTriangleArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
{
MathLib::Vector3 const u(a,c);
@@ -405,6 +434,15 @@ double scalarTriple(MathLib::Vector3 const& u, MathLib::Vector3 const& v, MathLi
return MathLib::scalarProduct(cross,w);
}
+double orientation3d(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
+{
+ MathLib::Vector3 const ap (a, p);
+ MathLib::Vector3 const bp (b, p);
+ MathLib::Vector3 const cp (c, p);
+ return scalarTriple(bp,cp,ap);
+}
+
bool dividedByPlane(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Point& c, const GeoLib::Point& d)
{
for (unsigned x=0; x<3; ++x)
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 65a241a7c30b31c64bf0efa2b606cc426dc8e83a..1fdc64904d1cfac8f8bb68dea0c7f06e7746fa5a 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -129,6 +129,21 @@ bool isPointInTriangle(GeoLib::Point const& p,
GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c,
double eps = 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(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
+ GeoLib::Point const& c, GeoLib::Point const& d, double eps = std::numeric_limits<double>::epsilon());
+
/**
* test for intersections of the line segments of the Polyline
* @param ply the polyline
@@ -173,7 +188,19 @@ GeoLib::Point* triangleLineIntersection(GeoLib::Point const& a, GeoLib::Point co
/// Calculates the scalar triple (u x v) . w
double scalarTriple(MathLib::Vector3 const& u, MathLib::Vector3 const& v, MathLib::Vector3 const& w);
-/**
+/**
+ * Checks if a point p is on the left or right side of a plane spanned by three points a, b, c.
+ * @param p point to test
+ * @param a first point on plane
+ * @param b second point on plane
+ * @param c third point on plane
+ * @return If the triangle abc is ordered counterclockwise when viewed from p, the method will return a negative value,
+ * otherwise it will return a positive value. If p is coplanar with abc, the function will return 0.
+ */
+double orientation3d(GeoLib::Point const& p,
+ GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c);
+
+/**
* Checks if a and b can be placed on a plane such that c and d lie on different sides of said plane.
* (In 2D space this checks if c and d are on different sides of a line through a and b.)
* @param a first point on plane
@@ -188,8 +215,7 @@ double scalarTriple(MathLib::Vector3 const& u, MathLib::Vector3 const& v, MathLi
/// Checks if the four given points are located on a plane.
bool isCoplanar(const GeoLib::Point& a, const GeoLib::Point& b,
const GeoLib::Point& c, const GeoLib::Point& d);
-
-
+
/**
* Method first computes the intersection points of line segements of GeoLib::Polyline objects
* (@see computeIntersectionPoints()) and pushes each intersection point in the GeoLib::PointVec
diff --git a/MeshLib/Elements/Element.cpp b/MeshLib/Elements/Element.cpp
index f63bbf5139db55ed782c0399f507bb5e3ac78017..276c19d8d85c0718abb15926ea28cdfe3f555257 100644
--- a/MeshLib/Elements/Element.cpp
+++ b/MeshLib/Elements/Element.cpp
@@ -170,4 +170,3 @@ bool Element::isBoundaryElement() const
}
}
-
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 6e2dd8003768da3f312185cd1e853ca5e7041db1..23435d802d59eb384ef047829c3d04afb81c4357 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -19,10 +19,13 @@
#include <limits>
#include <boost/optional.hpp>
+#include "GeoLib/Point.h"
+
#include "MeshLib/MeshEnums.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshQuality/ElementErrorCode.h"
+
namespace MeshLib {
class Node;
@@ -154,6 +157,14 @@ public:
/// Returns true if these two indeces form an edge and false otherwise
virtual bool isEdge(unsigned i, unsigned j) const = 0;
+ /**
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is not outside the element, false otherwise
+ */
+ virtual bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const = 0;
+
/**
* Tests if the element is geometrically valid.
*/
diff --git a/MeshLib/Elements/Face.h b/MeshLib/Elements/Face.h
index 6c767316171b038324f3d034ff2a6c119ed9576b..d9a5ad0babfb597da6c7960fb4adf0b36517229c 100644
--- a/MeshLib/Elements/Face.h
+++ b/MeshLib/Elements/Face.h
@@ -59,14 +59,6 @@ public:
/// Destructor
virtual ~Face();
- /**
- * Check if the 3d GeoLib::Point is inside of the element.
- * @param pnt the 3d GeoLib::Point object
- * @param eps tolerance for numerical algorithm used or computing the property
- * @return true if the point is inside the element, false otherwise
- */
- virtual bool isPntInside(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const = 0;
-
/**
* This method is pure virtual and is inherited from class @sa Element.
* It has to be implemented in the derived classes of class Face!
diff --git a/MeshLib/Elements/TemplateHex-impl.h b/MeshLib/Elements/TemplateHex-impl.h
index 3469c18b60ad7c4b98066b1e6d498d3e76dc9513..7a462b3357980099c836350cebbd780b19a9decd 100644
--- a/MeshLib/Elements/TemplateHex-impl.h
+++ b/MeshLib/Elements/TemplateHex-impl.h
@@ -139,6 +139,17 @@ bool TemplateHex<NNODES,CELLHEXTYPE>::isEdge(unsigned idx1, unsigned idx2) const
return false;
}
+template <unsigned NNODES, CellType CELLHEXTYPE>
+bool TemplateHex<NNODES,CELLHEXTYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
+{
+ 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));
+}
+
template <unsigned NNODES, CellType CELLHEXTYPE>
Element* TemplateHex<NNODES,CELLHEXTYPE>::clone() const
{
diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h
index 0fcb9bf951c8bfbc17ad95da2772c35510e99831..43231955e5ccd4d44aa4c7f1ecb47733cc6497f3 100644
--- a/MeshLib/Elements/TemplateHex.h
+++ b/MeshLib/Elements/TemplateHex.h
@@ -104,6 +104,14 @@ public:
/// Returns true if these two indices form an edge and false otherwise
bool isEdge(unsigned i, unsigned j) const;
+ /**
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is not outside the element, false otherwise
+ */
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+
/**
* Tests if the element is geometrically valid.
* @param check_zero_volume indicates if volume == 0 should be checked
diff --git a/MeshLib/Elements/TemplateLine-impl.h b/MeshLib/Elements/TemplateLine-impl.h
index 875cb8ecf29a59a215a515d04bb970ec04f899b7..3c3a1499c666c1319e3de5ebf296941c2efad36b 100644
--- a/MeshLib/Elements/TemplateLine-impl.h
+++ b/MeshLib/Elements/TemplateLine-impl.h
@@ -55,6 +55,15 @@ template <unsigned NNODES, CellType CELLLINETYPE>
TemplateLine<NNODES,CELLLINETYPE>::~TemplateLine()
{}
+template <unsigned NNODES, CellType CELLLINETYPE>
+bool TemplateLine<NNODES,CELLLINETYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
+{
+ double tmp;
+ double tmp_dst(0);
+ double const dist =MathLib::calcProjPntToLineAndDists(pnt.getCoords(), _nodes[0]->getCoords(), _nodes[1]->getCoords(), tmp, tmp_dst);
+ return (dist < eps);
+}
+
template <unsigned NNODES, CellType CELLLINETYPE>
ElementErrorCode TemplateLine<NNODES,CELLLINETYPE>::validate() const
{
diff --git a/MeshLib/Elements/TemplateLine.h b/MeshLib/Elements/TemplateLine.h
index b3c5bb29e246bdacf2ae45d34411f5d0b82b6971..3a3d661baea87eedf4d177037a983a9b2bbc0ca6 100644
--- a/MeshLib/Elements/TemplateLine.h
+++ b/MeshLib/Elements/TemplateLine.h
@@ -115,6 +115,14 @@ public:
return false;
}
+ /**
+ * Checks if a point is located on the line
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is located on the line, false otherwise
+ */
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+
/**
* Tests if the element is geometrically valid.
* @param check_zero_volume indicates if area == 0 should be checked
diff --git a/MeshLib/Elements/TemplatePrism-impl.h b/MeshLib/Elements/TemplatePrism-impl.h
index b6c00bf195b6ff47dca1fac979a14ed33e128080..b460682cd20d5014e2498d874771521e38b8046a 100644
--- a/MeshLib/Elements/TemplatePrism-impl.h
+++ b/MeshLib/Elements/TemplatePrism-impl.h
@@ -148,6 +148,14 @@ bool TemplatePrism<NNODES,CELLPRISMTYPE>::isEdge(unsigned idx1, unsigned idx2) c
return false;
}
+template <unsigned NNODES, CellType CELLPRISMTYPE>
+bool TemplatePrism<NNODES,CELLPRISMTYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
+{
+ 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));
+}
+
template <unsigned NNODES, CellType CELLPRISMTYPE>
Element* TemplatePrism<NNODES,CELLPRISMTYPE>::clone() const
{
diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h
index d91f295402421c6098b989d12c33028fdf1a2c59..ea98f9ce71fb009c14a54cf665b834bdc5bfeedb 100644
--- a/MeshLib/Elements/TemplatePrism.h
+++ b/MeshLib/Elements/TemplatePrism.h
@@ -102,6 +102,14 @@ public:
/// Returns true if these two indeces form an edge and false otherwise
bool isEdge(unsigned i, unsigned j) const;
+ /**
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is not outside the element, false otherwise
+ */
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+
/**
* Tests if the element is geometrically valid.
* @param check_zero_volume indicates if volume == 0 should be checked
diff --git a/MeshLib/Elements/TemplatePyramid-impl.h b/MeshLib/Elements/TemplatePyramid-impl.h
index d21b6b0031a5513e7078d582684acc5de8b4b7df..a1a3f6ee5413daf7e5dccb57056aad7dd6150109 100644
--- a/MeshLib/Elements/TemplatePyramid-impl.h
+++ b/MeshLib/Elements/TemplatePyramid-impl.h
@@ -150,6 +150,13 @@ bool TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::isEdge(unsigned idx1, unsigned idx
return false;
}
+template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
+bool TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
+{
+ 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));
+}
+
template <unsigned NNODES, CellType CELLPYRAMIDTYPE>
Element* TemplatePyramid<NNODES,CELLPYRAMIDTYPE>::clone() const
{
diff --git a/MeshLib/Elements/TemplatePyramid.h b/MeshLib/Elements/TemplatePyramid.h
index da8b34f776ea6e5fc6d8080e2c41c58521ed3de2..620075d88c404f532f2fc839c4965b50752da814 100644
--- a/MeshLib/Elements/TemplatePyramid.h
+++ b/MeshLib/Elements/TemplatePyramid.h
@@ -100,6 +100,14 @@ public:
/// Returns true if these two indeces form an edge and false otherwise
bool isEdge(unsigned i, unsigned j) const;
+ /**
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is not outside the element, false otherwise
+ */
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+
/**
* Tests if the element is geometrically valid.
* @param check_zero_volume indicates if volume == 0 should be checked
diff --git a/MeshLib/Elements/TemplateQuad-impl.h b/MeshLib/Elements/TemplateQuad-impl.h
index 5c6e7bb4d4a829cb7b76284e9926fbe288dd6aff..0c63f01fd3573fba812ba8de1437a29067feb8ce 100644
--- a/MeshLib/Elements/TemplateQuad-impl.h
+++ b/MeshLib/Elements/TemplateQuad-impl.h
@@ -106,10 +106,12 @@ bool TemplateQuad<NNODES,CELLQUADTYPE>::isEdge(unsigned idx1, unsigned idx2) con
}
template <unsigned NNODES, CellType CELLQUADTYPE>
-bool TemplateQuad<NNODES,CELLQUADTYPE>::isPntInside(GeoLib::Point const& pnt, double eps) const
+bool TemplateQuad<NNODES,CELLQUADTYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
{
+ bool a (GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps));
+ bool b ( GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
return (GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps) ||
- GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
+ GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[2], *_nodes[3], eps));
}
template <unsigned NNODES, CellType CELLQUADTYPE>
diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h
index 21c4fe4ed3e83dd42330fe245fa1bf0305a6a982..407c2b98056e693b2319d30436cd30b69e3d3cf4 100644
--- a/MeshLib/Elements/TemplateQuad.h
+++ b/MeshLib/Elements/TemplateQuad.h
@@ -90,12 +90,12 @@ public:
bool isEdge(unsigned i, unsigned j) const;
/**
- * Check if the 3d GeoLib::Point is inside of the quad element.
- * @param pnt the 3d GeoLib::Point object
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
* @param eps tolerance for numerical algorithm used or computing the property
- * @return true if the point is inside the element, false otherwise
+ * @return true if the point is not outside the element, false otherwise
*/
- virtual bool isPntInside(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
/**
* Tests if the element is geometrically valid.
diff --git a/MeshLib/Elements/TemplateTet-impl.h b/MeshLib/Elements/TemplateTet-impl.h
index 601f5199bcffa76a395aa59fb7e404dfed971bce..5d70c1778a6e9756ce27911acfc46312051bf594 100644
--- a/MeshLib/Elements/TemplateTet-impl.h
+++ b/MeshLib/Elements/TemplateTet-impl.h
@@ -128,6 +128,12 @@ bool TemplateTet<NNODES,CELLTETTYPE>::isEdge(unsigned idx1, unsigned idx2) const
return false;
}
+template <unsigned NNODES, CellType CELLTETTYPE>
+bool TemplateTet<NNODES,CELLTETTYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
+{
+ return GeoLib::isPointInTetrahedron(pnt, *_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3], eps);
+}
+
template <unsigned NNODES, CellType CELLTETTYPE>
Element* TemplateTet<NNODES,CELLTETTYPE>::clone() const
{
diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h
index b98c1e795802fbc65050cd3f573bdab93bc99cf1..7fee311b3484c277eaa9eb096630529e7cfb7d23 100644
--- a/MeshLib/Elements/TemplateTet.h
+++ b/MeshLib/Elements/TemplateTet.h
@@ -99,6 +99,14 @@ public:
/// Returns true if these two indeces form an edge and false otherwise
bool isEdge(unsigned i, unsigned j) const;
+ /**
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
+ * @param eps tolerance for numerical algorithm used or computing the property
+ * @return true if the point is not outside the element, false otherwise
+ */
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+
/**
* Tests if the element is geometrically valid.
* @param check_zero_volume indicates if volume == 0 should be checked
diff --git a/MeshLib/Elements/TemplateTri-impl.h b/MeshLib/Elements/TemplateTri-impl.h
index a4d2e0b735e29a4dc51fdb6cda2e83549ed7557a..04f2f581e27bd33c28c2bc6e62efeedce8a163db 100644
--- a/MeshLib/Elements/TemplateTri-impl.h
+++ b/MeshLib/Elements/TemplateTri-impl.h
@@ -88,7 +88,7 @@ bool TemplateTri<NNODES,CELLTRITYPE>::isEdge(unsigned idx1, unsigned idx2) const
}
template <unsigned NNODES, CellType CELLTRITYPE>
-bool TemplateTri<NNODES,CELLTRITYPE>::isPntInside(GeoLib::Point const& pnt, double eps) const
+bool TemplateTri<NNODES,CELLTRITYPE>::isPntInElement(GeoLib::Point const& pnt, double eps) const
{
return GeoLib::isPointInTriangle(pnt, *_nodes[0], *_nodes[1], *_nodes[2], eps);
}
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
index 465a8daed055a51d38ec5f365a8f126f5361db28..0ddaf4983489263af8efce6e282a867376b0413c 100644
--- a/MeshLib/Elements/TemplateTri.h
+++ b/MeshLib/Elements/TemplateTri.h
@@ -94,12 +94,12 @@ public:
bool isEdge(unsigned idx1, unsigned idx2) const;
/**
- * Check if the 3d GeoLib::Point is inside of the element.
- * @param pnt the 3d GeoLib::Point object
+ * Checks if a point is inside the element.
+ * @param pnt a 3D GeoLib::Point object
* @param eps tolerance for numerical algorithm used or computing the property
- * @return true if the point is inside the element, false otherwise
+ * @return true if the point is not outside the element, false otherwise
*/
- virtual bool isPntInside(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
+ bool isPntInElement(GeoLib::Point const& pnt, double eps = std::numeric_limits<double>::epsilon()) const;
/**
* Tests if the element is geometrically valid
diff --git a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
index 281829ac3fc8b0c19f58f7672eee8dcca305cda1..4f9ffafe00a60421aafcbb90b61e46805956b8ae 100644
--- a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
+++ b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
@@ -94,7 +94,7 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh(Mesh *dest_mes
for (size_t j(0); j<n_nodes_in_vec; j++) {
MeshLib::Node const*const j_th_node((*i_th_vec)[j]);
if (elem_aabb.containsPoint(*j_th_node)) {
- if (dynamic_cast<MeshLib::Face*>(dest_elements[k])->isPntInside(* dynamic_cast<GeoLib::Point const*>(j_th_node), 30)) {
+ if (dest_elements[k]->isPntInElement(*j_th_node, 30)) {
dest_properties[k] += interpolated_src_node_properties[(*i_th_vec)[j]->getID()];
cnt++;
}
diff --git a/Tests/MeshLib/TestPntInElement.cpp b/Tests/MeshLib/TestPntInElement.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb94fc0ec6bb65560dc55bb005f380422b403748
--- /dev/null
+++ b/Tests/MeshLib/TestPntInElement.cpp
@@ -0,0 +1,157 @@
+/**
+ * @file TestPntInElement.cpp
+ * @author Karsten Rink
+ * @date 2014-09-23
+ * @brief Tests for check if a point is located inside an element
+ *
+ * @copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include "gtest/gtest.h"
+
+#include "Mesh.h"
+#include "Node.h"
+#include "MeshEditing/MeshRevision.h"
+#include "Elements/Element.h"
+#include "Elements/Tri.h"
+#include "Elements/Quad.h"
+#include "Elements/Tet.h"
+#include "Elements/Hex.h"
+#include "Elements/Pyramid.h"
+#include "Elements/Prism.h"
+
+
+std::vector<MeshLib::Node*> createNodes()
+{
+ std::vector<MeshLib::Node*> nodes;
+ nodes.push_back(new MeshLib::Node(0,0,0));
+ nodes.push_back(new MeshLib::Node(1,0,0));
+ nodes.push_back(new MeshLib::Node(1,1,0));
+ nodes.push_back(new MeshLib::Node(0,1,0));
+ nodes.push_back(new MeshLib::Node(0,0,1));
+ nodes.push_back(new MeshLib::Node(1,0,1));
+ nodes.push_back(new MeshLib::Node(1,1,1));
+ nodes.push_back(new MeshLib::Node(0,1,1));
+ return nodes;
+}
+
+void deleteNodes(std::vector<MeshLib::Node*> &nodes)
+{
+ std::for_each(nodes.begin(), nodes.end(), [](MeshLib::Node* node){ delete node; });
+}
+
+TEST(IsPntInElement, Line)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 2> line_nodes = { nodes[0], nodes[4] };
+ MeshLib::Line line(line_nodes);
+ pnt = GeoLib::Point(0,0,0.7);
+ ASSERT_EQ(true, line.isPntInElement(pnt));
+ pnt = GeoLib::Point(0,0.1,0.7);
+ ASSERT_EQ(false, line.isPntInElement(pnt));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Tri)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 3> tri_nodes = { nodes[0], nodes[1], nodes[4] };
+ MeshLib::Tri tri(tri_nodes);
+
+ pnt = GeoLib::Point(0.1,0,0.1);
+ ASSERT_EQ(true, tri.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.9,0,0.7);
+ ASSERT_EQ(false, tri.isPntInElement(pnt));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Quad)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 4> quad_nodes = { nodes[0], nodes[1], nodes[5], nodes[4] };
+ MeshLib::Quad quad(quad_nodes);
+
+ pnt = GeoLib::Point(0.1,0,0.1);
+ ASSERT_EQ(true, quad.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.999,0,0.001);
+ ASSERT_EQ(true, quad.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.5,0.00001,1);
+ ASSERT_EQ(false, quad.isPntInElement(pnt));
+ ASSERT_EQ(true, quad.isPntInElement(pnt, 0.001));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Tet)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 4> tet_nodes = { nodes[0], nodes[1], nodes[2], nodes[4] };
+ MeshLib::Tet tet(tet_nodes);
+
+ pnt = GeoLib::Point(0.5,0.3,0.1);
+ ASSERT_EQ(true, tet.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.5,0.6,0.1);
+ ASSERT_EQ(false, tet.isPntInElement(pnt));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Pyramid)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 5> pyr_nodes;
+ std::copy(nodes.begin(), nodes.begin()+5, pyr_nodes.begin());
+ MeshLib::Pyramid pyr(pyr_nodes);
+
+ pnt = GeoLib::Point(0.5,0.00001,-0.000001);
+ ASSERT_EQ(false, pyr.isPntInElement(pnt));
+ ASSERT_EQ(true, pyr.isPntInElement(pnt, 0.0001));
+ pnt = GeoLib::Point(0.5,0.5,0.1);
+ ASSERT_EQ(true, pyr.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.5,0.5,0.51);
+ ASSERT_EQ(false, pyr.isPntInElement(pnt));
+ ASSERT_EQ(true, pyr.isPntInElement(pnt, 0.02));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Prism)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 6> prism_nodes = { nodes[0], nodes[1], nodes[2], nodes[4], nodes[5], nodes[6] };
+ MeshLib::Prism prism(prism_nodes);
+
+ pnt = GeoLib::Point(0.5,0.5,0.1);
+ ASSERT_EQ(true, prism.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.49,0.51,0.1);
+ ASSERT_EQ(false, prism.isPntInElement(pnt));
+ ASSERT_EQ(true, prism.isPntInElement(pnt, 0.03));
+ deleteNodes(nodes);
+}
+
+TEST(IsPntInElement, Hex)
+{
+ GeoLib::Point pnt;
+ std::vector<MeshLib::Node*> nodes (createNodes());
+ std::array<MeshLib::Node*, 8> hex_nodes;
+ std::copy(nodes.begin(), nodes.end(), hex_nodes.begin());
+ MeshLib::Hex hex(hex_nodes);
+
+ pnt = GeoLib::Point(0.99,0.99,0.99);
+ ASSERT_EQ(true, hex.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.99,0,0);
+ ASSERT_EQ(true, hex.isPntInElement(pnt));
+ pnt = GeoLib::Point(0.0, 0.0, 0.0);
+ ASSERT_EQ(true, hex.isPntInElement(pnt));
+ pnt = GeoLib::Point(1.01,0.99,0.99);
+ ASSERT_EQ(false, hex.isPntInElement(pnt));
+ ASSERT_EQ(true, hex.isPntInElement(pnt, 0.02));
+ deleteNodes(nodes);
+}