!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
