vlib.f90

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

 t4=(x4-x3)**2-(x3-x2)**2
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-2*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-2*(x3*y1-x1*y3)
 return
 end subroutine groupa

 subroutine groupb(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*(x2-x4)*(y4-y2)
 s2=s1
 s3=-s1/2
 s4=s3
 t1=x2*(y4-2*y2+y3)+x3*(y2-2*y3+y4)+x4*(y2-2*y4+y3)
 t2=t1
 t3=x2*(y2-y3)+x3*(y4-y2)+x4*(y3-y4)
 t4=t3
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-2*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-2*(x3*y1-x1*y3)
 return
 end subroutine groupb

 subroutine groupc(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=(y4-y2)*(2*y1-y3-y4)
 s2=(x4-x2)*(2*x1-x3-x4)
 s3=(y4-y2)*(y4-y1)
 s4=(x4-x2)*(x4-x1)
 t1=(y3-y1)*(2*y2-y3-y4)
 t2=(x3-x1)*(2*x2-x3-x4)
 t3=(y3-y1)*(y3-y2)
 t4=(x3-x1)*(x3-x2)
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-2*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-2*(x3*y1-x1*y3)
 return
 end subroutine groupc

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

 subroutine groupe(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=-(y4-y2)**2
 s2=-(x4-x2)**2
 s3=0
 s4=0
 t1=(y3+y1)*(y4+y2)-2*(y4-y2)**2-2*(y1*y3+y2*y4)
 t2=(x3+x1)*(x4+x2)-2*(x4-x2)**2-2*(x1*x3+x2*x4)
 t3=(y4-y2)*(y1-y2+y3-y4)
 t4=(x4-x2)*(x1-x2+x3-x4)
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-2*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-2*(x3*y1-x1*y3)
 return
 end subroutine groupe

 subroutine groupf(x1,x2,x3,x4,y1,y2,y3,y4,&
                   s1,s2,s3,s4,t1,t2,t3,t4)
 implicit none
 real,intent(in)::x1,x2,x3,x4,y1,y2,y3,y4
 real,intent(out)::s1,s2,s3,s4,t1,t2,t3,t4
 s1=(x4-x2)*(y4-y2)
 s2=s1
 s3=0
 s4=0
 t1=(x4-x2)*(y4-y2)+(x2-x1)*(y2-y3)+(x4-x1)*(y4-y3)
 t2=(y4-y2)*(x4-x2)+(y2-y1)*(x2-x3)+(y4-y1)*(x4-x3)
 t3=(x2-x1)*(y3-y2)+(x4-x1)*(y4-y3)
 t4=(y2-y1)*(x3-x2)+(y4-y1)*(x4-x3)
 return
 end subroutine groupf

 subroutine f1f2(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2)
 implicit none
 real,intent(in)::x1,x2,x3,x4,y1,y2,y3,y4
 real,intent(out)::f1,f2
 f1=(x1+x3)*(y4-y2)-(y1+y3)*(x4-x2)-2*(x2*y4-x4*y2)
 f2=(y2+y4)*(x3-x1)-(x2+x4)*(y3-y1)-2*(x3*y1-x1*y3)
 return
 end subroutine f1f2
subroutine bee4(coord,points,i,det,bee)

!  analytical version of the bee matrix for a 4-node quad      

 implicit none
 real,intent(in):: coord(:,:),points(:,:)
 real,intent(out):: det,bee(:,:)
 integer,intent(in):: i
 real:: x1,x2,x3,x4,y1,y2,y3,y4,xi,et,&
        x12,x13,x14,x23,x24,x34,y12,y13,y14,y23,y24,y34,&
        xy12,xy13,xy14,xy23,xy24,xy34,xifac,etfac,const,den,&
        dn1dx,dn2dx,dn3dx,dn4dx,dn1dy,dn2dy,dn3dy,dn4dy  

 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)
 xi=points(i,1); et=points(i,2)
 x12=x1-x2; x13=x1-x3; x14=x1-x4; x23=x2-x3; x24=x2-x4; x34=x3-x4
 y12=y1-y2; y13=y1-y3; y14=y1-y4; y23=y2-y3; y24=y2-y4; y34=y3-y4
 xy12=x1*y2-y1*x2; xy13=x1*y3-y1*x3; xy14=x1*y4-y1*x4; xy23=x2*y3-y2*x3
 xy24=x2*y4-y2*x4; xy34=x3*y4-y3*x4
 xifac= xy12-xy13+xy24-xy34; etfac= xy13-xy14-xy23+xy24
 const=-xy12+xy14-xy23-xy34; den=xifac*xi+etfac*et+const; det=0.125*den
 dn1dx=( y23*xi+y34*et-y24)/den; dn2dx=(-y14*xi-y34*et+y13)/den
 dn3dx=( y14*xi-y12*et+y24)/den; dn4dx=(-y23*xi+y12*et-y13)/den
 dn1dy=(-x23*xi-x34*et+x24)/den; dn2dy=( x14*xi+x34*et-x13)/den
 dn3dy=(-x14*xi+x12*et-x24)/den; dn4dy=( x23*xi-x12*et+x13)/den
 bee(1,1)=dn1dx; bee(1,2)=0.; bee(1,3)=dn2dx; bee(1,4)=0.
 bee(1,5)=dn3dx; bee(1,6)=0.; bee(1,7)=dn4dx; bee(1,8)=0.
 bee(2,1)=0.; bee(2,2)=dn1dy; bee(2,3)=0.; bee(2,4)=dn2dy
 bee(2,5)=0.; bee(2,6)=dn3dy; bee(2,7)=0.; bee(2,8)=dn4dy
 bee(3,1)=dn1dy; bee(3,2)=dn1dx; bee(3,3)=dn2dy; bee(3,4)=dn2dx
 bee(3,5)=dn3dy; bee(3,6)=dn3dx; bee(3,7)=dn4dy; bee(3,8)=dn4dx
return
end subroutine bee4           
subroutine seep4(coord,perm,kc)

!     this contains the analytical conductivity mstrix
!     for a four node quadrilateral

 implicit none
 real,intent(in):: coord(:,:),perm(:,:)
 real,intent(out):: kc(:,:)
 integer::i,j
 real:: x1,x2,x3,x4,y1,y2,y3,y4,y31,y42,x31,x42,x34,y34
 real:: a2,f1,f2,kx,ky,ht(4),it(4),temp,s1,s2,t1,t2,t3,gdiag,goppo
 
! input data
 
 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)
 kx=perm(1,1); ky=perm(2,2)
 
! procedure definitions
 
 y31=(y3-y1); y42=(y4-y2); x31=(x3-x1); x42=(x4-x2)
 x34=x3-x4; y34=y3-y4; ht(1)=x4-x1; ht(2)=x1-x2
 ht(3)=x2-x3; it(1)=y4-y1; it(2)=y1-y2; it(3)=y2-y3
 
 a2=y42*x31-y31*x42
 f1=2*y42*x34-2*x42*y34-a2; f2=2*y31*x34-2*x31*y34-a2
 s1=kx*y42**2+ky*x42**2; s2=kx*y34**2+ky*x34**2
 t1=kx*y42*y34+ky*x42*x34; t2=kx*y31*y34+ky*x31*x34; t3=kx*y42*y31+x31*x42*ky
 
! calculate parent terms
 
 gdiag=(-(2.*a2+f1)*s1)/(2.*(3.*a2**2-f1**2))
 goppo=-(s1*(2.*a2+f2)+2.*t1*(a2+f2)+2.*a2*s2)/(2.*(3.*a2**2-f2**2))
 kc(1,1)=gdiag+goppo
!
 gdiag=((2.*a2+f1)*t3-(a2+f1)*t1)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*t3+(a2+f2)*t2)/(2.*(3.*a2**2-f2**2))
 kc(2,1)=gdiag+goppo
