!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              MAKE_CUT    Nov. 2006                                           |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
                               u,v,w,temp,pressure,strain,e2dp,is_plastic,     &
                               nnode,f,lsf,nlsf,r0,s0,t0,rst,icut,ileaves,     &
                               eviscosity,vbounded,yield_ratio,frict_angle,    &
                               params,threadinfo,weightel,matnum)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this subroutine is a intermediary routine between build_system and make_matrix 
! to take into account the complex geometry of cut cells

! if we are in a non cut cell, make-matrix is called
! if we are in a cut cell but at a level that it smaller than levelmax, the
! cell is further cut and make_cut is recursively called
! if we are in a cut cell and level is equal to levelmax, we call make_matrix
! with material properties that have been interpolated from the
! various material properties contribnuting to the cut cell

! level       : current level in cut cell algorithm. It varies between 0 and levelmax.
! levelapprox : used to improve the postitive volutme calculation by further division
! mpe         : number of nodes per element (8)
! ndof        : number of degrees of freedom per node (3)
! ael         : computed finite element matrix
! bel         : computed rhs vector
! icon        : connectivity array for the current element
! xg,yg,zg    : global coordinate arrays of length nnode
! penalty     : penalty factor used to impose the boundary conditions
! tempscale   : temperature scaling parameter
! kfix        : bc array of length ndof*nnode (kfix=1 means the dof is fixed to the value stored in velo)
! mat         : material matrix for the nmat materials
! materialn   : contains the material number associated to each lsf
! dt          : time step length (only needed for the temperature calculations)
! u,v,w       : velocity array (obtained from previous time step or at least containing the proper velocity at the fixed dofs)
! nnode       : number of nodes
! f           : global rhs vector
! lsf         : global array of level set functions defining the surfaces
! nlsf        : number of lsfs
! r0,s0,t0    : bottom left back corner of the part of the element. (we are now integrating in local r,s,t coordinates)
! rst         : size of the part of the element we are integrating
! icut        : returned, 0 if homogeneous element - 1 if cut element
! ileaves     : current leaf number (useless except for debugging)
! matnum      : material number, stored for output

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

use threads      
use definitions

implicit none

type(parameters),intent(in) :: params
integer,intent(in) :: level
integer,intent(in) :: levelmax
integer,intent(in) :: ndof
double precision,intent(inout) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(inout) :: bel(params%mpe*ndof)
integer,intent(in) :: icon(params%mpe)
double precision,intent(in) :: x(nnode),y(nnode),z(nnode)
integer,intent(in) :: kfix(nnode*ndof)
type(material),intent(in) :: mat(0:params%nmat)
double precision,intent(in) :: u(nnode),v(nnode),w(nnode)
double precision,intent(in) :: temp(nnode)
double precision,intent(in) :: pressure
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: e2dp(nnode)
logical,intent(inout) :: is_plastic
integer,intent(in) :: nnode
double precision,intent(inout) :: f(nnode*ndof)
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(inout) :: eviscosity
logical,intent(inout) :: vbounded
double precision,intent(inout) :: yield_ratio
double precision,intent(inout) :: frict_angle
type(thread),intent(inout) :: threadinfo
double precision,intent(inout) :: weightel
integer,intent(inout) :: matnum

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

integer :: matel,matrule,i,j,k,err
double precision,allocatable :: vol_lsf(:)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision :: prod,vol_lsf0,eviscosityp,eps,weight,yield_ratio_temp
double precision :: xmean,ymean,zmean,frict_angle_temp
logical :: is_plastic_temp,vbounded_temp

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

eps=1.d-10
icut=0

if (level.eq.0) then
   ael=0.d0
   bel=0.d0
   eviscosity=0.d0
   yield_ratio=0.d0
   frict_angle=0.d0
   if (ndof.eq.3) weightel=0.d0
endif

matel=params%materialn(0)
matrule=params%matrule

