delaunay.f

关于网格自动生成delaunay算法的

**********************************************************************
* *                                                                * *
* *    2-dimensional Delaunay Triangulation Program                * *
* *                                                                * *
* *                                                                * *
**********************************************************************
      subroutine auto(ng,ntotal,nelm,x,ktj,mef,kbd,kcm,
     $     nodf,ibex,ibin,ibno,nex,nin,nob,nib,index,
     $     idm,mtj,jac,jnb,nei,ifix,list,istack,kv,ibr,iadres,
     $     map,angl,nedg,kstack)
      implicit real*8(a-h,o-z)
      dimension nodf(3,mef)
      dimension ibex(kbd),ibin(kbd),ibno(kbd,ktj)
      dimension index(kbd+1),idm(2*ktj+1)
      dimension mtj(2*ktj+1,3),jac(0:2*ktj+1,3)
      dimension x(ng,ktj+3)
      dimension jnb(ktj),nei(ktj,kcm),ifix(ktj),list(ktj)
      dimension istack(ktj),kv(2*ktj+1),ibr(2*ktj+1)
      dimension iadres(ktj+3),map(2*ktj+1)
      dimension angl(2*ktj+1),nedg(2*ktj+1),kstack(2*ktj+1)
      call model(nex,nin,ibex,ibin,ibno,nob,nib,index,
     $     node,x,nelm,mtj,jac,idm,ntotal,ktj,kbd,kcm,
     $     jnb,nei,ifix,list,istack,kv,ibr,iadres,map,angl,
     $     nedg,kstack)
      do i = 1,3
         do j = 1,nelm
            nodf(i,j) = mtj(j,i)        
         enddo
      enddo
      call nas(mtj,x,ktj,nelm,node)
      return
      end
*
**********************************************************************
*    Automatic Mesh Generation
**********************************************************************      
      subroutine  model(nex,nin,ibex,ibin,ibno,nob,nib,index,
     $     node,x,nelm,mtj,jac,idm,ntotal,ktj,kbd,kcm,
     $     jnb,nei,ifix,list,istack,kv,ibr,iadres,map,angl,
     $     nedg,kstack)
      implicit real*8(a-h,o-z)
      dimension ibex(kbd),ibin(kbd),ibno(kbd,ktj)
      dimension index(kbd+1),x(2,ktj+3)
      dimension jnb(ktj),nei(ktj,kcm),ifix(ktj),list(ktj)
      dimension istack(ktj),kv(2*ktj+1),ibr(2*ktj+1)
      dimension jac(0:2*ktj+1,3),idm(2*ktj+1),iadres(ktj+3)
      dimension mtj(2*ktj+1,3),map(2*ktj+1)
      dimension angl(2*ktj+1),nedg(2*ktj+1),kstack(2*ktj+1)
      do i=1,kbd+1
         index(i) = 0
      enddo
      kte=2*ktj+1
      node=0
      nelm=0
      do  i=1,kte
         idm(i) = 0
         do  j=1,3
            mtj(i,j)=0
            jac(i,j)=0
         enddo
      enddo
      do i=1,ktj
         jnb(i) = 0
         do j=1,kcm
            nei(i,j)=0
         enddo
      enddo
      do i=1,ktj
         ifix(i)=0
      enddo
      ntp = nob + nib
      xmin = x(1,1)
      xmax = xmin
      ymin = x(2,1)
      ymax = ymin
      do i = 2,ntp
         xmin = dmin1(xmin,x(1,i))
         xmax = dmax1(xmax,x(1,i))
         ymin = dmin1(ymin,x(2,i))
         ymax = dmax1(ymax,x(2,i))
      enddo
      rax = xmax-xmin
      ray = ymax-ymin
      dmax = dmax1(rax,ray)
      do i = 1,ntp
         x(1,i) = (x(1,i)-xmin)/dmax
         x(2,i) = (x(2,i)-ymin)/dmax
      enddo
      nelm = 1 
      ia = ktj+1
      ib = ktj+2
      ic = ktj+3
      do i = 1,3
      mtj(1,i) = ktj + i
      jac(1,i) = 0
      enddo
      x(1,ia) = -1.23d0
      x(2,ia) = -0.50d0
      x(1,ib) =  2.23d0
      x(2,ib) = -0.50d0
      x(1,ic) =  0.50d0
      x(2,ic) =  2.50d0
      call rough(nex,nin,ibex,ibin,ibno,ntp,nib,node,x,
     &     jnb,nei,nelm,mtj,jac,idm,list,iadres,istack,kv,
     &     ibr,map,kbd,kte,ktj,kcm)
      do i = 1,node
         ifix(i) =1
      enddo
      call fine(node,x,jnb,nei,ifix,nelm,mtj,
     &     jac,idm,istack,angl,nedg,kstack,map,kte,ktj,kcm,ntotal)
      if(nelm.ne.2*node+1) then
         write(7777,*) 'Error Message from subroutine model 1 !    '
         write(7777,*) 'incorrect number of elements'
         stop
      endif
      call laplas(node,ifix,mtj,x,jnb,nei,kte,ktj,
     &     kcm)
      call remove(nex,index,nelm,mtj,jac,idm,kbd,kte)
      call check(nelm,mtj,jac,kte)
      do i = 1,node
         x(1,i) = x(1,i) * dmax + xmin
         x(2,i) = x(2,i) * dmax + ymin
      enddo
      return 
      end
