move isPointOnLine

dune/iga/geometrykernel/higherorderalgorithms.hh → dune/iga/geometrykernel/algorithms.hh

@@ -178,4 +178,39 @@ namespace Dune::IGANEW::GeometryKernel {
     return dot(basis, localControlPointNet);
   }
 
+  /**
+ * @brief Checks if a point lies on a line defined by two other points with a specified tolerance.
+ *
+ * @details This function determines if a given point, specified by the local coordinates 'pt',
+ * lies on the line defined by two other points 'p1' and 'p2'. It uses the slope formula to check if
+ * the points are collinear and then verifies if the given point lies on the line within the specified tolerance.
+ *
+ * @tparam Coordinate Type representing the local coordinates of points.
+ * @tparam Tolerance Type representing the tolerance for comparing coordinates.
+ *
+ * @param p The point to check.
+ * @param linePoint0 First point defining the line.
+ * @param linePoint1 Second point defining the line.
+ * @param tol Tolerance for comparing coordinates. Default is 1e-6.
+ *
+ * @return bool - True if the point lies on the line, false otherwise.
+ */
+  template <typename Coordinate, std::floating_point Tolerance=double_t>
+  auto isPointOnLineSegment(const Coordinate& p, const Coordinate& linePoint0, const Coordinate& linePoint1,Tolerance tol =1e-6) -> bool {
+    const auto crossProductZ = [](const Coordinate& a, const Coordinate& b) {
+      return  a[0]*b[1] - a[1]*b[0];
+    };
+    if (std::abs(crossProductZ(linePoint1-linePoint0,p-linePoint0)) < tol) {
+      // The lines "linePoint1-linePoint0" and "p-linePoint0" are collinear, now check if p is within the line segment
+      auto withinRange = [](auto a, auto b, auto c) {
+        const auto& [left,right] = std::minmax(b, c);
+        return a >= left and a <= right;
+      };
+
+      if (withinRange(p[0], linePoint0[0], linePoint1[0]) and withinRange(p[1], linePoint0[1], linePoint1[1]))
+        return true;
+    }
+    return false;
+  }
+
 }  // namespace Dune::IGANEW::GeometryKernel

dune/iga/geometrykernel/nurbspatchgeometry.hh

@@ -18,7 +18,7 @@
 
 #include "dune/iga/hierarchicpatch/patchgridfwd.hh"
 #include <dune/iga/geometrykernel/geohelper.hh>
-#include <dune/iga/geometrykernel/higherorderalgorithms.hh>
+#include <dune/iga/geometrykernel/algorithms.hh>
 #include <dune/iga/geometrykernel/nurbspatchgeometrylocalview.hh>
 
 namespace Dune {

dune/iga/geometrykernel/nurbspatchgeometrylocalview.hh

@@ -12,7 +12,7 @@
 
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/iga/geometrykernel/higherorderalgorithms.hh>
+#include <dune/iga/geometrykernel/algorithms.hh>
 #include <dune/iga/hierarchicpatch/enums.hh>
 // #include <dune/iga/hierarchicpatch/hierachicpatchgridlocalgeometry.hh>
 // #include <dune/iga/hierarchicpatch/patchgridentity.hh>

dune/iga/test/testibrareader.cpp

@@ -36,7 +36,7 @@ auto testIbraReader() {
   for (auto& [name, min, max] : testCases) {
     for (int i = min; i <= max; i++) {
       gridFactory.insertJson(name, true, {i, i});
-      t.checkNoThrow([&]{auto grid = gridFactory.createGrid();});
+      t.checkNoThrow([&]{auto grid = gridFactory.createGrid();})<<"createGrid should not throw";
     }
   }
 

dune/iga/trimmer/defaulttrimmer/trimmingutils/clipelementrectangle.hh

@@ -117,7 +117,7 @@ namespace Dune::IGANEW::DefaultTrim::Util {
 
         for (auto e = 0; e < edgeLookUp.size(); ++e) {
           if (const auto& edgeIdx = edgeLookUp[e];
-              isPointOnLine(pt, eleGeo.corner(edgeIdx.front()), eleGeo.corner(edgeIdx.back()))
+              GeometryKernel::isPointOnLineSegment(pt, eleGeo.corner(edgeIdx.front()), eleGeo.corner(edgeIdx.back()))
               && !checkParallel(curve, e)) {
             result.addNewVertex(e, ptClipper, patchTrimData.getZValue(cI,0));
             break;

dune/iga/trimmer/defaulttrimmer/trimmingutils/cliputils.hh

@@ -24,29 +24,6 @@ namespace Dune::IGANEW::DefaultTrim::Util {
     return std::make_pair(it != eleRect.end(), std::ranges::distance(eleRect.begin(), it));
   }
 
-  auto isPointOnLine(const auto& pt, const auto& p1, const auto& p2) -> bool {
-    // Check if the points are collinear using the slope formula
-    auto x  = pt[0];
-    auto y  = pt[1];
-    auto x1 = p1[0];
-    auto y1 = p1[1];
-    auto x2 = p2[0];
-    auto y2 = p2[1];
-
-    // Check if the line is horizontal
-    if (std::abs(y1 - y2) < 1e-6) {
-      return std::abs(y - y1) < 1e-6 && ((x >= x1 && x <= x2) || (x >= x2 && x <= x1));
-    }
-    // Check if the line is vertical
-    if (std::abs(x1 - x2) < 1e-6) {
-      return std::abs(x - x1) < 1e-6 && ((y >= y1 && y <= y2) || (y >= y2 && y <= y1));
-    }
-
-    // Check if the given point lies on the line
-    const double slope = (y2 - y1) / (x2 - x1);
-    return std::abs((y - y1) - slope * (x - x1)) < 1e-6;
-  }
-
   struct ClippingResult {
     explicit ClippingResult(const std::vector<Clipper2Lib::PointD>& oldV) : originalVertices_(oldV) {}