¿ÉÖ´ÐгÌÐò Planar_8_Nodes.f90

!=============================================================
! ƽÃæ8½ÚµãµÈ²ÎÔªÍêÕû³ÌÐò
! Ç廪´óѧÍÁľ¹¤³Ìϵ£ºÂ½ÐÂÕ÷
£¡³ÌÐò˵Ã÷
!============================================================

module Elem_Rect8 ! °Ë½ÚµãµÈ²ÎÔª
implicit none
integer (kind(1)),parameter ::ikind=(kind(1))
integer (kind(1)),parameter ::rkind=(kind(0.d0))

type :: typ_Kcol
real(rkind),pointer :: Row(:)
end type typ_Kcol

type :: typ_GValue !×ÜÌå¿ØÖƱäÁ¿
integer(ikind) :: NNode, NElem, NLoad, NMat, NSupport
integer(ikind) :: NGlbDOF !ÕûÌå×ÔÓɶÈ×ÜÊý
integer(ikind) :: NGENS, NodeDOF,ElemNodeNo
integer(ikind) :: NInt
end type typ_GValue

type Typ_Node !¶¨Òå½ÚµãÀàÐÍ
real(rkind) :: coord(2) !½Úµã×ø±ê
integer(ikind) :: GDOF(2) !ÕûÌå×ÔÓɶȱàÂë
real(rkind) :: DISP(2) !½ÚµãλÒÆ
real(rkind) :: dDISP(2) !½ÚµãλÒÆÔöÁ¿
real(rkind) :: dForce(2) !½Úµã²»Æ½ºâÁ¦

end type typ_Node

!=============================================================================
Type typ_IntPoint !¶¨Òå»ý·Öµã²ÎÊý
real(rkind) :: EPS(3) !Ó¦±ä
real(rkind) :: SIG(3) !Ó¦Á¦
real(rkind) :: D(3,3) !±¾¹¹¾ØÕó
real(rkind) :: B(3,16) !¼¸ºÎ¾ØÕó
real(rkind) :: DETJ !Ñſ˱ÈÐÐÁÐʽ
end type Typ_IntPoint

type Typ_Rect8 !¶¨ÒåʵÌåµ¥Ôª
integer(ikind) :: NodeNo(8) !½Úµã±àºÅ
real(rkind) :: E !µ¯ÐÔÄ£Á¿
real(rkind) :: u !²´ËɱÈ
real(rkind) :: t !µ¥Ôªºñ¶È
real(rkind) :: EK(16,16) !µ¥Ôª¸Õ¶È¾ØÕó
type(typ_intpoint) :: IntP(9) !»ý·Öµã
!........................
end type Typ_Rect8

type Typ_Load ! ¶¨ÒåºÉÔØÀàÐÍ
integer(ikind) :: NodeNo !ºÉÔØ×÷Óýڵã±àºÅ
real(rkind) :: Value(2) !ºÉÔØ´óС
end type Typ_Load

type Typ_Support ! ¶¨ÒåÖ§×ùÀàÐÍ
integer(ikind) :: NodeNo !Ö§×ù½Úµã±àºÅ
integer(ikind) :: DOF !Ö§×ùÔ¼Êø×ÔÓɶÈ
real(rkind) :: Value !Ö§×ùλÒÆ´óС
end type Typ_Support

contains

subroutine Get_Elem_K(GValue,Elem,Node) ! ºËÐijÌÐò£¬µÃµ½µ¥ÔªµÄ¸Õ¶È¾ØÕó
implicit none
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
type(typ_Rect8) :: Elem(:)
real*8 :: InvJ(2,2) ! Ñſ˱ȾØÕó¼°ÆäÄæ¾ØÕó
real*8 :: Coord(8,2),IntPCoord(9,2),IntPwt(9)
real*8 :: dN(2,8) !½Úµã×ø±ê£¬»ý·Öµã×ø±ê£¬È¨º¯Êý£¬Ðκ¯Êýµ¼Êý
integer :: I,ElemNo

call Get_IntP_Prop(IntPCoord,IntPWt) ! ¼ÆËã»ý·ÖµãÐÅÏ¢

