!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! POSITIVE_VOLUME Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine compute_positive_volume (phi,vol,levelmax)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! finds the volume of a cube that is defined by the positive value of a lsf
! known at the nodes of the cube
! phi are the lsf values
! vol is the returned volume
! levelmax is the max level for accuracy (power of 2)
! this routine has been tested and provides accurate solutions in
! all ranges of phi with levelmax=3
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
double precision vol,phi(8)
integer levelmax
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
call volume_lsf (phi,vol,1,levelmax)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! VOLUME_LSF Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
recursive subroutine volume_lsf (phi,volp,level,levelmax)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! used by positive_volume
! internal routine
! should not be changed
! it computes the volume based on Marthijn's simple formula
! the algorithm has been improved by allowing
! multi-level estimates up to levelmax by simple octree
! division before the formula is used
! phi are the lsf values at the 8 nodes
! volp is volume where lsf is negative
! level is current level in octree division
! levelmax is maximum level of octree division
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
double precision phi(8)
double precision volp
integer level
integer levelmax
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision phip(8),rat
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
if (level.eq.1) then
volp=0.d0
endif
if (level.eq.levelmax) then
rat=(sum(phi)/sum(abs(phi))+1.d0)/2.d0
volp=volp+rat*(2.d0**(-level+1))**3
return
endif
if (phi(1).gt.0.d0 .and. phi(2).gt.0.d0 .and. &
phi(3).gt.0.d0 .and. phi(4).gt.0.d0 .and. &
phi(5).gt.0.d0 .and. phi(6).gt.0.d0 .and. &
phi(7).gt.0.d0 .and. phi(8).gt.0.d0) then
volp=volp+(2.d0**(-level+1))**3
return
endif
if (phi(1).lt.0.d0 .and. phi(2).lt.0.d0 .and. &
phi(3).lt.0.d0 .and. phi(4).lt.0.d0 .and. &
phi(5).lt.0.d0 .and. phi(6).lt.0.d0 .and. &
phi(7).lt.0.d0 .and. phi(8).lt.0.d0) then
return
endif
phip(1)=phi(1)
phip(2)=(phi(1)+phi(2))/2.d0
phip(3)=(phi(1)+phi(3))/2.d0
phip(4)=(phi(1)+phi(2)+phi(3)+phi(4))/4.d0
phip(5)=(phi(1)+phi(5))/2.d0
phip(6)=(phi(1)+phi(2)+phi(5)+phi(6))/4.d0
phip(7)=(phi(1)+phi(3)+phi(5)+phi(7))/4.d0
phip(8)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(2))/2.d0
phip(2)=phi(2)
phip(3)=(phi(1)+phi(2)+phi(3)+phi(4))/4.d0
phip(4)=(phi(2)+phi(4))/2.d0
phip(5)=(phi(1)+phi(2)+phi(5)+phi(6))/4.d0
phip(6)=(phi(2)+phi(6))/2.d0
phip(7)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(8)=(phi(2)+phi(4)+phi(6)+phi(8))/4.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(3))/2.d0
phip(2)=(phi(1)+phi(2)+phi(3)+phi(4))/4.d0
phip(3)=phi(3)
phip(4)=(phi(3)+phi(4))/2.d0
phip(5)=(phi(1)+phi(3)+phi(5)+phi(7))/4.d0
phip(6)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(7)=(phi(3)+phi(7))/2.d0
phip(8)=(phi(3)+phi(4)+phi(7)+phi(8))/4.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(2)+phi(3)+phi(4))/4.d0
phip(2)=(phi(2)+phi(4))/2.d0
phip(3)=(phi(3)+phi(4))/2.d0
phip(4)=phi(4)
phip(5)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(6)=(phi(2)+phi(4)+phi(6)+phi(8))/4.d0
phip(7)=(phi(3)+phi(4)+phi(7)+phi(8))/4.d0
phip(8)=(phi(4)+phi(8))/2.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(5))/2.d0
phip(2)=(phi(1)+phi(2)+phi(5)+phi(6))/4.d0
phip(3)=(phi(1)+phi(3)+phi(5)+phi(7))/4.d0
phip(4)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(5)=phi(5)
phip(6)=(phi(5)+phi(6))/2.d0
phip(7)=(phi(5)+phi(7))/2.d0
phip(8)=(phi(5)+phi(6)+phi(7)+phi(8))/4.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(2)+phi(5)+phi(6))/4.d0
phip(2)=(phi(2)+phi(6))/2.d0
phip(3)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(4)=(phi(2)+phi(4)+phi(6)+phi(8))/4.d0
phip(5)=(phi(5)+phi(6))/2.d0
phip(6)=phi(6)
phip(7)=(phi(5)+phi(6)+phi(7)+phi(8))/4.d0
phip(8)=(phi(6)+phi(8))/2.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(3)+phi(5)+phi(7))/4.d0
phip(2)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(3)=(phi(3)+phi(7))/2.d0
phip(4)=(phi(3)+phi(4)+phi(7)+phi(8))/4.d0
phip(5)=(phi(5)+phi(7))/2.d0
phip(6)=(phi(5)+phi(6)+phi(7)+phi(8))/4.d0
phip(7)=phi(7)
phip(8)=(phi(7)+phi(8))/2.d0
call volume_lsf (phip,volp,level+1,levelmax)
phip(1)=(phi(1)+phi(2)+phi(3)+phi(4)+phi(5)+phi(6)+phi(7)+phi(8))/8.d0
phip(2)=(phi(2)+phi(4)+phi(6)+phi(8))/4.d0
phip(3)=(phi(3)+phi(4)+phi(7)+phi(8))/4.d0
phip(4)=(phi(4)+phi(8))/2.d0
phip(5)=(phi(5)+phi(6)+phi(7)+phi(8))/4.d0
phip(6)=(phi(6)+phi(8))/2.d0
phip(7)=(phi(7)+phi(8))/2.d0
phip(8)=phi(8)
call volume_lsf (phip,volp,level+1,levelmax)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|