text2.f90

fortran程序

      case(3) ! three dimensional cases
       xi=points(i,1); eta=points(i,2); zeta=points(i,3)
       etam=1.-eta ;  xim=1.-xi  ;  zetam=1.-zeta
       etap=eta+1. ;  xip=xi+1.   ;  zetap=zeta+1.
       select case(nod)
        case(4)
         fun(1)=xi   ;   fun(2)= eta ;  fun(3)=zeta 
         fun(4)=1.-fun(1)-fun(2)-fun(3)
        case(8)
         fun=(/.125*xim*etam*zetam,.125*xim*etam*zetap,.125*xip*etam*zetap,&
               .125*xip*etam*zetam,.125*xim*etap*zetam,.125*xim*etap*zetap,&
               .125*xip*etap*zetap,.125*xip*etap*zetam/)
        case(14) !type 6 element
    x = points(i,1);  y = points(i,2);  z = points(i,3)
  fun(1)=((x*y+x*z+2.*x+y*z+2.*y+2.*z+2.)*(x-1.)*(y-1.)*(z-1.))/8.
  fun(2)=((x*y-x*z-2.*x+y*z+2.*y-2.*z-2.)*(x-1.)*(y+1.)*(z-1.))/8.
  fun(3)=((x*y+x*z+2.*x-y*z-2.*y-2.*z-2.)*(x+1.)*(y-1.)*(z-1.))/8.
  fun(4)=((x*y-x*z-2.*x-y*z-2.*y+2.*z+2.)*(x+1.)*(y+1.)*(z-1.))/8.
  fun(5)=-((x*y-x*z+2.*x-y*z+2.*y-2.*z+2.)*(x-1.)*(y-1.)*(z+1.))/8.
  fun(6)=-((x*y+x*z-2.*x-y*z+2.*y+2.*z-2.)*(x-1.)*(y+1.)*(z+1.))/8.
  fun(7)=-((x*y-x*z+2.*x+y*z-2.*y+2.*z-2.)*(x+1.)*(y-1.)*(z+1.))/8.
  fun(8)=-((x*y+x*z-2.*x+y*z-2.*y-2.*z+2.)*(x+1.)*(y+1.)*(z+1.))/8.
  fun(9)=-((x+1.)*(x-1.)*(y+1.)*(y-1.)*(z-1.))/2.
  fun(10)=((x+1.)*(x-1.)*(y+1.)*(y-1.)*(z+1.))/2.
  fun(11)=-((x+1.)*(x-1.)*(y-1.)*(z+1.)*(z-1.))/2.
  fun(12)=((x+1.)*(x-1.)*(y+1.)*(z+1.)*(z-1.))/2.
  fun(13)=-((x-1.)*(y+1.)*(y-1.)*(z+1.)*(z-1.))/2.
  fun(14)=((x+1.)*(y+1.)*(y-1.)*(z+1.)*(z-1.))/2.      
        case(20)
           xii=(/-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0/)
           etai=(/-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1/)
           zetai=(/-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1/)
           do l=1,20
            xi0=xi*xii(l); eta0=eta*etai(l); zeta0=zeta*zetai(l)
            if(l==4.or.l==8.or.l==16.or.l==20) then
              fun(l)=.25*(1.-xi*xi)*(1.+eta0)*(1.+zeta0)
            else if(l>=9.and.l<=12)then
              fun(l)=.25*(1.+xi0)*(1.-eta*eta)*(1.+zeta0)
            else if(l==2.or.l==6.or.l==14.or.l==18) then
              fun(l)=.25*(1.+xi0)*(1.+eta0)*(1.-zeta*zeta)
            else
              fun(l)=.125*(1.+xi0)*(1.+eta0)*(1.+zeta0)*(xi0+eta0+zeta0-2)
            end if
           end do
          case default
           print*,"wrong number of nodes in shape_fun"
        end select
      case default
        print*,"wrong number of dimensions in shape_fun"  
    end select
  return
 end subroutine shape_fun 
subroutine rod_km(km,ea,length)
  implicit none
  real,intent(in):: ea,length; real,intent(out)::km(:,:)
  km(1,1)=1. ; km(2,2)=1.; km(1,2)=-1. ; km(2,1)=-1.; km=km*ea/length
