!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! CHECK_DELAUNAY Apr. 2008 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine check_delaunay (params,surface,is,istep,string,ref_count)
use constants
use definitions
!use mpi
implicit none
include 'mpif.h'
!------------------------------------------------------------------------------|
!(((((((((((((((( 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 subroutine then updates the Delaunay triangulation around a set of points
! located on a surface in three dimensional space. A non Euclidian metrics is
! used that takes into account the divergence of normals attached to the points
! such that folding over of the surface is impossible
! The Delaunay triangumation is UPDATED, not reconstructed at each time step.
! First the delaunay triangulation is checked along each edge between two
! adjacent triangles using a generalized in-circle test based on the calculation
! of distances only; the test is based on an angle.
! where the test is positive, edge flipping takes place and all the arrays are
! permutated (icon, edge, and subsidiary arrays )
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
type (parameters) params
type (sheet) surface
integer is
integer istep
character*(*) string
integer ref_count
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer nsurfacen
integer ntriangle
integer nedge
integer iedge,jedge
integer nnmax,nflip,counter
integer ie,inode,kp,kpp
integer inodep,inodepp
integer iproc,nproc,ierr
integer i1,i2,j1,j2,k,k12,k21,n1,n2,t1,t2,nn,j,jj,err
integer,dimension(:) ,allocatable :: nedgepernode
integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
logical,dimension(:) ,allocatable :: icheck,icheckp
type (edge), dimension(:), allocatable :: ed
double precision dn1n2,dn1m1,dn1m2,dn2m1,dn2m2,dist1,dist2
double precision, external :: dist
character*72 :: shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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=err)
allocate (nedgepernode(nsurfacen),stat=err)
allocate (nodenodenumber(nnmax,nsurfacen),stat=err)
allocate (nodeedgenumber(nnmax,nsurfacen),stat=err)
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'
print *,'Node ',inode,' has ',nedgepernode(inode),'neighbors'
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) then
if (iproc.eq.0) write (8,*) 'iedge=',iedge,'nedge=',nedge
call stop_run ('pb_b in check_delaunay$')
endif
!--------------------------------------------------
! check state of triangulation with in-circle test
!--------------------------------------------------
allocate (icheck(nedge),stat=err)
allocate (icheckp(nedge),stat=err)
nflip=1
counter=0
do while (nflip > 0)
counter=counter+1
icheck=.false.
icheckp=.false.
do iedge=1+iproc,nedge,nproc
if ((ed(iedge)%t1.gt.0) .and. (ed(iedge)%t2.gt.0)) then
i1=ed(iedge)%n1
i2=ed(iedge)%n2
j1=ed(iedge)%m1
j2=ed(iedge)%m2
dn1n2=dist(surface%x(i1),surface%x(i2),surface%y(i1),surface%y(i2), &
surface%z(i1),surface%z(i2))
!dn1n2=dist(surface%x(i1),surface%x(i2),surface%y(i1),surface%y(i2), &
! surface%z(i1),surface%z(i2),params%distance_exponent)
dn1m1=dist(surface%x(i1),surface%x(j1),surface%y(i1),surface%y(j1), &
surface%z(i1),surface%z(j1))
!dn1m1=dist(surface%x(i1),surface%x(j1),surface%y(i1),surface%y(j1), &
! surface%z(i1),surface%z(j1),params%distance_exponent)
dn1m2=dist(surface%x(i1),surface%x(j2),surface%y(i1),surface%y(j2), &
surface%z(i1),surface%z(j2))
!dn1m2=dist(surface%x(i1),surface%x(j2),surface%y(i1),surface%y(j2), &
! surface%z(i1),surface%z(j2),params%distance_exponent)
dn2m1=dist(surface%x(i2),surface%x(j1),surface%y(i2),surface%y(j1), &
surface%z(i2),surface%z(j1))
!dn2m1=dist(surface%x(i2),surface%x(j1),surface%y(i2),surface%y(j1), &
! surface%z(i2),surface%z(j1),params%distance_exponent)
dn2m2=dist(surface%x(i2),surface%x(j2),surface%y(i2),surface%y(j2), &
surface%z(i2),surface%z(j2))
!dn2m2=dist(surface%x(i2),surface%x(j2),surface%y(i2),surface%y(j2), &
! surface%z(i2),surface%z(j2),params%distance_exponent)
dist1=(dn1m1**2+dn2m1**2-dn1n2**2)/(2.d0*dn1m1*dn2m1)
dist2=(dn1m2**2+dn2m2**2-dn1n2**2)/(2.d0*dn1m2*dn2m2)
dist1=acos(dist1)
dist2=acos(dist2)
icheckp(iedge) = (dist1+dist2.gt.pi) ! incircle test
endif
enddo
call mpi_allreduce (icheckp,icheck,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
nflip = count(icheck)
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i2,a,i5,a)') shift//'S.',is,':', nflip,' flips in delaunay2'
do iedge=1,nedge
if (icheck(iedge)) then
!--------------------------------------------------------
! finds where n2 is in icon(:,t1) and n1 is in icon(:,t2)
!--------------------------------------------------------
do k=1,3
if (surface%icon(k,ed(iedge)%t1).eq.ed(iedge)%n2) k12=k
if (surface%icon(k,ed(iedge)%t2).eq.ed(iedge)%n1) k21=k
enddo
!--------------------------------------------------------
! permutates nodal info in edge
!--------------------------------------------------------
n1=ed(iedge)%n1
n2=ed(iedge)%n2
ed(iedge)%n1=ed(iedge)%m2
ed(iedge)%n2=ed(iedge)%m1
ed(iedge)%m1=n1
ed(iedge)%m2=n2
!--------------------------------------------------------
! permutates nodal info in icon
!--------------------------------------------------------
t1=ed(iedge)%t1
t2=ed(iedge)%t2
surface%icon(k12,t1)=ed(iedge)%n1
surface%icon(k21,t2)=ed(iedge)%n2
!--------------------------------------------------------
! updates info in nedgepernode, nodenodenumber and nodeedgenumber
!--------------------------------------------------------
nn=ed(iedge)%m1
do j=1,nedgepernode(nn)
if (nodeedgenumber(j,nn).eq.iedge) then
do jj=j+1,nedgepernode(nn)
nodenodenumber(jj-1,nn)=nodenodenumber(jj,nn)
nodeedgenumber(jj-1,nn)=nodeedgenumber(jj,nn)
enddo
endif
enddo
nedgepernode(nn)=nedgepernode(nn)-1
nn=ed(iedge)%n1
nedgepernode(nn)=nedgepernode(nn)+1
nodenodenumber(nedgepernode(nn),nn)=ed(iedge)%n2
nodeedgenumber(nedgepernode(nn),nn)=iedge
!--------------------------------------------------------
! updates triangle info in adjacent edges
!--------------------------------------------------------
nn=ed(iedge)%n1
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%m1) then
ed(nodeedgenumber(j,nn))%t2=ed(iedge)%t1
ed(nodeedgenumber(j,nn))%m2=ed(iedge)%n2
endif
enddo
nn=ed(iedge)%m1
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%n1) then
ed(nodeedgenumber(j,nn))%t1=ed(iedge)%t1
ed(nodeedgenumber(j,nn))%m1=ed(iedge)%n2
endif
enddo
nn=ed(iedge)%n2
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%m2) then
ed(nodeedgenumber(j,nn))%t2=ed(iedge)%t2
ed(nodeedgenumber(j,nn))%m2=ed(iedge)%n1
endif
enddo
nn=ed(iedge)%m2
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%n2) then
ed(nodeedgenumber(j,nn))%t1=ed(iedge)%t2
ed(nodeedgenumber(j,nn))%m1=ed(iedge)%n1
endif
enddo
nn=ed(iedge)%n1
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%m2) then
ed(nodeedgenumber(j,nn))%m1=ed(iedge)%n2
endif
enddo
nn=ed(iedge)%m1
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%n2) then
ed(nodeedgenumber(j,nn))%m2=ed(iedge)%n1
endif
enddo
nn=ed(iedge)%n2
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%m1) then
ed(nodeedgenumber(j,nn))%m1=ed(iedge)%n1
endif
enddo
nn=ed(iedge)%m2
do j=1,nedgepernode(nn)
if (nodenodenumber(j,nn).eq.ed(iedge)%n1) then
ed(nodeedgenumber(j,nn))%m2=ed(iedge)%n2
endif
enddo
endif
enddo
if (params%debug .ge.2) call output_surf(surface,is,string,istep,ref_count)
end do
deallocate (icheck)
deallocate (icheckp)
deallocate (ed)
deallocate (nedgepernode)
deallocate (nodenodenumber)
deallocate (nodeedgenumber)
end subroutine
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
function dist (x1,x2,y1,y2,z1,z2)
!function dist (x1,x2,y1,y2,z1,z2,distance_exponent)
implicit none
double precision dist
double precision x1,y1,z1,x2,y2,z2
double precision xn1,yn1,zn1,xn2,yn2,zn2
!double precision distance_exponent
!prod=xxn(i1)*xxn(i2)+yyn(i1)*yyn(i2)+zzn(i1)*zzn(i2)
!dist=sqrt(x**2+y**2+z**2)/sqrt(prod)
!dist=sqrt(x**2+y**2+z**2)/prod
!dist=sqrt(x**2+y**2+z**2)/prod**distance_exponent
dist=sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|