geometry.f90

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

module geometry_lib
contains                                          
!----------------Node to freedom number conversion ---------------------------- 
 subroutine num_to_g(num,nf,g)
 !finds the g vector from num and nf
 implicit none
 integer,intent(in)::num(:),nf(:,:)  ; integer,intent(out):: g(:)
 integer::i,k,nod,nodof ; nod=ubound(num,1) ; nodof=ubound(nf,1)
  do i = 1 , nod
      k = i*nodof  ; g(k-nodof+1:k) = nf( : , num(i) )
  end do
 return
 end subroutine num_to_g   
!-------------------------------- Lines  --------------------------------------
subroutine geometry_2l(iel,ell,coord,num)
! node numbers, nodal coordinates and steering vectors for
! a line of (nonuniform) beam elements
 implicit none
  integer,intent(in)::iel; integer,intent(out)::num(:)
  real,intent(in)::ell; real,intent(out)::coord(:,:)
  num=(/iel,iel+1/)
  if(iel==1)then;coord(1,1)=.0; coord(2,1)=ell; else
    coord(1,1)=coord(2,1); coord(2,1)=coord(2,1)+ell; end if
 end subroutine geometry_2l 
!----------------------------- Triangles --------------------------------------
subroutine geometry_3tx(iel,nxe,aa,bb,coord,num)
!      this subroutine forms the coordinates and node vector
!      for a rectangular mesh of uniform 3-node triangles
!      counting in the x-direction    ; local numbering clockwise
  implicit none
  real,intent(in):: aa,bb; integer,intent(in):: iel,nxe
  real,intent(out) :: coord(:,:); integer,intent(out):: num(:)
  integer::ip,iq,jel
      jel= (2*nxe+iel-1)/(2*nxe)
      if(iel/2*2==iel)then; iq=2*jel; else; iq=2*jel-1; end if
      ip= (iel-2*nxe*(jel-1)+1)/2
      if(mod(iq,2)/=0) then
       num(1)=(nxe+1)*(iq-1)/2+ip  ;   num(3)=(nxe+1)*(iq+1)/2+ip
       num(2)=num(1)+1              ;   coord(1,1)=(ip-1)*aa
       coord(1,2)=-(iq-1)/2*bb      ;   coord(3,1)=(ip-1)*aa
       coord(2,2)=coord(1,2)        ;   coord(2,1)=ip*aa
       coord(3,2)=-(iq+1)/2*bb
      else
       num(1)=(nxe+1)*iq/2+ip+1     ;   num(3)=(nxe+1)*(iq-2)/2+ip+1
       num(2)=num(1)-1               ;   coord(1,1)=ip*aa
       coord(1,2)=-iq/2*bb           ;   coord(3,1)=ip*aa
       coord(3,2)=-(iq-2)/2*bb       ;   coord(2,1)=(ip-1)*aa
       coord(2,2)=coord(1,2)
      end if
     return
   end subroutine geometry_3tx
subroutine geometry_6tx(iel,nxe,aa,bb,coord,num)
!      this subroutine forms the coordinates and nodal vector
!      for a rectangular mesh of uniform 6-node triangles
!      counting in the x-direction ; local numbering clockwise
  implicit none
  real ,intent(in):: aa,bb; integer,intent(in):: iel,nxe
  real,intent(out) :: coord(:,:); integer,intent(out)::  num(:)
  integer::ip,iq,jel,i   
      jel= (2*nxe+iel-1)/(2*nxe)
      if(iel/2*2==iel)then; iq=2*jel; else; iq=2*jel-1; end if
      ip= (iel-2*nxe*(jel-1)+1)/2
      if(mod(iq,2)/=0) then
       num(1)=(iq-1)*(2*nxe+1)+2*ip-1     ;   num(2)=num(1)+1
       num(3)=num(1)+2                    ;   num(4)=num(2)+1 
       num(6)=(iq-1)*(2*nxe+1)+2*nxe+2*ip ;   num(5)=(iq+1)*(2*nxe+1)+2*ip-1
       coord(1,1)=(ip-1)*aa               ;   coord(1,2)=-(iq-1)/2*bb
       coord(5,1)=(ip-1)*aa               ;   coord(5,2)=-(iq+1)/2*bb
       coord(3,1)=ip*aa                   ;   coord(3,2)=coord(1,2)
      else
       num(1)=iq*(2*nxe+1)+2*ip+1       ; num(6)=(iq-2)*(2*nxe+1)+2*nxe+2*ip+2
       num(5)=(iq-2)*(2*nxe+1)+2*ip+1    ; num(4)=num(2)-1
       num(3)=num(1)-2                   ; num(2)=num(1)-1
       coord(1,1)=ip*aa                  ; coord(1,2)=-iq/2*bb
       coord(5,1)=ip*aa                  ; coord(5,2)=-(iq-2)/2*bb
       coord(3,1)=(ip-1)*aa              ; coord(3,2)=coord(1,2)
      end if
      do  i=1,2
      coord(2,i)=.5*(coord(1,i)+coord(3,i))
      coord(4,i)=.5*(coord(3,i)+coord(5,i))
      coord(6,i)=.5*(coord(5,i)+coord(1,i))
      end do
      return
   end subroutine geometry_6tx
  subroutine geometry_15tyv(iel,nye,width,depth,coord,num)