do ElemNo=1,size(Elem)
Elem(ElemNo)%EK=0.0
do I=1, GValue%ElemNodeNo; !µÃµ½½Úµã×ø±ê
coord(I,:)=Node(Elem(ElemNo)%NodeNo(i))%coord;
end do
Elem(ElemNo)%EK=0d0
DO I=1,size(Elem(ElemNo)%IntP)
! ¼ÆËã±¾¹¹¾ØÕó
call Get_D(Elem(ElemNo)%IntP(I)%D,Elem(ElemNo)%E,Elem(ElemNo)%u)
! ¼ÆËãÐκ¯Êý¶Ô¾Ö²¿×ø±êµ¼Êý
call Get_dN_dxi(dN,IntPCoord(i,1),IntPCoord(i,2))
InvJ=matmul(dN, Coord)
! µÃµ½Ñſ˱ÈÐÐÁÐʽ
Elem(ElemNo)%IntP(I)%DETJ=InvJ(1,1)*InvJ(2,2)-InvJ(1,2)*InvJ(2,1)
InvJ=matinv(InvJ); ! ¶ÔÑſ˱ÈÐÐÁÐʽÇóÄæ
dN = matmul(InvJ,dN) ! µÃµ½Ðκ¯Êý¶ÔÕûÌå×ø±êµÄµ¼Êý
call Get_B (Elem(ElemNo)%IntP(I)%B,dN) ! µÃµ½¼¸ºÎ¾ØÕó
Elem(ElemNo)%EK=Elem(ElemNo)%EK+&
matmul(matmul(transpose(Elem(ElemNo)%IntP(I)%B),&
Elem(ElemNo)%IntP(I)%D),Elem(ElemNo)%IntP(I)%B)*&
Elem(ElemNo)%IntP(I)%DetJ*IntPWt(i)*Elem(ElemNo)%t

end do
end do
return
end subroutine

subroutine Get_D(D,e,v) ! µÃµ½±¾¹¹¾ØÕó£¬Æ½ÃæÓ¦Á¦
implicit none
real*8 :: e,v
real*8 :: D (:,:)
D=0.d0
D(1,1)=1D0; D(2,2)=1d0
D(1,2)=v; D(2,1)=v
D(3,1:2)=0d0; D(1:2,3)=0d0;
D(3,3)=(1.d0-v)/2.d0;
D=E/(1.d0-v*v)*D;
return
end subroutine Get_D

subroutine Get_IntP_Prop(IntPCoord,IntPWt) ! µÃµ½¸ß˹»ý·ÖµãµÄ²ÎÊý
implicit none
real*8 ::IntPCoord(:,:),IntPWt(:) ! ¸ß˹»ý·Öµã×ø±ê£¬¸ß˹»ý·ÖµãȨÖØ
real*8 :: z,x(3),y(3),w(3)
integer :: i,j
z=.2d0*sqrt(15.d0)
x=(/-1.d0,0.d0,1.d0/); y=(/1.d0,0.d0,-1.d0/)
w=(/5.d0/9.d0,8.d0/9.d0,5.d0/9.d0/)
do i=1,3
do j=1,3
IntPCoord((i-1)*3+j,1)=x(j)*z
IntPCoord((i-1)*3+j,2)=y(i)*z
IntPWt((i-1)*3+j)=w(i)*w(j)
end do
end do

return
end subroutine Get_IntP_Prop

subroutine Get_dN_dxi(dN_dxi,xi,eta) ! µÃµ½Ðκ¯Êý¶Ô¾Ö²¿×ø±êϵµÄµ¼Êý
implicit none
real*8 :: dN_dxi(:,:), xi,eta
real*8:: x1, x2, x3, x4, x5 ,x6 ,x7 ,x8
x1=.25*(1.-eta);
x2=.25*(1.+eta);
x3=.25*(1.-xi);
x4=.25*(1.+xi)

dN_dxi(1,1)=x1*(2.*xi+eta)
dN_dxi(1,2)=-8.*x1*x2
dN_dxi(1,3)=x2*(2.*xi-eta)
dN_dxi(1,4)=-4.*x2*xi
dN_dxi(1,5)=x2*(2.*xi+eta)
dN_dxi(1,6)=8.*x2*x1
dN_dxi(1,7)=x1*(2.*xi-eta)
dN_dxi(1,8)=-4.*x1*xi
dN_dxi(2,1)=x3*(xi+2.*eta)
dN_dxi(2,2)=-4.*x3*eta
dN_dxi(2,3)=x3*(2.*eta-xi)
dN_dxi(2,4)=8.*x3*x4
dN_dxi(2,5)=x4*(xi+2.*eta)
dN_dxi(2,6)=-4.*x4*eta
dN_dxi(2,7)=x4*(2.*eta-xi)
dN_dxi(2,8)=-8.*x3*x4
return
end subroutine Get_dN_dxi

