delaunay.f

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

                  do k = 1,3
                     mkp(k) = mtj(ielm,k)
                     jkp(k) = jac(ielm,k)
                  enddo
                  ikp = idm(ielm)
                  do k = 1,3
                     kelm = jac(jelm,k)
                     if(kelm.ne.0) then
                        call edge(kelm,jelm,jac,kte,kedg)
                        jac(kelm,kedg) = ielm+kte+1
                     endif
                  enddo
                  do l = 1,3
                     lelm = jkp(l)
                     if(lelm.ne.0) then
                        call edge(lelm,ielm,jac,kte,ledg)
                        jac(lelm,ledg) = jelm+kte+1
                     endif
                  enddo
                  do k = 1,3
                     jac(ielm,k) = mod(jac(ielm,k),kte+1)
                     jac(jelm,k) = mod(jac(jelm,k),kte+1)
                     kelm = jac(ielm,k)
                     lelm = jac(jelm,k)
                     do l = 1,3
                        jac(kelm,l) = mod(jac(kelm,l),kte+1)
                        jac(lelm,l) = mod(jac(lelm,l),kte+1)
                     enddo
                  enddo
                  do k = 1,3
                     jkp(k) = jac(ielm,k)
                  enddo
                  do k = 1,3
                     mtj(ielm,k) = mtj(jelm,k)
                     jac(ielm,k) = jac(jelm,k)
                     mtj(jelm,k) = mkp(k)
                     jac(jelm,k) = jkp(k)
                  enddo
                  idm(ielm) = idm(jelm)
                  idm(jelm) = ikp
               endif
            endif
         enddo
         index(i+1) = index(i)+inelm
      enddo
      do i = 1,ielm
         do j = 1,3
            if(jac(i,j).gt.ielm) then
               jac(i,j) = 0
            endif
         enddo
      enddo
      do i = ielm+1,nelm
         do j = 1,3
            mtj(i,j) = 0
            jac(i,j) = 0
         enddo
      enddo
      nelm = ielm
      return
      end
*
**********************************************************************
*    Rough triangulation
********************************************************************** 
      subroutine 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)
      implicit real*8(a-h,o-z)
      dimension ibex(kbd),ibin(kbd),ibno(kbd,ktj)
      dimension x(2,ktj+3)
      dimension mtj(kte,3),jac(0:kte,3),idm(kte)
      dimension list(ktj),iadres(ktj+3),istack(ktj),jnb(ktj)
      dimension map(kte),nei(ktj,kcm),ibr(kte),kv(kte)
      do 10 i = 1,nex
         iz = 0
         nb = i
         np = ibex(nb)
         do 20 j = 1,np
            list(j) = ibno(nb,j)
 20      continue
         call bougen(iz,np,list,ntp,x,jnb,nei,nelm,mtj,
     &        jac,idm,iadres,istack,kv,ibr,map,node,
     &        kte,ktj,kcm)
         do 30 k = 1,nelm
            if(idm(k).eq.iz) then
               ma = iadres(mtj(k,1))
               mb = iadres(mtj(k,2))
               mc = iadres(mtj(k,3))
               ms = ma*mb*mc
               mp = (mb-ma)*(mc-mb)*(ma-mc)
               if((ms.ne.0).and.(mp.gt.0)) then
                  idm(k) = nb
               endif
            endif
 30      continue
 10   continue
      do 50 i = 1,nin
         nb = i
         np = ibin(nb)
         do 60 j = 1,np
            list(j) = node + j
 60      continue
         is = list(1)
         xs = x(1,is)
         ys = x(2,is)
         call locate(xs,ys,x,mtj,jac,nelm,kte,ktj,loc)
         iz = idm(loc)
         call bougen(iz,np,list,ntp,x,jnb,nei,nelm,mtj,
     &        jac,idm,iadres,istack,kv,ibr,map,node,
     &        kte,ktj,kcm)
         do 70 k = 1,nelm
            if(idm(k).eq.iz) then
               ma = iadres(mtj(k,1))
               mb = iadres(mtj(k,2))
               mc = iadres(mtj(k,3))
               ms = ma*mb*mc
               mp = (mb-ma)*(mc-mb)*(ma-mc)
               if((ms.ne.0).and.(mp.lt.0)) then
                  idm(k) = 0
               endif
            endif
 70      continue
 50   continue
      do i = 1,ktj+3
         iadres(i) = 0
      enddo
      do 100 i = 1,nib
         node = node+1
         xa = x(1,node)
         ya = x(2,node)
         call locate(xa,ya,x,mtj,jac,nelm,kte,ktj,it)
         iz = idm(it)
         call delaun(iz,node,node,ntp,x,jnb,nei,nelm,mtj,
     &        jac,idm,iadres,istack,kte,ktj,kcm)
 100  continue
      if(node.ne.ntp) then
         write(7777,*) 'Error Message from subroutine rough 1 ! '
         write(7777,*) 'incorrect number of nodes'
         stop
      endif
      return
      end
