geometry.f90

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

  return
 end subroutine geometry_8bxz                     
subroutine geometry_20bxz(iel,nxe,nze,aa,bb,cc,coord,num)
! nodal vector and nodal coordinates for boxes of 20-node
! bricks counting x-z planes in the y-direction
implicit none
  integer,intent(in)::iel,nxe,nze; real,intent(in)::aa,bb,cc
  real,intent(out)::coord(:,:); integer,intent(out)::num(:)
  integer::fac1,fac2,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
  fac1=((2*nxe+1)*(nze+1)+(2*nze+1)*(nxe+1))*(iq-1)
  fac2=((2*nxe+1)*(nze+1)+(2*nze+1)*(nxe+1))*iq
  num(1)=fac1+(3*nxe+2)*is+2*ip-1
  num(2)=fac1+(3*nxe+2)*is-nxe+ip-1; num(3)=num(1)-3*nxe-2
  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
  num(9)=fac2-(nxe+1)*(nze+1)+(nxe+1)*is+ip
  num(10)=num(9)-nxe-1; num(11)=num(10)+1; num(12)=num(9)+1
  num(13)=fac2+(3*nxe+2)*is+2*ip-1
  num(14)=fac2+(3*nxe+2)*is-nxe+ip-1
  num(15)=num(13)-3*nxe-2; num(16)=num(15)+1; num(17)=num(16)+1
  num(18)=num(14)+1; num(19)=num(13)+2; num(20)=num(13)+1 
  coord(1:3,1)=(ip-1)*aa; coord(9:10,1)=(ip-1)*aa; coord(13:15,1)=(ip-1)*aa
  coord(5:7,1)=ip*aa; coord(11:12,1)=ip*aa; coord(17:19,1)=ip*aa
  coord(4,1)=.5*(coord(3,1)+coord(5,1));coord(8,1)=.5*(coord(1,1)+coord(7,1))
  coord(16,1)=.5*(coord(15,1)+coord(17,1))
  coord(20,1)=.5*(coord(13,1)+coord(19,1))
  coord(1:8,2)=(iq-1)*bb; coord(13:20,2)=iq*bb
  coord(9,2)=.5*(coord(1,2)+coord(13,2))
  coord(10,2)=.5*(coord(3,2)+coord(15,2))
  coord(11,2)=.5*(coord(5,2)+coord(17,2))
  coord(12,2)=.5*(coord(7,2)+coord(19,2))
  coord(1,3)=-is*cc; coord(7:9,3)=-is*cc; coord(12:13,3)=-is*cc
  coord(19:20,3)=-is*cc; coord(3:5,3)=-(is-1)*cc
  coord(10:11,3)=-(is-1)*cc; coord(15:17,3)=-(is-1)*cc
  coord(2,3)=.5*(coord(1,3)+coord(3,3))
  coord(6,3)=.5*(coord(5,3)+coord(7,3))
  coord(14,3)=.5*(coord(13,3)+coord(15,3))
  coord(18,3)=.5*(coord(17,3)+coord(19,3))
 return
end subroutine geometry_20bxz
!--------------------Special purpose geometries  ------------------------------
 subroutine slope_geometry(iel,nye,top,bot,depth,coord,num)
