vlib.f90

I[1].M.Smith所著的《有限元方法编程》第三版Fortran程序

      r0(1,2)=1.; r0(2,1)=-cg; r0(3,3)=cg; r0(2,3)=sg; r0(3,1)=sg
     end if
     do i=1,3; do j=1,3 
      x=r0(i,j)
      do k=0,9,3; t(i+k,j+k)=x; end do
     end do; end do
     do i=1,12
      sum=0.
       do j=1,12
        sum=sum+t(i,j)*global(j)
       end do
       local(i)=sum
     end do
   end select
   return
end subroutine glob_to_loc                
subroutine kvdet(kv,n,iw,det,ksc)
!
!      this subroutine forms the determinant of the global stiffness
!      matrix stored as a vector (upper triangle)
!
 implicit none
 real,intent(inout)::kv(:)
 integer,intent(in)::n,iw
 real,intent(out)::det
 integer,intent(out)::ksc
 real::const,em
 integer::iwp2,l,ll,j,ma,ii,k,i
 iwp2=iw+2
 l=n-1
 do j=1,l
   if((iw-n+j) <= 0.0)ma=iwp2
   ma=ma-1
   const=1./kv(j)
   ii=j
   do k=2,ma
     ii=ii+1
     if(kv((k-1)*n+j) /= 0.0)then
       em=kv((k-1)*n+j)*const
       ll=0
       do i=k,ma
         ll=ll+1
         kv((ll-1)*n+ii)=kv((ll-1)*n+ii)-em*kv((i-1)*n+j)
       end do
     end if
   end do
 end do
 det=1.
 ksc=0
 do j=1,n
   det=det*kv(j)
   if(kv(j) < 0.)ksc=ksc+1
 end do
 return
end subroutine kvdet            
subroutine hinge(coord,holdr,action,react,prop,iel,etype,gamma)
!
!      this subroutine forms the end forces and moments to be
!      applied to a member if a joint has gone plastic
!
implicit none
!interface
!  subroutine glob_to_loc(local,global,gamma,coord)
!    real,intent(in)::global(:),gamma,coord(:,:)
!    real,intent(out)::local(:)
!  end subroutine glob_to_loc
!  subroutine loc_to_glob(local,global,gamma,coord)
!    real,intent(in)::local(:),gamma,coord(:,:)
!    real,intent(out)::global(:)
!  end subroutine loc_to_glob
!end interface
real,intent(in):: holdr(:,:),coord(:,:),action(:),prop(:,:),gamma(:)
real,intent(out):: react(:)
integer,intent(in)::etype(:),iel
real::ell,x1,x2,y1,y2,z1,z2,csch,snch,bm1,bm2,bm3,bm4,bm5,&
      s1,s2,s3,s4,s5,mpy,mpz,mpx
