delaunay.f

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

            ien(ir,jr2) = iv3
            jee(ir,jr1) = ia
            jee(ir,jr2) = il
            if(ia.ne.0) then
               call edge(ia,il,jee,lte,iedge)
               jee(ia,iedge) = ir
            endif
            if(ib.ne.0) then
               call edge(ib,ir,jee,lte,iedge)
               jee(ib,iedge) = il
            endif
            itop = itop+1
            jstack(itop) = ipush(il,maxstk,itop)
            go to 10
         endif
 20   continue
      go to 10
      endif
      return
      end
*
**********************************************************************
*   locate triangle which encloses point (xp,yp)
**********************************************************************  
      subroutine locate(xp,yp,x,mtj,jac,nelm,
     &     kte,ktj,itri)
      implicit real*8(a-h,o-z)  
      dimension x(2,ktj+3),mtj(kte,3),jac(0:kte,3)
      itri = nelm
 10   continue
      do i = 1,3
         ia = mtj(itri,i)
         ib = mtj(itri,mod(i,3)+1)
         if((x(2,ia)-yp)*(x(1,ib)-xp).gt.
     &        (x(1,ia)-xp)*(x(2,ib)-yp)) then
            itri = jac(itri,i)
            go to 10
         endif
      enddo
      return
      end
*
**********************************************************************
*    pick up points on boundary of given polygon
**********************************************************************  
      subroutine pick(iq,ip,iv,kv,mtj,jac,map,npa,npb,nsra,
     &     nsrb,kte,ltj)
      implicit real*8(a-h,o-z)      
      dimension kv(kte),mtj(kte,3),jac(0:kte,3),map(kte)
      dimension nsra(ltj),nsrb(ltj)
      npa = 1
      nsra(npa) = ip
      ivx = ivert(kv(1),ip,mtj,kte)
      npa = npa+1
      nsra(npa) = mtj(kv(1),mod(ivx,3)+1)
      do i = 2,iv-1
         jvx = ivert(kv(i),nsra(npa),mtj,kte)
         jelm = jac(kv(i),jvx)
         if(jelm.eq.0) then
            npa = npa + 1
            nsra(npa) = mtj(kv(i),mod(jvx,3)+1)
         elseif(map(jelm).eq.0) then
            npa = npa+1
            nsra(npa) = mtj(kv(i),mod(jvx,3)+1)
         endif
      enddo
      npa = npa + 1
      nsra(npa) = i
      npb = 1
      nsrb(npb) = iq
      ivx = ivert(kv(iv),iq,mtj,kte)
      npb = npb+1
      nsrb(npb) = mtj(kv(iv),mod(ivx,3)+1)
      do i = iv-1,2,-1
         jvx = ivert(kv(i),nsrb(npb),mtj,kte)
         jelm = jac(kv(i),jvx)
         if(jelm.eq.0) then
            npb = npb + 1
            nsrb(npb) = mtj(kv(i),mod(jvx,3)+1)
         elseif(map(jelm).eq.0) then
            npb = npb+1
            nsrb(npb) = mtj(kv(i),mod(jvx,3)+1)
         endif
      enddo
      npb = npb + 1
      nsrb(npb) = ip
      if(iv.ne.npa+npb-4) then
         write(7777,*) 'Error Message from subroutine pick 1 ! '
         write(7777,*) 'incorrect number of nodes'
         stop 
      endif
      return
      end