subroutine Get_B(B,dN_dx) ! µÃµ½¼¸ºÎ¾ØÕó
implicit none
real*8 :: B(:,:), dN_dx(:,:)
integer::k,l,m,n , ih,nod;
real*8 :: x,y,z
B=0.
do m=1,8
k=2*m; l=k-1; x=dN_dx(1,m); y=dN_dx(2,m)
B(1,l)=x; B(3,k)=x; B(2,k)=y; B(3,l)=y
end do
return
end subroutine Get_B

function matinv(A) result (B) ! ¼ÆËã2x2Äæ¾ØÕó
real(rkind) ,intent (in):: A(:,:)
real(rkind) , pointer::B(:,:)
real(rkind) :: x
integer :: N
N=size(A,dim=2)
allocate(B(N,N))
B(1,1)=A(2,2)
B(2,2)=A(1,1)
B(1,2)=-A(1,2)
B(2,1)=-A(2,1)
x=A(1,1)*A(2,2)-A(1,2)*A(2,1)
B=B/x
return
end function matinv

end module

module Data_Input ! Êý¾ÝÊäÈëÄ£¿é
use Elem_Rect8
implicit none

contains

subroutine DataInput(GValue,Node,Elem,Load,Support)
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)

call ReadGValue(GValue) !¶ÁÈëÕûÌå¿ØÖƱäÁ¿
call ReadNode(GValue,Node) !¶ÁÈë½Úµã×ø±êÐÅÏ¢
call ReadElem(GValue,Elem) !¶ÁÈ뵥ԪԭʼÐÅÏ¢
call ReadMaterial(GValue,Elem)
call ReadLoad(GValue,Load) !¶ÁÈëºÉÔØÐÅÏ¢
call ReadSupport(GValue,Support) !¶ÁÈëÖ§×ùÐÅÏ¢
return
end subroutine DataInput

subroutine ReadGValue(GValue)
type(typ_GValue) :: GValue
GValue%NGENS=3
GValue%NodeDOF=2
GValue%ElemNodeNo=8
GValue%NInt=3
return
end subroutine ReadGValue

subroutine ReadNode(GValue,Node) ! ¶ÁÈë½áµã×ø±ê
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
integer(ikind) :: I,J
open(55,file='node.txt')
read(55,*) GValue%NNode
GValue%NGlbDOF=2*GValue%NNode
allocate(Node(GValue%NNode))
do I=1, GValue%NNode
read(55,*) J,Node(I)%Coord(1:2)
end do
close(55)
write(*,*) " ½áµãÐÅÏ¢¶ÁÈëÍê±Ï"
return
end subroutine ReadNode

subroutine ReadElem(GValue,Elem) ! ¶ÁÈëµ¥ÔªÐÅÏ¢
type(typ_GValue) :: GValue
type(typ_Rect8),pointer :: Elem(:)
integer(ikind) :: I,J,K
integer(ikind) :: TransNode(8),N(8),N1(8) ! ½áµã×ø±ê±ä»»
TransNode=(/1,7,5,3,8, 6,4, 2/) ! ×¢Òâµ¥Ôª½Úµã±àºÅÆ¥Åä
open(55,file='Element.txt')
read(55,*) GValue%NElem
allocate(Elem(GValue%NElem))
do I=1,GValue%NElem
read(55,*) J, N
do J=1,size(N)
K=TransNode(J)
N1(K)=N(J)
end do
Elem(I)%NodeNo=N1
end do
close(55)
write(*,*) " µ¥ÔªÐÅÏ¢¶ÁÈëÍê±Ï"
return
end subroutine ReadElem

