Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "LineSegment.h"
      11             : 
      12             : #include "JsonIO.h"
      13             : 
      14             : #include "libmesh/plane.h"
      15             : #include "libmesh/vector_value.h"
      16             : 
      17         967 : LineSegment::LineSegment(const Point & p0, const Point & p1) : _p0(p0), _p1(p1) {}
      18             : 
      19             : bool
      20         479 : LineSegment::closest_point(const Point & p, bool clamp_to_segment, Point & closest_p) const
      21             : {
      22         479 :   Point p0_p = p - _p0;
      23         479 :   Point p0_p1 = _p1 - _p0;
      24         479 :   Real p0_p1_2 = p0_p1.norm_sq();
      25         479 :   Real perp = p0_p(0) * p0_p1(0) + p0_p(1) * p0_p1(1) + p0_p(2) * p0_p1(2);
      26         479 :   Real t = perp / p0_p1_2;
      27         479 :   bool on_segment = true;
      28             : 
      29         479 :   if (t < 0.0 || t > 1.0)
      30         136 :     on_segment = false;
      31             : 
      32         479 :   if (clamp_to_segment)
      33             :   {
      34         109 :     if (t < 0.0)
      35          16 :       t = 0.0;
      36          93 :     else if (t > 1.0)
      37          33 :       t = 1.0;
      38             :   }
      39             : 
      40         479 :   closest_p = _p0 + p0_p1 * t;
      41         479 :   return on_segment;
      42             : }
      43             : 
      44             : Point
      45         109 : LineSegment::closest_point(const Point & p) const
      46             : {
      47         109 :   Point closest_p;
      48         109 :   closest_point(p, true, closest_p);
      49         109 :   return closest_p;
      50             : }
      51             : 
      52             : bool
      53          33 : LineSegment::closest_normal_point(const Point & p, Point & closest_p) const
      54             : {
      55          33 :   return closest_point(p, false, closest_p);
      56             : }
      57             : 
      58             : bool
      59         337 : LineSegment::contains_point(const Point & p) const
      60             : {
      61         337 :   Point closest_p;
      62         337 :   return closest_point(p, false, closest_p) && closest_p.absolute_fuzzy_equals(p);
      63             : }
      64             : 
      65             : bool
      66         598 : LineSegment::intersect(const libMesh::Plane & pl, Point & intersect_p) const
      67             : {
      68             :   /**
      69             :    * There are three cases in 3D for intersection of a line and a plane
      70             :    * Case 1: The line is parallel to the plane - No intersection
      71             :    *         Numerator = non-zero
      72             :    *         Denominator = zero
      73             :    *
      74             :    * Case 2: The line is within the plane - Inf intersection
      75             :    *         Numerator = zero
      76             :    *         Denominator = zero
      77             :    *
      78             :    * Case 3: The line intersects the plane at a single point
      79             :    *         Denominator = non-zero
      80             :    */
      81             : 
      82         598 :   Point pl0 = pl.get_planar_point();
      83         598 :   RealVectorValue N = pl.unit_normal(_p0);
      84         598 :   RealVectorValue I = (_p1 - _p0).unit();
      85             : 
      86         598 :   Real numerator = (pl0 - _p0) * N;
      87         598 :   Real denominator = I * N;
      88             : 
      89             :   // The Line is parallel to the plane
      90         598 :   if (std::abs(denominator) < 1.e-10)
      91             :   {
      92             :     // The Line is on the plane
      93         451 :     if (std::abs(numerator) < 1.e-10)
      94             :     {
      95             :       // The solution is not unique so we'll just pick an end point for now
      96           5 :       intersect_p = _p0;
      97           5 :       return true;
      98             :     }
      99         446 :     return false;
     100             :   }
     101             : 
     102         147 :   Real d = numerator / denominator;
     103             : 
     104             :   // Make sure we haven't moved off the line segment!
     105         147 :   if (d + libMesh::TOLERANCE < 0 || d - libMesh::TOLERANCE > (_p1 - _p0).norm())
     106          48 :     return false;
     107             : 
     108          99 :   intersect_p = d * I + _p0;
     109             : 
     110          99 :   return true;
     111             : }
     112             : 
     113             : bool
     114         378 : LineSegment::intersect(const LineSegment & l, Point & intersect_p) const
     115             : {
     116             :   /**
     117             :    * First check for concurance:
     118             :    *
     119             :    *
     120             :    * | x1 y1 z1 1 |
     121             :    * | x2 y2 z2 1 | = (x3 - x1) * [(x2-x1) x (x4-x3)] = 0
     122             :    * | x3 y3 z3 1 |
     123             :    * | x4 y4 z4 1 |
     124             :    *
     125             :    *
     126             :    * Solve:
     127             :    *   x = _p0 + (_p1 - _p0)*s
     128             :    *   x = l.p0 + (l._p1 - l.p0)*t
     129             :    *
     130             :    *   where
     131             :    *   a = _p1 - _p0
     132             :    *   b = l._p1 - l._p0
     133             :    *   c = l._p0 - _p0
     134             :    *
     135             :    *   s = (c x b) * (a x b) / | a x b |^2
     136             :    */
     137         378 :   RealVectorValue a = _p1 - _p0;
     138         378 :   RealVectorValue b = l._p1 - l._p0;
     139         378 :   RealVectorValue c = l._p0 - _p0;
     140             : 
     141         378 :   RealVectorValue v = a.cross(b);
     142             : 
     143             :   // Check for parallel lines
     144         378 :   if (std::abs(v.norm()) < 1.e-10 && std::abs(c.cross(a).norm()) < 1.e-10)
     145             :   {
     146             :     // TODO: The lines are co-linear but more work is needed to determine and intersection point
     147             :     //       it could be the case that the line segments don't lie on the same line or overlap only
     148             :     //       a bit
     149           4 :     return true;
     150             :   }
     151             : 
     152             :   // Check that the lines are coplanar
     153         374 :   Real concur = c * (a.cross(b));
     154         374 :   if (std::abs(concur) > 1.e-10)
     155           8 :     return false;
     156             : 
     157         366 :   Real s = (c.cross(b) * v) / (v * v);
     158         366 :   Real t = (c.cross(a) * v) / (v * v);
     159             : 
     160             :   // if s and t are between 0 and 1 then the Line Segments intersect
     161             :   // TODO: We could handle other case of clamping to the end of Line
     162             :   //       Segements if we want to here
     163             : 
     164         366 :   if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
     165             :   {
     166         131 :     intersect_p = _p0 + s * a;
     167         131 :     return true;
     168             :   }
     169         235 :   return false;
     170             : 
     171             :   /**
     172             :    * Parameteric Equation of lines
     173             :    * _p0 + t(v0) = l._p0 + u(v1)
     174             :    *
     175             :    * Case 1: Parallel Lines
     176             :    *         v0 x v1 == 0
     177             :    *
     178             :    * Case 1a: Collinear Lines
     179             :    *         v0 x v1 == 0
     180             :    *         (l._p0 - _p0) x (_p1 - _p0) == 0
     181             :    *
     182             :    * Case 2: Intersecting Lines
     183             :    *         0 <= t <= 1
     184             :    *         0 <= u <= 1
     185             :    *
     186             :    *
     187             :    * Case 1: The lines do not intersect
     188             :    *         vleft cross vright = non-zero
     189             :    *
     190             :    * Case 2: The lines are co-linear
     191             :    *         vleft cross vright = zero
     192             :    *         vleft (Denominator) = zero
     193             :    *
     194             :    * Case 3: The line intersect at a single point
     195             :    *         vleft cross vright = zero
     196             :    *         vleft (Denominator) = non-zero
     197             :   RealVectorValue v0 = _p1 - _p0;
     198             :   RealVectorValue v1 = l._p1 - l._p0;
     199             :   RealVectorValue v2 = l._p0 - _p0;
     200             : 
     201             :   RealVectorValue vbot = v0.cross(v1);
     202             :   RealVectorValue vtop = v2.cross(v1);
     203             : 
     204             :   RealVectorValue crossed = vleft.cross(vright);
     205             : 
     206             :   // Case 1: No intersection
     207             :   if (std::abs(vleft.cross(vright).size()) > 1.e-10)
     208             :     return false;
     209             : 
     210             :   // Case 2: Co-linear (just return one of the end points)
     211             :   if (std::abs(vleft.size()) < 1.e-10)
     212             :   {
     213             :     intersect_p = _p0;
     214             :     return true;
     215             :   }
     216             : 
     217             :   // Case 3:
     218             : 
     219             :   //TODO: We could detect whether the Line Segments actually overlap
     220             :   //      instead of whether the Lines are co-linear
     221             : 
     222             :   Real a = vright.size()/vleft.size();
     223             :   intersect_p = _p0 + a*v0;
     224             :   return true;
     225             :      */
     226             : }
     227             : 
     228             : void
     229           0 : LineSegment::set(const Point & p0, const Point & p1)
     230             : {
     231           0 :   setStart(p0);
     232           0 :   setEnd(p1);
     233           0 : }
     234             : 
     235             : void
     236           0 : dataStore(std::ostream & stream, LineSegment & l, void * context)
     237             : {
     238           0 :   dataStore(stream, l.start(), context);
     239           0 :   dataStore(stream, l.end(), context);
     240           0 : }
     241             : 
     242             : void
     243           0 : dataLoad(std::istream & stream, LineSegment & l, void * context)
     244             : {
     245           0 :   Point p0;
     246           0 :   dataLoad(stream, p0, context);
     247           0 :   Point p1;
     248           0 :   dataLoad(stream, p1, context);
     249           0 :   l.set(p0, p1);
     250           0 : }
     251             : 
     252             : void
     253           0 : to_json(nlohmann::json & json, const LineSegment & l)
     254             : {
     255           0 :   to_json(json["start"], l.start());
     256           0 :   to_json(json["end"], l.end());
     257           0 : }