!      this subroutine forms the coordinates and node vector
!      for a rectangular mesh of nonuniform 15-node triangles
!      counting in the y-direction ; local numbering clockwise
 implicit none
 real,intent(in) :: width(:),depth(:) ; real,intent(out) :: coord(:,:)
 integer, intent(in) ::  iel,nye; integer,intent(out)::num(:)
 integer::ip,iq,jel,i , fac1,fac2   
   jel = (iel - 1)/nye; ip= jel+1; iq=iel-nye*jel
    if(mod(iq,2)/=0) then
      fac1=4*(2*nye+1)*(ip-1)+2*iq-1 ; num(1)=fac1;  num(12)=fac1+1
      num(11)=fac1+2 ;  num(10)=fac1+3 ; num(9)=fac1+4 ; num(8)=fac1+2*nye+4
      num(7)=fac1+4*nye+4 ;  num(6)=fac1+6*nye+4 ;  num(5)=fac1+8*nye+4
      num(4)=fac1+6*nye+3 ;  num(3)=fac1+4*nye+2 ;  num(2)=fac1+2*nye+1
      num(13)=fac1+2*nye+2 ;  num(15)=fac1+2*nye+3  ;  num(14)=fac1+4*nye+3
      coord(1,1)=width(ip) ;  coord(1,2)=depth((iq+1)/2)
      coord(9,1)=width(ip)   ;  coord(9,2)=depth((iq+3)/2)
      coord(5,1)=width(ip+1)   ;  coord(5,2)=depth((iq+1)/2)
    else
      fac2=4*(2*nye+1)*(ip-1)+2*iq+8*nye+5 ;  num(1)=fac2 ;  num(12)=fac2-1
      num(11)=fac2-2 ;  num(10)=fac2-3 ;  num(9)=fac2-4; num(8)=fac2-2*nye-4
      num(7)=fac2-4*nye-4 ; num(6)=fac2-6*nye-4  ;  num(5)=fac2-8*nye-4
      num(4)=fac2-6*nye-3 ; num(3)=fac2-4*nye-2 ; num(2)=fac2-2*nye-1
      num(13)=fac2-2*nye-2  ; num(15)=fac2-2*nye-3 ; num(14)=fac2-4*nye-3
      coord(1,1)=width(ip+1) ;  coord(1,2)=depth((iq+2)/2)
      coord(9,1)=width(ip+1) ;  coord(9,2)=depth(iq/2)
      coord(5,1)=width(ip)   ;   coord(5,2)=depth((iq+2)/2)
    end if
    do  i=1,2
      coord(3,i)=.5*(coord(1,i)+coord(5,i))
      coord(7,i)=.5*(coord(5,i)+coord(9,i))
      coord(11,i)=.5*(coord(9,i)+coord(1,i))
      coord(2,i)=.5*(coord(1,i)+coord(3,i))
      coord(4,i)=.5*(coord(3,i)+coord(5,i))
      coord(6,i)=.5*(coord(5,i)+coord(7,i))
      coord(8,i)=.5*(coord(7,i)+coord(9,i))
      coord(10,i)=.5*(coord(9,i)+coord(11,i))
      coord(12,i)=.5*(coord(11,i)+coord(1,i))
      coord(15,i)=.5*(coord(7,i)+coord(11,i))
      coord(14,i)=.5*(coord(3,i)+coord(7,i))
      coord(13,i)=.5*(coord(2,i)+coord(15,i))
    end do
  return
 end subroutine geometry_15tyv  