*
**********************************************************************
*    mesh generation 
*   domain surrounded by boundaries
********************************************************************** 
      subroutine bougen(iz,np,list,ntp,x,jnb,nei,nelm,mtj,jac,idm,
     &     iadres,istack,kv,ibr,map,node,kte,ktj,kcm)
      implicit real*8(a-h,o-z)
      dimension x(2,ktj+3),mtj(kte,3),jac(0:kte,3)
      dimension nei(ktj,kcm),list(ktj),iadres(ktj+3)
      dimension istack(ktj),kv(kte),ibr(kte),map(kte)
      dimension idm(kte),jnb(ktj)
      inp = 0 
      do i = 1,ktj+3
         iadres(i) = 0
      enddo 
      do i = 1,np
         iadres(list(i)) = i
      enddo
      is = list(1)
      if(is.gt.node) then 
         inp = inp+1
         call delaun(iz,is,is,ntp,x,jnb,nei,nelm,mtj,jac,
     &        idm,iadres,istack,kte,ktj,kcm)
      endif
      do 10 i = 1,np
         ip = list(mod(i,np)+1)
         iq = list(i)
         if((ip.gt.node).and.(i.ne.np)) then
            inp = inp+1
            call delaun(iz,ip,ip,ntp,x,jnb,nei,nelm,mtj,
     &           jac,idm,iadres,istack,kte,ktj,kcm)
         endif
         do j = 1,jnb(iq)
            do k = 1,3
               if(mtj(nei(iq,j),k).eq.ip) go to 10
             enddo
         enddo
         call search2(iz,ip,iq,jnb,nei,nelm,mtj,jac,idm,
     &        istack,iv,kv,iadres,ibr,kte,ktj,kcm)
         call poly(iq,ip,iv,kv,x,nelm,mtj,jac,jnb,nei,
     &        map,kte,ktj,kcm)
 10   continue
      node = node + inp
      return
      end
*
**********************************************************************
*   check for adjacent elements
**********************************************************************  
      subroutine check(nelm,mtj,jac,kte)
      implicit real*8(a-h,o-z)  
      dimension mtj(kte,3),jac(0:kte,3)
      do 10 i = 1,nelm
         do 10 j = 1,3
            ielm = i
            ia = mtj(i,j) 
            ib = mtj(i,mod(j,3)+1)
            jelm = jac(i,j)
            if(jelm.eq.0) go to 10
            do k = 1,3
               kelm = jac(jelm,k)
               if(kelm.eq.ielm) then
                  ja = mtj(jelm,k)
                  jb = mtj(jelm,mod(k,3)+1)
                  if((ia.eq.jb).and.(ib.eq.ja)) go to 10
                  write(7777,*) 'Error Message
     $                 from subroutine check 1 ! '
                  write(7777,*) 'jac is wrong'
                  stop
               endif
            enddo
            write(7777,*) 'Error Message from subroutine check 2 ! '
            write(7777,*) 'adjacent element not found'
            stop
 10   continue
      return
      end