!    this subroutine forms the coordinates and steering vector
!    for 8-node quadrilaterals in a 'slope' geometry
!    (numbering in the y-direction)
 implicit none
  real,intent(in):: top(:),bot(:),depth(:); real,intent(out)::coord(:,:)
  integer,intent(in):: iel,nye; integer,intent(out)::num(:)
  real :: fac1 , fac2  ; 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
   fac1=(bot(ip)-top(ip))/nye ;  fac2=(bot(ip+1)-top(ip+1))/nye
   coord(1,1)=top(ip)+iq*fac1;   coord(3,1)=top(ip)+(iq-1)*fac1
   coord(5,1)=top(ip+1)+(iq-1)*fac2; coord(7,1)=top(ip+1)+iq*fac2
   coord(2,1)=.5*(coord(1,1)+coord(3,1)); coord(6,1)=.5*(coord(5,1)+coord(7,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(8,2)=depth(iq+1); coord(7,2)=depth(iq+1)
   coord(3,2)=depth(iq) ;  coord(4,2)=depth(iq)  ; coord(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 slope_geometry    
subroutine geometry_freesurf(iel,nxe,fixed_seep,fixed_down,down,&
                             width,angs,surf,coord,num)
!             this subroutine forms the coordinates and steering vector
!             for 4-node quads numbering in the x-direction
!             (Laplace's equation, variable mesh, 1-freedom per node)
implicit none
real,intent(in)::width(:),angs(:),surf(:),down ;real,intent(out)::coord(:,:)
integer,intent(in)::iel,nxe,fixed_seep,fixed_down;integer,intent(out)::num(:)
real::angr(size(angs)),pi,b1,b2,tan1,tan2,fac1,fac2
integer::ip,iq
 pi=acos(-1.0); angr=angs*pi/180.    ; iq=(iel-1)/nxe+1; ip=iel-(iq-1)*nxe
 num(1)=iq*(nxe+1)+ip; num(2)=(iq-1)*(nxe+1)+ip
 num(3)=num(2)+1; num(4)=num(1)+1
 if(iq <= fixed_seep+1)then
 b1=(surf(ip)-down)/real(fixed_seep+1); b2=(surf(ip+1)-down)/real(fixed_seep+1)
 coord(1,2)=down+(fixed_seep+1-iq)*b1 ; coord(2,2)=down+(fixed_seep+2-iq)*b1
 coord(3,2)=down+(fixed_seep+2-iq)*b2 ; coord(4,2)=down+(fixed_seep+1-iq)*b2
 else
 b1=real(fixed_down+fixed_seep-iq)/real(fixed_down-1) 
 b2=real(fixed_down+fixed_seep-iq+1)/real(fixed_down-1) 
 coord(1,2)=down*b1;    coord(2,2)=down*b2
 coord(3,2)=coord(2,2); coord(4,2)=coord(1,2)
 end if
if(abs(angr(ip)-pi*.5) < 0.0001)then
  fac1=0.
else
 tan1=tan(angr(ip))      ; fac1=1.0/tan1
end if
if(abs(angr(ip+1)-pi*.5) < 0.0001)then
  fac2=0.
else
 tan2=tan(angr(ip+1))    ; fac2=1.0/tan2 
end if
 coord(1,1)=width(ip) +coord(1,2)*fac1 ; coord(2,1)=width(ip)  +coord(2,2)*fac1
 coord(3,1)=width(ip+1)+coord(3,2)*fac2 ;coord(4,1)=width(ip+1)+coord(4,2)*fac2
return
end subroutine geometry_freesurf                                               
!---------------------Construction and excavation -----------------------------
 subroutine fmglem(fnxe,fnye,lnxe,lnye,g_num,lifts)
 ! returns g_num for the mesh
 implicit none
 integer,intent(in)::fnxe,fnye,lnxe,lnye,lifts;integer,intent(out)::g_num(:,:)
 integer::ig,i,j,ii,ilast       ; ig = 0
 do i=1,fnye; do j=1,fnxe
  ig = ig + 1
  g_num(1,ig)=(j*2-1)+(i-1)*(fnxe+1)+(i-1)*(2*fnxe+1)
  g_num(8,ig) = g_num(1,ig)+1; g_num(7,ig) = g_num(1,ig)+2
  g_num(2,ig)=i*(2*fnxe+1)+j+(i-1)*(fnxe+1)
  g_num(6,ig)=g_num(2,ig)+1
  g_num(3,ig)=i*(2*fnxe+1)+i*(fnxe+1)+(j*2-1)
  g_num(4,ig)=g_num(3,ig)+1; g_num(5,ig)=g_num(3,ig)+2
 end do; end do
  ilast=g_num(5,ig)
  do ii = 1 , lifts - 1
    do i=1,lnye; do j=1,lnxe-(ii-1)*1
      ig = ig + 1
      g_num(1,ig)=ilast-(lnxe-ii+1)*2+(j-1)*2
      g_num(8,ig) = g_num(1,ig)+1; g_num(7,ig) = g_num(1,ig)+2
      g_num(2,ig)=ilast+j
      g_num(6,ig)=g_num(2,ig)+1
      g_num(3,ig)=ilast+(lnxe-ii+1)+1+(j*2-1)
      g_num(4,ig)=g_num(3,ig)+1; g_num(5,ig)=g_num(3,ig)+2
    end do; end do
    ilast=g_num(5,ig)
  end do
 return
 end subroutine fmglem
 subroutine fmcoem(g_num,g_coord,fwidth,fdepth,width,depth,           &
                   lnxe,lifts,fnxe,fnye,itype)
 ! returns g_coord for the mesh
 implicit none
 integer,intent(in)::g_num(:,:),lnxe,lifts,fnxe,fnye,itype
 real,intent(out)::g_coord(:,:);integer::ig,i,j,ii,lnye
 real,intent(in)::fwidth(:),fdepth(:),width(:),depth(:)
 lnye = 1; ig = 0
 do i=1,fnye; do j=1,fnxe
 ig = ig + 1
 g_coord(1,g_num(1,ig))=fwidth(j); g_coord(1,g_num(2,ig))=fwidth(j)
 g_coord(1,g_num(3,ig))=fwidth(j); g_coord(1,g_num(5,ig))=fwidth(j+1)
 g_coord(1,g_num(6,ig))=fwidth(j+1); g_coord(1,g_num(7,ig))=fwidth(j+1)
 g_coord(1,g_num(4,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(5,ig)))
 g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
 g_coord(2,g_num(1,ig))=fdepth(i); g_coord(2,g_num(8,ig))=fdepth(i)
 g_coord(2,g_num(7,ig))=fdepth(i); g_coord(2,g_num(3,ig))=fdepth(i+1)
 g_coord(2,g_num(4,ig))=fdepth(i+1); g_coord(2,g_num(5,ig))=fdepth(i+1)
 g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(3,ig)))
 g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(7,ig))+g_coord(2,g_num(5,ig)))
 end do; end do
 if(itype==1) then
 do ii = 1 , lifts - 1
  do i=1,lnye; do j=1,lnxe-(ii-1)*1
   ig = ig + 1
   if(j==1) then
    g_coord(1,g_num(1,ig))=width((ii-1)+j)
    g_coord(1,g_num(5,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(6,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(7,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(3,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(5,ig)))
    g_coord(1,g_num(2,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(1,ig)))
    g_coord(1,g_num(4,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(5,ig)))
    g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
    g_coord(2,g_num(1,ig))=depth((ii-1)+i)
    g_coord(2,g_num(8,ig))=depth((ii-1)+i)
    g_coord(2,g_num(7,ig))=depth((ii-1)+i)
    g_coord(2,g_num(5,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(3,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(5,ig)))
    g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(3,ig)))
    g_coord(2,g_num(4,ig))=0.5*(g_coord(2,g_num(3,ig))+g_coord(2,g_num(5,ig)))
    g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(5,ig))+g_coord(2,g_num(7,ig)))
   else
    g_coord(1,g_num(1,ig))=width((ii-1)+j)
    g_coord(1,g_num(2,ig))=width((ii-1)+j)
    g_coord(1,g_num(3,ig))=width((ii-1)+j)
    g_coord(1,g_num(5,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(6,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(7,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(4,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(5,ig)))
    g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
    g_coord(2,g_num(1,ig))=depth((ii-1)+i)
    g_coord(2,g_num(8,ig))=depth((ii-1)+i)
    g_coord(2,g_num(7,ig))=depth((ii-1)+i)
    g_coord(2,g_num(3,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(4,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(5,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(3,ig)))
    g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(5,ig))+g_coord(2,g_num(7,ig)))
  end if
  end do; end do
 end do
 else
 do ii = 1 , lifts - 1
  do i=1,lnye; do j=1,lnxe-(ii-1)*1
   ig = ig + 1
   if(j==1) then
    g_coord(1,g_num(1,ig))=width((ii-1)+j)
    g_coord(1,g_num(5,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(6,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(7,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(3,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(4,ig))=g_coord(1,g_num(5,ig))
    g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
    g_coord(1,g_num(2,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(5,ig)))
    g_coord(2,g_num(1,ig))=depth((ii-1)+i)
    g_coord(2,g_num(8,ig))=depth((ii-1)+i)
    g_coord(2,g_num(7,ig))=depth((ii-1)+i)
    g_coord(2,g_num(5,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(5,ig)))
    g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(5,ig))+g_coord(2,g_num(7,ig)))
    g_coord(2,g_num(3,ig))=g_coord(2,g_num(5,ig))
    g_coord(2,g_num(4,ig))=g_coord(2,g_num(5,ig))
   else
    g_coord(1,g_num(1,ig))=width((ii-1)+j)
    g_coord(1,g_num(2,ig))=width((ii-1)+j)
    g_coord(1,g_num(3,ig))=width((ii-1)+j)
    g_coord(1,g_num(5,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(6,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(7,ig))=width((ii-1)+j+1)
    g_coord(1,g_num(4,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(5,ig)))
    g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
    g_coord(2,g_num(1,ig))=depth((ii-1)+i)
    g_coord(2,g_num(8,ig))=depth((ii-1)+i)
    g_coord(2,g_num(7,ig))=depth((ii-1)+i)
    g_coord(2,g_num(3,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(4,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(5,ig))=depth((ii-1)+i+1)
    g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(3,ig)))
    g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(5,ig))+g_coord(2,g_num(7,ig)))
  end if
  end do; end do
 end do
 end if
 return
end subroutine fmcoem                       
 subroutine fmglob(nxe,nye,g_num)
 ! returns g_num for the mesh
 implicit none
 integer,intent(in)::nxe,nye;integer,intent(out)::g_num(:,:)
 integer::ig,i,j       ; ig = 0
 do i=1,nye; do j=1,nxe
  ig = ig + 1
  g_num(1,ig)=(j*2-1)+(i-1)*(nxe+1)+(i-1)*(2*nxe+1)
  g_num(8,ig) = g_num(1,ig)+1; g_num(7,ig) = g_num(1,ig)+2
  g_num(2,ig)=i*(2*nxe+1)+j+(i-1)*(nxe+1)
  g_num(6,ig)=g_num(2,ig)+1
  g_num(3,ig)=i*(2*nxe+1)+i*(nxe+1)+(j*2-1)
  g_num(4,ig)=g_num(3,ig)+1; g_num(5,ig)=g_num(3,ig)+2
 end do; end do
 return
 end subroutine fmglob
 subroutine fmcoco(g_num,g_coord,width,depth,nxe,nye)
 ! returns g_coord for the mesh
 implicit none
 integer,intent(in)::g_num(:,:),nxe,nye
 real,intent(out)::g_coord(:,:);integer::ig,i,j
 real,intent(in)::width(:),depth(:)
 ig = 0
 do i=1,nye; do j=1,nxe
 ig = ig + 1
 g_coord(1,g_num(1,ig))=width(j); g_coord(1,g_num(2,ig))=width(j)
 g_coord(1,g_num(3,ig))=width(j); g_coord(1,g_num(5,ig))=width(j+1)
 g_coord(1,g_num(6,ig))=width(j+1); g_coord(1,g_num(7,ig))=width(j+1)
 g_coord(1,g_num(4,ig))=0.5*(g_coord(1,g_num(3,ig))+g_coord(1,g_num(5,ig)))
 g_coord(1,g_num(8,ig))=0.5*(g_coord(1,g_num(1,ig))+g_coord(1,g_num(7,ig)))
 g_coord(2,g_num(1,ig))=depth(i); g_coord(2,g_num(8,ig))=depth(i)
 g_coord(2,g_num(7,ig))=depth(i); g_coord(2,g_num(3,ig))=depth(i+1)
 g_coord(2,g_num(4,ig))=depth(i+1); g_coord(2,g_num(5,ig))=depth(i+1)
 g_coord(2,g_num(2,ig))=0.5*(g_coord(2,g_num(1,ig))+g_coord(2,g_num(3,ig)))
 g_coord(2,g_num(6,ig))=0.5*(g_coord(2,g_num(7,ig))+g_coord(2,g_num(5,ig)))
 end do; end do
 return
end subroutine fmcoco 
end module geometry_lib