5. A point and the body

probe now inspects the point q = (0.8, 0.6, 0.8) of chapter 1: is it inside the tetrahedron, how far from each face, where does it fall on the slanted one? Then it checks points against lines and circles.

program probe
!< Tutorial, chapter 5: a point and the tetrahedron, distances, projections, collinear and concyclic points.
use penf, only : R8P
use vecfor, only : ex, ey, ez, vector
implicit none
type(vector) :: o          ! the origin
type(vector) :: p(4)       ! the vertices
integer      :: face(3,4)  ! the vertices of each face, counter-clockwise seen from outside
type(vector) :: q          ! the probe
type(vector) :: c          ! the centroid
type(vector) :: foot       ! the projection of the probe onto the slanted face
type(vector) :: d          ! a vectorial distance
type(vector) :: m          ! a point on an edge
type(vector) :: a          ! a point on the circle through the slanted face
real(R8P)    :: dist(4)    ! signed distances to the faces
integer      :: f          ! counter

p(1) = o ; p(2) = ex ; p(3) = ey ; p(4) = ez
face = reshape([1,3,2, 1,2,4, 1,4,3, 2,3,4], [3,4])
q = 0.8_R8P * ex + 0.6_R8P * ey + 0.8_R8P * ez
c = (p(1) + p(2) + p(3) + p(4)) / 4
do f = 1, 4
   dist(f) = q%distance_to_plane(pt1=p(face(1,f)), pt2=p(face(2,f)), pt3=p(face(3,f)))
enddo
print '(A,4F8.4)', 'q to the faces       ', dist
print '(A,L1)',    'q is inside:         ', all(dist <= 0._R8P)
do f = 1, 4
   dist(f) = c%distance_to_plane(pt1=p(face(1,f)), pt2=p(face(2,f)), pt3=p(face(3,f)))
enddo
print '(A,4F8.4)', 'centroid to the faces', dist
print '(A,L1)',    'centroid is inside:  ', all(dist <= 0._R8P)
d = q%distance_vectorial_to_plane(pt1=p(2), pt2=p(3), pt3=p(4))
foot = q%projection_onto_plane(pt1=p(2), pt2=p(3), pt3=p(4))
print '(A,3F8.4)', 'from the face to q   ', d%x, d%y, d%z
print '(A,3F8.4)', 'foot of q            ', foot%x, foot%y, foot%z
print '(A,3F8.4)', 'foot + distance      ', foot%x + d%x, foot%y + d%y, foot%z + d%z
print '(A,F8.4)', 'q to the line p2 p3  ', q%distance_to_line(pt1=p(2), pt2=p(3))
print '(A,F8.4)', 'q to the line p1 p4  ', q%distance_to_line(pt1=p(1), pt2=p(4))
m = (p(2) + p(3)) / 2                          ! the midpoint: exact in floating point
print '(A,L1)',    'midpoint on the line p2 p3:   ', m%is_collinear(pt1=p(2), pt2=p(3))
m = 0.3_R8P * p(2) + 0.7_R8P * p(3)            ! 0.3 and 0.7 are not exact in binary
print '(A,ES8.1)', 'distance of m from the line:  ', m%distance_to_line(pt1=p(2), pt2=p(3))
print '(A,L1)',    'm on the line p2 p3:          ', m%is_collinear(pt1=p(2), pt2=p(3))
print '(A,L1)',    'with a tolerance of 1e-12:    ', m%is_collinear(pt1=p(2), pt2=p(3), tolerance=1.e-12_R8P)
m = m + 1.e-9_R8P * ez
print '(A,L1)',    '1e-9 off, tolerance of 1e-12: ', m%is_collinear(pt1=p(2), pt2=p(3), tolerance=1.e-12_R8P)
a = 2 * (p(2) + p(3) + p(4)) / 3 - p(2)        ! opposite to p2, on the circle through p2, p3, p4
print '(A,L1)', 'a, p4, p2, p3 on a circle: ', a%is_concyclic(pt1=p(4), pt2=p(2), pt3=p(3))
print '(A,L1)', 'a, p2, p4, p3 on a circle: ', a%is_concyclic(pt1=p(2), pt2=p(4), pt3=p(3))
endprogram probe

$ probe
q to the faces        -0.8000 -0.6000 -0.8000  0.6928
q is inside:         F
centroid to the faces -0.2500 -0.2500 -0.2500 -0.1443
centroid is inside:  T
from the face to q     0.4000  0.4000  0.4000
foot of q              0.4000  0.2000  0.4000
foot + distance        0.8000  0.6000  0.8000
q to the line p2 p3    0.8485
q to the line p1 p4    1.0000
midpoint on the line p2 p3:   T
distance of m from the line:   2.0E-17
m on the line p2 p3:          F
with a tolerance of 1e-12:    T
1e-9 off, tolerance of 1e-12: F
a, p4, p2, p3 on a circle: T
a, p2, p4, p3 on a circle: F

Inside or outside

q%distance_to_plane(pt1=, pt2=, pt3=) is the signed distance of q from the plane through the three points: positive on the side the normal face_normal3(pt1, pt2, pt3) points to, negative on the other. With the faces listed as in chapter 4, every normal points outward, so a point is inside when all its distances are negative or zero. q is 0.6928 ($0.4\sqrt{3}$) outside the slanted face; the centroid is inside, a quarter from each face on the coordinate planes.

The foot of a point

q%distance_vectorial_to_plane(...) is the vector from the plane to q, along the normal: its length is the distance, its direction tells the side. q%projection_onto_plane(...) is the foot of q, the point of the plane nearest to it, (0.4, 0.2, 0.4); foot plus vectorial distance gives q back.

Lines

q%distance_to_line(pt1=, pt2=) is the distance of q from the line through two points, computed as $|(q-p_1)\times (q-p_2)|/|p_2-p_1|$: the line is infinite, not the segment between the points.

Collinear points and tolerances

m%is_collinear(pt1=, pt2=) asks whether m lies on the line through two points. Without tolerance the test is exact: the distance must be zero. The midpoint of p2 p3 passes, because its coordinates are exact in binary; the point 0.3 p2 + 0.7 p3 lies on the line in exact arithmetic, but 0.3 and 0.7 are not exact in binary, and it misses by $2\times {10}^{-17}$. With tolerance=1.e-12_R8P it passes, and a point 1e-9 off the line still fails. Pass a tolerance whenever the points are computed: scale it to the size of your geometry.

Concyclic points

a%is_concyclic(pt1=, pt2=, pt3=) asks whether four points lie on one circle, in this order around it: VecFor tests Ptolemy's theorem, $|a{p}_1||p_2{p}_3|+|p_1{p}_2||p_{3}a|=|a{p}_2||p_1{p}_3|$, which holds for the points of a circle taken in their order along it. a, the point of the circle through the slanted face opposite to p2, lies between p3 and p4: the order a, p4, p2, p3 passes, a, p2, p4, p3 does not. Its default tolerance is the smallest positive number, so computed points need a tolerance here too (see Are four points on one circle?).

What you learned

distance_to_plane (signed), distance_vectorial_to_plane, projection_onto_plane, distance_to_line, is_collinear and is_concyclic with their tolerance. Reference: Geometry.

Next: 6. Components: a slip wall.