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

subroutine embed_surface_in_octree (osolve,params,surface,is,istep,iter,       &
                                   threadinfo)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine embed a surface in a new octree, computes a lsf, and then  
! refines osolve accordingly and interpolates the lsf on osolve

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

use definitions
!use mpi
use threads

implicit none

include 'mpif.h'

type(octreesolve) osolve
type(parameters) params
type (octreelsf) olsf
type (sheet) surface
integer is
integer istep
integer iter
type (thread) threadinfo

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

integer,dimension(:),allocatable::levs
type (octreesolve) oint
integer ierr,iproc,nproc
integer ie,ilsf,k
double precision xxx,yyy,zzz
integer ioctree_number_of_elements
external ioctree_number_of_elements
character*72  :: shift


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

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

shift=' '

if(iproc.eq.0) write(8,*) 'embedding surface',is,'/',params%ns

! we want to create the octree olsf in which the surface under consideration will be imbedded.
! we first create a uniform octree at level leveluniform_oct

olsf%noctree=params%noctreemax
allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%octree', 'embed_surf...',size(olsf%octree),'int',+1)

call octree_init (olsf%octree,olsf%noctree)
call octree_create_uniform (olsf%octree,olsf%noctree,params%leveluniform_oct)

! we then improve the octree to best fit/represent the surface

call build_surface_octree (surface,olsf,params,threadinfo)

if (iproc.eq.0) write (8,*) olsf%nleaves,' leaves in lsf octree'

! we then calculate the  connectivity and node coordinates of the surface octree

