diff --git a/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp
index b29448a77d2513b564dad0a02a35bfeb20aad06d..7f99bc9f6072bed483c9baa62f13da3f5751bdec 100644
--- a/GeoLib/MinimalBoundingSphere.cpp
+++ b/GeoLib/MinimalBoundingSphere.cpp
@@ -182,7 +182,7 @@ MinimalBoundingSphere::recurseCalculation(
double MinimalBoundingSphere::pointDistanceSquared(MathLib::Point3d const& pnt) const
{
- return MathLib::sqrDist(_center.getCoords(), pnt.getCoords())-(_radius*_radius);
+ return MathLib::sqrDist(_center, pnt)-(_radius*_radius);
}
std::vector<MathLib::Point3d*>* MinimalBoundingSphere::getRandomSpherePoints(std::size_t n_points) const
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index 0c21c3cc374cf133d6899eadac6a0d4535281e04..cb617c4f0f9f0c3cc5c3c7377eb8f3e314362422 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -455,10 +455,10 @@ double Polyline::getDistanceAlongPolyline(const MathLib::Point3d& pnt,
// is the orthogonal projection of the j-th node to the
// line g(lambda) = _ply->getPoint(k) + lambda * (_ply->getPoint(k+1) - _ply->getPoint(k))
// at the k-th line segment of the polyline, i.e. 0 <= lambda <= 1?
- if (MathLib::calcProjPntToLineAndDists(pnt.getCoords(),
- (getPoint(k))->getCoords(), (getPoint(k + 1))->getCoords(),
- lambda, dist) <= epsilon_radius) {
-
+ if (MathLib::calcProjPntToLineAndDists(pnt, *getPoint(k),
+ *getPoint(k + 1), lambda,
+ dist) <= epsilon_radius)
+ {
double act_length_of_ply(getLength(k));
double seg_length (getLength(k+1)-getLength(k));
double lower_lambda (- epsilon_radius / seg_length);
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 84b7c83d4a2b934ffc34070fe07752af40fd8dd7..c9999f5a26337b29dab78bb7b77477d06cdb96eb 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -7,32 +7,39 @@
* http://www.opengeosys.org/project/license
*/
+#include "MathTools.h"
+
+#include <Eigen/Dense>
#include <cmath>
-#include "MathTools.h"
+#include "Point3d.h"
namespace MathLib
{
-double calcProjPntToLineAndDists(const double p[3], const double a[3],
- const double b[3], double &lambda, double &d0)
+double calcProjPntToLineAndDists(Point3d const& pp, Point3d const& pa,
+ Point3d const& pb, double& lambda, double& d0)
{
- // g (lambda) = a + lambda v, v = b-a
- double v[3] = {b[0] - a[0], b[1] - a[1], b[2] - a[2]};
- // orthogonal projection: (g(lambda)-p) * v = 0 => in order to compute lambda we define a help vector u
- double u[3] = {p[0] - a[0], p[1] - a[1], p[2] - a[2]};
- lambda = scalarProduct<double,3> (u, v) / scalarProduct<double,3> (v, v);
+ auto const a =
+ Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pa.getCoords()));
+ auto const b =
+ Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pb.getCoords()));
+ auto const p =
+ Eigen::Map<Eigen::Vector3d>(const_cast<double*>(pp.getCoords()));
+ // g(lambda) = a + lambda v, v = b-a
+ Eigen::Vector3d const v = b - a;
+
+ // orthogonal projection: (p - g(lambda))^T * v = 0
+ // <=> (a-p - lambda (b-a))^T * (b-a) = 0
+ // <=> (a-p)^T * (b-a) = lambda (b-a)^T ) (b-a)
+ lambda = (((p - a).transpose() * v) / v.squaredNorm())(0, 0);
// compute projected point
- double proj_pnt[3];
- for (std::size_t k(0); k < 3; k++)
- {
- proj_pnt[k] = a[k] + lambda * v[k];
- }
+ Eigen::Vector3d const proj_pnt = a + lambda * v;
- d0 = std::sqrt (sqrDist (proj_pnt, a));
+ d0 = (proj_pnt - a).norm();
- return std::sqrt (sqrDist (p, proj_pnt));
+ return (p - proj_pnt).norm();
}
double getAngle (const double p0[3], const double p1[3], const double p2[3])
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index 7925c269ea36648c624ffb43ea168e9f194910d0..361a1289c473298baaace868636b97c730add130 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -17,6 +17,9 @@
namespace MathLib
{
+template <typename T, std::size_t DIM> class TemplatePoint;
+using Point3d = MathLib::TemplatePoint<double, 3>;
+
/**
* standard inner product in R^N
* \param v0 array of type T representing the vector
@@ -65,24 +68,18 @@ inline T scalarProduct(T const* const v0, T const* const v1, int const n)
* \f$g(\lambda) = a + \lambda (b - a)\f$,
* the distance between p and the projected point
* and the distances between the projected point and the end
- * points a, b of the line
- * \param p the (mesh) point
- * \param a first point of line
- * \param b second point of line
+ * points pa, pb of the line
+ * \param pp the (mesh) point
+ * \param pa first point of line
+ * \param pb second point of line
* \param lambda the projected point described by the line equation above
* \param d0 distance to the line point a
- * \returns the distance between p and the orthogonal projection of p
+ * \returns the distance between pp and the orthogonal projection of pp
*/
-double calcProjPntToLineAndDists(const double p[3], const double a[3],
- const double b[3], double &lambda, double &d0);
-
-/** squared dist between double arrays p0 and p1 (size of arrays is 3) */
-inline
-double sqrDist(const double* p0, const double* p1)
-{
- const double v[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
- return scalarProduct<double,3>(v,v);
-}
+double calcProjPntToLineAndDists(MathLib::Point3d const& pp,
+ MathLib::Point3d const& pa,
+ MathLib::Point3d const& pb, double& lambda,
+ double& d0);
/**
* Let \f$p_0, p_1, p_2 \in R^3\f$. The function getAngle
diff --git a/MathLib/Point3d.h b/MathLib/Point3d.h
index 4f548648f75d4be7a112f89a490edf7e597f661c..d287253458da7dfdd87c9ce6e7ced4d1c1cb8de1 100644
--- a/MathLib/Point3d.h
+++ b/MathLib/Point3d.h
@@ -13,11 +13,11 @@
#pragma once
#include <limits>
+#include <Eigen/Dense>
#include "mathlib_export.h"
#include "TemplatePoint.h"
-#include "MathTools.h"
namespace MathLib
{
@@ -47,8 +47,9 @@ inline MathLib::Point3d operator*(MATRIX const& mat, MathLib::Point3d const& p)
inline
double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
{
- const double v[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
- return MathLib::scalarProduct<double,3>(v,v);
+ return (Eigen::Map<Eigen::Vector3d>(const_cast<double*>(p0.getCoords())) -
+ Eigen::Map<Eigen::Vector3d>(const_cast<double*>(p1.getCoords())))
+ .squaredNorm();
}
/// Computes the squared distance between the orthogonal projection of the two
diff --git a/MathLib/TemplatePoint.h b/MathLib/TemplatePoint.h
index 7f1fd9f2ea0815911c07ea6757941704b9e92d02..644d43706c90a3d0ce96a77fed154a646ac0aee1 100644
--- a/MathLib/TemplatePoint.h
+++ b/MathLib/TemplatePoint.h
@@ -171,17 +171,6 @@ bool lessEq(TemplatePoint<T, DIM> const& a, TemplatePoint<T, DIM> const& b,
return true;
}
-/** Distance between points p0 and p1 in the maximum norm. */
-template <typename T>
-T maxNormDist(const MathLib::TemplatePoint<T>* p0, const MathLib::TemplatePoint<T>* p1)
-{
- const T x = fabs((*p1)[0] - (*p0)[0]);
- const T y = fabs((*p1)[1] - (*p0)[1]);
- const T z = fabs((*p1)[2] - (*p0)[2]);
-
- return std::max(x, std::max(y, z));
-}
-
/** overload the output operator for class Point */
template <typename T, std::size_t DIM>
std::ostream& operator<< (std::ostream &os, const TemplatePoint<T,DIM> &p)
diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp
index 5afc280eb6d25fdb1086995b96e2d4cd1bb3cd11..7887ba457ad55f323ae60127d25cf01d8496fc86 100644
--- a/MeshGeoToolsLib/MeshNodeSearcher.cpp
+++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp
@@ -163,9 +163,7 @@ std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(
ss << "- bulk node: " << (*bulk_nodes[id])[0] << ", "
<< (*bulk_nodes[id])[1] << ", " << (*bulk_nodes[id])[2]
<< ", distance: "
- << std::sqrt(MathLib::sqrDist(bulk_nodes[id]->getCoords(),
- p.getCoords()))
- << "\n";
+ << std::sqrt(MathLib::sqrDist(*bulk_nodes[id], p)) << "\n";
}
OGS_FATAL(
"Found {:d} nodes in the mesh for point {:d} : ({:g}, {:g}, "
diff --git a/MeshLib/Elements/LineRule2.cpp b/MeshLib/Elements/LineRule2.cpp
index 37abacc47b7ee53aa2d1fc52365cd73883712260..373a8ef4b8b4146b72a8ad6b94469281638ef281 100644
--- a/MeshLib/Elements/LineRule2.cpp
+++ b/MeshLib/Elements/LineRule2.cpp
@@ -21,7 +21,7 @@ const unsigned LineRule2::edge_nodes[1][2] =
double LineRule2::computeVolume(Node const* const* _nodes)
{
- return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords()));
+ return sqrt(MathLib::sqrDist(*_nodes[0], *_nodes[1]));
}
bool LineRule2::isPntInElement(Node const* const* nodes,
@@ -30,8 +30,7 @@ bool LineRule2::isPntInElement(Node const* const* nodes,
double tmp;
double tmp_dst(0);
double const dist = MathLib::calcProjPntToLineAndDists(
- pnt.getCoords(), nodes[0]->getCoords(), nodes[1]->getCoords(), tmp,
- tmp_dst);
+ pnt, *nodes[0], *nodes[1], tmp, tmp_dst);
return (dist < eps);
}
diff --git a/MeshLib/MeshEditing/MeshRevision.cpp b/MeshLib/MeshEditing/MeshRevision.cpp
index 7783ff84722b3dc0f6292139ab6ae9eb461e6b5b..3815dc5ff6c3ebbcca96b6eca1a987383b7fab05 100644
--- a/MeshLib/MeshEditing/MeshRevision.cpp
+++ b/MeshLib/MeshEditing/MeshRevision.cpp
@@ -179,8 +179,7 @@ std::vector<std::size_t> MeshRevision::collapseNodeIndices(double eps) const
}
// calc distance
- if (MathLib::sqrDist(node->getCoords(),
- test_node->getCoords()) < sqr_eps)
+ if (MathLib::sqrDist(*node, *test_node) < sqr_eps)
{
id_map[test_node->getID()] = node->getID();
}