Merge branch 'MathLibCodeImprovements' into 'master'

MathLib code improvements

See merge request ogs/ogs!3289

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

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);

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])

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

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

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)

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}, "

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);
 }
 

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();
                 }