end subroutine rod_km 
subroutine axialkm(km,e,area,length)
  implicit none
  real,intent(in):: e,area,length; real,intent(out)::km(:,:)
  km(1,1)=1. ; km(2,2)=1.; km(1,2)=-1. ; km(2,1)=-1.; km=km*e*area/length
end subroutine axialkm    
  subroutine mocouf(phi,c,sigm,dsbar,theta,f)
!
!      this subroutine calculates the value of the yield function
!      for a Mohr-Coulomb material (phi in degrees)
!
 implicit none
  real,intent(in)::phi,c,sigm,dsbar,theta   ; real,intent(out)::f
  real::phir,snph,csph,csth,snth
  phir=phi*4.*atan(1.)/180.
  snph=sin(phir) ;  csph=cos(phir) ; csth=cos(theta);  snth=sin(theta)
  f=snph*sigm+dsbar*(csth/sqrt(3.)-snth*snph/3.)-c*csph
  return
  end  subroutine mocouf
  subroutine mocouq(psi,dsbar,theta,dq1,dq2,dq3)
!
!      this subroutine forms the derivatives of a Mohr-Coulomb
!      potential function with respect to the three invariants
!      psi in degrees
!
 implicit none
  real,intent(in)::psi,dsbar,theta; real,intent(out)::dq1,dq2,dq3
  real::psir,snth,snps,sq3,c1,csth,cs3th,tn3th,tnth     
  psir=psi*4.*atan(1.)/180. ;  snth=sin(theta) ;   snps=sin(psir)
  sq3=sqrt(3.)  ; dq1=snps
      if(abs(snth).gt..49)then
        c1=1.
        if(snth.lt.0.)c1=-1.
        dq2=(sq3*.5-c1*snps*.5/sq3)*sq3*.5/dsbar ;   dq3=0.
      else
        csth=cos(theta); cs3th=cos(3.*theta);tn3th=tan(3.*theta); tnth=snth/csth
        dq2=sq3*csth/dsbar*((1.+tnth*tn3th)+snps*(tn3th-tnth)/sq3)*.5
        dq3=1.5*(sq3*snth+snps*csth)/(cs3th*dsbar*dsbar)
      end if
  return
  end subroutine mocouq
  subroutine mocopl(phi,psi,e,v,stress,pl)
!
!      this subroutine forms the plastic stress/strain matrix
!      for a mohr-coulomb material  (phi,psi in degrees)
!
  implicit none
   real,intent(in)::stress(:),phi,psi,e,v  ;real,intent(out)::pl(:,:)
   integer::i,j;  real::row(4),col(4),sx,sy,txy,sz,pi,phir,psir,&
                            dx,dy,dz,d2,d3,th,snth,sig,rph,rps,cps,snps,sq3,&
                            cc,cph ,alp,ca,sa,dd ,snph,ee,s1,s2
   sx=stress(1);  sy=stress(2);  txy=stress(3) ;  sz=stress(4)
   pi=4.*atan(1.) ; phir=phi*pi/180.; psir=psi*pi/180.;  snph=sin(phir)
   snps=sin(psir)     ;  sq3=sqrt(3.) ;   cc=1.-2.*v
   dx=(2.*sx-sy-sz)/3.    ;   dy=(2.*sy-sz-sx)/3. ;  dz=(2.*sz-sx-sy)/3.
   d2=sqrt(-dx*dy-dy*dz-dz*dx+txy*txy) ;  d3=dx*dy*dz-dz*txy*txy
   th=-3.*sq3*d3/(2.*d2**3)
   if(th.gt.1.)th=1.   ;  if(th.lt.-1.)th=-1.
   th=asin(th)/3.   ;    snth=sin(th)
   if(abs(snth).gt..49)then
      sig=-1.
      if(snth.lt.0.)sig=1.
      rph=snph*(1.+v)/3.;  rps=snps*(1.+v)/3. ; cps=.25*sq3/d2*(1.+sig*snps/3.)
      cph=.25*sq3/d2*(1.+sig*snph/3.)
      col(1)=rph+cph*((1.-v)*dx+v*(dy+dz));col(2)=rph+cph*((1.-v)*dy+v*(dz+dx))
      col(3)=cph*cc*txy   ; col(4)=rph+cph*((1.-v)*dz+v*(dx+dy))
      row(1)=rps+cps*((1.-v)*dx+v*(dy+dz));row(2)=rps+cps*((1.-v)*dy+v*(dz+dx))
      row(3)=cps*cc*txy ;  row(4)=rps+cps*((1.-v)*dz+v*(dx+dy))
      ee=e/((1.+v)*cc*(rph*snps+2.*cph*cps*d2*d2*cc))
  else
      alp=atan(abs((sx-sy)/(2.*txy))) ; ca=cos(alp);  sa=sin(alp)
      dd=cc*sa ;  s1=1.  ;  s2=1.
      if((sx-sy).lt..0)s1=-1.
      if(txy.lt..0)s2=-1.
      col(1)=snph+s1*dd ; col(2)=snph-s1*dd;col(3)=s2*cc*ca; col(4)=2.*v*snph
      row(1)=snps+s1*dd ;row(2)=snps-s1*dd; row(3)=s2*cc*ca ; row(4)=2.*v*snps
      ee=e/(2.*(1.+v)*cc*(snph*snps+cc))
  end if
  do  i=1,4; do j=1,4 ; pl(i,j)=ee*row(i)*col(j) ; end do ; end do
  return
 end subroutine mocopl
