search.f

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

      data (sstg(i,3,3),i=1,3)
     $     /0.2000000000d0,0.8000000000d0,0.2000000000d0/
      data (sstg(i,4,3),i=1,3)
     $     /0.2000000000d0,0.2000000000d0,0.8000000000d0/
      data (sstg(i,5,3),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,6,3),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,7,3),i=1,3)
     $     /0.0000000000d0,0.0000000000d0,0.0000000000d0/
      data (sstg(i,1,4),i=1,3)
     $     /0.3333333333d0,0.3333333333d0,0.3333333333d0/
      data (sstg(i,2,4),i=1,3)
     $     /0.7974269853d0,0.1012865073d0,0.1012865073d0/
      data (sstg(i,3,4),i=1,3)
     $     /0.1012865073d0,0.7974269853d0,0.1012865073d0/
      data (sstg(i,4,4),i=1,3)
     $     /0.1012865073d0,0.1012865073d0,0.7974269853d0/
      data (sstg(i,5,4),i=1,3)
     $     /0.0597158717d0,0.4701420641d0,0.4701420641d0/
      data (sstg(i,6,4),i=1,3)
     $     /0.4701420641d0,0.0597158717d0,0.4701420641d0/
      data (sstg(i,7,4),i=1,3)
     $     /0.4701420641d0,0.4701420641d0,0.0597158717d0/
      data (wwtg(i,1),i=1,7)
     $     / 1.0000000000d0, 0.0000000000d0, 0.0000000000d0,
     $       0.0000000000d0, 0.0000000000d0, 0.0000000000d0,
     $       0.0000000000d0/
      data (wwtg(i,2),i=1,7)
     $     / 0.3333333333d0, 0.3333333333d0, 0.3333333333d0,
     $       0.0000000000d0, 0.0000000000d0, 0.0000000000d0,
     $       0.0000000000d0/
      data (wwtg(i,3),i=1,7)
     $     /-0.5625000000d0, 0.5208333333d0, 0.5208333333d0,
     $       0.5208333333d0, 0.0000000000d0, 0.0000000000d0,
     $       0.0000000000d0/
      data (wwtg(i,4),i=1,7)
     $     / 0.2250000000d0, 0.1259391805d0, 0.1259391805d0,
     $       0.1259391805d0, 0.1323941527d0, 0.1323941527d0,
     $       0.1323941527d0/
      end
*
**********************************************************************
*     making information table of Gauss points for boundary integral 
**********************************************************************
      subroutine mkitgb(mpoin,meb,ng,msrch,mtail,mnod,
     $     npoin,nadrs,x,nodb,isrch,ltail,
     $     ngbcls,nibcls,iflnod,nflnod,aldm,
     $     nbc,nbcbm,mgpb,mavnd,ngpb,
     $     xxgpb,wwgpb,dmaxb,ndgpb,idgpb)
      implicit real*8(a-h,o-z)
      integer    npoin
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  x(ng,mpoin+3)
      dimension  nodb(2,meb)
      dimension  iflnod(mnod), nflnod(0:10)
      dimension  xxgpb(ng,mgpb), wwgpb(mgpb), dmaxb(mgpb)
      dimension  ndgpb(0:mgpb),idgpb(mgpb*mavnd)
      dimension  nbcbm(3,meb)
      dimension  ndelm(2), xnode(2,2)
      dimension  xx(2)
      common/gaussa/ssg(10,10),wwg(10,10)
      nss= ngbcls
      igpb= 0
      ndgpb(0)= 0
      iptot= 0
      do 100 ibc= 1,nbc
         ieb= nbcbm(1,ibc)
         do 200 inode= 1,2
            ndelm(inode)= nodb(inode,ieb)
            do 210 ig= 1,ng
               xnode(ig,inode)= x(ig,ndelm(inode))
 210        continue
 200     continue
         aleng= sqrt((xnode(1,2)-xnode(1,1))*(xnode(1,2)-xnode(1,1))
     $              +(xnode(2,2)-xnode(2,1))*(xnode(2,2)-xnode(2,1)))
         wwj= 0.5d0*aleng
         do 120 iss= 1, nss
            igpb= igpb+1
            if(igpb.gt.mgpb) then
               write(7777,*) 
     $              'Error Message from subroutine mkitgb 1 !   '
               write(7777,*) 'igpb is greater than mgpb  '
               write(7777,*) 'mgpb is too small '
               write(7777,*) 'igpb =',igpb
               write(7777,*) 'mgpb =',mgpb
               stop
            endif
            do 350 ig= 1,ng
               xx(ig)= 0.5d0*(1.0d0-ssg(iss,ngbcls))*xnode(ig,1)
     $               + 0.5d0*(1.0d0+ssg(iss,ngbcls))*xnode(ig,2)
 350        continue
            www= wwj*wwg(iss,ngbcls)