xmean=0.d0
ymean=0.d0
zmean=0.d0
do i=1,params%mpe
  xmean=xmean+x(icon(i))
  ymean=ymean+y(icon(i))
  zmean=zmean+z(icon(i))*params%vex
enddo
xmean=xmean/params%mpe
ymean=ymean/params%mpe
zmean=zmean/params%mpe

if (.not.params%excl_vol) then
   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) then
            matrule=0
          elseif (params%dommat.gt.0) then
            if (vol_lsf(params%dommat).gt.params%domvol) then                   ! If the volume of dominant material within the element exceeds domvol,
              matrule=1                                                         ! use majority rule with the dominant material
              vol_lsf(params%dommat) = 1.d0
            endif
          endif
          if (matrule.eq.1 .or. matrule.eq.2) then                              ! Use majority or minority rule
            call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,  &
                               ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,&
                               w,temp,pressure,strain,e2dp,r0,s0,t0,rst,xmean, &
                               ymean,zmean,f,ael,bel,eviscosity,weightel,      &
                               yield_ratio,frict_angle,is_plastic,vbounded,    &
                               threadinfo)
          else                                                                  ! Use divFEM
            if (level.eq.levelmax) then
              call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,&
                                 ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,&
                                 v,w,temp,pressure,strain,e2dp,r0,s0,t0,rst,   &
                                 xmean,ymean,zmean,f,ael,bel,eviscosity,       &
                                 weightel,yield_ratio,frict_angle,is_plastic,  &
                                 vbounded,threadinfo)
            else
              call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode,   &
                                 icon,kfix,ileaves,matnum,lsf,r0,s0,t0,rst,x,y,&
                                 z,u,v,w,temp,pressure,strain,e2dp,eviscosity, &
                                 weightel,ael,bel,f,yield_ratio,frict_angle,   &
                                 is_plastic,vbounded,threadinfo)
            endif
          endif
          deallocate (vol_lsf)
        endif
      enddo
      if (prod.lt.0.d0 .and. icut.eq.0) matel=params%materialn(i)
   enddo
else
   !check whether this is a cut cell
   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) then
            matrule=0
          elseif (params%dommat.gt.0) then
            if (vol_lsf(params%dommat).gt.params%domvol) then                   ! If the volume of dominant material within the element exceeds domvol,
              matrule=1                                                         ! use majority rule with the dominant material
              vol_lsf(params%dommat) = 1.d0
            endif
          endif
          if (matrule.eq.1 .or. matrule.eq.2) then                              ! Use majority or minority rule
            call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,  &
                               ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,&
                               w,temp,pressure,strain,e2dp,r0,s0,t0,rst,xmean, &
                               ymean,zmean,f,ael,bel,eviscosity,weightel,      &
                               yield_ratio,frict_angle,is_plastic,vbounded,    &
                               threadinfo)
          else                                                                  ! Use divFEM
            if (level.eq.levelmax) then
              call set_elem_props(params,mat,matrule,ndof,nlsf,nnode,icon,kfix,&
                                 ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,&
                                 v,w,temp,pressure,strain,e2dp,r0,s0,t0,rst,   &
                                 xmean,ymean,zmean,f,ael,bel,eviscosity,       &
                                 weightel,yield_ratio,frict_angle,is_plastic,  &
                                 vbounded,threadinfo)
            else
              call divide_element(params,mat,level,levelmax,ndof,nlsf,nnode,   &
                                 icon,kfix,ileaves,matnum,lsf,r0,s0,t0,rst,x,y,&
                                 z,u,v,w,temp,pressure,strain,e2dp,eviscosity, &
                                 weightel,ael,bel,f,yield_ratio,frict_angle,   &
                                 is_plastic,vbounded,threadinfo)
            endif
          endif
          deallocate (vol_lsf)
        endif
      enddo
   enddo
   if (icut.eq.0) then
     !assign material to plain cell, since at that point we know that this is a plain cell
     do i=1,nlsf
       if (lsf(1,i).lt.0.d0) matel=params%materialn(i)
     enddo
   endif
endif