*
**********************************************************************
*   remove triangle l from adjacent lest of point n 
*              and decrement lest size
**********************************************************************  
      subroutine decr(n,l,jnb,nei,ktj,kcm)
      implicit real*8(a-h,o-z)  
      dimension jnb(ktj),nei(ktj,kcm)
      if(n.le.ktj) then
         inb = neibor(n,l,jnb,nei,ktj,kcm)
         do i = inb,jnb(n)-1
            nei(n,i) = nei(n,i+1)
         enddo
         nei(n,jnb(n)) = 0
         jnb(n) = jnb(n)-1
         if(jnb(n).eq.0) then
            write(7777,*) 'Error Message from subroutine decr 1 ! '
            write(7777,*) 'neighbor list underflow'
            stop
         endif
      endif
      return
      end
*
**********************************************************************
*    computation of distortion rate
**********************************************************************   
      subroutine distor(n,mtj,x,drate,jedg,kte,ktj)
      implicit real*8(a-h,o-z)
      dimension mtj(kte,3),x(2,ktj+3),dst(3)
      pi = 3.1415926535d0
      xa = x(1,mtj(n,1))
      ya = x(2,mtj(n,1))
      xb = x(1,mtj(n,2))
      yb = x(2,mtj(n,2))
      xc = x(1,mtj(n,3))
      yc = x(2,mtj(n,3))
      ang1 = theta(xc,yc,xa,ya,xb,yb)
      ang2 = theta(xa,ya,xb,yb,xc,yc)
      ang3 = pi-ang1-ang2
      dst(1) = ang1-pi/3.d0
      dst(2) = ang2-pi/3.d0
      dst(3) = ang3-pi/3.d0
      drate = dabs(dst(1))+dabs(dst(2))+dabs(dst(3))
      jedg = 1
      if(dst(2).gt.dst(jedg)) jedg = 2
      if(dst(3).gt.dst(jedg)) jedg = 3
      return
      end