*
**********************************************************************
*  subdivide given polygon into finte element
**********************************************************************  
      subroutine poly(iq,ip,iv,kv,x,nelm,mtj,jac,jnb,
     &     nei,map,kte,ktj,kcm)
      implicit real*8(a-h,o-z)  
      parameter(lte=100)
      dimension kv(kte),x(2,ktj+3),mtj(kte,3)
      dimension jnb(ktj),nei(ktj,kcm),map(kte),nsra(lte+2)
      dimension iena(lte,3),ienb(lte,3),jeea(lte,3)
      dimension ihen(2*lte+1,2),jhen(2*lte+1),iad(2*lte+1)
      dimension jac(0:kte,3),nsrb(lte+2),jeeb(lte,3)
      dimension jstack(lte+2)
      if(iv.eq.2) then
         call edge(kv(1),kv(2),jac,kte,ia)
         call edge(kv(2),kv(1),jac,kte,ja)
         ir = jac(kv(1),mod(ia,3)+1) 
         jr = jac(kv(2),mod(ja,3)+1) 
         mtj(kv(1),mod(ia,3)+1) = iq
         jac(kv(1),ia) = jr
         jac(kv(1),mod(ia,3)+1) = kv(2)
         mtj(kv(2),mod(ja,3)+1) = ip
         jac(kv(2),ja) = ir
         jac(kv(2),mod(ja,3)+1) = kv(1)
         if(ir.ne.0) then
            call edge(ir,kv(1),jac,kte,iedge)
            jac(ir,iedge) = kv(2)
         endif
         if(jr.ne.0) then
            call edge(jr,kv(2),jac,kte,iedge)
            jac(jr,iedge) = kv(1)
         endif
         iv1 = mtj(kv(1),ia)
         iv2 = mtj(kv(2),ja)
         call decr(iv1,kv(2),jnb,nei,ktj,kcm)
         call decr(iv2,kv(1),jnb,nei,ktj,kcm)
         call incr(iq,kv(1),jnb,nei,ktj,kcm)
         call incr(ip,kv(2),jnb,nei,ktj,kcm)
      else
         ltj = lte+2
         lhn = 2*lte+1
         npa = 0
         npb = 0
         do i =1,ltj
            nsra(i) = 0
            nsrb(i) = 0
         enddo
         do i = 1,nelm
            map(i) = 0 
         enddo
         do i = 1,iv
            map(kv(i)) = 1
         enddo
         do i = 1,iv
            ielm = kv(i)
            do j = 1,3
               ivx = mtj(ielm,j)
               call decr(ivx,ielm,jnb,nei,ktj,kcm)
            enddo
         enddo
         call pick(iq,ip,iv,kv,mtj,jac,map,npa,npb,nsra,nsrb,
     &        kte,ltj)
         call subdiv(npa,nsra,x,nta,iena,jeea,ihen,jhen,
     &        iad,jstack,ktj,lte,ltj,lhn)
         call subdiv(npb,nsrb,x,ntb,ienb,jeeb,ihen,jhen,
     &        iad,jstack,ktj,lte,ltj,lhn)
         if(iv.ne.nta+ntb) then
            write(7777,*) 'Error Message from subroutine poly 1 ! '
            write(7777,*) 'incorrect number of elements formed'
            stop
         endif
         do i = 1,iv
            do j = 1,3
               jac(kv(i),j) = 0
            enddo
         enddo
         do i = 1,nta
            ielm = kv(i)
            do j = 1,3
               mtj(ielm,j) = iena(i,j)
               if(jeea(i,j).ne.0) then
                  jac(ielm,j) = kv(jeea(i,j))
               endif
            enddo
         enddo
         do i = 1,ntb
            ielm = kv(nta+i)
            do j = 1,3
               mtj(ielm,j) = ienb(i,j)
               if(jeeb(i,j).ne.0) then
                  jac(ielm,j) = kv(nta+jeeb(i,j))
               endif
            enddo
         enddo
         do i = 1,iv
            ielm = kv(i)
            do j = 1,3
               ivx = mtj(ielm,j)
               call incr(ivx,ielm,jnb,nei,ktj,kcm)
            enddo
         enddo

         do i = 1,iv
            ielm = kv(i)
            do 10 j = 1,3
               jelm = jac(ielm,j)
               if(jelm.ne.0) go to 10
               ips = mtj(ielm,j)
               ipg = mtj(ielm,mod(j,3)+1)
               do k =1,jnb(ipg)
                  kelm = nei(ipg,k)
                  do l = 1,3
                     iva = mtj(kelm,l)
                     ivb = mtj(kelm,mod(l,3)+1)
                     if((iva.eq.ipg).and.(ivb.eq.ips)) then
                        jac(ielm,j) = kelm
                        jac(kelm,l) = ielm
                        go to 10
                     endif
                  enddo
               enddo
 10         continue
         enddo
      endif
      return
      end
