delaunay.f

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

                  ms = iadres(iv3)*iadres(iv1)
                  idf = iabs(iadres(iv3) - iadres(iv1))
                  if((idm(ib).eq.iz).and.((ms.eq.0).or.
     &                 (idf.ne.1))) then
                     itop = itop + 1
                     istack(itop) = ipush(ir,maxstk,itop)                 
                  endif
               endif
               if(ic.ne.0) then
                  call edge(ic,il,jac,kte,iedge)
                  jac(ic,iedge) = ir
               endif
            endif
            goto 50
         endif
 100  continue
      return
      end
*
**********************************************************************
*   find edge in triangle l which is adjacent to triangle k
**********************************************************************  
      subroutine edge(l,k,jac,mmm,iedge)
      implicit real*8(a-h,o-z)  
      dimension jac(0:mmm,3)
      do i = 1,3
         if(jac(l,i).eq.k) then
            iedge = i
            return
         endif
      enddo
      write(7777,*) 'Error Message from subroutine edge 1 ! '
      write(7777,*) 'elements not adjacent'
      stop
      end
*
**********************************************************************
*    fine triangulation
**********************************************************************   
      subroutine fine(node,x,jnb,nei,ifix,nelm,mtj,
     &     jac,idm,istack,angl,nedg,kstack,map,kte,ktj,kcm,ntotal)
      implicit real*8(a-h,o-z)
      dimension x(2,ktj+3),jnb(ktj),nei(ktj,kcm)
      dimension jac(0:kte,3),idm(kte),angl(kte),nedg(kte)
      dimension ifix(ktj),istack(ktj),map(kte),mtj(kte,3)
      dimension kstack(kte)
      if(node.gt.ntotal) then
         write(7777,*) 'Error Message from subroutine fine 1 ! '
         write(7777,*) 'total node is less than node'
         stop
      endif
      if(node.eq.ntotal) go to 180
      almt = 1.d0/dble(ntotal)
 20   almt = 0.5d0*almt
      neib = 0
      do 30 i = 1,kte
         angl(i) = 0.d0
         nedg(i) = 0
         kstack(i) = 0
 30   continue
      do 40 i = 1,nelm
         if(idm(i).ne.0) then
            ielm = i
            neib = neib + 1
            kstack(neib) = ielm
            call distor(ielm,mtj,x,drate,jedg,kte,ktj)
            s = area(ielm,mtj,x,ktj,kte,j1,j2,j3)
            angl(ielm) = drate
            nedg(ielm) = jedg
         endif
 40   continue
      do 50 i = 1,neib-1
         k = kstack(i)
         xx = angl(k)
         ii = i + 1
         l = i
         do 60 j = ii,neib
            kk = kstack(j)
            if(angl(kk).le.xx) goto 60
            xx = angl(kk)
            l = j
 60      continue
         mm = kstack(i)
         kstack(i) = kstack(l)
         kstack(l) = mm
 50   continue
 70   if(neib.eq.0) goto 20
      if(angl(kstack(1)).lt.1.d0) goto 20
      if(node.lt.ntotal) then
         ielm = kstack(1)
         ied = nedg(ielm)
         jelm = jac(ielm,ied)
         call edge(jelm,ielm,jac,kte,jed)
         iv1 = mtj(ielm,ied)
         iv2 = mtj(ielm,mod(ied,3)+1)
         node = node + 1
         x(1,node) = (x(1,iv1)+x(1,iv2))/2.d0
         x(2,node) = (x(2,iv1)+x(2,iv2))/2.d0
         call remesh(ielm,ied,jelm,jed,node,x,jnb,nei,
     &        nelm,mtj,jac,idm,istack,kte,ktj,kcm)
         if(idm(ielm).ne.idm(jelm)) then
            ifix(node) = 1
         else
            gx = 0.d0
            gy = 0.d0
            ar = 0.d0
            do 80 j = 1,jnb(node)
               kelm = nei(node,j)
               s = area(kelm,mtj,x,ktj,kte,j1,j2,j3)
               xc = (x(1,j1)+x(1,j2)+x(1,j3))/3.d0
               yc = (x(2,j1)+x(2,j2)+x(2,j3))/3.d0
               ar = ar+s
               gx = gx+s*xc
               gy = gy+s*yc
 80         continue
            cgrax = gx/ar
            cgray = gy/ar
            xkp = x(1,node)
            ykp = x(2,node)
            x(1,node) = cgrax
            x(2,node) = cgray
            do 90 j = 1,jnb(node)
               kelm = nei(node,j)
               s = area(kelm,mtj,x,ktj,kte,j1,j2,j3)
               if(s.le.0.1d0*almt) then
                  x(1,node) = xkp
                  x(2,node) = ykp
                  goto 100
               endif
 90         continue
 100        continue
         endif
         do 110 i = 1,nelm
            map(i) = 1
 110     continue
         do 120 i = 1,jnb(node)
            lelm = nei(node,i)
            map(lelm) = 0
            angl(lelm) = 0
            nedg(lelm) = 0 
 120     continue
         nn = 0
         do 130 i = 1,neib
            if(map(kstack(i)).eq.1) then
               nn = nn+1
               kstack(nn) = kstack(i)
            endif
 130     continue
         neib = nn
         do 140 i = 1,jnb(node)
            melm = nei(node,i)
            s = area(melm,mtj,x,ktj,kte,j1,j2,j3)
            if((s.gt.almt).and.(idm(melm).ne.0)) then
               call distor(melm,mtj,x,drate,jedg,kte,ktj)
               angl(melm) = drate
               nedg(melm) = jedg
               if(drate.gt.angl(kstack(1))) then
                  do 150 j = 1,neib
                     kstack(neib-j+2) = kstack(neib-j+1)
 150              continue
                  neib = neib+1
                  kstack(1) = melm
               else if(drate.lt.angl(kstack(neib))) then
                  neib = neib + 1
                  kstack(neib) = melm
               else
                  ls = 1
                  le = neib
 160              if(le-ls.ne.1) then
                     lm = (ls+le)/2
                     if(drate.gt.angl(kstack(lm))) then
                        le = lm
                     else
                        ls = lm
                     endif
                     goto 160
                  endif
                  do 170 k = neib,le,-1
                     kstack(k+1) = kstack(k)
 170              continue
                  neib = neib+1
                  kstack(le) = melm
               endif
            endif
 140     continue
         goto 70
      endif
 180  return
      end
