tstiff.f

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

**********************************************************************
*     Domain Integral for making Total Stiffness
*       RCM version
**********************************************************************
      subroutine domint(mpoin,ng,mnod,npoin,x,iflnod,
     $     ntype,props,alcm,kk,
     $     xi,bl,dldx,phi,bmat,ekmat,amxdia,
     $     mgpf,mavnd,ngpf,xxgpf,wwgpf,dmaxf,ndgpf,idgpf,
     $     msizl,mindx,stotl,index,idofh,idiag,kfx,
     $     ncrit,hdpara,epstn,strsg,iincs,gepst,epstnfpre)
      implicit real*8(a-h,o-z)
      integer    npoin
      dimension  x(ng,mpoin+3)
      dimension  iflnod(mnod)
      dimension  props(7)
      dimension  xi(ng,mnod)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  bmat(3,mnod*ng), ekmat(mnod*ng,mnod*ng)
      dimension  phi(2,mnod*ng)
      dimension  xx(2)
      dimension  aa(6,6),ai(6,6)
      dimension  dadx(6,6,2)
      dimension  xxgpf(ng,mgpf), wwgpf(mgpf), dmaxf(mgpf)
      dimension  ndgpf(0:mgpf),idgpf(mgpf*mavnd)
*
*---- for skyline method 
*
      dimension  stotl(msizl)
      integer    idofh(mpoin*ng),idiag(0:mpoin*ng)
      integer    index(mindx)
*
*---- for plastic-elastic problem
*
      dimension  kfx(mpoin)
      dimension  hdpara(20),ncrit(3)
      dimension  epstn(mgpf),strsg(4,mgpf),stres(4)
      dimension  epstnfpre(mgpf)
*
*-----plane strain and plane stress case
*
      nstre = 3
      nstr1 = nstre+1
      do 100 igpf= 1,ngpf
         do 210 ig= 1,ng
            xx(ig)= xxgpf(ig,igpf)
 210     continue
         www= wwgpf(igpf)
         ddmax= dmaxf(igpf)
         nnod= ndgpf(igpf)-ndgpf(igpf-1)
         do 220 inod= 1, nnod
            iflnod(inod)= idgpf(ndgpf(igpf-1)+inod)
 220     continue
         do 230 inod= 1,nnod
            ipoin= iflnod(inod)
            do 240 ig= 1,ng
               xi(ig,inod)= x(ig,ipoin) - xx(ig)
 240        continue
 230     continue
         dm= ddmax
         cm= alcm*dm
         if(dm.lt.1.0d-30) then
            write(7777,*) 'Error Message from subroutine domint 1 !    '
            write(7777,*) '     Diameter of infulence region is zero !'
            stop
         endif
         gepst= epstn(igpf)
         geppre= epstn(igpf)-epstnfpre(igpf)
         do  istr1= 1,nstr1
            stres(istr1)= strsg(istr1,igpf)
         enddo
         call estiff(ng,mnod,nnod,ntype,props,dm,cm,kk,
     $     nstre,nstr1,gepst,ncrit,hdpara,stres,iincs,
     $     xi,www,aa,ai,dadx,bl,dldx,bmat,ekmat,phi,geppre)
*
*---- store elemental stiffness to skyline storage matrix of 
*---- total system
*
         do 300 i = 1,nnod
            iflnod(i) = kfx(iflnod(i))
 300     continue
         call formts(msizl,mpoin,ng,stotl,index,idofh,idiag,
     $              mnod,nnod,iflnod,ekmat)
 100  continue
      amxdia= 0.0d0
      do 250 i= 1,idiag(npoin*ng)
         dia= abs(stotl(i))
         if(dia.gt.amxdia) then
            amxdia= dia
         endif
 250  continue
*
*---- writing penalty coefficient
*
      return
      end