subroutine ReadMaterial(GValue,Elem) ! ¶ÁÈë²ÄÁÏÐÅÏ¢
type(typ_GValue) :: GValue
type(typ_Rect8),pointer :: Elem(:)
integer (ikind) :: NElem
real(rkind) :: E, mu
integer(ikind) :: I,J
integer(ikind),allocatable :: K(:)
open(55,file='Material.txt')
read(55,*) E, mu
Elem(:)%E=E
Elem(:)%u=mu
Elem(:)%t=0.5
close(55)
write(*,*) " ²ÄÁÏÐÅÏ¢¶ÁÈëÍê±Ï"
return
end subroutine ReadMaterial

subroutine ReadLoad(GValue,Load) ! ¶ÁÈëºÉÔØÐÅÏ¢
type(typ_GValue) :: GValue
type(typ_Load),pointer :: Load(:)
integer(ikind) :: I,J
open(55,file='Load.txt')
read(55,*) GValue%NLoad
allocate(Load(GValue%NLoad))
do I=1, GValue%NLoad
read(55,*) J, Load(I)%NodeNo,Load(I)%Value(1:2)
end do
close(55)
write(*,*) " ºÉÔØÐÅÏ¢¶ÁÈëÍê±Ï"
return
end subroutine ReadLoad

subroutine ReadSupport(GValue,Support) ! ¶ÁÈëÖ§×ùÐÅÏ¢
type(typ_GValue) :: GValue
type(typ_Support),pointer :: Support(:)
integer(ikind) :: I,J
open(55,file='Support.txt')
read(55,*) GValue%NSupport
allocate(Support(GValue%NSupport))
do I=1, GValue%NSupport
read(55,*) J, Support(I)%NodeNo, Support(I)%DOF, Support(I)%Value
end do
close(55)
write(*,*) " Ö§×ùÐÅÏ¢¶ÁÈëÍê±"
return
end subroutine ReadSupport

end module

module Mat_solve ! ¾ØÕóÇó½âÄ£¿é
use Elem_Rect8
implicit none

integer :: NGlbDOF
contains

subroutine Matsolve(GValue,Elem,Node,Load, Support)
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
type(typ_Kcol),allocatable :: Kcol(:)
real(rkind),allocatable :: GLoad(:), GDisp(:)
real(rkind) :: Penatly
integer(ikind) :: BandWidth(2*GValue%NNode)
integer(ikind) :: ELocVec(16)
integer(ikind) :: I,J,K
integer(ikind) :: MinDOFNum

NGlbDOF=2*GValue%NNode
Penatly=1.0 ! ·£º¯Êý

allocate(GLoad(NGlbDOF))
allocate(GDisp(NGlbDOF))

!²éÕÒ´ø¿í
do I=1,NGlbDOF
BandWidth(I)=I
end do
do I=1,GValue%NElem
do J=1,size(Elem(I)%NodeNo)
ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1
ELocVec(J*2 )=2*Elem(I)%NodeNo(J)
end do
MinDOFNum=minval(ELocVec)
do J=1,size(ElocVec)
BandWidth(ELocVec(J))=min(MinDOFNum,BandWidth(ELocVec(J)))
end do
end do
write(*,*) " Íê³É´ø¿í²éÕÒ"
allocate(Kcol(NGlbDOF))
do I=1, NGlbDOF
allocate(Kcol(I)%Row(BandWidth(I):I))
Kcol(I)%Row=0.d0
end do
do I=1, GValue%NElem
do J=1,size(Elem(I)%NodeNo)
ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1
ELocVec(J*2 )=2*Elem(I)%NodeNo(J)
end do
do J=1,size(ElocVec)
do K=1,J
if(ELocVec(J)>ELocVec(K)) then
Kcol(ELocVec(J))%row(ELocVec(K))=&
Kcol(ELocVec(J))%row(ELocVec(K))+Elem(I)%EK(J,K)
else
Kcol(ELocVec(K))%row(ELocVec(J))=&
Kcol(ELocVec(K))%row(ELocVec(J))+Elem(I)%EK(J,K)
end if
end do
end do
Penatly=max(Penatly,maxval(abs(Elem(I)%EK)))
end do

write(*,*) " Íê³É×ܸռ¯³É"

GLoad=0.d0

do I=1, size(Load)
J=Load(I)%NodeNo
GLoad(J*2-1:J*2)=GLoad(J*2-1:J*2)+Load(I)%Value(:)
end do

