plast.f

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

**********************************************************************
*     This subroutine checks for convergence of the iteration process
**********************************************************************
      subroutine conver(nszf,asdis,bsdis,eload,tload,
     $     iiter,nchek,toler,pvalu)
      implicit real*8(a-h,o-z)
      dimension  asdis(nszf),bsdis(nszf)
      dimension  eload(nszf),tload(nszf)
      resid= 0.0d0
      retot= 0.0d0
      remax= 0.0d0
      nchek= 0
      nswch= 0
*
*---- Check boundary condition.
*
      do iszf= 1,nszf
         if(tload(iszf).ne.0.0d0) then
            nswch= 1
            goto 100
         endif
      enddo
 100  continue
      if(nswch.eq.1) then
         do iszf= 1,nszf
            if(tload(iszf).ne.0.0d0) then
               refor= tload(iszf) - eload(iszf)
               resid= resid + refor*refor
               retot= retot + tload(iszf)*tload(iszf)
               agash= abs(refor)
               if(agash.gt.remax) then
                  remax= agash
               endif
               eload(iszf)= refor
            else
               eload(iszf)= 0.0d0
            endif
         enddo
      else
         do iszf=1,nszf
            resid = resid + asdis(iszf) * asdis(iszf)
            retot = retot + bsdis(iszf) * bsdis(iszf)
            eload(iszf)= - eload(iszf)
         enddo
      endif
      resid= sqrt(resid)
      retot= sqrt(retot)
      if(retot.ne.0.0d0) then
         ratio= 100.0d0*resid/retot
      else
         write(7777,*) 'Error Message from subroutine conver 1 !      '
         write(7777,*) 'retot is equal to zero !'
         stop
      endif
      if(ratio.gt.toler) then
         nchek= 1
      endif
      if(iiter.gt.1.and.ratio.gt.pvalu) then
         nchek= 999
      endif
      pvalu= ratio
      write(6,30) iiter,nchek,ratio,remax
 30   format(5x,'It=',i4,3x,'Code=',i4,3x,'Sum R-St=',e13.5,3x,
     $     'Max R-St=',e13.5)
      return
      end
*
*********************************************************************
*     Incremnents the applied load
*********************************************************************
      subroutine increm(iincs,nbc,nszf,meb,mpoin,ng,dbcbm,epbc,
     $     rload,tload,eload,facto,tfact,toler,miter,noutp)
      implicit real*8(a-h,o-z)
      dimension  dbcbm(4,meb)
      dimension  epbc(4,meb)
      dimension  rload(mpoin*ng),tload(mpoin*ng),eload(mpoin*ng)
      write(6,900) iincs
 900  format(//5x,75('-')/,8x,'Increment Number= ',i5)
*
*-----reading incremental informations 
*
      read(1500,*) idmy,facto,toler,miter,noutp
*
*-----writing incremental informations 
*
      tfact= tfact+facto
      write(6,960) tfact,toler,miter,noutp
 960  format(8x,'Accumulated Increment Factor = ',f10.5/
     $       8x,'Convergence Tolerance        = ',f10.5/
     $       8x,'Maximum Number of Iterations = ',i5/
     $       8x,'Print No. = ',i5)
      do iszf= 1,nszf
         eload(iszf)= eload(iszf)+rload(iszf)*facto
         tload(iszf)= tload(iszf)+rload(iszf)*facto
      enddo
      do ibc= 1,nbc
         do i = 1,4
            epbc(i,ibc) = dbcbm(i,ibc)*facto
         enddo
      enddo
      return
      end
*
***********************************************************************
*     Regurate algorism of making stiffness matrix and solving matrix 
***********************************************************************
      subroutine algor(iincs,iiter,nalgo,kresl)
      implicit real*8(a-h,o-z)
      kresl= 2
      if(nalgo.eq.1.and.iincs.eq.1.and.iiter.eq.1) kresl= 1
      if(nalgo.eq.2) kresl= 1
      if(nalgo.eq.3.and.iiter.eq.1) kresl= 1
      if(nalgo.eq.4.and.iincs.eq.1.and.iiter.eq.1) kresl= 1
      if(nalgo.eq.4.and.iiter.eq.2) kresl= 1
      if(nalgo.eq.5.and.iiter.eq.1) kresl= 1
      if(nalgo.eq.5.and.iiter.eq.2) kresl= 1
      return
      end