geom_utils Namespace Reference

Functions
bool	isPointZero (const libMesh::Point &pt)
	Check whether a point is equal to zero. More...
libMesh::Point	unitVector (const libMesh::Point &pt, const std::string &name)
	Get the unit vector for a point parameter. More...
libMesh::Point	rotatePointAboutAxis (const libMesh::Point &p, const libMesh::Real angle, const libMesh::Point &axis)
	Rotate point about an axis. More...
libMesh::Real	minDistanceToPoints (const libMesh::Point &pt, const std::vector< libMesh::Point > &candidates, const unsigned int axis)
	Get the minimum distance from a point to another set of points, in the plane perpendicular to the specified axis. More...
std::vector< libMesh::Point >	polygonCorners (const unsigned int num_sides, const libMesh::Real radius, const unsigned int axis)
	Get the corner coordinates of a regular 2-D polygon, assuming a face of the polygon is parallel to the x0 axis. More...
std::pair< unsigned int, unsigned int >	projectedIndices (const unsigned int axis)
	Get the indices of the plane perpendicular to the specified axis. More...
libMesh::Point	projectPoint (const libMesh::Real x0, const libMesh::Real x1, const unsigned int axis)
	Given two coordinates, construct a point in the 2-D plane perpendicular to the specified axis. More...
libMesh::Point	projectedUnitNormal (libMesh::Point pt1, libMesh::Point pt2, const unsigned int axis)
	Get the unit normal vector between two points (which are first projected onto the plane perpendicular to the 'axis'), such that the cross product of the unit normal with the line from pt1 to pt2 has a positive 'axis' component. More...
libMesh::Real	distanceFromLine (const libMesh::Point &pt, const libMesh::Point &line0, const libMesh::Point &line1)
	Compute the distance from a 3-D line, provided in terms of two points on the line. More...
libMesh::Real	projectedDistanceFromLine (libMesh::Point pt, libMesh::Point line0, libMesh::Point line1, const unsigned int axis)
	Compute the distance from a 3-D line, provided in terms of two points on the line. More...
libMesh::Real	projectedLineHalfSpace (libMesh::Point pt1, libMesh::Point pt2, libMesh::Point pt3, const unsigned int axis)
	If positive, point is on the positive side of the half space (and vice versa). More...
bool	pointInPolygon (const libMesh::Point &point, const std::vector< libMesh::Point > &corners, const unsigned int axis)
	Whether a point is in 2-D a polygon in the plane perpendicular to the specified axis, given by corner points. More...
bool	pointOnEdge (const libMesh::Point &point, const std::vector< libMesh::Point > &corners, const unsigned int axis)
	Whether a point is on the edge of a 2-D polygon in the plane perpendicular to the specified axis, given its corner points. More...
std::vector< libMesh::Point >	boxCorners (const libMesh::BoundingBox &box, const libMesh::Real factor)
	Get corner points of a bounding box, with side length re-scaled. More...
bool	arePointsColinear (const Point &p1, const Point &p2, const Point &p3)
	Check if three points are colinear. More...
bool	segmentsIntersect (const Point &p1, const Point &p2, const Point &p3, const Point &p4)
	Check if the line segment p1-p2 intersects with line segment p3-p4 (only working in 2D (x-y plane)). More...
Real	pointSegmentDistanceSq (const Point &point, const Point &a, const Point &b)
	Compute the squared distance from a point to a 3-D line segment. More...
Real	pointTriangleDistanceSq (const Point &point, const Point &v0, const Point &v1, const Point &v2)
	Compute the squared distance from a point to a 3-D triangle. More...
Real	solidAngle (const Point &point, const Point &v0, const Point &v1, const Point &v2)
	Compute the signed solid angle subtended by one oriented triangle at the query point. More...