*
**********************************************************************
*   this program is sub-function
**********************************************************************
*   place item on ligo stack and increment stack size
**********************************************************************  
      function ipush(item,maxstk,itop)
      implicit real*8(a-h,o-z)  
      if(itop.gt.maxstk) then
         write(7777,*) 'Error Message from subroutine ipush 1 ! '
         write(7777,*) 'stack overflow'
         stop
      endif
      ipush = item
      return
      end
*
**********************************************************************
*  find triangle l in adjacent list of point n
**********************************************************************  
      function neibor(n,l,jnb,nei,ktj,kcm)
      implicit real*8(a-h,o-z)  
      dimension jnb(ktj),nei(ktj,kcm)
      do i = 1,jnb(n)
         if(nei(n,i).eq.l) then
            neibor = i
            return
         endif
      enddo
      write(7777,*) 'Error Message from subroutine neibor 1 ! '
      write(7777,*) 'elements not adjacent'
      stop
      end
*
**********************************************************************
*  find triangle l in adjacent list of point n
**********************************************************************  
      function ivert(l,k,mtj,mmm)
      implicit real*8(a-h,o-z)      
      dimension mtj(mmm,3)
      do i = 1,3
         if(mtj(l,i).eq.k) then
            ivert = i
            return
         endif
      enddo
      write(7777,*) 'Error Message from subroutine ivert 1 ! '
      write(7777,*) 'vertices not includes'
      stop
      end
*
**********************************************************************
*  computation of area
**********************************************************************  
      function area(ielm,mtj,x,ktj,kte,j1,j2,j3)
      implicit real*8(a-h,o-z)
      dimension mtj(kte,3),x(2,ktj+3)
      j1 = mtj(ielm,1)
      j2 = mtj(ielm,2)
      j3 = mtj(ielm,3)
      area = 0.5d0*(x(1,j1)*x(2,j2)+x(1,j2)*x(2,j3)+x(1,j3)*x(2,j1)
     &     -x(1,j1)*x(2,j3)-x(1,j2)*x(2,j1)-x(1,j3)*x(2,j2))
      return
      end
