!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              UPDATE_CLOUD_STRUCTURE    Feb. 2007                             |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine update_cloud_structure (cl,os,params,ninject,nremove,istep)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this file contains the operations done on a 3D cloud of particles
! to check the density of particles, accordingly add or remove particles
! and transfer the fields carried by the cloud (here the original
! position of the particles)

! First it initializes the fields to the current (thus original) position

! then we also interpolate the fields (original position) onto the
! corresponding fields of the octree solve in order for the information 
! to be available during matrix building operation

! cl is the cloud
! os is the octree solve
! npmin is the minimum number of particles in any leaf of the octree solve
! npmax is the maximum number of particles in the smallest leaves of the
!    octree solve
! ninject is number of new particles injected in the cloud
! nremove is number of particles removed from the cloud
! levelmax_oct is maximum level of octree solve

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

implicit none
      
type (cloud) cl
type (octreesolve) os
type (parameters) params
integer ninject
integer nremove
type (octreev) ov
integer istep

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

integer ierr,iproc,nproc,i,j,k,ic,jc,kk,ii
integer np,ne,nn,npnew,ie,npswap,ntest,ninjectp,kinject,kkinject,nchoice
integer,dimension(:),allocatable::loc,nloc,ntot,remove,lev,done
integer err,leaf,level,locc,ndist,nisolated,ntodo,nremovep
integer, dimension(:,:), allocatable :: imemswap

double precision x0,y0,z0,dxyz,distmax,dist,distmin,xx,yy,zz,h(8),ww,r,s,t
double precision,dimension(:),allocatable::xtest,ytest,ztest
double precision,dimension(:),allocatable::xinject,yinject,zinject
double precision,dimension(:),allocatable::straininject,x0inject,y0inject,z0inject
double precision,dimension(:),allocatable::strain,w
double precision,dimension(:,:),allocatable::memswap

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

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

np=cl%np
ne=os%nleaves
nn=os%nnode
ntest=10

! check number of particles to inject/remove to adjust the size of arrays
! in main program

allocate (ntot(ne),stat=err) ; if (err.ne.0) call stop_run ('error in allocating ntot in update_cloud$')
allocate (lev(ne),stat=err) ; if (err.ne.0) call stop_run ('error in allocating lev in update_cloud$')

ntot=0
do i=1,np
   call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
                          leaf,level,locc,x0,y0,z0,dxyz)
   ntot(leaf)=ntot(leaf)+1
enddo

call  octree_find_element_level (os%octree,os%noctree,lev,ne)

ninject=0
nremove=0
do ie=1,ne
   if (ntot(ie).lt.params%npmin) ninject=ninject+(params%npmin-ntot(ie))
   if (ntot(ie).gt.params%npmax .and. lev(ie).ge.params%levelmax_oct) &
      nremove=nremove+(ntot(ie)-params%npmax)
enddo

deallocate (ntot)
deallocate (lev)

if (ninject.gt.0.and.iproc.eq.0) write (8,*) 'Injecting ',ninject,' particles'
if (nremove.gt.0.and.iproc.eq.0) write (8,*) 'Removing ',nremove,' particles'

npswap=cl%np
if (npswap.ne.0) then
   allocate (memswap(10,npswap),stat=err) ; if (err.ne.0) call stop_run ('error in alloc memswap in main$')
   allocate (imemswap(1,npswap),stat=err) ; if (err.ne.0) call stop_run ('error in alloc imemswap in main$')
   memswap(1,:)=cl%x
   memswap(2,:)=cl%y
   memswap(3,:)=cl%z
   memswap(4,:)=cl%x0
   memswap(5,:)=cl%y0
   memswap(6,:)=cl%z0
   memswap(7,:)=cl%strain
   memswap(8,:)=cl%lsf0
   memswap(9,:)=cl%temp
   memswap(10,:)=cl%press
   imemswap(1,:)=cl%tag
endif

deallocate (cl%x)
deallocate (cl%y)
deallocate (cl%z)
deallocate (cl%x0)
deallocate (cl%y0)
deallocate (cl%z0)
deallocate (cl%strain)
deallocate (cl%temp)
deallocate (cl%press)
deallocate (cl%lsf0)
deallocate (cl%tag)

npnew=cl%np+ninject-nremove
npnew=max(npnew,npswap)

