skysub.f

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

********************************************************************** 
*     making skyline index                                           *
**********************************************************************
      subroutine skysub(mpoin,mnod,mgpf,mgpb,mavnd,msizl,ng,nszf,idofh,
     $     iflnod,ngpf,ngpb,ndgpf,idgpf,ndgpb,idgpb,index,lhigh,
     $     idiag,kfx)
      implicit real*8(a-h,o-z)
      integer idofh(mpoin*ng),index(mnod*ng),lhigh(mnod*ng)
      dimension  iflnod(mnod),ndgpf(0:mgpf),idgpf(mgpf*mavnd)
      dimension  kfx(mpoin),idiag(0:nszf)
      dimension  ndgpb(0:mgpb),idgpb(mgpb*mavnd)
      call zeroi(nszf,idofh)
*---- for domain
      call skyhi(mpoin,mnod,mgpf,mavnd,ng,ngpf,
     $     idofh,iflnod,ndgpf,idgpf,index,lhigh,kfx)
*---- for boundary
      call skyhi(mpoin,mnod,mgpb,mavnd,ng,ngpb,
     $     idofh,iflnod,ndgpb,idgpb,index,lhigh,kfx)
      call skydia(mpoin,ng,msizl,nszf,idofh,idiag)
      return
      end
*
********************************************************************** 
*     find column hights of system equations                         *
*     skyline store mode                                             *
**********************************************************************
      subroutine skyhi(mpoin,mnod,mgpf,mavnd,ng,ngpf,
     $     idofh,iflnod,ndgpf,idgpf,index,lhigh,kfx)
      implicit real*8(a-h,o-z)
      integer idofh(mpoin*ng),index(mnod*ng),lhigh(mnod*ng)
      dimension  iflnod(mnod),ndgpf(0:mgpf),idgpf(mgpf*mavnd)
      dimension  kfx(mpoin)
*
*     loop over gaussian point
*
      do 20 igpf=1,ngpf
         nnod= ndgpf(igpf)-ndgpf(igpf-1)
         nelfre= nnod*ng
         do 220 inod= 1, nnod
            iflnod(inod)= kfx(idgpf(ndgpf(igpf-1)+inod))
 220     continue
         call indxel(ng,mnod,nnod,iflnod,index)
*
*     find element column heights
*
         call elhigh(mnod,ng,nelfre,index,lhigh)
*
*     compare with current maximums
*
         do 10 j=1,nelfre
            ndx= index(j)
            if(ndx.ge.1) then
               if(idofh(ndx).lt.lhigh(j)) then
                  idofh(ndx)= lhigh(j)
               endif
            endif
 10      continue
 20   continue
      return
      end
**********************************************************************
*     use column height to find diagonal coefficients for            *
*     symmetric skyline storage mode                                 *
**********************************************************************
      subroutine skydia(mpoin,ng,msizl,nszf,idofh,idiag)
      implicit real*8(a-h,o-z)
	integer idofh(mpoin*ng),idiag(0:nszf)
      ipoin= 0
      idiag(0)= 0
      do 10 i= 1,nszf
         ipoin= ipoin + idofh(i)
         idiag(i)= ipoin
 10   continue
      ifullm= (nszf*nszf-nszf)/2 + nszf
      write(6,100) ifullm, ipoin, ipoin, msizl
      if(ipoin.gt.msizl) then
         write(7777,*) 'Error Message from subroutine skydia 1 !      '
         write(7777,*) 'Maximum size of skyline storage is too small !'
         write(7777,*) 'Please change the parameter msizl larger than',
     $        ipoin
         write(7777,*) 'in main program and re-compile !!'
         stop
      endif
      if(ipoin*2-nszf.gt.nszf*nszf) then
         write(7777,*) 'Error Message from subroutine skydia 2 !      '
         write(7777,*) 'skyline storage size of matrix is larger !'
         stop
      endif
 100  format(//5x,75('*'),
     $        /8x,'Message from Skydia',
     $        /8x,'Half matrix size of whole stiffness',i10,
     $        /8x,'Skyline storage size of matrix     ',i10,
     $       //8x,'Parameter msizl must be more than  ',i10,
     $        /8x,'Current size of parameter msizl    ',i10,
     $        /5x,75('*'))
      return
      end
*
********************************************************************** 
*     determine degrees of freedom numbers of infuluence region      *
**********************************************************************
      subroutine indxel(ng,mnod,nnod,iflnod,index)
      implicit real*8(a-h,o-z)
	integer iflnod(mnod),index(mnod*ng)
      do 20 k= 1,nnod
         index(k*2-1) = (iflnod(k)-1)*ng+1
         index(k*2) = (iflnod(k)-1)*ng+2        
 20   continue
      return
      end
*
********************************************************************** 
*     find system column heights of an infuluence region             *
**********************************************************************
      subroutine elhigh(mnod,ng,nfree,index,lhigh)
      implicit real*8(a-h,o-z)
      integer index(mnod*ng),lhigh(mnod*ng)
      min= index(1)
      do 10 i= 2,nfree
         lhigh(i)= 0
         ndx= index(i)
         if(ndx.ge.1) then
            if(ndx.lt.min) then
               min=ndx
            endif
         endif
 10   continue
      min= min-1
      do 20 i= 1,nfree
         ndx= index(i)
         if(ndx.ge.1) then
            lhigh(i)= ndx - min
         endif
 20   continue
      return
      end