search.f

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

**********************************************************************
*     Making node search table
**********************************************************************
      subroutine mkstbl(mef,msrch,npoin,nef,
     $     nodf,nadrs,isrch)
      implicit real*8(a-h,o-z)
      integer    npoin,nef
      dimension  nodf(3,mef)
      dimension  isrch(2,msrch)
      nnode= 3
      do 100 i= 1, npoin
         isrch(1,i)= i
         isrch(2,i)= 0
 100  continue
      iadress= npoin+1 
      do 200 ief= 1,nef
         do 210 i= 1,nnode
            itop= nodf(i,ief)
            do 220 j= 1,nnode
               itail= nodf(j,ief)
               ifrag= itop
 250           continue
               if(isrch(1,ifrag).eq.itail) then
                  goto 260
               endif
               if(isrch(2,ifrag).eq.0) then
                  isrch(2,ifrag)= iadress
                  isrch(1,iadress)= itail
                  isrch(2,iadress)= 0
                  iadress= iadress+1
                  if(iadress.gt.msrch) then
                     write(6,900) iadress, msrch
                  endif
                  goto 260
               else
                  ifrag= isrch(2,ifrag)
               endif
               goto 250
 260           continue
 220        continue
 210     continue
 200  continue
      nadrs= iadress-1
 900  format(/5x,'** Message from subrouthine mkstbl **',
     $       /5x,'   iadress is larger than msrch !       ',
     $       /5x,'   iadress =',i10,
     $       /5x,'   msrch=   ',i10)
      return
      end
*
**********************************************************************
*     searching nodes from node search table
**********************************************************************
      subroutine search(msrch,npoin,ntgpt,nadrs,isrch,mtail,ltail)
      implicit real*8(a-h,o-z)
      integer    npoin, ntgpt
      dimension  isrch(2,msrch), ltail(mtail)
      if(ntgpt.gt.npoin) then
         write(7777,*) 'Error Message from subroutine search 1 !      '
         write(7777,*) 'number of a target point is out of the range ! '
         write(7777,*) 'ntgpt =',ntgpt
         write(7777,*) 'npoin =',npoin
         stop
      endif
      ntail= 0
      ifrag= ntgpt
 250  continue
      if(ifrag.gt.nadrs) then
         write(7777,*) 'Error Message from subroutine search 2 !      '
         write(7777,*) 'seaching has been failed !       '
         write(7777,*) 'ifrag =',ifrag
         write(7777,*) 'nadrs =',nadrs
         stop
      endif
      ntail= ntail+1
      ltail(ntail+1)= isrch(1,ifrag)
      if(isrch(2,ifrag).eq.0) then
         goto 260
      else
         ifrag= isrch(2,ifrag)
      endif
      goto 250
 260  continue
      ltail(1)= ntail
 900  format(/5x,'** Message from subrouthine search **',
     $       /5x,'   seaching has been failed !       ',
     $       /5x,'   ifrag = ',i10,
     $       /5x,'   nadrs = ',i10)
 910  format(/5x,'** Message from subrouthine search **',
     $       /5x,'   number of a target point is out of the range ! ',
     $       /5x,'   ntgpt = ',i10,
     $       /5x,'   npoin = ',i10)
      return
      end
*
**********************************************************************
*     making information table of Gauss points for domain integral 
**********************************************************************
      subroutine mkitgf(mpoin,mef,ng,msrch,mtail,mnod,
     $     npoin,nef,nadrs,x,nodf,isrch,ltail,
     $     ngfcls,nifcls,iflnod,nflnod,aldm,
     $     mgpf,mavnd,ngpf,xxgpf,wwgpf,dmaxf,ndgpf,idgpf)
      implicit real*8(a-h,o-z)
      integer    npoin
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  x(ng,mpoin+3)
      dimension  nodf(3,mef)
      dimension  iflnod(mnod), nflnod(0:10)
      dimension  xxgpf(ng,mgpf), wwgpf(mgpf), dmaxf(mgpf)
      dimension  ndgpf(0:mgpf),idgpf(mgpf*mavnd)
      dimension  ndelm(3), xnode(2,3), xx(2)
      common /gnntri/ nnss(4)
      common /gustri/ sstg(3,7,4),wwtg(7,4)
      nss= nnss(ngfcls)
      igpf= 0
      ndgpf(0)= 0
      iptot= 0
      do 100 ief= 1,nef
         do 110 inode= 1,3
            ndelm(inode)= nodf(inode,ief)
            do 200 ig= 1,ng
               xnode(ig,inode)= x(ig,ndelm(inode))
 200        continue
 110     continue
         wwj= (xnode(1,2)-xnode(1,1))*(xnode(2,3)-xnode(2,1))
     $       -(xnode(2,2)-xnode(2,1))*(xnode(1,3)-xnode(1,1))
         do 120 iss= 1, nss
            igpf= igpf+1
            if(igpf.gt.mgpf) then
               write(7777,*) 'Error Message from
     $              subroutine mkitgf 1 !   '
               write(7777,*) 'igpf is greater than mgpf  '
               write(7777,*) 'mgpf is too small '
               write(7777,*) 'igpf =',igpf
               write(7777,*) 'mgpf =',mgpf
               stop
            endif
            do 210 ig= 1,ng
               xx(ig)= 0.0d0
 210        continue
            do 130 inode= 1,3
               do 220 ig= 1,ng
                  xx(ig)= xx(ig) + 
     $                 sstg(inode,iss,ngfcls)*xnode(ig,inode)
 220           continue
 130        continue
            www= 0.5d0*wwtg(iss,ngfcls)*wwj
         call getndf(mpoin,ng,msrch,mtail,mnod,
     $        npoin,nadrs,isrch,ltail,ndelm,
     $        nifcls,nnod,iflnod,nflnod,x,xx,aldm,ddmax)
            iptot= iptot+nnod
            ndgpf(igpf)= iptot
            if(iptot.gt.mgpf*mavnd) then
               write(7777,*) 
     $              'Error Message from subroutine mkitgf 2 !   '
               write(7777,*) 
     $              'iptot is greater than mgpf*mavnd for idgpf  '
               write(7777,*) 'mgpf or mavnd may be too small  '
               write(7777,*) 'iptot =',iptot
               write(7777,*) 'mgpf*mavnd =',mgpf*mavnd
               stop
            endif
            do 180 i= 1,nnod
               idgpf(ndgpf(igpf-1)+i)= iflnod(i)
 180        continue
            do 240 ig= 1,ng
               xxgpf(ig,igpf)= xx(ig)
 240        continue
            wwgpf(igpf)= www
            dmaxf(igpf)= ddmax
 120     continue
 100  continue
      ngpf= igpf
      return
      end