real::global(12),local(12),total(12),gam
integer::ndim
ndim=ubound(coord,2)
bm1=0.; bm2=0.; bm3=0.; bm4=0.; bm5=0.
total(:)=holdr(:,iel)
select case(ndim)
  case(1)
    mpy=prop(2,etype(iel))
    ell=coord(2,1)-coord(1,1)
    s1=total(2)+action(2); s2=total(4)+action(4)
  if(abs(s1) > mpy)then
    if(s1 >  0.)bm1= mpy-s1
    if(s1 <= 0.)bm1=-mpy-s1
  end if
  if(abs(s2) > mpy)then
    if(s2 >  0.)bm2= mpy-s2
    if(s2 <= 0.)bm2=-mpy-s2
  end if
    react(1)= (bm1+bm2)/ell; react(2)=bm1
    react(3)=-react(1);      react(4)=bm2
  case(2)
    mpy=prop(3,etype(iel))  
    x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2)
    ell=sqrt((y2-y1)**2+(x2-x1)**2)
    csch=(x2-x1)/ell; snch=(y2-y1)/ell
    s1=total(3)+action(3); s2=total(6)+action(6)
  if(abs(s1) > mpy)then
    if(s1 >  0.)bm1= mpy-s1
    if(s1 <= 0.)bm1=-mpy-s1
  end if
  if(abs(s2) > mpy)then
    if(s2 >  0.)bm2= mpy-s2
    if(s2 <= 0.)bm2=-mpy-s2
  end if
    react(1)=-(bm1+bm2)*snch/ell; react(2)= (bm1+bm2)*csch/ell; react(3)=bm1
    react(4)=-react(1); react(5)=-react(2); react(6)=bm2
  case(3)
    gam=gamma(iel)
    mpy=prop(5,etype(iel))
    mpz=prop(6,etype(iel))
    mpx=prop(7,etype(iel))
    x1=coord(1,1); y1=coord(1,2); z1=coord(1,3)
    x2=coord(2,1); y2=coord(2,2); z2=coord(2,3)
    ell=sqrt((z2-z1)**2+(y2-y1)**2+(x2-x1)**2)
    global=total+action
    call glob_to_loc(local,global,gam,coord)
    global=0.0
    s1=local(5); s2=local(11)
    if(abs(s1) > mpy)then
      if(s1 >  0.)bm1= mpy-s1
      if(s1 <= 0.)bm1=-mpy-s1
    end if
      if(abs(s2) > mpy)then
      if(s2 >  0.)bm2= mpy-s2
      if(s2 <= 0.)bm2=-mpy-s2
    end if
      local( 3)=-(bm1+bm2)/ell
      local( 9)=-local(3)
      local( 5)= bm1
      local(11)= bm2
    s3=local(6); s4=local(12)
    if(abs(s3) > mpz)then
      if(s3 >  0.)bm1= mpz-s3
      if(s3 <= 0.)bm1=-mpz-s3
    end if
      if(abs(s4) > mpy)then
      if(s4 >  0.)bm2= mpz-s4
      if(s4 <= 0.)bm2=-mpz-s4
    end if
      local( 2)=(bm3+bm4)/ell
      local( 8)=-local(2)
      local( 6)= bm3
      local(12)= bm4
      s5=local(4)
      if(abs(s5) > mpx)then
        if(s5 >  0.)global(4)= mpx-s5
        if(s5 <= 0.)global(4)=-mpx-s5
      end if
      local(10)=-local(4)
      call loc_to_glob(local,react,gam,coord)
  end select
return
end subroutine hinge                           
 subroutine rigid_jointed(km,prop,gamma,etype,iel,coord) 
