!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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,i8,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,i8,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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|