! Check for overriding material transitions
!write (*,*) 'Calling mat_trans for an uncut element'
call mat_trans(params,mat,matel,xmean,ymean,zmean)
!write (*,*) 'Finished calling mat_trans for an uncut element'
matnum=matel

if (icut.eq.0) then
  call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,mat(matel)%viscosity,&
                   mat(matel)%density,mat(matel)%penalty,mat(matel)%expon,     &
                   mat(matel)%activationenergy,mat(matel)%expansion,           &
                   mat(matel)%diffusivity,mat(matel)%heat,                     &
                   mat(matel)%plasticity_type,                                 &
                   mat(matel)%plasticity_ss_type_coh,                          &
                   mat(matel)%plasticity_ss_type_phi,                          &
                   mat(matel)%plasticity_parameters,mat(matel)%fviscosity,u,v, &
                   w,temp,pressure,strain,e2dp,is_plastic_temp,nnode,f,r0,s0,  &
                   t0,rst,ileaves,eviscosityp,vbounded_temp,yield_ratio_temp,  &
                   frict_angle_temp,threadinfo,weight)

  ael=ael+aelp/(8.d0**level)
  bel=bel+belp/(8.d0**level)
  eviscosity=eviscosity+eviscosityp/(8.d0**level)   
  is_plastic=(is_plastic.or.is_plastic_temp)
  vbounded=(vbounded.or.vbounded_temp)
  yield_ratio=yield_ratio+yield_ratio_temp/(8.d0**level)
  frict_angle=frict_angle+frict_angle_temp/(8.d0**level)
  if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
endif

end subroutine make_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=nlsf-1,1,-1
    vol_lsf(i)=max(vol_lsf(i),vol_lsf(i+1))
  enddo
endif

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,ndof,nlsf,nnode,icon,kfix,       &
                          ileaves,level,matnum,vol_lsf,vol_lsf0,x,y,z,u,v,w,   &
                          temp,pressure,strain,e2dp,r0,s0,t0,rst,xmean,ymean,  &
                          zmean,f,ael,bel,eviscosity,weightel,yield_ratio,     &
                          frict_angle,is_plastic,vbounded,threadinfo)

use threads
use definitions

implicit none

type(parameters),intent(in) :: params
type(material),intent(in)   :: mat(0:params%nmat)
integer,intent(in) :: matrule
integer,intent(in) :: ndof
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: kfix(nnode*ndof)
integer,intent(in) :: ileaves
integer,intent(in) :: level
integer,intent(out) :: matnum
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) :: pressure
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: e2dp(nnode)
double precision,intent(in) :: r0,s0,t0,rst
double precision,intent(in) :: xmean,ymean,zmean
double precision,intent(inout) :: f(nnode*ndof)
double precision,intent(out) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(out) :: bel(params%mpe*ndof)
double precision,intent(out) :: eviscosity
double precision,intent(out) :: weightel
double precision,intent(out) :: yield_ratio
double precision,intent(out) :: frict_angle
logical,intent(out) :: is_plastic
logical,intent(out) :: vbounded
type(thread),intent(inout)  :: threadinfo

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

double precision :: eps,viscosity,density,penal,expon,activ,expan,diffusivity
double precision :: heat,eviscosityp,weight,yield_ratio_temp,frict_angle_temp
double precision :: plasticity_parameters(9)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision :: fvisc
integer :: matel,i
character(len=8) :: plasticity_type
character(len=16) :: plasticity_ss_type_coh,plasticity_ss_type_phi
logical :: is_plastic_temp,vbounded_temp

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

eps=1.d-10

matnum=0