*
**********************************************************************
*     getting nodes in the influence region
**********************************************************************
      subroutine getndf(mpoin,ng,msrch,mtail,mnod,
     $     npoin,nadrs,isrch,ltail,ndelm,
     $     nifcls,nnod,iflnod,nflnod,x,xx,aldm,ddmax)
      implicit real*8(a-h,o-z)
      integer    npoin
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  x(ng,mpoin+3)
      dimension  xx(2)
      dimension  ndelm(3)
      dimension  iflnod(mnod), nflnod(0:10)
*
*---- search nodes under influence
*
      nflnod(0)= 0
      nflnod(1)= 3
      do 230 inode= 1,3
         iflnod(inode)= ndelm(inode)
 230  continue
      dd1= sqrt((x(1,ndelm(1))-xx(1))**2+
     $     (x(2,ndelm(1))-xx(2))**2 )
      dd2= sqrt((x(1,ndelm(2))-xx(1))**2+
     $     (x(2,ndelm(2))-xx(2))**2 )
      dd3= sqrt((x(1,ndelm(3))-xx(1))**2+
     $     (x(2,ndelm(3))-xx(2))**2 )
      ddmax= aldm*max(dd1,dd2,dd3)
      ip= 3
      if(nifcls.ge.2) then
         jfrag= 0
         iiclas= 1
         do 140 while(iiclas.le.nifcls-1.and.jfrag.eq.0)
            if(nflnod(iiclas).gt.nflnod(iiclas-1)) then
               do 150 i= nflnod(iiclas-1)+1, nflnod(iiclas)
                  ntgpt= iflnod(i)
                  call search(msrch,npoin,ntgpt,nadrs,
     $                 isrch,mtail,ltail)
                  if(ltail(1).ge.1) then
                     do 160 j= 2, ltail(1)+1
                        ifrag= 1
                        do 170 k=1, ip
                           if(iflnod(k).eq.ltail(j)) then
                              ifrag= 0
                           endif
 170                    continue
                        if(ifrag.eq.1) then
                           ip= ip+1
                           if(ip.gt.mnod) then
                              write(7777,*)
     $                             'Error Message from subroutine 
     $                             getndf 1 !'
                              write(7777,*) 
     $                             'ip is greater than mnod for iflnod'
                              write(7777,*) 
     $                             'Too large nifcls or too small mnod'
                              write(7777,*) 'ip =',ip
                              write(7777,*) 'mnod =',mnod
                              stop
                           endif
                           iflnod(ip)= ltail(j)
                        endif
 160                 continue
                  endif
 150           continue
            endif
            nflnod(iiclas+1)= ip
            ip= nflnod(iiclas)
            jfrag= 1
            if(nflnod(iiclas+1).gt.nflnod(iiclas)) then
               do 180 i= nflnod(iiclas)+1, nflnod(iiclas+1)
                  ipoin= iflnod(i)
                  xi1 = x(1,ipoin)
                  xi2 = x(2,ipoin)
                  dd= sqrt( (xx(1)-xi1)**2+
     $                 (xx(2)-xi2)**2  )
                  if(dd.le.ddmax) then
                     ip= ip+1
                     iflnod(ip)= ipoin
                     jfrag= 0
                  endif
 180           continue
            endif
            nflnod(iiclas+1)= ip
            iiclas= iiclas+1
 140     continue
      endif
      nnod= ip
      return
      end
*
**********************************************************************
*     constants for gauss integration formula (triangle area)
**********************************************************************
      block data gaust
      implicit real*8(a-h,o-z)
      common /gnntri/ nnss(4)
      common /gustri/ sstg(3,7,4),wwtg(7,4)
      data nnss/1,3,4,7/
      data (sstg(i,1,1),i=1,3)
     $     /0.3333333333d0,0.3333333333d0,0.3333333333d0/
      data (sstg(i,2,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,3,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,4,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,5,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,6,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,7,1),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,1,2),i=1,3)
     $     /0.5000000000d0,0.5000000000d0,0.0000000000d0/
      data (sstg(i,2,2),i=1,3)
     $     /0.0000000000d0,0.5000000000d0,0.5000000000d0/
      data (sstg(i,3,2),i=1,3)
     $     /0.5000000000d0,0.0000000000d0,0.5000000000d0/
      data (sstg(i,4,2),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,5,2),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,6,2),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,7,2),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,1,3),i=1,3)
     $     /0.3333333333d0,0.3333333333d0,0.3333333333d0/
      data (sstg(i,2,3),i=1,3)
     $     /0.8000000000d0,0.2000000000d0,0.2000000000d0/