refine_surface.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              REFINE_SURFACE    May 2008                                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
      
subroutine refine_surface (params,surface,threadinfo,nadd,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 threads      
use definitions

implicit none

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

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

integer naddp
integer nnmax
type (edge), dimension(:), allocatable :: ed
type (edge), dimension(:), allocatable :: edswap
integer,dimension(:)  ,allocatable :: nedgepernode
integer,dimension(:)  ,allocatable :: refine_list
integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
logical,dimension(:),allocatable::refine
logical,dimension(:),allocatable::refinep
character*72  :: shift
integer ie,j,jedge,kp,kpp,k,inode,inodep,inodepp,iadd
integer ntriangle
integer nedge
integer nedgen,nsurfacen
double precision dist1,distmax,prod
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

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

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
         
shift=' '

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

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

      
allocate (ed(nedge),stat=threadinfo%err)                       ; call heap (threadinfo,'ed','refine_surface',size(ed),'int',+1)
allocate (nedgepernode(nsurfacen),stat=threadinfo%err)
call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',+1)
allocate (nodenodenumber(nnmax,nsurfacen),stat=threadinfo%err)
call heap (threadinfo,'nodenodenumber','refine_surface',size(nodenodenumber),'int',+1)
allocate (nodeedgenumber(nnmax,nsurfacen),stat=threadinfo%err)
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

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

refine=.false.
refinep=.false.

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), &
              surface%xn(i1),surface%xn(i2),surface%yn(i1),surface%yn(i2),surface%zn(i1),surface%zn(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 ) 
enddo

call mpi_allreduce (refinep,refine,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)


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

nadd=count(refine)

allocate(refine_list(nadd))

iadd=0
naddp=0
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
enddo

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

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 (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 (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$')

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

deallocate(refine_list)

call heap (threadinfo,'ed','refine_surface',size(ed),'int',-1) ; deallocate (ed)         
call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',-1) ; deallocate (nedgepernode)
call heap (threadinfo,'nodenodenumber','refine_surface',size(nodenodenumber),'int',-1) ; deallocate (nodenodenumber)
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|