!---------------------- Quadrilaterals ----------------------------------------
subroutine geometry_4qx(iel,nxe,aa,bb,coord,num)
! coordinates and nodal vectors for equal four node quad
! elements, numbering in x 
implicit none
integer,intent(in)::iel,nxe; real,intent(in)::aa,bb
real,intent(out)::coord(:,:); integer,intent(out)::num(:)
integer :: ip,iq    ; iq=(iel-1)/nxe+1; ip=iel-(iq-1)*nxe
   num=(/iq*(nxe+1)+ip,(iq-1)*(nxe+1)+ip,   &
         (iq-1)*(nxe+1)+ip+1, iq*(nxe+1)+ip+1/)
   coord(1:2,1)=(ip-1)*aa; coord(3:4,1)=ip*aa
   coord(1:4:3,2)=-iq*bb; coord(2:3,2)=-(iq-1)*bb
 return
end subroutine geometry_4qx
   subroutine geometry_4qy(iel,nye,aa,bb,coord,num)
   ! rectangles of equal 4-node quads numbered in y
   implicit none
   integer,intent(in)::iel,nye; real,intent(in)::aa,bb
   real,intent(out)::coord(:,:);integer,intent(out)::num(:)
       num=(/iel+(iel-1)/nye+1,iel+(iel-1)/nye,iel+(iel-1)/nye+nye+1,  &
             iel+(iel-1)/nye+nye+2/)
       coord(1:2,1)= aa*((iel-1)/nye); coord(3:4,1)=aa*((iel-1)/nye+1)
       coord(1:4:3,2)=-(iel-((iel-1)/nye)*nye)*bb; coord(2:3,2)=coord(1,2)+bb
     return
   end subroutine geometry_4qy              
subroutine geometry_4qyv(iel,nye,width,depth,coord,num)
!  coordinates and steering vector for a variable rectangular
!  mesh of 4-node quad elements numbering in the y-direction
  implicit none
  real,intent(in)::width(:),depth(:); integer,intent(in)::iel,nye
  real,intent(out)::coord(:,:); integer,intent(out):: num(:)
  integer:: ip,iq; ip=(iel-1)/nye+1; iq=iel-(ip-1)*nye
  num(1)=(ip-1)*(nye+1)+iq+1; num(2)=num(1)-1
  num(3)=ip*(nye+1)+iq;num(4)= num(3) + 1
  coord(1:2,1)=width(ip); coord(3:4,1)=width(ip+1)
  coord(1,2)=depth(iq+1); coord(2:3,2)=depth(iq); coord(4,2)=coord(1,2)
 return
end subroutine geometry_4qyv
subroutine geometry_8qx(iel,nxe,aa,bb,coord,num)
! coordinates and steering vector for a rectangular mesh of
! equal  8-node  elements  numbering in x
implicit none
 real,intent(out)::coord(:,:); integer,intent(out)::num(:)
 integer,intent(in)::iel,nxe; real,intent(in)::aa,bb
 integer:: ip,iq ; iq=(iel-1)/nxe+1; ip=iel-(iq-1)*nxe
 num(1)=iq*(3*nxe+2)+2*ip-1; num(2)=iq*(3*nxe+2)+ip-nxe-1
 num(3)=(iq-1)*(3*nxe+2)+2*ip-1; num(4)=num(3)+1
 num(5)=num(4)+1; num(6)=num(2)+1; num(7)=num(1)+2; num(8)=num(1)+1
 coord(1:3,1)=aa*(ip-1); coord(5:7,1)=aa*ip
 coord(4,1)=.5*(coord(3,1)+coord(5,1))
 coord(8,1)=.5*(coord(7,1)+coord(1,1))
 coord(1,2)=-bb*iq; coord(7:8,2)=-bb*iq
 coord(3:5,2)=-bb*(iq-1); coord(2,2)=.5*(coord(1,2)+coord(3,2))
 coord(6,2)=.5*(coord(5,2)+coord(7,2))  