select case (matrule)
case (1)                                                                        ! Assign material properties of volumetric majority material
  matel=params%materialn(sum(maxloc(vol_lsf)))
  ! Check for overriding material transitions
  !write (*,*) 'Calling mat_trans for an cut element (majority rule)'
  call mat_trans(params,mat,matel,xmean,ymean,zmean)
  !write (*,*) 'Finished calling mat_trans for an cut element (majority rule)'
  matnum=matel
  viscosity=mat(matel)%viscosity
  density=mat(matel)%density
  penal=mat(matel)%penalty
  expon=mat(matel)%expon
  activ=mat(matel)%activationenergy
  expan=mat(matel)%expansion
  diffusivity=mat(matel)%diffusivity
  heat=mat(matel)%heat
  plasticity_type=mat(matel)%plasticity_type
  plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
  plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
  plasticity_parameters=mat(matel)%plasticity_parameters
  fvisc=mat(matel)%fviscosity

case(2)                                                                         ! Assign material properties of volumetric minority material
  matel=params%materialn(sum(minloc(vol_lsf)))
  ! Check for overriding material transitions
  !write (*,*) 'Calling mat_trans for an cut element (minority rule)'
  call mat_trans(params,mat,matel,xmean,ymean,zmean)
  !write (*,*) 'Finished calling mat_trans for an cut element (minority rule)'
  matnum=matel
  viscosity=mat(matel)%viscosity
  density=mat(matel)%density
  penal=mat(matel)%penalty
  expon=mat(matel)%expon
  activ=mat(matel)%activationenergy
  expan=mat(matel)%expansion
  diffusivity=mat(matel)%diffusivity
  heat=mat(matel)%heat
  plasticity_type=mat(matel)%plasticity_type
  plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
  plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
  plasticity_parameters=mat(matel)%plasticity_parameters
  fvisc=mat(matel)%fviscosity

case default
  if (matrule.ne.0) write(*,*) 'Invalid matrule value, using divFEM'
  viscosity=0.d0
  density=0.d0
  penal=0.d0
  expon=0.d0
  activ=0.d0
  expan=0.d0
  diffusivity=0.d0
  heat=0.d0
  fvisc=0.d0
  do i=1,nlsf
    matel=params%materialn(i)
    ! Check for overriding material transitions
    !write (*,*) 'Calling mat_trans for an cut element (divFEM)'
    call mat_trans(params,mat,matel,xmean,ymean,zmean)
    !write (*,*) 'Finished calling mat_trans for an cut element (divFEM)'
    if (vol_lsf(i).ge.0.015625-eps) matnum=matnum+matel
    viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
    density=density+vol_lsf(i)*mat(matel)%density
    penal=penal+vol_lsf(i)*mat(matel)%penalty
    expon=expon+vol_lsf(i)*mat(matel)%expon
    activ=activ+vol_lsf(i)*mat(matel)%activationenergy
    expan=expan+vol_lsf(i)*mat(matel)%expansion
    diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
    heat=heat+vol_lsf(i)*mat(matel)%heat
    fvisc=fvisc+vol_lsf(i)*mat(matel)%fviscosity
  enddo
  matel=params%materialn(0)
  if (vol_lsf0.ge.0.015625-eps) matnum=matnum+100
  viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
  density=density+vol_lsf0*mat(matel)%density
  penal=penal+vol_lsf0*mat(matel)%penalty
  expon=expon+vol_lsf0*mat(matel)%expon
  activ=activ+vol_lsf0*mat(matel)%activationenergy
  expan=expan+vol_lsf0*mat(matel)%expansion
  heat=heat+vol_lsf0*mat(matel)%heat
  fvisc=fvisc+vol_lsf0*mat(matel)%fviscosity
  diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
  diffusivity=1.d0/diffusivity                                                  ! Note that some properties add geometrically not algebraically
  matel=params%materialn(sum(maxloc(vol_lsf)))
  plasticity_type=mat(matel)%plasticity_type
  plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
  plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
  plasticity_parameters=mat(matel)%plasticity_parameters

end select

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

call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity,density,     &
                 penal,expon,activ,expan,diffusivity,heat,plasticity_type,     &
                 plasticity_ss_type_coh,plasticity_ss_type_phi,                &
                 plasticity_parameters,fvisc,u,v,w,temp,pressure,strain,e2dp,  &
                 is_plastic_temp,nnode,f,r0,s0,t0,rst,ileaves,eviscosityp,     &
                 vbounded_temp,yield_ratio_temp,frict_angle_temp,threadinfo,   &
                 weight)

