text2.f90

fortran程序

       40      s = 1.0/sqrt(1.0+b*b);  c = b*s; c2 = c*c ;s2 = s*s  ; cs = c*s
               u = c2*a(j-1,1) - 2.0*cs*a(j-1,2) + s2*a(j,1)
               u1 = s2*a(j-1,1) + 2.0*cs*a(j-1,2) + c2*a(j,1)
               a(j-1,2) = cs*(a(j-1,1)-a(j,1)) + (c2-s2)*a(j-1,2)
               a(j-1,1) = u        ;   a(j,1) = u1 ;     j2 = j - 2
               do  l=iugl,j2
                  jl = j - l      ;   u = c*a(l,jl) - s*a(l,jl+1)
                  a(l,jl+1) = s*a(l,jl) + c*a(l,jl+1) ;   a(l,jl) = u
               end do
               jm = j - iw
               if (j/=kr) a(jm-1,iw+1) = c*a(jm-1,iw+1) - s*g
               maxl = iw - 1     ; if (n-j<iw-1) maxl = n - j
               if (maxl>0) then
                do  l=1,maxl
                  u = c*a(j-1,l+2) - s*a(j,l+1)
                  a(j,l+1) = s*a(j-1,l+2) + c*a(j,l+1) ;  a(j-1,l+2) = u
                end do
               end if
            if (j+iw>n) go to 120
               g = -s*a(j,iw+1)   ;   a(j,iw+1) = c*a(j,iw+1)
     120    continue
     140    continue
     160 continue
   end if
    e(1) = 0.0
      d(1:n) = a(1:n,1)  ;   if (2>n) go to 240
      do  i=2,n    ;    e(i) = a(i-1,2)  ;   end do
    240  e2 = e*e
   return
 end subroutine bandred
 subroutine bisect(d,e,acheps,ifail)
!
!      this subroutine finds the eigenvalues of a tridiagonal
!      matrix, given with its diagonal elements in the array d(n) and
!      its subdiagonal elements in the last n - 1 stores of the
!      array e(n), using ql transformations. the eigenvalues are
!      overwritten on the diagonal elements in the array d in
!      ascending order. the subroutine will fail if any one
!      eigenvalue takes more than 30 iterations.
!
 implicit none
    real,intent(in)::acheps; real,intent(in out)::d(0:),e(0:)
    integer,intent(in out)::ifail  
      integer ::m, isave, n, i, l, j, i1, m1, ii
      real    :: b, f, h, g, p, r, c, s
      isave = ifail    ;  n = ubound(d,1)
      if (n==1) go to 40
      do  i=2,n ;  e(i-1) = e(i)  ; end do
    40 e(n) = 0.0
      b = 0.0    ;   f = 0.0
      do 340 l=1,n
         j = 0    ;   h = acheps*(abs(d(l))+abs(e(l)))
         if (b<h) b = h
!     look for small sub diagonal element
         do  m=l,n ;  if (abs(e(m))<=b) exit   ; end do
         if (m==l) go to 260
    100  if (j==30) go to 360  ;  j = j + 1
!     form shift
         g = d(l)    ;    h = d(l+1) - g
         if (abs(h)>=abs(e(l))) go to 120
         p = h*0.5/e(l) ;  r = sqrt(p*p+1.0) ;  h = p + r
         if (p<0.0) h = p - r  ;  d(l) = e(l)/h ;   go to 140
    120  p = 2.0*e(l)/h
         r = sqrt(p*p+1.0)  ;   d(l) = e(l)*p/(1.0+r)
    140  h = g - d(l)
         i1 = l + 1    ;  if (i1>n) go to 180
         do i=i1,n  ;   d(i) = d(i) - h   ;  end do
   180   f = f + h
!     ql transformation
         p = d(m)  ;  c = 1.0 ;   s = 0.0 ;  m1 = m - 1
         do 240 ii=l,m1
            i = m1 - ii + l  ;  g = c*e(i) ;   h = c*p
            if (abs(p)<abs(e(i))) go to 200
            c = e(i)/p    ;  r = sqrt(c*c+1.0)  ;   e(i+1) = s*p*r
            s = c/r  ;  c = 1.0/r  ;    go to 220
    200     c = p/e(i)  ; r = sqrt(c*c+1.0);   e(i+1) = s*e(i)*r
            s = 1.0/r   ;    c = c/r
    220     p = c*d(i) - s*g  ;  d(i+1) = h + s*(c*g+s*d(i))
   240    continue
         e(l)= s*p  ;   d(l)= c*p   ;  if (abs(e(l))>b) go to 100
   260    p = d(l) + f
