SUBROUTINE hexcalc(elements,nodes,vol,vtotal,centroids,nelms,nnodes)
    IMPLICIT NONE
    REAL*8,INTENT(IN)::nodes(nnodes,3)
    INTEGER*8,INTENT(IN)::elements(:,:)
    REAL*8,INTENT(INOUT)::vol(:,:),vtotal(nelms),centroids(:,:)
    INTEGER,INTENT(IN)::nelms,nnodes
    REAL*8::c,p,vrat
    REAL*8,ALLOCATABLE::r(:),s(:),t(:),sgn(:,:),dndeta(:,:,:),dxdeta(:,:),x(:,:),xicent(:)
    INTEGER*8::a,i,j,k,n

    ALLOCATE(r(8))
    ALLOCATE(s(8))
    ALLOCATE(t(8))
    ALLOCATE(sgn(8,3))
    ALLOCATE(dndeta(8,8,3))
    ALLOCATE(dxdeta(3,3))
    ALLOCATE(x(8,3))
    ALLOCATE(xicent(3))
    
    sgn = RESHAPE((/ -1, 1, 1, -1, -1, 1, 1, -1, &
        -1,-1,1,1,-1,-1,1,1, &
        -1,-1,-1,-1,1,1,1,1 /),(/ 8,3 /))

    c = 1.0/8.0
    p = 1.0/SQRT(3.0)
    DO a=1,8
        r(a) = sgn(a,1)*p
        s(a) = sgn(a,2)*p
        t(a) = sgn(a,3)*p
        DO i=1,8
            dndeta(a,i,1) = c*sgn(i,1)*(1.0+sgn(i,2)*s(a))*(1.0+sgn(i,3)*t(a))
            dndeta(a,i,2) = c*sgn(i,2)*(1.0+sgn(i,1)*r(a))*(1.0+sgn(i,3)*t(a))
            dndeta(a,i,3) = c*sgn(i,3)*(1.0+sgn(i,1)*r(a))*(1.0+sgn(i,2)*s(a))
        END DO
    END DO

    DO n = 1,nelms
        DO i = 1,8
            x(i,:) = nodes(elements(n,i),:)
        END DO
        DO a=1,8
            DO i=1,3
                DO j=1,3
                    dxdeta(i,j) = 0.0
                    DO k=1,8
                        dxdeta(i,j) = dxdeta(i,j) + x(k,i)*dndeta(a,k,j)
                    END DO
                END DO
            END DO
            vol(n,a) = dxdeta(1,1)*(dxdeta(2,2)*dxdeta(3,3)-dxdeta(2,3)*dxdeta(3,2))- &
            dxdeta(1,2)*(dxdeta(2,1)*dxdeta(3,3)-dxdeta(2,3)*dxdeta(3,1))+ &
            dxdeta(1,3)*(dxdeta(2,1)*dxdeta(3,2)-dxdeta(2,2)*dxdeta(3,1))
        END DO
        vtotal(n) = SUM(vol(n,:))
        xicent = (/ 0.0, 0.0, 0.0 /)
        DO i= 1,8
            vrat = vol(n,i)/vtotal(n)
            xicent(1) = xicent(1) + r(i)*vrat
            xicent(2) = xicent(2) + s(i)*vrat
            xicent(3) = xicent(3) + t(i)*vrat
        END DO
        DO j= 1,3
            DO i= 1,8
                centroids(n,j) = centroids(n,j) + x(i,j)*c*(1.0+sgn(i,1)*xicent(1))*&
                    (1.0+sgn(i,2)*xicent(2))*(1.0+sgn(i,3)*xicent(3))
            END DO
        END DO
    END DO
    DEALLOCATE(r,s,t,sgn,dndeta,dxdeta,x,xicent)
    RETURN
END SUBROUTINE hexcalc