tstiff.f

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

            call invar(ncrit,props,stres,
     $           devia,sint3,steff,theta,varj2,yield)
            call yieldf(ncrit,props,nstr1,
     $           devia,steff,theta,varj2,avect)
            call flowpl(ncrit,props,ntype,hdpara,nstr1,
     $           gepst,avect,abeta,dvect)
            do istre= 1,nstre
               do jstre= 1,nstre
                  dmat(istre,jstre)= dmat(istre,jstre)
     $                 -abeta*dvect(istre)*dvect(jstre)
               enddo
            enddo
         endif
      do 200 i= 1,nnod*ng
         do 210 j= 1,nnod*ng
            ekmat(i,j)= 0.0d0
 210     continue
 200  continue
      do 100 i= 1,nnod*ng
         do 110 j= 1,nnod*ng
            do 120 k= 1,nstrs
               do 130 l= 1,nstrs
                  ekmat(i,j)= ekmat(i,j) + 
     $                 bmat(k,i)*dmat(k,l)*bmat(l,j)
 130           continue
 120        continue
 110     continue
 100  continue
      do 300 i= 1,nnod*ng
         do 310 j= 1,nnod*ng
            ekmat(i,j)= ekmat(i,j)*www
 310     continue
 300  continue
      return
      end
*
**********************************************************************
*     Making Boundary Elemental Stiffness Matrix
**********************************************************************
      subroutine bstiff(ng,mnod,nnod,ntype,props,dm,cm,kk,
     $     uu,xi,www,aa,ai,dadx,bl,dldx,phi,
     $     nstre,nstr1,gepst,ncrit,hdpara,stres,iincs,
     $     smat,ekmat,eforc,geppre)
      implicit real*8(a-h,o-z)
      dimension  props(7)
      dimension  uu(2)
      dimension  xi(ng,mnod)
      dimension  aa(6,6),ai(6,6)
      dimension  dadx(6,6,2)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  phi(2,mnod*ng)
      dimension  smat(2,2)
      dimension  ekmat(mnod*ng,mnod*ng), eforc(mnod*ng)
      dimension  dmat(4,4)
      dimension  hdpara(20),ncrit(3)
      dimension  devia(4),stres(4),dvect(4),avect(4)
      md = 3
      call makeaa(ng,mnod,md,nnod,xi,dm,cm,
     $     kk,aa,bl,ai,dadx,dldx)
      call itfunc(ng,mnod,md,nnod,bl,ai,phi)
      call dmatrx(ntype,props,nstrs,dmat)
         if(iincs.ne.1.and.gepst.ne.0.0d0
     $     .and.geppre.ne.0.0d0) then
            call invar(ncrit,props,stres,
     $           devia,sint3,steff,theta,varj2,yield)
            call yieldf(ncrit,props,nstr1,
     $           devia,steff,theta,varj2,avect)
            call flowpl(ncrit,props,ntype,hdpara,nstr1,gepst,
     $           avect,abeta,dvect)
            do istre= 1,nstre
               do jstre= 1,nstre
                  dmat(istre,jstre)= dmat(istre,jstre)
     $                 -abeta*dvect(istre)*dvect(jstre)
               enddo
            enddo
         endif
      do 200 i= 1,nnod*ng
         do 210 j= 1,nnod*ng
            ekmat(i,j)= 0.0d0
 210     continue
         eforc(i)= 0.0d0
 200  continue
      do 100 i= 1,2
         do 110 j= 1,nnod*ng
            do 120 k= 1,2
               do 130 l= 1,nnod*ng
                     ekmat(j,l)= ekmat(j,l) +
     $                    phi(i,j)*smat(i,k)*phi(k,l)
 130           continue
 120        continue
 110     continue
 100  continue
      do 400 i= 1,nnod*ng
         do 410 j= 1,nnod*ng
            ekmat(i,j)= ekmat(i,j)*www
 410     continue
 400  continue
      do 500 i= 1,2
         do 510 j= 1,nnod*ng
            do 520 k= 1,2
               eforc(j)= eforc(j) + phi(i,j)*smat(i,k)*uu(k)
 520        continue
 510     continue
 500  continue
      do 580 i= 1,nnod*ng
            eforc(i)= eforc(i)*www
 580  continue
      return
      end