*
**********************************************************************
*  computation of angle
**********************************************************************  
      function theta(x0,y0,x1,y1,x2,y2)
      implicit real*8(a-h,o-z)
      pi = 3.141592653589793d0
      error = 1.0d-15
      xa = x1-x0
      ya = y1-y0
      xb = x2-x0
      yb = y2-y0
      prdin = xa*xb+ya*yb
      prdex = xa*yb-xb*ya
      if(dabs(prdin).lt.error) then
         theta = pi/2.d0
      else
         theta = datan(prdex/prdin)
         if(theta.lt.0) then
            theta = theta+pi
         endif
         if(theta.gt.pi) then
            theta = theta-pi
         endif
      endif
      return
      end
*
**********************************************************************
*   place triangle l on adjacent list of point n
*        and increment list size
**********************************************************************  
      subroutine incr(n,l,jnb,nei,ktj,kcm)
      implicit real*8(a-h,o-z)  
      dimension jnb(ktj),nei(ktj,kcm)
      if(n.le.ktj) then
         jnb(n) = jnb(n)+1
         if(jnb(n).gt.kcm) then
            write(7777,*) 'Error Message from subroutine incr 1 ! '
            write(7777,*) 'neighbor list overflow'
            stop
         endif
         nei(n,jnb(n)) = l
      endif
      return
      end
*
**********************************************************************
*    laplacian method
**********************************************************************   
      subroutine laplas(node,ifix,mtj,x,jnb,nei,
     &     kte,ktj,kcm)
      implicit real*8(a-h,o-z)
      dimension ifix(ktj),mtj(kte,3),x(2,ktj+3)
      dimension nei(ktj,kcm),jnb(ktj)
      itera = 5
      do 10 it = 1,itera
         do 20 i = 1,node
            if(ifix(i).eq.0) then
               gx = 0.d0
               gy = 0.d0
               ar = 0.d0
               do 30 j = 1,jnb(i)
                  ielm = nei(i,j)
                  j1 = mtj(ielm,1)
                  j2 = mtj(ielm,2)
                  j3 = mtj(ielm,3)
                  s = area(ielm,mtj,x,ktj,kte,j1,j2,j3)
                  xc = (x(1,j1)+x(1,j2)+x(1,j3))/3.d0
                  yc = (x(2,j1)+x(2,j2)+x(2,j3))/3.d0
                  ar = ar+s
                  gx = gx+s*xc
                  gy = gy+s*yc
 30            continue
               cgrax = gx/ar
               cgray = gy/ar
               x(1,i) = cgrax
               x(2,i) = cgray
            endif
 20      continue
 10   continue
      return
      end
*
**********************************************************************
*  apply lawson's swapping algorithm
**********************************************************************  
      subroutine lawson(nte,ien,jee,npl,x,jstack,ktj,lte,ltj)
      implicit real*8(a-h,o-z)  
      dimension ien(lte,3),jee(lte,3),x(2,ktj+3)
      dimension jstack(ltj)
      itop = 0 
      maxstk = npl
      ncount = 0 
      do i = 1,nte
         ielm = 1
         itop = itop+1
         jstack(itop) = ipush(ielm,maxstk,itop)
      enddo
 10   if(itop.gt.0) then
         ncount = ncount + 1
         if(ncount.ge.lte) then
            write(7777,*) 'Error Message from subroutine lawson 1 ! '
            write(7777,*) 'non-convergence'
            stop
         endif
         il = jstack(itop)
         itop = itop - 1
         do 20 j = 1,3
            jl1 = j
            jl2 = mod(jl1,3) + 1
            jl3 = mod(jl2,3) + 1
            ir = jee(il,jl1)
            if(ir.eq.0) go to 20
            iv1 = ien(il,jl1)
            iv2 = ien(il,jl2)
            iv3 = ien(il,jl3)
            xx = x(1,ien(il,jl3))
            yy = x(2,ien(il,jl3))
            call edge(ir,il,jee,lte,jr1)
            jr2 = mod(jr1,3) + 1
            jr3 = mod(jr2,3) + 1
            iv4 = ien(ir,jr3)
            call swap(x(1,iv2),x(2,iv2),x(1,iv1),x(2,iv1),x(1,iv4),
     &           x(2,iv4),xx,yy,iswap)
            if(iswap.eq.1) then
            ia = jee(il,jl2)
            ib = jee(ir,jr2)
            ien(il,jl2) = iv4
            jee(il,jl1) = ib
            jee(il,jl2) = ir