4. Faces, area, volume
A solver sees a body through its faces: their normals, their areas, and the volume they enclose. probe computes all of them for its tetrahedron.
program probe
!< Tutorial, chapter 4: the faces, their normals and areas, the surface and the volume of the tetrahedron.
use penf, only : R8P
use vecfor, only : ex, ey, ez, face_normal3, vector
implicit none
type(vector) :: o ! the origin
type(vector) :: p(4) ! the vertices
integer :: face(3,4) ! the vertices of each face
type(vector) :: n ! the normal of a face
type(vector) :: c ! the centroid of a face
real(R8P) :: area ! the surface of the tetrahedron
real(R8P) :: volume ! the volume of the tetrahedron
integer :: f ! counter
p(1) = o ; p(2) = ex ; p(3) = ey ; p(4) = ez
face(:,1) = [1, 3, 2] ! counter-clockwise seen from outside: the normal points outward
face(:,2) = [1, 2, 4]
face(:,3) = [1, 4, 3]
face(:,4) = [2, 3, 4]
area = 0._R8P
volume = 0._R8P
do f = 1, 4
n = face_normal3(pt1=p(face(1,f)), pt2=p(face(2,f)), pt3=p(face(3,f))) ! its length is the face area
c = (p(face(1,f)) + p(face(2,f)) + p(face(3,f))) / 3
area = area + n%normL2()
volume = volume + (c .dot. n) / 3 ! divergence theorem
print '(A,I0,A,3F8.4,A,F8.4)', 'face ', f, ': normal', n%x, n%y, n%z, ', area', n%normL2()
enddo
print '(A,F8.4)', 'surface ', area
print '(A,F8.4,A,F8.4)', 'volume ', volume, ', 1/6 =', 1._R8P / 6
n = face_normal3(pt1=p(2), pt2=p(3), pt3=p(4), norm='y')
call n%printf(prefix='unit normal of face 4 = ')
n = face_normal3(pt1=p(2), pt2=p(4), pt3=p(3), norm='y')
call n%printf(prefix='p3 and p4 swapped = ')
endprogram probe$ probe
face 1: normal 0.0000 0.0000 -0.5000, area 0.5000
face 2: normal 0.0000 -0.5000 0.0000, area 0.5000
face 3: normal -0.5000 0.0000 0.0000, area 0.5000
face 4: normal 0.5000 0.5000 0.5000, area 0.8660
surface 2.3660
volume 0.1667, 1/6 = 0.1667
unit normal of face 4 = +0.5773502691896258, +0.5773502691896258, +0.5773502691896258
p3 and p4 swapped = -0.5773502691896258, -0.5773502691896258, -0.5773502691896258The normal of a face
face_normal3(pt1=, pt2=, pt3=) is the normal of the triangle through the three points: norm='y' (any character does) the normal is a unit vector; swapping two points reverses it, as the last line shows.
face_normal4(pt1=, pt2=, pt3=, pt4=) does the same for a quadrilateral, from its diagonals:
Outward normals
face lists the vertices of each face counter-clockwise seen from outside, so every normal points out of the body: (0, 0, -0.5) for the face on the plane z = 0, (0.5, 0.5, 0.5) for the slanted one. A mesh generator usually guarantees such an order; when it does not, compare each normal with the vector from the centroid of the body to the centroid of the face, and swap two vertices where their dot product is negative.
Surface and volume
The surface is the sum of the lengths of the normals. The volume comes from the divergence theorem: for a closed surface of flat faces, probe finds 0.1667, the volume of the corner tetrahedron, 1/6. The same sum gives the volume of any closed mesh of triangles, and a negative volume reveals faces listed clockwise.
What you learned
face_normal3 and face_normal4, with and without norm; areas from normals; outward normals from the order of the vertices; the volume of a body by the divergence theorem. Reference: Geometry.
Next: 5. A point and the body.