frmain.f

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

**********************************************************************
* *                                                                * *
* *   This code is 2-dimensional elast-plastic analysis by EFGM    * *
* *                                                                * *
**********************************************************************
      program efgm_plast
      implicit real*8(a-h,o-z)
      parameter(mpoin=3000,ng=2)
      parameter(mef=mpoin*2+1,meb=300)
      parameter(msrch=mpoin*10,mtail=100,mnod=1000)
      parameter(mgpf=mef*7,mgpb=meb*7,mavnd=30)
      parameter(aldm= 1.8d0, bldm= 2.5d0, alcm= 0.25d0, kk=1)
      parameter(coplty=1.0d8)
      parameter(ngfcls= 1, nifcls= 4)
      parameter(ngbcls= 4, nibcls= 4)
      parameter(kbd=10, kcm=400)
      parameter(ircm=1)
      parameter(msizl=400000,mindx=mnod*ng)
      integer  npoin,ntype
      integer  nef,neb
      dimension  x(ng,mpoin+3)
      dimension  props(7)
      dimension  nodf(3,mef), nodb(2,meb)
      dimension  nbcbm(3,meb), dbcbm(4,meb)
*
*---- dimensions for searching nodes under influence region
*
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  iflnod(mnod), nflnod(0:10)
      dimension  xxgpf(ng,mgpf), wwgpf(mgpf), dmaxf(mgpf)
      dimension  ndgpf(0:mgpf),  idgpf(mgpf*mavnd)
      dimension  xxgpb(ng,mgpb), wwgpb(mgpb), dmaxb(mgpb)
      dimension  ndgpb(0:mgpb),  idgpb(mgpb*mavnd)
*
*---- dimensions for making elemental stiffness matrix
*
      dimension  xi(ng,mnod)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  bmat(3,mnod*ng)
      dimension  ekmat(mnod*ng,mnod*ng),eforc(mnod*ng)
      dimension  phi(2,mnod*ng)
*
*---- dimensions for delaunay triangulation
*
      dimension  ibex(kbd),ibin(kbd),ibno(kbd,mpoin)
      dimension  ind(kbd+1),idm(2*mpoin+1)
      dimension  mtj(2*mpoin+1,3),jac(0:2*mpoin+1,3)
      dimension  jnb(mpoin),nei(mpoin,kcm),ifix(mpoin),list(mpoin)
      dimension  istack(mpoin),kv(2*mpoin+1),ibr(2*mpoin+1)
      dimension  iadres(mpoin+3),map(2*mpoin+1)
      dimension  angl(2*mpoin+1),nedg(2*mpoin+1),kstack(2*mpoin+1)
*
*---- dimensions for calculating stress
*
      dimension  disp(mnod*ng)
*
*---- dimensions for renumbering table of nodes by RCM
*
      dimension kfx(mpoin), kfdx(mpoin)
*
*     kfx(): Original No. ---> Renumbered No.
*     kfdx(): Renumbered No. ---> Original No.
*
*---- dimensions for working area for Reversed Cuthill-Mckee (RCM) method
*     
      dimension  ma(0:mpoin),ir(mpoin),iw(mpoin)
*
*---- dimensions for stiffness matrix (skyline storage)
*
      dimension  stotl(msizl)
      dimension  tu(mpoin*ng), pload(mpoin*ng)
*
*---- dimensions for skyline method 
*
      integer index(mindx),lhigh(mindx)
      integer idofh(mpoin*ng),idiag(0:mpoin*ng)
*
*      index(i): system dof(degree of freedom) no. of EFGM dof i
*      lhigh(i): row width of EFGM dof i for triangle matrix
*      idofh(i): row width of system dof i for triangle matrix
*      idiag(i): location of diagonal of i-th equation for matrix
*
*---- dimensions for working area for solver
*
      dimension ww(mpoin*ng)
*
*---- dimensions for plastic-elastic problem
*
      dimension  epbc(4,meb)
      dimension  ncrit(3)