!
!      this subroutine forms the stiffness matrix of a
!      general beam/column element(1-d, 2-d or 3-d)
!
 implicit none
   real,intent(in)::gamma(:),coord(:,:),prop(:,:) 
   integer,intent(in)::etype(:),iel
   real,intent(out):: km(:,:)
   integer::ndim,i,j,k
   real::ell,x1,x2,y1,y2,z1,z2,c,s,e1,e2,e3,e4,pi,xl,yl,zl,cg,sg,den
   real::ea,ei,eiy,eiz,gj
   real::a1,a2,a3,a4,a5,a6,a7,a8,sum,gamrad,x
   real::t(12,12),tt(12,12),cc(12,12),r0(3,3)
   ndim=ubound(coord,2)
   select case(ndim)
    case(1)
     ei=prop(1,etype(iel)); ell=coord(2,1)-coord(1,1)
     km(1,1)=12.*ei/(ell*ell*ell) ;km(3,3)=km(1,1)
     km(1,2)=6.*ei/(ell*ell) ;  km(2,1)=km(1,2) ; km(1,4)=km(1,2)
     km(4,1)=km(1,4) ;  km(1,3)=-km(1,1) ; km(3,1)=km(1,3) ; km(3,4)=-km(1,2)
     km(4,3)=km(3,4) ; km(2,3)=km(3,4) ;  km(3,2)=km(2,3); km(2,2)=4.*ei/ell
     km(4,4)=km(2,2) ;km(2,4)=2.*ei/ell ; km(4,2)=km(2,4)
    case(2)
     ea=prop(1,etype(iel)); ei=prop(2,etype(iel))
     x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2)
     ell=sqrt((y2-y1)**2+(x2-x1)**2)
     c=(x2-x1)/ell; s=(y2-y1)/ell
     e1=ea/ell; e2=12.*ei/(ell*ell*ell); e3=ei/ell; e4=6.*ei/(ell*ell)
     km(1,1)=c*c*e1+s*s*e2; km(4,4)=km(1,1); km(1,2)=s*c*(e1-e2)
     km(2,1)=km(1,2); km(4,5)=km(1,2); km(5,4)=km(4,5); km(1,3)=-s*e4
     km(3,1)=km(1,3); km(1,6)=km(1,3); km(6,1)=km(1,6); km(3,4)=s*e4 
     km(4,3)=km(3,4); km(4,6)=km(3,4); km(6,4)=km(4,6); km(1,4)=-km(1,1) 
     km(4,1)=km(1,4); km(1,5)=s*c*(-e1+e2); km(5,1)=km(1,5); km(2,4)=km(1,5)
     km(4,2)=km(2,4); km(2,2)=s*s*e1+c*c*e2; km(5,5)=km(2,2); km(2,5)=-km(2,2)
     km(5,2)=km(2,5); km(2,3)=c*e4; km(3,2)=km(2,3); km(2,6)=km(2,3)
     km(6,2)=km(2,6); km(3,3)=4.*e3; km(6,6)=km(3,3); km(3,5)=-c*e4
     km(5,3)=km(3,5); km(5,6)=km(3,5); km(6,5)=km(5,6); km(3,6)=2.*e3
     km(6,3)=km(3,6)
    case(3)
     ea =prop(1,etype(iel)); eiy=prop(2,etype(iel))
     eiz=prop(3,etype(iel)); gj =prop(4,etype(iel))
     x1=coord(1,1); y1=coord(1,2); z1=coord(1,3)
     x2=coord(2,1); y2=coord(2,2); z2=coord(2,3)
     xl=x2-x1; yl=y2-y1; zl=z2-z1; ell=sqrt(xl*xl+yl*yl+zl*zl)
     km=0.0; t=0.0; tt=0.0
     a1=ea/ell; a2=12.*eiz/(ell*ell*ell); a3=12.*eiy/(ell*ell*ell)
     a4=6.*eiz/(ell*ell); a5=6.*eiy/(ell*ell); a6=4.*eiz/ell
     a7=4.*eiy/ell; a8=gj/ell
     km(1,1)=a1; km(7,7)=a1; km(1,7)=-a1; km(7,1)=-a1; km(2,2)=a2; km(8,8)=a2
     km(2,8)=-a2; km(8,2)=-a2; km(3,3)=a3; km(9,9)=a3; km(3,9)=-a3; km(9,3)=-a3
     km(4,4)=a8; km(10,10)=a8; km(4,10)=-a8; km(10,4)=-a8; km(5,5)=a7
     km(11,11)=a7; km(5,11)=.5*a7; km(11,5)=.5*a7; km(6,6)=a6; km(12,12)=a6
     km(6,12)=.5*a6; km(12,6)=.5*a6; km(2,6)=a4; km(6,2)=a4; km(2,12)=a4
     km(12,2)=a4; km(6,8)=-a4; km(8,6)=-a4; km(8,12)=-a4; km(12,8)=-a4
     km(5,9)=a5; km(9,5)=a5; km(9,11)=a5; km(11,9)=a5; km(3,5)=-a5
     km(5,3)=-a5; km(3,11)=-a5; km(11,3)=-a5
     pi=acos(-1.); gamrad=gamma(iel)*pi/180.; cg=cos(gamrad); sg=sin(gamrad)
     den=ell*sqrt(xl*xl+zl*zl)
     if(den /= 0.0)then
      r0(1,1)=xl/ell; r0(1,2)=yl/ell; r0(1,3)=zl/ell
      r0(2,1)=(-xl*yl*cg-ell*zl*sg)/den; r0(2,2)=den*cg/(ell*ell)
      r0(2,3)=(-yl*zl*cg+ell*xl*sg)/den; r0(3,1)=(xl*yl*sg-ell*zl*cg)/den
      r0(3,2)=-den*sg/(ell*ell); r0(3,3)=(yl*zl*sg+ell*xl*cg)/den
     else
      r0(1,1)=0.; r0(1,3)=0.; r0(2,2)=0.; r0(3,2)=0.
      r0(1,2)=1.; r0(2,1)=-cg; r0(3,3)=cg; r0(2,3)=sg; r0(3,1)=sg
     end if
     do i=1,3; do j=1,3 
      x=r0(i,j)
      do k=0,9,3; t(i+k,j+k)=x; tt(j+k,i+k)=x; end do
     end do; end do
     do i=1,12; do j=1,12
      sum=0.0
      do k=1,12; sum=sum+km(i,k)*t(k,j); end do
      cc(i,j)=sum
     end do; end do
     do i=1,12; do j=1,12
      sum=0.
      do k=1,12; sum=sum+tt(i,k)*cc(k,j); end do
      km(i,j)=sum
     end do; end do
  end select
 return
 end subroutine rigid_jointed                
 subroutine pin_jointed(km,ea,coord)
