!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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),intent(in) :: params
integer,intent(in) :: level
integer,intent(in) :: levelmax
integer,intent(in) :: icon(params%mpe)
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(in) :: materialn(0:nlsf)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(inout) :: pressure
double precision,intent(in) :: strain(nnode)
integer,intent(in) :: nnode
double precision,intent(in) :: lsf(params%mpe,nlsf)
integer,intent(in) :: nlsf
double precision,intent(in) :: r0,s0,t0,rst
integer,intent(out) :: icut
integer,intent(in) :: ileaves
double precision,intent(in) :: eviscosity

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

integer :: matel,matrule,i,j,k,err
double precision :: prod,vol_lsf0,eps,pressurep
double precision,allocatable :: vol_lsf(:)

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

eps=1.d-10
icut=0

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

matel=materialn(0)
matrule=params%matrule

do i=1,nlsf
   prod=lsf(1,i)
   do k=2,params%mpe
      if (prod*lsf(k,i).le.0.d0 .and. icut.eq.0) then
        icut=1
        allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in make_cut$')
        call get_mat_volume(params,nlsf,lsf,vol_lsf,vol_lsf0)
        if (vol_lsf0.gt.eps .or. level.gt.0) matrule=0
        if (matrule.eq.1 .or. matrule.eq.2) then                                ! Use majority or minority rule
          call set_elem_props(params,mat,matrule,nlsf,nnode,icon,level,        &
                             materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp,      &
                             pressure,strain,r0,s0,t0,rst,eviscosity)
        else                                                                    ! Use divFEM
          if (level.eq.levelmax) then
          call set_elem_props(params,mat,matrule,nlsf,nnode,icon,level,        &
                             materialn,vol_lsf,vol_lsf0,x,y,z,u,v,w,temp,      &
                             pressure,strain,r0,s0,t0,rst,eviscosity)
          else
            call divide_element(params,mat,level,levelmax,nlsf,nnode,icon,     &
                               ileaves,materialn,lsf,r0,s0,t0,rst,x,y,z,u,v,w, &
                               temp,pressure,strain,eviscosity)
          endif
        endif
        deallocate (vol_lsf)
      endif
   enddo
   if (prod.lt.0.d0 .and. icut.eq.0) then
     matel=materialn(i)
   endif
enddo

if (icut.eq.0) then
  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)
endif

end subroutine pressure_cut

!-------------------------------------------------------------------------------

subroutine get_mat_volume (params,nlsf,lsf,vol_lsf,vol_lsf0)

use definitions

implicit none

type(parameters),intent(in) :: params
integer,intent(in) :: nlsf
double precision,intent(in)  :: lsf(params%mpe,nlsf)
double precision,intent(out) :: vol_lsf(nlsf)
double precision,intent(out) :: vol_lsf0

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

integer :: i

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

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
endif

! this is a little fix that is necessary due to the finite precision with which
! we compute the positive volumes and could lead to a contributing volume being
! either marginally negative or greater than 1
do i=1,nlsf
  vol_lsf(i)=max(vol_lsf(i),0.d0)
  vol_lsf(i)=min(vol_lsf(i),1.d0)
enddo
vol_lsf0=1.d0-sum(vol_lsf)

end subroutine get_mat_volume

!-------------------------------------------------------------------------------

subroutine set_elem_props (params,mat,matrule,nlsf,nnode,icon,level,materialn, &
                          vol_lsf,vol_lsf0,x,y,z,u,v,w,temp,pressure,strain,r0,&
                          s0,t0,rst,eviscosity)

use definitions

implicit none

type(parameters),intent(in) :: params
type(material),intent(in)   :: mat(0:params%nmat)
integer,intent(in) :: matrule
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: level
integer,intent(in) :: materialn(0:nlsf)
double precision,intent(in) :: vol_lsf(nlsf)
double precision,intent(in) :: vol_lsf0
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: r0,s0,t0,rst
double precision,intent(in) :: eviscosity
double precision,intent(out) :: pressure

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

integer :: matel,i
double precision :: eps,viscosity,penal,expon,pressurep
double precision :: plasticity_parameters(9)
character(len=8) :: plasticity_type

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

eps=1.d-10

select case (matrule)
case (1)                                                                        ! Assign material properties of volumetric majority material
  matel=materialn(sum(maxloc(vol_lsf)))
  viscosity=mat(matel)%viscosity
  penal=mat(matel)%penalty
  expon=mat(matel)%expon
  plasticity_type=mat(matel)%plasticity_type
  plasticity_parameters=mat(matel)%plasticity_parameters

case(2)                                                                         ! Assign material properties of volumetric minority material
  matel=materialn(sum(minloc(vol_lsf)))
  viscosity=mat(matel)%viscosity
  penal=mat(matel)%penalty
  expon=mat(matel)%expon
  plasticity_type=mat(matel)%plasticity_type
  plasticity_parameters=mat(matel)%plasticity_parameters

case default
  if (matrule.ne.0) write(*,*) 'Invalid matrule value, using divFEM'
  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
  matel=materialn(0)
  viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
  penal=penal+vol_lsf0*mat(matel)%penalty
  expon=expon+vol_lsf0*mat(matel)%expon
  matel=materialn(sum(maxloc(vol_lsf)))
  plasticity_type=mat(matel)%plasticity_type
  plasticity_parameters=mat(matel)%plasticity_parameters

end select

if (maxval(vol_lsf).lt.eps) then
  matel=materialn(0)
  plasticity_type=mat(matel)%plasticity_type
  plasticity_parameters=mat(matel)%plasticity_parameters
endif

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)

end subroutine set_elem_props

!-------------------------------------------------------------------------------

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

use definitions

implicit none

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

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

integer :: err,i,j,k,ii,jj,kk,levelp,jcut
double precision,allocatable :: lsfp(:,:)
double precision :: r(params%mpe),s(params%mpe),t(params%mpe),h(params%mpe)
double precision :: r0p,s0p,t0p,rstp

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

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)

end subroutine divide_element

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------