!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              PRESSURE_CUT    Nov. 2006    I                                  |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

recursive subroutine pressure_cut (params,level,levelmax,icon,   &
                                   x,y,z,mat,materialn,                   &
                                   u,v,w,temp,pressure,strain,nnode,           &
                                   lsf,nlsf,r0,s0,t0,rst,icut,ileaves,eviscosity)
                                   

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this is an intermediary routine between find_pressure and make_pressure
! in the same way as make_cut is between build_system and make_matrix

! level is current level in cut cell algortihm
! it varies between 0 and levelmax
! icon is connectivity array for the current element
! x,y,z are the global coordinate arrays of length nnode
! mat is the material matrix for the nmat materials
! materialn contains the material number associated to each lsf
! u,v,w is the velocity array (obtained from previous time step
! or at least containing the proper velocity at the fixed dofs)
! temp, pressure and strain are temperature, pressure and strain
! nnode is number of nodes
! lsf is global array of level set functions defining the surfaces
! nlsf is number of lsfs
! r0,s0,t0 are the bottom left back corner of the part of the element
!   we are now integrating (in local r,s,t coordinates)
! rst is the size of the part of the element we are integrating
! icut is returned (0 if homogeneous element - 1 if cut element)

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|
    
use definitions

implicit none

type (parameters) params
integer level
integer levelmax
integer icon(params%mpe)
double precision x(nnode),y(nnode),z(nnode)
type (material) mat(0:params%nmat)
integer materialn(0:nlsf)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
integer nnode
double precision lsf(params%mpe,nlsf)
integer nlsf
double precision r0,s0,t0,rst
integer icut
integer ileaves
double precision eviscosity

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer err,i,j,k,ii,jj,kk,levelp,matel,jcut
double precision r(params%mpe),s(params%mpe),t(params%mpe)
double precision h(params%mpe),vol_lsf0
double precision r0p,s0p,t0p,rstp,dt
double precision viscosity,density,penal,expon,diffusivity,heat
double precision pressurep,prod,volmax
double precision,dimension(:,:),allocatable::lsfp
double precision,dimension(:),allocatable::vol_lsf
character (len=8) plasticity_type
double precision plasticity_parameters(9)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

if (level.eq.0) then
   pressure=0.d0
endif

matel=materialn(0)

do i=1,nlsf
   prod=lsf(1,i)
   do k=2,params%mpe
      if (prod*lsf(k,i).le.0.d0) goto 222
   enddo
   if (prod.lt.0.d0) then
      matel=materialn(i)
   endif
enddo

call make_pressure (params,icon,x,y,z, &
                    mat(matel)%viscosity, &
                    mat(matel)%penalty,mat(matel)%expon, &
                    u,v,w,temp,pressurep,strain,nnode, &
                    r0,s0,t0,rst,mat(matel)%plasticity_type,&
                    mat(matel)%plasticity_parameters,eviscosity)

pressure=pressure+pressurep/(8.d0**level)


icut=0

return

222   continue

icut=1

! when we get to the bottom of the division
! we use an approximate algorithm

if (level.eq.levelmax) then
   allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in make_cut$')
   do i=1,nlsf
      call compute_positive_volume (lsf(1,i),vol_lsf(i),params%levelapprox)
   enddo
   vol_lsf=1.d0-vol_lsf

   if (.not.params%excl_vol) then
      do i=1,nlsf-1
         vol_lsf(i)=vol_lsf(i)-vol_lsf(i+1)
      enddo
   end if

    do i=1,nlsf
       vol_lsf(i)=max(vol_lsf(i),0.d0)
       vol_lsf(i)=min(vol_lsf(i),1.d0)
    enddo

   viscosity=0.d0
   penal=0.d0
   expon=0.d0
   do i=1,nlsf
      matel=materialn(i)
      viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
      penal=penal+vol_lsf(i)*mat(matel)%penalty
      expon=expon+vol_lsf(i)*mat(matel)%expon
   enddo
   vol_lsf0=1.d0-sum(vol_lsf)
   matel=materialn(0)
   viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
   penal=penal+vol_lsf0*mat(matel)%penalty
   expon=expon+vol_lsf0*mat(matel)%expon

   volmax=0.d0
   matel=params%materialn(0)
   plasticity_type=mat(matel)%plasticity_type
   plasticity_parameters=mat(matel)%plasticity_parameters
   do i=1,nlsf
      if (vol_lsf(i).gt.volmax) then
         volmax=vol_lsf(i)
         matel=params%materialn(i)
         plasticity_type=mat(matel)%plasticity_type
         plasticity_parameters=mat(matel)%plasticity_parameters
      endif
   enddo


   call make_pressure (params,icon,x,y,z, &
                       viscosity,penal,expon, &
                       u,v,w,temp,pressurep,strain,nnode, &
                       r0,s0,t0,rst,plasticity_type,plasticity_parameters,eviscosity)
   pressure=pressure+pressurep/(8.d0**level)
   deallocate (vol_lsf)
   return
endif

! If we are not at the bottom level, we keep dividing

allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in make_cut$')

do kk=1,2
   t(1:4)=-1.d0+float(kk-1)
   t(5:8)=float(kk-1)
   do jj=1,2
      s(1:2)=-1.d0+float(jj-1)
      s(3:4)=float(jj-1)
      s(5:6)=-1.d0+float(jj-1)
      s(7:8)=float(jj-1)
      do ii=1,2
         r(1)=-1.d0+float(ii-1)
         r(2)=float(ii-1)
         r(3)=-1.d0+float(ii-1)
         r(4)=float(ii-1)
         r(5)=-1.d0+float(ii-1)
         r(6)=float(ii-1)
         r(7)=-1.d0+float(ii-1)
         r(8)=float(ii-1)
         do k=1,8
            h(1)=(1.d0-r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
            h(2)=(1.d0+r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
            h(3)=(1.d0-r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
            h(4)=(1.d0+r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
            h(5)=(1.d0-r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
            h(6)=(1.d0+r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
            h(7)=(1.d0-r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
            h(8)=(1.d0+r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
            lsfp(k,:)=0.d0
            do i=1,nlsf
               do j=1,8
                  lsfp(k,i)=lsfp(k,i)+h(j)*lsf(j,i)
               enddo
            enddo
         enddo
         r0p=r0+rst*(r(1)+1.d0)/2.d0
         s0p=s0+rst*(s(1)+1.d0)/2.d0
         t0p=t0+rst*(t(1)+1.d0)/2.d0
         rstp=rst/2.d0
         levelp=level+1
         call pressure_cut (params,levelp,levelmax, &
                            icon,x,y,z, &
                            mat,materialn, &
                            u,v,w,temp,pressure,strain,nnode, &
                            lsfp,nlsf,r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity)
      enddo
   enddo
enddo

deallocate (lsfp)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
