jinteg.f

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

************************************************************************
*     J integral
************************************************************************
      subroutine jinteg(mpoin,mef,msrch,mtail,mnod,npoin,nef,ng,nadrs,
     $     nodf,nifcls,aldm,dm,cm,kk,alcm,x,jac,isrch,ltail,iflnod,
     $     xi,bl,dldx,bmat,bmat2,phi,disp,tu,sts_node,eps_node,wn,
     $     ntip,ndiv,ndivr,ngauss,r,resj,tstar)
      implicit real*8(a-h,o-z)
      dimension  x(ng,mpoin+3),xi(ng,mnod)
      dimension  nodf(3,mef),jac(0:2*mpoin+1,3)
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  iflnod(mnod)
      dimension  bl(6,mnod),dldx(6,mnod,ng)
      dimension  bmat(3,mnod*ng),bmat2(3,mnod*ng)
      dimension  disp(mnod*ng),tu(mpoin*ng)
      dimension  sts_node(4,mpoin),eps_node(4,mpoin)
      dimension  phi(2,mnod*ng),wn(mpoin)
*
*---- Set coordinate for center of Circle using J-integral
*
      xj=x(1,ntip)
      yj=x(2,ntip)
      call jcalc(mpoin,mef,msrch,mtail,mnod,npoin,nef,ng,nadrs,nodf,
     $     x,jac,isrch,ltail,iflnod,xi,bl,dldx,bmat,bmat2,phi,
     $     disp,tu,nifcls,aldm,dm,cm,kk,alcm,sts_node,eps_node,xj,yj,
     $     ndiv,ngauss,r,wn,resj)
      call tcalc(mpoin,mef,msrch,mtail,mnod,npoin,nef,ng,nadrs,nodf,
     $     x,jac,isrch,ltail,iflnod,xi,bl,dldx,bmat,phi,
     $     nifcls,aldm,dm,cm,kk,alcm,sts_node,xj,yj,ndivr,
     $     ndiv,r,resj,eps_node,wn,tstar)
      return
      end
*
************************************************************************
*     Read data for J integral
************************************************************************
      subroutine jin(jdata,rj)
      implicit real*8(a-h,o-z)
      dimension  jdata(5),rj(10)
*
*---- read input data for fracture mechanics parameters
*
      read(1600,*) (jdata(i),i=1,5)
      do i= 1,jdata(1)
         read(1600,*) rj(i)
      enddo
      return
      end
*
************************************************************************
*     calculating J integral
************************************************************************
      subroutine jcalc(mpoin,mef,msrch,mtail,mnod,npoin,nef,ng,nadrs,
     $     nodf,x,jac,isrch,ltail,iflnod,xi,bl,dldx,bmat,bmat2,phi,
     $     disp,tu,nifcls,aldm,dm,cm,kk,alcm,sts_node,eps_node,xj,yj,
     $     ndiv,ngauss,r,wn,resj)
      implicit real*8(a-h,o-z)
      parameter(pi=3.141592653589793d0)
      dimension  x(ng,mpoin+3),xi(ng,mnod),xx(2)
      dimension  nodf(3,mef),jac(0:2*mpoin+1,3)
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  ndelm(3),iflnod(mnod), nflnod(0:10)
      dimension  aa(6,6), ai(6,6), dadx(6,6,2)
      dimension  bl(6,mnod),dldx(6,mnod,ng)
      dimension  bmat(3,mnod*ng),bmat2(3,mnod*ng)
      dimension  sts(4),disp(mnod*ng),tu(mpoin*ng)
      dimension  t(2),ep(4),ep2(4)
      dimension  sts_node(4,mpoin),eps_node(4,mpoin)
      dimension  phi(2,mnod*ng)
      dimension  wn(mpoin)
      common/gaussa/ssg(10,10),wwg(10,10)
      md = 3
      resj = 0.0d0
      ang = pi / dble(ndiv)
      subc = ang/2.0d0
      do idiv = 1,ndiv
         addc = (ang*(2.0d0*dble(idiv)-1.0d0))/2.0d0
         do igauss = 1,ngauss
            cta = addc - subc * ssg(igauss,ngauss)
