check_delaunay.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              CHECK_DELAUNAY    Apr. 2008                                     |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine check_delaunay (params,surface,is,istep,string,ref_count)

use definitions
use constants 

!------------------------------------------------------------------------------|
!(((((((((((((((( 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  ))))))))))))))))))))
!------------------------------------------------------------------------------|

implicit none

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

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

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=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'
      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), &
                    surface%xn(i1),surface%xn(i2),surface%yn(i1),surface%yn(i2),surface%zn(i1),surface%zn(i2), &
                    params%distance_exponent)
         dn1m1=dist(surface%x (i1),surface%x (j1),surface%y (i1),surface%y (j1),surface%z (i1),surface%z (j1), &
                    surface%xn(i1),surface%xn(j1),surface%yn(i1),surface%yn(j1),surface%zn(i1),surface%zn(j1), &
                    params%distance_exponent)
         dn1m2=dist(surface%x (i1),surface%x (j2),surface%y (i1),surface%y (j2),surface%z (i1),surface%z (j2), &
                    surface%xn(i1),surface%xn(j2),surface%yn(i1),surface%yn(j2),surface%zn(i1),surface%zn(j2), &
                    params%distance_exponent)
         dn2m1=dist(surface%x (i2),surface%x (j1),surface%y (i2),surface%y (j1),surface%z (i2),surface%z (j1), &
                    surface%xn(i2),surface%xn(j1),surface%yn(i2),surface%yn(j1),surface%zn(i2),surface%zn(j1), &
                    params%distance_exponent)
         dn2m2=dist(surface%x (i2),surface%x (j2),surface%y (i2),surface%y (j2),surface%z (i2),surface%z (j2), &
                    surface%xn(i2),surface%xn(j2),surface%yn(i2),surface%yn(j2),surface%zn(i2),surface%zn(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,xn1,xn2,yn1,yn2,zn1,zn2,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

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