subroutine formkv(bk,km,g,n)
!global stiffness matrix stored as a vector (upper triangle)
implicit none
real,intent(in)::km(:,:);real,intent(out)::bk(:)
integer,intent(in)::g(:),n
integer::idof,i,j,icd,ival
idof=size(km,1)
     do i=1,idof
        if(g(i)/=0) then
           do j=1,idof
              if(g(j)/=0) then
                 icd=g(j)-g(i)+1
                 if(icd-1>=0) then
                    ival=n*(icd-1)+g(i)      
                    bk(ival)=bk(ival)+km(i,j)
                 end if
               end if
            end do
         end if
     end do
return
end subroutine formkv
subroutine fsparv(bk,km,g,kdiag)
! assembly of element matrices into skyline global matrix
implicit none
real,intent(in)::km(:,:); integer,intent(in)::g(:),kdiag(:)
real,intent(out)::bk(:) ;  integer::i,idof,k,j,iw,ival
 idof=ubound(g,1)
   do i=1,idof
      k=g(i)
      if(k/=0) then
         do j=1,idof
            if(g(j)/=0) then
               iw=k-g(j)
               if(iw>=0) then
                   ival=kdiag(k)-iw
                   bk(ival)=bk(ival)+km(i,j) 
                end if
            end if
         end do
      end if
   end do
 return
end subroutine fsparv
subroutine banred(bk,n)
! gaussian reduction on a vector stored as an upper triangle
implicit none
real,intent(in out)::bk(:);integer,intent(in)::n
integer::i,il1,kbl,j,ij,nkb,m,ni,nj,iw ; real::sum
 iw = ubound(bk,1)/n-1
       do i=2,n
          il1=i-1;kbl=il1+iw+1
          if(kbl-n>0)kbl=n
          do j=i,kbl
             ij=(j-i)*n+i;sum=bk(ij);nkb=j-iw
             if(nkb<=0)nkb=1
             if(nkb-il1<=0)then
                do m=nkb,il1
                   ni=(i-m)*n+m ; nj=(j-m)*n+m
                   sum=sum-bk(ni)*bk(nj)/bk(m) 
                end do
             end if
             bk(ij)=sum
           end do
       end do
return
end subroutine banred   
subroutine bacsub(bk,loads)
! performs the complete gaussian backsubstitution
implicit none
real,intent(in)::bk(:);real,intent(in out)::loads(0:)
integer::nkb,k,i,jn,jj,i1,n,iw;real::sum
n = ubound(loads,1); iw = ubound(bk,1)/n - 1
loads(1)=loads(1)/bk(1)
   do i=2,n
      sum=loads(i);i1=i-1 ; nkb=i-iw
      if(nkb<=0)nkb=1
      do k=nkb,i1
         jn=(i-k)*n+k;sum=sum-bk(jn)*loads(k)
      end do
      loads(i)=sum/bk(i)
   end do
   do jj=2,n
      i=n-jj+1;sum=.0;i1=i+1;nkb=i+iw
      if(nkb-n>0)nkb=n
      do k=i1,nkb
           jn=(k-i)*n+i  ; sum=sum+bk(jn)*loads(k)
      end do
      loads(i)=loads(i)-sum/bk(i)
   end do
