!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              REFINE_SURFACE    May 2008                                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
      
subroutine refine_surface(params,surface,surface0,threadinfo,nadd,nrem,is,&
           istep,ref_count) 

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine builds an edge array between a set of particles on a surface.
! It uses the delaunay triangulation and then steps through the triangles
! to build the edge information.
! ed is the computed edge array
! nedge is the number of edges
! nedgepernode,nodenodenumber and nodenodenumber contain the list of edges 
! that start from each node; for a given node i
! their number is nedgepernode(i), the edge number in the list of edges is
! nodeedgenumber(j,i) and the node at the end of the edge is
! nodenodenumber(j,i) for j=1,nedgepernode(i)

! this routine builds the boolean 'refine' array of length nedge (number of edges 
! in the triangulation connecting the particles) that contains the list of 
! edges to be refined.  

! surface is the sheet/surface to be refined
! it will contain the new number of particles
! ed is the computed edge array
! refine is the integer array determining the edges to be refined
! nedge is the number of edges
! nadd is the number of edges to be changed/split
! nedgepernode, nodenodenumber,nodeedgenumber are arrays containing edge info
! nnmax is the maximum number of nn

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

implicit none

include 'mpif.h'

type (parameters) params
type (sheet) surface 
type (sheet) surface0 
type (thread) threadinfo
integer nadd,nrem
integer is,istep
integer ref_count

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

integer naddp,nremp
integer nnmax
type (edge), dimension(:), allocatable :: ed
type (edge), dimension(:), allocatable :: edswap
integer,dimension(:)  ,allocatable :: nedgepernode
integer,dimension(:)  ,allocatable :: refine_list,remove_list
integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
logical,dimension(:),allocatable::refine,refinep,remove,removep
character*72  :: shift
integer ie,j,jedge,kp,kpp,k,inode,inodep,inodepp,iadd,irem
integer ntriangle
integer nedge
integer nedgen,nsurfacen
double precision dist1,distmax,prod,distmin,eps
double precision, external :: dist


integer err,ierr,iproc,nproc,i,iedge,i1,i2
integer nnodeint,nelemint,nedgeint,nelemmax,nnodemax,nedgemax
double precision,dimension(:,:),allocatable::memswap
integer,dimension(:,:),allocatable::iconswap
double precision xm,ym,zm,xyzn

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

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
         
shift=' '
eps=1.d-8

!------------------------
! build edge information 
!------------------------

nsurfacen=surface%nsurface
ntriangle=surface%nt
nedge=nsurfacen+ntriangle-1
if (surface%closed.eq.1) nedge=nedge-1
nnmax=12 ! nnmax is maximum number of neighbours in triangulation (should be around 6 on average)