bool	isPointZero (const Point &pt)
Point	unitVector (const Point &pt, const std::string &name)
Real	minDistanceToPoints (const Point &pt, const std::vector< Point > &candidates, const unsigned int axis)
Point	projectPoint (const Real x0, const Real x1, const unsigned int axis)
Real	projectedLineHalfSpace (Point pt1, Point pt2, Point pt3, const unsigned int axis)
bool	pointInPolygon (const Point &point, const std::vector< Point > &corners, const unsigned int axis)
bool	pointOnEdge (const Point &point, const std::vector< Point > &corners, const unsigned int axis)
Point	projectedUnitNormal (Point pt1, Point pt2, const unsigned int axis)
Real	distanceFromLine (const Point &pt, const Point &line0, const Point &line1)
Real	projectedDistanceFromLine (Point pt, Point line0, Point line1, const unsigned int axis)
std::vector< Point >	polygonCorners (const unsigned int num_sides, const Real radius, const unsigned int axis)
Point	rotatePointAboutAxis (const Point &p, const Real angle, const Point &axis)
std::vector< Point >	boxCorners (const libMesh::BoundingBox &box, const Real factor)

Function Documentation

arePointsColinear()

bool geom_utils::arePointsColinear	(	const Point &	p1,
		const Point &	p2,
		const Point &	p3
	)

Check if three points are colinear.

Parameters

p1	First point
p2	Second point
p3	Third point

Returns

true if the three points are colinear, false otherwise

Definition at line 285 of file GeometryUtils.C.

Referenced by BoundaryLayerUtils::collectExteriorVertexPointsFromMesh(), and segmentsIntersect().

{
  const Point v1 = p2 - p1;
  const Point v2 = p3 - p1;
  const Point cross_prod = v1.cross(v2);

  return MooseUtils::absoluteFuzzyEqual(cross_prod.norm(), 0.0);
}

boxCorners() [1/2]

std::vector<libMesh::Point> geom_utils::boxCorners	(	const libMesh::BoundingBox &	box,
		const libMesh::Real	factor
	)

Get corner points of a bounding box, with side length re-scaled.

Parameters

[in]	box	bounding box to start from
[in]	factor	by which to multiply the bounding box side

boxCorners() [2/2]

std::vector<Point> geom_utils::boxCorners	(	const libMesh::BoundingBox &	box,
		const Real	factor
	)

Definition at line 261 of file GeometryUtils.C.

{
  Point diff = (box.max() - box.min()) / 2.0;
  const Point origin = box.min() + diff;

  // Rescale side length of box by specified factor
  diff *= factor;

  // Vectors for sides of box
  const Point dx(2.0 * diff(0), 0.0, 0.0);
  const Point dy(0.0, 2.0 * diff(1), 0.0);
  const Point dz(0.0, 0.0, 2.0 * diff(2));

  std::vector<Point> verts(8, origin - diff);
  const unsigned int pts_per_dim = 2;
  for (const auto z : make_range(pts_per_dim))
    for (const auto y : make_range(pts_per_dim))
      for (const auto x : make_range(pts_per_dim))
        verts[pts_per_dim * pts_per_dim * z + pts_per_dim * y + x] += x * dx + y * dy + z * dz;

  return verts;
}

distanceFromLine() [1/2]

libMesh::Real geom_utils::distanceFromLine	(	const libMesh::Point &	pt,
		const libMesh::Point &	line0,
		const libMesh::Point &	line1
	)

Compute the distance from a 3-D line, provided in terms of two points on the line.

Parameters

[in]	pt	point of interest
[in]	line0	first point on line
[in]	line1	second point on line

Returns

distance from line

Referenced by projectedDistanceFromLine().

distanceFromLine() [2/2]

Real geom_utils::distanceFromLine	(	const Point &	pt,
		const Point &	line0,
		const Point &	line1
	)

Definition at line 192 of file GeometryUtils.C.

{
  const Point a = pt - line0;
  const Point b = pt - line1;
  const Point c = line1 - line0;

  return (a.cross(b).norm()) / c.norm();
}

isPointZero() [1/2]

bool geom_utils::isPointZero

(

const Point &

pt

)

Definition at line 21 of file GeometryUtils.C.

{
  const Point zero(0.0, 0.0, 0.0);
  return pt.absolute_fuzzy_equals(zero);
}

isPointZero() [2/2]

bool geom_utils::isPointZero

(

const libMesh::Point &

pt

)

Check whether a point is equal to zero.

Parameters

[in]

pt

point

Returns

whether point is equal to the zero point

Referenced by unitVector().

minDistanceToPoints() [1/2]

