smooth.f

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

**********************************************************************
*     This subroutine calculate nodal data.
*     eps,stress,W
**********************************************************************
      subroutine 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,aldm,dm,cm,alcm,kk,xi,
     $     eps_node,sts_node,wn)
      implicit real*8(a-h,o-z)
      dimension  x(ng,mpoin+3)
      dimension  nodf(3,mef),jac(0:2*mpoin+1,3)
      dimension  isrch(2,msrch), ltail(mtail)
      dimension  iflnod(mnod)
      dimension  tu(mpoin*ng)
      dimension  bl(6,mnod), dldx(6,mnod,ng)
      dimension  xi(ng,mnod),bmat(3,mnod*ng)
      dimension  xxgpf(ng,mgpf),strsgf(4,mgpf),epsgf(4,mgpf)
      dimension  disp(mnod*ng)
      dimension  phi(2,mnod*ng)
      dimension  eps_node(4,mpoin),sts_node(4,mpoin)
      dimension  wn(mpoin)
      dimension  eps(4),sts(4)
      do i=1,npoin
         xt = x(1,i)
         yt = x(2,i)
         call stress(xt,yt,mpoin,ng,mef,nef,mnod,mgpf,ngpf,
     $        npoin,nodf,iflnod,x,aldm,dm,cm,kk,alcm,xi,dldx,
     $        bl,phi,jac,xxgpf,strsgf,epsgf,sts,eps)
         do j=1,4
            eps(j) = eps(j) - eps_node(j,i)
            sts(j) = 0.5d0 * (sts(j) + sts_node(j,i))
         enddo
         w = 0.0d0
         do j=1,4
            w = w + sts(j) * eps(j)
         enddo
         do j=1,4
            eps_node(j,i) = eps_node(j,i) + eps(j)
            sts_node(j,i) = 2.0d0 * sts(j) - sts_node(j,i)
         enddo
         wn(i) = wn(i) + w
      enddo
      return
      end
*
**********************************************************************
*     Calculate stress on nodes.
**********************************************************************
      subroutine stress(xt,yt,mpoin,ng,mef,nef,mnod,mgpf,ngpf,
     $     npoin,nodf,iflnod,x,aldm,dm,cm,kk,alcm,xi,dldx,bl,phi,
     $     jac,xxgpf,strsgf,epsgf,sts,eps)
      implicit real*8(a-h,o-z)
      dimension  jac(0:2*mpoin+1,3)
      dimension  nodf(3,mef)
      dimension  xxgpf(ng,mgpf)
      dimension  x(ng,mpoin),xx(2),xi(ng,mnod)
      dimension  ndelm(3),iflnod(mnod)
      dimension  aa(6,6),ai(6,6), bl(6,mnod)
      dimension  dadx(6,6,2), dldx(6,mnod,ng)
      dimension  strsgf(4,mgpf),sts(4)
      dimension  epsgf(4,mgpf),eps(4)
      dimension  phi(2,mnod*ng)
      md = 3
      xx(1) = xt
      xx(2) = yt
      call locate2(xx(1),xx(2),x,nodf,jac,mpoin,mef,
     $     nef,npoin,itri)
      do inode= 1,3
         ndelm(inode)= nodf(inode,itri)
      enddo
      bldm= aldm*2.0d0/3.0d0
      dd1= sqrt((x(1,ndelm(1))-xx(1))**2+
     $     (x(2,ndelm(1))-xx(2))**2 )
      dd2= sqrt((x(1,ndelm(2))-xx(1))**2+
     $     (x(2,ndelm(2))-xx(2))**2 )
      dd3= sqrt((x(1,ndelm(3))-xx(1))**2+
     $     (x(2,ndelm(3))-xx(2))**2 )
      ddmax= bldm*max(dd1,dd2,dd3)
      nnod = 0
      do igpf=1,ngpf
         dd= sqrt((xxgpf(1,igpf)-xx(1))**2+
     $        (xxgpf(2,igpf)-xx(2))**2 )
         if(dd.le.ddmax) then
            nnod = nnod+1
            if(nnod.gt.mnod) then
               write(7777,*) 
     $              'Error Message from subroutine stress 1 !    '
               write(7777,*) '     mnod is smaller than nnod !'
               write(7777,*) 'mnod=',mnod,' nnod=',nnod
               stop
            endif
            iflnod(nnod) = igpf
         endif
      enddo
      do i = 1,nnod
         xi(1,i) = xxgpf(1,iflnod(i)) - xx(1)
         xi(2,i) = xxgpf(2,iflnod(i)) - xx(2)
      enddo
      dm = ddmax
      cm = alcm*dm
      if(dm.lt.1.0d-30) then
         write(7777,*) 'Error Message from subroutine stress 2 !    '
         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)
      do i=1,4
         sts(i) = 0.0d0
         eps(i) = 0.0d0
      enddo
      do i = 1,4
         do inod= 1,nnod
            sts(i) = sts(i) + phi(1,inod*2-1)*strsgf(i,iflnod(inod))
            eps(i) = eps(i) + phi(1,inod*2-1)*epsgf(i,iflnod(inod))
         enddo
      enddo
      return
      end
*
**********************************************************************
*     locate triangle which encloses point (xp,yp)
*     Using for arbitrary domain (include hollow)
**********************************************************************  
      subroutine locate2(xp,yp,x,nodf,jac,mpoin,mef,nef,npoin,itri)
      implicit real*8(a-h,o-z)  
      dimension x(2,mpoin+3),nodf(3,mef),jac(0:mef,3)
      itri = 1
 10   continue
      do i = 1,3
         ia = nodf(i,itri)
         ib = nodf(mod(i,3)+1,itri)
         if((x(2,ia)-yp)*(x(1,ib)-xp).gt.
     &        (x(1,ia)-xp)*(x(2,ib)-yp)) then
            itri = jac(itri,i)
            if(itri.eq.0) then
               goto 20
            else
               go to 10
            endif
         endif
      enddo
 20   continue
      if(itri.ne.0.and.itri.le.npoin) then
         return
      else
         do n=1,nef
            ia = nodf(1,n)
            ib = nodf(2,n)
            ic = nodf(3,n)
            if(((x(2,ia)-yp)*(x(1,ib)-xp).le.
     $           (x(1,ia)-xp)*(x(2,ib)-yp))) then
               if(((x(2,ib)-yp)*(x(1,ic)-xp).le.
     $              (x(1,ib)-xp)*(x(2,ic)-yp))) then
                  if(((x(2,ic)-yp)*(x(1,ia)-xp).le.
     $                 (x(1,ic)-xp)*(x(2,ia)-yp))) then
                     itri = n
                     return
                  endif
               endif
            endif
         enddo
      endif
      write(7777,*) 'Error Message from subroutine locate2 1 !    '
      write(7777,600) xp,yp
      stop
 600  format(8x,'x=',e15.7,8x,'y=',e15.7)
      end