*
**********************************************************************
*   search for treangules to be modified
**********************************************************************  
      subroutine search2(iz,ip,iq,jnb,nei,nelm,mtj,jac,idm,
     &     istack,iv,kv,iadres,ibr,kte,ktj,kcm)
      implicit real*8(a-h,o-z)  
      dimension mtj(kte,3),jac(0:kte,3),idm(kte),jnb(ktj)
      dimension istack(ktj),kv(kte),iadres(ktj+3),ibr(kte)
      dimension nei(ktj,kcm)
      iv = 0
      mstk = 0 
      nbr = 0
      do i = 1,nelm
         kv(i) = 0
         ibr(i) = 0
      enddo
      do i = 1,ktj
         istack(i) = 0
      enddo 
      do i = 1,jnb(iq)
         nbr = nbr + 1
         ibr(nei(iq,i)) = nbr
      enddo
      do 10 i = 1,jnb(iq)
         ielm = nei(iq,i)
         j = ivert(ielm,iq,mtj,kte)
         ja = mod(j,3) + 1
         jb = mod(ja,3) + 1
         ia = mtj(ielm,ja)
         ib = mtj(ielm,jb)
         idf = iabs(iadres(ia) - iadres(ib))
         ms = iadres(ia)*iadres(ib)
         if((idf.eq.1).and.(ms.ne.0)) go to 10
         jelm = jac(ielm,ja)
         if(jelm.eq.0) go to 10
         if(idm(jelm).ne.iz) go to 10
         nbr = nbr + 1
         ibr(jelm) = nbr
         if(mtj(jelm,1).eq.ip) go to 40
         if(mtj(jelm,2).eq.ip) go to 40
         if(mtj(jelm,3).eq.ip) go to 40
         mstk = mstk + 1
         istack(mstk) = jelm
 10   continue
 20   if(mstk.eq.0) then
         write(7777,*) 'Error Message from subroutine search2 1 ! '
         write(7777,*) 'istack is empty'
         stop
      endif
      ielm = istack(1)
      mstk = mstk - 1
      do i = 1,mstk
         istack(i) = istack(i+1)
      enddo
      istack(mstk+1) = 0
      do 30 j =1,3
         ja = mod(j,3) + 1
         ia = mtj(ielm,j)
         ib = mtj(ielm,ja)
         idf = iabs(iadres(ia)-iadres(ib))
         ms = iadres(ia)*iadres(ib)
         if((idf.eq.1).and.(ms.ne.0)) go to 30
         jelm = jac(ielm,j)
         if(jelm.eq.0) go to 30
         if(ibr(jelm).ne.0) go to 30
         if(idm(jelm).ne.iz) go to 30
         nbr = nbr + 1
         ibr(jelm) = nbr
         if(mtj(jelm,1).eq.ip) go to 40
         if(mtj(jelm,2).eq.ip) go to 40
         if(mtj(jelm,3).eq.ip) go to 40
         mstk = mstk + 1
         istack(mstk) = jelm
 30   continue
      go to 20
 40   iv = iv + 1
      kv(iv) = jelm
      if(ibr(jelm).le.jnb(iq)) go to 60
      min = kte + 1
      do 50 j = 1,3
         jr = jac(jelm,j)
         if(jr.eq.0) go to 50
         if(ibr(jr).eq.0) go to 50
         ja = mod(j,3) + 1
         ia = mtj(jelm,j)
         ib = mtj(jelm,ja)
         idf = iabs(iadres(ia) - iadres(ib))
         ms = iadres(ia) * iadres(ib)
         if((idf.eq.1).and.(ms.ne.0)) go to 50
         if(ibr(jr).lt.min) then
            kelm = jr
            min = ibr(jr)
         endif
 50      continue
         if(min.eq.kte+1) then
            write(7777,*) 'Error Message from subroutine search2 2 ! '
            write(7777,*) 'branch is closed'
            stop
         endif
         jelm = kelm
         go to 40
 60      continue
         return
         end