Real geom_utils::minDistanceToPoints	(	const Point &	pt,
		const std::vector< Point > &	candidates,
		const unsigned int	axis
	)

Definition at line 37 of file GeometryUtils.C.

{
  const auto idx = projectedIndices(axis);

  Real distance = std::numeric_limits<Real>::max();
  for (const auto & c : candidates)
  {
    const Real dx = c(idx.first) - pt(idx.first);
    const Real dy = c(idx.second) - pt(idx.second);
    const Real d = std::sqrt(dx * dx + dy * dy);
    distance = std::min(d, distance);
  }

  return distance;
}

minDistanceToPoints() [2/2]

libMesh::Real geom_utils::minDistanceToPoints	(	const libMesh::Point &	pt,
		const std::vector< libMesh::Point > &	candidates,
		const unsigned int	axis
	)

Get the minimum distance from a point to another set of points, in the plane perpendicular to the specified axis.

Parameters

[in]	pt	point
[in]	candidates	set of points we will find the nearest of
[in]	axis	axis perpendicular to the plane of the polygon

Returns

minimum distance to the provided points

pointInPolygon() [1/2]

bool geom_utils::pointInPolygon	(	const Point &	point,
		const std::vector< Point > &	corners,
		const unsigned int	axis
	)

Definition at line 82 of file GeometryUtils.C.

{
  const auto n_pts = corners.size();

  std::vector<bool> negative_half_space;
  std::vector<bool> positive_half_space;
  for (const auto i : index_range(corners))
  {
    const int next = (i == n_pts - 1) ? 0 : i + 1;
    const auto half = projectedLineHalfSpace(point, corners[i], corners[next], axis);
    negative_half_space.push_back(half < 0);
    positive_half_space.push_back(half > 0);
  }

  const bool negative = std::find(negative_half_space.begin(), negative_half_space.end(), true) !=
                        negative_half_space.end();
  const bool positive = std::find(positive_half_space.begin(), positive_half_space.end(), true) !=
                        positive_half_space.end();

  const bool in_polygon = !(negative && positive);
  if (in_polygon)
    return true;

  if (pointOnEdge(point, corners, axis))
    return true;

  return false;
}

pointInPolygon() [2/2]

bool geom_utils::pointInPolygon	(	const libMesh::Point &	point,
		const std::vector< libMesh::Point > &	corners,
		const unsigned int	axis
	)

Whether a point is in 2-D a polygon in the plane perpendicular to the specified axis, given by corner points.

Parameters

[in]	point	point of interest
[in]	corners	corner points of polygon
[in]	axis	axis perpendicular to the plane of the polygon

Returns

whether point is inside the polygon

pointOnEdge() [1/2]

bool geom_utils::pointOnEdge	(	const Point &	point,
		const std::vector< Point > &	corners,
		const unsigned int	axis
	)

Definition at line 112 of file GeometryUtils.C.

{
  const auto n_pts = corners.size();
  const auto idx = projectedIndices(axis);

  constexpr Real tol = 1e-8;
  for (const auto i : index_range(corners))
  {
    const int next = (i == n_pts - 1) ? 0 : i + 1;
    const auto & pt1 = corners[i];
    const auto & pt2 = corners[next];
    const bool close_to_line = projectedDistanceFromLine(point, pt1, pt2, axis) < tol;

    // we can stop early if we know we're not close to the line
    if (!close_to_line)
      continue;

    // check that the point is "between" the two points; TODO: first pass
    // we can just compare x and y coordinates
    const bool between_points = (point(idx.first) >= std::min(pt1(idx.first), pt2(idx.first))) &&
                                (point(idx.first) <= std::max(pt1(idx.first), pt2(idx.first))) &&
                                (point(idx.second) >= std::min(pt1(idx.second), pt2(idx.second))) &&
                                (point(idx.second) <= std::max(pt1(idx.second), pt2(idx.second)));

    // point needs to be close to the line AND "between" the two points
    if (close_to_line && between_points)
      return true;
  }

  return false;
}

pointOnEdge() [2/2]

bool geom_utils::pointOnEdge	(	const libMesh::Point &	point,
		const std::vector< libMesh::Point > &	corners,
		const unsigned int	axis
	)