return
end subroutine geometry_8qx
subroutine geometry_8qy(iel,nye,aa,bb,coord,num)
!  coordinates and steering vector for a constant rectangular
!  mesh of 8-node quad elements numbering in the y-direction
  implicit none
  real,intent(in):: aa,bb ; integer,intent(in):: iel,nye
  real,intent(out)::coord(:,:); integer,intent(out):: num(:)
  integer:: ip,iq; ip=(iel-1)/nye+1; iq=iel-(ip-1)*nye
  num(1)=(ip-1)*(3*nye+2)+2*iq+1; num(2)=num(1)-1; num(3)=num(1)-2
  num(4)=(ip-1)*(3*nye+2)+2*nye+iq+1;num(5)=ip*(3*nye+2)+2*iq-1
  num(6)=num(5)+1; num(7)=num(5)+2; num(8)=num(4)+1
  coord(1:3,1)=(ip-1)*aa; coord(5:7,1)=ip*aa
  coord(4,1)=.5*(coord(3,1)+coord(5,1))
  coord(8,1)=.5*(coord(7,1)+coord(1,1))
  coord(1,2)=-iq*bb; coord(7:8,2)=-iq*bb; coord(3:5,2)=-(iq-1)*bb
  coord(2,2)=.5*(coord(1,2)+coord(3,2))
  coord(6,2)=.5*(coord(5,2)+coord(7,2))
 return
end subroutine geometry_8qy
subroutine geometry_8qxv(iel,nxe,width,depth,coord,num)
! nodal coordinates and node vector for a variable mesh of
! 8-node quadrilaterals numbering in the x-direction
implicit none
  integer,intent(in)::iel,nxe;real,intent(in)::width(:),depth(:)
  real,intent(out)::coord(:,:); integer,intent(out)::num(:)                
  integer::ip,iq; iq=(iel-1)/nxe+1; ip=iel-(iq-1)*nxe
  num(1)=iq*(3*nxe+2)+2*ip-1; num(2)=iq*(3*nxe+2)+ip-nxe-1
  num(3)=(iq-1)*(3*nxe+2)+2*ip-1; num(4)=num(3)+1;num(5)=num(4)+1
  num(6)=num(2)+1; num(7)=num(1)+2; num(8)=num(1)+1
  coord(1:3,1)=width(ip); coord(5:7,1)=width(ip+1)
  coord(4,1)=.5*(coord(3,1)+coord(5,1));coord(8,1)=.5*(coord(7,1)+coord(1,1))
  coord(1,2)=depth(iq+1); coord(7:8,2)=depth(iq+1); coord(3:5,2)=depth(iq)
  coord(2,2)=.5*(coord(1,2)+coord(3,2));coord(6,2)=.5*(coord(5,2)+coord(7,2))
 return
end subroutine geometry_8qxv    
subroutine geometry_8qyv(iel,nye,width,depth,coord,num)
!  coordinates and steering vector for a variable rectangular
!  mesh of 8-node quad elements numbering in the y-direction
  implicit none
  real,intent(in)::width(:),depth(:); integer,intent(in)::iel,nye
  real,intent(out)::coord(:,:); integer,intent(out)::num(:)
  integer::ip,iq; ip=(iel-1)/nye+1; iq=iel-(ip-1)*nye
  num(1)=(ip-1)*(3*nye+2)+2*iq+1; num(2)=num(1)-1; num(3)=num(1)-2
  num(4)=(ip-1)*(3*nye+2)+2*nye+iq+1;num(5)=ip*(3*nye+2)+2*iq-1
  num(6)=num(5)+1; num(7)=num(5)+2; num(8)=num(4)+1
  coord(1:3,1)=width(ip); coord(5:7,1)=width(ip+1)
  coord(4,1)=.5*(coord(3,1)+coord(5,1))
  coord(8,1)=.5*(coord(7,1)+coord(1,1))
  coord(1,2)=depth(iq+1); coord(7:8,2)=depth(iq+1); coord(3:5,2)=depth(iq)
  coord(2,2)=.5*(coord(1,2)+coord(3,2))
  coord(6,2)=.5*(coord(5,2)+coord(7,2))
 return
