I really like the idea of using the shape functions. It will take some algebra, but you should be able to do it easily by hand with the element shape functions. For instance, the 8 shape functions are N(1), N(2)...N(8).
x=sum(X(i)*N(i))=X(1)*N(1)+....X(8)*N(8),..X(i) are the
x coordinates of the nodes.
y=sum(Y(i)*N(i)), z=sum(Z(i)*N(i)). The shape functions N(i) are a function of the parent element coordinates, xi, eta, zeta.
Volume=integral(1*dx*dy*dz)=integral(J*dxi*deta*dzeta), where J is the Jacobian, the determinant of a 3 x 3 matrix, first row are the derivatives dx/dxi, dy/dxi, and dz/dxi, second row are the derivates with respect to eta, third row zeta. The derivatives are easy to compute:
dx/dxi=sum(X(i)*dN(i)/dxi)=X(1)*dN(1)/dxi+X(2)*dN(2)/dxi+...X(8)*dN(8)/dxi, same for the other 8 derivatives.
Once you compute those derivatives, compute the Jacobian, the integral(J*dxi*deta*dzeta) you can do by hand I think. Even if not, you may be able to write a small VBA in Excel to compute the integral numerically.