*
      dimension  hdpara(20)
      dimension  effstf(mgpf),epstnf(mgpf),strsgf(4,mgpf),epsgf(4,mgpf)
      dimension  epstnfpre(mgpf),epstnbpre(mgpb)
      dimension  effstb(mgpb),epstnb(mgpb),strsgb(4,mgpb)
      dimension  eload(mpoin*ng),asdis(mpoin*ng),bsdis(mpoin*ng)
      dimension  rload(mpoin*ng),tload(mpoin*ng)
*
      character  header*70
*
*---- dimensions for post process for J integral and T*
*
      dimension  eps_node(4,mpoin),sts_node(4,mpoin)
      dimension  wn(mpoin)
      dimension  bmat2(3,mnod*ng)
      dimension  jdata(5),rj(10)
      character  ch*10
      ij=0
*
*---- open files statement
*
      open(5,file='efgmstdin.txt')
      open(1500,file='efgmplast.txt')
      open(1600,file='efgmjread.txt',status='old',
     $     err=100,form='formatted',access='sequential')
      ij=1
 100  continue
      open(7777,file='efgmerror.txt')
*
*---- Read input data
*
      call gdata1(ng,mpoin,kbd,meb,ntype,nalgo,props,nbc,
     $     nbcbm,dbcbm,x,nex,nin,ibex,ibin,ibno,nob,nib,npoin,
     $     nszf,ncrit,hdpara,nincs,neb,nodb)
      if(ij.eq.1) then
         call jin(jdata,rj)
         do j= 1,jdata(1)
            write(ch,'(i6,4h.txt)') 10000+j
            ch(1:4)= 'JOUT'
            open(1000+j,file=ch)
         enddo
      endif
*
*---- Delaunay triangulation
*
      call auto(ng,npoin,nef,x,mpoin,mef,kbd,kcm,
     $     nodf,ibex,ibin,ibno,nex,nin,nob,nib,ind,
     $     idm,mtj,jac,jnb,nei,ifix,list,istack,kv,ibr,
     $     iadres,map,angl,nedg,kstack)
*
*---- Making search table
*
      call mkstbl(mef,msrch,npoin,nef,nodf,nadrs,isrch)
*
*---- Renumbering no. of nodes by reversed Cuthill-Mckee method 
*
      call rcm(mpoin,msrch,mtail,npoin,nadrs,isrch,
     $     ircm,ltail,ma,ir,iw,kfx,kfdx)
*
*---- Making list of nodes under influence region (in domain)
*
      call mkitgf(mpoin,mef,ng,msrch,mtail,mnod,
     $     npoin,nef,nadrs,x,nodf,isrch,ltail,
     $     ngfcls,nifcls,iflnod,nflnod,aldm,
     $     mgpf,mavnd,ngpf,xxgpf,wwgpf,dmaxf,ndgpf,idgpf)
*
*---- Making list of nodes under influence region (on boundary)
*
      call 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)
*
*---- Prepare skyline storage address
*
      call skysub(mpoin,mnod,mgpf,mgpb,mavnd,msizl,ng,nszf,idofh,
     $     iflnod,ngpf,ngpb,ndgpf,idgpf,ndgpb,idgpb,index,lhigh,
     $     idiag,kfx)
*
*---- Making total load vector
*
      call 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)
*
*---- Loop over each increment
*
      do 10000 iincs= 1,nincs
         call zerod(nszf,bsdis)
*
*---- Read data for current increment
*
         call increm(iincs,nbc,nszf,meb,mpoin,ng,dbcbm,epbc,
     $        rload,tload,eload,facto,tfact,toler,miter,noutp)
         header='     ***** incremental external force *****'
         call ldout(meb,nbc,nbcbm,epbc,header)
*     
*---- Loop over each iteration
*
         do 20000 iiter= 1,miter
            call algor(iincs,iiter,nalgo,kresl)
*---- Set displacement to zero in iteration loop.
            if(iiter.ne.1) then
               do ibc= 1,nbc
                  do i = 1,4
                     epbc(i,ibc) = 0.0d0
                  enddo
               enddo
            endif
