Compute the volume of a hexahedron, given a list of points
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in), | dimension(1:3, 1:8) | :: | p_list |
function vol_hexahedron(p_list)
!< Compute the volume of a hexahedron, given a list of points
! The points are arranged in a specific order. For the
! element i,j,k, the order of nodes required are:
! i,j,k
! i+1, j, k
! i+1, j+1, k
! i, j+1, k
! i, j, k+1
! i+1, j, k+1
! i+1, j+1, k+1
! i, j+1, k+1
!
! The hexahedron is to be split into 5 tetrahedrons.
! Source: Split hex into 5 tetrahedron:
! No assumptions about planarity seem to be made. All cuts
! were made with a plane containing only 3 vertices at a time.
! https://ieeexplore.ieee.org/ieee_pilot/articles/06/ttg2009061587/assets/img/article_1/fig_6/large.gif
!
! The indices of the 5 split tetrahedra can be visualised from
! the above link. But since the volume of each tetrahedron
! depends on the determinant calculated, it is IMPERATIVE to
! ensure that a "correct" order is followed for the 4 points.
!
! The logic to get the "correct" order is explained as
! follows (Refer wiki article on parallelepiped):
! The determinant is taken of a matrix of pi - p4, i = 1, 2, 3.
! Graphically it denotes the sides with p4 as common vertex,
! with direction outward from p4, i.e., directed from
! p4 to pi, i = 1, 2, 3
! Hence, if you ensure that cross(p1-p4, p2-p4) is along
! p3-p4, then the determinant will be positive.
!
! From the above link, a set of 5 tetrahedra was obtained.
! Each tetrahedra has 4 points, and in the function calls
! below, care was taken to ensure that the order is observed
! while passing parameters into the vol_tetrahedron function
!-----------------------------------------------------------
implicit none
real(wp), dimension(1:3, 1:8), intent(in) :: p_list
real(wp) :: vol_hexahedron
real(wp) :: vol_hexahedron1
vol_hexahedron1 = 0.
vol_hexahedron1 = vol_hexahedron1 + &
vol_tetrahedron(p_list(:,1), p_list(:,5), &
p_list(:,8), p_list(:,6))
vol_hexahedron1 = vol_hexahedron1 + &
vol_tetrahedron(p_list(:,7), p_list(:,8), &
p_list(:,6), p_list(:,3))
vol_hexahedron1 = vol_hexahedron1 + &
vol_tetrahedron(p_list(:,8), p_list(:,4), &
p_list(:,1), p_list(:,3))
vol_hexahedron1 = vol_hexahedron1 + &
vol_tetrahedron(p_list(:,6), p_list(:,1), &
p_list(:,3), p_list(:,8))
vol_hexahedron1 = vol_hexahedron1 + &
vol_tetrahedron(p_list(:,1), p_list(:,2), &
p_list(:,6), p_list(:,3))
vol_hexahedron = 0.
vol_hexahedron = vol_hexahedron + &
vol_tetrahedron(p_list(:,2), p_list(:,6), &
p_list(:,5), p_list(:,7))
vol_hexahedron = vol_hexahedron + &
vol_tetrahedron(p_list(:,8), p_list(:,5), &
p_list(:,7), p_list(:,4))
vol_hexahedron = vol_hexahedron + &
vol_tetrahedron(p_list(:,5), p_list(:,1), &
p_list(:,2), p_list(:,4))
vol_hexahedron = vol_hexahedron + &
vol_tetrahedron(p_list(:,7), p_list(:,2), &
p_list(:,4), p_list(:,5))
vol_hexahedron = vol_hexahedron + &
vol_tetrahedron(p_list(:,2), p_list(:,3), &
p_list(:,7), p_list(:,4))
vol_hexahedron = max(vol_hexahedron,vol_hexahedron1)
end function vol_hexahedron