*
**********************************************************************
*    delaunay triangulation
**********************************************************************   
      subroutine delaun(iz,is,ig,ntp,x,jnb,nei,nelm,mtj,
     &     jac,idm,iadres,istack,kte,ktj,kcm)
      implicit real*8(a-h,o-z)
      dimension x(2,ktj+3),jnb(ktj),nei(ktj,kcm)
      dimension jac(0:kte,3),idm(kte),iadres(ktj+3)
      dimension mtj(kte,3),istack(ktj)
      itop = 0
      maxstk = ntp
      do 100 i = is,ig
         ip = i
         xp = x(1,ip)
         yp = x(2,ip)
         call locate(xp,yp,x,mtj,jac,nelm,kte,ktj,it)
         if(idm(it).ne.iz) then
            write(7777,*) 'Error Message from subroutine delaun 1 ! '
            write(7777,*) 'new point exists outside'
            stop
         endif
         ia = jac(it,1)
         ib = jac(it,2)
         ic = jac(it,3)
         iv1 = mtj(it,1)
         iv2 = mtj(it,2)
         iv3 = mtj(it,3)
         mtj(it,1) = ip
         mtj(it,2) = iv1
         mtj(it,3) = iv2
         jac(it,1) = nelm + 2
         jac(it,2) = ia
         jac(it,3) = nelm + 1
         nelm = nelm + 1
         idm(nelm) = iz
         mtj(nelm,1) = ip
         mtj(nelm,2) = iv2
         mtj(nelm,3) = iv3
         jac(nelm,1) = it
         jac(nelm,2) = ib
         jac(nelm,3) = nelm + 1
         nelm = nelm + 1
         idm(nelm) = iz
         mtj(nelm,1) = ip
         mtj(nelm,2) = iv3
         mtj(nelm,3) = iv1
         jac(nelm,1) = nelm - 1
         jac(nelm,2) = ic
         jac(nelm,3) = it
         call incr(iv1,nelm,jnb,nei,ktj,kcm)
         call incr(iv2,nelm-1,jnb,nei,ktj,kcm)
         if(iv3.le.ktj) then
            nei(iv3,neibor(iv3,it,jnb,nei,ktj,kcm)) = nelm - 1
         endif
         call incr(iv3,nelm,jnb,nei,ktj,kcm) 
         jnb(ip) = 3
         nei(ip,1) = it
         nei(ip,2) = nelm - 1
         nei(ip,3) = nelm
         if(ia.ne.0) then
            ms = iadres(iv1) * iadres(iv2)
            idf = iabs(iadres(iv1) - iadres(iv2))
            if((idm(ia).eq.iz).and.((ms.eq.0).or.(idf.ne.1))) 
     &           then
               itop = itop + 1
               istack(itop) = ipush(it,maxstk,itop)
            endif
         endif
         if(ib.ne.0) then
            call edge(ib,it,jac,kte,iedge)
            jac(ib,iedge) = nelm - 1
            ms = iadres(iv2)*iadres(iv3)
            idf = iabs(iadres(iv2) - iadres(iv3))
            if((idm(ib).eq.iz).and.((ms.eq.0).or.
     &           (idf.ne.1))) then
               itop = itop + 1
               istack(itop) = ipush(nelm-1,maxstk,itop)              
            endif
         endif
         if(ic.ne.0) then
            call edge(ic,it,jac,kte,iedge)
            jac(ic,iedge) = nelm
            ms = iadres(iv3)*iadres(iv1)
            idf = iabs(iadres(iv3)-iadres(iv1))
            if((idm(ic).eq.iz).and.((ms.eq.0).or.
     &           (idf.ne.1))) then
               itop = itop +1
               istack(itop) = ipush(nelm,maxstk,itop)
            endif
         endif
 50      if(itop.gt.0) then
            il = istack(itop)
            itop = itop-1
            ir = jac(il,2)
            call edge(ir,il,jac,kte,ierl)
            iera = mod(ierl,3) + 1
            ierb = mod(iera,3) + 1
            iv1 = mtj(ir,ierl)
            iv2 = mtj(ir,iera)
            iv3 = mtj(ir,ierb)
            call swap(x(1,iv1),x(2,iv1),x(1,iv2),x(2,iv2),x(1,iv3),
     &           x(2,iv3),xp,yp,iswap)
            if(iswap.eq.1) then
               ia = jac(ir,iera)
               ib = jac(ir,ierb)
               ic = jac(il,3)
               mtj(il,3) = iv3
               jac(il,2) = ia
               jac(il,3) = ir
               mtj(ir,1) = ip
               mtj(ir,2) = iv3
               mtj(ir,3) = iv1
               jac(ir,1) = il
               jac(ir,2) = ib
               jac(ir,3) = ic
               call decr(iv1,il,jnb,nei,ktj,kcm)
               call decr(iv2,ir,jnb,nei,ktj,kcm)
               call incr(ip,ir,jnb,nei,ktj,kcm)
               call incr(iv3,il,jnb,nei,ktj,kcm)
               if(ia.ne.0) then
                  call edge(ia,ir,jac,kte,iedge)
                  jac(ia,iedge) = il
                  ms = iadres(iv2)*iadres(iv3)
                  idf = iabs(iadres(iv2) - iadres(iv3))
                  if((idm(ia).eq.iz).and.((ms.eq.0).or.
     &                 (idf.ne.1))) then
                     itop = itop + 1
                     istack(itop) = ipush(il,maxstk,itop)                   
                  endif
               endif
               if(ib.ne.0) then