text2.f90

fortran程序

 module new_library
contains                                  
 subroutine formblock_k(bigk,km,g,g_min)
 ! forms complete n*n matrix or submatrix
 implicit none
 real,intent(in)::km(:,:); integer,intent(in):: g(:),g_min
 real,intent(out)::bigk(:,:)
 integer::i,j,ndof;       ndof=ubound(km,1)
 do i = 1 , ndof
  if(g(i)/=0) then
   do j = 1 , ndof
    if(g(j)/=0) then
     bigk(g(i)-g_min+1,g(j)-g_min+1)=bigk(g(i)-g_min+1,g(j)-g_min+1) + km(i,j)
    end if
   end do
  end if
 end do
return
end subroutine formblock_k
subroutine sparin_gauss(kv,kdiag)
! Gaussian factorisation of a skyline matrix
implicit none
real,intent(out)::kv(:) ; integer,intent(in)::kdiag(:)
real::num,den,fac; integer::n,ii,i,j,k,l,kk,l1,l2,l3; n = ubound(kdiag,1)
 do j = 1 , n-1
    den = kv(kdiag(j))
    ii = 0 
    do i = j+1 , n
    ii = ii + 1 ;   l = kdiag(i) - ii
     if(l-kdiag(i-1)>.0) then
       num = kv(l)   ; fac = num/den ; kk = -1
       do k = i , n
          kk = kk + 1; l1=kdiag(i+kk)-kk;  l2=l1-ii; l3=kdiag(i+kk-1)
          if(l2-l3>.0) then
             kv(l1) = kv(l1) - fac*kv(l2)
          end if
       end do 
     end if
    end do
 end do
 return
end subroutine sparin_gauss
subroutine spabac_gauss(kv,loads,kdiag)
! Gaussian back-substitution on a skyline matrix
implicit none
real,intent(in)::kv(:); real,intent(inout)::loads(0:)
integer,intent(in)::kdiag(:)
real::num,den,fac,asum;integer::i,j,l,n,ii,jj,l1,l2; n=ubound(kdiag,1)
 do j = 1 , n-1
    den = kv(kdiag(j))  ;  ii = 0
    do i = j+1 , n
       ii = ii + 1      ; l = kdiag(i) - ii
       if(l-kdiag(i-1)>.0) then
          num = kv(l)   ; fac = num/den  ; loads(i)=loads(i)-fac*loads(j)
       end if
    end do
 end do
 loads(n) = loads(n)/kv(kdiag(n))
 do i = n-1 , 1 , -1
    jj = 0     ;  asum = .0
    do j = i+1 , n
       jj = jj + 1 ; l1 = kdiag(i+jj)-jj ;  l2 = kdiag(i+jj-1)
       if(l1-l2>.0) then
          asum = asum + kv(l1) * loads(j)
       end if
    end do
    loads(i) = (loads(i) - asum)/kv(kdiag(i))
 end do
 return
end subroutine spabac_gauss  
     subroutine readbf(nf,nodof,nrb)
  !----- blocks of input of nf data  -----------------------
      implicit none
      integer,intent(out)::nf(:,:); integer,intent(in)::nodof,nrb
      integer::nfd(nodof),i,j,l,n,if,it,is,nn;nn=ubound(nf,2)
      nf=1
      do l=1,nrb
       read(10,*)if,it,is,nfd
       do i=if,min(nn,it),is 
         nf(:,i) = nf(:,i) * nfd(:) 
       end do
      end do
      n=0
      do j = 1,nn  ;    do i= 1,nodof
          if(nf(i,j) /= 0) then
           n=n+1     ;   nf(i,j)=n
          end if
      end do       ;    end do
     return
     end subroutine readbf                