allocate (cl%x(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%x in main$')
allocate (cl%y(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%y in main$')
allocate (cl%z(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%z in main$')
allocate (cl%x0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%x0 in main$')
allocate (cl%y0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%y0 in main$')
allocate (cl%z0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%z0 in main$')
allocate (cl%strain(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%strain in main$')
allocate (cl%temp(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%temp in main$')
allocate (cl%press(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%press in main$')
allocate (cl%lsf0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%lsf0 in main$')
allocate (cl%tag(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%tag in main$')
        
if (npswap.ne.0) then
   cl%x(1:npswap)=memswap(1,:)
   cl%y(1:npswap)=memswap(2,:)
   cl%z(1:npswap)=memswap(3,:)
   cl%x0(1:npswap)=memswap(4,:)
   cl%y0(1:npswap)=memswap(5,:)
   cl%z0(1:npswap)=memswap(6,:)
   cl%strain(1:npswap)=memswap(7,:)
   cl%lsf0(1:npswap)=memswap(8,:)
   cl%temp(1:npswap)=memswap(9,:)
   cl%press(1:npswap)=memswap(10,:)
   cl%tag(1:npswap)=imemswap(1,:)
   deallocate (memswap)
   deallocate (imemswap)
else
   cl%strain=0.d0
   cl%lsf0=0.d0
   cl%temp=0.d0
   cl%press=0.d0
   cl%tag=1
endif

! first build location array
! ntot(ie) is number of particles in leaf ie
! nloc(ie) is location in array loc of where the particles belonging to ie are stored
! loc(nloc(ie):nloc(ie)+ntot(ie)-1) is list of particles belonging to ie

allocate (loc(np),stat=err) ; if (err.ne.0) call stop_run ('error in alloc loc in update_cloud$')
allocate (nloc(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc nloc in update_cloud$')
allocate (ntot(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ntot in update_cloud$')
allocate (lev(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc lev in update_cloud$')
allocate (remove(np),stat=err) ; if (err.ne.0) call stop_run ('error in alloc remove in update_cloud$')

ntot=0
do i=1,np
   call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
                          leaf,level,locc,x0,y0,z0,dxyz)
   ntot(leaf)=ntot(leaf)+1
enddo

nloc(1)=1
do ie=2,ne
   nloc(ie)=nloc(ie-1)+ntot(ie-1)
enddo

ntot=0
loc=0
do i=1,np
   call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
                          leaf,level,locc,x0,y0,z0,dxyz)
   loc(nloc(leaf)+ntot(leaf))=i
   ntot(leaf)=ntot(leaf)+1
enddo

call  octree_find_element_level (os%octree,os%noctree,lev,ne)

! first check for injection

ninjectp=ninject
ninject=0
do ie=1,ne
   if (ntot(ie).lt.params%npmin) then
      ninject=ninject+(params%npmin-ntot(ie))
   endif
enddo

allocate (xinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc xinject in update_cloud$')
allocate (yinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc yinject in update_cloud$')
allocate (zinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc zinject in update_cloud$')
allocate (straininject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc straininject in update_cloud$')
allocate (x0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc x0inject in update_cloud$')
allocate (y0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc y0inject in update_cloud$')
allocate (z0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc z0inject in update_cloud$')


if (ninject.ne.ninjectp) call stop_run ('error 1 in update_cloud$')

kinject=0
do ie=1,ne
   if (ntot(ie).lt.params%npmin) then
      kkinject=kinject+1
      do ii=1,params%npmin-ntot(ie)
         allocate (xtest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc xtest in update_cloud$')
         allocate (ytest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ytest in update_cloud$')
         allocate (ztest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ztest in update_cloud$')
         call random_number (xtest)
         call random_number (ytest)
         call random_number (ztest)
         xtest=os%x(os%icon(1,ie))+xtest*(os%x(os%icon(2,ie))-os%x(os%icon(1,ie)))
         ytest=os%y(os%icon(1,ie))+ytest*(os%y(os%icon(3,ie))-os%y(os%icon(1,ie)))
         ztest=os%z(os%icon(1,ie))+ztest*(os%z(os%icon(5,ie))-os%z(os%icon(1,ie)))
         distmax=0.d0
         nchoice=0
         do i=1,ntest
            dist=0.d0
            ndist=0
            do k=1,ntot(ie)
               dist=dist+(xtest(i)-cl%x(loc(nloc(ie)+k-1)))**2+ &
                         (ytest(i)-cl%y(loc(nloc(ie)+k-1)))**2+ &
                         (ztest(i)-cl%z(loc(nloc(ie)+k-1)))**2
               ndist=ndist+1
            enddo
            do k=1,8
               dist=dist+(xtest(i)-os%x(os%icon(k,ie)))**2+ &
                         (ytest(i)-os%y(os%icon(k,ie)))**2+ &
                         (ztest(i)-os%z(os%icon(k,ie)))**2
               ndist=ndist+1
            enddo
            do k=kkinject,kinject
               dist=dist+(xtest(i)-xinject(k))**2+ &
                         (ytest(i)-yinject(k))**2+ &
                         (ztest(i)-zinject(k))**2
               ndist=ndist+1
            enddo
            dist=dist/ndist
            if (dist.gt.distmax) then
               distmax=dist
               nchoice=i
            endif
         enddo
         if (nchoice.eq.0) call stop_run ('error in finding inject part loc$')
         ! found the random particle that is the furthest away from the others and the
         ! 8 nodes of the leaf and previously injected ones; now we inject a particle
         kinject=kinject+1
         xinject(kinject)=xtest(nchoice)
         yinject(kinject)=ytest(nchoice)
         zinject(kinject)=ztest(nchoice)
         deallocate (xtest)
         deallocate (ytest)
         deallocate (ztest)
      enddo
   endif
enddo

if (kinject.ne.ninjectp) call stop_run ('error 2 in update_cloud$')

! second check for removal

remove=0
do ie=1,ne
   if (ntot(ie).gt.params%npmax.and.lev(ie).ge.params%levelmax_oct) then
      do ii=1,ntot(ie)-params%npmax
         distmin=1.d2
         nchoice=0
         do i=1,ntot(ie)
            if (remove(loc(nloc(ie)+i-1)).eq.0) then
               dist=0.d0
               ndist=0
               xx=cl%x(loc(nloc(ie)+i-1))
               yy=cl%y(loc(nloc(ie)+i-1))
               zz=cl%z(loc(nloc(ie)+i-1))
               do k=1,ntot(ie)
                  if (remove(loc(nloc(ie)+k-1)).eq.0 .and. k.ne.i) then
                     dist=dist+(xx-cl%x(loc(nloc(ie)+k-1)))**2+ &
                               (yy-cl%y(loc(nloc(ie)+k-1)))**2+ &
                               (zz-cl%z(loc(nloc(ie)+k-1)))**2
                     ndist=ndist+1
                  endif
               enddo
               do k=1,8
                  dist=dist+(xx-os%x(os%icon(k,ie)))**2+ &
                            (yy-os%y(os%icon(k,ie)))**2+ &
                            (zz-os%z(os%icon(k,ie)))**2
                  ndist=ndist+1
               enddo
               dist=dist/ndist
               if (dist.lt.distmin) then
                  distmin=dist
                  nchoice=loc(nloc(ie)+i-1)
               endif
            endif
         enddo
         if (nchoice.eq.0) call stop_run ('error in finding inject part loc$')
         ! found the particle that is the closest to the others and the
         ! 8 nodes of the leaf; we remove it
         remove(nchoice)=1
      enddo
   endif
enddo




! now interpolate onto the nodes of the octree solve

os%strain=0.d0
allocate (strain(nn),stat=err) ; if (err.ne.0) call stop_run ('error in alloc strain in update_cloud$')
allocate (w(nn),stat=err) ; if (err.ne.0) call stop_run ('error in alloc /w in update_cloud')

strain=0.d0
w=0.d0
do i=1,np
   call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
                          leaf,level,locc,x0,y0,z0,dxyz)
   r=(cl%x(i)-x0)/dxyz*2.d0-1.d0
   s=(cl%y(i)-y0)/dxyz*2.d0-1.d0
   t=(cl%z(i)-z0)/dxyz*2.d0-1.d0
   h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
   h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
   h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
   h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
   h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
   h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
   h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
   h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
   do k=1,8
      ic=os%icon(k,leaf)
      ww=h(k)/dxyz**3
      w(ic)=w(ic)+ww
      strain(ic)=strain(ic)+ww*cl%strain(i)
   enddo
enddo

if (np.gt.0) then
   allocate (done(os%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc done in update_cloud$')
   nisolated=0
   done=0
   do i=1,os%nnode
      if (w(i).gt.0.d0) then
         os%strain(i)=strain(i)/w(i)
      else
         done(i)=1
         nisolated=nisolated+1
      endif
   enddo
   if (iproc.eq.0) write (8,*) 'There were ',nisolated,' octree nodes'

111   continue

   strain=0.d0
   w=0.d0
   do ie=1,os%nleaves
      do k=1,8
         ic=os%icon(k,ie)
         if (done(ic).eq.1) then
            do kk=1,8
               jc=os%icon(kk,ie)
               if (done(jc).eq.0) then
                  ww=1.d0/dxyz**3
                  w(ic)=w(ic)+ww
                  strain(ic)=strain(ic)+os%strain(jc)*ww
               endif
            enddo
         endif
      enddo
   enddo

   ntodo=0
   do i=1,os%nnode
      if (done(i).eq.1) then
         if (w(i).gt.0.d0) then
            done(i)=0
            os%strain(i)=strain(i)/w(i)
         else
            ntodo=ntodo+1
         endif
      endif
   enddo

   if (iproc.eq.0) write (8,*) ntodo,' Isolated left '

   if (ntodo.ne.0) goto 111

   deallocate (done)

endif

deallocate (strain)
deallocate (w)

! now calculate the field value at the new (injected points)

straininject=0.d0
x0inject=0.d0
y0inject=0.d0
z0inject=0.d0
do i=1,ninject
   call octree_find_leaf (os%octree,os%noctree,xinject(i),yinject(i),zinject(i), &
                          leaf,level,locc,x0,y0,z0,dxyz)
   r=(xinject(i)-x0)/dxyz*2.d0-1.d0
   s=(yinject(i)-y0)/dxyz*2.d0-1.d0
   t=(zinject(i)-z0)/dxyz*2.d0-1.d0
   h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
   h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
   h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
   h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
   h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
   h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
   h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
   h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
   do k=1,8
      straininject(i)=straininject(i)+os%strain(os%icon(k,leaf))*h(k)
      x0inject(i)=x0inject(i)+os%x(os%icon(k,leaf))*h(k)
      y0inject(i)=y0inject(i)+os%y(os%icon(k,leaf))*h(k)
      z0inject(i)=z0inject(i)+os%z(os%icon(k,leaf))*h(k)
   enddo
enddo

! modify the particle cloud for removal/injection of particles

nremovep=0
do i=1,np
   if (remove(i).eq.1) then
      nremovep=nremovep+1
   else
      cl%x(i-nremovep)=cl%x(i)
      cl%y(i-nremovep)=cl%y(i)
      cl%z(i-nremovep)=cl%z(i)
      cl%x0(i-nremovep)=cl%x0(i)
      cl%y0(i-nremovep)=cl%y0(i)
      cl%z0(i-nremovep)=cl%z0(i)
      cl%strain(i-nremovep)=cl%strain(i)
      cl%lsf0(i-nremovep)=cl%lsf0(i)
      cl%temp(i-nremovep)=cl%temp(i)
      cl%press(i-nremovep)=cl%press(i)
      cl%tag(i-nremovep)=cl%tag(i)

   endif
enddo

if (nremovep.ne.nremove) call stop_run ('error in update_cloud/remove$')

np=np-nremovep

do i=1,ninject
   np=np+1
   cl%x(np)=xinject(i)
   cl%y(np)=yinject(i)
   cl%z(np)=zinject(i)
   cl%x0(np)=x0inject(i)
   cl%y0(np)=y0inject(i)
   cl%z0(np)=z0inject(i)
   cl%strain(np)=straininject(i)
   cl%lsf0(np)=0.d0
   cl%temp(np)=0.d0
   cl%press(np)=0.d0
   cl%ntag=cl%ntag+1
   cl%tag(np)=cl%ntag
enddo

! in first step, we set the original particle positions to their
! current position

if (cl%np.eq.0) then
   do i=1,np
      cl%x0(i)=cl%x(i)
      cl%y0(i)=cl%y(i)
      cl%z0(i)=cl%z(i)
      cl%strain(i)=0.d0
      cl%temp(i)=0.d0
      cl%press(i)=0.d0
      cl%lsf0(i)=0.d0
      cl%tag(i)=i
   enddo
   cl%ntag=np
endif

cl%np=np

deallocate (xinject)
deallocate (yinject)
deallocate (zinject)
deallocate (straininject)
deallocate (x0inject,y0inject,z0inject)
deallocate (loc)
deallocate (nloc)
deallocate (ntot)
deallocate (remove)
deallocate (lev)

if (iproc.eq.0) then
   write (8,*) 'Strain on osolve ',minval(os%strain),maxval(os%strain)
   write (8,*) cl%np,' particles in 3D cloud'
end if

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