!
 gdiag=(s1*a2)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*s1+(a2+f2)*(2.*t1-t3)+2.*a2*(s2-t2))/(2.*(3.*a2**2-f2**2))
 kc(3,1)=gdiag+goppo
 
! calculate elements using parent terms and transforms
 
 temp=y31; y31=y42; y42=-temp; temp=x31; 
 x31=x42; x42=-temp; x34=ht(1); y34=it(1)
!
 f1=2*y42*x34-2*x42*y34-a2; f2=2*y31*x34-2*x31*y34-a2
 s1=kx*y42**2+ky*x42**2; s2=kx*y34**2+ky*x34**2
 t1=kx*y42*y34+ky*x42*x34; t2=kx*y31*y34+ky*x31*x34; t3=kx*y42*y31+x31*x42*ky
!
 gdiag=(-(2.*a2+f1)*s1)/(2.*(3.*a2**2-f1**2))
 goppo=-(s1*(2.*a2+f2)+2.*t1*(a2+f2)+2.*a2*s2)/(2.*(3.*a2**2-f2**2))
 kc(2,2)=gdiag+goppo
!
 gdiag=((2.*a2+f1)*t3-(a2+f1)*t1)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*t3+(a2+f2)*t2)/(2.*(3.*a2**2-f2**2))
 kc(3,2)=gdiag+goppo
!
 gdiag=(s1*a2)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*s1+(a2+f2)*(2.*t1-t3)+2.*a2*(s2-t2))/(2.*(3.*a2**2-f2**2))
 kc(4,2)=gdiag+goppo
!
 temp=y31; y31=y42; y42=-temp; temp=x31
 x31=x42; x42=-temp; x34=ht(2); y34=it(2)