subroutine vmflow(stress,dsbar,vmfl)
! Forms the von Mises flow vector
implicit none
real,intent(in)::stress(:),dsbar; real,intent(out)::vmfl(:)
real::sigm; sigm=(stress(1)+stress(2)+stress(4))/3.
vmfl(1)=stress(1)-sigm;vmfl(2)=stress(2)-sigm
vmfl(3)=stress(3)*2. ; vmfl(4)=stress(4)-sigm
vmfl = vmfl*1.5/dsbar
return
end subroutine vmflow
subroutine formaa(flow,rmat,modee)
! Modification to dee matrix
implicit none
real,intent(in)::flow(:),rmat(:,:); real,intent(out)::modee(:,:)
real::flowt(1,4),flowa(4,1)   ; flowt(1,:) = flow ; flowa(:,1) = flow
modee = matmul(matmul(matmul(rmat,flowa),flowt),rmat)
return
end subroutine formaa
subroutine fmrmat(vmfl,dsbar,dlam,dee,temp,rmat)
! Forms the r matrix
implicit none
real,intent(in)::vmfl(:),dsbar,dlam,dee(:,:),temp(:,:)
real,intent(out)::rmat(:,:); real::acat(4,4),acatc(4,4),qmat(4,4)
integer::i,j  ; real:: con
do i=1,4;do j=1,4;acat(i,j)=vmfl(i)*vmfl(j);end do;end do
acat = (temp-acat)/dsbar; acatc = matmul(dee,acat); qmat = acatc*dlam
do i=1,4; qmat(i,i)=qmat(i,i)+1.; end do
 do i=1,4
  con = qmat(i,i); qmat(i,i) = 1.
  qmat(i,:) = qmat(i,:)/con
  do j=1,4
   if(j/=i) then
    con = qmat(j,i); qmat(j,i)=0.0
    qmat(j,:) = qmat(j,:) - qmat(i,:)*con
   end if
  end do
 end do
 rmat = matmul(qmat,dee)
return
end subroutine fmrmat
subroutine fmacat(vmfl,temp,acat)
! Intermediate step
implicit none
real,intent(in)::vmfl(:),temp(:,:);real,intent(out)::acat(:,:);integer::i,j
do i=1,4; do j=1,4; acat(i,j)=vmfl(i)*vmfl(j); end do; end do
acat = temp - acat
return 
end subroutine fmacat
 subroutine formke(km,kp,c,ke,theta)
 ! ke matrix for incremental Biot
 implicit none
    real,intent(in)::km(:,:),kp(:,:),c(:,:),theta;real,intent(out)::ke(:,:)
    integer::ndof; ndof = ubound(km,1)
    ke(:ndof, :ndof)=km ;    ke(:ndof,ndof+1:)=c
    ke(ndof+1:,:ndof)=transpose(c) ; ke(ndof+1:,ndof+1:) = -theta * kp
  return
 end subroutine formke

  subroutine bmat_nonaxi(bee,radius,coord,deriv,fun,iflag,lth)
!
!      this subroutine forms the strain-displacement matrix for
!      axisymmetric solids subjected to non-axisymmetric loading
!
 implicit none
  real,intent(in)::deriv(:,:),fun(:),coord(:,:)
  real,intent(out):: bee(:,:),radius; integer,intent(in)::iflag,lth
  integer::nod,k,l,m,n; nod = ubound(deriv , 2 ) ; bee = .0
  radius = sum(fun * coord(:,1))
    do  m=1,nod
      n=3*m      ;    k=n-1      ;    l=k-1
      bee(1,l)=deriv(1,m);   bee(2,k)=deriv(2,m);   bee(3,l)=fun(m)/radius
      bee(3,n)=iflag*lth*bee(3,l) ; bee(4,l)=deriv(2,m); bee(4,k)=deriv(1,m)
      bee(5,k)=-iflag*lth*fun(m)/radius ;   bee(5,n)=deriv(2,m)
      bee(6,l)=bee(5,k) ; bee(6,n)=deriv(1,m)-fun(m)/radius
    end do
   return
  end subroutine bmat_nonaxi