*
*---- set coordinates
*
            xx(1) = xj + r * dcos(cta)
            xx(2) = yj + r * dsin(cta)
*
*---- calculate stress and bmatrix
*
            call locate2(xx(1),xx(2),x,nodf,jac,mpoin,mef,
     $           nef,npoin,itri)
            do inode= 1,3
               ndelm(inode)= nodf(inode,itri)
            enddo
            call getndf(mpoin,ng,msrch,mtail,mnod,
     $           npoin,nadrs,isrch,ltail,ndelm,
     $           nifcls,nnod,iflnod,nflnod,x,xx,aldm,ddmax)
            do i = 1,nnod
               xi(1,i) = x(1,iflnod(i)) - xx(1)
               xi(2,i) = x(2,iflnod(i)) - xx(2)
            enddo
            dm = ddmax
            cm = alcm*dm
            if(dm.lt.1.0d-30) then
               write(7777,*) 'Error Message from subroutine jcalc 1 !'
               write(7777,*) 'Diameter of infulence region is zero !'
               stop
            endif
            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)
            do i = 1,4
*---- Zero clear
               ep(i) = 0.0d0
               sts(i) = 0.0d0
*---- Calculate eps,stress
               do inod= 1,nnod
                  ep(i) = ep(i) + phi(1,inod*2-1)*
     $                 eps_node(i,iflnod(inod))
                  sts(i) = sts(i) + phi(1,inod*2-1)*
     $                 sts_node(i,iflnod(inod))
               enddo
            enddo
*     
*---- Calculate du/dy , dv/dx
*
            do i = 1,nnod
               disp(i*2-1) = tu(iflnod(i)*2-1)
               disp(i*2) = tu(iflnod(i)*2)
            enddo
*     
            do inod = 1,nnod
               bmat2(1,(inod-1)*2+1)= bmat(2,(inod-1)*2+2)
               bmat2(2,(inod-1)*2+2)= bmat(1,(inod-1)*2+1)
               bmat2(3,(inod-1)*2+1)= bmat2(2,(inod-1)*2+2)
               bmat2(3,(inod-1)*2+2)= bmat2(1,(inod-1)*2+1)
            enddo
            do j= 1,3
*---- Zero clear
               ep2(j) = 0.0d0
               do inod= 1,nnod*ng
                  ep2(j)= ep2(j) + bmat2(j,inod)*disp(inod)
               enddo
            enddo
*
*---- Calculate strain energy density
*
            w = 0.0d0
            do inod= 1,nnod
               w = w + phi(1,inod*2-1) * wn(iflnod(inod))
            enddo
*
*---- Calculate J-integral
*
            t(1)= sts(1)*dcos(cta)+sts(3)*dsin(cta)
            t(2)= sts(2)*dsin(cta)+sts(3)*dcos(cta)
            func = w * dcos(cta) - (t(1)*ep(1)+t(2)*ep2(2))
            www = 0.5d0 * r * ang * wwg(igauss,ngauss)
            resj = resj + func * www
         enddo
      enddo
      resj = resj * 2.0d0
      return
      end
*
************************************************************************
*     Calculate T*
************************************************************************
      subroutine tcalc(mpoin,mef,msrch,mtail,mnod,npoin,nef,ng,
     $     nadrs,nodf,x,jac,isrch,ltail,iflnod,xi,bl,dldx,bmat,
     $     phi,nifcls,aldm,dm,cm,kk,alcm,sts_node,
     $     xj,yj,ndivr,ndiv,r,resj,eps_node,wn,tstar)
      implicit real*8(a-h,o-z)
      parameter (pi=3.141592653589793d0)
      parameter (gline=0.0d0)
      dimension  x(ng,mpoin+3),xi(ng,mnod),xx(2)
      dimension  nodf(3,mef),jac(0:2*mpoin+1,3)
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  ndelm(3),iflnod(mnod), nflnod(0:10)
      dimension  aa(6,6), ai(6,6), dadx(6,6,2)
      dimension  bl(6,mnod),dldx(6,mnod,ng)
      dimension  bmat(3,mnod*ng)
      dimension  sts_node(4,mpoin),eps_node(4,mpoin)
      dimension  phi(2,mnod*ng)
      dimension  sts(4),dqeps(4)
      dimension  wn(mpoin)
      common /integ/xq(4,2),ww(4)
      md=3
      cang= (pi-gline) / dble(ndiv)
      rang= r/dble(ndivr)
      ttt=0.0d0
      do 10 idivr=1,ndivr