allocate (ed(nedge),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'ed','refine_surface',size(ed),'int',+1)
allocate (nedgepernode(nsurfacen),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',+1)
allocate (nodenodenumber(nnmax,nsurfacen),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'nodenodenumber','refine_surface',size(nodenodenumber),'int',+1)
allocate (nodeedgenumber(nnmax,nsurfacen),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'nodeedgenumber','refine_surface',size(nodeedgenumber),'int',+1)

ed%t1=0
ed%t2=0
nedgepernode=0
iedge=0

do ie=1,ntriangle
  do k=1,3
    inode=surface%icon(k,ie)
    kp=mod(k,3)+1
    kpp=mod(k+1,3)+1
    inodep=surface%icon(kp,ie)
    inodepp=surface%icon(kpp,ie)
    do j=1,nedgepernode(inodep)
      if (nodenodenumber(j,inodep).eq.inode) then
        jedge=nodeedgenumber(j,inodep)
        ed(jedge)%m2=inodepp
        ed(jedge)%t2=ie
        goto 111
      endif
    enddo
    iedge=iedge+1
    if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
    nedgepernode(inode)=nedgepernode(inode)+1
    if (nedgepernode(inode).gt.nnmax) then
      if (iproc.eq.0) write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
      call stop_run ('too many neighbours$')
    endif
    nodenodenumber(nedgepernode(inode),inode)=inodep
    nodeedgenumber(nedgepernode(inode),inode)=iedge
    ed(iedge)%n1=inode
    ed(iedge)%n2=inodep
    ed(iedge)%m1=inodepp
    ed(iedge)%t1=ie
    111 continue
  enddo
enddo

if (iedge/=nedge) call stop_run ('pb_b in refine_surface$')

!------------------------------------
! find edges that need to be refined
!------------------------------------

!distmax=1.d0/2.d0**surface%levelt*sqrt(2.d0)*surface%stretch
distmax=1.d0/(2.d0**surface%levelt+1)*sqrt(2.d0)*surface%stretch   !opla
distmin=1.d0/(2.d0**(surface%levelt+4)+1)*sqrt(2.d0)

allocate (refine(nedge),stat=err)
allocate (refinep(nedge),stat=err)
refine=.false.
refinep=.false.

if (params%remove_surf_pts) then
  allocate (remove(nedge),stat=err)
  allocate (removep(nedge),stat=err)
  remove=.false.
  removep=.false.
endif

do iedge=1+iproc,nedge,nproc
   i1=ed(iedge)%n1
   i2=ed(iedge)%n2
   dist1=dist(surface%x(i1),surface%x(i2),surface%y(i1),surface%y(i2),         &
             surface%z(i1),surface%z(i2))
   !dist1=dist(surface%x(i1),surface%x(i2),surface%y(i1),surface%y(i2),         &
   !          surface%z(i1),surface%z(i2),params%distance_exponent)
   prod=surface%xn(i1)*surface%xn(i2)+surface%yn(i1)*surface%yn(i2)+surface%zn(i1)*surface%zn(i2)
   !if (dist1.gt.distmax .or. prod.lt.prodmin) refinep(iedge)=1
   refinep(iedge)=(dist1.gt.distmax) 
   if (params%remove_surf_pts) removep(iedge)=(dist1.gt.eps .and. dist1.le.distmin)
   !if (dist1.le.distmin) print *,'dist1 (',dist1,') less than distmin'
enddo

call mpi_allreduce (refinep,refine,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
if (params%remove_surf_pts) call mpi_allreduce (removep,remove,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)

!-----------------------
! compute nadd and naddp
!-----------------------

nadd=count(refine)
allocate(refine_list(nadd))
iadd=0
naddp=0

if (params%remove_surf_pts) then
  nrem=count(remove)
  allocate(remove_list(nrem))
  irem=0
  nremp=0
else
  nrem=0
endif

do iedge=1,nedge
  if (refine(iedge)) then
     iadd=iadd+1
     if (ed(iedge)%t1.eq.0 .or. ed(iedge)%t2.eq.0) naddp=naddp+1
     refine_list(iadd)=iedge
  endif
  if (params%remove_surf_pts) then
    if (remove(iedge)) then
      irem=irem+1
      if (ed(iedge)%t1.eq.0 .or. ed(iedge)%t2.eq.0) nremp=nremp+1
      !remove_list(irem)=iedge
      remove_list(irem)=ed(iedge)%n1
    endif
  endif
enddo

! Modify remove_list to account for removal of nodes from surface
if (params%remove_surf_pts) then
  if (nrem.gt.0) then
    do i=1,nrem
      do j=i,nrem
        if (remove_list(i).lt.remove_list(j)) remove_list(j)=remove_list(j)-1
      enddo
    enddo
  endif
  if (irem/=nrem) call stop_run ('pb_b in refine_surface$')
  deallocate (remove)
  deallocate (removep)
endif

if (iadd/=nadd) call stop_run ('pb_b in refine_surface$')

deallocate (refine)
deallocate (refinep)

if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i2,a,i4,a)') shift//'S.',is,':', nadd,' added ptcls in refine_surface'

!-----------------------------------
! resizing ed (allocate/deallocate)
!-----------------------------------

nedgen=(surface%nsurface+nadd)+(surface%nt+(nadd-naddp)*2+naddp)-1
if (surface%closed.eq.1) nedgen=nedgen-1

! Not sure this is correct...
!nedgen=(surface%nsurface+(nadd-nrem))+(surface%nt+(((nadd-naddp)*2+naddp)-((nrem-nremp)*2+nremp))-1

allocate (edswap(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main$')
edswap(1:nedge)=ed(1:nedge)
deallocate (ed)
allocate (ed(nedgen),stat=err) ; if (err.ne.0) call stop_run ('Error alloc edswap in main$')
ed(1:nedge)=edswap(1:nedge)
deallocate (edswap)

!----------------------------------
! prepare memory for the refinment
!----------------------------------

nnodemax=surface%nsurface+nadd
nelemmax=surface%nt+(nadd-naddp)*2+naddp
nedgemax=nnodemax+nelemmax-1
if (surface%closed.eq.1) nedgemax=nedgemax-1

if (nedgemax.ne.nedgen) call stop_run ('error in counting number of edges in refine_surface$')

allocate (memswap(surface%nsurface,8),stat=err) ; if (err.ne.0) call stop_run('Error alloc memswap in refine_surface$')
allocate (iconswap(3,surface%nt),stat=err) ; if (err.ne.0) call stop_run('Error alloc iconswap in refine_surface$')
      
memswap(:,1)=surface%x
memswap(:,2)=surface%y
memswap(:,3)=surface%z
memswap(:,4)=surface%xn
memswap(:,5)=surface%yn
memswap(:,6)=surface%zn
memswap(:,7)=surface%r
memswap(:,8)=surface%s
iconswap(1:3,1:surface%nt)=surface%icon(1:3,1:surface%nt)

deallocate (surface%x)
deallocate (surface%y)
deallocate (surface%z)
deallocate (surface%xn)
deallocate (surface%yn)
deallocate (surface%zn)
deallocate (surface%r)
deallocate (surface%s)
deallocate (surface%icon)

allocate (surface%x(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in refine_surface$')
allocate (surface%y(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in refine_surface$')
allocate (surface%z(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in refine_surface$')
allocate (surface%xn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xn in refine_surface$')
allocate (surface%yn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc yn in refine_surface$')
allocate (surface%zn(nnodemax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zn in refine_surface$')
allocate (surface%r(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc r in refine_surface$')
allocate (surface%s(nnodemax), stat=err) ; if (err.ne.0) call stop_run ('Error alloc s in refine_surface$')
allocate (surface%icon(3,nelemmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc icon in refine_surface$')
      
surface%x(1:surface%nsurface)=memswap(:,1)
surface%y(1:surface%nsurface)=memswap(:,2)
surface%z(1:surface%nsurface)=memswap(:,3)
surface%xn(1:surface%nsurface)=memswap(:,4)
surface%yn(1:surface%nsurface)=memswap(:,5)
surface%zn(1:surface%nsurface)=memswap(:,6)
surface%r(1:surface%nsurface)=memswap(:,7)
surface%s(1:surface%nsurface)=memswap(:,8)
surface%icon(1:3,1:surface%nt)=iconswap(1:3,1:surface%nt)

deallocate (memswap)
deallocate (iconswap)

!---------------------
! updates the surface
!---------------------

nnodeint=surface%nsurface
nelemint=surface%nt
nedgeint=nedge

do i=1,nadd
   iedge=refine_list(i)
   i1=ed(iedge)%n1
   i2=ed(iedge)%n2
   call middle (xm,ym,zm,i1,i2,surface%x,surface%y,surface%z, &
                surface%xn,surface%yn,surface%zn,nnodeint)
   nnodeint=nnodeint+1
   surface%x(nnodeint)=xm
   surface%y(nnodeint)=ym
   surface%z(nnodeint)=zm
   surface%r(nnodeint)=(surface%r(i1)+surface%r(i2))/2.d0
   surface%s(nnodeint)=(surface%s(i1)+surface%s(i2))/2.d0
   xm=(surface%xn(i1)+surface%xn(i2))/2.d0
   ym=(surface%yn(i1)+surface%yn(i2))/2.d0
   zm=(surface%zn(i1)+surface%zn(i2))/2.d0
   xyzn=sqrt(xm**2+ym**2+zm**2)
   surface%xn(nnodeint)=xm/xyzn
   surface%yn(nnodeint)=ym/xyzn
   surface%zn(nnodeint)=zm/xyzn
   call update_icon (surface%icon,nelemmax,ed,nedgemax,iedge,nnodeint,nelemint,nedgeint)
enddo

if (iproc.eq.0 .and. params%debug>=1 .and. params%remove_surf_pts) write(*,'(a,i2,a,i4,a)') shift//'S.',is,':', nrem,' removed ptcls in refine_surface'

do i=1,size(ed%n1)
   i1=ed(i)%n1
   i2=ed(i)%n2
   if (i1 == i2 .and. iproc == 0) then
     print *,'*******************************************************************'
     print *,'Nodes for edge ',i,' have the same spatial coordinates'
     print *,'The edges nodes are are node i1: ',i1,' and i2: ',i2
     print *,'x1:',surface%x(i1),' y1: ',surface%y(i1),' z1: ',surface%z(i1)
     print *,'x2:',surface%x(i2),' y2: ',surface%y(i2),' z2: ',surface%z(i2)
     print *,'*******************************************************************'
   endif
enddo

surface%nsurface=nnodemax
surface%nt=nelemmax

if (params%remove_surf_pts) then
  do i=1,nrem
    i1=remove_list(i)
    call remove_point(i1,surface,0,surface0)
  end do
endif

if (.not.params%normaladvect) then
   call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
                         surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if

if (nnodeint.ne.nnodemax) call stop_run ('Error counting nodes in refine_surface$')
if (nelemint.ne.nelemmax) call stop_run ('Error counting elements in refine_surface$')

if (nrem.gt.0) then
  nnodemax=surface%nsurface
  nelemmax=surface%nt
  nedgemax=nnodemax+nelemmax-1
endif

surface%nsurface=nnodemax
surface%nt=nelemmax
nedge=nedgemax

deallocate(refine_list)
if (params%remove_surf_pts) deallocate(remove_list)

if (params%debug.gt.1) call heap (threadinfo,'ed','refine_surface',size(ed),'int',-1)
deallocate (ed)         
if (params%debug.gt.1) call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',-1)
deallocate (nedgepernode)
if (params%debug.gt.1) call heap (threadinfo,'nodenodenumber','refine_surface',size(nodenodenumber),'int',-1)
deallocate (nodenodenumber)
if (params%debug.gt.1) call heap (threadinfo,'nodeedgenumber','refine_surface',size(nodeedgenumber),'int',-1)
deallocate (nodeedgenumber)

if (params%debug .ge.2) call output_surf(surface,is,'after_refine',istep,ref_count) 

end subroutine  

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

subroutine middle (xm,ym,zm,i1,i2,xx,yy,zz,xxn,yyn,zzn,nnode)
  
implicit none

double precision xx(nnode),yy(nnode),zz(nnode)
double precision xxn(nnode),yyn(nnode),zzn(nnode)
double precision xm,ym,zm,xnm,ynm,znm,alpha,dist12,prod,xyzn
integer i1,i2,nnode
  
xm=(xx(i2)+xx(i1))/2.d0
ym=(yy(i2)+yy(i1))/2.d0
zm=(zz(i2)+zz(i1))/2.d0
  
return
  
! what follows is an attempt at defining the middle of the segment
! outside of the euclidian line
!xnm=(xxn(i2)+xxn(i1))/2.d0
!ynm=(yyn(i2)+yyn(i1))/2.d0
!znm=(zzn(i2)+zzn(i1))/2.d0
!xyzn=sqrt(xnm**2+ynm**2+znm**2)
!xnm=xnm/xyzn
!ynm=ynm/xyzn
!znm=znm/xyzn
!dist12=sqrt((xx(i2)-xx(i1))**2+(yy(i2)-yy(i1))**2+(zz(i2)-zz(i1))**2)
!prod=xxn(i1)*xxn(i2)+yyn(i1)*yyn(i2)+zzn(i1)*zzn(i2)
!alpha=dist12*(1./prod**2-1.)/4.d0
!xm=xm+xnm*alpha
!ym=ym+ynm*alpha
!zm=zm+ynm*alpha
  
end subroutine

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
  
subroutine update_icon (icon,nelemmax,ed,nedgemax,iedge,nnodeint,nelemint,nedgeint)
  
use definitions
  
implicit none
  
integer icon(3,nelemmax),nelemmax,nedge,iedge,nnodeint,nelemint
integer jedge,nedgemax,nedgeint
type(edge) ed(nedgemax)
integer k,k12,k21
integer n1,n2,m1,m2,t1,t2,t1n,t2n
  
nedge=nedgeint

n1=ed(iedge)%n1
n2=ed(iedge)%n2
m1=ed(iedge)%m1
m2=ed(iedge)%m2
t1=ed(iedge)%t1
t2=ed(iedge)%t2

if (t1.ne.0) then
   do k=1,3
      if (icon(k,t1).eq.n1) k12=k
   enddo
   icon(k12,t1)=nnodeint
   nelemint=nelemint+1
   t1n=nelemint
   icon(1,nelemint)=m1
   icon(2,nelemint)=n1
   icon(3,nelemint)=nnodeint
endif

if (t2.ne.0) then
   do k=1,3
      if (icon(k,t2).eq.n2) k21=k
   enddo
   icon(k21,t2)=nnodeint
   nelemint=nelemint+1
   t2n=nelemint
   icon(1,nelemint)=m2
   icon(2,nelemint)=n2
   icon(3,nelemint)=nnodeint
endif

ed(iedge)%n2=nnodeint
if (t1.ne.0) then
   ed(iedge)%t1=t1n
else
   ed(iedge)%t1=0
endif

if (t1.ne.0) then
   nedgeint=nedgeint+1
   ed(nedgeint)%n1=m1
   ed(nedgeint)%n2=nnodeint
   ed(nedgeint)%m1=n2
   ed(nedgeint)%m2=n1
   ed(nedgeint)%t1=t1
   ed(nedgeint)%t2=t1n
endif

nedgeint=nedgeint+1
ed(nedgeint)%n1=n2
ed(nedgeint)%n2=nnodeint

if (t1.ne.0) then
   ed(nedgeint)%m2=m1
   ed(nedgeint)%t2=t1
else
   ed(nedgeint)%t2=0
endif

if (t2.ne.0) then
   ed(nedgeint)%m1=m2
   ed(nedgeint)%t1=t2n
else
   ed(nedgeint)%t1=0
endif

if (t2.ne.0) then
   nedgeint=nedgeint+1
   ed(nedgeint)%n1=m2
   ed(nedgeint)%n2=nnodeint
   ed(nedgeint)%m1=n1
   ed(nedgeint)%m2=n2
   ed(nedgeint)%t1=t2
   ed(nedgeint)%t2=t2n
endif

do jedge=1,nedge
   if (t1.ne.0) then
      if (ed(jedge)%t1.eq.t1) then
         if (ed(jedge)%n1.eq.m1) then
            ed(jedge)%m1=nnodeint
            ed(jedge)%t1=t1n
         endif
         if (ed(jedge)%n1.eq.n2) ed(jedge)%m1=nnodeint
      endif
      if (ed(jedge)%t2.eq.t1) then
         if (ed(jedge)%n1.eq.n1) then
            ed(jedge)%m2=nnodeint
            ed(jedge)%t2=t1n
         endif
         if (ed(jedge)%n1.eq.m1) ed(jedge)%m2=nnodeint
      endif
   endif
   if (t2.ne.0) then
      if (ed(jedge)%t1.eq.t2) then
         if (ed(jedge)%n1.eq.m2) then
            ed(jedge)%m1=nnodeint
            ed(jedge)%t1=t2n
         endif
         if (ed(jedge)%n1.eq.n1) ed(jedge)%m1=nnodeint
      endif
      if (ed(jedge)%t2.eq.t2) then
         if (ed(jedge)%n1.eq.n2) then
            ed(jedge)%m2=nnodeint
            ed(jedge)%t2=t2n
         endif
         if (ed(jedge)%n1.eq.m2) ed(jedge)%m2=nnodeint
      endif
   endif
enddo

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