return
end subroutine bacsub                                                       
subroutine print_vector(vec,n,channel)
! prints out first n terms of a vector in e format to channel x 
implicit none
real,intent(in out)::vec(0:);integer,intent(in)::n,channel
integer::i
write(channel,'(1x,6g12.4)')(vec(i),i=1,n)
return
end subroutine print_vector 
subroutine print_array(array,m,n,channel)
! prints out first m rows and n columns of 'array' to channel x
implicit none
real,intent(in out)::array(:,:);integer,intent(in)::m,n,channel
integer::i,j
   do i=1,m
      write(channel,'(1x,6g12.4)')(array(i,j),j=1,n)
   end do
end subroutine print_array  
subroutine cholin(kb)
! Choleski reduction on kb(l,iw+1) stored as a lower triangle
implicit none
real,intent(in out)::kb(:,:);integer::i,j,k,l,ia,ib,n,iw;real::x
n=ubound(kb,1); iw=ubound(kb,2)-1
   do i=1,n
      x=.0
      do j=1,iw; x=x+kb(i,j)**2; end do
      kb(i,iw+1)=sqrt(kb(i,iw+1)-x)
      do k=1,iw
         x=.0
         if(i+k<=n) then
           if(k/=iw) then
             do l=iw-k,1,-1
                x=x+kb(i+k,l)*kb(i,l+k)
             end do
           end if
           ia=i+k; ib=iw-k+1
           kb(ia,ib)=(kb(ia,ib)-x)/kb(i,iw+1)
         end if
      end do
   end do
 return
end subroutine cholin  
subroutine chobac(kb,loads)
!Choleski back-substitution
implicit none
real,intent(in)::kb(:,:);real,intent(in out)::loads(0:)
integer::iw,n,i,j,k,l,m; real::x
n=size(kb,1); iw=size(kb,2)-1
loads(1)=loads(1)/kb(1,iw+1)
    do i=2,n
       x=.0;k=1
       if(i<=iw+1)k=iw-i+2
       do j=k,iw; x=x+kb(i,j)*loads(i+j-iw-1); end do
       loads(i)=(loads(i)-x)/kb(i,iw+1)
    end do
    loads(n)=loads(n)/kb(n,iw+1)
    do i=n-1,1,-1
       x=0.0; l=i+iw
       if(i>n-iw)l=n;  m=i+1
       do j=m,l; x=x+kb(j,iw+i-j+1)*loads(j); end do
       loads(i)=(loads(i)-x)/kb(i,iw+1)
    end do
  return
end subroutine chobac
subroutine sparin(a,kdiag)
! Choleski factorisation of variable bandwidth matrix a
! stored as a vector and overwritten
implicit none
real,intent(in out)::a(:);integer,intent(in)::kdiag(:)
integer::n,i,ki,l,kj,j,ll,m,k; real::x
 n=ubound(kdiag,1)  ; a(1)=sqrt(a(1))
 do i=2,n
    ki=kdiag(i)-i;  l=kdiag(i-1)-ki+1
    do j=l,i
       x=a(ki+j);  kj=kdiag(j)-j
       if(j/=1) then
          ll=kdiag(j-1)-kj+1; ll=max0(l,ll)
          if(ll/=j) then
              m=j-1
              do k=ll,m ; x=x-a(ki+k)*a(kj+k) ; end do
          end if
       end if
       a(ki+j)=x/a(kj+j)
    end do
    a(ki+i)=sqrt(x)
 end do
 return
end subroutine sparin
subroutine spabac(a,b,kdiag)
! Choleski forward and backward substitution combined
! variable bandwidth factorised matrix a stored as a vector 
implicit none
real,intent(in)::a(:);real,intent(in out)::b(0:);integer,intent(in)::kdiag(:)
integer::n,i,ki,l,m,j,it,k; real::x
n=ubound(kdiag,1)
 b(1)=b(1)/a(1)
  do i=2,n
     ki=kdiag(i)-i;  l=kdiag(i-1)-ki+1 ; x=b(i)
     if(l/=i) then
        m=i-1
        do j=l,m ; x=x-a(ki+j)*b(j); end do
     end if
     b(i)=x/a(ki+i)
  end do
  do it=2,n
     i=n+2-it; ki=kdiag(i)-i; x=b(i)/a(ki+i); b(i)=x; l=kdiag(i-1)-ki+1
     if(l/=i) then