*
*----    search nodes under influence
*
            call getndb(mpoin,ng,msrch,mtail,mnod,
     $           npoin,nadrs,isrch,ltail,ndelm,
     $           nibcls,nnod,iflnod,nflnod,x,xx,aldm,ddmax)
*
*---- store coordinates, influence nodes, wait of a gauss point.
*
            iptot= iptot+nnod
            ndgpb(igpb)= iptot
            if(iptot.gt.mgpb*mavnd) then
               write(7777,*) 
     $              'Error Message from subroutine mkitgb 2 !   '
               write(7777,*) 
     $              'iptot is greater than mgpb*mavnd for idgpb  '
               write(7777,*) 'mgpb or mavnd may be too small  '
               write(7777,*) 'iptot =',iptot
               write(7777,*) 'mgpb*mavnd =',mgpb*mavnd
               stop
            endif
            do 180 i= 1,nnod
               idgpb(ndgpb(igpb-1)+i)= iflnod(i)
 180        continue
             do 240 ig= 1,ng
               xxgpb(ig,igpb)= xx(ig)
 240        continue
            wwgpb(igpb)= www
            dmaxb(igpb)= ddmax
 120     continue
 100  continue
      ngpb= igpb
      return
      end
*
**********************************************************************
*     getting nodes in the influence region
**********************************************************************
      subroutine getndb(mpoin,ng,msrch,mtail,mnod,
     $              npoin,nadrs,isrch,ltail,ndelm,
     $              nibcls,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(2)
      dimension  iflnod(mnod), nflnod(0:10)
*
*---- search nodes under influence
*
      nflnod(0)= 0
      nflnod(1)= 2
*
*---- set ddmax 
*
      d = 0.0d0
      dm = 1.0d30
      do i = 1,2
         ntgpt = ndelm(i)
         call search(msrch,npoin,ntgpt,nadrs,
     $        isrch,mtail,ltail)
         do j=2,ltail(1)+1
            if(ltail(j).ne.ndelm(1).and.ltail(j).ne.ndelm(2)) then
               d = sqrt((x(1,ltail(j))-xx(1))**2+
     $              (x(2,ltail(j))-xx(2))**2 )
               if(d.lt.dm) then
                  ipp = ltail(j)
                  dm = d
               endif
            endif
         enddo
      enddo
      iflnod(1)= ndelm(1)
      iflnod(2)= ndelm(2)
      iflnod(3)= ipp
      dd12= sqrt((x(1,ndelm(1))-xx(1))**2+
     $     (x(2,ndelm(1))-xx(2))**2 )
      dd23= sqrt((x(1,ndelm(2))-xx(1))**2+
     $     (x(2,ndelm(2))-xx(2))**2 )
      ddmax= aldm*max(dd12,dd23,dm)
      ip= 3
      if(nibcls.ge.2) then
         jfrag= 0
         iiclas= 1
         do 140 while(iiclas.le.nibcls-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 
     $                             getndb 1 !'
                              write(7777,*) 
     $                             'ip is greater than mnod for iflnod'
                              write(7777,*) 
     $                             'Too large nibcls 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 (line integral)
**********************************************************************
      block data gauss
      implicit real*8(a-h,o-z)
      common/gaussa/ssg(10,10),wwg(10,10)
      data ssg(1,1) / 0.0d0/
      data (ssg(i,2),i=1,2)/
     $ -0.577350269189626d0, 0.577350269189626d0/
      data (ssg(i,3),i=1,3)/
     $ -0.774596669241483d0, 0.000000000000000d0, 0.774596669241483d0/
      data (ssg(i,4),i=1,4)/
     $ -0.861136311594953d0,-0.339981043584856d0, 0.339981043584856d0,
     $  0.861136311594953d0/
      data (ssg(i,5),i=1,5)/
     $ -0.906179845938664d0,-0.538469310105683d0, 0.000000000000000d0,
     $  0.538469310105683d0, 0.906179845938664d0/
      data (ssg(i,6),i=1,6)/
     $ -0.932469514203152d0,-0.661209386466265d0,-0.238619186083197d0,
     $  0.238619186083197d0, 0.661209386466265d0, 0.932469514203152d0/
      data (ssg(i,7),i=1,7)/
     $ -0.949107912342759d0,-0.741531185599394d0,-0.405845151377397d0,
     $  0.000000000000000d0, 0.405845151377397d0, 0.741531185599394d0,
     $  0.949107912342759d0/
      data (ssg(i,8),i=1,8)/
     $ -0.960289856497536d0,-0.796666477413627d0,-0.525532409916329d0,
     $ -0.183434642495650d0, 0.183434642495650d0, 0.525532409916329d0,
     $  0.796666477413627d0, 0.960289856497536d0/
      data (ssg(i,9),i=1,9)/
     $ -0.968160239507626d0,-0.836031107326636d0,-0.613371432700590d0,
     $ -0.324253423403809d0, 0.000000000000000d0, 0.324253423403809d0,
     $  0.613371432700590d0, 0.836031107326636d0, 0.968160239507626d0/
      data (ssg(i,10),i=1,10)/
     $ -0.973906528517172d0,-0.865063366688985d0,-0.679409568299024d0,
     $ -0.433395394129247d0,-0.148874338981631d0, 0.148874338981631d0,
     $  0.433395394129247d0, 0.679409568299024d0, 0.865063366688985d0,
     $  0.973906528517172d0/
      data wwg(1,1)/ 2.0d0/
      data (wwg(i,2),i=1,2)/
     $  1.000000000000000d0, 1.000000000000000d0/
      data (wwg(i,3),i=1,3)/
     $  0.555555555555556d0, 0.888888888888889d0, 0.555555555555556d0/
      data (wwg(i,4),i=1,4)/
     $  0.347854845137454d0, 0.652145154862546d0, 0.652145154862546d0,
     $  0.347854845137454d0/
      data (wwg(i,5),i=1,5)/
     $  0.236926885056189d0, 0.478628670499366d0, 0.568888888888889d0,
     $  0.478628670499366d0, 0.236926885056189d0/
      data (wwg(i,6),i=1,6)/
     $  0.171324492379170d0, 0.360761573048139d0, 0.467913934572691d0,
     $  0.467913934572691d0, 0.360761573048139d0, 0.171324492379170d0/
      data (wwg(i,7),i=1,7)/
     $  0.129484966168870d0, 0.279705391489277d0, 0.381830050505119d0,
     $  0.417959183673469d0, 0.381830050505119d0, 0.279705391489277d0,
     $  0.129484966168870d0/
      data (wwg(i,8),i=1,8)/
     $  0.101228536290376d0, 0.222381034453374d0, 0.313706645877887d0,
     $  0.362683783378362d0, 0.362683783378362d0, 0.313706645877887d0,
     $  0.222381034453374d0, 0.101228536290376d0/
      data (wwg(i,9),i=1,9)/
     $  0.081274388361574d0, 0.180648160694857d0, 0.260610696402935d0,
     $  0.312347077040003d0, 0.330239355001260d0, 0.312347077040003d0,
     $  0.260610696402935d0, 0.180648160694857d0, 0.081274388361574d0/
      data (wwg(i,10),i=1,10)/
     $  0.066671344308688d0, 0.149451349150581d0, 0.219086362515982d0,
     $  0.269266719309996d0, 0.295524224714753d0, 0.295524224714753d0,
     $  0.269266719309996d0, 0.219086362515982d0, 0.149451349150581d0,
     $  0.066671344308688d0/
      end