*
**********************************************************************
*     This subroutine evaluates the consistent nodal forces
**********************************************************************
      subroutine loadps(mpoin,meb,ng,mnod,x,nodb,ngbcls,iflnod,
     $     alcm,kk,nibcls,xi,bl,dldx,phi,eforc,rload,nbc,nbcbm,
     $     dbcbm,mgpb,mavnd,dmaxb,ndgpb,idgpb,kfx)
      implicit real*8(a-h,o-z)
      dimension  x(ng,mpoin)
      dimension  nodb(2,meb)
      dimension  iflnod(mnod)
      dimension  nbcbm(3,meb), dbcbm(4,meb)
      dimension  rload(mpoin*ng)
      dimension  xi(ng,mnod)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  phi(2,mnod*ng)
      dimension  eforc(mnod*ng)
      dimension  dmaxb(mgpb)
      dimension  ndgpb(0:mgpb),idgpb(mgpb*mavnd)
      dimension  ndelm(2), xnode(2,2), tnode(2,2)
      dimension  xx(2), tt(2)
      dimension  aa(6,6),ai(6,6)
      dimension  dadx(6,6,2)
      dimension   kfx(mpoin)
      common/gaussa/ssg(10,10),wwg(10,10)
      nss= ngbcls
      igpb= 0
      do 100 ibc= 1,nbc
         ieb= nbcbm(1,ibc)
         do inode= 1,2
            ndelm(inode)= nodb(inode,ieb)
            do ig= 1,ng
               xnode(ig,inode)= x(ig,ndelm(inode))
            enddo
         enddo
         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 ig= 1,ng
            if(nbcbm(1+ig,ibc).eq.0) then
               tnode(ig,1)= 0.0d0
               tnode(ig,2)= 0.0d0
            else
               tnode(ig,1)= dbcbm(ig*2-1,ibc)
               tnode(ig,2)= dbcbm(ig*2,ibc)
            endif
         enddo
         do 120 iss= 1, nss
            igpb= igpb + 1
            do 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)
               tt(ig)= 0.5d0*(1.0d0-ssg(iss,ngbcls))*tnode(ig,1)
     $              + 0.5d0*(1.0d0+ssg(iss,ngbcls))*tnode(ig,2)
            enddo
            www= wwj*wwg(iss,ngbcls)
            ddmax= dmaxb(igpb)
            nnod= ndgpb(igpb) - ndgpb(igpb-1)
            do inod= 1, nnod
               iflnod(inod)= idgpb(ndgpb(igpb-1)+inod)
            enddo
            do inod= 1, nnod
               do ig= 1,2
                  xi(ig,inod)= x(ig,iflnod(inod)) - xx(ig)
               enddo
            enddo
            dm= ddmax
            cm= alcm*dm
            if(dm.lt.1.0d-30) then
               write(7777,*)
     $              'Error Message from subroutine loadps 1 !    '
               write(7777,*)
     $              '     Diameter of infulence region is zero !'
               stop
            endif
            md = 3
            call makeaa(ng,mnod,md,nnod,xi,dm,cm,
     $           kk,aa,bl,ai,dadx,dldx)
            call itfunc(ng,mnod,md,nnod,bl,ai,phi)
            do i= 1,nnod*ng
               eforc(i)= 0.0d0
            enddo
            do i= 1,nnod*ng
               do j= 1,ng
                  eforc(i)= eforc(i)+www*phi(j,i)*tt(j)
               enddo
            enddo
            do i = 1,nnod
               iflnod(i) = kfx(iflnod(i))
            enddo
            call formtf(ng,mpoin,mnod,nnod,iflnod,
     $           eforc,rload)
 120     continue
 100  continue
      return
      end
*
**********************************************************************
*     Forming total stiffness matrix
**********************************************************************
      subroutine formts(msizl,mpoin,ng,stotl,index,idofh,idiag,
     $     mnod,nnod,iflnod,ekmat)
      implicit real*8(a-h,o-z)
      dimension ekmat(mnod*ng,mnod*ng)
      dimension stotl(msizl)
      dimension idofh(mpoin*ng),idiag(0:mpoin*ng)
      dimension iflnod(mnod),index(mnod*ng)
      call indxel(ng,mnod,nnod,iflnod,index)
      ncn = nnod*ng
      do 350 ii= 1,ncn
         i= index(ii)
         do 330 jj= 1,ncn
            j= index(jj)
            if(i.ge.j) then
               ill= idiag(i)-i+j
               illmax= idiag(i)
               illmin= idiag(i)-idofh(i)+1
               if(ill.lt.illmin.or.ill.gt.illmax) then
                  write(7777,*) 'Error Message from
     $                 subroutine formts 1 !'
                  write(6,9100) ii,jj,i,j,illmax,illmin,ill,idiag(i)
                  stop
               endif
               stotl(ill)= stotl(ill)+ekmat(ii,jj)
            endif
 330     continue
 350  continue
      return
 9100 format(/2x,'irregular skyline storage has been tried !!!'/
     $     2x,'(ii,jj)=        ',2i5/
     $     2x,'( i, j)=        ',2i5/
     $     2x,'(illmax,illmin)=  ',2i10/
     $     2x,'ill=            ',i10/
     $     2x,'idiag(i)=       ',i10/)
      end
*
**********************************************************************
*     Forming total forc vector
*       RCM version
**********************************************************************
      subroutine formtf(ng,mpoin,mnod,nnod,iflnod,eforc,tforc)
      implicit real*8(a-h,o-z)
      dimension  tforc(mpoin*ng)
      dimension  iflnod(mnod)
      dimension  eforc(mnod*ng)
      do 300 inod= 1,nnod
         ii= iflnod(inod)
         do 200 ig= 1,ng
            tforc((ii-1)*ng+ig) =
     $           tforc((ii-1)*ng+ig) + eforc((inod-1)*ng+ig)
 200     continue
 300  continue
      return
      end