Whether a point is on the edge of a 2-D polygon in the plane perpendicular to the specified axis, given its corner points.

Parameters

[in]	point	point of interest
[in]	corners	corner points of polygon
[in]	axis	axis perpendicular to the plane of the polygon

Returns

whether point is on edge of polygon

Referenced by pointInPolygon().

pointSegmentDistanceSq()

Real geom_utils::pointSegmentDistanceSq	(	const Point &	point,
		const Point &	a,
		const Point &	b
	)

Compute the squared distance from a point to a 3-D line segment.

Parameters

[in]	point	point of interest
[in]	a	first segment endpoint
[in]	b	second segment endpoint

Returns

squared distance from the point to the segment

Definition at line 351 of file GeometryUtils.C.

Referenced by TriangleManifold::rayIntersectsTriangle().

{
  const Point ab = b - a;
  const auto length_sq = ab.norm_sq();
  if (length_sq <= std::numeric_limits<Real>::epsilon())
    return (point - a).norm_sq();

  const auto t = std::clamp(((point - a) * ab) / length_sq, 0.0, 1.0);
  const Point projection = a + t * ab;
  return (point - projection).norm_sq();
}

pointTriangleDistanceSq()

Real geom_utils::pointTriangleDistanceSq	(	const Point &	point,
		const Point &	v0,
		const Point &	v1,
		const Point &	v2
	)

Compute the squared distance from a point to a 3-D triangle.

Parameters

[in]	point	point of interest
[in]	v0	first triangle vertex
[in]	v1	second triangle vertex
[in]	v2	third triangle vertex

Returns

squared distance from the point to the triangle

Definition at line 364 of file GeometryUtils.C.

{
  const Point ab = v1 - v0;
  const Point ac = v2 - v0;
  const Point ap = point - v0;
  const Real d1 = ab * ap;
  const Real d2 = ac * ap;
  if (d1 <= 0.0 && d2 <= 0.0)
    return (point - v0).norm_sq();

  const Point bp = point - v1;
  const Real d3 = ab * bp;
  const Real d4 = ac * bp;
  if (d3 >= 0.0 && d4 <= d3)
    return (point - v1).norm_sq();

  const Real vc = d1 * d4 - d3 * d2;
  if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
  {
    const Real v = d1 / (d1 - d3);
    const Point projection = v0 + v * ab;
    return (point - projection).norm_sq();
  }

  const Point cp = point - v2;
  const Real d5 = ab * cp;
  const Real d6 = ac * cp;
  if (d6 >= 0.0 && d5 <= d6)
    return (point - v2).norm_sq();

  const Real vb = d5 * d2 - d1 * d6;
  if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
  {
    const Real w = d2 / (d2 - d6);
    const Point projection = v0 + w * ac;
    return (point - projection).norm_sq();
  }

  const Real va = d3 * d6 - d5 * d4;
  if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
  {
    const Point bc = v2 - v1;
    const Real w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
    const Point projection = v1 + w * bc;
    return (point - projection).norm_sq();
  }

  const Real denom = 1.0 / (va + vb + vc);
  const Real v = vb * denom;
  const Real w = vc * denom;
  const Point projection = v0 + ab * v + ac * w;
  return (point - projection).norm_sq();
}

polygonCorners() [1/2]

std::vector<libMesh::Point> geom_utils::polygonCorners	(	const unsigned int	num_sides,
		const libMesh::Real	radius,
		const unsigned int	axis
	)

Get the corner coordinates of a regular 2-D polygon, assuming a face of the polygon is parallel to the x0 axis.

Parameters

[in]	num_sides	number of sides to polygon
[in]	radius	distance from polygon center to a corner
[in]	axis	axis perpendicular to the plane of the polygon

Returns

corner coordinates

polygonCorners() [2/2]

std::vector<Point> geom_utils::polygonCorners	(	const unsigned int	num_sides,
		const Real	radius,
		const unsigned int	axis
	)

Definition at line 213 of file GeometryUtils.C.

