LineSegment.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "LineSegment.h"

#include "JsonIO.h"

#include "libmesh/plane.h"
#include "libmesh/vector_value.h"

using libMesh::Point;
using libMesh::Real;
using libMesh::RealVectorValue;

LineSegment::LineSegment(const Point & p0, const Point & p1) : _p0(p0), _p1(p1) {}

bool
LineSegment::closest_point(const Point & p, bool clamp_to_segment, Point & closest_p) const
{
  Point p0_p = p - _p0;
  Point p0_p1 = _p1 - _p0;
  Real p0_p1_2 = p0_p1.norm_sq();
  Real perp = p0_p(0) * p0_p1(0) + p0_p(1) * p0_p1(1) + p0_p(2) * p0_p1(2);
  Real t = perp / p0_p1_2;
  bool on_segment = true;

  if (t < 0.0 || t > 1.0)
    on_segment = false;

  if (clamp_to_segment)
  {
    if (t < 0.0)
      t = 0.0;
    else if (t > 1.0)
      t = 1.0;
  }

  closest_p = _p0 + p0_p1 * t;
  return on_segment;
}

Point
LineSegment::closest_point(const Point & p) const
{
  Point closest_p;
  closest_point(p, true, closest_p);
  return closest_p;
}

bool
LineSegment::closest_normal_point(const Point & p, Point & closest_p) const
{
  return closest_point(p, false, closest_p);
}

bool
LineSegment::contains_point(const Point & p) const
{
  Point closest_p;
  return closest_point(p, false, closest_p) && closest_p.absolute_fuzzy_equals(p);
}

bool
LineSegment::intersect(const libMesh::Plane & pl, Point & intersect_p) const
{
  Point pl0 = pl.get_planar_point();
  RealVectorValue N = pl.unit_normal(_p0);
  RealVectorValue I = (_p1 - _p0).unit();

  Real numerator = (pl0 - _p0) * N;
  Real denominator = I * N;

  // The Line is parallel to the plane
  if (std::abs(denominator) < 1.e-10)
  {
    // The Line is on the plane
    if (std::abs(numerator) < 1.e-10)
    {
      // The solution is not unique so we'll just pick an end point for now
      intersect_p = _p0;
      return true;
    }
    return false;
  }

  Real d = numerator / denominator;

  // Make sure we haven't moved off the line segment!
  if (d + libMesh::TOLERANCE < 0 || d - libMesh::TOLERANCE > (_p1 - _p0).norm())
    return false;

  intersect_p = d * I + _p0;

  return true;
}

bool
LineSegment::intersect(const LineSegment & l, Point & intersect_p) const
{
  RealVectorValue a = _p1 - _p0;
  RealVectorValue b = l._p1 - l._p0;
  RealVectorValue c = l._p0 - _p0;

  RealVectorValue v = a.cross(b);

  const auto tol = 1.e-10;

  // Parallel lines check
  if (v.norm() < tol)
  {
    // Parallel but not collinear: no intersection
    if (c.cross(a).norm() >= tol)
      return false;

    // Collinear: segments overlap if any endpoint lies on the other segment
    const bool overlap = this->contains_point(l._p0) || this->contains_point(l._p1) ||
                         l.contains_point(_p0) || l.contains_point(_p1);

    if (!overlap)
      return false;

    // Pick a deterministic intersection point on the overlap
    if (this->contains_point(l._p0))
      intersect_p = l._p0;
    else if (this->contains_point(l._p1))
      intersect_p = l._p1;
    else if (l.contains_point(_p0))
      intersect_p = _p0;
    else
      intersect_p = _p1;

    return true;
  }

  // Check that the lines are coplanar
  Real concur = c * (a.cross(b));
  if (std::abs(concur) > 1.e-10)
    return false;

  Real s = (c.cross(b) * v) / (v * v);
  Real t = (c.cross(a) * v) / (v * v);

  // if s and t are between 0 and 1 then the Line Segments intersect
  // TODO: We could handle other case of clamping to the end of Line
  //       Segements if we want to here

  if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
  {
    intersect_p = _p0 + s * a;
    return true;
  }
  return false;

}

bool
LineSegment::intersect(const LineSegment & line_segment) const
{
  Point p;
  return intersect(line_segment, p);
}

void
LineSegment::set(const Point & p0, const Point & p1)
{
  setStart(p0);
  setEnd(p1);
}

void
dataStore(std::ostream & stream, LineSegment & l, void * context)
{
  dataStore(stream, l.start(), context);
  dataStore(stream, l.end(), context);
}

void
dataLoad(std::istream & stream, LineSegment & l, void * context)
{
  Point p0;
  dataLoad(stream, p0, context);
  Point p1;
  dataLoad(stream, p1, context);
  l.set(p0, p1);
}

void
to_json(nlohmann::json & json, const LineSegment & l)
{
  to_json(json["start"], l.start());
  to_json(json["end"], l.end());
}

Point
LineSegment::normal() const
{
  const auto tangent = _p1 - _p0;

  // Rotate 90 degrees counter-clockwise (2D)
  Point n(-tangent(1), tangent(0), 0.0);
  n /= n.norm();

  return n;
}

Ball
LineSegment::computeBoundingBall() const
{
  const Point center = 0.5 * (_p0 + _p1);
  const Real radius = 0.5 * (_p1 - _p0).norm();
  return Ball(center, radius);
}