is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level) 
yield_ratio=yield_ratio+yield_ratio_temp/(8.d0**level)
frict_angle=frict_angle+frict_angle_temp/(8.d0**level)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)

end subroutine set_elem_props

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

subroutine divide_element(params,mat,level,levelmax,ndof,nlsf,nnode,icon,kfix, &
                         ileaves,matnum,lsf,r0,s0,t0,rst,x,y,z,u,v,w,temp,     &
                         pressure,strain,e2dp,eviscosity,weightel,ael,bel,f,   &
                         yield_ratio,frict_angle,is_plastic,vbounded,threadinfo)

use threads      
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) :: ndof
integer,intent(in) :: nlsf
integer,intent(in) :: nnode
integer,intent(in) :: icon(params%mpe)
integer,intent(in) :: kfix(nnode*ndof)
integer,intent(in) :: ileaves
integer,intent(inout) :: matnum
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) :: pressure
double precision,intent(in) :: strain(nnode)
double precision,intent(in) :: e2dp(nnode)
double precision,intent(inout) :: eviscosity
double precision,intent(inout) :: weightel
double precision,intent(inout) :: ael(params%mpe*ndof,params%mpe*ndof)
double precision,intent(inout) :: bel(params%mpe*ndof)
double precision,intent(inout) :: f(nnode*ndof)
double precision,intent(inout) :: yield_ratio
double precision,intent(inout) :: frict_angle
logical,intent(inout) :: is_plastic
logical,intent(inout) :: vbounded
type(thread),intent(inout)  :: threadinfo

!------------------------------------------------------------------------------|
!(((((((((((((((( 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 make_cut (levelp,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,u,v,w,&
                       temp,pressure,strain,e2dp,is_plastic,nnode,f,lsfp,nlsf, &
                       r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,      &
                       yield_ratio,frict_angle,params,threadinfo,weightel,     &
                       matnum)
      enddo
   enddo
enddo

deallocate (lsfp)

end subroutine divide_element

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

subroutine mat_trans (params,mat,matel,xmean,ymean,zmean)

use definitions

implicit none

type(parameters),intent(in) :: params
type(material),intent(in) :: mat(0:params%nmat)
integer,intent(inout) :: matel
double precision,intent(in) :: xmean,ymean,zmean

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

logical :: checkmattrans

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

checkmattrans=.true.
!write (*,*) 'checkmattrans at start of call: ',checkmattrans
do while (checkmattrans)
  checkmattrans=.false.
  !write (*,*) 'checkmattrans at start of loop: ',checkmattrans
  if (mat(matel)%mattrans(1).gt.0.d0 .and. xmean.lt.mat(matel)%mattrans(1)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for xmean lt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(1)
  endif
  if (mat(matel)%mattrans(2).gt.0.d0 .and. xmean.gt.mat(matel)%mattrans(2)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for xmean gt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(2)
  endif
  if (mat(matel)%mattrans(3).gt.0.d0 .and. ymean.lt.mat(matel)%mattrans(3)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for ymean lt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(3)
  endif
  if (mat(matel)%mattrans(4).gt.0.d0 .and. ymean.gt.mat(matel)%mattrans(4)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for ymean gt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(4)
  endif
  if (mat(matel)%mattrans(5).gt.0.d0 .and. zmean.lt.mat(matel)%mattrans(5)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for zmean lt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(5)
  endif
  if (mat(matel)%mattrans(6).gt.0.d0 .and. zmean.gt.mat(matel)%mattrans(6)) then
    checkmattrans=.true.
    !write (*,*) 'checkmattrans for zmean gt mattrans: ',checkmattrans
    matel=mat(matel)%transnum(6)
  endif
  !write (*,*) 'checkmattrans at end of loop: ',checkmattrans
enddo

!write (*,*) 'checkmattrans at end of call: ',checkmattrans

end subroutine mat_trans

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