subroutine ecmat(ecm,fun,ndof,nodof)
implicit none
 real,intent(in)::fun(:); real,intent(out)::ecm(:,:)
 integer,intent(in):: nodof,ndof
 integer:: nod,i,j; real::nt(ndof,nodof),tn(nodof,ndof)
 nod = ndof/nodof; nt = .0; tn = .0
  do i = 1 , nod ; do j = 1 , nodof
     nt((i-1)*nodof+j,j) = fun(i); tn(j,(i-1)*nodof+j) = fun(i)
  end do; end do
  ecm = matmul( nt , tn )
 return
end subroutine ecmat
subroutine vmpl(e,v,stress,pl)
! plastic matrix for a von-Mises material
implicit none
real,intent(in)::e,v,stress(:); real,intent(out)::pl(:,:)
real::sx,sy,txy,sz,dsbar,ee,term(4); integer::i,j,nst; nst = ubound(stress,1)
  sx=stress(1); sy=stress(2); txy=stress(3); sz=stress(4)
  dsbar=sqrt((sx-sy)**2+(sy-sz)**2+(sz-sx)**2+6.*txy**2)/sqrt(2.)
  ee=1.5*e/((1.+v)*dsbar*dsbar)
  term(1)=(2.*sx-sy-sz)/3.; term(2)=(2.*sy-sz-sx)/3.
  term(3)=txy             ; term(4)=(2.*sz-sx-sy)/3.
  do i=1,nst;do j=1,nst;pl(i,j)=term(i)*term(j)*ee;pl(j,i)=pl(i,j);end do;end do
return
end subroutine vmpl 
 subroutine formupv(ke,c11,c12,c21,c23,c32)
!
!      this subroutine forms the unsymmetrical stiffness matrix
!      for the u-v-p version of the Navier-Stokes equations
!
 implicit none
  real,intent(in):: c11(:,:),c21(:,:),c23(:,:),c32(:,:),c12(:,:)
  real,intent(out):: ke(:,:)  ; integer::nod,nodf,ntot
   nod = ubound(c11,1); nodf = ubound(c21,1) ;ntot=nod+nodf+nod
   ke(1:nod,1:nod)=c11  ;  ke(1:nod,nod+1:nod+nodf)=c12
   ke(nod+1:nod+nodf,1:nod)=c21  ; ke(nod+1:nod+nodf,nod+nodf+1:ntot)=c23
   ke(nod+nodf+1:ntot,nod+1:nod+nodf)=c32 
   ke(nod+nodf+1:ntot,nod+nodf+1:ntot)=c11
  return
 end subroutine formupv 
 subroutine gauss_band(pb,work)
!
!      this subroutine performs gaussian reduction of an
!      unsymmetric banded matrix 'pb' .  array 'work'
!      used as working space.
!
implicit none
  real,intent(in out):: pb(:,:),work(:,:)   
  integer::n,iwp1,iq,iqp,iwp11,i,j,k,l,ip,k1  ;   real::s
  n=ubound(pb,1); iwp1=(ubound(pb,2)-1)/2+1
  iq=2*iwp1-1  ;  iqp=iwp1 ;  iwp11=iwp1-1
   do 1 i=1,iwp11
     do 1 j=1,iq
      if(j>=iwp1+i)goto 2
       pb(i,j)=pb(i,j+iwp1-i)  
      goto 1
   2   pb(i,j)=0.
      pb(n-i+1,j)=0.
   1 continue
      do 3 k=1,n
      l=k+iwp1-1
      if(l>n)l=n  ;  ip=0  ;    s=1.e-10
     do i=k,l
      if(abs(pb(i,1))<=s) cycle
      s=abs(pb(i,1))    ;   ip=i
     end do
      if(ip==0)goto 5   ;    if(k==n)goto 11
      work(iwp1,k)=ip  ;     iqp=iqp-1 ;   j=iwp1+ip-k
      if(iqp<j)iqp=j    ;   if(j==iwp1)goto 6
    do  j=1,iqp
      s=pb(k,j)
      pb(k,j)=pb(ip,j)
      pb(ip,j)=s
    end do
    6 k1=k+1
      do 8 i=k1,l
      s=pb(i,1)/pb(k,1)
      do 9 j=2,iq
      if(j>iqp)goto 10
      pb(i,j-1)=pb(i,j)-s*pb(k,j)
      goto 9
   10 pb(i,j-1)=pb(i,j)
    9 continue
      pb(i,iq)=0.      ;    work(i-k,k)=s
    8 continue
    3 continue
    5 write(6,'(" singular")')
   11 return
 end subroutine gauss_band
 subroutine solve_band(pb,copy,ans)