{
  std::vector<Point> corners;
  const Real theta = 2.0 * M_PI / num_sides;
  const Real first_angle = M_PI / 2.0 - theta / 2.0;

  for (const auto i : make_range(num_sides))
  {
    const Real angle = first_angle + i * theta;
    const Real x = radius * cos(angle);
    const Real y = radius * sin(angle);

    corners.push_back(projectPoint(x, y, axis));
  }

  return corners;
}

projectedDistanceFromLine() [1/2]

libMesh::Real geom_utils::projectedDistanceFromLine	(	libMesh::Point	pt,
		libMesh::Point	line0,
		libMesh::Point	line1,
		const unsigned int	axis
	)

Compute the distance from a 3-D line, provided in terms of two points on the line.

Both the input point and the points on the line are projected into the 2-d plane perpendicular to the specified axis.

Parameters

[in]	pt	point of interest
[in]	line0	first point on line
[in]	line1	second point on line
[in]	axis	axis index (0 = x, 1 = y, 2 = z) perpendicular to the projection plane

Returns

distance from line

Referenced by pointOnEdge().

projectedDistanceFromLine() [2/2]

Real geom_utils::projectedDistanceFromLine	(	Point	pt,
		Point	line0,
		Point	line1,
		const unsigned int	axis
	)

Definition at line 202 of file GeometryUtils.C.

{
  // project all the points to the plane perpendicular to the axis
  pt(axis) = 0.0;
  line0(axis) = 0.0;
  line1(axis) = 0.0;

  return distanceFromLine(pt, line0, line1);
}

projectedIndices()

std::pair< unsigned int, unsigned int > geom_utils::projectedIndices

(

const unsigned int

axis

)

Get the indices of the plane perpendicular to the specified axis.

For example, if the axis is the y-axis (1), then this will return (0, 2), indicating that the coordinates of a general 3-D point once projected onto the x-z plane can be obtained with the 0 and 2 indices.

Parameters

[in]

axis

axis perpendicular to the projection plane

Returns

indices of coordinates on plane

Definition at line 145 of file GeometryUtils.C.

Referenced by minDistanceToPoints(), pointOnEdge(), projectedLineHalfSpace(), projectedUnitNormal(), and projectPoint().

{
  std::pair<unsigned int, unsigned int> indices;

  if (axis == 0)
  {
    indices.first = 1;
    indices.second = 2;
  }
  else if (axis == 1)
  {
    indices.first = 0;
    indices.second = 2;
  }
  else
  {
    indices.first = 0;
    indices.second = 1;
  }

  return indices;
}

projectedLineHalfSpace() [1/2]

Real geom_utils::projectedLineHalfSpace	(	Point	pt1,
		Point	pt2,
		Point	pt3,
		const unsigned int	axis
	)

Definition at line 68 of file GeometryUtils.C.

{
  // project points onto plane perpendicular to axis
  pt1(axis) = 0.0;
  pt2(axis) = 0.0;
  pt3(axis) = 0.0;

  const auto i = projectedIndices(axis);

  return (pt1(i.first) - pt3(i.first)) * (pt2(i.second) - pt3(i.second)) -
         (pt2(i.first) - pt3(i.first)) * (pt1(i.second) - pt3(i.second));
}

projectedLineHalfSpace() [2/2]

libMesh::Real geom_utils::projectedLineHalfSpace	(	libMesh::Point	pt1,
		libMesh::Point	pt2,
		libMesh::Point	pt3,
		const unsigned int	axis
	)

If positive, point is on the positive side of the half space (and vice versa).

Because a 3-D line does not have a negative or positive "side," you must provide the 'axis' perpendicular to the plane into which the point and line are first projected.

Parameters

[in]	pt1	point of interest
[in]	pt2	one end point of line
[in]	pt3	other end point of line
[in]	axis	axis perpendicular to plane onto which point and line are first projected

Returns

half space of line

Referenced by pointInPolygon().

projectedUnitNormal() [1/2]

libMesh::Point geom_utils::projectedUnitNormal	(	libMesh::Point	pt1,
		libMesh::Point	pt2,
		const unsigned int	axis
	)

Get the unit normal vector between two points (which are first projected onto the plane perpendicular to the 'axis'), such that the cross product of the unit normal with the line from pt1 to pt2 has a positive 'axis' component.

Parameters