do I=1,GValue%NSupport
J=2*(Support(I)%NodeNo-1)+Support(I)%DOF
GLoad(J)=GLoad(J)+Support(I)%Value*Penatly*1.0D8
Kcol(J)%Row(J)=KCol(J)%Row(J)+Penatly*1.0d8
end do
write(*,*) " Íê³ÉÖ§×ù¼¯³É"

call BandSolv(Kcol,Gload,GDisp)

call Get_Node_dDisp(GValue,Node,GDisp)
do I=1,NGlbDOF
deallocate(Kcol(I)%Row)
end do
deallocate(Kcol)
deallocate(GDisp)
deallocate(GLoad)
return
end subroutine

subroutine Get_Node_dDisp(GValue,Node,GDisp) ! ´ÓÕûÌåλÒÆÏòÁ¿µÃµ½½áµãλÒÆÏòÁ¿
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
real(rkind) :: GDisp(:)
integer(ikind) :: I,J
do I=1, GValue%NNode
Node(I)%dDisp(1:2)=GDisp(I*2-1:I*2)
Node(I)%Disp=Node(I)%Disp+Node(I)%dDisp
end do
return
end subroutine Get_Node_dDisp

!---------------------------------------
subroutine BandSolv(Kcol,GLoad,diag)
!---------------------------------------
type (typ_Kcol) :: Kcol(:);
real(rkind) :: GLoad(:),diag(:);
integer :: row1,ncol,row,j,ie
real(rkind) :: s
ncol=NGlbDOF
diag(1:ncol)=(/(Kcol(j)%row(j),j=1,ncol)/)
do j=2,ncol
row1=lbound(Kcol(j)%row,1)
do ie=row1,j-1
row=max(row1,lbound(Kcol(ie)%row,1))
s=sum(diag(row:ie-1)*Kcol(ie)%row(row:ie-1)*Kcol(j)%row(row:ie-1))
Kcol(j)%row(ie)=(Kcol(j)%row(ie)-s)/diag(ie)
end do
s=sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)
diag(j)=diag(j)-s
end do

do ie=2,ncol
row1=lbound(Kcol(ie)%row,dim=1)
GLoad(ie)=GLoad(ie)-sum(Kcol(ie)%row(row1:ie-1)*GLoad(row1:ie-1))
end do
GLoad(:)=GLoad(:)/diag(:)

do j=ncol,2,-1
row1=lbound(Kcol(j)%row,dim=1)
GLoad(row1:j-1)=GLoad(row1:j-1)-GLoad(j)*Kcol(j)%row(row1:j-1)
end do
diag(:)=GLoad(:)

return;
end subroutine
end module

module Data_Output ! Êý¾ÝÊä³öÄ£¿é
use Elem_Rect8
implicit none

contains

subroutine DataOutput(GValue,Elem,Node,Load, Support) ! Êä³ö½ÚµãλÒÆ
type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)
integer i
open (55,file='out.txt')
do i=1,size(Node)
write(55,'(I5, 2G14.5)') i, Node(i)%Disp(1:2)
end do
close(55)
return
end subroutine

end module

Program Main ! Ö÷³ÌÐò
use Elem_Rect8
use Data_Input
use Mat_Solve
use Data_Output
implicit none

type(typ_GValue) :: GValue
type(typ_Node),pointer :: Node(:)
type(typ_Rect8),pointer :: Elem(:)
type(typ_Load),pointer :: Load(:)
type(typ_Support),pointer :: Support(:)

write(*,*) "¿ªÊ¼¶ÁÈ¡Êý¾Ý"
call DataInput(GValue,Node,Elem,Load,Support) ! ¶ÁÈëÊý¾Ý
write(*,*) "¿ªÊ¼¼ÆËãµ¥Ôª¸Õ¶È¾ØÕó"
call Get_Elem_K(GValue,Elem,Node) ! µÃµ½µ¥Ôª¸Õ¶È¾ØÕó
write(*,*) "¿ªÊ¼Çó½âÕûÌå·½³Ì"
call MatSolve(GValue,Elem,Node,Load, Support) ! Çó½âÕûÌåƽºâ·½³Ì
write(*,*) "¿ªÊ¼Êä³ö½á¹û"
call DataOutput(GValue,Elem,Node,Load, Support) ! Êä³ö½á¹û
stop
end program