!
!      this subroutine forms the stiffness matrix of a
!      general rod element(1-d, 2-d or 3-d)
!
 implicit none
   real,intent(in)::ea,coord(:,:); real,intent(out):: km(:,:)
   integer::ndim ,i,j
   real::eaol,ell,cs,sn,x1,x2,y1,y2,z1,z2,a,b,c,d,e,f,xl,yl,zl
   ndim=ubound(coord,2)
   select case(ndim)
    case(1)
      ell=coord(2,1)-coord(1,1)
      eaol=ea/ell
      km(1,1)=eaol;  km(1,2)=-eaol
      km(2,1)=-eaol; km(2,2)=eaol  
    case(2)
      x1=coord(1,1); y1=coord(1,2)
      x2=coord(2,1); y2=coord(2,2)
      ell=sqrt((y2-y1)**2+(x2-x1)**2)
      cs=(x2-x1)/ell; sn=(y2-y1)/ell
      a=cs*cs; b=sn*sn; c=cs*sn
      km(1,1)=a; km(3,3)=a; km(1,3)=-a; km(3,1)=-a; km(2,2)=b; km(4,4)=b
      km(2,4)=-b; km(4,2)=-b; km(1,2)=c; km(2,1)=c; km(3,4)=c; km(4,3)=c
      km(1,4)=-c; km(4,1)=-c; km(2,3)=-c; km(3,2)=-c
      do i=1,4; do j=1,4; km(i,j)=km(i,j)*ea/ell; end do; end do
    case(3)
      x1=coord(1,1); y1=coord(1,2); z1=coord(1,3)
      x2=coord(2,1); y2=coord(2,2); z2=coord(2,3)
      xl=x2-x1; yl=y2-y1; zl=z2-z1
      ell=sqrt(xl*xl+yl*yl+zl*zl)
      xl=xl/ell; yl=yl/ell; zl=zl/ell
      a=xl*xl; b=yl*yl; c=zl*zl; d=xl*yl; e=yl*zl; f=zl*xl
      km(1,1)=a; km(4,4)=a; km(2,2)=b; km(5,5)=b; km(3,3)=c; km(6,6)=c
      km(1,2)=d; km(2,1)=d; km(4,5)=d; km(5,4)=d; km(2,3)=e; km(3,2)=e
      km(5,6)=e; km(6,5)=e; km(1,3)=f; km(3,1)=f; km(4,6)=f; km(6,4)=f
      km(1,4)=-a; km(4,1)=-a; km(2,5)=-b; km(5,2)=-b; km(3,6)=-c; km(6,3)=-c
      km(1,5)=-d; km(5,1)=-d; km(2,4)=-d; km(4,2)=-d; km(2,6)=-e; km(6,2)=-e
      km(3,5)=-e; km(5,3)=-e; km(1,6)=-f; km(6,1)=-f; km(3,4)=-f; km(4,3)=-f
      do i=1,6; do j=1,6; km(i,j)=km(i,j)*ea/ell; end do; end do
   end select
 return
 end subroutine pin_jointed       
 subroutine glob_to_axial(axial,global,coord)
!
!      this subroutine transforms the global end reactions
!      into an axial force for rod elements (2-d or 3-d)
!
 implicit none
 real,intent(in)::global(:),coord(:,:)
 real,intent(out)::axial
 real::add,ell
 integer::ndim,i
 ndim=ubound(coord,2)
 add = 0.
 do i = 1 , ndim
   add = add + (coord(2,i)-coord(1,i))**2
 end do
 ell=sqrt(add)
 axial = 0.
 do i = 1 , ndim
   axial = axial + (coord(2,i)-coord(1,i))/ell*global(ndim+i)
 end do
 return
 end subroutine glob_to_axial                       
end module vlib