[in]	pt1	first point for line
[in]	pt2	second point for line
[in]	axis	project points onto plane perpendicular to this axis

Returns

unit normal

projectedUnitNormal() [2/2]

Point geom_utils::projectedUnitNormal	(	Point	pt1,
		Point	pt2,
		const unsigned int	axis
	)

Definition at line 169 of file GeometryUtils.C.

{
  // project the points to the plane perpendicular to the axis
  pt1(axis) = 0.0;
  pt2(axis) = 0.0;

  const auto i = projectedIndices(axis);

  const Real dx = pt2(i.first) - pt1(i.first);
  const Real dy = pt2(i.second) - pt1(i.second);

  const Point normal = projectPoint(dy, -dx, axis);
  const Point gap_line = pt2 - pt1;

  const auto cross_product = gap_line.cross(normal);

  if (cross_product(axis) > 0)
    return normal.unit();
  else
    return projectPoint(-dy, dx, axis).unit();
}

projectPoint() [1/2]

Point geom_utils::projectPoint	(	const Real	x0,
		const Real	x1,
		const unsigned int	axis
	)

Definition at line 56 of file GeometryUtils.C.

{
  const auto i = projectedIndices(axis);
  Point point;
  point(i.first) = x0;
  point(i.second) = x1;
  point(axis) = 0.0;

  return point;
}

projectPoint() [2/2]

libMesh::Point geom_utils::projectPoint	(	const libMesh::Real	x0,
		const libMesh::Real	x1,
		const unsigned int	axis
	)

Given two coordinates, construct a point in the 2-D plane perpendicular to the specified axis.

Parameters

[in]	x0	first coordinate
[in]	x1	second coordinate
[in]	axis	axis perpendicular to the projection plane

Returns

point

Referenced by polygonCorners(), and projectedUnitNormal().

rotatePointAboutAxis() [1/2]

libMesh::Point geom_utils::rotatePointAboutAxis	(	const libMesh::Point &	p,
		const libMesh::Real	angle,
		const libMesh::Point &	axis
	)

Rotate point about an axis.

Parameters

[in]	p	point
[in]	angle	angle to rotate (radians)
[in]	axis	axis expressed as vector

Returns

rotated point

Referenced by TransformedPositions::initialize().

rotatePointAboutAxis() [2/2]

Point geom_utils::rotatePointAboutAxis	(	const Point &	p,
		const Real	angle,
		const Point &	axis
	)

Definition at line 232 of file GeometryUtils.C.

{
  const Real cos_theta = cos(angle);
  const Real sin_theta = sin(angle);

  Point pt;
  const Real xy = axis(0) * axis(1);
  const Real xz = axis(0) * axis(2);
  const Real yz = axis(1) * axis(2);

  const Point x_op(cos_theta + axis(0) * axis(0) * (1.0 - cos_theta),
                   xy * (1.0 - cos_theta) - axis(2) * sin_theta,
                   xz * (1.0 - cos_theta) + axis(1) * sin_theta);

  const Point y_op(xy * (1.0 - cos_theta) + axis(2) * sin_theta,
                   cos_theta + axis(1) * axis(1) * (1.0 - cos_theta),
                   yz * (1.0 - cos_theta) - axis(0) * sin_theta);

  const Point z_op(xz * (1.0 - cos_theta) - axis(1) * sin_theta,
                   yz * (1.0 - cos_theta) + axis(0) * sin_theta,
                   cos_theta + axis(2) * axis(2) * (1.0 - cos_theta));

  pt(0) = x_op * p;
  pt(1) = y_op * p;
  pt(2) = z_op * p;
  return pt;
}

segmentsIntersect()

bool geom_utils::segmentsIntersect	(	const Point &	p1,
		const Point &	p2,
		const Point &	p3,
		const Point &	p4
	)

Check if the line segment p1-p2 intersects with line segment p3-p4 (only working in 2D (x-y plane)).

Parameters

p1	First point of first line segment
p2	Second point of first line segment
p3	First point of second line segment
p4	Second point of second line segment

Returns

true if the two line segments intersect, false otherwise

Definition at line 295 of file GeometryUtils.C.

Referenced by BoundaryLayerUtils::generateOffsetPolyline().

