vlib.f90

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

module vlib
contains                                         
 subroutine analy4(kma,coord,e,v)
!
! analytical stiffness of a 4-node quadrilateral element
! plane strain formulation
!
 implicit none
 real,intent(in):: coord(:,:),e,v
 real,intent(out):: kma(:,:)
 real:: e1,e2,g,x1,x2,x3,x4,y1,y2,y3,y4,a2,a2st3,alph,beta
 real:: c11,c21,c31,c41,c51,c61,c71,c81,c22,c32,c42,c52,c62,c72,c82
 real:: c33,c43,c53,c63,c73,c83,c44,c54,c64,c74,c84,c55,c65,c75,c85
 real:: c66,c76,c86,c77,c87,c88,f1,f2,s1,s2,s3,s4,t1,t2,t3,t4
 integer:: i,j
!
 e1=e*(1-v)/(1+v)/(1-2*v)
 e2=v*e1/(1-v)
 g=e/2/(1+v)      

 x1=coord(1,1)
 x2=coord(2,1)
 x3=coord(3,1)
 x4=coord(4,1)
 y1=coord(1,2)
 y2=coord(2,2)
 y3=coord(3,2)
 y4=coord(4,2)

 a2=(x4-x2)*(y3-y1)-(x3-x1)*(y4-y2)
 a2st3=3*a2*a2

 call groupa(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c11=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(1,1)=c11

 call groupa(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c33=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(3,3)=c33

 call groupa(x3,x4,x1,x2,y3,y4,y1,y2,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c55=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(5,5)=c55

 call groupa(x4,x1,x2,x3,y4,y1,y2,y3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c77=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,7)=c77

 call groupa(y4,y3,y2,y1,x4,x3,x2,x1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c88=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,8)=c88

 call groupa(y1,y4,y3,y2,x1,x4,x3,x2,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c22=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(2,2)=c22

 call groupa(y2,y1,y4,y3,x2,x1,x4,x3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c44=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(4,4)=c44

 call groupa(y3,y2,y1,y4,x3,x2,x1,x4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c66=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,6)=c66
 
 call groupb(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c21=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(2,1)=c21

 call groupb(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c43=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(4,3)=c43

 call groupb(x3,x4,x1,x2,y3,y4,y1,y2,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c65=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,5)=c65

 call groupb(x4,x1,x2,x3,y4,y1,y2,y3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c87=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,7)=c87

 call groupc(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c31=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(3,1)=c31

 call groupc(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c53=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(5,3)=c53

 call groupc(x3,x4,x1,x2,y3,y4,y1,y2,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c75=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,5)=c75

 call groupc(x4,x1,x2,x3,y4,y1,y2,y3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c71=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,1)=c71
 
 call groupc(y4,y3,y2,y1,x4,x3,x2,x1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c86=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,6)=c86

 call groupc(y1,y4,y3,y2,x1,x4,x3,x2,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c82=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,2)=c82

 call groupc(y2,y1,y4,y3,x2,x1,x4,x3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c42=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(4,2)=c42

 call groupc(y3,y2,y1,y4,x3,x2,x1,x4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c64=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,4)=c64

 call groupd(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c41=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(4,1)=c41

 call groupd(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c63=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,3)=c63

 call groupd(x3,x4,x1,x2,y3,y4,y1,y2,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x3,x4,x1,x2,y3,y4,y1,y2,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c85=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,5)=c85

 call groupd(x4,x1,x2,x3,y4,y1,y2,y3,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x4,x1,x2,x3,y4,y1,y2,y3,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c72=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,2)=c72

 call groupd(y1,y2,y3,y4,x1,x2,x3,x4,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c32=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(3,2)=c32

 call groupd(y2,y3,y4,y1,x2,x3,x4,x1,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c54=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(5,4)=c54

 call groupd(y3,y4,y1,y2,x3,x4,x1,x2,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x3,x4,x1,x2,y3,y4,y1,y2,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c76=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,6)=c76

 call groupd(y4,y1,y2,y3,x4,x1,x2,x3,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x4,x1,x2,x3,y4,y1,y2,y3,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c81=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,1)=c81

 call groupe(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c51=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(5,1)=c51

 call groupe(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c73=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,3)=c73

 call groupe(y2,y1,y4,y3,x2,x1,x4,x3,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c84=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,4)=c84

 call groupe(y3,y2,y1,y4,x3,x2,x1,x4,&
             s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 alph=a2*(s1*e1+s2*g)+f1*(s3*e1+s4*g)
 beta=a2*(t1*e1+t2*g)+f2*(t3*e1+t4*g)
 c62=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,2)=c62

 call groupf(x1,x2,x3,x4,y1,y2,y3,y4,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c61=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(6,1)=c61

 call groupf(x2,x3,x4,x1,y2,y3,y4,y1,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c83=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(8,3)=c83

 call groupf(y1,y2,y3,y4,x1,x2,x3,x4,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c52=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(5,2)=c52

 call groupf(y2,y3,y4,y1,x2,x3,x4,x1,&
             s1,s2,s3,s4,t1,t2,t3,t4)
 call f1f2(x2,x3,x4,x1,y2,y3,y4,y1,f1,f2)
 alph=a2*(s1*e2+s2*g)+f1*(s3*e2+s4*g)
 beta=a2*(t1*e2+t2*g)+f2*(t3*e2+t4*g)
 c74=(alph/(a2st3-f1**2)+beta/(a2st3-f2**2))*0.5
 kma(7,4)=c74
 do i=1,8; do j=i+1,8; kma(i,j)=kma(j,i); end do; end do

 return 
 end  subroutine analy4

 subroutine groupa(x1,x2,x3,x4,y1,y2,y3,y4,&
                   s1,s2,s3,s4,t1,t2,t3,t4,f1,f2)
 implicit none
 real,intent(in)::x1,x2,x3,x4,y1,y2,y3,y4
 real,intent(out)::s1,s2,s3,s4,t1,t2,t3,t4,f1,f2
 s1=2*(y4-y2)**2
 s2=2*(x4-x2)**2
 s3=-s1/2
 s4=-s2/2
 t1=(y2-y3)**2+(y3-y4)**2+(y4-y2)**2
 t2=(x2-x3)**2+(x3-x4)**2+(x4-x2)**2
 t3=(y4-y3)**2-(y3-y2)**2