!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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,is_plastic,     &
                           nnode,f,lsf,nlsf,r0,s0,t0,rst,icut,ileaves,         &
                           eviscosity,vbounded,params,threadinfo, &
                           weightel)

!------------------------------------------------------------------------------|
!(((((((((((((((( 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 materail 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)


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

use threads      
use definitions

implicit none

type(parameters) params
integer level
integer levelmax
integer ndof
double precision ael(params%mpe*ndof,params%mpe*ndof)
double precision bel(params%mpe*ndof)
integer icon(params%mpe)
double precision x(nnode),y(nnode),z(nnode)
integer kfix(nnode*ndof)
type (material) mat(0:params%nmat)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
logical is_plastic
integer nnode
double precision f(nnode*ndof)
double precision lsf(params%mpe,nlsf)
integer nlsf
double precision r0,s0,t0,rst
integer icut
integer ileaves
double precision eviscosity
logical vbounded
type (thread) threadinfo
double precision weightel

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

double precision volmax
double precision r(params%mpe),s(params%mpe),t(params%mpe)
double precision aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision h(params%mpe),vol_lsf0,prod
double precision r0p,s0p,t0p,rstp
double precision viscosity,density,penal,expon,diffusivity,heat,activ,expan
character (len=8) plasticity_type
double precision plasticity_parameters(9)
double precision,dimension(:,:),allocatable::lsfp
double precision,dimension(:),allocatable::vol_lsf
integer i,j,k,ii,jj,kk,jcut,levelp, err
integer matel
double precision eviscosityp
logical  is_plastic_temp,vbounded_temp
double precision weight
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

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

matel=params%materialn(0)

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) goto 222
      enddo
      if (prod.lt.0.d0) then
         matel=params%materialn(i)
      endif
   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) goto 222
      enddo
   end do
   !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)
   end do

   !if (lsf(1,1).lt.0) matel=params%materialn(1)
   !if (lsf(1,2).lt.0) matel=params%materialn(2)
end if

!=====[end new stuff]=====

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_parameters, &
                  u,v,w,temp,pressure,strain,is_plastic_temp,nnode, &
                  f,r0,s0,t0,rst,ileaves,eviscosityp,vbounded_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)
!if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
weightel=weightel+weight/(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

    ! 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
    viscosity=0.d0
    density=0.d0
    penal=0.d0
    expon=0.d0
    activ=0.d0
    expan=0.d0
    diffusivity=0.d0
    heat=0.d0
    do i=1,nlsf
       matel=params%materialn(i)
       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
    enddo
    vol_lsf0=1.d0-sum(vol_lsf)
    matel=params%materialn(0)
    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
    diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
    diffusivity=1.d0/diffusivity! note that some properties add geometrically not algebraically

    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

!    if (plasticity_type/='No') then
!       print *,vol_lsf
!       call stop_run ('pb$')
!    end if

    call make_matrix (params,ndof,aelp,belp,icon,x,y,z, &
                        kfix,viscosity,density, &
                        penal,expon, &
                        activ,expan, &
                        diffusivity,heat,plasticity_type, &
                        plasticity_parameters, &
                        u,v,w,temp,pressure,strain,is_plastic_temp,nnode, &
                        f,r0,s0,t0,rst,ileaves,eviscosityp,vbounded_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) 
    !if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
    weightel=weightel+weight/(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 make_cut (levelp,levelmax, &
                           ndof,ael,bel,icon,x,y,z, &
                           kfix,mat, &
                           u,v,w,temp,pressure,strain,is_plastic,nnode, &
                           f,lsfp,nlsf,r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,params,threadinfo, &
                           weightel)
      enddo
   enddo
enddo

deallocate (lsfp)

return

end

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