*
**********************************************************************
*  subdivide given poltgon using by the modified-delaunay
**********************************************************************  
      subroutine subdiv(npl,nsr,x,nte,ien,jee,ihen,jhen,
     &     iad,jstack,ktj,lte,ltj,lhn)
      implicit real*8(a-h,o-z)  
      dimension nsr(ltj),x(2,ktj+3),ien(lte,3)
      dimension ihen(lhn,2),jhen(lhn),iad(lhn),jstack(ltj)
      dimension jee(lte,3)
      nte = 0 
      nnpl = npl
      do i = 1,lte
         do j = 1,3
            ien(i,j) = 0
            jee(i,j) = 0
         enddo
      enddo
      icheck = 0
 10   if(nnpl.ge.3) then
         icheck = icheck + 1
         if (icheck.ge.10000) then
            write(7777,*) 'Error Message from subroutine subdiv 1 ! '
            write(7777,*) 'rooped in computing ien'
            stop
         endif
         nbs = 1
 20      if(nbs.le.nnpl-1) then
            ia = nsr(nbs)
            ib = nsr(nbs+1)
            ic = nsr(mod(nbs+1,nnpl)+1)
            xa = x(1,ib)-x(1,ia)
            ya = x(2,ib)-x(2,ia)
            xb = x(1,ic)-x(1,ia)
            yb = x(2,ic)-x(2,ia)
            see = xa*yb-xb*ya
            if(see.gt.1.0d-15) then
               nte = nte+1
               ien(nte,1) = ia
               ien(nte,2) = ib
               ien(nte,3) = ic
               nnpl = nnpl-1
               do i = nbs+1,nnpl
                  nsr(i) = nsr(i+1) 
               enddo
            endif
            nbs = nbs+1
            go to 20
         endif
         go to 10
      endif
      ix = 0
      do i = 1,lhn
         ihen(i,1) = 0
         ihen(i,2) = 0
         jhen(i) = 0
         iad(i) = 0
      enddo
      do i = 1,nte
         do 30 j = 1,3
            ia = ien(i,j)
            ib = ien(i,mod(j,3)+1)
            do k = 1,ix
               if((ihen(k,1).eq.ib).and.(ihen(k,2).eq.ia)) then
                  jee(i,j) = jhen(k)
                  jee(jhen(k),iad(k)) = i
                  go to 30
               endif
            enddo
            ix = ix+1
            jhen(ix) = i
            iad(ix) = j
            ihen(ix,1) = ia 
            ihen(ix,2) = ib
 30      continue
      enddo
      call lawson(nte,ien,jee,npl,x,jstack,ktj,lte,ltj)
      return
      end
*
**********************************************************************
*   check if point lies inside the circumcircle
**********************************************************************  
      subroutine swap(x1,y1,x2,y2,x3,y3,xp,yp,iswap)
      implicit real*8(a-h,o-z)  
      x13 = x1-x3
      y13 = y1-y3
      x23 = x2-x3
      y23 = y2-y3
      x1p = x1-xp
      y1p = y1-yp
      x2p = x2-xp
      y2p = y2-yp
      cosa = x13*x23+y13*y23
      cosb = x2p*x1p+y1p*y2p
      if((cosa.ge.0.d0).and.(cosb.ge.0.d0)) then
         iswap = 0
      elseif((cosa.lt.0.d0).and.(cosb.lt.0.d0)) then
         iswap = 1
      else
         sina = x13*y23-x23*y13
         sinb = x2p*y1p-x1p*y2p
         if((sina*cosb+sinb*cosa).lt.0.d0) then
            iswap = 1
         else
            iswap = 0
         endif
      endif
      return
      end  
*
************************************************************
*     This program is used for translating a format of     *
*       dlaunay triangulation to that of nasmarc           *
************************************************************
      subroutine nas(mtj,x,ktj,nelm,node)
      implicit real*8(a-h,o-z)
      parameter(ii=1) 
      dimension mtj(2*ktj+1,4)
      dimension x(2,ktj+3)
      open(9,file='MESH.NAS',status='unknown')
      do i = 1,node
         write(9,6000) i,x(1,i),x(2,i)
      write(9,6010)
      enddo
      do i = 1,nelm
         write(9,6100) i,ii,(mtj(i,j),j=1,3)
      enddo
      write(9,6200)
      close(9)
 6000 format('GRID*           ',i8,'               0',f16.10,f17.10)
 6010 format('*                     0.               0                ')
 6100 format('CTRIA3    ',i6,4i8)
 6200 format('ENDDATA')
      end