*
**********************************************************************
*     boundary integral
*       RCM version
**********************************************************************
      subroutine bunint(mpoin,meb,ng,mnod,x,nodb,
     $     ngbcls,iflnod,ntype,props,alcm,kk,nibcls,
     $     xi,bl,dldx,phi,ekmat,eforc,
     $     tforc,nbc,nbcbm,dbcbm,pnlty,
     $     mgpb,mavnd,dmaxb,ndgpb,idgpb,
     $     msizl,mindx,stotl,index,idofh,idiag,kfx,
     $     ncrit,hdpara,epstnb,strsgb,iincs,gepst,
     $     epstnbpre)
      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  props(7)
      dimension  tforc(mpoin*ng)
      dimension  xi(ng,mnod)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  phi(2,mnod*ng)
      dimension  ekmat(mnod*ng,mnod*ng),eforc(mnod*ng)
      dimension  dmaxb(mgpb)
      dimension  ndgpb(0:mgpb),idgpb(mgpb*mavnd)
      dimension  ndelm(2), xnode(2,2), unode(2,2)
      dimension  xx(2), uu(2)
      dimension  aa(6,6),ai(6,6)
      dimension  dadx(6,6,2)
      dimension  smat(2,2)
      dimension  stotl(msizl)
      integer idofh(mpoin*ng),idiag(0:mpoin*ng)
      integer index(mindx)
      dimension kfx(mpoin)
*
*---- for plastic-elastic problem
*
      dimension  hdpara(20),ncrit(3)
      dimension  epstnb(mgpb),strsgb(4,mgpb),stres(4)
      dimension  epstnbpre(mgpb)
      common/gaussa/ssg(10,10),wwg(10,10)
      nstre = 3
      nstr1 = nstre+1
      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
*     
*---- set unode and tnode
*     
         do ig= 1,ng
            if(nbcbm(1+ig,ibc).eq.0) then
               unode(ig,1)= dbcbm(ig*2-1,ibc)
               unode(ig,2)= dbcbm(ig*2,ibc)
            else
               unode(ig,1)= 0.0d0
               unode(ig,2)= 0.0d0
            endif
         enddo
*     
*---- set smat
*     
         do i= 1,2
            do j= 1,2
               smat(i,j)= 0.0d0
            enddo
         enddo
         if(nbcbm(2,ibc).eq.0) then
            smat(1,1)= 1.0d0*pnlty
         endif
         if(nbcbm(3,ibc).eq.0) then
            smat(2,2)= 1.0d0*pnlty
         endif
         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)
               uu(ig)= 0.5d0*(1.0d0-ssg(iss,ngbcls))*unode(ig,1)
     $              + 0.5d0*(1.0d0+ssg(iss,ngbcls))*unode(ig,2)
            enddo
            www= wwj*wwg(iss,ngbcls)
*     
*---- set local group
*     
            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
*     
*---- make ekmat
*     
            dm= ddmax
            cm= alcm*dm
            if(dm.lt.1.0d-30) then
               write(7777,*)
     $              'Error Message from subroutine bunint 1 !    '
               write(7777,*)
     $              '     Diameter of infulence region is zero !'
               stop
            endif
            gepst= epstnb(igpb)
            geppre= epstnb(igpb) - epstnbpre(igpb)
            do  istr1= 1,nstr1
               stres(istr1)= strsgb(istr1,igpb)
            enddo
            call 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)
            do i = 1,nnod
               iflnod(i) = kfx(iflnod(i))
            enddo
            call formts(msizl,mpoin,ng,stotl,index,idofh,idiag,
     $           mnod,nnod,iflnod,ekmat)
            call formtf(ng,mpoin,mnod,nnod,iflnod,
     $           eforc,tforc)
 120     continue
 100  continue
      return
      end
*     
**********************************************************************
*     Making Elemental Stiffness Matrix
**********************************************************************
      subroutine estiff(ng,mnod,nnod,ntype,props,dm,cm,kk,
     $     nstre,nstr1,gepst,ncrit,hdpara,stres,iincs,
     $     xi,www,aa,ai,dadx,bl,dldx,bmat,ekmat,phi,geppre)
      implicit real*8(a-h,o-z)
      dimension  props(7)
      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  bmat(3,mnod*ng), ekmat(mnod*ng,mnod*ng)
      dimension  phi(2,mnod*ng)
      dimension  dmat(4,4)
*
*---- for plastic-elastic problem
*
      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 bmatrx(ng,mnod,md,nnod,bl,ai,dadx,dldx,bmat)
      call dmatrx(ntype,props,nstrs,dmat)
*
*---- making plastic stiffness matrix (D sup p)
*
         if(iincs.ne.1.and.gepst.ne.0.0d0.
     $     and.geppre.ne.0.0d0) then