*     
*---- Domain integral for making stiffness matrix
*     
            if(kresl.eq.1) then
               call zerod(msizl,stotl)
               call 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,epstnf,strsgf,iincs,
     $              gepst,epstnfpre)
*     
*---- Boundary integral for making stiffness matrix
*     
               pnlty= coplty*amxdia
               call bunint(mpoin,meb,ng,mnod,x,nodb,ngbcls,
     $              iflnod,ntype,props,alcm,kk,nibcls,xi,bl,
     $              dldx,phi,ekmat,eforc,eload,nbc,nbcbm,
     $              epbc,pnlty,mgpb,mavnd,dmaxb,ndgpb,idgpb,
     $              msizl,mindx,stotl,index,idofh,idiag,kfx,
     $              ncrit,hdpara,epstnb,strsgb,iincs,gepst,
     $              epstnbpre)
               call skydec(msizl,nszf,stotl,idiag,ww,eps,ind)
            endif
*     
*---- Skyline solver with modified cholesky decomposition
*
            call skysol(msizl,nszf,stotl,idiag,eload)
            call renode(mpoin,ng,npoin,nszf,kfdx,eload,pload)
*     
*---- Calculate residual forces (for domain)
*     
            call domresidu(mpoin,mnod,mgpf,mavnd,ng,nszf,
     $           ntype,ngpf,x,tu,asdis,bsdis,eload,effstf,disp,iflnod,
     $           xxgpf,wwgpf,dmaxf,ndgpf,idgpf,dldx,bl,bmat,eforc,
     $           alcm,kk,xi,kfx,epstnf,strsgf,epsgf,props,ncrit,hdpara,
     $           epstnfpre)
            call bunresidu(mpoin,mnod,mgpb,mavnd,ng,
     $           ntype,ngpb,x,effstb,disp,iflnod,xxgpb,
     $           asdis,dmaxb,ndgpb,idgpb,dldx,bl,bmat,
     $           alcm,kk,xi,epstnb,strsgb,props,ncrit,
     $           hdpara,epstnbpre)
*
*---- Check for convergence
*
            call conver(nszf,asdis,bsdis,eload,tload,
     $           iiter,nchek,toler,pvalu)
*     
*---- If solution has converged stop iterating and output results
*     
            if(nchek.eq.0) then
               goto 600
            elseif(nchek.eq.999) then
               write(7777,*) 'Error Message from subroutine main 1 !'
               write(7777,*) 'nchek is 999 !!!!!!'
               stop
            endif
20000    continue
 600     continue
*     
*---- Calculate nodal data.
*     
         call smooth(mpoin,mef,mnod,mgpf,msrch,
     $        mtail,ng,nef,npoin,ngpf,nadrs,x,nodf,jac,tu,
     $        disp,iflnod,isrch,ltail,nifcls,xxgpf,dldx,bl,
     $        bmat,phi,strsgf,epsgf,bldm,dm,cm,alcm,kk,xi,
     $        eps_node,sts_node,wn)
 60      format(i5,6e15.7)
         if(noutp.eq.1.or.noutp.eq.2.or.noutp.eq.3) then
*
*---- Output results for step if required
*     
            call output(ng,mpoin,npoin,mgpf,ngpf,iincs,noutp,x,xxgpf,
     $           tu,sts_node,eps_node,strsgf,epsgf,effstf,epstnf)
         endif
         if(ij.eq.1) then
*=====================================================================
*     Post Process for Fracture Mechanics
*
*     1. J-integral for plast-elastic problem
*     2. T* for plast-elastic problem
*
*=====================================================================
            call postfr(mpoin,ng,x,mef,msrch,mtail,mnod,npoin,nef,
     $           nadrs,nodf,jac,isrch,ltail,iflnod,xi,bl,dldx,bmat,
     $           bmat2,phi,ntype,nifcls,aldm,dm,cm,kk,alcm,tu,disp,
     $           sts_node,eps_node,wn,jdata,rj)
         endif
10000 continue
      close(1500)
      close(1600)
      close(7777)
      if(ij.eq.1) then
         do j= 1,jdata(1)
            close(1000+j)
         enddo
      endif
      stop
      end