*
**********************************************************************
*    subdivide given domain by new data point
**********************************************************************   
      subroutine remesh(ielm,ied,jelm,jed,node,x,jnb,
     &     nei,nelm,mtj,jac,idm,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),istack(ktj),mtj(kte,3)
      itop = 0
      maxstk = node
      ip = node
      xp = x(1,ip)
      yp = x(2,ip)
      ibj = mod(ied,3) + 1
      jbj = mod(jed,3) + 1
      ia = jac(ielm,ibj) 
      ib = jac(ielm,mod(ibj,3)+1)
      ic = jac(jelm,jbj)
      id = jac(jelm,mod(jbj,3)+1)
      iv1 = mtj(ielm,ibj)
      iv2 = mtj(ielm,mod(ibj,3)+1)
      iv3 = mtj(jelm,jbj)
      iv4 = mtj(jelm,mod(jbj,3)+1)
      mtj(ielm,1) = ip
      mtj(ielm,2) = iv1
      mtj(ielm,3) = iv2
      jac(ielm,1) = nelm+2
      jac(ielm,2) = ia
      jac(ielm,3) = nelm+1
      nelm = nelm+1
      idm(nelm) = idm(ielm)
      mtj(nelm,1) = ip
      mtj(nelm,2) = iv2
      mtj(nelm,3) = iv3
      jac(nelm,1) = ielm
      jac(nelm,2) = ib
      jac(nelm,3) = jelm
      mtj(jelm,1) = ip
      mtj(jelm,2) = iv3
      mtj(jelm,3) = iv4
      jac(jelm,1) = nelm
      jac(jelm,2) = ic
      jac(jelm,3) = nelm+1
      nelm = nelm+1
      idm(nelm) = idm(jelm)
      mtj(nelm,1) = ip
      mtj(nelm,2) = iv4
      mtj(nelm,3) = iv1
      jac(nelm,1) = jelm
      jac(nelm,2) = id
      jac(nelm,3) = ielm
      if(iv1.le.ktj) then
         nei(iv1,neibor(iv1,jelm,jnb,nei,ktj,kcm)) = nelm
      endif
      call incr(iv2,nelm-1,jnb,nei,ktj,kcm)
      if(iv3.le.ktj) then
         nei(iv3,neibor(iv3,ielm,jnb,nei,ktj,kcm)) = nelm-1
      endif
      call incr(iv4,nelm,jnb,nei,ktj,kcm)
      jnb(ip) = 4
      nei(ip,1) = ielm
      nei(ip,2) = nelm-1
      nei(ip,3) = jelm
      nei(ip,4) = nelm
      if(ia.ne.0) then
         if(idm(ia).eq.idm(ielm)) then
            itop = itop + 1
            istack(itop) = ipush(ielm,maxstk,itop)
         endif
      endif
      if(ib.ne.0) then
         call edge(ib,ielm,jac,kte,iedge)
         jac(ib,iedge) = nelm-1
         if(idm(ib).eq.idm(nelm-1)) then
            itop = itop + 1
            istack(itop) = ipush(nelm-1,maxstk,itop)
         endif
      endif
      if(ic.ne.0) then
         if(idm(ic).eq.idm(jelm)) then
            itop = itop +1
            istack(itop) = ipush(jelm,maxstk,itop)
         endif
      endif
      if(id.ne.0) then
         call edge(id,jelm,jac,kte,iedge)
         jac(id,iedge) = nelm
         if(idm(id).eq.idm(nelm)) then
            itop = itop + 1
            istack(itop) = ipush(nelm,maxstk,itop)
         endif
      endif
 10   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
               if(idm(ia).eq.idm(il)) then
                  itop = itop + 1
                  istack(itop) = ipush(il,maxstk,itop)
               endif
            endif
            if(ib.ne.0) then
               if(idm(ib).eq.idm(ir)) 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
         go to 10
      endif
      return
      end
*
**********************************************************************
*   remove alltriangles outside of boundary
**********************************************************************  
      subroutine remove(nex,index,nelm,mtj,jac,idm,kbd,kte)
      implicit real*8(a-h,o-z)  
      dimension index(kbd+1),mtj(kte,3),jac(0:kte,3),idm(kte)
      dimension mkp(3),jkp(3)
      ielm = 0
      index(1) = 1
      do i = 1,nex
         iz = i
         inelm = 0
         do j = index(i),nelm
            if(idm(j).eq.iz) then
               ielm = ielm+1
               jelm = j
               inelm = inelm+1
               if(ielm.ne.jelm) then