*---- loop over r direction
         r2= rang*dble(idivr)
         r1= r2-rang
*
         do 100 idivc= 1,ndiv
            c2= cang*dble(idivc)
            c1= c2-cang
            do 1000 igauss= 1,4
*---- loop over gauss integral (4 gaussian point only)
               r = 0.5d0 * (r1+r2 - (r1-r2)*xq(igauss,1))
               c = 0.5d0 * (c1+c2 - (c1-c2)*xq(igauss,2))
               xx(1)= xj + r * dcos(c)
               xx(2)= yj + r * dsin(c)
*
*---- calculate stress and bmatrix
*
               call locate2(xx(1),xx(2),x,nodf,jac,mpoin,
     $              mef,nef,npoin,itri)
               do inode= 1,3
                  ndelm(inode)= nodf(inode,itri)
               enddo
               call getndf(mpoin,ng,msrch,mtail,mnod,
     $              npoin,nadrs,isrch,ltail,ndelm,
     $              nifcls,nnod,iflnod,nflnod,x,xx,
     $              aldm,ddmax)
               do i = 1,nnod
                  xi(1,i) = x(1,iflnod(i)) - xx(1)
                  xi(2,i) = x(2,iflnod(i)) - xx(2)
               enddo
               dm = ddmax
               cm = alcm*dm
               if(dm.lt.1.0d-30) then
                  write(7777,*) 
     $                 'Error Message from subroutine tcalc 1 !     '
                  write(7777,*)
     $                 '     Diameter of infulence region is zero !'
                  stop
               endif
               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)
*---- Zero clear
               do i = 1,4
                  sts(i) = 0.0d0
                  dqeps(i) = 0.0d0
                  do inod= 1,nnod
*---- Calculate stress
                     sts(i) = sts(i) + phi(1,inod*2-1)
     $                    * sts_node(i,iflnod(inod))
                     dqeps(i) = dqeps(i) + bmat(1,inod*2-1)
     $                    * eps_node(i,iflnod(inod))
                  enddo
               enddo
*     
*---- Calculate dw/dx
*     
               dqwn = 0.0d0
               do inod= 1,nnod
                  dqwn = dqwn + bmat(1,inod*2-1)
     $                 * wn(iflnod(inod))
               enddo
*     
*---- Calculate T*
*     
               tfunc = dqwn
               do i=1,4
                  tfunc = tfunc - sts(i) * dqeps(i)
               enddo
               wwj = 0.25d0 * r * cang * rang
               www = wwj * ww(igauss)
               ttt = ttt + tfunc * www
 1000       continue
*---- end of gauss integral loop
 100     continue
*---- end of theta-element loop
 10   continue
*---- end of r-element loop
      tstar = resj - 2.0d0 * ttt
      return
      end
*
**********************************************************************
* integral location and weight for Gauss integration
**********************************************************************
      block data gaussrec
      implicit real*8(a-h,o-z)
      common /integ/xq(4,2),ww(4)
      data (xq(i,1),i=1,4)
     $     /-0.577350269189626d0,0.577350269189626d0,
     $     0.577350269189626d0,-0.577350269189626d0/
      data (xq(i,2),i=1,4)      
     $     /-0.577350269189626d0,-0.577350269189626d0,
     $     0.577350269189626d0,0.577350269189626d0/
      data (ww(i),i=1,4) /1.0d0,1.0d0,1.0d0,1.0d0/
      end