!
!      this subroutine performs the gaussian back-substitution
!      on the reduced unsymmetric band matrix 'pb'.
!
 implicit none
   real,intent(in):: pb(:,:),copy(:,:); real,intent(out)::ans(0:)
   integer::iwp1,n,n1,i,iv,l,iq,iv1     ;  real ::s 
   iwp1= (ubound(pb,2)-1)/2+1; n=ubound(pb,1); iq=2*iwp1-1  ;   n1=n-1
   do  iv=1,n1
     i=int(copy(iwp1,iv)+.5)
     if(i/=iv) then
      s=ans(iv) ;   ans(iv)=ans(i) ;    ans(i)=s
     end if
     l=iv+iwp1-1
     if(l>n)l=n   ;    iv1=iv+1
     do  i=iv1,l
      ans(i)=ans(i)-copy(i-iv,iv)*ans(iv)
     end do
   end do
   ans(n)=ans(n)/pb(n,1)
      iv=n-1
    6 s=ans(iv)
      l=iq
      if(iv+l-1>n)l=n-iv+1
     do  i=2,l
      s=s-pb(iv,i)*ans(iv+i-1)   ;      ans(iv)=s/pb(iv,1)
     end do
      iv=iv-1
      if(iv/=0)goto 6
     return
   end subroutine solve_band                              
 subroutine formku(ku,km,g)
!
!      this subroutine assembles element matrices into symmetrical
!      global matrix(stored as an upper rectangle)
!
      real,intent(in):: km(:,:) ; real,intent(out)::ku(:,:)
      integer,intent(in):: g(:) ;integer::i,j,icd,ndof; ndof=ubound(km,1)
      do  i=1,ndof
      if(g(i)/=0) then
        do  j=1,ndof
          if(g(j)/=0) then
           icd=g(j)-g(i)+1
           if(icd>=1) ku(g(i),icd)=ku(g(i),icd)+km(i,j)
          end if
        end do
      end if
     end do
   return
 end subroutine formku
 subroutine bandred(a,d,e,e2)
!
!      this subroutine transforms a real symmetric band matrix a,
!      of order n and band width iw,to tridiagonal form by an appropriate
!      sequence of jacobi rotations. during the transformation the
!      property of the band matrix is maintained. the method yields
!      a tridiagonal matrix, the diagonal elements of which are in
!      d(n) and off-diagonal elements in e(n).
!
 implicit none  
    real,intent(in out)::a(:,:); real,intent(out)::d(0:),e(0:),e2(0:)    
      integer::  iw, n2, n, k, maxr, irr, ir, kr, j, jm, iugl, j2,    &
                 l, jl, maxl, i
      real ::    g, b, s, c, c2, s2, cs, u, u1
       n=ubound(a,1) ; iw = ubound(a,2)-1
      n2 = n - 2
      if (n2>=1) then
      do 160 k=1,n2
         maxr = iw      ;     if (n-k<iw) maxr = n - k
         do 140 irr=2,maxr
            ir = 2 + maxr - irr  ;     kr = k + ir
            do 120 j=kr,n,iw
               if (j==kr) go to 20  ;  if (g==0.0) go to 140
               jm = j - iw   ; b = -a(jm-1,iw+1)/g   ; iugl = j - iw
               go to 40
       20      if (a(k,ir+1)==0.0) go to 140
               b = -a(k,ir)/a(k,ir+1) ;  iugl = k