!
 f1=2*y42*x34-2*x42*y34-a2; f2=2*y31*x34-2*x31*y34-a2
 s1=kx*y42**2+ky*x42**2; s2=kx*y34**2+ky*x34**2
 t1=kx*y42*y34+ky*x42*x34; t2=kx*y31*y34+ky*x31*x34; t3=kx*y42*y31+x31*x42*ky
!
 gdiag=(-(2.*a2+f1)*s1)/(2.*(3.*a2**2-f1**2))
 goppo=-(s1*(2.*a2+f2)+2.*t1*(a2+f2)+2.*a2*s2)/(2.*(3.*a2**2-f2**2))
 kc(3,3)=gdiag+goppo
!
 gdiag=((2.*a2+f1)*t3-(a2+f1)*t1)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*t3+(a2+f2)*t2)/(2.*(3.*a2**2-f2**2))
 kc(4,3)=gdiag+goppo
!
 temp=y31; y31=y42; y42=-temp; temp=x31
 x31=x42; x42=-temp; x34=ht(3); y34=it(3)
!
 f1=2*y42*x34-2*x42*y34-a2; f2=2*y31*x34-2*x31*y34-a2
 s1=kx*y42**2+ky*x42**2; s2=kx*y34**2+ky*x34**2
 t1=kx*y42*y34+ky*x42*x34; t2=kx*y31*y34+ky*x31*x34; t3=kx*y42*y31+x31*x42*ky
!
 gdiag=(-(2.*a2+f1)*s1)/(2.*(3.*a2**2-f1**2))
 goppo=-(s1*(2.*a2+f2)+2.*t1*(a2+f2)+2.*a2*s2)/(2.*(3.*a2**2-f2**2))
 kc(4,4)=gdiag+goppo
!
 gdiag=((2.*a2+f1)*t3-(a2+f1)*t1)/(2.*(3.*a2**2-f1**2))
 goppo=((2.*a2+f2)*t3+(a2+f2)*t2)/(2.*(3.*a2**2-f2**2))
 kc(4,1)=gdiag+goppo
 
! mirror matrix about main diagonal
 
 do i=2,4; do j=1,i-1; kc(j,i)=kc(i,j); enddo; enddo

 return
 end subroutine seep4         
subroutine loc_to_glob(local,global,gamma,coord)
!
!      this subroutine transforms the local end reactions and
!      moments into the global system (3-d)
!
implicit none
real,intent(in)::local(:),gamma,coord(:,:)
real,intent(out)::global(:)
real::t(12,12),r0(3,3),x1,x2,y1,y2,z1,z2,xl,yl,zl,&
      pi,gamrad,cg,sg,den,ell,x,sum
integer::i,j,k
     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)
     t=0.0
     pi=acos(-1.); gamrad=gamma*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(2,1)=yl/ell; r0(3,1)=zl/ell
      r0(1,2)=(-xl*yl*cg-ell*zl*sg)/den; r0(2,2)=den*cg/(ell*ell)
      r0(3,2)=(-yl*zl*cg+ell*xl*sg)/den; r0(1,3)=(xl*yl*sg-ell*zl*cg)/den
      r0(2,3)=-den*sg/(ell*ell); r0(3,3)=(yl*zl*sg+ell*xl*cg)/den
     else
      r0(1,1)=0.; r0(3,1)=0.; r0(2,2)=0.; r0(2,3)=0.
      r0(2,1)=1.; r0(1,2)=-cg; r0(3,3)=cg; r0(3,2)=sg; r0(1,3)=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)*local(j)
       end do
      global(i)=sum
     end do
      return
end subroutine loc_to_glob    
subroutine glob_to_loc(local,global,gamma,coord)
!
!      this subroutine transforms the global end reactions and
!      moments into the local system (2-d, 3-d)
!
implicit none
real,intent(in)::global(:),gamma,coord(:,:)
real,intent(out)::local(:)
real::t(12,12),r0(3,3),x1,x2,y1,y2,z1,z2,xl,yl,zl,&
      pi,gamrad,cg,sg,den,ell,x,sum
integer::i,j,k,ndim
ndim=ubound(coord,2)
select case(ndim)
  case(2)
      x1=coord(1,1); y1=coord(1,2)
      x2=coord(2,1); y2=coord(2,2)
      ell=sqrt((x2-x1)**2+(y2-y1)**2)
      cg=(x2-x1)/ell
      sg=(y2-y1)/ell
      local(1)=cg*global(1)+sg*global(2)
      local(2)=cg*global(2)-sg*global(1)
      local(3)=global(3)
      local(4)=cg*global(4)+sg*global(5)
      local(5)=cg*global(5)-sg*global(4)
      local(6)=global(6)
  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)
     t=0.0
     pi=acos(-1.); gamrad=gamma*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.