{
  mooseAssert(
      MooseUtils::absoluteFuzzyEqual(p1(2), 0.0) && MooseUtils::absoluteFuzzyEqual(p2(2), 0.0) &&
          MooseUtils::absoluteFuzzyEqual(p3(2), 0.0) && MooseUtils::absoluteFuzzyEqual(p4(2), 0.0),
      "segmentsIntersect only works in 2D (x-y plane)");

  mooseAssert(MooseUtils::absoluteFuzzyGreaterThan((p1 - p2).norm(), 0.0) &&
                  MooseUtils::absoluteFuzzyGreaterThan((p3 - p4).norm(), 0.0),
              "Zero length segments are not allowed in segmentsIntersect");

  const Real a1 = p2(1) - p1(1);
  const Real b1 = p1(0) - p2(0);
  const Real c1 = p2(0) * p1(1) - p1(0) * p2(1);

  const Real a2 = p4(1) - p3(1);
  const Real b2 = p3(0) - p4(0);
  const Real c2 = p4(0) * p3(1) - p3(0) * p4(1);

  const Real denom = a1 * b2 - a2 * b1;
  Point intersection_pt;
  // Parallel case
  if (MooseUtils::absoluteFuzzyEqual(denom, 0.0))
  {
    if (arePointsColinear(p1, p2, p3))
    {
      // for colinear segments, we construct a "virtual" intersection point using weighted average
      // In that case, the virtual point will always lie outside of both segments unless they
      // overlap
      const Point p12 = (p1 + p2) / 2.0;
      const Point p34 = (p3 + p4) / 2.0;
      const Real dist12 = (p1 - p2).norm();
      const Real dist34 = (p3 - p4).norm();
      intersection_pt = p12 + (p34 - p12) * (dist12 / (dist12 + dist34));
    }
    else
      return false;
  }
  else
  {
    intersection_pt = Point((b1 * c2 - b2 * c1) / denom, (a2 * c1 - a1 * c2) / denom, 0.0);
  }

  const Real ratio_p1p2 = (intersection_pt - p1) * (p2 - p1) / ((p2 - p1).norm_sq());
  const Real ratio_p3p4 = (intersection_pt - p3) * (p4 - p3) / ((p4 - p3).norm_sq());

  if (MooseUtils::absoluteFuzzyGreaterEqual(ratio_p1p2, 0.0) &&
      MooseUtils::absoluteFuzzyLessEqual(ratio_p1p2, 1.0) &&
      MooseUtils::absoluteFuzzyGreaterEqual(ratio_p3p4, 0.0) &&
      MooseUtils::absoluteFuzzyLessEqual(ratio_p3p4, 1.0))
    return true;
  else
    return false;
}

solidAngle()

Real geom_utils::solidAngle	(	const Point &	point,
		const Point &	v0,
		const Point &	v1,
		const Point &	v2
	)

Compute the signed solid angle subtended by one oriented triangle at the query point.

Parameters

[in]	point	query point
[in]	v0	first triangle vertex
[in]	v1	second triangle vertex
[in]	v2	third triangle vertex

Returns

signed solid angle in steradians

Definition at line 419 of file GeometryUtils.C.

Referenced by TriangleManifold::containsBySolidAngle().

{
  const Point a = v0 - point;
  const Point b = v1 - point;
  const Point c = v2 - point;

  const Real la = a.norm();
  const Real lb = b.norm();
  const Real lc = c.norm();

  const Real numerator = a * (b.cross(c));
  const Real denominator = la * lb * lc + (a * b) * lc + (b * c) * la + (c * a) * lb;

  return 2.0 * std::atan2(numerator, denominator);
}

unitVector() [1/2]

Point geom_utils::unitVector	(	const Point &	pt,
		const std::string &	name
	)

Definition at line 28 of file GeometryUtils.C.

{
  if (isPointZero(pt))
    mooseError("'" + name + "' cannot have zero norm!");

  return pt.unit();
}

unitVector() [2/2]

libMesh::Point geom_utils::unitVector	(	const libMesh::Point &	pt,
		const std::string &	name
	)

Get the unit vector for a point parameter.

Parameters

[in]	pt	point
[in]	name	name of the parameter