!     order eigenvalue
         if (l==1) go to 300
         do  ii=2,l  ;   i = l - ii + 2
            if (p>=d(i-1)) go to 320  ;    d(i) = d(i-1)
         end do
   300    i = 1      ; 320    d(i) = p
   340 continue
      ifail = 0     ;    return
  360 ifail=1       ;    return
 end subroutine bisect
 subroutine formlump(diag,emm,g)
!
!      this subroutine forms the global mass matrix as vector diag
!
  implicit none
   real,intent(in)::emm(:,:)  ; real,intent(out)::diag(0:)
   integer,intent(in)::g(:); integer::i,ndof; ndof=ubound(emm,1)
     do  i=1,ndof;   diag(g(i))=diag(g(i))+emm(i,i) ; end do
    return
 end subroutine formlump   
   function bandwidth(g) result(nband)
   ! finds the element bandwidth from g
   implicit none     ; integer :: nband
   integer,intent(in)::g(:) 
      nband= maxval(g,1,g>0)-minval(g,1,g>0)
   end function bandwidth
   subroutine formnf(nf)
   ! reform nf
   implicit none
   integer,intent(in out)::nf(:,:)
    integer:: i,j,m
    m=0
    do j=1,ubound(nf,2)
       do i=1,ubound(nf,1)
          if(nf(i,j)/=0) then
             m=m+1; nf(i,j)=m
          end if
       end do
    end do
    return
  end subroutine formnf 
  subroutine invert(matrix)
  ! invert a small square matrix onto itself
  implicit none
  real,intent(in out)::matrix(:,:)
  integer::i,k,n; real::con  ; n= ubound(matrix,1)
  do k=1,n
     con=matrix(k,k); matrix(k,k)=1.
     matrix(k,:)=matrix(k,:)/con
     do i=1,n
        if(i/=k) then
           con=matrix(i,k); matrix(i,k)=0.0
           matrix(i,:)=matrix(i,:) - matrix(k,:)*con
        end if
     end do
  end do
  return
 end subroutine invert
 function determinant (jac) result(det)
 ! returns the determinant of a 1x1 2x2 3x3 jacobian matrix 
 implicit none    ; real :: det
 real,intent(in)::jac(:,:); integer:: it ; it = ubound(jac,1)  
 select case (it)
   case (1)
    det=1.0
   case (2)
    det=jac(1,1)*jac(2,2) - jac(1,2) * jac(2,1)
   case (3)
    det= jac(1,1)*(jac(2,2) * jac(3,3) -jac(3,2) * jac(2,3))
    det= det-jac(1,2)*(jac(2,1)*jac(3,3)-jac(3,1)*jac(2,3))
    det= det+jac(1,3)*(jac(2,1)*jac(3,2)-jac(3,1)*jac(2,2))
   case default
    print*,' wrong dimension for jacobian matrix'
 end select
 return
 end function determinant
 subroutine deemat(dee,e,v)
 ! returns the elastic dee matrix for given ih
 ! ih=3,plane strain; =4,axisymmetry or plane strain elastoplasticity
 ! =6 , three dimensional
 implicit none
   real,intent(in)::e,v; real,intent(out)::dee(:,:)
 ! local variables
   real::v1,v2,c,vv; integer :: i,ih;  dee=0.0  ; ih = ubound(dee,1)
         v1 = 1. - v; c = e/((1.+v)*(1.-2.*v))
   select case (ih)
          case(3)
             dee(1,1)=v1*c; dee(2,2)=v1*c; dee(1,2)=v*c; dee(2,1)=v*c
             dee(3,3)=.5*c*(1.-2.*v)
          case(4)
             dee(1,1)=v1*c; dee(2,2)=v1*c; dee(4,4)=v1*c
             dee(3,3)=.5*c*(1.-2.*v) ; dee(1,2)=v*c; dee(2,1)=v*c
             dee(1,4)=v*c; dee(4,1)=v*c; dee(2,4)=v*c; dee(4,2)=v*c
          case(6)
             v2=v/(1.-v); vv=(1.-2.*v)/(1.-v)*.5
             do i=1,3; dee(i,i)=1.;end do; do i=4,6; dee(i,i)=vv; end do
             dee(1,2)=v2; dee(2,1)=v2; dee(1,3)=v2; dee(3,1)=v2
             dee(2,3)=v2; dee(3,2)=v2
             dee = dee*e/(2.*(1.+v)*vv)
          case default
             print*,'wrong size for dee matrix'
   end select
 return
 end subroutine deemat    
 subroutine beemat(bee,deriv)
 ! bee matrix for 2-d elasticity or elastoplasticity (ih=3 or 4 respectively)
 ! or for 3-d (ih = 6)
 implicit none
 real,intent(in)::deriv(:,:);  real,intent(out)::bee(:,:)
 ! local variables
 integer::k,l,m,n , ih,nod; real::x,y,z
 bee=0. ; ih = ubound(bee,1); nod = ubound(deriv,2)
     select case (ih)
       case(3,4)
        do m=1,nod
           k=2*m; l=k-1; x=deriv(1,m); y=deriv(2,m)
           bee(1,l)=x; bee(3,k)=x; bee(2,k)=y; bee(3,l)=y
        end do
       case(6)
        do m=1,nod
           n=3*m;  k=n-1; l=k-1
           x=deriv(1,m); y=deriv(2,m); z=deriv(3,m)
           bee(1,l)=x; bee(4,k)=x; bee(6,n)=x
           bee(2,k)=y; bee(4,l)=y; bee(5,n)=y
           bee(3,n)=z; bee(5,k)=z; bee(6,l)=z
        end do
       case default
        print*,'wrong dimension for nst in bee matrix'        
      end select   
  return
 end subroutine beemat
 subroutine bmataxi(bee,radius,coord,deriv,fun)
 ! b matrix for axisymmetry
 real,intent(in)::deriv(:,:),fun(:),coord(:,:);real,intent(out)::bee(:,:),radius
 integer::nod ,k,l,m; real :: x,y
  radius = sum(fun * coord(:,1))  ; nod = ubound(deriv , 2) ; bee = .0
  do m = 1 , nod
   k=2*m; l = k-1 ; x = deriv(1,m); bee(1,l) = x; bee(3 , k) = x
   y = deriv(2,m); bee(2,k)=y; bee(3,l) = y; bee(4,l)=fun(m)/radius
  end do
  return
 end subroutine bmataxi

 subroutine sample(element,s,wt)
 ! returns the local coordinates of the integrating points
 implicit none 
  real,intent(out)::s(:,:),wt(:)  ; character(*),intent(in):: element
  integer::nip ;  real:: root3, r15 , w(3),v(9),b,c   
  root3 = 1./sqrt(3.)   ;  r15 = .2*sqrt(15.) 
  nip = ubound( s , 1 ) 
         w = (/5./9.,8./9.,5./9./); v=(/5./9.*w,8./9.*w,5./9.*w/)
     select case (element)
            case('line')
             select case(nip)
              case(1)
               s(1,1)=0.  ;  wt(1)=2.   
              case(2)
               s(1,1)=root3  ; s(2,1)=-s(1,1)  ;  wt(1)=1.  ; wt(2)=1.
              case(3)
               s(1,1)=r15 ; s(2,1)=.0     ; s(3,1)=-s(1,1)
               wt = w
              case(4)
               s(1,1)=.861136311594053  ; s(2,1)=.339981043584856  
               s(3,1)=-s(2,1)  ; s(4,1)=-s(1,1)
               wt(1)=.347854845137454 ; wt(2)=.652145154862546 
               wt(3)=wt(2) ; wt(4)=wt(1)
              case(5)
               s(1,1)=.906179845938664 ; s(2,1)=.538469310105683  
               s(3,1)=.0 ; s(4,1)=-s(2,1) ; s(5,1)=-s(1,1)
               wt(1)=.236926885056189 ; wt(2)=.478628670499366
               wt(3)=.568888888888889 ; wt(4)=wt(2) ; wt(5)=wt(1)
              case(6)
               s(1,1)=.932469514203152 ; s(2,1)=.661209386466265 
               s(3,1)=.238619186083197
               s(4,1)=-s(3,1) ; s(5,1)=-s(2,1) ; s(6,1)=-s(1,1)
               wt(1)=.171324492379170 ; wt(2)=.360761573048139 
               wt(3)=.467913934572691
               wt(4)=wt(3); wt(5)=wt(2) ; wt(6)=wt(1)
                    case default
                     print*,"wrong number of integrating points for a line"
             end select
            case('triangle') 
             select case(nip)
              case(1)   ! for triangles weights multiplied by .5
                s(1,1)=1./3.  ; s(1,2)=1./3.  ;  wt(1)= .5
              case(3)   
               s(1,1)=.5 ;  s(1,2)=.5 ;  s(2,1)=.5  
               s(2,2)=0.;  s(3,1)=0.  ;  s(3,2)=.5
               wt(1)=1./3.  ;  wt(2)=wt(1) ; wt(3)=wt(1)   ; wt = .5*wt
              case(6)
 s(1,1)=.816847572980459  ; s(1,2)=.091576213509771
 s(2,1)=s(1,2);  s(2,2)=s(1,1) ;  s(3,1)=s(1,2); s(3,2)=s(1,2)
 s(4,1)=.108103018168070 ;  s(4,2)=.445948490915965
 s(5,1)=s(4,2) ;   s(5,2)=s(4,1) ;  s(6,1)=s(4,2)  ; s(6,2)=s(4,2)
 wt(1)=.109951743655322 ;   wt(2)=wt(1)  ;   wt(3)=wt(1)
 wt(4)=.223381589678011 ;   wt(5)=wt(4)  ;   wt(6)=wt(4)    ; wt = .5*wt
              case(7)
 s(1,1)=1./3. ; s(1,2)=1./3.;s(2,1)=.797426985353087 ;s(2,2)=.101286507323456
 s(3,1)=s(2,2) ;  s(3,2)=s(2,1) ; s(4,1)=s(2,2) ;  s(4,2)=s(2,2)
 s(5,1)=.470142064105115 ;   s(5,2)=.059715871789770
 s(6,1)=s(5,2) ; s(6,2)=s(5,1);  s(7,1)=s(5,1);  s(7,2)=s(5,1)
 wt(1)=.225 ; wt(2)=.125939180544827 ;  wt(3)=wt(2);  wt(4)=wt(2)
 wt(5)=.132394152788506;  wt(6)=wt(5)      ;  wt(7)=wt(5)     ;wt = .5*wt
              case(12)
 s(1,1)=.873821971016996 ; s(1,2)=.063089014491502
 s(2,1)=s(1,2) ;  s(2,2)=s(1,1);  s(3,1)=s(1,2) ;  s(3,2)=s(1,2)
 s(4,1)=.501426509658179 ;  s(4,2)=.249286745170910
 s(5,1)=s(4,2); s(5,2)=s(4,1)   ;  s(6,1)=s(4,2) ;  s(6,2)=s(4,2)
 s(7,1)=.636502499121399 ;      s(7,2)=.310352451033785
 s(8,1)=s(7,1) ;  s(8,2)=.053145049844816 ;  s(9,1)=s(7,2) ; s(9,2)=s(7,1)
 s(10,1)=s(7,2) ; s(10,2)=s(8,2) ; s(11,1)=s(8,2);   s(11,2)=s(7,1)
 s(12,1)=s(8,2) ;  s(12,2)=s(7,2)
 wt(1)=.050844906370207 ; wt(2)=wt(1); wt(3)=wt(1)
 wt(4)=.116786275726379 ; wt(5)=wt(4); wt(6)=wt(4)
 wt(7)=.082851075618374 ; wt(8:12)=wt(7)           ; wt = .5*wt
              case(16)
 s(1,1)=1./3. ;  s(1,2)=1./3.  ;  s(2,1)=.658861384496478
 s(2,2)=.170569307751761 ; s(3,1)=s(2,2)   ;  s(3,2)=s(2,1)
 s(4,1)=s(2,2)  ; s(4,2)=s(2,2)
 s(5,1)=.898905543365938 ; s(5,2)=.050547228317031
 s(6,1)=s(5,2);  s(6,2)=s(5,1) ; s(7,1)=s(5,2)  ;  s(7,2)=s(5,2)
 s(8,1)=.081414823414554; s(8,2)=.459292588292723
 s(9,1)=s(8,2)  ;  s(9,2)=s(8,1);  s(10,1)=s(8,2) ;  s(10,2)=s(8,2)
 s(11,1)=.008394777409958; s(11,2)=.263112829634638
 s(12,1)=s(11,1)    ;  s(12,2)=.728492392955404
 s(13,1)=s(11,2) ;   s(13,2)=s(11,1)  ;  s(14,1)=s(11,2); s(14,2)=s(12,2)
 s(15,1)=s(12,2) ;  s(15,2)=s(11,1) ;  s(16,1)=s(12,2) ;  s(16,2)=s(11,2)
 wt(1)=.144315607677787 ; wt(2)=.103217370534718 ; wt(3)=wt(2); wt(4)=wt(2)
 wt(5)=.032458497623198 ; wt(6)=wt(5)   ;  wt(7)=wt(5)
 wt(8)=.095091634267284 ; wt(9)=wt(8)   ;  wt(10)=wt(8)
 wt(11)=.027230314174435 ; wt(12:16) = wt(11)  ;     wt = .5*wt
              case default
                  print*,"wrong number of integrating points for a triangle"
             end select
            case ('quadrilateral')
             select case (nip)
              case(1)
                s(1,1) = .0 ; wt(1) = 4.
              case(4)
                s(1,1)=-root3; s(1,2)= root3
                s(2,1)= root3; s(2,2)= root3
                s(3,1)=-root3; s(3,2)=-root3
                s(4,1)= root3; s(4,2)=-root3
                wt = 1.0
              case(9)
                s(1:7:3,1) = -r15; s(2:8:3,1) = .0
                s(3:9:3,1) =  r15; s(1:3,2)   = r15
                s(4:6,2)   =  .0 ; s(7:9,2)   =-r15
                     wt= v
              case default
                print*,"wrong number of integrating points for a quadrilateral"
            end select
            case('tetrahedron')    
             select case(nip)