diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 334434e499d844520a8dd93e88c5f010a9a3bb3a..e7dcf9ab60d151431227fb8c56b9bf961ebe41de 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -314,48 +314,6 @@ std::unique_ptr<GeoLib::Point> triangleLineIntersection(
u * a[2] + v * b[2] + w * c[2])};
}
-bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
- const MathLib::Point3d& c, const MathLib::Point3d& d)
-{
- for (unsigned x=0; x<3; ++x)
- {
- const unsigned y=(x+1)%3;
- const double abc = (b[x] - a[x])*(c[y] - a[y]) - (b[y] - a[y])*(c[x] - a[x]);
- const double abd = (b[x] - a[x])*(d[y] - a[y]) - (b[y] - a[y])*(d[x] - a[x]);
-
- if ((abc>0 && abd<0) || (abc<0 && abd>0))
- return true;
- }
- return false;
-}
-
-bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
- const MathLib::Point3d& c, const MathLib::Point3d& d)
-{
- const MathLib::Vector3 ab(a,b);
- const MathLib::Vector3 ac(a,c);
- const MathLib::Vector3 ad(a,d);
-
- if (ab.getSqrLength() < pow(std::numeric_limits<double>::epsilon(),2) ||
- ac.getSqrLength() < pow(std::numeric_limits<double>::epsilon(),2) ||
- ad.getSqrLength() < pow(std::numeric_limits<double>::epsilon(),2)) {
- return true;
- }
-
- // In exact arithmetic <ac*ad^T, ab> should be zero
- // if all four points are coplanar.
- const double sqr_scalar_triple(pow(MathLib::scalarProduct(MathLib::crossProduct(ac,ad), ab),2));
- // Due to evaluating the above numerically some cancellation or rounding
- // can occur. For this reason a normalisation factor is introduced.
- const double normalisation_factor =
- (ab.getSqrLength() * ac.getSqrLength() * ad.getSqrLength());
-
- // tolerance 1e-11 is choosen such that
- // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-6) are considered as coplanar
- // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-5) are considered as not coplanar
- return (sqr_scalar_triple/normalisation_factor < 1e-11);
-}
-
void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec,
std::vector<GeoLib::Polyline*> & plys)
{
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 83bac893dd2d963925ed8846d255d18fd6bcfed9..460b9bd5223e439f30bbaf6261daff76af45b964 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -227,22 +227,6 @@ std::unique_ptr<GeoLib::Point> triangleLineIntersection(
MathLib::Point3d const& c, MathLib::Point3d const& p,
MathLib::Point3d const& q);
-/**
- * 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
- * @param b second point on plane
- * @param c first point to test
- * @param d second point to test
- * @return true, if such a plane can be found, false otherwise
- */
- bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
- const MathLib::Point3d& c, const MathLib::Point3d& d);
-
- /// Checks if the four given points are located on a plane.
- bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
- const MathLib::Point3d& c, const MathLib::Point3d & 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/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp
index e243ad5f1171f01962ce829c36e8f56de3604fc3..2482792c8ad045d90d367b3623f1913e58b64818 100644
--- a/GeoLib/MinimalBoundingSphere.cpp
+++ b/GeoLib/MinimalBoundingSphere.cpp
@@ -17,6 +17,7 @@
#include <ctime>
#include "MathLib/Point3d.h"
+#include "MathLib/GeometricBasics.h"
#include "MathLib/Vector3.h"
#include "AnalyticalGeometry.h"
@@ -86,7 +87,7 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p,
MathLib::Vector3 const b(p, r);
MathLib::Vector3 const c(p, s);
- if (!GeoLib::isCoplanar(p, q, r, s))
+ if (!MathLib::isCoplanar(p, q, r, s))
{
double const denom = 2.0 * MathLib::scalarTriple(a,b,c);
MathLib::Vector3 const o = (scalarProduct(c,c) * crossProduct(a,b)
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index c00836551df11e2f21e1238230596e7d88bf40eb..5f0df0f0482bf4ccd8151ad1fdc6e96e9f380a8e 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -19,6 +19,7 @@
#include "Polyline.h"
#include "AnalyticalGeometry.h"
+#include "MathLib/GeometricBasics.h"
namespace GeoLib
{
@@ -208,7 +209,7 @@ bool Polyline::isCoplanar() const
GeoLib::Point const& p2 (*this->getPoint(2));
for (std::size_t i=3; i<n_points; ++i)
{
- if (!GeoLib::isCoplanar(p0, p1, p2, *this->getPoint(i)))
+ if (!MathLib::isCoplanar(p0, p1, p2, *this->getPoint(i)))
{
DBUG ("Point %d is not coplanar to the first three points of the line.", i);
return false;
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index 11e2fcf4e25cebb751a64cd231ab055b6336b30c..a90c5bf86f08bc3d02e6b080681273dd758dbb48 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -169,4 +169,52 @@ bool barycentricPointInTriangle(MathLib::Point3d const& p,
return true;
}
+bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
+ const MathLib::Point3d& c, const MathLib::Point3d& d)
+{
+ for (unsigned x = 0; x < 3; ++x)
+ {
+ const unsigned y = (x + 1) % 3;
+ const double abc =
+ (b[x] - a[x]) * (c[y] - a[y]) - (b[y] - a[y]) * (c[x] - a[x]);
+ const double abd =
+ (b[x] - a[x]) * (d[y] - a[y]) - (b[y] - a[y]) * (d[x] - a[x]);
+
+ if ((abc > 0 && abd < 0) || (abc < 0 && abd > 0))
+ return true;
+ }
+ return false;
+}
+
+bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
+ const MathLib::Point3d& c, const MathLib::Point3d& d)
+{
+ const MathLib::Vector3 ab(a, b);
+ const MathLib::Vector3 ac(a, c);
+ const MathLib::Vector3 ad(a, d);
+
+ if (ab.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2) ||
+ ac.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2) ||
+ ad.getSqrLength() < pow(std::numeric_limits<double>::epsilon(), 2))
+ {
+ return true;
+ }
+
+ // In exact arithmetic <ac*ad^T, ab> should be zero
+ // if all four points are coplanar.
+ const double sqr_scalar_triple(
+ pow(MathLib::scalarProduct(MathLib::crossProduct(ac, ad), ab), 2));
+ // Due to evaluating the above numerically some cancellation or rounding
+ // can occur. For this reason a normalisation factor is introduced.
+ const double normalisation_factor =
+ (ab.getSqrLength() * ac.getSqrLength() * ad.getSqrLength());
+
+ // tolerance 1e-11 is choosen such that
+ // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-6) are considered as
+ // coplanar
+ // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-5) are considered as not
+ // coplanar
+ return (sqr_scalar_triple / normalisation_factor < 1e-11);
+}
+
} // end namespace MathLib
diff --git a/MathLib/GeometricBasics.h b/MathLib/GeometricBasics.h
index d1c45460ad380e541fc2eda10816b137d77de7ea..80a5783e559250089acf8c8acbc23676d2bca63b 100644
--- a/MathLib/GeometricBasics.h
+++ b/MathLib/GeometricBasics.h
@@ -154,6 +154,23 @@ bool barycentricPointInTriangle(
double eps_pnt_out_of_plane = std::numeric_limits<float>::epsilon(),
double eps_pnt_out_of_tri = std::numeric_limits<float>::epsilon());
+/**
+ * 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
+ * @param b second point on plane
+ * @param c first point to test
+ * @param d second point to test
+ * @return true, if such a plane can be found, false otherwise
+ */
+bool dividedByPlane(const MathLib::Point3d& a, const MathLib::Point3d& b,
+ const MathLib::Point3d& c, const MathLib::Point3d& d);
+
+/// Checks if the four given points are located on a plane.
+bool isCoplanar(const MathLib::Point3d& a, const MathLib::Point3d& b,
+ const MathLib::Point3d& c, const MathLib::Point3d& d);
+
} // end namespace MathLib
#endif
diff --git a/MeshLib/Elements/QuadRule4.cpp b/MeshLib/Elements/QuadRule4.cpp
index d370adc593c608e819105708c6644428225400e0..b010d5ba0949dae8de3b702131be512a30168ce7 100644
--- a/MeshLib/Elements/QuadRule4.cpp
+++ b/MeshLib/Elements/QuadRule4.cpp
@@ -11,8 +11,6 @@
#include "logog/include/logog.hpp"
-#include "GeoLib/AnalyticalGeometry.h"
-
#include "MathLib/GeometricBasics.h"
#include "MeshLib/Node.h"
@@ -57,14 +55,19 @@ unsigned QuadRule4::identifyFace(Node const* const* _nodes, Node* nodes[3])
ElementErrorCode QuadRule4::validate(const Element* e)
{
ElementErrorCode error_code;
- error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
+ error_code[ElementErrorFlag::ZeroVolume] = e->hasZeroVolume();
Node const* const* _nodes = e->getNodes();
- error_code[ElementErrorFlag::NonCoplanar] = (!GeoLib::isCoplanar(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]));
- // for collapsed quads (i.e. reduced to a line) this test might result "false" as all four points are actually located on a line.
+ error_code[ElementErrorFlag::NonCoplanar] =
+ (!MathLib::isCoplanar(*_nodes[0], *_nodes[1], *_nodes[2], *_nodes[3]));
+ // for collapsed quads (i.e. reduced to a line) this test might result
+ // "false" as all four points are actually located on a line.
if (!error_code[ElementErrorFlag::ZeroVolume])
- error_code[ElementErrorFlag::NonConvex] = (!(GeoLib::dividedByPlane(*_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&
- GeoLib::dividedByPlane(*_nodes[1], *_nodes[3], *_nodes[0], *_nodes[2])));
- error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder();
+ error_code[ElementErrorFlag::NonConvex] =
+ (!(MathLib::dividedByPlane(
+ *_nodes[0], *_nodes[2], *_nodes[1], *_nodes[3]) &&
+ MathLib::dividedByPlane(
+ *_nodes[1], *_nodes[3], *_nodes[0], *_nodes[2])));
+ error_code[ElementErrorFlag::NodeOrder] = !e->testElementNodeOrder();
return error_code;
}
diff --git a/MeshLib/MeshEditing/MeshRevision.cpp b/MeshLib/MeshEditing/MeshRevision.cpp
index 97106b291922ca163e5e78e3ce8358012501be44..9ab05f9d075848847bef2bd828aa4df1596840ff 100644
--- a/MeshLib/MeshEditing/MeshRevision.cpp
+++ b/MeshLib/MeshEditing/MeshRevision.cpp
@@ -16,11 +16,11 @@
#include <numeric>
-// ThirdParty/logog
#include "logog/include/logog.hpp"
#include "GeoLib/Grid.h"
#include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/GeometricBasics.h"
#include "MeshLib/Elements/Elements.h"
#include "MeshLib/Mesh.h"
@@ -657,10 +657,10 @@ unsigned MeshRevision::reducePrism(MeshLib::Element const*const org_elem,
addTetrahedron(i_offset, j_offset, k_offset, i);
const unsigned l =
- (GeoLib::isCoplanar(*org_elem->getNode(i_offset),
- *org_elem->getNode(k_offset),
- *org_elem->getNode(i),
- *org_elem->getNode(k)))
+ (MathLib::isCoplanar(*org_elem->getNode(i_offset),
+ *org_elem->getNode(k_offset),
+ *org_elem->getNode(i),
+ *org_elem->getNode(k)))
? j
: i;
const unsigned l_offset = (i>2) ? l-3 : l+3;
@@ -754,7 +754,8 @@ MeshLib::Element* MeshRevision::constructFourNodeElement(
}
// test if quad or tet
- const bool isQuad (GeoLib::isCoplanar(*new_nodes[0], *new_nodes[1], *new_nodes[2], *new_nodes[3]));
+ const bool isQuad(MathLib::isCoplanar(*new_nodes[0], *new_nodes[1],
+ *new_nodes[2], *new_nodes[3]));
if (isQuad && min_elem_dim < 3)
{
MeshLib::Element* elem (new MeshLib::Quad(new_nodes));
diff --git a/Tests/GeoLib/TestDividedByPlane.cpp b/Tests/GeoLib/TestDividedByPlane.cpp
deleted file mode 100644
index 3e6121dc7636e2052b94f3e4548704e43f2151db..0000000000000000000000000000000000000000
--- a/Tests/GeoLib/TestDividedByPlane.cpp
+++ /dev/null
@@ -1,49 +0,0 @@
-/**
- * \file TestDividedByPlane.cpp
- * \author Karsten Rink
- * \date 2014-02-26
- *
- * \copyright
- * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- *
- */
-
-#include "gtest/gtest.h"
-
-#include "GeoLib/AnalyticalGeometry.h"
-#include "GeoLib/Point.h"
-
-TEST(GeoLib, TestDividedByPlane)
-{
- // xy plane
- GeoLib::Point a(0,0,0);
- GeoLib::Point b(1,0,0);
- GeoLib::Point c(0,1,0);
- GeoLib::Point d(1,1,0);
-
- bool result = GeoLib::dividedByPlane(a, d, b, c);
- ASSERT_TRUE(result);
-
- result = GeoLib::dividedByPlane(b, c, a, d);
- ASSERT_TRUE(result);
-
- d = GeoLib::Point(0.1, 0.1, 0);
- result = GeoLib::dividedByPlane(b, c, a, d);
- ASSERT_FALSE(result);
-
- // xz plane
- c = GeoLib::Point(0, 0, 1);
- d = GeoLib::Point(1, 0, 1);
- result = GeoLib::dividedByPlane(a, d, b, c);
- ASSERT_TRUE(result);
-
- // yz plane
- b = GeoLib::Point(0, 1, 0);
- d = GeoLib::Point(0, 1, 1);
- result = GeoLib::dividedByPlane(a, d, b, c);
- ASSERT_TRUE(result);
-}
diff --git a/Tests/GeoLib/TestPointsOnAPlane.cpp b/Tests/GeoLib/TestPointsOnAPlane.cpp
index a2d1dbcb2d872de1e263bd106da7fb33a0391c4f..2c1a9e05f0656284877b9f164cbe710ec62673be 100644
--- a/Tests/GeoLib/TestPointsOnAPlane.cpp
+++ b/Tests/GeoLib/TestPointsOnAPlane.cpp
@@ -15,9 +15,9 @@
#include <ctime>
#include <random>
-#include "gtest/gtest.h"
+#include <gtest/gtest.h>
-#include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/GeometricBasics.h"
#include "GeoLib/Point.h"
void testAllPossibilities(
@@ -27,33 +27,33 @@ void testAllPossibilities(
GeoLib::Point const& d,
bool expected)
{
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, b, c, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, b, d, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, c, b, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, c, d, b));
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, d, b, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(a, d, c, b));
-
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, a, c, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, a, d, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, c, a, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, c, d, a));
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, d, a, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(b, d, c, a));
-
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, b, a, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, b, d, a));
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, a, b, d));
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, a, d, b));
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, d, b, a));
- ASSERT_EQ(expected, GeoLib::isCoplanar(c, d, a, b));
-
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, b, c, a));
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, b, a, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, c, b, a));
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, c, a, b));
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, a, b, c));
- ASSERT_EQ(expected, GeoLib::isCoplanar(d, a, c, b));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, b, c, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, b, d, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, c, b, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, c, d, b));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, d, b, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(a, d, c, b));
+
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, a, c, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, a, d, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, c, a, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, c, d, a));
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, d, a, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(b, d, c, a));
+
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, b, a, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, b, d, a));
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, a, b, d));
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, a, d, b));
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, d, b, a));
+ ASSERT_EQ(expected, MathLib::isCoplanar(c, d, a, b));
+
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, b, c, a));
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, b, a, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, c, b, a));
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, c, a, b));
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, a, b, c));
+ ASSERT_EQ(expected, MathLib::isCoplanar(d, a, c, b));
}
TEST(GeoLib, TestPointsOnAPlane)
diff --git a/Tests/MathLib/TestDividedByPlane.cpp b/Tests/MathLib/TestDividedByPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..94325588bb8c4b046a5af2d03e8628a074050c9c
--- /dev/null
+++ b/Tests/MathLib/TestDividedByPlane.cpp
@@ -0,0 +1,49 @@
+/**
+ * \file TestDividedByPlane.cpp
+ * \author Karsten Rink
+ * \date 2014-02-26
+ *
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ *
+ */
+
+#include <gtest/gtest.h>
+
+#include "MathLib/GeometricBasics.h"
+#include "MathLib/Point3d.h"
+
+TEST(MathLib, TestDividedByPlane)
+{
+ // xy plane
+ MathLib::Point3d a(std::array<double,3>{{0, 0, 0}});
+ MathLib::Point3d b(std::array<double,3>{{1, 0, 0}});
+ MathLib::Point3d c(std::array<double,3>{{0, 1, 0}});
+ MathLib::Point3d d(std::array<double,3>{{1, 1, 0}});
+
+ bool result = MathLib::dividedByPlane(a, d, b, c);
+ ASSERT_TRUE(result);
+
+ result = MathLib::dividedByPlane(b, c, a, d);
+ ASSERT_TRUE(result);
+
+ d = MathLib::Point3d(std::array<double,3>{{0.1, 0.1, 0}});
+ result = MathLib::dividedByPlane(b, c, a, d);
+ ASSERT_FALSE(result);
+
+ // xz plane
+ c = MathLib::Point3d(std::array<double,3>{{0, 0, 1}});
+ d = MathLib::Point3d(std::array<double,3>{{1, 0, 1}});
+ result = MathLib::dividedByPlane(a, d, b, c);
+ ASSERT_TRUE(result);
+
+ // yz plane
+ b = MathLib::Point3d(std::array<double,3>{{0, 1, 0}});
+ d = MathLib::Point3d(std::array<double,3>{{0, 1, 1}});
+ result = MathLib::dividedByPlane(a, d, b, c);
+ ASSERT_TRUE(result);
+}