olsf%nnode=olsf%nleaves*3
allocate (olsf%icon(8,olsf%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',+1)
   
call octree_find_node_connectivity (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%x,olsf%y,olsf%z,olsf%nnode)
! olsf%nnode has been changed by octree_find_node_connectivity, so re-size x,y,z
call octreelsf_shrink_xyz(olsf, threadinfo, params)

if (params%debug.gt.1) then
  write(threadinfo%Logunit,*) '[same] embed_surface_in_octree: olsf%nnode ',olsf%nnode
  call flush(threadinfo%Logunit) 
endif

allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf', 'embed_surf...',size(olsf%lsf),'dp',+1)

! we then calculate the lsf on the olsf octree

call calculate_lsf (olsf%lsf,olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,  &
                   olsf%nnode,surface%x,surface%y,surface%z,surface%xn,        &
                   surface%yn,surface%zn,surface%icon,surface%nt,              &
                   surface%nsurface,params%levelmax_oct,params%normaladvect)

if (params%debug>=2) call output_octree_lsf (olsf,surface,is,istep,iter)


!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! we create and allocate the oint octree that will receive the current osolve octree.
! this will allow us to resize osolve so that the new surface can be embedded in it.

oint%noctree=osolve%noctree
oint%nleaves=osolve%nleaves
oint%nnode=osolve%nnode
oint%nlsf=osolve%nlsf
! write(threadinfo%mem_heap_unit,*) 'oint noctree,nleaves,nnode,nlsf are', oint%noctree,oint%nleaves,oint%nnode,oint%nlsf

allocate (oint%octree(params%noctreemax),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%octree', 'embed_surf...',size(oint%octree),'int',+1)
allocate (oint%icon(8,oint%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%icon', 'embed_surf...',size(oint%icon),'int',+1)
allocate (oint%lsf(oint%nnode,oint%nlsf),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%lsf', 'embed_surf...',size(oint%lsf),'dp',+1)
allocate (oint%x(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%x', 'embed_surf...',size(oint%x),'dp',+1)
allocate (oint%y(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%y', 'embed_surf...',size(oint%y),'dp',+1)
allocate (oint%z(oint%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'oint%z', 'embed_surf...',size(oint%z),'dp',+1)

oint%octree=osolve%octree
oint%icon=osolve%icon
if (size(oint%lsf).ne.size(osolve%lsf)) then 
   write(*,*) 'oint%lsf = osolve%lsf; sizes differ', size(oint%lsf), size(osolve%lsf)
endif
oint%lsf=osolve%lsf
oint%x=osolve%x
oint%y=osolve%y
oint%z=osolve%z

! having replicated osolve, we can deallocate it
if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',-1)
deallocate (osolve%icon)
if (params%debug.gt.1) call heap (threadinfo,'osolve%x','embed_surf...',size(osolve%x),'dp',-1)
deallocate (osolve%x)    
if (params%debug.gt.1) call heap (threadinfo,'osolve%y','embed_surf...',size(osolve%y),'dp',-1)
deallocate (osolve%y) 
if (params%debug.gt.1) call heap (threadinfo,'osolve%z','embed_surf...',size(osolve%z),'dp',-1)
deallocate (osolve%z)
if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf','embed_surf...',size(osolve%lsf),'dp',-1)
deallocate (osolve%lsf) 

! we reinitialise it
call octree_init (osolve%octree,osolve%noctree)

osolve%octree=oint%octree

! we modify the osolve%octree by embedding the olsf in it
allocate (levs(olsf%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',+1)
call octree_find_element_level (olsf%octree,olsf%noctree,levs,olsf%nleaves)

do ie=1,olsf%nleaves
   xxx=sum(olsf%x(olsf%icon(:,ie)))/8.d0
   yyy=sum(olsf%y(olsf%icon(:,ie)))/8.d0
   zzz=sum(olsf%z(olsf%icon(:,ie)))/8.d0
   call octree_create_from_particles (osolve%octree,osolve%noctree,xxx,yyy,zzz,1,levs(ie))
enddo

if (params%debug.gt.1) call heap (threadinfo,'levs', 'embed_surf...',size(levs),'dp',-1)
deallocate (levs)


! we smoothen osolve
call octree_smoothen (osolve%octree,osolve%noctree)


! the number of leaves in osolve is computed by means of an octree function, and the 
! number of nodes is over estimated
osolve%nleaves = ioctree_number_of_elements(osolve%octree,osolve%noctree)
osolve%nnode   = osolve%nleaves*3
osolve%nlsf    = params%ns

allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','embed_surf...',size(osolve%icon),'int',+1)
allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%x',   'embed_surf...',size(osolve%x),'dp',+1)
allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%y',   'embed_surf...',size(osolve%y),'dp',+1)
allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%z',   'embed_surf...',size(osolve%z),'dp',+1)


! we build the icon array of osolve. the routine also outputs the real number of nodes osolve%nnode
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
                                    osolve%icon,osolve%nleaves,   &
                                    osolve%x,osolve%y,osolve%z,osolve%nnode)
! osolve%nnode has been changed by octree_find_node_connectivity, so  re-size x,y,z
call octreesolve_shrink_xyz(osolve, threadinfo, params)

! now that osolve%nnode is known we can allocate osolve%lsf
allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf', 'embed_surf...',size(osolve%lsf),'dp',+1)

! if we are dealing with the last surface, nodes are renumbered
if (params%renumber_nodes .and. is.eq.params%ns) then
   call octree_renumber_nodes (osolve%icon,osolve%nleaves, &
                               osolve%x,osolve%y,osolve%z,osolve%nnode)
end if

! the lsf values of the previous is-1 surfaces are known on the oint octree which is no more
! the same as osolve, so we have to interpolates these lsfs on final octree osolve


do ilsf=1,is-1
   do k=1,osolve%nnode

      call octree_interpolate (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf),oint%nnode, &
                               osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))

!      call octree_interpolate3 (oint%octree,oint%noctree,oint%icon,oint%nleaves,oint%lsf(1,ilsf),          &
!                                oint%x(1:oint%nnode),oint%y(1:oint%nnode),oint%z(1:oint%nnode),oint%nnode, &
!                               osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,ilsf))

   enddo
enddo


! finally the lsf corresponding to the current surface also has to be interpolated onto osolve
! since it was computed on a dedicated olsf 

do k=1,osolve%nnode
   call octree_interpolate (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,olsf%nnode, &
                           osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
!   call octree_interpolate3 (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,olsf%x(1:olsf%nnode),olsf%y(1:olsf%nnode),olsf%z(1:olsf%nnode),olsf%nnode, &
!                            osolve%x(k),osolve%y(k),osolve%z(k),osolve%lsf(k,is))
enddo

if (params%debug.gt.1) call heap (threadinfo,'olsf%x','embed_surf...',size(olsf%x),'dp',-1)
deallocate (olsf%x)
if (params%debug.gt.1) call heap (threadinfo,'olsf%y','embed_surf...',size(olsf%y),'dp',-1)
deallocate (olsf%y)
if (params%debug.gt.1) call heap (threadinfo,'olsf%z','embed_surf...',size(olsf%z),'dp',-1)
deallocate (olsf%z)
if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','embed_surf...',size(olsf%lsf),'dp',-1)
deallocate (olsf%lsf)
if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','embed_surf...',size(olsf%icon),'int',-1)
deallocate (olsf%icon)
if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','embed_surf...',size(olsf%octree),'int',-1)
deallocate (olsf%octree)

if (params%debug.gt.1) call heap (threadinfo,'oint%x','embed_surf...',size(oint%x),'dp',-1)
deallocate (oint%x)
if (params%debug.gt.1) call heap (threadinfo,'oint%y','embed_surf...',size(oint%y),'dp',-1)
deallocate (oint%y)
if (params%debug.gt.1) call heap (threadinfo,'oint%z','embed_surf...',size(oint%z),'dp',-1)
deallocate (oint%z)
if (params%debug.gt.1) call heap (threadinfo,'oint%lsf','embed_surf...',size(oint%lsf),'dp',-1)
deallocate (oint%lsf)
if (params%debug.gt.1) call heap (threadinfo,'oint%icon','embed_surf...',size(oint%icon),'int',-1)
deallocate (oint%icon)
if (params%debug.gt.1) call heap (threadinfo,'oint%octree','embed_surf...',size(oint%octree),'int',-1)
deallocate (oint%octree)

if (params%debug>=1 .and. iproc==0) then
   write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%nleaves
end if

if (params%debug.gt.1) then
  write(threadinfo%Logunit,*) '[same] build_surface_octree: osolve%nleaves ', osolve%nleaves
  call flush(threadinfo%Logunit)
endif

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