end subroutine geometry_8qyv    
 subroutine geometry_9qx(iel,nxe,aa,bb,coord,num)
!      this subroutine forms the coordinates and steering vector
!      for equal 9-node Lagrangian quads counting in x-direction
 implicit none
 real,intent(out)::coord(:,:); integer,intent(out)::num(:)
 integer,intent(in)::iel,nxe; real,intent(in)::aa,bb
 integer:: ip,iq ;iq=(iel-1)/nxe+1;ip=iel-(iq-1)*nxe
   num(1)=iq*(4*nxe+2)+2*ip-1 ; num(2)=iq*(4*nxe+2)+2*ip-nxe-4
   num(3)= (iq-1)*(4*nxe+2)+2*ip-1 ;   num(4)=num(3)+1
   num(5)=num(4)+1; num(6)=num(2)+2 ;  num(7)=num(1)+2
   num(8)=num(1)+1      ;   num(9)=num(2)+1
   coord(1,1)=(ip-1)*aa  ;  coord(3,1)=(ip-1)*aa   ; coord(5,1)=ip*aa
   coord(7,1)=ip*aa ;   coord(1,2)=-iq*bb ;  coord(3,2)=-(iq-1)*bb
   coord(5,2)=-(iq-1)*bb    ;  coord(7,2)=-iq*bb
   coord(2,1)=.5*(coord(1,1)+coord(3,1)); coord(2,2)=.5*(coord(1,2)+coord(3,2))
   coord(4,1)=.5*(coord(3,1)+coord(5,1)); coord(4,2)=.5*(coord(3,2)+coord(5,2))
   coord(6,1)=.5*(coord(5,1)+coord(7,1)); coord(6,2)=.5*(coord(5,2)+coord(7,2))
   coord(8,1)=.5*(coord(1,1)+coord(7,1)); coord(8,2)=.5*(coord(1,2)+coord(7,2))
   coord(9,1)=.5*(coord(2,1)+coord(6,1)); coord(9,2)=.5*(coord(4,2)+coord(8,2))
  return
 end subroutine geometry_9qx  
!-----------------------Hexahedra "Bricks" ------------------------------------
 subroutine geometry_8bxz(iel,nxe,nze,aa,bb,cc,coord,num)
!      this subroutine forms the coordinates and nodal vector
!      for boxes of 8-node brick elements counting x-z planes in y-direction
 implicit none
   integer,intent(in)::iel,nxe,nze;integer,intent(out)::num(:)
   real,intent(in)::aa,bb,cc; real,intent(out)::coord(:,:)
   integer::ip,iq,is,iplane
   iq=(iel-1)/(nxe*nze)+1 ; iplane = iel -(iq-1)*nxe*nze
   is=(iplane-1)/nxe+1; ip = iplane-(is-1)*nxe   
   num(1)=(iq-1)*(nxe+1)*(nze+1)+is*(nxe+1)+ip ;  num(2)=num(1)-nxe-1
   num(3)=num(2)+1 ;   num(4)=num(1)+1 ; num(5)=num(1)+(nxe+1)*(nze+1)
   num(6)=num(5)-nxe-1 ;   num(7)=num(6)+1    ;   num(8)=num(5)+1
   coord(1,1)=(ip-1)*aa ; coord(2,1)=(ip-1)*aa ;  coord(5,1)=(ip-1)*aa
   coord(6,1)=(ip-1)*aa ; coord(3,1)=ip*aa ;   coord(4,1)=ip*aa
   coord(7,1)=ip*aa   ;   coord(8,1)=ip*aa
   coord(1,2)=(iq-1)*bb  ;   coord(2,2)=(iq-1)*bb ; coord(3,2)=(iq-1)*bb
   coord(4,2)=(iq-1)*bb  ;   coord(5,2)=iq*bb  ;    coord(6,2)=iq*bb
   coord(7,2)=iq*bb ;    coord(8,2)=iq*bb ;   coord(1,3)=-is*cc
   coord(4,3)=-is*cc  ;   coord(5,3)=-is*cc ;   coord(8,3)=-is*cc
   coord(2,3)=-(is-1)*cc ; coord(3,3)=-(is-1)*cc  ;  coord(6,3)=-(is-1)*cc
   coord(7,3)=-(is-1)*cc